From e6c43c2a299d71ff3c1517945ab0c59ed210e38f Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 13 Jul 2023 12:01:02 -0400 Subject: [PATCH 1/2] load_last_mh_history_file: make sure number of blocks and parameters matches loaded files --- matlab/CutSample.m | 7 ++--- matlab/GetAllPosteriorDraws.m | 27 +++++++++-------- matlab/GetPosteriorMeanVariance.m | 2 +- matlab/GetPosteriorParametersStatistics.m | 29 +++++++++---------- matlab/compute_mh_covariance_matrix.m | 17 ++++------- .../convergence_diagnostics/McMCDiagnostics.m | 21 +++----------- matlab/generate_trace_plots.m | 4 +-- matlab/get_name_of_the_last_mh_file.m | 4 +-- matlab/marginal_density.m | 9 +++--- matlab/metropolis_draw.m | 4 +-- matlab/mh_autocorrelation_function.m | 8 +++-- matlab/posterior_sampler.m | 4 +-- matlab/posterior_sampler_initialization.m | 4 +-- matlab/prior_posterior_statistics.m | 9 +++--- matlab/selec_posterior_draws.m | 2 +- matlab/trace_plot.m | 5 ++-- 16 files changed, 68 insertions(+), 88 deletions(-) diff --git a/matlab/CutSample.m b/matlab/CutSample.m index 409f65cbf..e29f010d5 100644 --- a/matlab/CutSample.m +++ b/matlab/CutSample.m @@ -16,7 +16,7 @@ function CutSample(M_, options_, estim_params_) % SPECIAL REQUIREMENTS % none -% Copyright © 2005-2017 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -33,8 +33,6 @@ function CutSample(M_, options_, estim_params_) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx; - % Get the path to the metropolis files. MetropolisFolder = CheckPath('metropolis',M_.dname); @@ -42,7 +40,8 @@ MetropolisFolder = CheckPath('metropolis',M_.dname); ModelName = M_.fname; % Load the last mh-history file. -load_last_mh_history_file(MetropolisFolder, M_.fname); +record=load_last_mh_history_file(MetropolisFolder, M_.fname); +npar=size(record.LastParameters,2); % Get the list of files where the mcmc draw are saved. mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]); diff --git a/matlab/GetAllPosteriorDraws.m b/matlab/GetAllPosteriorDraws.m index bc6a21c37..f8cde8ef2 100644 --- a/matlab/GetAllPosteriorDraws.m +++ b/matlab/GetAllPosteriorDraws.m @@ -1,22 +1,23 @@ -function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, blck) - -% function Draws = GetAllPosteriorDraws(column,FirstMhFile,FirstLine,TotalNumberOfMhFile,NumberOfDraws) +function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck) +% function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck) % Gets all posterior draws % % INPUTS -% column: column +% column: column of desired parameter in draw matrix % FirstMhFile: first mh file % FirstLine: first line % TotalNumberOfMhFile: total number of mh file % NumberOfDraws: number of draws - +% nblcks: total number of blocks +% blck: desired block to read +% % OUTPUTS % Draws: draws from posterior distribution % % SPECIAL REQUIREMENTS % none -% Copyright © 2005-2017 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -33,19 +34,17 @@ function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumbe % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ options_ - -nblck = options_.mh_nblck; +global M_ iline = FirstLine; linee = 1; DirectoryName = CheckPath('metropolis',M_.dname); -if nblck>1 && nargin<6 - Draws = zeros(NumberOfDraws*nblck,1); +if nblcks>1 && nargin<7 + Draws = zeros(NumberOfDraws*nblcks,1); iline0=iline; if column>0 - for blck = 1:nblck + for blck = 1:nblcks iline=iline0; for file = FirstMhFile:TotalNumberOfMhFile load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2') @@ -56,7 +55,7 @@ if nblck>1 && nargin<6 end end else - for blck = 1:nblck + for blck = 1:nblcks iline=iline0; for file = FirstMhFile:TotalNumberOfMhFile load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2') @@ -68,7 +67,7 @@ if nblck>1 && nargin<6 end end else - if nblck==1 + if nblcks==1 blck=1; end if column>0 diff --git a/matlab/GetPosteriorMeanVariance.m b/matlab/GetPosteriorMeanVariance.m index fa095c6aa..562f39486 100644 --- a/matlab/GetPosteriorMeanVariance.m +++ b/matlab/GetPosteriorMeanVariance.m @@ -20,7 +20,7 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop) MetropolisFolder = CheckPath('metropolis',M.dname); FileName = M.fname; BaseName = [MetropolisFolder filesep FileName]; -load_last_mh_history_file(MetropolisFolder, FileName); +record=load_last_mh_history_file(MetropolisFolder, FileName); NbrDraws = sum(record.MhDraws(:,1)); NbrFiles = sum(record.MhDraws(:,2)); NbrBlocks = record.Nblck; diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index e33fae487..dee97ac76 100644 --- a/matlab/GetPosteriorParametersStatistics.m +++ b/matlab/GetPosteriorParametersStatistics.m @@ -16,7 +16,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay % SPECIAL REQUIREMENTS % None. -% Copyright © 2006-2018 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -38,19 +38,17 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay %end TeX = options_.TeX; -nblck = options_.mh_nblck; nvx = estim_params_.nvx; nvn = estim_params_.nvn; ncx = estim_params_.ncx; ncn = estim_params_.ncn; np = estim_params_.np ; -nx = nvx+nvn+ncx+ncn+np; MetropolisFolder = CheckPath('metropolis',M_.dname); OutputFolder = CheckPath('Output',M_.dname); FileName = M_.fname; -load_last_mh_history_file(MetropolisFolder,FileName); +record=load_last_mh_history_file(MetropolisFolder,FileName); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; @@ -58,6 +56,7 @@ TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); FirstMhFile = record.KeepedDraws.FirstMhFile; NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); +mh_nblck = size(record.LastParameters,1); clear record; header_width = row_header_width(M_, estim_params_, bayestopt_); @@ -76,7 +75,7 @@ catch disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean)) end -num_draws=NumberOfDraws*options_.mh_nblck; +num_draws=NumberOfDraws*mh_nblck; hpd_draws = round((1-options_.mh_conf_sig)*num_draws); if hpd_draws<2 @@ -97,7 +96,7 @@ if np ip = nvx+nvn+ncx+ncn+1; for i=1:np if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); name = bayestopt_.name{ip}; oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density); @@ -106,7 +105,7 @@ if np name = bayestopt_.name{ip}; [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type); catch - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); name = bayestopt_.name{ip}; oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density); @@ -140,7 +139,7 @@ if nvx disp(tit2) for i=1:nvx if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); k = estim_params_.var_exo(i,1); name = M_.exo_names{k}; @@ -152,7 +151,7 @@ if nvx name = M_.exo_names{k}; [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type); catch - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ... posterior_moments(Draws, 1, options_.mh_conf_sig); k = estim_params_.var_exo(i,1); @@ -184,7 +183,7 @@ if nvn ip = nvx+1; for i=1:nvn if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density); @@ -193,7 +192,7 @@ if nvn name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type); catch - Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws); + Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig); name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)}; oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density); @@ -223,7 +222,7 @@ if ncx ip = nvx+nvn+1; for i=1:ncx if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig); k1 = estim_params_.corrx(i,1); k2 = estim_params_.corrx(i,2); @@ -240,7 +239,7 @@ if ncx NAME = sprintf('%s_%s', M_.exo_names{k1}, M_.exo_names{k2}); [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type); catch - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); k1 = estim_params_.corrx(i,1); k2 = estim_params_.corrx(i,2); @@ -273,7 +272,7 @@ if ncn ip = nvx+nvn+ncx+1; for i=1:ncn if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) - Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws); + Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); k1 = estim_params_.corrn(i,1); k2 = estim_params_.corrn(i,2); @@ -288,7 +287,7 @@ if ncn NAME = sprintf('%s_%s', M_.endo_names{k1}, M_.endo_names{k2}); [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type); catch - Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig); k1 = estim_params_.corrn(i,1); k2 = estim_params_.corrn(i,2); diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m index b4435190a..c8f96e270 100644 --- a/matlab/compute_mh_covariance_matrix.m +++ b/matlab/compute_mh_covariance_matrix.m @@ -15,7 +15,7 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at % SPECIAL REQUIREMENTS % None. -% Copyright © 2006-2017 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -32,27 +32,20 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ options_ estim_params_ bayestopt_ - - -n = estim_params_.np + ... - estim_params_.nvn+ ... - estim_params_.ncx+ ... - estim_params_.ncn+ ... - estim_params_.nvx; - -nblck = options_.mh_nblck; +global M_ bayestopt_ MetropolisFolder = CheckPath('metropolis',M_.dname); ModelName = M_.fname; BaseName = [MetropolisFolder filesep ModelName]; -load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, ModelName); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); +[nblck, n] = size(record.LastParameters); + posterior_kernel_at_the_mode = -Inf; posterior_mean = zeros(n,1); posterior_mode = NaN(n,1); diff --git a/matlab/convergence_diagnostics/McMCDiagnostics.m b/matlab/convergence_diagnostics/McMCDiagnostics.m index 16f58ed85..1d0da5d61 100644 --- a/matlab/convergence_diagnostics/McMCDiagnostics.m +++ b/matlab/convergence_diagnostics/McMCDiagnostics.m @@ -16,7 +16,7 @@ function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) % PARALLEL CONTEXT % See the comment in posterior_sampler.m funtion. -% Copyright © 2005-2018 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -38,25 +38,12 @@ MetropolisFolder = CheckPath('metropolis',M_.dname); ModelName = M_.fname; TeX = options_.TeX; -nblck = options_.mh_nblck; -npar = estim_params_.nvx; -npar = npar + estim_params_.nvn; -npar = npar + estim_params_.ncx; -npar = npar + estim_params_.ncn; -npar = npar + estim_params_.np ; -MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); -load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, ModelName); NumberOfMcFilesPerBlock = record.LastFileNumber; +[nblck, npar] = size(record.LastParameters); npardisp = options_.convergence.brooksgelman.plotrows; -% Check that the declared number of blocks is consistent with informations saved in mh-history files. -if ~isequal(nblck,record.Nblck) - disp(['Estimation::mcmc::diagnostics: The number of MCMC chains you declared (' num2str(nblck) ') is inconsistent with the information available in the mh-history files (' num2str(record.Nblck) ' chains)!']) - disp([' I reset the number of MCMC chains to ' num2str(record.Nblck) '.']) - nblck = record.Nblck; -end - % check if all the mh files are available. issue_an_error_message = 0; for b = 1:nblck @@ -91,7 +78,7 @@ for jj = 1:npar par_name_temp = get_the_name(jj, options_.TeX, M_, estim_params_, options_); param_name = vertcat(param_name, par_name_temp); end - Draws = GetAllPosteriorDraws(jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + Draws = GetAllPosteriorDraws(jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck); Draws = reshape(Draws, [NumberOfDraws nblck]); Nc = min(1000, NumberOfDraws/2); for ll = 1:nblck diff --git a/matlab/generate_trace_plots.m b/matlab/generate_trace_plots.m index 2fcbda4f3..a241cc076 100644 --- a/matlab/generate_trace_plots.m +++ b/matlab/generate_trace_plots.m @@ -11,7 +11,7 @@ function generate_trace_plots(chain_number) % SPECIAL REQUIREMENTS % none -% Copyright © 2016-2018 Dynare Team +% Copyright © 2016-2023 Dynare Team % % This file is part of Dynare. % @@ -33,7 +33,7 @@ global M_ options_ estim_params_ % Get informations about the posterior draws: MetropolisFolder = CheckPath('metropolis', M_.dname); -load_last_mh_history_file(MetropolisFolder, M_.fname); +record=load_last_mh_history_file(MetropolisFolder, M_.fname); if chain_number>record.Nblck error('generate_trace_plots:: chain number is bigger than existing number of chains') end diff --git a/matlab/get_name_of_the_last_mh_file.m b/matlab/get_name_of_the_last_mh_file.m index 5d49a780b..6c9ccca24 100644 --- a/matlab/get_name_of_the_last_mh_file.m +++ b/matlab/get_name_of_the_last_mh_file.m @@ -11,7 +11,7 @@ function [mhname,info] = get_name_of_the_last_mh_file(M_) % file. Otherwise info is equal to zero (a likely reason for this is % that the mcmc simulations were not completed). -% Copyright © 2008-2017 Dynare Team +% Copyright © 2008-2023 Dynare Team % % This file is part of Dynare. % @@ -35,7 +35,7 @@ MetropolisFolder = CheckPath('metropolis',M_.dname); ModelName = M_.fname; BaseName = [MetropolisFolder filesep ModelName]; -load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, ModelName); mh_number = record.LastFileNumber ; bk_number = record.Nblck ; diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m index 8bca13e66..8eb56bed4 100644 --- a/matlab/marginal_density.m +++ b/matlab/marginal_density.m @@ -15,7 +15,7 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay % SPECIAL REQUIREMENTS % none -% Copyright © 2005-2018 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -33,14 +33,13 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay % along with Dynare. If not, see . -npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx; -nblck = options_.mh_nblck; - MetropolisFolder = CheckPath('metropolis',M_.dname); ModelName = M_.fname; BaseName = [MetropolisFolder filesep ModelName]; -load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, ModelName); + +[nblck, npar] = size(record.LastParameters); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine; diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m index 533d79ad3..b5dbe38e4 100644 --- a/matlab/metropolis_draw.m +++ b/matlab/metropolis_draw.m @@ -24,7 +24,7 @@ function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params % % Requires CutSample to be run before in order to set up mh_history-file -% Copyright © 2003-2017 Dynare Team +% Copyright © 2003-2023 Dynare Team % % This file is part of Dynare. % @@ -59,7 +59,7 @@ if init FileName = M_.fname; BaseName = [MetropolisFolder filesep FileName]; %load mh_history-file with info on what to load - load_last_mh_history_file(MetropolisFolder, FileName); + record=load_last_mh_history_file(MetropolisFolder, FileName); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); diff --git a/matlab/mh_autocorrelation_function.m b/matlab/mh_autocorrelation_function.m index 16524a190..8c07143bd 100644 --- a/matlab/mh_autocorrelation_function.m +++ b/matlab/mh_autocorrelation_function.m @@ -18,7 +18,7 @@ function mh_autocorrelation_function(options_,M_,estim_params_,type,blck,name1,n % % SPECIAL REQUIREMENTS -% Copyright © 2003-2017 Dynare Team +% Copyright © 2003-2023 Dynare Team % % This file is part of Dynare. % @@ -48,17 +48,19 @@ end % Get informations about the posterior draws: MetropolisFolder = CheckPath('metropolis',M_.dname); -load_last_mh_history_file(MetropolisFolder, M_.fname); +record=load_last_mh_history_file(MetropolisFolder, M_.fname); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); +nblck = size(record.LastParameters,1); + clear record; % Get all the posterior draws: -PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, blck); +PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck, blck); % Compute the autocorrelation function: [autocov,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size); diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m index 2942a6129..d16ff7540 100644 --- a/matlab/posterior_sampler.m +++ b/matlab/posterior_sampler.m @@ -37,7 +37,7 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun % Then the comments write here can be used for all the other pairs of % parallel functions and also for management functions. -% Copyright © 2006-2017 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -62,7 +62,7 @@ vv = sampler_options.invhess; InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2); % Load last mh history file -load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, ModelName); % Only for test parallel results!!! diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index 60a909d2c..9a8190269 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -340,7 +340,7 @@ if ~options_.load_mh_file && ~options_.mh_recover elseif options_.load_mh_file && ~options_.mh_recover % Here we consider previous mh files (previous mh did not crash). disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...') - load_last_mh_history_file(MetropolisFolder, ModelName); + record=load_last_mh_history_file(MetropolisFolder, ModelName); if ~isnan(record.MCMCConcludedSuccessfully) && ~record.MCMCConcludedSuccessfully error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.') end @@ -400,7 +400,7 @@ elseif options_.load_mh_file && ~options_.mh_recover elseif options_.mh_recover % The previous metropolis-hastings crashed before the end! I try to recover the saved draws... disp('Estimation::mcmc: Recover mode!') - load_last_mh_history_file(MetropolisFolder, ModelName); + record=load_last_mh_history_file(MetropolisFolder, ModelName); NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains. options_.mh_nblck = NumberOfBlocks; diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index e93387afb..5f2b61db5 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -19,7 +19,7 @@ function prior_posterior_statistics(type,dataset,dataset_info) % See the comments in the posterior_sampler.m funtion. -% Copyright © 2005-2020 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -209,7 +209,8 @@ localVars.DirectoryName = DirectoryName; if strcmpi(type,'posterior') BaseName = [DirectoryName filesep M_.fname]; - load_last_mh_history_file(DirectoryName, M_.fname); + record=load_last_mh_history_file(DirectoryName, M_.fname); + [nblck, npar] = size(record.LastParameters); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); @@ -219,9 +220,9 @@ if strcmpi(type,'posterior') mh_nblck = options_.mh_nblck; if B==NumberOfDraws*mh_nblck % we load all retained MH runs ! - logpost=GetAllPosteriorDraws(0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + logpost=GetAllPosteriorDraws(0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck); for column=1:npar - x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); + x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck); end else logpost=NaN(B,1); diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m index 0a10dab82..0b5dadcc5 100644 --- a/matlab/selec_posterior_draws.m +++ b/matlab/selec_posterior_draws.m @@ -67,7 +67,7 @@ ModelName = M_.fname; BaseName = [MetropolisFolder filesep ModelName]; % Get informations about the mcmc: -load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, ModelName); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); diff --git a/matlab/trace_plot.m b/matlab/trace_plot.m index af4a6b876..8deecb817 100644 --- a/matlab/trace_plot.m +++ b/matlab/trace_plot.m @@ -53,16 +53,17 @@ end % Get informations about the posterior draws: MetropolisFolder = CheckPath('metropolis',M_.dname); -load_last_mh_history_file(MetropolisFolder, M_.fname); +record=load_last_mh_history_file(MetropolisFolder, M_.fname); FirstMhFile = 1; FirstLine = 1; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); +[mh_nblck] = size(record.LastParameters,2); clear record; % Get all the posterior draws: -PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, blck); +PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck); % Plot the posterior draws: From 8532d6abd73b86b1c8ff3540e041cbb3c6c6275d Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 13 Jul 2023 14:05:52 -0400 Subject: [PATCH 2/2] Cosmetic changes to various estimation routines --- matlab/GetPosteriorMeanVariance.m | 22 +++++++++++--- matlab/GetPosteriorParametersStatistics.m | 24 +++++++-------- .../convergence_diagnostics/McMCDiagnostics.m | 18 +++++------- matlab/get_name_of_the_last_mh_file.m | 2 +- matlab/marginal_density.m | 17 +++++------ matlab/metropolis_draw.m | 4 +-- matlab/posterior_sampler.m | 29 +++++++++---------- matlab/posterior_sampler_initialization.m | 9 +----- matlab/prior_posterior_statistics.m | 10 +++---- 9 files changed, 65 insertions(+), 70 deletions(-) diff --git a/matlab/GetPosteriorMeanVariance.m b/matlab/GetPosteriorMeanVariance.m index 562f39486..5a2c44af2 100644 --- a/matlab/GetPosteriorMeanVariance.m +++ b/matlab/GetPosteriorMeanVariance.m @@ -1,6 +1,20 @@ -function [mean,variance] = GetPosteriorMeanVariance(M,drop) +function [mean,variance] = GetPosteriorMeanVariance(M_,drop) +%function [mean,variance] = GetPosteriorMeanVariance(M,drop) +% Computes the posterior mean and variance +% (+updates of oo_ & TeX output). +% +% INPUTS +% M_ [structure] Dynare model structure +% drop [double] share of draws to drop +% +% OUTPUTS +% mean [double] mean parameter vector +% variance [double] variance +% +% SPECIAL REQUIREMENTS +% None. -% Copyright © 2012-2017 Dynare Team +% Copyright © 2012-2023 Dynare Team % % This file is part of Dynare. % @@ -17,8 +31,8 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -MetropolisFolder = CheckPath('metropolis',M.dname); -FileName = M.fname; +MetropolisFolder = CheckPath('metropolis',M_.dname); +FileName = M_.fname; BaseName = [MetropolisFolder filesep FileName]; record=load_last_mh_history_file(MetropolisFolder, FileName); NbrDraws = sum(record.MhDraws(:,1)); diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index dee97ac76..b7713b66d 100644 --- a/matlab/GetPosteriorParametersStatistics.m +++ b/matlab/GetPosteriorParametersStatistics.m @@ -1,17 +1,18 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames) +% function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames) % This function prints and saves posterior estimates after the mcmc % (+updates of oo_ & TeX output). % % INPUTS -% estim_params_ [structure] -% M_ [structure] -% options_ [structure] -% bayestopt_ [structure] -% oo_ [structure] -% pnames [cell] cell of strings, names of the prior shapes available +% estim_params_ [structure] Dynare estimation parameter structure +% M_ [structure] Dynare model structure +% options_ [structure] Dynare options structure +% bayestopt_ [structure] Dynare structure describing priors +% oo_ [structure] Dynare results structure +% pnames [cell] cell of strings, names of the prior shapes available % % OUTPUTS -% oo_ [structure] +% oo_ [structure] Dynare results structure % % SPECIAL REQUIREMENTS % None. @@ -50,7 +51,6 @@ FileName = M_.fname; record=load_last_mh_history_file(MetropolisFolder,FileName); -FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); @@ -68,12 +68,10 @@ skipline(2) disp('ESTIMATION RESULTS') skipline() -try - disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean)) -catch - [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_); - disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean)) +if ~isfield(oo_,'MarginalDensity') || ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean') + [~,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_); end +fprintf('\nLog data density is %f.\n', oo_.MarginalDensity.ModifiedHarmonicMean); num_draws=NumberOfDraws*mh_nblck; hpd_draws = round((1-options_.mh_conf_sig)*num_draws); diff --git a/matlab/convergence_diagnostics/McMCDiagnostics.m b/matlab/convergence_diagnostics/McMCDiagnostics.m index 1d0da5d61..cbe02ce8d 100644 --- a/matlab/convergence_diagnostics/McMCDiagnostics.m +++ b/matlab/convergence_diagnostics/McMCDiagnostics.m @@ -1,11 +1,12 @@ function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) -% function McMCDiagnostics +% function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) % Computes convergence tests % % INPUTS -% options_ [structure] -% estim_params_ [structure] -% M_ [structure] +% options_ [structure] Dynare options structure +% estim_params_ [structure] Dynare estimation parameter structure +% M_ [structure] Dynare model structure +% oo_ [structure] Dynare results structure % % OUTPUTS % oo_ [structure] @@ -58,7 +59,6 @@ if issue_an_error_message end % compute inefficiency factor -FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); @@ -68,6 +68,7 @@ NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws param_name = {}; param_name_tex = {}; +Ifac=NaN(nblck,npar); for jj = 1:npar if options_.TeX [par_name_temp, par_name_tex_temp] = get_the_name(jj, options_.TeX, M_,estim_params_, options_); @@ -107,8 +108,6 @@ update_last_mh_history_file(MetropolisFolder, ModelName, record); PastDraws = sum(record.MhDraws,1); -LastFileNumber = PastDraws(2); -LastLineNumber = record.MhDraws(end,3); NumberOfDraws = PastDraws(1); if NumberOfDraws<=2000 @@ -217,8 +216,6 @@ ALPHA = 0.2; % increase too much with the number time = 1:NumberOfDraws; xx = Origin:StepSize:NumberOfDraws; NumberOfLines = length(xx); -tmp = zeros(NumberOfDraws*nblck,3); -UDIAG = zeros(NumberOfLines,6,npar); if NumberOfDraws < Origin disp('Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.') @@ -256,7 +253,6 @@ if isnumeric(options_.parallel) clear fout % Parallel execution! else - ModelName = ModelName; if ~isempty(M_.bvar) ModelName = [ModelName '_bvar']; end @@ -391,7 +387,7 @@ if reste fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1)); fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,'\n'); - fprintf(fidTeX,'% End Of TeX file.'); + fprintf(fidTeX,'%% End Of TeX file.'); fclose(fidTeX); end end % if reste > 0 diff --git a/matlab/get_name_of_the_last_mh_file.m b/matlab/get_name_of_the_last_mh_file.m index 6c9ccca24..1ddf7a060 100644 --- a/matlab/get_name_of_the_last_mh_file.m +++ b/matlab/get_name_of_the_last_mh_file.m @@ -1,4 +1,5 @@ function [mhname,info] = get_name_of_the_last_mh_file(M_) +% function [mhname,info] = get_name_of_the_last_mh_file(M_) % This function returns the name of the last mh file and test if the metropolis was completed. % % INPUTS @@ -28,7 +29,6 @@ function [mhname,info] = get_name_of_the_last_mh_file(M_) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -mhname = []; info = 1; MetropolisFolder = CheckPath('metropolis',M_.dname); diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m index 8eb56bed4..48014ba61 100644 --- a/matlab/marginal_density.m +++ b/matlab/marginal_density.m @@ -1,16 +1,16 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_) -% function marginal = marginal_density() +% function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_) % Computes the marginal density % % INPUTS -% options_ [structure] -% estim_params_ [structure] -% M_ [structure] -% oo_ [structure] +% options_ [structure] Dynare options structure +% estim_params_ [structure] Dynare estimation parameter structure +% M_ [structure] Dynare model structure +% oo_ [structure] Dynare results structure % % OUTPUTS -% marginal: [double] marginal density (modified harmonic mean) -% oo_ [structure] +% marginal: [double] marginal density (modified harmonic mean) +% oo_ [structure] Dynare results structure % % SPECIAL REQUIREMENTS % none @@ -42,10 +42,9 @@ record=load_last_mh_history_file(MetropolisFolder, ModelName); [nblck, npar] = size(record.LastParameters); FirstMhFile = record.KeepedDraws.FirstMhFile; -FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine; +FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); -MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws); fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... '); diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m index b5dbe38e4..690dcac5d 100644 --- a/matlab/metropolis_draw.m +++ b/matlab/metropolis_draw.m @@ -1,5 +1,5 @@ function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_) -% function [xparams, logpost]=metropolis_draw(init) +% function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_) % Builds draws from metropolis % % INPUTS: @@ -62,8 +62,6 @@ if init record=load_last_mh_history_file(MetropolisFolder, FileName); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; - TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); - LastMhFile = TotalNumberOfMhFiles; TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ? diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m index d16ff7540..fbc395895 100644 --- a/matlab/posterior_sampler.m +++ b/matlab/posterior_sampler.m @@ -3,21 +3,20 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun % Random Walk Metropolis-Hastings algorithm. % % INPUTS -% o TargetFun [char] string specifying the name of the objective -% function (posterior kernel). -% o ProposalFun [char] string specifying the name of the proposal -% density -% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). -% o sampler_options structure -% .invhess [double] (p*p) matrix, posterior covariance matrix (at the mode). -% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. -% o dataset_ data structure -% o dataset_info dataset info structure -% o options_ options structure -% o M_ model structure -% o estim_params_ estimated parameters structure -% o bayestopt_ estimation options structure -% o oo_ outputs structure +% o TargetFun [char] string specifying the name of the objective +% function (posterior kernel). +% o ProposalFun [char] string specifying the name of the proposal +% density +% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). +% o sampler_options structure +% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. +% o dataset_ [structure] data structure +% o dataset_info [structure] dataset info structure +% o options_ [structure] options structure +% o M_ [structure] model structure +% o estim_params_ [structure] estimated parameters structure +% o bayestopt_ [structure] prior specification structure +% o oo_ [structure] output structure % % SPECIAL REQUIREMENTS % None. diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index 9a8190269..2ac7f6963 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -58,16 +58,9 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa %Initialize outputs ix2 = []; ilogpo2 = []; -ModelName = []; -MetropolisFolder = []; FirstBlock = []; FirstLine = []; -npar = []; -NumberOfBlocks = []; -nruns = []; NewFile = []; -MAX_nruns = []; -d = []; ModelName = M_.fname; if ~isempty(M_.bvar) @@ -171,7 +164,7 @@ if ~options_.load_mh_file && ~options_.mh_recover else bayestopt0 = load([PreviousFolder0 filesep 'prior' filesep 'definition.mat']); end - [common_parameters,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name); + [~,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name); new_estimated_parameters = ~ismember(bayestopt_.name,bayestopt0.bayestopt_.name); ix2 = zeros(NumberOfBlocks,npar); ilogpo2 = zeros(NumberOfBlocks,1); diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index 5f2b61db5..51ab6f02d 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -1,6 +1,5 @@ function prior_posterior_statistics(type,dataset,dataset_info) - -% function prior_posterior_statistics(type,dataset) +% function prior_posterior_statistics(type,dataset,dataset_info) % Computes Monte Carlo filter smoother and forecasts % % INPUTS @@ -8,6 +7,7 @@ function prior_posterior_statistics(type,dataset,dataset_info) % prior % gsa % dataset: data structure +% dataset_info: dataset structure % % OUTPUTS % none @@ -208,13 +208,11 @@ localVars.ifil=ifil; localVars.DirectoryName = DirectoryName; if strcmpi(type,'posterior') - BaseName = [DirectoryName filesep M_.fname]; record=load_last_mh_history_file(DirectoryName, M_.fname); [nblck, npar] = size(record.LastParameters); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); - LastMhFile = TotalNumberOfMhFiles; TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); mh_nblck = options_.mh_nblck; @@ -318,9 +316,9 @@ if strcmpi(type,'gsa') return end -if ~isnumeric(options_.parallel), +if ~isnumeric(options_.parallel) leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen; - if options_.parallel_info.leaveSlaveOpen == 0, + if options_.parallel_info.leaveSlaveOpen == 0 % Commenting for testing!!! options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls) end