%%--------------------------------------------------------------------------
%%  nk_policy_3e.m Basic New Keynesian Model
%%  Used to construct figures 8.2 and 8.3 of third edition of
%%  Monetary Theory and Policy
%%
%%  Based on programs by Paul Soderlind (nepso@hhs.se), 20 June 2000.
%%----------------------------------------------------------------------------*/
clear all
%%-----------------------------------------------------------------------*/

% choice of policy solution and shock
pc = 3; % commitment = 1, discretion = 2, both = 3
% choice of shock for impulse response plots
sk = 2; % natural rate shock = 1, cost shock = 2;

% Parameter values
lambda_x = .25;  %weight on employment variance in simple rules relative to inflation
beta = 0.99;
sig = 1/.16;       % coefficient of relative risk aversion
kappa = 0.05;      % coeff. on x in inflation adj. 
sdeps = 1;         % s.d. of cost shock innovation
sdrn = 1;          % s.d. of natural real rate innovation
rhorn = 0;         % serial correlation of natural real rate
rhoe =  0;         % serial correlation of cost shock

% Weights in loss function
lambda_x = .25;
lambda_i = 0;

disp('New Keynesian model: Chapter 8, Monetary Theory and Policy, 3rd ed.')
disp('weight on output and intereest rate variables in social loss function')
[lambda_x lambda_i]

%%
% Setting up model
n1 = 2; % rn e -- the two shocks (the natural real rate and the cost shock) 
n2 = 2; % pi, x -- inflation and the output gap

% housekeeping
index_pi = 3;
index_x = 4;
n = n1 + n2;
AA = zeros(n,n);
AB = zeros(n,n);
BB = zeros(n,1);
%%
%%
%% equations for rn, e, pi, x
% for rn
AA(1,1) = 1;
AB(1,1) = rhorn;

% for e
AA(2,2) = 1;
AB(2,2) = rhoe;

% NK Phillips curve
AA(3,3) = beta;
AB(3,2) = -1;
AB(3,3) = 1;
AB(3,4) = -kappa;

% IS equation
AA(4,3) = sig;
AA(4,4) = 1;
AB(4,1) = -1;
AB(4,4) = 1;

BB(4,1) = sig;
%%
invAA = inv(AA);
A = invAA*AB;
B = invAA*BB;

%%Loss function
%%  The period loss function is
%%
%%  x'Qx + 2x'Uu + u'Ru.
%%
Q  = zeros(size(A));
U  = zeros(size(A,1),1); 
R  = 0; 
Q(3,3) = 1;
Q(4,4) = lambda_x;

SigmaXX = zeros(n1,n1);
SigmaXX(1,1) = sdrn^2;
SigmaXX(2,2) = sdeps^2;

Tbig   = 16;                       %%periods to simulate in VAR*/
tid = linspace(1,Tbig,Tbig)';
tz = tid - tid;
   
%%-----------------------------------------------------------------------*/

pcc = pc;
if pc == 3;
    pcc = 3;
    pc  = 1;
end

if pc == 1;
    % chooses shock for impulse responses
    if sk == 1;
        Shock0 = [sdrn; zeros(n-1,1)];   
    else;
        if sk == 2;
            Shock0 = [0; sdeps; zeros(n-2,1)]; 
        end
    end
      
%% COMMITMENT*/
[ M_Commit,C_Commit ] = ComItAlg( A,B,Q,R,U,beta,n1,n2,1.0 );

k      = VarImpResp( M_Commit,Shock0,Tbig ); %%VAR of x1(t),p2(t)*/

lambda = k*C_Commit';               
uu     = lambda(:,n2+1:n2+1);       %%control -- i(t)*/
x      = [ lambda(1:Tbig,1:n2) ];   %%x(t) = [x1(t),x2(t)]*/

y_Commit = [ x(:,1) x(:,2)]; % pi, x

figure(1);
plot( tid,y_Commit,tid,tz,'k' );
axis tight
    if sk == 1;
        title('Commitment: impulse response to a one std natural rate shock');
    else;
        title('Commitment: impulse response to a one std cost shock');
    end

legend('inflation','output gap');
  
end ;    % end of pc = 1

if pcc == 3;
    pc = 2;
end


%----------------------------------------------
if pc == 2;
       
      if sk == 1;
          Shock0 = [sdrn; zeros(n1-1,1)]; 
      else;
          if sk == 2;
              Shock0 = [0; sdeps; zeros(n1-2,1)];  
          end
      end
   
%% ---------------------------------------------      
%% DISCRETION */

[M_Disc,C_Disc,VS,F ] = DiscAlg( A,B,Q,R,U,beta,n1,n2,eye(n1),zeros(n2,n1),...
                                 [1E-1;1E-5],0,1,0,1,1,100 );

x1 = zeros(Tbig,n1);
x2 = zeros(Tbig,n2);
x  = zeros(Tbig,n);
x1 = VarImpResp( M_Disc,Shock0,Tbig );  %%VAR of x1(t)*/
x2 = x1*C_Disc';                        %%x2(t)*/
x  = [x1 x2];                           %%x(t) = [x1(t),x2(t)]*/
uu = -x1*F';                            %%control -- i(t)*/
y_Disc = [ x2(1:Tbig,1:2) ]; % [pi x]

tid = linspace(1,Tbig,Tbig)';

figure(2);
plot( tid,y_Disc,tid,tz,'k' );
axis tight
if sk == 1;
    title('Discretion: impulse response to a one std natural rate shock');
else;
    title('Discretion: impulse response to a one std cost shock');
end
legend('inflation','output gap');

end; % end of pc = 2
   
%%-----------------------------------------------------------------------*/
figure (3)
plot(tid,tz,'k',tid,y_Commit(:,2),tid,y_Disc(:,2),':')
axis([0 16 -.5 .5])
title('Fig 8.2: Output gap: commitment (solid) and discretion (dotted)');
figure (4)
plot(tid,tz,'k',tid,y_Commit(:,1),tid,y_Disc(:,1),':')
axis([0 16 -.5 1.5])
title('Fig 8.3: Inflation: commitment (solid) and discretion (dotted)');

