% autor programu: Wouter den Haan (UvA)
% laskave zapujceni a interpretace: Pavel Herber (ESF)
% modifikace: Miroslav Hlousek (ESF)
% pomoc pri odladeni modelu: Daniel Nemec (ESF)

% !!!! running this program requires the updated version of disp_dr.m
% (which writes the policy functions to the matrix decision, exactly
% the way you see it on the screen). Alternatively, you can upload them
% from the files created by Dynare, but I find this easier especially for
% second-order


% Vse najdete na adrese 
% http://www1.feb.uva.nl/mint/wdenhaan/numerical/dynareprograms.htm
% 
% nutne nejdriv prekopirovat  disp_dr.m v dynaru, pripadne mit v adresari,
% ve kterem pracujete
% - dynare na zacatku nesmaze pamet
% - dynare ulozi do matice policy funkce


T = 100; % number of periods for IRFs

% matice pro vsechny irf
capirf=zeros(T,30);
conirf=zeros(T,30);
outirf=zeros(T,30);
tfpirf=zeros(T,30);
hourirf=zeros(T,30);
rateirf=zeros(T,30);
invirf=zeros(T,30);


% indicate which column of decision corresponds with which variable
% note that for each application you have to check that these column
% indications are correct

iout = 1
icon =2;
icap =3;
itfp = 4;
ihour = 5;
irate = 6;
iinv = 7;

alpha = .36
beta  = .99;
delta = .025;
psi = 1.72;
%rho = .95; 
sigma = .0035;

for i =  1:50
    
    rho = (i-1)*0.02;

save parameterfile alpha beta delta psi rho sigma
%this saves parameter values into parameterfile.mat
%which is then read by modelcloglinear2

dynare rbc_loglin_loop 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 of the policy functions
% 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

out    = zeros(T,1);
cap    = zeros(T,1);
con    = zeros(T,1);
tfp   = zeros(T,1);	
hour = zeros(T,1);
rate =  zeros(T,1); 
inv  =  zeros(T,1); 

shocks = zeros(T,1);
shocks(2) = sigma;    % one standard deviation shock
tfp(1)= 0;


for t = 2:T	
 tfp(t) = rho*tfp(t-1)+shocks(t);

 out(t) =decision(1,iout) + decision(2,iout)*cap(t-1) ...	
           +decision(3,iout)*tfp(t-1) + decision(4,iout)*shocks(t);
      
 cap(t) =decision(1,icap)+decision(2,icap)*cap(t-1) ...	
           +decision(3,icap)*tfp(t-1) + decision(4,icap)*shocks(t);
       
 con(t) =decision(1,icon)+decision(2,icon)*cap(t-1) ...	
           +decision(3,icon)*tfp(t-1) + decision(4,icon)*shocks(t);
 
 hour(t) =decision(1,ihour)+decision(2,ihour)*cap(t-1) ...	
           +decision(3,ihour)*tfp(t-1) + decision(4,ihour)*shocks(t);
       
 rate(t) =decision(1,irate)+decision(2,irate)*cap(t-1) ...	
           +decision(3,irate)*tfp(t-1) + decision(4,irate)*shocks(t);

 inv(t) =decision(1,iinv)+decision(2,iinv)*cap(t-1) ...	
           +decision(3,iinv)*tfp(t-1) + decision(4,iinv)*shocks(t);

end	

 outirf(:,i) = out;
 capirf(:,i) = cap;
 conirf(:,i) = con;
 hourirf(:,i) = hour;
 rateirf(:,i) = rate;
 invirf(:,i) = inv;
 
       
        plot(capirf(1:40,1:i-1))
        hold on
        plot(capirf(1:40,i),'linewidth',2.2,'color',[0 0 0])
        hold off
        title('irf for capital (different value of \rho )')
        
%         plot(outirf(1:40,1:i-1))
%         hold on
%         plot(outirf(1:40,i),'linewidth',2.2,'color',[0 0 0])
%         hold off
%         title('irf for output (different value of \rho )')

%         plot(hourirf(1:T,1:i-1))
%         hold on
%         plot(hourirf(1:T,i),'linewidth',2.2,'color',[0 0 0])
%         hold off
%         title('irf for hours worked (different value of \rho )')
%         
%         plot(invirf(1:T,1:i-1))
%         hold on
%         plot(invirf(1:T,i),'linewidth',2.2,'color',[0 0 0])
%         hold off
%         title('irf for investment (different value of \rho )')
end    



