Cholesky decomposition: only add to diagonal if really necessary

Closes #1891

matlab/PosteriorIRF_core1.m

@@ -177,11 +177,10 @@ while fpar<B
             error(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ') with sample ' type])
         end
     end
-    SS = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
-    SS = transpose(chol(SS));
+    SS = get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
     irf_shocks_indx = getIrfShocksIndx(M_, options_);
     for i=irf_shocks_indx
-        if SS(i,i) > 1e-13
+        if SS(i,i) > 5e-7
             if options_.order>1 && options_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
                 y=irf(M_,options_,dr,SS(:,i)./SS(i,i)/100, options_.irf, options_.drop,options_.replic,options_.order);
             else

matlab/backward/backward_model_forecast.m

@@ -85,8 +85,7 @@ end
 
 % Get the covariance matrix of the shocks.
 if withuncertainty
-    Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
-    sigma = transpose(chol(Sigma));
+    sigma = get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
 end
 
 % Compute forecast without shock

matlab/backward/backward_model_irf.m

@@ -146,10 +146,8 @@ end
 if ~deterministicshockflag
     if nnz(M_.Sigma_e)
         % Add ϵ>0 on the diagonal, so that the Cholesky won't fail
-        % if a shock has zero variance
-        Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
-        % Factorize Sigma (C is such that C*C' == Sigma)
-        C = chol(Sigma, 'lower');
+        % if a shock has zero variance and factorize Sigma (C is such that C*C' == Sigma)
+        C = get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
     else
         error('You did not specify the size of the shocks!')
     end

matlab/default_option_values.m

@@ -67,6 +67,7 @@ options_.mode_check.number_of_points = 20;
 options_.mode_check.nolik = false;
 
 options_.huge_number = 1e7;
+options_.add_tiny_number_to_cholesky=1e-14;
 
 % Default number of threads for parallelized mex files.
 options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = num_procs;

matlab/get_lower_cholesky_covariance.m

@@ -0,0 +1,92 @@
+function chol_sigma=get_lower_cholesky_covariance(Sigma_e,add_tiny_number_to_cholesky)
+% function chol_sigma=get_lower_cholesky_covariance(Sigma_e)
+% Computes the lower triangular Cholesky decomposition of a covariance matrix,
+% working around zero entries on the diagonal and perfect correlation
+%
+% INPUTS
+%   Sigma_e         [double]      covariance matrix
+%
+% OUTPUTS
+%   chol_sigma      [cell]        Cholesky factor 
+%
+% ALGORITHM
+%   Add small value to diagonal to break perfect correlation
+%
+% SPECIAL REQUIREMENTS.
+%   None.
+%
+% Copyright © 2023 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare 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.
+%
+% Dynare 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 Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+if nargin<2
+    add_tiny_number_to_cholesky=1e-14;
+end
+std_deviation=sqrt(diag(Sigma_e));
+non_zero_indices=find(std_deviation~=0); %find non-zero shocks;
+try    
+    chol_sigma=zeros(size(Sigma_e));
+    chol_sigma(non_zero_indices,non_zero_indices)=chol(Sigma_e(non_zero_indices,non_zero_indices),'lower');
+catch
+    % cases with perfect correlation
+    fprintf('Non-positive definite covariance matrix encountered. Using add_tiny_number_to_cholesky one the diagonal.\n')
+    chol_sigma=zeros(size(Sigma_e));
+    chol_sigma(non_zero_indices,non_zero_indices)=chol(Sigma_e(non_zero_indices,non_zero_indices)+add_tiny_number_to_cholesky*eye(length(non_zero_indices)),'lower');
+    % correlation=diag(std_deviation(non_zero_indices))\Sigma_e(non_zero_indices,non_zero_indices)/diag(std_deviation(non_zero_indices));
+end
+
+return % --*-- Unit tests --*--
+
+%@test:1
+
+Sigma_e=diag(4*ones(3,1));
+Sigma_e(2,2)=0;
+chol_1=get_lower_cholesky_covariance(Sigma_e);
+if max(max(abs(chol_1-diag([2,0,2]))))>eps 
+    t(1)=false;
+else
+    t(1)=true;
+end
+
+Sigma_e=ones(3,3);
+chol_2=get_lower_cholesky_covariance(Sigma_e,1e-14);
+chol_3=get_lower_cholesky_covariance(Sigma_e+1e-14*eye(3),1e-14);
+if max(max(abs(chol_2-chol_3)))>eps || any(any(triu(chol_3,1)))
+    t(2)=false;
+else
+    t(2)=true;
+end
+
+Sigma_e=ones(3,3);
+Sigma_e(2,:)=0;
+Sigma_e(:,2)=0;
+chol_4=get_lower_cholesky_covariance(Sigma_e,1e-14);
+if chol_4(2,2)~=0 || any(any(triu(chol_4,1)))
+    t(3)=false;
+else
+    t(3)=true;
+end   
+
+Sigma_e=[4 0.5 0; 0.5 9 0; 0 0 16];
+chol_5=get_lower_cholesky_covariance(Sigma_e,1e-14);
+if any(any(triu(chol_5,1))) %should be lower triangular
+    t(4)=false;
+else
+    t(4)=true;
+end
+
+T = all(t);
+%@eof:1

matlab/stoch_simul.m

@@ -231,15 +231,14 @@ if options_.irf
         fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
         fprintf(fidTeX,' \n');
     end
-    SS=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
-    cs = transpose(chol(SS));
+    cs=get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
     tit = M_.exo_names;
     if TeX
         titTeX = M_.exo_names_tex;
     end
     irf_shocks_indx = getIrfShocksIndx(M_, options_);
     for i=irf_shocks_indx
-        if SS(i,i) > 1e-13
+        if cs(i,i) > 5e-7            
             if PI_PCL_solver
                 y=PCL_Part_info_irf (0, PCL_varobs, i_var, M_, oo_.dr, options_.irf, i);
             else

matlab/th_autocovariances.m

@@ -172,8 +172,7 @@ if options_.hp_filter == 0 && ~options_.bandpass.indicator
             Gamma_y{nar+2} = ones(nvar,1);
         else
             Gamma_y{nar+2} = NaN(nvar,M_.exo_nbr);
-            SS=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
-            cs = chol(SS)';
+            cs = get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
             b1 = ghu1;
             b1 = b1*cs;
             b2 = ghu(iky,:);
@@ -250,8 +249,7 @@ else% ==> Theoretical filters.
             Gamma_y{nar+2} = ones(nvar,1);
         else
             Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
-            SS = M_.Sigma_e+1e-14*eye(M_.exo_nbr); %make sure Covariance matrix is positive definite
-            cs = chol(SS)';
+            cs = get_lower_cholesky_covariance(M_.Sigma_e); %make sure Covariance matrix is positive definite
             SS = cs*cs';
             b1 = ghu1;
             b2 = ghu(iky,:);