% riskaversion.m
%
% this program solves the standard growth model sequentially for different values of riskaversion
% it compares the simulated series and IRFs
% each time a first-order log approximation is used

isimul = 0; % plot simulation pics if 1
iIRF = 1;   % plot IRF pics if 1

T=200;  % observations used in simulations
T2=100; % number of periods for IRFs

capstore=zeros(T,10);
constore=zeros(T,10);
capirf=zeros(T2,10);
conirf=zeros(T2,10);

%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;

for i = 1:10
    
nu = (i-1)*0.05;
save wouterrocks nu
%this saves nu into wouterrocks (which is then read by by Dynare)

dynare modelcriskaversion noclearall
% the command noclearall means that variables defined here (like T) are not
% deleted out of memory
% the file modelcriskaversion.mod makes it possible to read different
% values for riskaversion 

load dynarerocks
% this loads the matrix "decision" with coefficients into memory

% to simulate series we use random shocks
% set is reset in each iteration to ensure we get the same draw each time
%initialize stuff
cap    = zeros(T,1);
con    = zeros(T,1);
prod   = zeros(T,1);	
randn('state',20080701);
shocks = sig*randn(T,1); 
prod(1)= 0;


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(:,i)=exp(cap);
constore(:,i)=exp(con);

if isimul == 1
%plot(capstore(:,1:i))
%pause
plot(constore(2:end,1:i))
pause
plot(constore(91:110,1:i))
pause
end

% to calculate IRFs for linear policy function we set shocks (except for t=2) equal to zero
% set is reset in each iteration to ensure we get the same draw each time
%initialize stuff
cap    = zeros(T2,1);
con    = zeros(T2,1);
prod   = zeros(T2,1);	
randn('state',20080701);
shocks = zeros(T,1);
shocks(2) = sig; % a one standard deviation shock
prod(1)= 0;

capstst = decision(1,iiicap); 	
cap(1)  = capstst; % initial value is steady state level
constst = exp(alpha*capstst)-delta*exp(capstst); %consumption level steady state
constst = log(constst);
con(1)  = constst; % this is only needed for IRFs

for t = 2:T2	
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	


% without the shock in period 2 the values would have remained equal 
% to the steady state
% IRF is the percentage difference between series with and without shock
% since series are in logs the difference is the percentage difference

capirf(:,i)=cap-capstst;
conirf(:,i)=con-constst;

if iIRF ==1
plot(conirf(1:20,1:i))
pause
plot(capirf(1:20,1:i))
pause
end

end
