% Author: Olivier Loisel (CREST). Date: June 2020.
% This .m file produces Figure 1 in the paper "The Implementation of Stabilization Policy" (O. Loisel, June 2020).
% The names of the objects that appear also in the paper (scalars, vectors, matrices) correspond to their notations in the paper.
% Figure 1 in the paper is obtained when nn is set to 1000. The run time may then be relatively long (on my computer: about 13 minutes for mm=0 and 48 minutes for mm=1).
% When nn is set to 100 (as below), the run time is about 100 times shorter, but the figure produced is then less accurate than Figure 1 in the paper.
% The papers cited in this file but not in the paper are referenced at the end of this file.

nn = 100; % number of alternative values considered for theta and rho (so that the number of points considered in the plane of each sub-figure of Figure 1 is nn^2)
mm = 0; % resolution method chosen: mm=0 for the method of Blanchard and Kahn (1980) and Klein (2000), mm=1 for the method of Sims (2001);
        % choosing mm=1 requires to download first the gensys.m file, which is available at http://sims.princeton.edu/yftp/gensys/

beta = 0.99; % discount factor
theta = 0.67; % probability of not being allowed to change one's price

for cc=0:1,
    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)
    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);
    MM=ones(nn);
    for jj=1:nn,
        rho_eta = -1 + 2*jj/(nn+1); % auto-regressive parameter for disturbance eta
        rho_u = rho_eta; % auto-regressive parameter for disturbance u
        for ii=1:nn,
            theta_eta = (-1+2*ii/(nn+1))/(1-abs(-1+2*ii/(nn+1))); % moving-average parameter for disturbance eta
            theta_u = theta_eta; % moving-average parameter for disturbance u
            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);
            
            if (mm==0)
                % Write the system in Blanchard and Kahn's (1980) form M1*E_t(X_{t+1}) = M0*X_t + exogenous term,
                % where X_t = [pi_{t-1}; ...; pi_{t-p}; y_{t-1}; ...; y_{t-p}; i_{t-1}; ...; i_{t-p}; pi_t; y_t]
                p = length(G)-1;
                M1 = [zeros(1,3*p), beta, 0;
                    zeros(1,2*p), -1/sigma, zeros(1,p-1), 1/sigma, 1;
                    zeros(1,2*p), G(p+1), zeros(1,p+1);
                    eye(2*p), zeros(2*p,p+2);
                    zeros(p-1,2*p+1), eye(p-1), zeros(p-1,2)];
                M0 = [zeros(1,3*p), 1, -kappa;
                    zeros(1,3*p+1), 1;
                    fliplr(F_pi), fliplr(F_y), -fliplr(G(1:p)), zeros(1,2);
                    zeros(1,3*p), 1, 0;
                    eye(p-1), zeros(p-1,2*p+3);
                    zeros(1,3*p+1), 1;
                    zeros(p-1,p), eye(p-1), zeros(p-1,p+3);
                    zeros(p-1,2*p), eye(p-1), zeros(p-1,3)];
                BK = M1\M0;
                V = sort(abs(eig(BK)),'descend');
                % Determine whether Blanchard and Kahn's (1980) root-counting condition is met.
                if (V(2)<1)
                    MM(jj,ii)=0.7; % In this case, there is at least multiplicity.
                elseif (V(3)>1)
                    MM(jj,ii)=0.4; % In this case, there is at least non-robustness.
                else
                    % When Blanchard and Kahn's (1980) root-counting condition is met, determine whether their no-decoupling condition is met.
                    % Due to the well known numerical instability of the Jordan decomposition, I follow Klein (2000) and use the generalized Schur decomposition.
                    [AA,BB,QQ,ZZ] = qz(M1,M0);
                    [AAS,BBS,QQS,ZZS] = ordqz(AA,BB,QQ,ZZ,'udo');
                    if (rank([eye(3*p), zeros(3*p,2)]*ZZS*[eye(3*p); zeros(2,3*p)])~=3*p)
                        MM(jj,ii)=0; % In this case, there is multiplicity and non-robustness.
                    end
                end
            end
            
            if (mm==1)
                % Write the system in Sims' (2001) form g0*Z_t = g1*Z_{t-1} + pi*eta_t + exogenous term,
                % where Z_t = [pi_t; ...; pi_{t-p+1}; y_t; ...; y_{t-p+1}; i_t; ...; i_{t-p+1}] and eta_t is an expectational error
                g0 = eye(3*p);
                g0(1,1) = beta;
                g0(p+1,1) = 1/sigma;
                g0(2*p+1,2*p+1) = G(p+1);
                g1 = [1, zeros(1,p-1), -kappa, zeros(1,2*p-1);
                    eye(p-1), zeros(p-1,2*p+1);
                    zeros(1,p), 1, zeros(1,p-1), 1/sigma, zeros(1,p-1);
                    zeros(p-1,p), eye(p-1), zeros(p-1,p+1);
                    fliplr(F_pi), fliplr(F_y), -fliplr(G(1:p));
                    zeros(p-1,2*p), eye(p-1), zeros(p-1,1)];
                c = zeros(3*p,1);
                psi = zeros(3*p,1);
                pi = [1, 0;
                    zeros(p-1,2);
                    0, 1;
                    zeros(2*p-1,2)];
                [G1,C,impact,fmat,fwt,ywt,gev,eu,loose]=gensys(g0,g1,c,psi,pi);
                if (eu==[1;0])
                    MM(jj,ii)=0.7; % In this case, there is multiplicity.
                elseif (eu==[0;1])
                    MM(jj,ii)=0.4; % In this case, there is non-robustness.
                elseif (eu==[0;0])
                    MM(jj,ii)=0; % In this case, there is multiplicity and non-robustness.
                end
            end
            
        end
    end
    if cc==0
        MMG=MM;
    end
    if cc==1
        MMW=MM;
    end
end

f=figure('Name','Figure 1');

subplot(1,2,1);
contourf(MMG,[0, 0.4, 0.7, 1]);
colormap(gray);
title('Gal\''i''s calibration','Interpreter','latex');
caxis([-1, 1]);
xlabel('$\theta(1+\vert\theta\vert)^{-1}$','Interpreter','latex','FontSize',10);
set(gca,'XTickLabel',{'-1';'-1/2';'0';'1/2';'1'},'XTick',[1, nn/4, nn/2, 3*nn/4, nn],'Xgrid','on','FontSize',9);
ylabel('$\rho$','Interpreter','latex','FontSize',10);
set(gca,'YTickLabel',{'-1';'0';'1'},'YTick',[1, nn/2, nn],'Ygrid','on','FontSize',9);
axis square;

subplot(1,2,2);
contourf(MMW,[0, 0.4, 0.7, 1]);
colormap(gray);
title('Woodford''s calibration','Interpreter','latex');
caxis([-1, 1]);
xlabel('$\theta(1+\vert\theta\vert)^{-1}$','Interpreter','latex','FontSize',10);
set(gca,'XTickLabel',{'-1';'-1/2';'0';'1/2';'1'},'XTick',[1, nn/4, nn/2, 3*nn/4, nn],'Xgrid','on','FontSize',9);
ylabel('$\rho$','Interpreter','latex','FontSize',10);
set(gca,'YTickLabel',{'-1';'0';'1'},'YTick',[1, nn/2, nn],'Ygrid','on','FontSize',9);
axis square;

print(f,'-dpdf','fig1');

% Additional Reference:
% Klein, P. (2000), ''Using the Generalized Schur Form to Solve a Multivariate Linear Rational Expectations Model,'' Journal of Economic Dynamics and Control, 24, 1405-1423.