% simul.m
%
% this program solves the standard growth model in four different ways (in
% levels (1st & 2nd order) and in loglevels (1st & 2nd order). 
%
% for each solution it simulates capital and consumption levels and stores
% them in a matrix and finally plots them
%

T=1000;
capstore=zeros(T,4);
constore=zeros(T,4);

%indicate which column of decision corresponds with which variable
%note that for each application you have to check that these column
%indications are correct. Here the labels are the same for each program
iiicap=1;
iiicons=3;

%initialize stuff
%this is the same for each numerical solution
cap    = zeros(T,1);
con    = zeros(T,1);
prod   = zeros(T,1);	
randn('state',20080701);
sig = 0.01; % this value of sigma has to equal the one given in the Dynare programs (for 2nd-order)
shocks = sig*randn(T,1); 
prod(1)= 0;


dynare modelclevel noclearall
% the command noclearall means that variables defined here (like T) are not
% deleted out of memory
load dynarerocks
% this loads the matrix "decision" with coefficients into memory

capstst = decision(1,iiicap); 	
cap(1)  = capstst; % initial value is steady state level
for t = 2:T	
prod(t)=rho*prod(t-1)+shocks(t);
cap(t) =decision(1,iiicap)+decision(2,iiicap)*(cap(t-1)-capstst) ...	
           +decision(3,iiicap)*prod(t-1) + decision(4,iiicap)*shocks(t);
con(t) =decision(1,iiicons)+decision(2,iiicons)*(cap(t-1)-capstst) ...	
           +decision(3,iiicons)*prod(t-1) + decision(4,iiicons)*shocks(t);
end	

capstore(:,1)=cap;
constore(:,1)=con;

dynare modelcloglinear noclearall
% now we get a policy rule for the LOG of capital and the LOG of consumption
% note that productivity is always in logs
load dynarerocks
% this loads the matrix "decision" with coefficients into memory

capstst = decision(1,iiicap); 	
cap(1)  = capstst; % initial value is steady state level
for t = 2:T	
prod(t)=rho*prod(t-1)+shocks(t);
cap(t) =decision(1,iiicap)+decision(2,iiicap)*(cap(t-1)-capstst) ...	
           +decision(3,iiicap)*prod(t-1) + decision(4,iiicap)*shocks(t);
con(t) =decision(1,iiicons)+decision(2,iiicons)*(cap(t-1)-capstst) ...	
           +decision(3,iiicons)*prod(t-1) + decision(4,iiicons)*shocks(t);
end	

capstore(:,2)=exp(cap);
constore(:,2)=exp(con);
% since this solution is generating logs I have to exp it to get the level

dynare modelclevel2ndorder noclearall
% the command noclearall means that variables defined here (like T) are not
% deleted out of memory
load dynarerocks
% this loads the matrix "decision" with coefficients into memory

capstst = decision(1,iiicap)-decision(2,iiicap); %recall that second coef gives correction from st.state 	
cap(1)  = capstst; % initial value is steady state level
for t = 2:T	
prod(t)=rho*prod(t-1)+shocks(t);
cap(t)=decision(1,iiicap) ...
      +decision(3,iiicap)*(cap(t-1)-capstst) ...
      +decision(4,iiicap)*prod(t-1) ...
      +decision(5,iiicap)*shocks(t) ...
      +decision(6,iiicap)*(cap(t-1)-capstst)^2 ...
      +decision(7,iiicap)*(cap(t-1)-capstst)*prod(t-1) ...
      +decision(8,iiicap)*prod(t-1)^2 ...
      +decision(9,iiicap)*shocks(t)^2 ...
      +decision(10,iiicap)*(cap(t-1)-capstst)*shocks(t) ...
      +decision(11,iiicap)*prod(t-1)*shocks(t);
con(t)=decision( 1,iiicons) ...
      +decision( 3,iiicons)*(cap(t-1)-capstst) ...
      +decision( 4,iiicons)*prod(t-1) ...
      +decision( 5,iiicons)*shocks(t) ...
      +decision( 6,iiicons)*(cap(t-1)-capstst)^2 ...
      +decision( 7,iiicons)*(cap(t-1)-capstst)*prod(t-1) ...
      +decision( 8,iiicons)*prod(t-1)^2 ...
      +decision( 9,iiicons)*shocks(t)^2 ...
      +decision(10,iiicons)*(cap(t-1)-capstst)*shocks(t) ...
      +decision(11,iiicons)*prod(t-1)*shocks(t);
end	

capstore(:,3)=cap;
constore(:,3)=con;

dynare modelclog2ndorder noclearall
% the command noclearall means that variables defined here (like T) are not
% deleted out of memory
load dynarerocks
% this loads the matrix "decision" with coefficients into memory

capstst = decision(1,iiicap)-decision(2,iiicap); %recall that second coef gives correction from st.state 	
cap(1)  = capstst; % initial value is steady state level
for t = 2:T	
prod(t)=rho*prod(t-1)+shocks(t);
cap(t)=decision(1,iiicap) ...
      +decision(3,iiicap)*(cap(t-1)-capstst) ...
      +decision(4,iiicap)*prod(t-1) ...
      +decision(5,iiicap)*shocks(t) ...
      +decision(6,iiicap)*(cap(t-1)-capstst)^2 ...
      +decision(7,iiicap)*(cap(t-1)-capstst)*prod(t-1) ...
      +decision(8,iiicap)*prod(t-1)^2 ...
      +decision(9,iiicap)*shocks(t)^2 ...
      +decision(10,iiicap)*(cap(t-1)-capstst)*shocks(t) ...
      +decision(11,iiicap)*prod(t-1)*shocks(t);
con(t)=decision( 1,iiicons) ...
      +decision( 3,iiicons)*(cap(t-1)-capstst) ...
      +decision( 4,iiicons)*prod(t-1) ...
      +decision( 5,iiicons)*shocks(t) ...
      +decision( 6,iiicons)*(cap(t-1)-capstst)^2 ...
      +decision( 7,iiicons)*(cap(t-1)-capstst)*prod(t-1) ...
      +decision( 8,iiicons)*prod(t-1)^2 ...
      +decision( 9,iiicons)*shocks(t)^2 ...
      +decision(10,iiicons)*(cap(t-1)-capstst)*shocks(t) ...
      +decision(11,iiicons)*prod(t-1)*shocks(t);
end	

capstore(:,4)=exp(cap);
constore(:,4)=exp(con);
% since this solution is generating logs I have to exp it to get the level

plot(capstore)
pause
plot(constore(2:end,:))
pause
plot(capstore(101:200,:))
pause
plot(constore(101:200,:))

