% PHILLIPS.M
clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data construction:

% Averaging CPI by quarter, then compute annualized geowth rate
% as in CG (2015)
load infl_q.txt; infl=infl_q*4;

% 1960Q1-2020Q1, MSC, one year inflation expectations
load msc_q.txt;
inflexp=msc_q;

% FRED, 1960Q1-2020Q1, CBP short-run unemployment rate, NSA
load urcbo.txt

% 1960Q1-2020Q1, unemployment rate, ...
load ur.txt; ur=ur(:,2);

% Unemployment rate gap
gap=ur-urcbo;

% Aggregate inflation expectations w/o component driven by mominal gasoline 
% price shock from SVAR counterfactual, monthly, 1982.7-2020.4, converted
% to quarterly inflation expectations by averaging (consistent with how MSC
% constructs its quarterly data)
load counterfactual; 
for i=1:151
    inflexpnogas(i,1)=mean(counterfactual(i*3-2:i*3));
end;
inflexpnogas=[ones(90,1)*NaN; inflexpnogas];

% 1981Q3-2020Q1, SPF 1 ys inflation forecast
load spf.txt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit Phillips curve 1981Q3-2007Q3 (timing and specification as in CG 2015, 
% panel C in Figure 6, p, 218
y1=infl(87:end,1)-inflexp(87:end,1);
y2=infl(87:end,1)-spf(1:end,1);
    %X1=[ones(length(y1),1) ur(87:end,1) inflexp(87:end,1)]; 
    %X2=[ones(length(y2),1) ur(87:end,1) spf(1:end,1)];
X1=[ones(length(y1),1) gap(87:end,1) inflexp(87:end,1)]; 
X2=[ones(length(y2),1) gap(87:end,1) spf(1:end,1)];
[beta1,~,~,~,~,~,~,R2]=ols(y1(1:end-50,1),X1(1:end-50,1:2)); yhat1=X1*[beta1;1];
[beta2,~,~,~,~,~,~,R2]=ols(y2(1:end-50,1),X2(1:end-50,1:2)); yhat2=X2*[beta2;1];
    %yhatcf1=[ones(155,1) ur(87:end,1) inflexpnogas(87:end,1)]*[beta1;1];
yhatcf1=[ones(155,1) gap(87:end,1) inflexpnogas(87:end,1)]*[beta1;1];
time=1981+3/4:1/4:2020+1/4;
y=infl(87:end,1);
plot(time,y,'k-',time,yhat1,'r:',time,yhatcf1,'k:',time,yhat2,'k--','linewidth',3)
legend('Actual Inflation','PC MSC','PC Counterfactual MSC','PC SPF')
set(gca,'XTick',[2008 2009 2010 2011 2012 2013])
axis([2008 2013+1/4 -10 7])
ylabel('Percent, Annualized')

% Table 1
X=[y-yhat1 y-yhat2 y-yhatcf1 yhat1-yhatcf1] 
X=X(106:127,:);
table=[mean(abs(X(1:5,2))) mean(abs(X(6:13,2)))  mean(abs(X(14:22,2)))  mean(abs(X(6:22,2))) mean(abs(X(1:22,2))); 
       mean(abs(X(1:5,1))) mean(abs(X(6:13,1)))  mean(abs(X(14:22,1)))  mean(abs(X(6:22,1))) mean(abs(X(1:22,1)));
       mean(abs(X(1:5,3))) mean(abs(X(6:13,3)))  mean(abs(X(14:22,3)))  mean(abs(X(6:22,3))) mean(abs(X(1:22,3)))]

 [table(2,1)-table(1,1) table(2,2)-table(1,2) table(2,3)-table(1,3) table(2,4)-table(1,4); 
  table(2,1)-table(3,1) table(2,2)-table(3,2) table(2,3)-table(3,3) table(2,4)-table(3,4)]

