% autor programu: Wouter den Haan (UvA)
% laskave zapujceni a interpretace: Pavel Herber (ESF, KBC Bank, PPF Vietnam Finance)
% modifikace: Miroslav Hloušek (ESF)
% pomoc pri odladeni modelu: Daniel Němec (ESF)

% Vše najdete na adrese 
% http://www.wouterdenhaan.com/numerical/dynareprograms.htm
% nutne mít soubor disp_dr.m v pracovním adresáři: 
% - dynare na začátku nesmaže pamět
% - dynare uloží do matice policy function

set(0,'defaultlinelinewidth',1.5); % tloušťka čáry

% indexovani proměnných (který sloupec matice decision náleží příslušné proměnné
% Pro každou aplikace je nutné zkontrolovat, zda je pořadí odpovídá

iout = 1;       %výstup 
icon =2;        %spotřeba
icap =3;        %kapitál
ihour = 4;      %odpracované hodiny
irate = 5;      %nájemní cena kapitálu
itfp = 6;       %šok v technologii
iinv = 7;       %investice

dynare indivisible_rbc_plus.mod noclearall  % spustí mod file a nevymaže paměť

load dynarerocks    % nahraje matici  "decision" s koeficienty rozhodovacího pravidla


%% vypocet impulsnich odezev z (linearniho) rozhodovaciho pravidla

% předdefinovani promennych 
T = 100;     % definování počtu období 

out    = zeros(T,1);    % (y)
con    = zeros(T,1);    % (c)
cap    = zeros(T,1);    % (k)
hour  = zeros(T,1);     % (h)
rate =  zeros(T,1);     % (r)
tfp   = zeros(T,1);     % (z)
inv   = zeros(T,1);     % (i)


shocks = zeros(T,1);    % epsilon
%  shocks(2) = sigma;      % % pro irf, je šok pouze ve druhém období, velikost jedné směrodatné odchylky

randn('state',20151203);
shocks = randn(T,1)*sigma;

% Správně bychom měli odečitat s.s. hodnotu [viz promenna out(t)], ale ta nulová proto dále vynechano
% podobne s.s. u kazde promenne (na zacatku) je take 0

capss = decision(1,3);
tfpss = decision(1,6);

for t = 2:T	
 tfp(t) = rho*tfp(t-1)+shocks(t);   % nejprve cyklus pro sok

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

% bez šoku v období 2 by proměnné zůstaly na svých steady statových hodnotách
% IRF je v procentních odchylkách od steady statu


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

figure
plot(tt,cap(tt),tt,con(tt),tt,out(tt),tt,hour(tt),tt,rate(tt),tt, tfp(tt),tt,inv(tt));
legend('capital','consumption','output','hours','interest rate','tfp','investment')
title('Simulace')
ylabel('odchylka od 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('Impulzní odezva')
% ylabel('odchylka od  s.s. v %')

% 
sdi.y = std(out)*100;
sdi.h = std(hour)*100;
sdi.h_y = sdi.h/sdi.y;

disp('volatilita (std): y  h  h/y')
disp([sdi.y sdi.h sdi.h_y])
 
