%Hopefully 2nd order IRFs for SimpleGrowth.mod
clear
randn('state',1);
dynare SimpleGrowth.mod

T=41; %Length of IRFs (+1)
sigma_e=0.07; %Shock size
eee=zeros(T,1); %Reserving space
eee(1,1)=sigma_e; %Shock is the first element

% The shock series will be the same for each solution. 
randn('state',1);
shocks = sigma_e*randn(T,1);
shocks1 = [shocks(1,1)+sigma_e;shocks(2:end,1)];

load dynarerocks
X=decision;
%Reserving space:
y = zeros(T,1);
c = zeros(T,1);
k = zeros(T,1);
z = zeros(T,1);  
%For Dynare steady states use:
            y(1,1)    = X(1,1)-X(2,1);
            c(1,1)    = X(1,2)-X(2,2);
            k(1,1)    = X(1,3)-X(2,3);
            z(1,1)    = X(1,4)-X(2,4);
            y1(1,1)    = X(1,1)-X(2,1);
            c1(1,1)    = X(1,2)-X(2,2);
            k1(1,1)    = X(1,3)-X(2,3);
            z1(1,1)    = X(1,4)-X(2,4);

            
%Recursion; NOTE: The first element of series contains steady states!
        for i=2:T
            %S=[1,0,k(i-1)-k(1,1),z(i-1)-z(1,1),e(i-1),(k(i-1)-k(1,1))*(k(i-1)-k(1,1)),(z(i-1)-z(1,1))*(k(i-1)-k(1,1)),(z(i-1)-z(1,1))*(z(i-1)-z(1,1)),e(i-1)*e(i-1),(k(i-1)-k(1,1))*e(i-1),(z(i-1)-z(1,1))*e(i-1)];
            S =[1,0,k(i-1)-k(1,1),z(i-1)-z(1,1),shocks(i-1),(k(i-1)-k(1,1))*(k(i-1)-k(1,1)),(z(i-1)-z(1,1))*(k(i-1)-k(1,1)),(z(i-1)-z(1,1))*(z(i-1)-z(1,1)),shocks(i-1)*shocks(i-1),(k(i-1)-k(1,1))*shocks(i-1),(z(i-1)-z(1,1))*shocks(i-1)];
            S1=[1,0,k1(i-1)-k1(1,1),z1(i-1)-z1(1,1),shocks1(i-1),(k1(i-1)-k1(1,1))*(k1(i-1)-k1(1,1)),(z1(i-1)-z1(1,1))*(k1(i-1)-k1(1,1)),(z1(i-1)-z1(1,1))*(z1(i-1)-z1(1,1)),shocks1(i-1)*shocks1(i-1),(k1(i-1)-k1(1,1))*shocks1(i-1),(z1(i-1)-z1(1,1))*shocks1(i-1)];

            y(i,1)     = S*X(:,1);
            c(i,1)     = S*X(:,2);
            k(i,1)     = S*X(:,3);
            z(i,1)     = S*X(:,4);
            y1(i,1)    = S1*X(:,1);
            c1(i,1)    = S1*X(:,2);
            k1(i,1)    = S1*X(:,3);
            z1(i,1)    = S1*X(:,4);
        end

    cIRF(:,1) = c1-c;
    yIRF(:,1) = y1-y;
    kIRF(:,1) = k1-k;
    zIRF(:,1) = z1-z;

% figure
% subplot(2,2,1),plot((yIRF(2:end,:))),title('y')
% subplot(2,2,2),plot((cIRF(2:end,:))),title('c')
% subplot(2,2,3),plot((kIRF(2:end,:))),title('k')
% subplot(2,2,4),plot((zIRF(2:end,:))),title('z')

figure
subplot(2,2,1),plot((y_e-yIRF(2:end,:))),title('Difference: Dynare and own IRF for y')
subplot(2,2,2),plot((c_e-cIRF(2:end,:))),title('Difference: Dynare and own IRF for c')
subplot(2,2,3),plot((k_e-kIRF(2:end,:))),title('Difference: Dynare and own IRF for k')
subplot(2,2,4),plot((z_e-zIRF(2:end,:))),title('Difference: Dynare and own IRF for z')