%MAKEPLOTS 
%
%  This file produces some of the plots in 
%
%     The New Keynesian Model is Not Yet Useful for Policy Analysis
%     by V.V. Chari, Patrick Kehoe, and Ellen McGrattan
%
%  The first part of the code is the same as in bca.m
%

%  Ellen McGrattan, 5-19-08

%
% Load dynare matrices dr_. oo_. from:
%
load  usmodel_bca_results.mat


% Vars in X and u, where X = AX(-1) + Bu(-1), u= ea,eb,eg,em,epinf,eqs,ew
%   1. astate
%   2. bstate
%   3. epinfma
%   4. ewma
%   5. gstate
%   6. kp
%   7. kpf
%   8. mstate
%   9. qstate
%  10. r
%  11. spstate
%  12. swstate
%  13. y
%  14. yf
%  15. c
%  16. cf
%  17. inve
%  18. invef
%  19. pinf
%  20. w
%
%  NOTE: to figure out order of variables, look at oo_.SmoothedVariables 
%        a,b,dc,dinve,dw,dy,g,k,kf,labobs,mc,ms,pinfobs,qs,robs,rrf,
%         spinf,sw,wf,zcap,zcapf,astate,bstate,epinfma,ewma,gstate,
%         kp,kpf,mstate,qstate,r,spstate,swstate,y,yf,c,cf,inve,
%         invef,pinf,w,lab,labf,pk,pkf,rk,rkf
%  

A      = dr_.ghx(22:41,:);
B      = dr_.ghu(22:41,:);

% Vars in Y = CX(-1) + Du(-1):
%   c,inve,y,labobs,w,pinfobs,robs
%
C      = dr_.ghx([36,38,34,10,41,13,15],:);
D      = dr_.ghu([36,38,34,10,41,13,15],:);

Ysm    = [oo_.SmoothedVariables.c, ...
          oo_.SmoothedVariables.inve, ...
          oo_.SmoothedVariables.y, ...
          oo_.SmoothedVariables.labobs, ...
          oo_.SmoothedVariables.w, ...
          oo_.SmoothedVariables.pinfobs, ...
          oo_.SmoothedVariables.robs];
Xsm    = [oo_.SmoothedVariables.astate, ...
          oo_.SmoothedVariables.bstate, ...
          oo_.SmoothedVariables.epinfma, ...
          oo_.SmoothedVariables.ewma, ...
          oo_.SmoothedVariables.gstate, ...
          oo_.SmoothedVariables.kp, ...
          oo_.SmoothedVariables.kpf, ...
          oo_.SmoothedVariables.mstate, ...
          oo_.SmoothedVariables.qstate, ...
          oo_.SmoothedVariables.r, ...
          oo_.SmoothedVariables.spstate, ...
          oo_.SmoothedVariables.swstate, ...
          oo_.SmoothedVariables.y, ...
          oo_.SmoothedVariables.yf, ...
          oo_.SmoothedVariables.c, ...
          oo_.SmoothedVariables.cf, ...
          oo_.SmoothedVariables.inve, ...
          oo_.SmoothedVariables.invef, ...
          oo_.SmoothedVariables.pinf, ...
          oo_.SmoothedVariables.w];
usm    = [oo_.SmoothedShocks.ea, ...
          oo_.SmoothedShocks.eb, ...
          oo_.SmoothedShocks.eg, ...
          oo_.SmoothedShocks.em, ...
          oo_.SmoothedShocks.epinf, ...
          oo_.SmoothedShocks.eqs, ...
          oo_.SmoothedShocks.ew];
usma   = 0*usm; usma(:,1)=usm(:,1);
usmb   = 0*usm; usmb(:,2)=usm(:,2);
usmg   = 0*usm; usmg(:,3)=usm(:,3);
usmm   = 0*usm; usmm(:,4)=usm(:,4);
usmp   = 0*usm; usmp(:,5)=usm(:,5);
usmq   = 0*usm; usmq(:,6)=usm(:,6);
usmw   = 0*usm; usmw(:,7)=usm(:,7);

[Y,X]   = dlsim(A,B,C,D,usm,Xsm(1,:)');
      % = Ya+Yb+Yg+Ym+Yq+Yp+Yw-6*Yn;
[Ya,Xa] = dlsim(A,B,C,D,usma,Xsm(1,:)');
[Yb,Xb] = dlsim(A,B,C,D,usmb,Xsm(1,:)');
[Yg,Xg] = dlsim(A,B,C,D,usmg,Xsm(1,:)');
[Ym,Xm] = dlsim(A,B,C,D,usmm,Xsm(1,:)');
[Yp,Xp] = dlsim(A,B,C,D,usmp,Xsm(1,:)');
[Yq,Xq] = dlsim(A,B,C,D,usmq,Xsm(1,:)');
[Yw,Xw] = dlsim(A,B,C,D,usmw,Xsm(1,:)');
[Yn,Xn] = dlsim(A,B,C,D,0*usm,Xsm(1,:)');

%
% Some statistics for the model
%
sdu    = [oo_.posterior_std.shocks_std.ea;
          oo_.posterior_std.shocks_std.eb;
          oo_.posterior_std.shocks_std.eg;
          oo_.posterior_std.shocks_std.em;
          oo_.posterior_std.shocks_std.epinf;
          oo_.posterior_std.shocks_std.eqs;
          oo_.posterior_std.shocks_std.ew];
Sig    = diag(sdu.^2);
VX     = edlyap(A,B*Sig*B');
VY     = C*VX*C'+D*Sig*D';
VX     = diag(VX);

%
% As a check on dynare, compare oo_.SmoothedShocks to u_chk
% and oo_.SmoothedVariables to X_chk
%

u_chk  = 0*usm;
X(1,:) = Xsm(1,:);
X0     = X(1,:)';
for i=1:159;
  u_chk(i,:) = (inv(D)*(Ysm(i+1,:)'-C*X0))';
  X(i+1,:)   = (A*X0+B*u_chk(i,:)')';
  X0         = X(i+1,:)';
end;
X_chk  = X;

%
% For the business cycle accounting, choose an initial state vector
% (Xsm(1,:) is initial smoothed vector from dynare) 
%
X0     = 0*Xsm(1,:)';
u      = 0*usm;
X(1,:) = X0;
for i=1:159;
  u(i,:)     = (inv(D)*(Ysm(i+1,:)'-C*X(i,:)'))';
  X(i+1,:)   = (A*X(i,:)'+B*u(i,:)')';
end;


%
% Load matrices based on one-shock-at-a-time versions of dynare code
%
loadmats

%
% Construct observations for all shocks and subsets of shocks
%

[Ybca,Xbca]   = dlsim(A,B,C,D,u,X0);
            % = Ybcaa+Ybcab+Ybcag+Ybcam+Ybcaq+Ybcap+Ybcaw
[Ybcaa,Xbcaa] = dlsim(Aa,Ba,Ca,Da,u,X0);
[Ybcab,Xbcab] = dlsim(Ab,Bb,Cb,Db,u,X0);
[Ybcag,Xbcag] = dlsim(Ag,Bg,Cg,Dg,u,X0);
[Ybcam,Xbcam] = dlsim(Am,Bm,Cm,Dm,u,X0);
[Ybcap,Xbcap] = dlsim(Ap,Bp,Cp,Dp,u,X0);
[Ybcaq,Xbcaq] = dlsim(Aq,Bq,Cq,Dq,u,X0);
[Ybcaw,Xbcaw] = dlsim(Aw,Bw,Cw,Dw,u,X0);

t=ptime(1965.2,length(oo_.SmoothedVariables.c),4);

load actual.out
tact = actual(:,1);
Yact = actual(:,2:8);

figure(1)
plot(tact,Yact(:,3), ...
     t,Ybcab(:,3)+Ybcag(:,3)+Ybcap(:,3)+Ybcaw(:,3))
     %t,Ybcaa(:,3)+Ybcam(:,3)+Ybcaq(:,3))
set(gca,'Fontsize',14)
title('Movements in Output','Fontsize',14)
legend('Actual','Predicted due to b,g,p,w','Location','Southeast')
v=axis;
axis([1965,2005,v(3),v(4)])

figure(2)
plot(tact,Yact(:,4), ...
     t,Ybcab(:,4)+Ybcag(:,4)+Ybcap(:,4)+Ybcaw(:,4))
     %t,Ybcaa(:,4)+Ybcam(:,4)+Ybcaq(:,4))
set(gca,'Fontsize',14)
title('Movements in Hours','Fontsize',14)
legend('Actual','Predicted due to b,g,p,w','Location','Southeast')
v=axis;
axis([1965,2005,v(3),v(4)])

figure(3)
plot(tact,Yact(:,6), ...
     t,Ybcab(:,6)+Ybcag(:,6)+Ybcap(:,6)+Ybcaw(:,6))
     %t,Ybcaa(:,6)+Ybcam(:,6)+Ybcaq(:,6))
set(gca,'Fontsize',14)
title('Movements in Inflation','Fontsize',14)
legend('Actual','Predicted due to b,g,p,w','Location','Northeast')
v=axis;
axis([1965,2005,v(3),v(4)])


fig2act  = [tact,Yact(:,[3,4,6])];
fig2pred = [t,Ybcab(:,3)+Ybcag(:,3)+Ybcap(:,3)+Ybcaw(:,3), ...
              Ybcab(:,4)+Ybcag(:,4)+Ybcap(:,4)+Ybcaw(:,4), ...
              Ybcab(:,6)+Ybcag(:,6)+Ybcap(:,6)+Ybcaw(:,6)];

