%% Program for solving a basic new Keynesian model with a simple Taylor
%% rule
% 
%
% pi(t) = beta E_t pi(t+1) + g y(t) + epi(t)
% y(t)  = E_t y(t+1) - br(R(t) - E_t pi(t+1)] + ey(t)
% R(t) = gpi pi(t) + gy y(t)) + er(t)
% epi(t) = rhopi pi(t-1 ) + upi(t)
% ey(t) = rhoy ey(t-1) + uy(t)
% er(t) = rhor er(t-1) + ur(t)
% 

%% Housekeeping 
% variables listed in following order:
% x1 = [ep(t+1) ey(t+1) ei(t+1) R(t)]
% x2 = [E_t y(t+1) E_t pi(t+1)]
% z = [x1 x2]';
clear all
n1 = 4;
n2 = 2;
n = n1+n2;
A0 = eye(n);
A1 = zeros(n);
B1 = zeros(n,1);
SigmaXX = zeros(n,n);
Tbig = 20;        % number of periods for impulse responses
gra = 'joint   '; % or 'separate'

%% Parameters
p.beta = 0.99;
p.g = 0.0081;
p.br = 0.0048;
p.gpi = 1.5;
p.gy = 0.5/4;

% Serial correlation parameters for disturbances
p.rhopi = 0;
p.rhoy = 0;
p.rhor = 0;

% Innovation variances (in percent)
cp = 0.5923;
cy = 0.4162;
ci = 0.9918;
p.vp = cp^2;
p.vy = cy^2;
p.vr = ci^2;

%% Matrices for model
% A0*E_z(t+1) = A1*z(t) + e(t+1)
% 
% E_z(t+1) = inv(A0)*A1*z(t) + inv(A0)*e(t+1)

% Exogenous shocks
A1(1,1) = p.rhopi;
A1(2,2) = p.rhoy;
A1(3,3) = p.rhor;

% Policy equation
A1(4,3) = 1;
A1(4,5) = p.gy;
A1(4,6) = p.gpi;

% Expectational IS relationship
A0(5,5) = 1;
A0(5,6) = p.br;
A0(5,4) = -p.br;
A1(5,2) = -1;
A1(5,5) = 1;

% Inflation equation
A0(6,6) = p.beta;
A1(6,1) = -1;
A1(6,5) = -p.g;
A1(6,6) = 1;

SigmaXX(1,1) = p.vp;
SigmaXX(2,2) = p.vy;
SigmaXX(3,3) = p.vr;

A = A0\A1;

%% Solving the model using basicsol.m
%  Output:    M      -   n1xn1 matrix, x1(t+1) = M*x1(t) + e(t+1)
%             C          n2xn1 matrix, x2(t) = C*x1(t)

x10 = zeros(n1,1);
[M,C] = basicsol(A,n1,n2,SigmaXX,x10,1.0 );

disp('M')
M
disp('C')
C

%% Constructing the impulse responses
% Output gap is variable 9, inflation is variable 10.
% Nominal interest rate is variable 5

x = 1:Tbig-1;
for i = 1:3;  % loops through each of the three shocks (cost, demand, policy)
shk = i; % 1 = inflation shock, 2 = demand shock, 3 = policy shock
shocks = zeros(n1,1);
shocks(shk,1) = (SigmaXX(shk,shk))^.5; % shock is one std. dev. of shock # shk
IMPx1 = VarImpResp(M,shocks,Tbig);
IMPx2 = IMPx1*C';
if shk == 1;
    IMP1 = [IMPx1 IMPx2];
    if gra == 'separate'
    figure (1)    
    plot(x,IMPx1(1:Tbig-1,5),x,IMPx1(1:Tbig-1,6),'.',x,IMPx1(2:Tbig,4),':')
    title('Impulse response of output gap and inflation to inflation shock')
    legend('output gap', 'inflation','interest rate')
    end
else
    if shk == 2;
        IMP2 = [IMPx1 IMPx2];
        if gra == 'separate'
        figure (2)
        plot(x,IMPx1(1:Tbig-1,5),x,IMPx1(1:Tbig-1,6),'.',x,IMPx1(2:Tbig,4),':')
        title('Impulse response of output gap and inflation to demand shock')
        legend('output gap', 'inflation','interest rate')
        end
    else
        IMP3 = [IMPx1 IMPx2];
        if gra == 'separate'
        figure (3)
        plot(x,IMPx1(1:Tbig-1,5),x,IMPx1(1:Tbig-1,6),'.',x,IMPx1(2:Tbig,4),':')
        title('Impulse response of output gap and inflation to policy shock')
        legend('output gap', 'inflation','interest rate')
        end
    end
end
end; % for i statement

if gra == 'joint   ';
figure (4)
x = 1:Tbig;
subplot(3,1,1)
plot(x,IMP1(1:Tbig,5),x,IMP1(1:Tbig,6),':');  
legend('output gap', 'inflation')
title('Responses to inflation shock')
subplot(3,1,2)
plot(x,IMP2(1:Tbig,5),x,IMP2(1:Tbig,6),':');  
legend('output gap', 'inflation')
title('Responses to demand shock')
subplot(3,1,3)
plot(x,IMP3(1:Tbig,5),x,IMP3(1:Tbig,6),':'); 
legend('output gap', 'inflation')
title('Responses to policy shock')
end






