% Solving New Keynesian Models

clear all

% Parameter values
beta    = 0.99;
ksi     = 0.25;
gamma   = 2;
nu      = 2;
phi_y   = 0.5;
phi_pi  = 1.5;
rho_a   = 0.95;
rho_eta = 0.95;
rho_ki  = 0.95;

%%%%% The sticky prices model %%%%%
g0  = zeros(9,9);
g1  = zeros(9,9);
c   = zeros(9,1);
psi = zeros(9,3);
pi  = zeros(9,2);

% Aggregate supply
g0(1,3) = nu + gamma;
g0(1,5) = -1;
g0(1,7) = -(1+nu);
g0(1,8) = -1;

% IS curve
g0(2,2) = 1/gamma;
g0(2,3) = -1;
g0(2,4) = 1;
g0(2,6) = -1/gamma;
g0(2,8) = (1-rho_eta)/gamma;

% Interest policy rule
g0(3,1) = phi_pi;
g0(3,3) = phi_y;
g0(3,6) = -1;
g0(3,9) = 1;

% Phillips curve
g0(4,1) = -1;
g0(4,2) = beta;
g0(4,5) = ksi*(1-beta*(1-ksi))/(1-ksi);

% Productivity
g0(5,7)  = 1;
g1(5,7)  = rho_a;
psi(5,1) = 1;

% Demand shocks
g0(6,8)  = 1;
g1(6,8)  = rho_eta;
psi(6,2) = 1;

% Monetary policy shocks
g0(7,9)  = 1;
g1(7,9)  = rho_ki;
psi(7,3) = 1;

% Expectational errors
g0(8,1) = 1;    g1(8,2) = 1;    pi(8,1) = 1;    % Inflation
g0(9,3) = 1;    g1(9,4) = 1;    pi(9,2) = 1;    % Output


% Solving the model
[G1,C,impact,fmat,fwt,ywt,gev,eu] = gensys(g0,g1,c,psi,pi);
eu

% Impulse response functions
% Monetary shock
vars_sp(:,1) = impact*[0; 0; 1];
for i=2:101
    vars_sp(:,i) = G1*vars_sp(:,i-1);
end



%%%%% The sticky information model %%%%%
global N m n A A0 A1tilde Btilde BETA OMEGA THETA PSI
n  = 8;    % Number of endogenous variables
m  = 3;    % Number of exogenous variables
N  = 40;   % Number of lagged expectations
if N==0 & ksi<1
    display('Inconsistency in the parameters of the sticky information model: ksi<1 while N=0.');
    return
end
A0 = zeros(n,n);
A1 = zeros(n,n);
A2 = zeros(n,n);
if N>0
    GAMMA = zeros(n,n,N);
end
BETA = zeros(n,m);

% Aggregate supply
A1(1,2) = nu + gamma;
A1(1,3) = -1;
A1(1,5) = -(1+nu);
A1(1,6) = -1;

% IS curve
A0(2,1) = 1/gamma;
A0(2,2) = 1;
A1(2,2) = -1;
A1(2,4) = -1/gamma;
A1(2,6) = (1-rho_eta)/gamma;

% Interest policy rule
A1(3,1) = phi_pi;
A1(3,2) = phi_y;
A1(3,4) = -1;
A1(3,7) = 1;

% Phillips curve
A1(4,1) = -1;
A1(4,3) = ksi/(1-ksi);
if N>0
    for i=1:N
        GAMMA(4,1,i) = ksi*((1-ksi)^(i-1));
        GAMMA(4,8,i) = ksi*((1-ksi)^(i-1));
    end
end

% Productivity
A1(5,5)   = -1;
A2(5,5)   = rho_a;
BETA(5,1) = 1;

% Demand shocks
A1(6,6)   = -1;
A2(6,6)   = rho_eta;
BETA(6,2) = 1;

% Monetary policy shocks
A1(7,7)   = -1;
A2(7,7)   = rho_ki;
BETA(7,3) = 1;

% Real marginal cost - first difference
A1(8,3) = 1;
A1(8,8) = -1;
A2(8,3) = -1;

% Sims' matrices
if N>0
    A1tilde = A1 + sum(GAMMA,3);
else
    A1tilde = A1;
end
g0  = [A0         A1tilde;
       zeros(n,n) eye(n)];
g1  = [zeros(n,n) -A2;
       eye(n)     zeros(n,n)];
c   = zeros(2*n,1);
psi = zeros(2*n,N*m);
pi  = [zeros(n,n); eye(n)];


% Solving for A
[G1,C,impact,fmat,fwt,ywt,gev,eu] = gensys(g0,g1,c,psi,pi);
eu

G11 = G1(1:n,1:n);
G12 = G1(1:n,n+1:2*n);
G21 = G1(n+1:2*n,1:n);
G22 = G1(n+1:2*n,n+1:2*n);
A = G22;

% Check solution for A
tol = 1e-12;
Apol = A0*A*A + A1tilde*A + A2;
max_err = max(max(abs(Apol)));
if max_err > tol
    fprintf('Invalid solution for A. pol(A) does not equal zero. Max abs error: %6.4e',max_err);
    fprintf(' \n')
    return
elseif max(abs(eig(A))) >= 1
    fprintf('Invalid solution for A. pol(A) equals zero, but solution is unstable.');
    fprintf('Max eig: %6.4e',max(abs(eig(A))));
    fprintf('Max abs error: %6.4e',max_err);
    fprintf(' \n')
    return
else
    fprintf('Valid solution for A!  pol(A) equals zero. Max abs error: %6.4e',max_err);
    fprintf(' \n')
end


% Solving for B
% if N>0 then Bflex is used as an initial guess
Bflex = -inv(A0*A+A1tilde)*BETA;

% Law of motion for innovations
if N==0
    THETA = zeros(m,m);
else
    THETA = zeros(m*N,m*N);
    THETA(m+1:m*N,1:m*(N-1))=eye(m*(N-1));
end

if N==0
    B = Bflex;
else
    Btilde = [BETA zeros(n,m*(N-1))];
    for i=1:N
        OMEGA(:,(i-1)*n+1:i*n) = sum(GAMMA(:,:,i:N),3);
    end
    B_guess = [Bflex zeros(n,m*(N-1))];
    tic
    B = solveB_le(B_guess);
    toc
end

% Check solution for B
if N==0
    Beq = (A0*A+A1tilde)*B + A0*B*THETA + BETA;
else
    Beq = (A0*A+A1tilde)*B + A0*B*THETA + Btilde - OMEGA*PSI;
end
max_err = max(max(abs(Beq)));
if max_err > tol
    fprintf('Invalid solution for B. Beq does not equal zero. Max abs error: %6.4e',max_err);
    fprintf(' \n')
    return
else
    fprintf('Valid solution for B!  Beq equals zero.    Max abs error: %6.4e',max_err);
    fprintf(' \n')
end


% Impulse response functions
% Monetary shock
if N==0, NN=1; else, NN=N; end
shocks = [kron(eye(NN),[0; 0; 1]) zeros(NN*m,101)];
vars_si(:,1) = B*shocks(:,1);
for i=2:101
    vars_si(:,i) = A*vars_si(:,i-1) + B*shocks(:,i);
end

% Inflation
figure(1)
line('Xdata',[0:60],'Ydata',vars_sp(1,1:61), ...
     'Color','b','LineWidth',1.5,'Marker','d','MarkerFaceColor','w','MarkerSize',6)
line('Xdata',[0:60],'Ydata',vars_si(1,1:61), ...
     'Color','r','LineWidth',1.5,'Marker','o','MarkerFaceColor','w','MarkerSize',6)
axis([0 60 -1.6 0])
box off
set(gca,'Ygrid','on','XTick',[0:4:60],'YTick',[-1.6:0.2:0],'FontSize',20)
xlabel('Quarters','FontSize',16);
ylabel('Percentage Points','FontSize',16);
legend('Sticky Prices', ...
       'Sticky Information','Location','SouthEast')