function [FEVDs,LRVD,status,RiccatiResults,UniquenessCheck] = DSGEVarianceDecompLevelsFcn(H,F,B0,R,AccMat,h,DAMaxIter,DAConvValue,RicMaxIter,RicConvValue)
% DSGEVarianceDecompLevelsFcn: Calculates the forecast error variance decomposition given that the matlab function dare
%                        exists.
%
% USAGE:
%
%       [FEVDs,LRVD,status,RiccatiResults,UniquenessCheck] = DSGEVarianceDecompLevelsFcn(H,F,B0,R,AccMat,h,DAMaxIter,DAConvValue,RicMaxIter,RicConvValue)
%
% 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 B0*B0' is the covariance matrix for the shocks to the
%                          state equation.
%
%                  R (nxn) matrix with the covariance for the measurement errors. May also be the scalar
%                          0.
%
%                  AccMat (nxn) diagonal 0-1 matrix. Only those diagonal elements that that 1
%                          will be accumulated in the
%
%                  h (positive integer) giving the forecast horizon.
%
%                  DAMaxIter (integer) which determines the maximum number of iterations that
%                          the doubling algorithm for the Lyapunov aquations can use.
%
%                  DAConvValue (integer) which determined the tolerance level for the convergence of the
%                          doubling algorithms for the Lyapunov equations.
%
%                  RicMaxIter (integer) which determined the maximum number of iterations that the
%                          Riccati equation solver can use.
%
%                  RicConvValue (integer) which determines the tolerance level for the convergence of
%                          the Riccati equations solver.
%
% REQUIRED OUTPUT: FEVDs (nx(n+q+1)xh) matrix with the forecast error variance decompositions. For each h, the nx(n+q+1) matrix
%                          is ordered such the the measurement errors are located in the first n columns, and the q economic
%                          shocks in the following columns. The remaining column captures possible signal extraction error
%                          variances. The observed variables are determined through the row numbers.
%
%                  LRVD (nx(n+q)) matrix with the long run forecast error variance decomposition.
%
%                  status (integer) takes the value 0 if the solution method converged and 1 otherwise. It takes
%                          on the value 2 in the special case that infinite or NaN's are found in P1 during the
%                          iterations.
%
%                  RiccatiResults (1x2 vector), with the number of iterations and the TestValue for the Riccati
%                          solver in the columns. The TestValue refers to the degree of convergence.
%
%                  UniquenessCheck (scalar), with largest eigenvalues of (F-K*H'), where K is the asymptote of the
%                          Kalman gain matrix, i.e., K = F*P1*H*inv(H'*P1*H+R) for the full forecast error
%                          variance.
%
%
%                       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 Europea Central Bank.
%
%                       First version: April 21, 2008. 
%                        This version: January 23, 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:
%
% * 22-04-2008: Finished the function.
%
% * 23-04-2008: If P0 with P1 set to Q is sufficiently close to zero, then P1 = Q is used.
%
% * 23-05-2008: Updated the documentation.
%
% * 23-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

FEVDs = [];
LRVD = [];
status = 0;
RiccatiResults = [1 0];
UniquenessCheck = 1;
%
% Use the doubling algorithm to solve for an initial value for P1, given by the unconditional
% covariance matrix of the state variables.
%
[r n] = size(H);
q = size(B0,2);
Q = B0*B0';
P0 = Q - (Q*H*inv((H'*Q*H)+R)*H'*Q);
[P_init,stat] = DoublingAlgorithmLyapunov(F,Q,DAMaxIter*100,10^(-(DAConvValue+1)));
%
if max(max(abs(P0)))>1.0e-10;
   [P1,status,NumIter,TestValue] = RiccatiSolver(F,H,Q,R,P_init,RicMaxIter*1000,10^(-(RicConvValue+1)));
   RiccatiResults = [NumIter TestValue];
else;
   P1 = Q;
end;
if status==0;
   K = F*P1*H*inv((H'*P1*H)+R);
   UniquenessCheck = max(abs(eig(F-K*H')));
   %
   % initialize the matrix with forecast error variance decompositions.
   %
   FEVDs = zeros(n,q+n+1,h);
   LRVD = zeros(n,q+n);
   %
   % compute the P1,...,Ph matrices collected into PMatrix
   %
   PMatrix = zeros(r,r,max(h,200));
   PMatrix(:,:,1) = P1;
   for i=2:max(h,200);
      PMatrix(:,:,i) = (F*PMatrix(:,:,i-1)*F')+Q;
   end;
   %
   % First compute the variances
   %
   FEVars = zeros(n,n,h);
   FEVars(:,:,1) = (H'*P1*H)+R;
   LRVars = H'*B0;
   FjProd = zeros(r,r);
   Fj = F;
   for i=2:h;
      FjProd = F*(PMatrix(:,:,i-1)+FjProd);
      FEVars(:,:,i) = (H'*PMatrix(:,:,i)*H)+R+(AccMat*(FEVars(:,:,i-1)+(H'*(FjProd+FjProd')*H))*AccMat);
      LRVars = LRVars+(H'*Fj*B0);
      Fj = Fj*F;
   end;
   %
   % check if we need to add additional terms to LRVars
   %
   if h<200;
      for i=h+1:200;
         LRVars = LRVars+(H'*Fj*B0);
         Fj = Fj*F;
      end;
   end;
   LRVars = (AccMat*diag(diag((LRVars*LRVars')+R)))+((eye(n)-AccMat)*diag(diag((H'*P_init*H)+R)));
   %
   % Start with the measurement errors. Now we apply B0 = 0, so that Ph = 0 for all h
   %
   for i=1:h;
      FEVDs(:,1:n,i) = inv(diag(diag(FEVars(:,:,i))))*diag(diag((i*AccMat*R)+((eye(n)-AccMat)*R)));
   end;
   LRVD(:,1:n) = inv(LRVars)*diag(diag(R));
   %
   % Now we turn to the shocks
   %
   for j=1:q;
      B = zeros(size(B0));
      B(:,j) = B0(:,j);
      Q = B*B';
      P0 = Q - (Q*H*inv(H'*Q*H)*H'*Q);
      [P_init,stat] = DoublingAlgorithmLyapunov(F,Q,DAMaxIter*100,10^(-(DAConvValue+1)));
      %
      if max(max(abs(P0)))>1.0e-10;
         [P1,status,NumIter,TestValue] = RiccatiSolver(F,H,Q,0,P_init,RicMaxIter*1000,10^(-(RicConvValue+1)));
      else;
         P1 = Q;
      end;
      if status==0;
         %
         % compute the P1,...,Ph matrices collected into PMatrix
         %
         PMatrix = zeros(r,r,max(h,200));
         PMatrix(:,:,1) = P1;
         for i=2:max(h,200);
            PMatrix(:,:,i) = (F*PMatrix(:,:,i-1)*F')+Q;
         end;
         %
         % First compute the variances
         %
         FEVarj = zeros(n,n,h);
         FEVarj(:,:,1) = H'*P1*H;
         FEVDs(:,j+n,1) = diag(FEVarj(:,:,1))./diag(FEVars(:,:,1));
         LRVarj = H'*B;
         FjProd = zeros(r,r);
         Fj = F;
         for i=2:h;
            FjProd = F*(PMatrix(:,:,i-1)+FjProd);
            FEVarj(:,:,i) = (H'*PMatrix(:,:,i)*H)+(AccMat*(FEVarj(:,:,i-1)+(H'*(FjProd+FjProd')*H))*AccMat);
            FEVDs(:,j+n,i) = diag(FEVarj(:,:,i))./diag(FEVars(:,:,i));
            LRVarj = LRVarj+(H'*Fj*B);
            Fj = Fj*F;
         end;
         %
         % check if we need to add additional terms to LRVars
         %
         if h<200;
            for i=h+1:200;
               LRVarj = LRVarj+(H'*Fj*B);
               Fj = Fj*F;
            end;
         end;
         LRVarj = (AccMat*diag(diag(LRVarj*LRVarj')))+((eye(n)-AccMat)*diag(diag(H'*P_init*H)));
         LRVD(:,n+j) = inv(LRVars)*diag(LRVarj);
      else;
         FEVDs = [];
         LRVD = [];
         break;
      end;
   end;
   %
   % now compute the signal extraction error quickly, i.e. as the remainder.
   %
   if status==0;
      for i=1:h;
         FEVDs(:,q+n+1,i) = ones(n,1)-sum(FEVDs(:,:,i),2);
      end;
   end;
end;

%
% end of DSGEVarianceDecompLevelsFcn.m
%
