function IRStructure = DSGEImpulseResponseFcn(H,F,B0,h)
% DSGEImpulseResponseFcn: Calculates the impulse response functions for Y, the observed variables, and
%                         Xi, the state variables, from period 0 to h for all possible shocks
%
% USAGE:
%
%       IRStructure = DSGEImpulseResponseFcn(H,F,B0,h)
%
% REQUIRED INPUT: H (rxn) matrix where n is the number of observed variables and r the number of
%                         state variables. This matrix is part of the measurement equation of the
%                         state-space representation for the Kalman filter.
%
%                 F (rxr) matrix from the state equation.
%
%                 B0 (rxq) matrix, where Bo*B0' is the covariance matrix for the shocks to the
%                         state equation.
%
%                 h (positive integer) giving the horizon for the impulse responses.
%
% REQUIRED OUTPUT: IRStructure (structure) with fields Ksi (rxqxh+1) and Y (nxqxh+1). The responses in the
%                         invidual variables are located in the rows, while the columns give the shocks.
%                         The third dimension is the response horizon.
%
%
%                       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: August 14, 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-08-2006: Fixed a bug using Fprevious instead of FPrevious in loop!
%
% * 31-10-2006: Made sure that all small values, i.e., less than 1.0e-12, are replaced by 0.
%
% * 26-02-2007: Changed the output into a structure with fields "Ksi" and "Y", both being 3D matrices.
%
% * 23-08-2007: New and faster method for setting small B0 elements to zero.
%
% * 13-11-2007: Updated the documentation.
%
% * 23-05-2008: Updated the documentation.
%
% * 22-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% Fix all B0 elements
%
[r,q] = size(B0);
n = size(H,2);
IRStructure.Ksi = zeros(r,q,h+1);
IRStructure.Y = zeros(n,q,h+1);
%
% set small responses to zero
%
Bpos = (B0>1.0e-12);
Bneg = (B0<1.0e-12);
B0 = (B0 .* Bpos)+(B0 .* Bneg);
%
% initial responses
%
IRStructure.Ksi(:,:,1) = B0;
IRStructure.Y(:,:,1) = H'*B0;
FPrevious = eye(r);
for i=1:h;
   FPrevious = F*FPrevious;
   NewIR = FPrevious*B0;
   %
   % set small responses to zero
   %
   for i1=1:r;
      for j1=1:q;
         if abs(NewIR(i1,j1))<10.e-12;
            NewIR(i1,j1) = 0;
         end;
      end; 
   end;
   IRStructure.Ksi(:,:,i+1) = NewIR;
   IRStructure.Y(:,:,i+1) = H'*NewIR;
end;

%
% end of DSGEImpulseResponseFcn.m
%
