%
%

figuresdir = '../20070130favarpostfomc' ;

ipind = 1 ;
unempind = 16 ;
cpiind = 1 ;
cpixind = 2 ;
pceind = 3 ;
pcexind = 4 ;
poilind = 19 ;
ffrind = 1 ;
dollarind = 10 ;

HPloose = 5e5 ;
HPtight = 5e7 ;
HPlooselt = 3e7 ;
HPtightlt = 3e8 ;

if strcmp(diffspecification,'no') ;
% 1. Output gap and Unemployment gap
figure(1) ;
plot(date(1:T),data(burn+1:T+burn,unempind),'LineWidth',1) ;
hold on ;
axis([1965 2010 3 11]) ;
shaderec ;
h1 = plot(date(1:T),hpfit(data(burn+1:T+burn,unempind),HPtight),'k--','LineWidth',1) ;
h2 = plot(date(1:T),hpfit(data(burn+1:T+burn,unempind),HPtightlt),'r:','LineWidth',1) ;
h3 = plot(date(1:T),data(burn+1:T+burn,unempind),'LineWidth',1) ;
set(gca,'FontSize',12) ;
hold off ;
legend([h3,h1,h2],'unemployment rate','shorter-term trend','longer-term trend') ;
ylabel('percent') ;
print('-dpng',strcat(figuresdir,'/figure1a.png')) ;

figure(2) ;
plot(date(1:T),gdpm(1:T),'LineWidth',1) ;
hold on ;
axis([1965 2010 3000 12000]) ;
shaderec ;
h1 = plot(date(1:T),exp(trendgdpm(1:T)),'k--','LineWidth',1) ;
h2 = plot(date(1:T-4),exp(trendgdpltm(1:T-4)),'r:','LineWidth',1) ;
h3 = plot(date(1:T),gdpm(1:T),'LineWidth',1) ;
set(gca,'FontSize',12) ;
hold off ;
legend([h3,h1,h2],'GDP','shorter-term trend','longer-term trend','Location','SouthEast') ;
ylabel('Billions of Chained 2000 $') ;
print('-dpng',strcat(figuresdir,'/figure1b.png')) ;

figure(3) ;
unempgap = datastd(unempind) * eafactors * ealoadings(:,unempind) ;
plot(date(1:T),100*gdpgapm(1:T),date(1:T),unempgap(1:T),'--r') ;
axis([1965 2010 -4 4]) ;
set(gca,'FontSize',12) ;
legend('GDP gap','unemployment gap','Location','SouthEast') ;
ylabel('percent') ;
print('-dpng',strcat(figuresdir,'/figure1c.png')) ;


% 1b. Okun's Law
okun1 = 100 * std(gdpgapm(1:T)) /std(unempgap(1:T)) ;

Dunempgap = diff(unempgap(1:T)) ;
Dgdpgapm = diff(gdpgapm(1:T)) ;
D4unempgap = unempgap(5:T) - unempgap(1:T-4) ;
D4gdpgapm = gdpgapm(5:T) - gdpgapm(1:T-4) ;

D12gdp = lgdpm(13:T) - lgdpm(1:T-12) ;
D12unemp = data(burn+13:burn+T,unempind) - data(burn+1:burn+T-12,unempind) ;

okun1a = ols(100*D12gdp,[ones(T-12,1),D12unemp]) ;
okun2 = ols(100*gdpgapm(1:T),[ones(size(gdpgapm(1:T))),unempgap(1:T)]) ;
okun3 = ols(100*Dgdpgapm,[ones(size(Dgdpgapm)),Dunempgap]) ;
okun4 = ols(100*D4gdpgapm,[ones(size(D4gdpgapm)),D4unempgap]) ;
okun5 = ols(100*lag(Dgdpgapm),[ones(size(Dgdpgapm)),Dunempgap]) ;


% 2. Inflation factor "Phillips Curve" coefficients:
coeffmat = reshape(B(1:54,3),9,6) ;
recursedata = [prfactors(:,1),trafactors(:,2),eafactors(:,1),eafactors(:,2),...
               prfactors(:,[3 2]),intfactors(:,1),intfactors(:,2),trafactors(:,1)] ; 
figure(4) ;
irfs = VARrecirf(recursedata,6,[],50) ;
shocksd = diag(squeeze(irfs(1,:,:))) ;

nperiods = size(irfs,1) ;
figure(5) ; % CPI impulse responses
 plot([0:nperiods-1],irfs(:,1,1) *prloadings(1,cpiind) *datastd(72+cpiind) *100 ...
                    /(shocksd(1) *prloadings(1,cpiind) *datastd(72+cpiind) *100)) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,0,1.1]) ;
 print('-dpng',strcat(figuresdir,'/figure2a.png')) ;
figure(6) ;
 plot([0:nperiods-1],irfs(:,1,3) *prloadings(1,cpiind) *datastd(72+cpiind) *100 ...
                    /(shocksd(3) *ealoadings(1,unempind) *datastd(unempind))) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,-.5,.1]) ;
 print('-dpng',strcat(figuresdir,'/figure2b.png')) ;
figure(7) ;
 plot([0:nperiods-1],irfs(:,1,6) *prloadings(1,cpiind) *datastd(72+cpiind) *100 ...
                    /(shocksd(6) *prloadings(2,poilind) *datastd(72+poilind) *100)) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,-.01,.012]) ;
 print('-dpng',strcat(figuresdir,'/figure2c.png')) ;
figure(8) ;
 plot([0:nperiods-1],irfs(:,1,9) *prloadings(1,cpiind) *datastd(72+cpiind) *100 ...
                    /(shocksd(9) *traloadings(1,dollarind) *datastd(107+dollarind))) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,-.01,.01]) ;
 print('-dpng',strcat(figuresdir,'/figure2d.png')) ;

figure(5) ; % unemployment rate impulse responses
 plot([0:nperiods-1],irfs(:,3,1) *ealoadings(1,unempind) *datastd(unempind) ...
                    /(shocksd(1) *prloadings(1,cpiind) *datastd(72+cpiind) *100)) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
axis([0,50,-.1,.4]) ;
 print('-dpng',strcat(figuresdir,'/figure2e.png')) ;
figure(6) ;
 plot([0:nperiods-1],irfs(:,3,3) *ealoadings(1,unempind) *datastd(unempind) ...
                    /(shocksd(3) *ealoadings(1,unempind) *datastd(unempind))) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,-.6,1.4]) ;
 print('-dpng',strcat(figuresdir,'/figure2f.png')) ;
figure(7) ;
 plot([0:nperiods-1],irfs(:,3,6) *ealoadings(1,unempind) *datastd(unempind)  ...
                    /(shocksd(6) *prloadings(2,poilind) *datastd(72+poilind) *100)) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,-.01,.01]) ;
 print('-dpng',strcat(figuresdir,'/figure2g.png')) ;
figure(8) ;
 plot([0:nperiods-1],irfs(:,3,9) *ealoadings(1,unempind) *datastd(unempind)  ...
                    /(shocksd(9) *traloadings(1,dollarind) *datastd(107+dollarind))) ;
 hold on; plot([1:nperiods],zeros(nperiods),':k'); hold off ;
 axis([0,50,-.015,.01]) ;
 print('-dpng',strcat(figuresdir,'/figure2h.png')) ;

end ; % if difference specification

% 4. Inflation forecast breakdown:
% Long-run average inflation rates:
lrcpi = 100 * ((1+robustmean(diff(log(prfinal(1:T,cpiind)))))^12 -1) ;
lrcpix = 100 * ((1+robustmean(diff(log(prfinal(1:T,cpixind)))))^12 -1) ;
lrpce = 100 * ((1+robustmean(diff(log(prfinal(1:T,pceind)))))^12 -1) ;
lrpcex = 100 * ((1+robustmean(diff(log(prfinal(1:T,pcexind)))))^12 -1) ;
lrpcex85 = 100 * ((1+robustmean(diff(log(prfinal(237:T,pcexind)))))^12 -1) ;

%
% Plot some counterfactuals for the current forecast, setting all factors to their
%  historical average and then setting each latent factor in turn equal to
%  its recent values:
%

figure(9) ;
piyoy = 100*[repmat(NaN,12,1);(prfinal(13:end,pcexind)-prfinal(1:end-12,pcexind))./prfinal(1:end-12,pcexind)] ;
pi3m =  100*((1+[repmat(NaN,3,1);(prfinal(4:end,pcexind)-prfinal(1:end-3,pcexind))./prfinal(1:end-3,pcexind)]).^4 -1) ;
T1 = lastobs(prindexes(pcexind)) ;
h1 = plot(date(1:T1),pi3m(1:T1,1),'LineWidth',1) ;
h = patch([2000,2000,2010,2010]',[1,2,2,1]',[.9,.9,.9]) ;
set(h,'EdgeColor','none') ;
hold on ;
h2 = plot(date(1:T1),piyoy(1:T1),'r','LineWidth',1) ;
plot(date,pi3m(:,1),'--','LineWidth',1) ;
plot(date,piyoy,'r--','LineWidth',1) ;
ylabel('percent') ;
axis([2005,2010,0,4]) ;
shaderec(date(T)-.03/12, date(T)+.03/12, [.3,.3,.3]) ;
set(gca,'FontSize',12) ;
hold off ;
legend([h1,h2],'3-month change','12-month change','Location','SouthEast') ;
print('-dpng',strcat(figuresdir,'/figure3.png')) ;



n = size(vardata,2) ;
for i = 1:n ;
  vardata1 = repmat(mean(vardata(1:T-1,:)),varlags,1) ;
  vardata1(1:varlags,i) = vardata(T-varlags+1:T,i) ;
  cfactfcast(:,:,i) = [zeros(1,n); VARfcast(vardata1,varlags,fcasthorizon+1,B)] ;
end ;
figure(10) ;
plot(date,zeros(T+fcasthorizon,1),'k','LineWidth',1) ; % Baseline forecast
axis([2006,2010,-1,.5]) ;
hold on ;
plot([1964 2011],zeros(2,1),':k') ;
shaderec(date(T)-.02/12, date(T)+.02/12, [.3,.3,.3]) ;
set(gca,'FontSize',12) ;
h1 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,1)),'--r','LineWidth',1) ; % better
h2 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,2)),'-.','Color',[.9,.4,.4],'LineWidth',1) ; % much worse
h3 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,3)),'--','Color',[1 .6 0],'LineWidth',1) ;
h4 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,4)),'-.y','LineWidth',1) ; % better
h5 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,5)),'--g','LineWidth',1) ; % better
h6 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,6)),'-.','Color',[0,.6,0],'LineWidth',1) ; % better
h7 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,7)),'--b','LineWidth',1) ; % better
h8 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,8)),'-.','Color',[.67,0,1],'LineWidth',1) ; % much worse
h9 = plot(date(T-1:T+fcasthorizon),squeeze(cfactfcast(:,3,9)),'--','Color',[.4,.4,.4],'LineWidth',1) ; % much worse
hold off ;
legend([h1,h2,h3,h4,h5,h6,h7,h8,h9],...
    'output','leading ind','inflation','oil price','commod pr',...
    'level','slope','exch rate','trade bal','Location','SouthWest') ;

print('-dpng',strcat(figuresdir,'/figure4.png')) ;

  


