% This Matlab file (dated September 2020) produces Figures 1 and G.1-G.5 in the paper "Pegging the Interest Rate on Bank Reserves:
% A Resolution of New Keynesian Puzzles and Paradoxes" (September 2020).

clear all;

% Calibration of the parameters of the basic NK model (Gali, 2008, chapter 3)
beta = 0.99; % quarterly discount factor
sigmatilde = 1; % relative risk aversion coefficient
epsilon = 6; % elasticity of substitution between differentiated goods
thetastar = 2/3; % benchmark probability of not being allowed to change one's price in a given quarter
eta = 1; % inverse of Frisch labor supply elasticity
alpha = 2/3; % labor elasticity of output

% Calibration of some other parameters
im = -0.0025; % change in quarterly interest rate on money (in percentage points per quarter)
Tmax = 21; % maximal horizon for forward guidance + 1
T = 1:Tmax;
g = 0.01; % change in government purchases expressed as a fraction of steady-state output
Tgmax = 21; % maximal horizon for government purchases + 1
Tg = 2:Tgmax;
Nmax = 5; % number of positive price-stickiness degrees considered

for cc=1:4
    if cc==1 % Zero steady-state government purchases, middle-of-the-range semi-elasticity (Gali, 2008, chapter 3)
        sigma = sigmatilde;
        deltag = 0;
        chiy = 1;
        chii = 4;
    end
    if cc==2 % Zero steady-state government purchases, low semi-elasticity
        sigma = sigmatilde;
        deltag = 0;
        chiy = 1;
        chii = 0.4;
    end
    if cc==3 % Zero steady-state government purchases, high semi-elasticity
        sigma = sigmatilde;
        deltag = 0;
        chiy = 1;
        chii = 33;
    end
    if cc==4 % Positive steady-state government purchases, middle-of-the-range semi-elasticity
        govery = 0.3; % steady-state ratio of government purchases to output
        sigma = sigmatilde/(1-govery);
        deltag = sigma / (sigma+(eta+1-alpha)/alpha) ;
        chiy = 1;
        chii = 4;
    end
    
    % Positive price-stickiness degrees
    for ii=1:Nmax
        N(ii) = 2^(ii-1);
    end
    theta = N*thetastar/(2^(Nmax-1));
    kappa = ((1-theta).*(1-beta*theta)./theta) * (alpha/(alpha+(1-alpha)*epsilon)) * (sigma+(eta+1-alpha)/alpha);
    P = [ones(Nmax,1)  -(2+1/beta+kappa.'/(beta*sigma)+chiy/(sigma*chii))  1+2/beta+(1+1/chii)*kappa.'/(beta*sigma)+(1+1/beta)*chiy/(sigma*chii)  -(1/beta+chiy/(beta*sigma*chii))*ones(Nmax,1)];
    for ii=1:Nmax
        R(:,ii) = sort(roots(P(ii,:)));
    end
    rho = R(1,:);
    omega1 = R(2,:);
    omega2 = R(3,:);
    xi1 = beta*(omega1+rho-1)-1;
    xi2 = beta*(omega2+rho-1)-1;
    for ii=1:Tmax
        pi(:,ii) = -kappa./(beta*sigma*(omega2-omega1)).*(omega1.^(-T(ii))-omega2.^(-T(ii)))*im;
        y(:,ii) = 1./(beta*sigma*(omega2-omega1)).*(xi1.*omega1.^(-T(ii))-xi2.*omega2.^(-T(ii)))*im;
    end
    inflation = 400*pi;
    output = 100*y;
    xig1 = (1-deltag).*(omega1-1) + chiy*deltag/(sigma*chii) ;
    xig2 = (1-deltag).*(omega2-1) + chiy*deltag/(sigma*chii) ;
    for ii=1:Tgmax-1
        pig(:,ii) = -kappa./(beta*(omega2-omega1)).*(xig1.*omega1.^(-Tg(ii))-xig2.*omega2.^(-Tg(ii)))*g;
        yg(:,ii) = 1./(beta.*(omega2-omega1)).*(xi1.*xig1.*omega1.^(-Tg(ii))-xi2.*xig2.*omega2.^(-Tg(ii)))*g;
    end
    inflationg = 400*pig;
    outputg = 100*yg;
    
    % Perfect price flexibility
    piflex = -(1+1/chii).^(-T)*im;
    inflationflex = 400*piflex;
    outputflex = zeros(1,Tmax);
    pigflex = -(chiy*deltag/chii + (1-deltag)*sigma/chii)*(1+1/chii).^(-Tg)*g;
    inflationgflex = 400*pigflex;
    outputgflex = zeros(1,Tgmax-1);

    % Basic NK model
    rhob = (1+beta+kappa./sigma-((1+beta+kappa./sigma).^2-4*beta).^(1/2))/(2*beta);
    omegab = (1+beta+kappa./sigma+((1+beta+kappa./sigma).^2-4*beta).^(1/2))/(2*beta);
    for ii=1:Tmax
        pink(:,ii) = -kappa*im./(beta*sigma*(omegab-rhob)).*(rhob.^(-T(ii))-omegab.^(-T(ii)));
        ynk(:,ii) = -kappa*im./(beta*sigma^2*(omegab-rhob)).*(rhob.^(-T(ii)+1)./(1-rhob)-omegab.^(-T(ii)+1)./(omegab-1));
    end
    inflationnk = 400*pink;
    outputnk = 100*ynk;
    outputnkflex = zeros(1,Tmax);
    for ii=1:Tgmax-1
        pinkg(:,ii) = kappa*(1-deltag)*g./(beta*(omegab-rhob)).*((1-rhob).*rhob.^(-Tg(ii))+(omegab-1).*omegab.^(-Tg(ii)));
        ynkg(:,ii) = kappa*(1-deltag)*g./(beta*sigma*(omegab-rhob)).*(rhob.^(-Tg(ii)+1)-omegab.^(-Tg(ii)+1));
    end
    inflationnkg = 400*pinkg;
    outputnkg = 100*ynkg;
    outputnkgflex = zeros(1,Tgmax-1);
    
    if cc==1
        inflation1 = inflation;
        output1 = output;
        inflationg1 = inflationg;
        outputg1 = outputg;
        inflationflex1 = inflationflex;
        outputflex1 = outputflex;
        inflationgflex1 = inflationgflex;
        outputgflex1 = outputgflex;
        inflationnk1 = inflationnk;
        outputnk1 = outputnk;
        inflationnkg1 = inflationnkg;
        outputnkg1 = outputnkg;
        outputnkflex1 = outputnkflex;
        outputnkgflex1 = outputnkgflex;
    end
    if cc==2
        inflation2 = inflation;
        output2 = output;
        inflationg2 = inflationg;
        outputg2 = outputg;
        inflationflex2 = inflationflex;
        outputflex2 = outputflex;
        inflationgflex2 = inflationgflex;
        outputgflex2 = outputgflex;
        inflationnk2 = inflationnk;
        outputnk2 = outputnk;
        inflationnkg2 = inflationnkg;
        outputnkg2 = outputnkg;
        outputnkflex2 = outputnkflex;
        outputnkgflex2 = outputnkgflex;
    end
    if cc==3
        inflation3 = inflation;
        output3 = output;
        inflationg3 = inflationg;
        outputg3 = outputg;
        inflationflex3 = inflationflex;
        outputflex3 = outputflex;
        inflationgflex3 = inflationgflex;
        outputgflex3 = outputgflex;
        inflationnk3 = inflationnk;
        outputnk3 = outputnk;
        inflationnkg3 = inflationnkg;
        outputnkg3 = outputnkg;
        outputnkflex3 = outputnkflex;
        outputnkgflex3 = outputnkgflex;
    end
    if cc==4
        inflation4 = inflation;
        output4 = output;
        inflationg4 = inflationg;
        outputg4 = outputg;
        inflationflex4 = inflationflex;
        outputflex4 = outputflex;
        inflationgflex4 = inflationgflex;
        outputgflex4 = outputgflex;
        inflationnk4 = inflationnk;
        outputnk4 = outputnk;
        inflationnkg4 = inflationnkg;
        outputnkg4 = outputnkg;
        outputnkflex4 = outputnkflex;
        outputnkgflex4 = outputnkgflex;
    end
end


% Figures

Hor = T-ones(1,Tmax);
Horg = Tg-ones(1,Tgmax-1);

f1=figure('Name','Effect of a policy-rate cut at date t+k on inflation at date t','PaperUnits','inches','PaperPosition',[-0.8 0 9.6 3],'PaperSize',[8 3]);
subplot(1,2,1);
plot(Hor,inflationflex1,'-ko',Hor,inflation1(1,:),':k',Hor,inflation1(2,:),'-.k',Hor,inflation1(3,:),'--k',Hor,inflation1(4,:),'-k',Hor,inflation1(5,:),'-k*');
title('Simple model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');
legend({'$\theta =0$','$\theta = \theta^\ast/16$','$\theta = \theta^\ast/8$','$\theta =\theta^\ast/4$','$\theta =\theta^\ast/2$','$\theta =\theta^\ast\equiv 2/3$'},'Interpreter','latex','FontSize',10,'Location','northeast');
subplot(1,2,2);
plot(Hor,inflationnk1(5,:),'-k*',Hor,inflationnk1(4,:),'-k',Hor,inflationnk1(3,:),'--k',Hor,inflationnk1(2,:),'-.k',Hor,inflationnk1(1,:),':k');
ylim([0 inflationnk1(5,Tmax)])
title('Basic NK model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');

fG1=figure('Name','Effect of a policy-rate cut at date t+k on output at date t','PaperUnits','inches','PaperPosition',[-0.71 0 9.6 3],'PaperSize',[8.05 3]);
subplot(1,2,1);
plot(Hor,output1(5,:),'-k*',Hor,output1(4,:),'-k',Hor,output1(3,:),'--k',Hor,output1(2,:),'-.k',Hor,output1(1,:),':k',Hor,outputflex1,'-ko');
title('Simple model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');
legend({'$\theta =\theta^\ast\equiv 2/3$','$\theta =\theta^\ast/2$','$\theta =\theta^\ast/4$','$\theta = \theta^\ast/8$','$\theta = \theta^\ast/16$','$\theta =0$'},'Interpreter','latex','FontSize',10,'Location','northeast');
subplot(1,2,2);
plot(Hor,outputnk1(5,:),'-k*',Hor,outputnk1(4,:),'-k',Hor,outputnk1(3,:),'--k',Hor,outputnk1(2,:),'-.k',Hor,outputnk1(1,:),':k',Hor,outputnkflex1,'-ko');
ylim([0 outputnk1(5,Tmax)])
title('Basic NK model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');

fG2=figure('Name','Effect of government purchases at date t+k on inflation at date t','PaperUnits','inches','PaperPosition',[-0.71 0 9.6 3],'PaperSize',[8.05 3]);
subplot(1,2,1);
plot(Horg,inflationg4(5,:),'-k*',Horg,inflationg4(4,:),'-k',Horg,inflationg4(3,:),'--k',Horg,inflationg4(2,:),'-.k',Horg,inflationg4(1,:),':k',Horg,inflationgflex4,'-ko');
title('Simple model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');
legend({'$\theta =\theta^\ast\equiv 2/3$','$\theta =\theta^\ast/2$','$\theta =\theta^\ast/4$','$\theta = \theta^\ast/8$','$\theta = \theta^\ast/16$','$\theta = 0$'},'Interpreter','latex','FontSize',10,'Location','southeast');
subplot(1,2,2);
plot(Horg,inflationnkg4(5,:),'-k*',Horg,inflationnkg4(4,:),'-k',Horg,inflationnkg4(3,:),'--k',Horg,inflationnkg4(2,:),'-.k',Horg,inflationnkg4(1,:),':k');
ylim([0 inflationnkg4(5,Tgmax-1)])
title('Basic NK model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');

fG3=figure('Name','Effect of government purchases at date t+k on output at date t','PaperUnits','inches','PaperPosition',[-0.67 0 9.6 3],'PaperSize',[8.1 3]);
subplot(1,2,1);
plot(Horg,outputgflex4,'-ko',Horg,outputg4(1,:),':k',Horg,outputg4(2,:),'-.k',Horg,outputg4(3,:),'--k',Horg,outputg4(4,:),'-k',Horg,outputg4(5,:),'-k*');
title('Simple model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');
legend({'$\theta =0$','$\theta = \theta^\ast/16$','$\theta = \theta^\ast/8$','$\theta =\theta^\ast/4$','$\theta =\theta^\ast/2$','$\theta =\theta^\ast\equiv 2/3$'},'Interpreter','latex','FontSize',10,'Location','southeast');
subplot(1,2,2);
plot(Horg,outputnkg4(5,:),'-k*',Horg,outputnkg4(4,:),'-k',Horg,outputnkg4(3,:),'--k',Horg,outputnkg4(2,:),'-.k',Horg,outputnkg4(1,:),':k',Horg,outputnkgflex4,'-ko');
ylim([0 outputnkg4(5,Tgmax-1)])
title('Basic NK model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');

fG4=figure('Name','Effect of a policy-rate cut at date t+k on inflation at date t for alternative values of chii','PaperUnits','inches','PaperPosition',[-0.67 0 9.6 3],'PaperSize',[8.1 3]);
subplot(1,3,1);
plot(Hor,inflationflex2,'-ko',Hor,inflation2(1,:),':k',Hor,inflation2(2,:),'-.k',Hor,inflation2(3,:),'--k',Hor,inflation2(4,:),'-k',Hor,inflation2(5,:),'-k*');
title('Simple model with $\chi_i=0.4$','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');
legend({'$\theta =0$','$\theta = \theta^\ast/16$','$\theta = \theta^\ast/8$','$\theta =\theta^\ast/4$','$\theta =\theta^\ast/2$','$\theta =\theta^\ast\equiv 2/3$'},'Interpreter','latex','FontSize',10,'Location','northeast');
subplot(1,3,2);
plot(Hor,inflationflex3(:),'-ko',Hor,inflation3(1,:),':k',Hor,inflation3(2,:),'-.k',Hor,inflation3(3,:),'--k',Hor,inflation3(4,:),'-k',Hor,inflation3(5,:),'-k*');
title('Simple model with $\chi_i=33$','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');
subplot(1,3,3);
plot(Hor,inflationnk1(5,:),'-k*',Hor,inflationnk1(4,:),'-k',Hor,inflationnk1(3,:),'--k',Hor,inflationnk1(2,:),'-.k',Hor,inflationnk1(1,:),':k');
ylim([0 inflationnk1(5,Tmax)])
title('Basic NK model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points per annum','Interpreter','latex');

fG5=figure('Name','Effect of a policy-rate cut at date t+k on output at date t for alternative values of chii','PaperUnits','inches','PaperPosition',[-0.67 0 9.6 3],'PaperSize',[8.1 3]);
subplot(1,3,1);
plot(Hor,output2(5,:),'-k*',Hor,output2(4,:),'-k',Hor,output2(3,:),'--k',Hor,output2(2,:),'-.k',Hor,output2(1,:),':k',Hor,outputflex2,'-ko');
title('Simple model with $\chi_i=0.4$','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');
legend({'$\theta =\theta^\ast\equiv 2/3$','$\theta =\theta^\ast/2$','$\theta =\theta^\ast/4$','$\theta = \theta^\ast/8$','$\theta = \theta^\ast/16$','$\theta =0$'},'Interpreter','latex','FontSize',10,'Location','northeast');
subplot(1,3,2);
plot(Hor,output3(5,:),'-k*',Hor,output3(4,:),'-k',Hor,output3(3,:),'--k',Hor,output3(2,:),'-.k',Hor,output3(1,:),':k',Hor,outputflex3,'-ko');
title('Simple model with $\chi_i=33$','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');
subplot(1,3,3);
plot(Hor,outputnk1(5,:),'-k*',Hor,outputnk1(4,:),'-k',Hor,outputnk1(3,:),'--k',Hor,outputnk1(2,:),'-.k',Hor,outputnk1(1,:),':k',Hor,outputnkflex1,'-ko');
ylim([0 outputnk1(5,Tmax)])
title('Basic NK model','Interpreter','latex');
xlabel('$k$ in quarters','Interpreter','latex');
ylabel('Percentage points','Interpreter','latex');



print(f1,'-dpdf','fig1');
print(fG1,'-dpdf','figG1');
print(fG2,'-dpdf','figG2');
print(fG3,'-dpdf','figG3');
print(fG4,'-dpdf','figG4');
print(fG5,'-dpdf','figG5');
