close all
clear all
clc


dynare model_withebpp_3yvesting noclearall 
load dynarerocks


X=decision;


%--------------------------------------------------------------------------
%  Computing impulse response functions
%--------------------------------------------------------------------------
 
s1  = 0.08;           %Standard deviation of the shock
s2  = 0.01;
Ti = 24;              %Length of the series

%Shock series:
ez    = zeros(Ti,1);
ev    = zeros(Ti,1);

ev(2) = s1;
%ez(2) = s2;

% defining variables
% Reserving space:
k = zeros(Ti,1);
n = zeros(Ti,1);
w = zeros(Ti,1);
e = zeros(Ti,1);
cash = zeros(Ti,1);
theta=zeros(Ti,1);
c=zeros(Ti,1);
z = zeros(Ti,1);
varphi = zeros(Ti,1);

%Initial values:
k(1,1) = X(1,2);
n(1,1) = X(1,3);
w(1,1) = X(1,4);
theta(1,1) = X(1,8);
e(1,1) = X(1,9);
cash(1,1) = X(1,10);
z(1,1) = X(1,11);
varphi(1,1) = X(1,12);

%Recursion that computes simulated data using the shock series and the
%policy functions:

for t = 2:Ti
    k(t,1)     = X(1,2) + X(2,2)*(n(t-1)-n(1)) + X(3,2)*(k(t-1)-k(1)) +  X(4,2)*(w(t-1)-w(1)) + X(5,2)*(theta(t-1)-theta(1)) + X(6,2)*(z(t-1)-z(1)) + X(7,2)*(theta(t-2)-theta(1)) + X(8,2)*(w(t-2)-w(1)) + X(9,2)*(n(t-2)-n(1)) + X(10,2)*(theta(t-3)-theta(1)) + X(11,2)*(w(t-3)-w(1)) + X(12,2)*(n(t-3)-n(1)) + X(13,2)*(c(t-1)-c(1)) + X(14,2)*(varphi(t-1)-varphi(1)) +X(15,2)*ez(t) + X(16,2)*ev(t);
    n(t,1)     = X(1,3) + X(2,3)*(n(t-1)-n(1)) + X(3,3)*(k(t-1)-k(1)) +  X(4,3)*(w(t-1)-w(1)) + X(5,3)*(theta(t-1)-theta(1)) + X(6,3)*(z(t-1)-z(1)) + X(7,3)*(theta(t-2)-theta(1)) + X(8,3)*(w(t-2)-w(1)) + X(9,3)*(n(t-2)-n(1)) + X(10,3)*(theta(t-3)-theta(1)) + X(11,3)*(w(t-3)-w(1)) + X(12,3)*(n(t-3)-n(1)) + X(13,3)*(c(t-1)-c(1)) + X(14,3)*(varphi(t-1)-varphi(1)) +X(15,3)*ez(t) + X(16,3)*ev(t);
    w(t,1)     = X(1,4) + X(2,4)*(n(t-1)-n(1)) + X(3,4)*(k(t-1)-k(1)) +  X(4,4)*(w(t-1)-w(1)) + X(5,4)*(theta(t-1)-theta(1)) + X(6,4)*(z(t-1)-z(1)) + X(7,4)*(theta(t-2)-theta(1)) + X(8,4)*(w(t-2)-w(1)) + X(9,4)*(n(t-2)-n(1)) + X(10,4)*(theta(t-3)-theta(1)) + X(11,4)*(w(t-3)-w(1)) + X(12,4)*(n(t-3)-n(1)) + X(13,4)*(c(t-1)-c(1)) + X(14,4)*(varphi(t-1)-varphi(1)) +X(15,4)*ez(t) + X(16,4)*ev(t);
    theta(t,1) = X(1,8) + X(2,8)*(n(t-1)-n(1)) + X(3,8)*(k(t-1)-k(1)) +  X(4,8)*(w(t-1)-w(1)) + X(5,8)*(theta(t-1)-theta(1)) + X(6,8)*(z(t-1)-z(1)) + X(7,8)*(theta(t-2)-theta(1)) + X(8,8)*(w(t-2)-w(1)) + X(9,8)*(n(t-2)-n(1)) + X(10,8)*(theta(t-3)-theta(1)) + X(11,8)*(w(t-3)-w(1)) + X(12,8)*(n(t-3)-n(1)) + X(13,8)*(c(t-1)-c(1)) + X(14,8)*(varphi(t-1)-varphi(1)) +X(15,8)*ez(t) + X(16,8)*ev(t);
    c(t,1)     = X(1,7) + X(2,7)*(n(t-1)-n(1)) + X(3,7)*(k(t-1)-k(1)) +  X(4,7)*(w(t-1)-w(1)) + X(5,7)*(theta(t-1)-theta(1)) + X(6,7)*(z(t-1)-z(1)) + X(7,7)*(theta(t-2)-theta(1)) + X(8,7)*(w(t-2)-w(1)) + X(9,7)*(n(t-2)-n(1)) + X(10,7)*(theta(t-3)-theta(1)) + X(11,7)*(w(t-3)-w(1)) + X(12,7)*(n(t-3)-n(1)) + X(13,7)*(c(t-1)-c(1)) + X(14,7)*(varphi(t-1)-varphi(1)) +X(15,7)*ez(t) + X(16,7)*ev(t);
    e(t,1)     = X(1,9) + X(2,9)*(n(t-1)-n(1)) + X(3,9)*(k(t-1)-k(1)) +  X(4,9)*(w(t-1)-w(1)) + X(5,9)*(theta(t-1)-theta(1)) + X(6,9)*(z(t-1)-z(1)) + X(7,9)*(theta(t-2)-theta(1)) + X(8,9)*(w(t-2)-w(1)) + X(9,9)*(n(t-2)-n(1)) + X(10,9)*(theta(t-3)-theta(1)) + X(11,9)*(w(t-3)-w(1)) + X(12,9)*(n(t-3)-n(1)) + X(13,9)*(c(t-1)-c(1)) + X(14,9)*(varphi(t-1)-varphi(1)) +X(15,9)*ez(t) + X(16,9)*ev(t);
    cash(t,1)  = X(1,10)+ X(2,10)*(n(t-1)-n(1))+ X(3,10)*(k(t-1)-k(1))+  X(4,10)*(w(t-1)-w(1))+ X(5,10)*(theta(t-1)-theta(1))+ X(6,10)*(z(t-1)-z(1))+ X(7,10)*(theta(t-2)-theta(1))+ X(8,10)*(w(t-2)-w(1))+ X(9,10)*(n(t-2)-n(1))+ X(10,10)*(theta(t-3)-theta(1))+ X(11,10)*(w(t-3)-w(1))+ X(12,10)*(n(t-3)-n(1))+ X(13,10)*(c(t-1)-c(1))+ X(14,10)*(varphi(t-1)-varphi(1))+X(15,10)*ez(t)+ X(16,10)*ev(t);
    z(t,1)     = X(1,11)+ X(2,11)*(n(t-1)-n(1))+ X(3,11)*(k(t-1)-k(1))+  X(4,11)*(w(t-1)-w(1))+ X(5,11)*(theta(t-1)-theta(1))+ X(6,11)*(z(t-1)-z(1))+ X(7,11)*(theta(t-2)-theta(1))+ X(8,11)*(w(t-2)-w(1))+ X(9,11)*(n(t-2)-n(1))+ X(10,11)*(theta(t-3)-theta(1))+ X(11,11)*(w(t-3)-w(1))+ X(12,11)*(n(t-3)-n(1))+ X(13,11)*(c(t-1)-c(1))+ X(14,11)*(varphi(t-1)-varphi(1))+X(15,11)*ez(t)+ X(16,11)*ev(t);
    varphi(t,1)= X(1,12)+ X(2,12)*(n(t-1)-n(1))+ X(3,12)*(k(t-1)-k(1))+  X(4,12)*(w(t-1)-w(1))+ X(5,12)*(theta(t-1)-theta(1))+ X(6,12)*(z(t-1)-z(1))+ X(7,12)*(theta(t-2)-theta(1))+ X(8,12)*(w(t-2)-w(1))+ X(9,12)*(n(t-2)-n(1))+ X(10,12)*(theta(t-3)-theta(1))+ X(11,12)*(w(t-3)-w(1))+ X(12,12)*(n(t-3)-n(1))+ X(13,12)*(c(t-1)-c(1))+ X(14,12)*(varphi(t-1)-varphi(1))+X(15,12)*ez(t)+ X(16,12)*ev(t);
    
end

% Define impulse responses as %deviations from the steady state 

figure
title('XXX')
subplot(3,2,1)
plot(1:Ti-1,((n(2:end) - n(1))/n(1))*100)
title('Employment')
subplot(3,2,2)
plot(1:Ti-1,((k(2:end) - k(1))/k(1))*100)
title('Capital')
subplot(3,2,3)
plot(1:Ti-1,((cash(2:end) - cash(1))/cash(1))*100)
title('CW')
subplot(3,2,4)
plot(1:Ti-1,((w(2:end) - w(1))/w(1))*100)
title('Wage')
%annotation('textbox','string',{'IRFs'},'FontWeight','bold', ...
%    'FontSize',12,'FitBoxToText','on','Position', ...
%   [0.45 0.9903 0 0]);


