%%--------------------------------------------------------------------------
%%  prg_nk_roots.m Basic new Keynesian model, solves for eigenvalues
%
% Calculates min. egenvalue as function of policy coefficient on inflation
% Graphs result
% Demonstrates the Taylor Principle
% Can use rule based on either inflation or expected inflation.
% (see also prg_nk_roots2.m for rule that includes the output gap)
%%----------------------------------------------------------------------------*/
clear all
%%-----------------------------------------------------------------------*/
% Choose form of simple rule
% if rule == 0, i = delta*pi(t)
% if rule == 1, i = delta*Et pi(t+1)
rule = 1;

%% model parameters
beta = 0.99;
sig = 1;
w = 0.5;
kappa = (1-w)*(1-w*beta)/w;      % coeff. on x in inflation adj. 
sde = 1;                         % s.d. of cost shock innovation
sdu = 1;                         % s.d. of demand innovation
sdv = 1;                         % s.d. of policy shock

gamma = 0.5; % coefficient on interest rate in NKPC

N = 20; % number of values of delta to evaluate
M = zeros(20,2); % will be used to store eigenvalues
d = zeros(N,1);
%% Calculating eigenvalues
% delta is the coefficient on inflation in the policy rule
if rule == 0;
    for i = 1:N+1; %  will cycle through N values for delta
        delta = .1*(i-1); % delta runs from 0 to .1*N
        d(i) = delta;
        % read in matrices for model
        % Z = [pi x]'
        % AA*E(t)Z(t+1) = BB*Z(t) + U(t+1)
        AA = [1 1/sig; 0 beta]; %[beta 0; 1/sig 1];
        BB = [1 delta/sig; -kappa 1-delta*gamma]; % [1  -kappa; delta/sig 1];
        % E(t)Z(t+1) = A*Z(t) + V(t+1)
        A = inv(AA)*BB;
        M(i,:) = abs(eig(A))'; % each row of M gives the two eigenvalues of A for a given delta
    end
else
    for i = 1:N+1; %  will cycle through N values for delta
        delta = .5*(i-1); % delta runs from 0 to .1*N
        d(i) = delta;
        % read in matrices for model
        % Z = [pi x]'
        % AA*E(t)Z(t+1) = BB*Z(t) + U(t+1)
        AA = [1 1/sig; 0 beta]; %[beta 0; 1/sig 1];
        BB = [1 delta/sig; -kappa 1-delta*gamma]; % [1  -kappa; delta/sig 1];
        % E(t)Z(t+1) = A*Z(t) + V(t+1)
        A = inv(AA)*BB;
        M(i,:) = abs(eig(A))'; % each row of M gives the two eigenvalues of A for a given delta
    end
end
    
minEig = min(M');  % min(M') gives min value of each row of M, stores as column vector minEig
maxEig = max(M');

%% Graphing results
tid = 1:N+1; tid = 0.1*(tid' -1);  % creates series equal to values of delta to use as axis
ti0 = ones(N+1,1); % col. of ones
plot(d,minEig,d,ti0,d,maxEig);
title('Minimum eigenvalue as function of \delta')
ylabel('Minimum eigenvalue')
if rule == 0;
    xlabel('Coefficient on inflation: \delta');
else
    xlabel('Coefficient on expected future inflation: \delta');
end
