% Author: Olivier Loisel (CREST). Date: June 2020.
% This Matlab file checks Figure 1 and the Supplementary Appendix S.1 in the paper "The Implementation of Stabilization Policy" (O. Loisel, June 2020).
% It calls two Dynare files: ISP_Fig1_check_ramsey.mod and ISP_Fig1_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_Fig1.m.
% The inputs are cc (referring to the sub-figure of Figure 1 considered), theta and rho (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 (32).
% 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]
theta = 0.5; % real-number value to be chosen in (-infinity,+infinity)
rho = 0.5; % real-number value to be chosen in(-1,1)

theta_eta = theta; % moving-average parameter for disturbance eta
theta_u = theta; % moving-average parameter for disturbance u
rho_eta = rho; % auto-regressive parameter for disturbance eta
rho_u = rho; % auto-regressive parameter for disturbance u
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','theta_u','theta_eta','sigma','kappa','lambda');

disp(' ');
disp('Step 1: compute the IRFs using the ramsey_policy command of Dynare');
dynare ISP_Fig1_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_Fig1_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 (32) 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_Fig1.m
mu = ((beta*lambda+kappa^2+lambda)-sqrt((beta*lambda+kappa^2+lambda)^2-4*beta*lambda^2))/(2*beta*lambda);
M = [beta*kappa^2+kappa*lambda/sigma+kappa^2+kappa^3/sigma+lambda, beta*kappa*(kappa-lambda/sigma), -beta^2*lambda;
    (beta*kappa+kappa^2/sigma+lambda/sigma)*kappa, beta*kappa*(kappa-lambda/sigma), beta*(-beta*lambda*mu+kappa^2+lambda);
    1, -(1+beta+kappa/sigma), beta];
A = inv(M)*[0; 0; kappa];
a_0 = A(1);
a_1 = A(2);
a_2 = A(3);
b = lambda*rho_u*(rho_u+theta_u)*(1-rho_u)/(beta*lambda*rho_u^2-(beta*lambda+kappa^2+lambda)*rho_u+lambda);
B = inv(M)*[lambda*(1+beta*(rho_u+theta_u)+kappa/sigma); lambda*(beta*rho_u*(rho_u+theta_u)-beta^2*(mu-rho_u)*b+kappa/sigma); 1-(rho_u+theta_u)];
b_0 = B(1);
b_1 = B(2);
b_2 = B(3);
TZ_0 = [a_0, b_0; (a_0-beta*a_1)/kappa, (b_0-beta*b_1-1)/kappa];
TZ_1 = [a_1, b_1; (a_1-beta*a_2)/kappa, (b_1-beta*b_2-(rho_u+theta_u))/kappa];
TZ_u = [0, b; 0, ((1-beta*rho_u)*b-rho_u*(rho_u+theta_u))/kappa];
TZ_mu = [a_2, b_2-b; (1-beta*mu)*a_2/kappa, (1-beta*mu)*(b_2-b)/kappa];
Ti_1 = [((1+beta-beta*mu+kappa/sigma)*a_2-a_1+kappa*(rho_eta+theta_eta))/(kappa/sigma), ((1+beta-beta*mu+kappa/sigma)*b_2-b_1+beta*(mu-rho_u)*b+(rho_u+theta_u)*(1-rho_u))/(kappa/sigma)];
Ti_eta = [rho_eta*(rho_eta+theta_eta)*sigma, 0];
Ti_u = [0, b*rho_u-(1-rho_u)*((1-beta*rho_u)*b-rho_u*(rho_u+theta_u))*sigma/kappa];
Ti_mu = [-(kappa*sigma-lambda)*a_2*mu/lambda, -(kappa*sigma-lambda)*(b_2-b)*mu/lambda];

TZ11 = [(rho_u*mu*TZ_1(1,1)-mu*TZ_u(1,1)-rho_u*TZ_mu(1,1)) (rho_u*mu*TZ_0(1,1)-(rho_u+mu)*TZ_1(1,1)+TZ_u(1,1)+TZ_mu(1,1)) (-(rho_u+mu)*TZ_0(1,1)+TZ_1(1,1)) TZ_0(1,1)];
TZ12 = [(rho_u*mu*TZ_1(1,2)-mu*TZ_u(1,2)-rho_u*TZ_mu(1,2)) (rho_u*mu*TZ_0(1,2)-(rho_u+mu)*TZ_1(1,2)+TZ_u(1,2)+TZ_mu(1,2)) (-(rho_u+mu)*TZ_0(1,2)+TZ_1(1,2)) TZ_0(1,2)];
TZ21 = [(rho_u*mu*TZ_1(2,1)-mu*TZ_u(2,1)-rho_u*TZ_mu(2,1)) (rho_u*mu*TZ_0(2,1)-(rho_u+mu)*TZ_1(2,1)+TZ_u(2,1)+TZ_mu(2,1)) (-(rho_u+mu)*TZ_0(2,1)+TZ_1(2,1)) TZ_0(2,1)];
TZ22 = [(rho_u*mu*TZ_1(2,2)-mu*TZ_u(2,2)-rho_u*TZ_mu(2,2)) (rho_u*mu*TZ_0(2,2)-(rho_u+mu)*TZ_1(2,2)+TZ_u(2,2)+TZ_mu(2,2)) (-(rho_u+mu)*TZ_0(2,2)+TZ_1(2,2)) TZ_0(2,2)];
Ti11 = [(-rho_eta*rho_u*mu*Ti_1(1,1)+rho_u*mu*Ti_eta(1,1)+rho_eta*mu*Ti_u(1,1)+rho_eta*rho_u*Ti_mu(1,1)) ((rho_eta*rho_u+rho_eta*mu+rho_u*mu)*Ti_1(1,1)-(rho_u+mu)*Ti_eta(1,1)-(rho_eta+mu)*Ti_u(1,1)-(rho_eta+rho_u)*Ti_mu(1,1)) (-(rho_eta+rho_u+mu)*Ti_1(1,1)+Ti_eta(1,1)+Ti_u(1,1)+Ti_mu(1,1)) Ti_1(1,1)];
Ti12 = [(-rho_eta*rho_u*mu*Ti_1(1,2)+rho_u*mu*Ti_eta(1,2)+rho_eta*mu*Ti_u(1,2)+rho_eta*rho_u*Ti_mu(1,2)) ((rho_eta*rho_u+rho_eta*mu+rho_u*mu)*Ti_1(1,2)-(rho_u+mu)*Ti_eta(1,2)-(rho_eta+mu)*Ti_u(1,2)-(rho_eta+rho_u)*Ti_mu(1,2)) (-(rho_eta+rho_u+mu)*Ti_1(1,2)+Ti_eta(1,2)+Ti_u(1,2)+Ti_mu(1,2)) Ti_1(1,2)];

G = conv([-rho_eta 1],conv(TZ11,TZ22)-conv(TZ12,TZ21));
F_pi = conv(Ti11,TZ22)-conv(Ti12,TZ21);
F_y = -conv(Ti11,TZ12)+conv(Ti12,TZ11);
% End of the part of the code that is copied and pasted from the file ISP_Fig1.m

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

save('ISP_Fig1_check.mat','F_pi_ext','F_y_ext','G_ext');
dynare ISP_Fig1_check_rule;
disp(' ');

if (isfield(oo_,'irfs'))
    disp('Blanchard and Kahn''s (1980) root-counting and no-decoupling conditions are met under the rule (32)');
    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_Fig1_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(' ');