%VDEC
%  This file computes forecast error variance decompositions for the model
%  in Smets-Wouters (AER 2007)
%

%  Ellen McGrattan, 5-20-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],:);

sdu    = [oo_.posterior_mode.shocks_std.ea;
          oo_.posterior_mode.shocks_std.eb;
          oo_.posterior_mode.shocks_std.eg;
          oo_.posterior_mode.shocks_std.em;
          oo_.posterior_mode.shocks_std.epinf;
          oo_.posterior_mode.shocks_std.eqs;
          oo_.posterior_mode.shocks_std.ew];
Sig    = diag(sdu.^2);
sSig   = diag(sdu);


%
% Calculate the variance decomposition of the j-step ahead prediction errors:
%
%      v(j) = E[ Z(t+j)-E[Z(t+j)| X(t) ] [ Z(t+j)-E[Z(t+j)| X(t) ]'
%
% j in J, for the system
%
%         X(t+1) = A X(t) + B u(t)
%          Z(t)  = C X(t) + D u(t)    (note: Z(t)=Y(t+1) above)
%          E u(t) u(t)' = Sig
%
% For this system:
%
%      v(0) = D*Sig*D'
%      v(1) = D*Sig*D' + C*F*F'*C' 
%      v(2) = D*Sig*D' + C*F*F'*C' + C*A*F*F'*C'
%       :
%      v(j) = v(j-1) + C*A^(j-1)*F*F'*A^(j-1)'*C'
%
% where F=B*Sig^(1/2).  To calculate the j-step ahead variance due to the 
% kth component of u(t), k in K, we compute
%
%      vk(j) = E[ Z(t+j)-E[Z(t+j) | u(k,t),u(k,t-1),...]]*
%                [Z(t+j)-E[Z(t+j) | u(k,t),u(k,t-1),...]]'
%
%            = vk(j-1) + C*A^(j-1)*F*i(k)*F'*A^(j-1)'*C'
%
% where vk(0) = D*Sig^(1/2)*i(k)*Sig^(1/2)*D' and i(k) is a matrix
% of zeros with element (k,k) equal to 1.

J     = input('What terminal horizon for forecast error decomposition? ');
i1    = zeros(7); i1(1,1)=1;
i2    = zeros(7); i2(2,2)=1;
i3    = zeros(7); i3(3,3)=1;
i4    = zeros(7); i4(4,4)=1;
i5    = zeros(7); i5(5,5)=1;
i6    = zeros(7); i6(6,6)=1;
i7    = zeros(7); i7(7,7)=1;
F     = B*sSig;
vj    = D*Sig*D';
vj1   = D*sSig*i1*sSig*D';
vj2   = D*sSig*i2*sSig*D';
vj3   = D*sSig*i3*sSig*D';
vj4   = D*sSig*i4*sSig*D';
vj5   = D*sSig*i5*sSig*D';
vj6   = D*sSig*i6*sSig*D';
vj7   = D*sSig*i7*sSig*D';

Vj    = zeros(J,8);
Vj1   = zeros(J,8);
Vj2   = zeros(J,8);
Vj3   = zeros(J,8);
Vj4   = zeros(J,8);
Vj5   = zeros(J,8);
Vj6   = zeros(J,8);
Vj7   = zeros(J,8);
for j=1:J;
  AI       = A^(j-1);
  vj       = vj   + C*AI*F*F'*AI'*C';
  vj1      = vj1  + C*AI*F*i1*F'*AI'*C';
  vj2      = vj2  + C*AI*F*i2*F'*AI'*C';
  vj3      = vj3  + C*AI*F*i3*F'*AI'*C';
  vj4      = vj4  + C*AI*F*i4*F'*AI'*C';
  vj5      = vj5  + C*AI*F*i5*F'*AI'*C';
  vj6      = vj6  + C*AI*F*i6*F'*AI'*C';
  vj7      = vj7  + C*AI*F*i7*F'*AI'*C';
  Vj(j,:)  = [j,diag(vj)'];
  Vj1(j,:) = [j,diag(vj1)'];
  Vj2(j,:) = [j,diag(vj2)'];
  Vj3(j,:) = [j,diag(vj3)'];
  Vj4(j,:) = [j,diag(vj4)'];
  Vj5(j,:) = [j,diag(vj5)'];
  Vj6(j,:) = [j,diag(vj6)'];
  Vj7(j,:) = [j,diag(vj7)'];
end
tab        = [Vj1(J,2:8)./Vj(J,2:8)
              Vj2(J,2:8)./Vj(J,2:8)
              Vj3(J,2:8)./Vj(J,2:8)
              Vj4(J,2:8)./Vj(J,2:8)
              Vj5(J,2:8)./Vj(J,2:8)
              Vj6(J,2:8)./Vj(J,2:8)
              Vj7(J,2:8)./Vj(J,2:8)]'*100;
tab        = [tab,sum(tab(:,[2,3,5,7])')'];

disp(sprintf('Forecast Error Decomposition at Horizon %g',J))
disp('   Observables\Shocks:  a     b     g     m     p     q     w  b+g+p+w')
disp(sprintf( ...
     '     Consumption     %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(1,:)))
disp(sprintf( ...
     '     Investment      %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(2,:)))
disp(sprintf( ...
     '     Output          %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(3,:)))
disp(sprintf( ...
     '     Hours           %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(4,:)))
disp(sprintf( ...
     '     Wages           %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(5,:)))
disp(sprintf( ...
     '     Inflation       %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(6,:)))
disp(sprintf( ...
     '     Interest Rate   %6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f%6.2f',tab(7,:)))
disp(' ')
disp('   Note: Decompositions are reported for variables in levels ')
disp('         Smets-Wouters report results for differenced variables')
