function [FEVDs,LRVD,status,RiccatiResults,UniquenessCheck] = DSGEVarianceDecompFcn(H,F,B0,R,h,DAMaxIter,DAConvValue,RicMaxIter,RicConvValue)
% DSGEVarianceDecompFcn: Calculates the forecast error variance decomposition given that the matlab function dare
%                        exists.
%
% USAGE:
%
%       [FEVDs,LRVD,status,RiccatiResults,UniquenessCheck] = DSGEVarianceDecompFcn(H,F,B0,R,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.
%
%                  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 European Central Bank.
%
%                       First version: October 23, 2006.
%                        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:
%
% * 25-10-2006: Changes the convergence value for the Riccati equation solver from 1.0e-8 to 1.0e-5.
%
% * 31-10-2006: Made sure that the value of DSGEModel.DAToleranceValue is sent in the appropriate form to the
%               "DoublingAlgorithmLyapunov" function. Furthermore, the function "RiccatiSolver" now takes the
%               model selected number of iterations and tolerance level into account. Long run variance
%               decompositions were also added and provided as output. The input CurrINI has been removed in
%               the call to the function. Furthermore, the DSGEModel structure is no longer given as input.
%               Instead, the function takes the 4 field values of this structure that it previously required
%               as the 4 variables "DAMaxIter", "DAConvValue", "RicMaxIter", and
%               "RicConvValue".
%
% * 20-11-2006: The status from the Riccati solver function is given as output.
%
% * 21-11-2006: The function also checks for the uniqueness of P1. "UniquenessCheck" variable therefore added
%               to the output list.
%
% * 22-11-2006: Added one column to the FEVDs output matrix and one row to RiccatiResults.
%
% * 24-11-2006: Made sure the last run of the Riccati solver is only considered if all previous runs were
%               successful. Moreover, if this run is unsuccessful, then FEVDs and LRVD are set to empty
%               matrices.
%
% * 13-11-2007: Updated the documentation.
%
% * 21-04-2008: Updated the documnetation.
%
% * 22-04-2008: Using B*B' = P1 when the measurement errors are set to zero. UniquenessCheck is now a scalar, while
%               RiccatiResults is a vector.
%
% * 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;
%
[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);
   LRVars = diag(diag((H'*P_init*H)+R));
   %
   % First compute the variances
   %
   FEVars = zeros(n,h);
   FEVars(:,1) = diag((H'*P1*H)+R);
   Pi = P1;
   for i=2:h;
      Pi = (F*Pi*F')+(B0*B0');
      FEVars(:,i) = diag((H'*Pi*H)+R);
   end;
   %
   % 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(FEVars(:,i)))*diag(diag(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';
      [P_init,stat] = DoublingAlgorithmLyapunov(F,Q,DAMaxIter*100,10^(-(DAConvValue+1)));
      LRVD(:,n+j) = inv(LRVars)*diag(H'*P_init*H);
      %
      P0 = Q - (Q*H*inv(H'*Q*H)*H'*Q);
      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 = B*B';
      end;
      %
      if status==0;
         FEVDs(:,j+n,1) = inv(diag(FEVars(:,1)))*diag((H'*P1*H));
         Pi = P1;
         for i=2:h;
            Pi = (F*Pi*F')+(B*B');
            FEVDs(:,j+n,i) = inv(diag(FEVars(:,i)))*diag((H'*Pi*H));
         end;
      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 DSGEVarianceDecompFcn.m
%
