transit1dx.m (2711B)
1 %%% A function representing a system of 6 ODEs for the time evolution of the model. Specifically, I think it gives the changes in the time-dependent variables corresponding to the current values and step-size 2 3 4 %%% 07/28/19 5 %%% The system of 6 differential equations for the life & growth model 6 %%% x=[s,ell,sigma,dlta,y,z, N] y==gA z==gB 7 %%% Used in Transition.m 8 9 function dx=transit1dx(t,x,alpha,epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,nbar, withShock); 10 11 s=x(1); 12 ell=x(2); 13 sigma=x(3); 14 dlta=x(4); 15 y=x(5); 16 z=x(6); 17 N=x(7); 18 19 20 if withShock==1 21 shockYear=860; 22 23 if (t>=shockYear & t<(shockYear+30)) 24 nbar=0.02; 25 end 26 27 %%% if (t>=(shockYear+20) & t<(shockYear+40)) 28 %%% nbar=0; 29 %%% end 30 end 31 32 if withShock==2 33 shockYear=840; 34 35 if (t>=shockYear & t<(shockYear+40)) 36 nbar=0.02; 37 end 38 39 if (t>=(shockYear+40) & t<(shockYear+80)) 40 nbar=0; 41 end 42 43 %%% if (t>=(shockYear+20) & t<(shockYear+40)) 44 %%% nbar=0; 45 %%% end 46 end 47 48 49 %% Terms that show up in many equations 50 ll=ell./(1-ell); 51 ebl=(epsilon./beta)./ll; 52 sigsig=sigma./(1-sigma); 53 ss=s./(1-s); 54 55 AoverB=(z./y).^(1/(1-phi)).*ss.^(lambda/(1-phi)); 56 coverh=(AoverB.^alpha) .*ll; 57 c=(((dlta./dltabar).*(coverh.^(-beta))).^(1/(epsilon-beta)))./N; 58 utilde=ubar.*(c.^(gamma-1))+1/(1-gamma); % u(c)/u'(c)c 59 uoverv =ll.*beta.*dlta.*utilde+epsilon.*dlta.*utilde; 60 61 %% shat (¿s hat?) first because it makes everything else easier 62 shat=alpha.*z.*(lambda/(1-lambda)).*(1-ell)./sigsig -alpha.*y.*(lambda/(1-lambda)).*(ell)./ss./sigsig; 63 64 %% Terms special to ellhat and sigmahat solution 65 thetaellblah=(1-ell).*(1+ebl); 66 thetaell= thetaellblah./(1 + (gamma - 1 + epsilon + beta.*ll).*thetaellblah); 67 68 thetasig=(1-sigma)./(1-(1-sigma).*lambda - epsilon.*sigma + beta.*sigma); 69 70 71 omegaell=(gamma-1-beta+epsilon).*sigsig; 72 omegasig=-(epsilon + beta.*ll + ll); 73 Acirc = (beta-epsilon)*nbar+(1 -gamma-epsilon).*alpha.*y+beta.*alpha.*z-rho-dlta + uoverv; 74 Bcirc=(1-lambda).*shat.*ss+(lambda+beta-epsilon).*nbar + alpha.*beta.*z-alpha.*epsilon.*y+ uoverv - alpha.*lambda.*z.*(1-ell)./(1-s)./sigsig; 75 76 77 ellhat=thetaell.*(Acirc+thetasig.*omegaell.*Bcirc)./(1-omegaell.*omegasig.*thetasig); 78 sigmahat=thetasig.*(Bcirc+omegasig.*ellhat); 79 80 dltahat=(epsilon-beta).*(nbar - sigmahat.*sigsig) + epsilon.*alpha.*y - beta.*alpha.*z + ellhat.*(epsilon + beta.*ll); 81 82 yhat=lambda.*(nbar+shat+sigmahat)-(1-phi).*y; 83 zhat=lambda.*(nbar-shat.*ss+sigmahat)-(1-phi).*z; 84 85 %% From growth rates back to changes... 86 ds=shat.*s; 87 dell=ellhat.*ell; 88 dsigma=sigmahat.*sigma; 89 ddlta=dltahat.*dlta; 90 dy=yhat.*y; 91 dz=zhat.*z; 92 93 dN=nbar.*N; 94 95 dx=[ds dell dsigma ddlta dy dz dN]'; 96 97 %%% keyboard 98 %%%% Unclear what "keyboard" here refers to. 99 end