%sticky_2
% does a simulation of the effects of monetary policy with k period price
% stckiness

clear all
close all

% baseline parameters
Bbar = 100;
Pbar = 100; 
r = 0.05; % delta in the paper 
beta = 1/(1+r); 
sbar = Bbar/Pbar*(1-beta);
shifter = 0.05; % shifts lines up and down so you can see them


% intervention
ichange = 0.01; 
Delta = (1+r+ichange)/(1+r);
rho = 0.9; % interest rate decay factor
Deltat = 1+(Delta-1)*rho.^(0:6)';
cpDeltat = cumprod(Deltat); 
sp = 1; % no fiscal shock 
create_time_series; % script to create time series

% --------------------
% plot pure monetary 
% -------------------

T = [0:7]';

figure;
hold on;
plot(T,100*it-5+shifter,'-vb','linewidth',2);
plot(T,(Pt-100)+2*shifter,'-ro','linewidth',2);
plot(T,(ct-1)*100+3*shifter,'--vk','linewidth',2);
plot(T,100*rt-5+4*shifter,'--sg','linewidth',2);
legend('i','P','c','r-5',4);
text(5.5,4,'P','color','r','fontsize',16);
text(3.5,2.2,'c','color','k','horizontalalignment','l','fontsize',16);
text(3.2,-2,'r-5','color','g','horizontalalignment','l','fontsize',16);
text(5.5,1.3,'i','color','b','fontsize',16);
xlabel('Time')
ylabel('Percent')
print -depsc2 sticky_2_plot.eps
print -dpng sticky_2_plot.png

%-----------------------
% 2.  Pure fiscal. Dropped from paper but fun to look at 
% ----------------------

Deltat = ones(7,1);
cpDeltat = cumprod(Deltat); 
sp = 1.02; % poure fiscal shock 
create_time_series; 

% plot
T = [0:7]';

% p = beta^n u'(t+n)/u'(t) so y = p^(-1/n) = 1/beta * u'(t)/u'(t+n)^(1/n)
longrates = 1/beta*[ [ 1; (uct(2:end)./[uct(3:end);1])]...
                     [ 1; (uct(2:end)./[uct(4:end);1;1]).^(1/2)]...
                     [ 1; (uct(2:end)./[uct(5:end);1;1;1]).^(1/3)]];


figure;
hold on;
plot(T,100*it-5,'-vb','linewidth',2);
plot(T,(Pt-100),'-ro','linewidth',2);
plot(T,(ct-1)*100,'--vk','linewidth',2);
plot(T,100*rt-5,'--sg','linewidth',2);
plot(T,100*(longrates-1)-5,'--sg','linewidth',2);
legend('i','P','c','r-5',4);
text(5.5,4,'P','color','r','fontsize',16);
text(3.5,2.2,'c','color','k','horizontalalignment','l','fontsize',16);
text(3.2,-2,'r-5','color','g','horizontalalignment','l','fontsize',16);
text(5.5,1.3,'i','color','b','fontsize',16);
xlabel('time')
ylabel('percent')
print -depsc2 sticky_2_plot_with_fiscal.eps
print -dpng sticky_2_plot_with_fiscal.png


% ---------------------
% mixed fiscal-monetary 
% ---------------------

figure;

subplot(2,2,1)
% first no fiscal policy

ichange = 0.01; 
x = 0.00; % final fall in rates
rho = 1; % interest rate decay factor
max = 2; % turn off interest rate change

% the program allows for monetary policy to eventually reverse, if you want
% to create eventual disinflation. For that, set a < 1. The paper ended up
% not using all the generality allowed here
%  Delta t = a + b*rho^(t-1). 
%  Want impact rise of 0.01, so (1+r)*Delta_1 = (1+r+0.01) so (1+r)*(a+b)=(1+r+0.01).
%  Want final fall of 0.0x, so (1+r)*a = (1+r-0.0x); 
% thus, a = (1+r-0.0x)/(1+r), and 
% (1+r-0.0x)/(1+r) + b = (1+r+0.01)/(1+r)
%  b = (1+r+0.01 - (1+r-0.0x))/(1+r)=(0.0x+0.01)/(1+r)

a = (1+r-x)/(1+r);
b = (x+ichange)/(1+r); 
Deltat = a + b*rho.^(0:6)'; % Deltat goes from 1 on, other series have a time 0 
Deltat(max+1:end)=1; 

cpDeltat = cumprod(Deltat); 
sp = 1; % fiscal shock 
create_time_series; % script to create time series

hold on;
plot(T,100*it-5,'-vb','linewidth',2);
plot(T,(Pt-100)+shifter,'-ro','linewidth',2);
plot(T,(ct-1)*100+2*shifter,'--vk','linewidth',2);
plot(T,100*rt-5+3*shifter,'--sg','linewidth',2);
legend('i','P','c','r-5',4);
xlabel('Time')
ylabel('Percent')


subplot(2,2,2)
% now, just fiscal policy

ichange = 0.00; 
x = 0.00; % final fall in rates
rho = 1; % interest rate decay factor
max = 2; % turn off interest rate change


a = (1+r-x)/(1+r);
b = (x+ichange)/(1+r); 
Deltat = a + b*rho.^(0:6)'; % Deltat goes from 1 on, other series have a time 0 
Deltat(max+1:end)=1; 

cpDeltat = cumprod(Deltat); 
sp = 1.03; % fiscal shock 
create_time_series; % script to create time series

hold on;
plot(T,100*it-5,'-vb','linewidth',2);
plot(T,(Pt-100)+shifter,'-ro','linewidth',2);
plot(T,(ct-1)*100+2*shifter,'--vk','linewidth',2);
plot(T,100*rt-5+3*shifter,'--sg','linewidth',2);
%axis([-3 3 70 120]);  
legend('i','P','c','r-5',1);
xlabel('Time')
ylabel('Percent')



subplot(2,2,3)
% now, fiscal policy too

ichange = 0.01; 
x = 0.00; % final fall in rates
rho = 1; % interest rate decay factor
max = 2; % turn off interest rate change

a = (1+r-x)/(1+r);
b = (x+ichange)/(1+r); 
Deltat = a + b*rho.^(0:6)'; % Deltat goes from 1 on, other series have a time 0 
Deltat(max+1:end)=1; 

cpDeltat = cumprod(Deltat); 
sp = 1.02; % fiscal shock 
create_time_series; % script to create time series

hold on;
plot(T,100*it-5,'-vb','linewidth',2);
plot(T,(Pt-100)+shifter,'-ro','linewidth',2);
plot(T,(ct-1)*100+2*shifter,'--vk','linewidth',2);
plot(T,100*rt-5+3*shifter,'--sg','linewidth',2);
%axis([-3 3 70 120]);  
legend('i','P','c','r-5',4);
xlabel('Time')
ylabel('Percent')


subplot(2,2,4)
% now, more fiscal policy too

ichange = 0.01; 
x = 0.00; % final fall in rates
rho = 1; % interest rate decay factor
max = 2; % turn off interest rate change
a = (1+r-x)/(1+r);
b = (x+ichange)/(1+r); 
Deltat = a + b*rho.^(0:6)'; % Deltat goes from 1 on, other series have a time 0 
Deltat(max+1:end)=1; 

cpDeltat = cumprod(Deltat); 
sp = 1.03; % fiscal shock 
create_time_series; % script to create time series

hold on;
plot(T,100*it-5,'-vb','linewidth',2);
plot(T,(Pt-100)+shifter,'-ro','linewidth',2);
plot(T,(ct-1)*100+2*shifter,'--vk','linewidth',2);
plot(T,100*rt-5+3*shifter,'--sg','linewidth',2);
%axis([-3 3 70 120]);  
legend('i','P','c','r-5',4);
xlabel('Time')
ylabel('Percent')

print -depsc2 sticky_2_plot_with_fiscal_and_m.eps
print -dpng sticky_2_plot_with_fiscal_and_m.png


stop;
% script create_time_series. Cut out and save as a separate file. 
% create time series
      %  before      shock                                          k                     |   free
      %  0            1            2                   3               4                  5                      6      
it =  [ 1+r     (1+r)*Deltat(1)  (1+r)*Deltat(2)   (1+r)*Deltat(3)   (1+r)*Deltat(4)       (1+r)*Deltat(5)      (1+r)*Deltat(6)      (1+r)*Deltat(7)   ]'-1;
uct = [ 1        1*sp             1/cpDeltat(1)*sp 1/cpDeltat(2)*sp   1/cpDeltat(3)*sp      1                   1                   1                 ]';
rt =  [ 1+r     (1+r)*Deltat(1)  (1+r)*Deltat(2)   (1+r)*Deltat(3)   (1+r)/cpDeltat(3)*sp  1+r                  1+r                 1+r               ]'-1;
Pt =  [ Pbar     Pbar             Pbar              Pbar              Pbar                 cpDeltat(4)*Pbar/sp  cpDeltat(5)*Pbar/sp cpDeltat(6)*Pbar/sp  ]';
pit = [ 1        1                1                 1                 1                    cpDeltat(4)/sp       Deltat(5)           Deltat(6)         ]'-1;
Bt =  [ 1        cpDeltat(1)      cpDeltat(2)       cpDeltat(3)       cpDeltat(4)          cpDeltat(5)          cpDeltat(6)         cpDeltat(7)       ]';


ct = 1./uct; % log utility; 

