function HistDecompStructure = DSGEHistDecompFcn(Y,X,A,H,F,B0,Ksitt,KsitT,etatT,wtT,h)
% DSGEHistDecompFcn: Calculates historical decompositions of the h-steps ahead forecast errors.
%
% USAGE:
%
%       HistDecompStructure = DSGEHistDecompFcn(Y,X,A,H,F,B0,Ksitt,KsitT,etatT,wtT,h)
%
% REQUIRED INPUT:  Y (nxT matrix), with the observed variables.
%
%                  X (kxT matrix), with the exogenous variables.
%
%                  A (kxn matrix), with the coefficients on X in the measurement equation of the state-space
%                                  representation of the log-linearized DSGE model.
%
%                  H (rxnxT matrix), with the coefficients on the state variables Ksi in the measurement equation.
%
%                  F (rxr matrix), with the state transition parameters in the state equation of the state-space
%                                  representation of the log-linearized DSGE model.
%
%                  B0 (rxq matrix), with the parameters on the economic shocks, eta, in the state equation.
%
%                  Ksitt (rxT matrix), with the filtered state variables.
%
%                  KsitT (rxT matrix), with the smoothed state variables.
%
%                  etatT (qxT matrix), with the smoothed estimates of the economic shocks in the state equation.
%
%                  wtT (nxT matrix), with the smoothed estimates of the non-zero measurement errors in the
%                                   measurement equation. This is an empty matrix if there are no measurement
%                                   errors in the model.
%
%                  h (positive integer), with the forecast horizon.
%
% REQUIRED OUTPUT: HistDecompStructure (structure), with the fields "epstth", an nx(T-h) matrix with
%                                   h-steps ahead forecast errors, "wthT", an nx(T-h) or empty matrix
%                                   with the smoothed measurement errors, "KsiErr", an nx(T-h) matrix
%                                   with the state projection error influence, "etathT", and nx(T-h)Xq
%                                   matrix with structural error terms, and "KeepVar" a vector with
%                                   indices that gives the non-zero columns of B0.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2009 European Central Bank.
%
%                       First version: November 10, 2006.
%                        This version: January 22, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 15-11-2006: Finished the function now that StateSmoother has been fixed.
%
% * 16-11-2006: Renamed the function.
%
% * 13-11-2007: Updated the documentation.
%
% * 23-05-2008: Updated the documentation.
%
% * 09-07-2008: Added code to allow for time-varing measurement matrix H.
%
% * 22-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

HistDecompStructure = [];
%
% determine some dimensions
%
[n,T] = size(Y);
k = min(size(X));
[B,KeepVar] = RemoveRedundantColumns(B0);
B0 = B;
[r,q] = size(B0);
HistDecompStructure.KeepVar = KeepVar;

%
% Start by calculating the h-steps ahead forecasts of the state variables
%
Fh = F^h;
Ksitth = Fh*Ksitt(:,1:T-h);
%
% Now we compute the forecast errors
%
if length(size(H))==2;
   if k>0;
      HistDecompStructure.epstth = Y(:,h+1:T)-(A'*X(:,h+1:T))-(H'*Ksitth);
   else;
      HistDecompStructure.epstth = Y(:,h+1:T)-(H'*Ksitth);
   end;
   %
   % now for the measurement errors
   %
   if isempty(wtT)==0;
      HistDecompStructure.wthT = wtT(:,h+1:T);
   else;
      HistDecompStructure.wthT = [];
   end;
   %
   % now for the state variables
   %
   HistDecompStructure.KsiErr = H'*Fh*(KsitT(:,1:T-h)-Ksitt(:,1:T-h));
   HistDecompStructure.etathT = zeros(n,T-h,q);
   for j=1:q;
      for i=0:h-1;
         %
         % add the economic errors of shock j for horizon i
         %
         HistDecompStructure.etathT(:,:,j) = HistDecompStructure.etathT(:,:,j)+(H'*(F^i)*B0(:,j)*etatT(j,h+1-i:T-i));
      end;
   end;
else;
   if k>0;
      HistDecompStructure.epstth = zeros(n,T-h);
      for t=h+1:T;
         HistDecompStructure.epstth(:,t-h) = Y(:,t)-(A'*X(:,t))-(H(:,:,t)'*Ksitth(:,t-h));
      end;
   else;
      for t=h+1:T;
         HistDecompStructure.epstth(:,t-h) = Y(:,t)-(H(:,:,t)'*Ksitth(:,t-h));
      end;
   end;
   %
   % now for the measurement errors
   %
   if isempty(wtT)==0;
      HistDecompStructure.wthT = wtT(:,h+1:T);
   else;
      HistDecompStructure.wthT = [];
   end;
   %
   % now for the state variables
   %
   HistDecompStructure.KsiErr = zeros(n,T-h);
   for t=h+1:T;
      HistDecompStructure.KsiErr(:,t-h) = H(:,:,t)'*Fh*(KsitT(:,t-h)-Ksitt(:,t-h));
   end;
   HistDecompStructure.etathT = zeros(n,T-h,q);
   for j=1:q;
      for i=0:h-1;
         %
         % add the economic errors of shock j for horizon i
         %
         for t=h+1:T;
            HistDecompStructure.etathT(:,t-h,j) = HistDecompStructure.etathT(:,t-h,j)+(H(:,:,t)'*(F^i)*B0(:,j)*etatT(j,t-i));
         end;
      end;
   end;
end;

%
% end of DSGEHistDecompFcn.m
%
