From d386bb9f76cc20b5d879ca961a3f2b3c8713fa71 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 23 Jun 2023 09:57:17 -0400 Subject: [PATCH] Cholesky decomposition: only add to diagonal if really necessary Closes #1891 --- matlab/PosteriorIRF_core1.m | 5 +- matlab/backward/backward_model_forecast.m | 3 +- matlab/backward/backward_model_irf.m | 6 +- matlab/default_option_values.m | 1 + matlab/get_lower_cholesky_covariance.m | 92 +++++++++++++++++++++++ matlab/stoch_simul.m | 5 +- matlab/th_autocovariances.m | 6 +- 7 files changed, 102 insertions(+), 16 deletions(-) create mode 100644 matlab/get_lower_cholesky_covariance.m diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index 15d9fd238..a0a47dd5a 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -177,11 +177,10 @@ while fpar 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 diff --git a/matlab/backward/backward_model_forecast.m b/matlab/backward/backward_model_forecast.m index 066df7757..dd032ca42 100644 --- a/matlab/backward/backward_model_forecast.m +++ b/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 diff --git a/matlab/backward/backward_model_irf.m b/matlab/backward/backward_model_irf.m index c2accdee0..8e36457f4 100644 --- a/matlab/backward/backward_model_irf.m +++ b/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 diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m index 0c3616b91..66cec2f3f 100644 --- a/matlab/default_option_values.m +++ b/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; diff --git a/matlab/get_lower_cholesky_covariance.m b/matlab/get_lower_cholesky_covariance.m new file mode 100644 index 000000000..05a85564b --- /dev/null +++ b/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 . + +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 diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index f4dbdb731..742b19f7e 100644 --- a/matlab/stoch_simul.m +++ b/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 diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m index 57f978194..eb65105bb 100644 --- a/matlab/th_autocovariances.m +++ b/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,:);