% ch8fig1_2.m
%
% Interest and Prices by M. Woodford
% ==================================
%
% Generates Figures 8.1 - 8.2
% requires file vardat.mat
%
% Computes weights of Target Criteria 1 and 2
% developed in "Optimal Inflation Targeting Rules" by
% M. Giannoni and M. Woodford.
%
% Marc Giannoni, 1/18/03; latest revision 4/17/03
% Similar to tc.m done for "Optimal Inflation Targeting Rules"

% ****************************************************************
% [0] Required Data file
% ****************************************************************

% addpath c:\Documents\Research\opint\NBER\estim\var
load vardat.mat  


% ****************************************************************
% [1] Estimated and implied parameters
% ****************************************************************

format compact

disp('Calibrated parameters:')
b = .99    % beta
ci = 0     % chi_i = chi*eta_i/(1-b*eta)
wp = 1/3   % omega_p
ap = .66
aw = .66
phi = 4/3

disp(' ')
disp('Estimated parameters:')
vp = 0.74839
xip = 0.0020355 
xiw = 0.0041512 
ww  = 19.5507 
eta = 1
gp = 1
gw = 1

disp(' ')
disp('Implied parameters:')
kp = xip*wp
w = wp+ww
nu = ww/phi
tp = ((1-ap)*(1-ap*b)/ap-xip)/(xip*wp)  % theta_p
tw = ((1-aw)*(1-aw*b)/aw-xiw)/(xiw*nu)  % theta_w

mup = tp/(tp-1)
muw = tw/(tw-1)

% Loss function:
disp(' ')
disp('Weights in loss function:')
b3 = (w+vp*(1+b*eta^2))/(b*vp); % b2*eta
ed =b/2*(b3+sqrt(b3^2-(4*eta^2)/b)); % eta/d
d = eta/ed;
d0 = ed*vp/2;

denom = tp/xip + tw/(phi*xiw);

Lp = (tp/xip)/denom
Lw = (tw/(phi*xiw))/denom
Lx = (2*d0)/denom
disp(['Lx*16 = ' num2str(Lx*16)])

Li = 0 % lambda_i

format

% ****************************************************************
% [2] Optimal Target Criterion
% ****************************************************************

xout = 1; %input('Output gap: (1) x  or (2) Y ? ');

J = 100; %input('Horizon = ');  % horizon
istar = 0; % i*
xstar = 0; % x*

% Construct polynomials
% AL = (s-b)*(1-s)+(xip+xiw)*s
AL = [-1 (1+b+xip+xiw) -b];
rAL = roots(AL);
mu1 = min(diag(inv(diag(rAL))));    % smallest of inverse of roots of AL
mu2 = max(diag(inv(diag(rAL))));    % largest of inverse of roots of AL

B_1 = vp*(eta*b-1)*(1-eta*1)-ww*1;
B_mu1inv = vp*(eta*b-1/mu1)*(1-eta/mu1)-ww/mu1;

% alpha1 and alpha2
a1 = xiw*(xip*B_1-kp)/(xip+xiw);
a2 = (xiw*mu1^2*B_mu1inv+kp*mu1)/(xip+xiw);

% C(L)
c0 = eta*b^2;
c1 = -(b+eta*b+eta^2*b^2);
c2 = 1+eta*b*(1+eta);
c3 = -eta;

% F(L)=f0+f1*L
f0 = xiw*xip*vp*eta*b/(xip+xiw);
f1 = -xiw*(xip*ww+xip*vp*(1-eta*b+eta^2*b)+kp)/(xip+xiw);

% G(L)=g0+g1*L
g0 = xiw*vp*eta*b/(xip+xiw);
g1 = -(xiw*ww+xiw*vp*(1-eta*b*mu1+eta^2*b)-kp)/(xip+xiw);

% Coefficients of f(t), g(t), e(t), and h(t)

% Let f(t) = fc + fpit E(t)[PI] + fwt E(t)[PIW] + fi E(t)[R]
% compute parameters of f(t-2), f(t-1)
fc = -Li*istar/(1-b);
% Param of f(t-2)
fpitm2 = zeros(1,J+7);
fwtm2 = zeros(1,J+7);
fitm2 = zeros(1,J+7);

fpitm2(2) = -b*(Lp+Lw); fpitm2(3) = b*Lp;
fwtm2(3) = b*Lw;
for j=0:J+5
    fitm2(2+j) = Li*b^j;
end

% Param of f(t-1)
fpitm1 = [0 fpitm2(1,1:end-1)];
fwtm1 = [0 fwtm2(1,1:end-1)];
fitm1 = [0 fitm2(1,1:end-1)];

% Param of f(t)
fpit = [0 fpitm1(1,1:end-1)];
fwt = [0 fwtm1(1,1:end-1)];
fit = [0 fitm1(1,1:end-1)];

% Let g(t) = gc + gpit E(t)[PI] + gwt E(t)[PIW] + gi E(t)[R]
% compute parameters of g(t-2), g(t-1)
gc = Li*istar*xip*mu1/(1-1/mu2);
% Param of g(t-2)
gpitm2 = zeros(1,J+7);
gwtm2 = zeros(1,J+7);
gitm2 = zeros(1,J+7);

gpitm2(2) = -b*mu1*(xiw*Lw-xip*Lp); gpitm2(3) = -b*mu1*(xip*Lp+(xiw*Lw-xip*Lp)*(1/mu2-b));
gwtm2(3) = b*mu1*xiw*Lw;
gitm2(2) = -Li*mu1*xip; gitm2(3) = -Li*xip*mu1/mu2;
for j=0:J+3
    gpitm2(4+j) = -b*mu1*(xip*Lp+(xiw*Lw-xip*Lp)/mu2)*(1/mu2-b)*(1/mu2)^j;
    gwtm2(4+j) = b*mu1*xiw*Lw*(1/mu2-b)*(1/mu2)^j;
    gitm2(4+j) = -Li*xip*mu1/(mu2^2)*(1/mu2)^j;
end

% Param of g(t-1)
gpitm1 = [0 gpitm2(1,1:end-1)];
gwtm1 = [0 gwtm2(1,1:end-1)];
gitm1 = [0 gitm2(1,1:end-1)];

% Param of g(t)
gpit = [0 gpitm1(1,1:end-1)];
gwt = [0 gwtm1(1,1:end-1)];
git = [0 gitm1(1,1:end-1)];

% Compute e(t)
ec = -b*Lx*(1-b*d)*xstar-vp*Li*(c0+c1+c2+c3)*istar;
% Param of e(t-2)
extm2 = zeros(1,J+7);
eitm2 = zeros(1,J+7);

extm2(1,2:4)= b*Lx*[-d, (1+b*d^2), -b*d];
eitm2(1,1:4)= vp*Li*[c3, c2, c1, c0];

% Compute param of h(t-2)
hpitm2 = -(f0*fpitm1+f1*fpitm2+g0*gpitm1+g1*gpitm2);
hwtm2 = -(f0*fwtm1+f1*fwtm2+g0*gwtm1+g1*gwtm2);
hxtm2 = extm2;
hitm2 = eitm2-(f0*fitm1+f1*fitm2+g0*gitm1+g1*gitm2);

% Param of h(t-1)
hpitm1 = [0 hpitm2(1,1:end-1)];
hwtm1 = [0 hwtm2(1,1:end-1)];
hxtm1 = [0 hxtm2(1,1:end-1)];
hitm1 = [0 hitm2(1,1:end-1)];
% Param of h(t)
hpit = [0 hpitm1(1,1:end-1)];
hwt = [0 hwtm1(1,1:end-1)];
hxt = [0 hxtm1(1,1:end-1)];
hit = [0 hitm1(1,1:end-1)];
% Param of h(t+1)
hpitp1 = [0 hpit(1,1:end-1)];
hwtp1 = [0 hwt(1,1:end-1)];
hxtp1 = [0 hxt(1,1:end-1)];
hitp1 = [0 hit(1,1:end-1)];

% Policy rule 11 (case Li=0)
% --------------
Pc = -a1*(1-mu1)*fc;
Ppi0 = hpitp1-a1*fpit-a2*gpit;
Pw0 = hwtp1-a1*fwt-a2*gwt;
Px0 = hxtp1;
% Pi0 = hit;

% Coeff wrt E(t-1)
Ppi1 = -(1+mu1)*hpit+a1*mu1*fpitm1+a2*gpitm1;
Pw1 = -(1+mu1)*hwt+a1*mu1*fwtm1+a2*gwtm1;
Px1 = -(1+mu1)*hxt;
% Pi1 = -((1+mu1)*hitm1+a1*fitm1+a2*gitm1);

% Coeff wrt E(t-2)
Ppi2 = mu1*hpitm1;
Pw2 = mu1*hwtm1;
Px2 = mu1*hxtm1;
% Pi2 = mu1*hitm2+a1*mu1*fitm2+a2*gitm2;


% ****************************************************************
% [3] Target criterion 1: in terms of inflation and level of real wage
% ****************************************************************

% LHS
phipi0 = sum(Ppi0(1,5:J+5)+Pw0(1,5:J+5));
apisTC1 = [Ppi0(1,5:J+5)+Pw0(1,5:J+5)]/phipi0;
phiw0 = sum(Pw0(1,5:J+5)-Pw0(1,6:J+6));
awsTC1 = [Pw0(1,5:J+5)-Pw0(1,6:J+6)]/phiw0;
phix0 = sum(Px0(1,5:end));
axsTC1 = Px0(1,5:end)/phix0;

% RHS
phipi1 = -sum([Ppi0(1,4)+Ppi1(1,4)+Pw1(1,4), Ppi1(1,5:J+5)+Pw1(1,5:J+5)]);
api1TC1 = -[Ppi0(1,4)+Ppi1(1,4)+Pw1(1,4), Ppi1(1,5:J+5)+Pw1(1,5:J+5)]/phipi1;
phiw1 = -sum([-Pw0(1,5)+Pw1(1,4)-Pw1(1,5), Pw1(1,5:J+5)-Pw1(1,6:J+6)]);
aw1TC1 = -[-Pw0(1,5)+Pw1(1,4)-Pw1(1,5), Pw1(1,5:J+5)-Pw1(1,6:J+6)]/phiw1;
phix1 = -sum(Px1(1,4:end));
ax1TC1 = -Px1(1,4:end)/phix1;

phipi2 = -sum([Ppi1(1,3)+Ppi2(1,3)+Pw2(1,3), Ppi2(1,4:J+4)+Pw2(1,4:J+4)]);
api2TC1 = -[Ppi1(1,3)+Ppi2(1,3)+Pw2(1,3), Ppi2(1,4:J+4)+Pw2(1,4:J+4)]/phipi2;
phiw2 = -sum([-Pw1(1,4)+Pw2(1,3)-Pw2(1,4), Pw2(1,4:J+4)-Pw2(1,5:J+5)]);
aw2TC1 = -[-Pw1(1,4)+Pw2(1,3)-Pw2(1,4), Pw2(1,4:J+4)-Pw2(1,5:J+5)]/phiw2;
phix2 = -sum(Px2(1,3:end));
ax2TC1 = -Px2(1,3:end)/phix2;

% normalization: phipis=1
norTC1 = phipi0;
phipiTC1 = phipi0/norTC1;
phiwTC1 = phiw0/norTC1;
phixTC1 = phix0/norTC1;

thpiTC1 = phipi1/phipi0-1;
thwTC1 = (phiw1-phiw0)/phipi0;
thxTC1 = (phix1-phix0)/phipi0;

disp(' ')
disp('Coefficients of Optimal Target Criterion (1) (long-term):')
disp('Representation in terms of pi(t) and w(t)')
disp(' ')
disp('    phipiTC1   phiwTC1*4    phixTC1*4')
disp([    phipiTC1   phiwTC1*4    phixTC1*4])
disp(' ')
disp('    thpiTC1    thwTC1*4     thxTC1*4')
disp([    thpiTC1    thwTC1*4     thxTC1*4])
disp(' ')

% Generate figure of az_i
Jpl = 6;
xax = 1:Jpl;
% Standard position of subplots given by vector 
%         [left bottom width height]
% 331:    0.1300    0.7012    0.2128    0.2238
% 332:    0.4111    0.7012    0.2128    0.2238
% 333:    0.6922    0.7012    0.2128    0.2238
% 334:    0.1300    0.4056    0.2128    0.2238
% 335:    0.4111    0.4056    0.2128    0.2238
% 336:    0.6922    0.4056    0.2128    0.2238
% 337:    0.1300    0.1100    0.2128    0.2238
% 338:    0.4111    0.1100    0.2128    0.2238
% 339:    0.6922    0.1100    0.2128    0.2238
 
figure('Name', 'Weights in Target criterion 1 (in terms of inflation and level of real wage)')
% subplot(331)
subplot('Position',[0.1300    0.72    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, apisTC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{*\pi}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
% subplot(332)
subplot('Position',[0.4111    0.72    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, awsTC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{*{\itw}}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
% subplot(333)
subplot('Position',[0.6922    0.72    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, axsTC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{*{\itx}}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')

% subplot(334)
subplot('Position',[0.1300    0.39    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, api1TC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{\pi1}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
% subplot(335)
subplot('Position',[0.4111    0.39    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, aw1TC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{{\itw}1}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
% subplot(336)
subplot('Position',[0.6922    0.39    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, ax1TC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{{\itx}1}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')

% subplot(337)
subplot('Position',[0.1300    0.06    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, api2TC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{\pi2}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
% subplot(338)
subplot('Position',[0.4111    0.06    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, aw2TC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{{\itw}2}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
% subplot(339)
subplot('Position',[0.6922    0.06    0.2128    0.2]) 
    plot(xax, zeros(1,Jpl), 'k:', xax, ax2TC1(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{{\itx}2}_{\itk}', 'FontName', 'Times')
%     get(gca,'Position')
print figalphaTC1.eps -deps
print figalphaTC1.ps -dps

% ****************************************************************
% [5] Target Criterion 2 (Short-run)
% ****************************************************************

% Compute coefficients of: 
% TC2pit * E(t)[PI] + TC2wt * E(t)[PIW] = TC2pit * E(t-1)[PI] + TC2wt * E(t-1)[PIW]
TC2pit = xip*fpit(1,5:J+5) + gpit(1,5:J+5);
TC2wt = xip*fwt(1,5:J+5) + gwt(1,5:J+5);

% Normalize coefficients
phipiTC2 = sum(TC2pit);
apiTC2 = TC2pit/phipiTC2;
phiwTC2 = sum(TC2wt);
awTC2 = TC2wt/phiwTC2;

norTC2 = phipiTC2 + phiwTC2;
phipiTC2 = phipiTC2/norTC2;
phiwTC2 = phiwTC2/norTC2;

% Alternative representation: phipiTC2b Ft(pi) + phiwTC2b Ft(w)
apiTC2b = phipiTC2*apiTC2+phiwTC2*awTC2; % sum(apiTC2b)=phipiTC2+phiwTC2=1
phipiTC2b = (phipiTC2+phiwTC2);
phiwTC2b = sum([phiwTC2*(awTC2(1:end)-[awTC2(2:end) 0])]');
awTC2b = phiwTC2*(awTC2(1:end)-[awTC2(2:end) 0])/phiwTC2b;

% Mean lead
mleadpi = sum(apiTC2b.*[1:size(apiTC2b,2)]);
mleadw = sum(awTC2b.*[1:size(awTC2b,2)]);

disp(' ')
disp('Coefficients of Optimal Target Criterion 2 (Short run)')
disp('Representation in terms of pi(t) and w(t)')
disp(' ')
disp('    phipiTC2b  phiwTC2b*4')
disp([phipiTC2b   phiwTC2b*4])
disp(' ')
disp('Mean lead of weights on projections:')
disp(' ')
disp('   mleadpi   mleadw')
disp([   mleadpi   mleadw])

Jpl = 6;
xax1 = 1:Jpl;

figure('Name', 'Weights in Target criterion 2 (in terms of pi(t) and w(t))')
    subplot(121)
        plot(xax1, zeros(1,Jpl), 'k:', xax1, apiTC2b(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{\pi}_{\itk}', 'FontName', 'Times')
    subplot(122)
        plot(xax1, zeros(1,Jpl), 'k:', xax1, awTC2b(1:Jpl), 'k-'), set(gca,'XLim', [1 Jpl]), title('\alpha^{\itw}_{\itk}', 'FontName', 'Times')
print figalphaTC2.eps -deps
print figalphaTC2.ps -dps