% Author: Olivier Loisel (CREST). Date: June 2020.
% This Matlab file checks Figure 2 and the Supplementary Appendix S.2 in the paper "The Implementation of Stabilization Policy" (O. Loisel, June 2020).
% It calls two Dynare files: ISP_Fig2_check_ramsey.mod and ISP_Fig2_check_rule.mod (which are consistent at least with version 4.3.0 of Dynare).
% The names of the objects that appear also in the paper (scalars, vectors, matrices) correspond to their notations in the paper.
% A large part of this program (specified in the code) is copied and pasted from the file ISP_Fig2.m.
% The inputs are cc (referring to the sub-figure of Figure 2 considered), rho_eta and rho_u (referring to the specific point considered in the plane of this sub-figure).
% This file determines whether Blanchard and Kahn's (1980) root-counting and no-decoupling conditions are met under the rule (41).
% When there are met, it checks that the corresponding impulse-reponse functions (IRFs) coincide with those obtained using the ramsey_policy command of Dynare.
% Running this file for the first time may generate an error message of type "Undefined function or variable", but things work well from the second time.

cc = 0; % integer value to be chosen equal to 0 [for Gali's (2015, Chapter 3) calibration] or equal to 1 [for Woodford's (2003, Chapter 4) calibration]
rho_eta = 0.5; % auto-regressive parameter for disturbance eta (real-number value to be chosen in (0,1))
rho_u = 0.5; % auto-regressive parameter for disturbance u (real-number value to be chosen in (0,1))

beta = 0.99; % discount factor
theta = 0.67; % probability of not being allowed to change one's price
if cc==0 % Gali's (2015, Chapter 3) calibration
    alpha = 0.67; % labor elasticity of output
    epsilon = 6; % steady-state elasticity of substitution between differentiated goods
    sigma = 1; % inverse of intertemporal elasticity of substitution
    chi = 1; % inverse of Frisch labor supply elasticity
    kappa = (1-theta)*(1-beta*theta)/theta*(alpha*sigma+chi+1-alpha)/(alpha+(1-alpha)*epsilon); % slope of the Phillips curve (when there is a single kind of labor)
end
if cc==1 % Woodford's (2003, Chapter 4) calibration
    alpha = 0.75;
    epsilon = 7.88;
    sigma = 0.16;
    chi = 0.11;
    kappa = (1-theta)*(1-beta*theta)/theta*(alpha*sigma+chi+1-alpha)/(alpha+(chi+1-alpha)*epsilon); % slope of the Phillips curve (when there are differentiated kinds of labor)
end
lambda = kappa/epsilon; % relative weight on output-gap-stabilization objective (for quarterly inflation rates)
save('Param.mat','beta','theta','rho_eta','rho_u','sigma','kappa','lambda');

disp(' ');
disp('Step 1: compute the IRFs using the ramsey_policy command of Dynare');
dynare ISP_Fig2_check_ramsey;
IRF_ramsey_pi_epseta = oo_.irfs.pi_epseta ;
IRF_ramsey_pi_epsu = oo_.irfs.pi_epsu ;
IRF_ramsey_y_epseta = oo_.irfs.y_epseta ;
IRF_ramsey_y_epsu = oo_.irfs.y_epsu ;
IRF_ramsey_i_epseta = oo_.irfs.i_epseta ;
IRF_ramsey_i_epsu = oo_.irfs.i_epsu ;
save('ISP_Fig2_check_ramsey.mat','IRF_ramsey_pi_epseta','IRF_ramsey_pi_epsu','IRF_ramsey_y_epseta','IRF_ramsey_y_epsu','IRF_ramsey_i_epseta','IRF_ramsey_i_epsu');

disp(' ');
disp('Step 2: determine whether Blanchard and Kahn''s (1980) conditions are met under the rule (41) and, if they are, compute the corresponding IRFs');
load('Param.mat');

% Start of the part of the code that is copied and pasted from the file ISP_Fig2.m
mu = ((beta*lambda+kappa^2+lambda)-sqrt((beta*lambda+kappa^2+lambda)^2-4*beta*lambda^2))/(2*beta*lambda);
TZSW11 = [0];
TZSW12 = [(-beta*rho_u^2*mu^2/(1-beta*rho_u*mu)) (rho_u*mu/(1-beta*rho_u*mu)-rho_u-mu) 1];
TZSW21 = [(rho_u*mu) (-(rho_u+mu)) 1];
TZSW22 = [(-kappa*rho_u*mu/(lambda*(1-beta*rho_u*mu))) 0];
TiSW11 = [(sigma*rho_eta*rho_u*mu) (-sigma*rho_eta*(rho_u+mu)) (sigma*rho_eta)];
TiSW12 = [((lambda-kappa*sigma)*rho_eta*rho_u^2*mu^2/(lambda*(1-beta*rho_u*mu))) (-(lambda-kappa*sigma)*rho_u*mu*(rho_u*mu-rho_eta*(1-rho_u-mu))/(lambda*(1-beta*rho_u*mu))) (-(lambda-kappa*sigma)*rho_u*mu*(1-rho_u-mu)/(lambda*(1-beta*rho_u*mu)))];
[G,RG] = deconv(conv([-rho_eta 1],-conv(TZSW12,TZSW21)),[rho_u*mu (-rho_u-mu) 1]);
[F_pi,RF_pi] = deconv([0 conv(TiSW11,TZSW22)]-conv(TiSW12,TZSW21),[rho_u*mu (-rho_u-mu) 1]);
[F_y,RF_y] = deconv(-conv(TiSW11,TZSW12),[rho_u*mu (-rho_u-mu) 1]);
% End of the part of the code that is copied and pasted from the file ISP_Fig2.m

G_ext = G;
F_pi_ext = F_pi;
F_y_ext = F_y;

save('ISP_Fig2_check.mat','F_pi_ext','F_y_ext','G_ext');
dynare ISP_Fig2_check_rule;
disp(' ');

if (isfield(oo_,'irfs'))
    disp('Blanchard and Kahn''s (1980) root-counting and no-decoupling conditions are met under the rule (41)');
    disp(' ');
    disp('Step 3: compute the distance between the two sets of IRFs');
    disp(' ');
    IRF_rule_pi_epseta = oo_.irfs.pi_epseta ;
    IRF_rule_pi_epsu = oo_.irfs.pi_epsu ;
    IRF_rule_y_epseta = oo_.irfs.y_epseta ;
    IRF_rule_y_epsu = oo_.irfs.y_epsu ;
    IRF_rule_i_epseta = oo_.irfs.i_epseta ;
    IRF_rule_i_epsu = oo_.irfs.i_epsu ;
    load('ISP_Fig2_check_ramsey.mat');
    Check = max([norm(IRF_ramsey_pi_epseta-IRF_rule_pi_epseta); norm(IRF_ramsey_pi_epsu-IRF_rule_pi_epsu); norm(IRF_ramsey_y_epseta-IRF_rule_y_epseta); ...
    norm(IRF_ramsey_y_epsu-IRF_rule_y_epsu); norm(IRF_ramsey_i_epseta-IRF_rule_i_epseta); norm(IRF_ramsey_i_epsu-IRF_rule_i_epsu)]);
    disp(['distance = ', num2str(Check)]);
elseif (error == 3)
    disp('Blanchard and Kahn''s (1980) root-counting condition is not met: non-robustness');
elseif (error == 4)
    disp('Blanchard and Kahn''s (1980) root-counting condition is not met: multiplicity');
elseif (error == 5)
    disp('Blanchard and Kahn''s (1980) root-counting condition is met, but their no-decoupling condition is not: multiplicity and non-robustness');
end
disp(' ');