From c356db4531a5e8786daa4a25b274de7a5ce34942 Mon Sep 17 00:00:00 2001 From: Willi Mutschler Date: Fri, 1 Sep 2023 16:13:23 +0200 Subject: [PATCH] Factorize estimation: set_mcmc_jumping_covariance --- matlab/dynare_estimation_1.m | 33 +-------- .../estimation/set_mcmc_jumping_covariance.m | 67 +++++++++++++++++++ 2 files changed, 68 insertions(+), 32 deletions(-) create mode 100644 matlab/utilities/estimation/set_mcmc_jumping_covariance.m diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index acd2f0135..f59fd7100 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -376,38 +376,7 @@ if np > 0 save([M_.dname filesep 'Output' filesep M_.fname '_params.mat'],'pindx'); end -switch options_.MCMC_jumping_covariance - case 'hessian' %Baseline - %do nothing and use hessian from mode_compute - case 'prior_variance' %Use prior variance - if any(isinf(bayestopt_.p2)) - error('Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.') - else - hh = diag(1./(bayestopt_.p2.^2)); - end - hsd = sqrt(diag(hh)); - invhess = inv(hh./(hsd*hsd'))./(hsd*hsd'); - case 'identity_matrix' %Use identity - invhess = eye(nx); - otherwise %user specified matrix in file - try - load(options_.MCMC_jumping_covariance,'jumping_covariance') - hh=jumping_covariance; - catch - error(['No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat']) - end - [nrow, ncol]=size(hh); - if ~isequal(nrow,ncol) && ~isequal(nrow,nx) %check if square and right size - error(['jumping_covariance matrix must be square and have ',num2str(nx),' rows and columns']) - end - try %check for positive definiteness - chol(hh); - hsd = sqrt(diag(hh)); - invhess = inv(hh./(hsd*hsd'))./(hsd*hsd'); - catch - error(['Specified jumping_covariance is not positive definite']) - end -end +invhess = set_mcmc_jumping_covariance(invhess, nx, options_.MCMC_jumping_covariance, bayestopt_, 'dynare_estimation_1'); if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... (any(bayestopt_.pshape >0 ) && options_.load_mh_file) %% not ML estimation diff --git a/matlab/utilities/estimation/set_mcmc_jumping_covariance.m b/matlab/utilities/estimation/set_mcmc_jumping_covariance.m new file mode 100644 index 000000000..b0cb99d49 --- /dev/null +++ b/matlab/utilities/estimation/set_mcmc_jumping_covariance.m @@ -0,0 +1,67 @@ +function invhess = set_mcmc_jumping_covariance(invhess, xparam_nbr, MCMC_jumping_covariance, bayestopt_, stringForErrors) +% function invhess = set_mcmc_jumping_covariance(invhess, xparam_nbr, MCMC_jumping_covariance, bayestopt_, stringForErrors) +% ------------------------------------------------------------------------- +% sets the jumping covariance matrix for the MCMC algorithm +% ========================================================================= +% INPUTS +% o invhess: [matrix] already computed inverse of the hessian matrix +% o xparam_nbr: [integer] number of estimated parameters +% o MCMC_jumping_covariance: [string] name of option or file setting the jumping covariance matrix +% o bayestopt_: [struct] information on priors +% o stringForErrors: [string] string to be used in error messages +% ------------------------------------------------------------------------- +% OUTPUTS +% o invhess: [matrix] jumping covariance matrix +% ------------------------------------------------------------------------- +% This function is called by +% o dynare_estimation_1.m +% ------------------------------------------------------------------------- +% 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 . +% ========================================================================= + +switch MCMC_jumping_covariance + case 'hessian' % do nothing and use hessian from previous mode optimization + case 'prior_variance' % use prior variance + if any(isinf(bayestopt_.p2)) + error('%s: Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.',stringForErrors); + else + hess = diag(1./(bayestopt_.p2.^2)); + end + hsd = sqrt(diag(hess)); + invhess = inv(hess./(hsd*hsd'))./(hsd*hsd'); + case 'identity_matrix' % use identity + invhess = eye(xparam_nbr); + otherwise % user specified matrix in file + try + load(MCMC_jumping_covariance,'jumping_covariance') + hess = jumping_covariance; + catch + error(['%s: No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat!'],stringForErrors); + end + [nrow, ncol] = size(hess); + if ~isequal(nrow,ncol) && ~isequal(nrow,xparam_nbr) % check if square and right size + error(['%s: ''jumping_covariance'' matrix (loaded from ',options_.MCMC_jumping_covariance,'.mat) must be square and have ',num2str(xparam_nbr),' rows and columns!'],stringForErrors); + end + try % check for positive definiteness + chol(hess); + hsd = sqrt(diag(hess)); + invhess = inv(hess./(hsd*hsd'))./(hsd*hsd'); + catch + error('%s: Specified ''MCMC_jumping_covariance'' is not positive definite!',stringForErrors); + end +end \ No newline at end of file