% IRF.m
%
% calculates IRFs for 2nd-order program
% there are different ways to calculate IRFs for non-linear models.
% this program takes as given the history of an economy from t=1 to t=T,
% and then asks the question how different things would have been if the
% innovation at t=Tirf would have been one standard deviation higher. The
% IRF is the difference between these two paths
%
% an alternative would to take some initial conditions (e.g. those at
% t=Tirf) and then generate many benchmark histories and for each of these
% benchmark histories resimulate the economy but with a one standard
% deviation higher innovation in t=Tirf. The IRF is then the AVERAGE of
% each of the differences
%
% One can think of the second as IRFs that are unconditional on future 
% draws, whereas the first is an IRF conditional not only on initial values 
% but also on the particular realization of future draws.

T=1000;

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

%
% first do a standard simulation
%

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	


% calculate IRFs when the initial conditions are like they are at t=Tirf
% Tirf must be less than T - 200
% shocks considered are in s (times standard deviation)

Tirf = 501;
s = [1; 5; -1; -5];


%initialize stuff
capalt    = zeros(T,1);
conalt    = zeros(T,1);
prodalt   = zeros(T,1);	
capirf    = zeros(200,4);
conirf    = zeros(200,4);
prodirf   = zeros(200,4);
capirfs   = zeros(200,4);
conirfs   = zeros(200,4);

for i = 1:4
    
% calculate what shocks would have been like if there would have been
% an ADDITIONAL shock at t = Tirf

shocksalt=shocks;
shocksalt(Tirf)=shocks(Tirf)+sig*s(i);

% recalculate what the series would have been
% you could start at t=2 (but we are only interested in what happens after
% Tirf)

capalt(Tirf-1)  = cap(Tirf-1); % same initial values
prodalt(Tirf-1) = prod(Tirf-1);

for t = Tirf:T	
prodalt(t)=rho*prodalt(t-1)+shocksalt(t);
capalt(t)=decision(1,iiicap) ...
      +decision(3,iiicap)*(capalt(t-1)-capstst) ...
      +decision(4,iiicap)*prodalt(t-1) ...
      +decision(5,iiicap)*shocksalt(t) ...
      +decision(6,iiicap)*(capalt(t-1)-capstst)^2 ...
      +decision(7,iiicap)*(capalt(t-1)-capstst)*prodalt(t-1) ...
      +decision(8,iiicap)*prodalt(t-1)^2 ...
      +decision(9,iiicap)*shocksalt(t)^2 ...
      +decision(10,iiicap)*(capalt(t-1)-capstst)*shocksalt(t) ...
      +decision(11,iiicap)*prodalt(t-1)*shocksalt(t);
conalt(t)=decision( 1,iiicons) ...
      +decision( 3,iiicons)*(capalt(t-1)-capstst) ...
      +decision( 4,iiicons)*prodalt(t-1) ...
      +decision( 5,iiicons)*shocksalt(t) ...
      +decision( 6,iiicons)*(capalt(t-1)-capstst)^2 ...
      +decision( 7,iiicons)*(capalt(t-1)-capstst)*prodalt(t-1) ...
      +decision( 8,iiicons)*prodalt(t-1)^2 ...
      +decision( 9,iiicons)*shocksalt(t)^2 ...
      +decision(10,iiicons)*(capalt(t-1)-capstst)*shocksalt(t) ...
      +decision(11,iiicons)*prodalt(t-1)*shocksalt(t);
end	

% define IRFs (IRF = perct. difference between what happens without and
% with the additional shock)
% since these variables are logs you get percentage difference simply by
% substraction

capirf(:,i)  =capalt(Tirf:Tirf+199) - cap(Tirf:Tirf+199);
conirf(:,i)  =conalt(Tirf:Tirf+199) - con(Tirf:Tirf+199);
prodirf(:,i) =prodalt(Tirf:Tirf+199)- prod(Tirf:Tirf+199);

% scale the IRFs to make them comparable

capirfs(:,i) = capirf(:,i)/s(i);
conirfs(:,i) = conirf(:,i)/s(i);

end


plot(capirf(1:150,:))
pause
plot(conirf(1:150,:))
pause
plot(capirfs(1:200,:))
pause
plot(conirfs(1:200,:))

% ! there are not much nonlinearities in this model so the scaled IRFs are
% virtually the same

