% 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 = 1000; % number of periods for IRFs

% napiste, ktery sloupec matice decision nalezi prislusne promenne
% Pozor, pro kazdou aplikaci musite zkontrolovat, ze tyto sloupce jsou
% spravne

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

dynare indivisible_rbc_plus noclearall
%dynare basic_rbc_plus 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 calculate IRFs for linear policy function we set shocks (except for t=2) equal to zero
% initialize stuff

%% vypocet impulsnich odezev z (linearniho) rozhodovaciho pravidla
% preddefinovani promennych 

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

randn('state',20101111);  % nastaveni generatoru nahodnych cisel
shocks = sigma*randn(T,1); % sok ve vsech obdobich
%shocks(2) = sigma;    % puvodne sok jen v jednom obdobi


% cyklus pro sok
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	

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


%% vykresleni
tt = (1:T)';

figure
plot(tt,cap(tt),tt,con(tt),tt,out(tt),tt,hour(tt),tt,tfp(tt),tt,inv(tt));
legend('capital','consumption','output','hours','tfp','investment')
title('impulse response function')
ylabel('deviation from s.s.')


% zkracene irf, sok v obdobi 0, vyjadrene v procentech (nikoliv v
% desetinnem vyjadreni)

% ttt = (2:T)';
% time = (0:T-2)';
% figure
% plot(time,cap(ttt)*100,time,con(ttt)*100,time,out(ttt)*100,time,tfp(ttt)*100)
% legend('capital','consumption','output','tfp')
% title('impulse response function')
% ylabel('deviation from s.s. in %')

stdev.y = std(out);
stdev.h = std(hour);
stdev.yh = stdev.h/stdev.y;

