solvetransition.m (3730B)
1 %%% A function that computes the path for the system over the specified time interval. It finds the "best" s_0 and l_0 given the other parameters provided, and returns the transition corresponding to that. It uses getels0.m to get the best s_0 and l_0 2 3 %%% Basic routine to solve for transition dynamics of the model. 4 %%% Used in Transition.m. Uses transit1dx.m as a subroutine. 5 %%% Show solution if ShowResults==1 6 7 function [t,x,chat,hhat,gdpgrowth,shat,ellhat, dltahat, sigmahat]=solvetransition(dltabar,ubar,epsilon,beta,gamma,dlta0, Nend,alpha,lambda,phi,rho,nbar,T,tstep,ShowResults,withShock) 8 9 10 %% Define variables 11 [sstar, ellstar, sigmastar, dltastar, ystar, zstar, gs, gc, gh, gdelta] = getsteadystate(dltabar,ubar,epsilon,beta,gamma,dlta0,Nend,alpha,lambda,phi,rho,nbar); 12 13 tdv=0; 14 %%% tinydev=1e-5; % s ell sigma dlta y z 15 %%% tdv=tinydev; 16 17 sguess=.005; 18 ellguess=.01; 19 sellguess=[sguess ellguess]; 20 sell0=getells0(sellguess,dlta0,Nend,alpha,epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,zstar,ystar,sigmastar,nbar,gs,T); 21 sstart=sell0(1); 22 ellstart=sell0(2); 23 24 xT=[sstart ellstart sigmastar dlta0 ystar-tdv zstar+tdv Nend]'; 25 26 %% ODE : Do the solving 27 %%% Unclear where exactly the solving happens. 28 tspan=[T:-tstep:0]; 29 options = odeset('Stats','off','RelTol',1e-6); 30 31 [t,x] = ode23t(@(t,x) transit1dx(t,x,alpha,epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,nbar,withShock),tspan,xT,options); 32 33 %% Recover the key variables 34 s=x(:,1); 35 ell=x(:,2); 36 sigma=x(:,3); 37 dlta=x(:,4); 38 y=x(:,5); 39 z=x(:,6); 40 N=x(:,7); 41 42 %% Recover terms that show up in many equations 43 ll=ell./(1-ell); 44 sigsig=sigma./(1-sigma); 45 ss=s./(1-s); 46 47 %% utilde=ubar*c.^(gamma-1)+1/(1-gamma); % u(c)/u'(c)c = value lifeyear / c 48 49 %%% For loop. What exactly is happening here? 50 51 for i=1:length(t); 52 dx(:,i)=transit1dx(t(i),x(i,:),alpha,epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,nbar,withShock); 53 end; 54 55 56 %% Recover more variables 57 58 dx=dx'; %%% [s ell sigma dlta y z] 59 shat=dx(:,1)./s; 60 ellhat=dx(:,2)./ell; 61 sigmahat=dx(:,3)./sigma; 62 dltahat=dx(:,4)./dlta; 63 yhat=dx(:,5)./y; 64 zhat=dx(:,6)./z; 65 chat=alpha*y+ellhat-sigmahat.*sigsig; 66 hhat=alpha*z-ellhat.*ll-sigmahat.*sigsig; 67 %%% vtilde=1./(beta*ll.*dlta); 68 gdpgrowth=(1-ell).*hhat+ell.*chat; 69 70 AoverB=(z./y).^(1/(1-phi)).*ss.^(lambda/(1-phi)); 71 coverh=(AoverB.^alpha) .*ll; 72 c=(((dlta./dltabar).*(coverh.^(-beta))).^(1/(epsilon-beta)))./N; 73 utilde=ubar.*(c.^(gamma-1))+1/(1-gamma); 74 %%% u(c)/u'(c)c 75 uoverv =ll.*beta.*dlta.*utilde+epsilon.*dlta.*utilde; 76 77 gvt = rho - uoverv + dlta; 78 vtilde = utilde ./(rho - dlta + gvt); 79 dltavtilde = dlta.*vtilde; 80 81 %% Pretty print 82 if ShowResults; 83 %%% Show data 84 fmt='%6.0f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.1e %8.4f %8.4f %8.4f %8.4f %8.4f %8.2f %8.2f %8.4f'; 85 cshow(' ',[t x [chat hhat gdpgrowth]*100 shat ellhat utilde vtilde dltavtilde],fmt,'Time s ell sigma dlta gA gB N chat hhat gdphat shat ellhat utilde vtilde deltavtilde'); 86 87 %%% Checks on initial values -- for close to SS 88 disp ' '; 89 90 disp 'Checks on initial values to ensure we are close to SS:'; 91 %%% fprintf('beta*dlta*utilde*ell/(1-ell) at T should be close to rho=%4.2f: %8.6f\n',[rho beta*dlta(1)*utilde(1)*ell(1)/(1-ell(1))]); 92 fprintf('gs=%6.4f shat=%6.4f ellhat=%6.4f\n',[gs shat(1) ellhat(1)]); 93 fprintf('sigma=%6.4f sigmastar=%6.4f sigmahat=%6.4f sigmahatstar=%6.4f \n',[sigma(1) sigmastar sigmahat(1) 0]); 94 fprintf('y=%6.4f ystar=%6.4f z=%6.4f zstar=%6.4f\n',[y(1) ystar z(1) zstar]); 95 fprintf('yhat=%6.4f yhatstar=0 zhat=%6.4f zhatstar=0\n',[yhat(1) zhat(1)]); 96 fprintf('dltahat=%6.4f dltahatstar=%6.4f\n',[dltahat(1) gdelta]); 97 98 end 99 100 %%% Matlab doesn't seem to use the return keyword. Instead, return variables are defined with the function definition.