stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

load_last_mh_history_file: make sure number of blocks and parameters matches loaded files

kalman-mex
Johannes Pfeifer 2023-07-13 12:01:02 -04:00
parent b5a74b96e5
commit e6c43c2a29
16 changed files with 68 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
matlab/CutSample.m
View File

@ -16,7 +16,7 @@ function CutSample(M_, options_, estim_params_)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2005-2017 Dynare Team % Copyright © 2005-2023 Dynare Team
% %
% This file is part of Dynare. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
% Get the path to the metropolis files. % Get the path to the metropolis files.
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
@ -42,7 +40,8 @@ MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname; ModelName = M_.fname;
% Load the last mh-history file. % 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. % Get the list of files where the mcmc draw are saved.
mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]); mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]);

27
matlab/GetAllPosteriorDraws.m
View File

@ -1,22 +1,23 @@
function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, blck) function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
% function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
% function Draws = GetAllPosteriorDraws(column,FirstMhFile,FirstLine,TotalNumberOfMhFile,NumberOfDraws)
% Gets all posterior draws % Gets all posterior draws
% %
% INPUTS % INPUTS
% column: column % column: column of desired parameter in draw matrix
% FirstMhFile: first mh file % FirstMhFile: first mh file
% FirstLine: first line % FirstLine: first line
% TotalNumberOfMhFile: total number of mh file % TotalNumberOfMhFile: total number of mh file
% NumberOfDraws: number of draws % NumberOfDraws: number of draws
% nblcks: total number of blocks
% blck: desired block to read
%
% OUTPUTS % OUTPUTS
% Draws: draws from posterior distribution % Draws: draws from posterior distribution
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2005-2017 Dynare Team % Copyright © 2005-2023 Dynare Team
% %
% This file is part of Dynare. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ options_ global M_
nblck = options_.mh_nblck;
iline = FirstLine; iline = FirstLine;
linee = 1; linee = 1;
DirectoryName = CheckPath('metropolis',M_.dname); DirectoryName = CheckPath('metropolis',M_.dname);
if nblck>1 && nargin<6 if nblcks>1 && nargin<7
Draws = zeros(NumberOfDraws*nblck,1); Draws = zeros(NumberOfDraws*nblcks,1);
iline0=iline; iline0=iline;
if column>0 if column>0
for blck = 1:nblck for blck = 1:nblcks
iline=iline0; iline=iline0;
for file = FirstMhFile:TotalNumberOfMhFile for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2') load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
@ -56,7 +55,7 @@ if nblck>1 && nargin<6
end end
end end
else else
for blck = 1:nblck for blck = 1:nblcks
iline=iline0; iline=iline0;
for file = FirstMhFile:TotalNumberOfMhFile for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2') load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
@ -68,7 +67,7 @@ if nblck>1 && nargin<6
end end
end end
else else
if nblck==1 if nblcks==1
blck=1; blck=1;
end end
if column>0 if column>0

2
matlab/GetPosteriorMeanVariance.m
View File

@ -20,7 +20,7 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop)
MetropolisFolder = CheckPath('metropolis',M.dname); MetropolisFolder = CheckPath('metropolis',M.dname);
FileName = M.fname; FileName = M.fname;
BaseName = [MetropolisFolder filesep FileName]; BaseName = [MetropolisFolder filesep FileName];
load_last_mh_history_file(MetropolisFolder, FileName); record=load_last_mh_history_file(MetropolisFolder, FileName);
NbrDraws = sum(record.MhDraws(:,1)); NbrDraws = sum(record.MhDraws(:,1));
NbrFiles = sum(record.MhDraws(:,2)); NbrFiles = sum(record.MhDraws(:,2));
NbrBlocks = record.Nblck; NbrBlocks = record.Nblck;

29
matlab/GetPosteriorParametersStatistics.m
View File

@ -16,7 +16,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None. % None.
% Copyright © 2006-2018 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -38,19 +38,17 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
%end %end
TeX = options_.TeX; TeX = options_.TeX;
nblck = options_.mh_nblck;
nvx = estim_params_.nvx; nvx = estim_params_.nvx;
nvn = estim_params_.nvn; nvn = estim_params_.nvn;
ncx = estim_params_.ncx; ncx = estim_params_.ncx;
ncn = estim_params_.ncn; ncn = estim_params_.ncn;
np = estim_params_.np ; np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
OutputFolder = CheckPath('Output',M_.dname); OutputFolder = CheckPath('Output',M_.dname);
FileName = M_.fname; FileName = M_.fname;
load_last_mh_history_file(MetropolisFolder,FileName); record=load_last_mh_history_file(MetropolisFolder,FileName);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
@ -58,6 +56,7 @@ TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
mh_nblck = size(record.LastParameters,1);
clear record; clear record;
header_width = row_header_width(M_, estim_params_, bayestopt_); header_width = row_header_width(M_, estim_params_, bayestopt_);
@ -76,7 +75,7 @@ catch
disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean)) disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean))
end end
num_draws=NumberOfDraws*options_.mh_nblck; num_draws=NumberOfDraws*mh_nblck;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws); hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
if hpd_draws<2 if hpd_draws<2
@ -97,7 +96,7 @@ if np
ip = nvx+nvn+ncx+ncn+1; ip = nvx+nvn+ncx+ncn+1;
for i=1:np for i=1:np
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) 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); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
name = bayestopt_.name{ip}; name = bayestopt_.name{ip};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density); 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}; name = bayestopt_.name{ip};
[post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type); [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
catch 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); [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
name = bayestopt_.name{ip}; name = bayestopt_.name{ip};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density); oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@ -140,7 +139,7 @@ if nvx
disp(tit2) disp(tit2)
for i=1:nvx for i=1:nvx
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) 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); [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); k = estim_params_.var_exo(i,1);
name = M_.exo_names{k}; name = M_.exo_names{k};
@ -152,7 +151,7 @@ if nvx
name = M_.exo_names{k}; name = M_.exo_names{k};
[post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type); [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
catch 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] = ... [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws, 1, options_.mh_conf_sig); posterior_moments(Draws, 1, options_.mh_conf_sig);
k = estim_params_.var_exo(i,1); k = estim_params_.var_exo(i,1);
@ -184,7 +183,7 @@ if nvn
ip = nvx+1; ip = nvx+1;
for i=1:nvn for i=1:nvn
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) 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); [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)}; 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); 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)}; name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
[post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type); [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
catch 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); [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)}; 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); 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; ip = nvx+nvn+1;
for i=1:ncx for i=1:ncx
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) 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); [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); k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2); k2 = estim_params_.corrx(i,2);
@ -240,7 +239,7 @@ if ncx
NAME = sprintf('%s_%s', M_.exo_names{k1}, M_.exo_names{k2}); NAME = sprintf('%s_%s', M_.exo_names{k1}, M_.exo_names{k2});
[post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type); [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
catch 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); [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); k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2); k2 = estim_params_.corrx(i,2);
@ -273,7 +272,7 @@ if ncn
ip = nvx+nvn+ncx+1; ip = nvx+nvn+ncx+1;
for i=1:ncn for i=1:ncn
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) 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); [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); k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2); k2 = estim_params_.corrn(i,2);
@ -288,7 +287,7 @@ if ncn
NAME = sprintf('%s_%s', M_.endo_names{k1}, M_.endo_names{k2}); NAME = sprintf('%s_%s', M_.endo_names{k1}, M_.endo_names{k2});
[post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type); [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
catch 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); [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); k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2); k2 = estim_params_.corrn(i,2);

17
matlab/compute_mh_covariance_matrix.m
View File

@ -15,7 +15,7 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None. % None.
% Copyright © 2006-2017 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ options_ estim_params_ bayestopt_ global M_ bayestopt_
n = estim_params_.np + ...
estim_params_.nvn+ ...
estim_params_.ncx+ ...
estim_params_.ncn+ ...
estim_params_.nvx;
nblck = options_.mh_nblck;
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname; ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName]; BaseName = [MetropolisFolder filesep ModelName];
load_last_mh_history_file(MetropolisFolder, ModelName); record=load_last_mh_history_file(MetropolisFolder, ModelName);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
[nblck, n] = size(record.LastParameters);
posterior_kernel_at_the_mode = -Inf; posterior_kernel_at_the_mode = -Inf;
posterior_mean = zeros(n,1); posterior_mean = zeros(n,1);
posterior_mode = NaN(n,1); posterior_mode = NaN(n,1);

21
matlab/convergence_diagnostics/McMCDiagnostics.m
View File

@ -16,7 +16,7 @@ function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
% PARALLEL CONTEXT % PARALLEL CONTEXT
% See the comment in posterior_sampler.m funtion. % See the comment in posterior_sampler.m funtion.
% Copyright © 2005-2018 Dynare Team % Copyright © 2005-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -38,25 +38,12 @@ MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname; ModelName = M_.fname;
TeX = options_.TeX; 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; NumberOfMcFilesPerBlock = record.LastFileNumber;
[nblck, npar] = size(record.LastParameters);
npardisp = options_.convergence.brooksgelman.plotrows; 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. % check if all the mh files are available.
issue_an_error_message = 0; issue_an_error_message = 0;
for b = 1:nblck 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_); par_name_temp = get_the_name(jj, options_.TeX, M_, estim_params_, options_);
param_name = vertcat(param_name, par_name_temp); param_name = vertcat(param_name, par_name_temp);
end end
Draws = GetAllPosteriorDraws(jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); Draws = GetAllPosteriorDraws(jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
Draws = reshape(Draws, [NumberOfDraws nblck]); Draws = reshape(Draws, [NumberOfDraws nblck]);
Nc = min(1000, NumberOfDraws/2); Nc = min(1000, NumberOfDraws/2);
for ll = 1:nblck for ll = 1:nblck

4
matlab/generate_trace_plots.m
View File

@ -11,7 +11,7 @@ function generate_trace_plots(chain_number)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2016-2018 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -33,7 +33,7 @@ global M_ options_ estim_params_
% Get informations about the posterior draws: % Get informations about the posterior draws:
MetropolisFolder = CheckPath('metropolis', M_.dname); 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 if chain_number>record.Nblck
error('generate_trace_plots:: chain number is bigger than existing number of chains') error('generate_trace_plots:: chain number is bigger than existing number of chains')
end end

4
matlab/get_name_of_the_last_mh_file.m
View File

@ -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 % file. Otherwise info is equal to zero (a likely reason for this is
% that the mcmc simulations were not completed). % that the mcmc simulations were not completed).
% Copyright © 2008-2017 Dynare Team % Copyright © 2008-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -35,7 +35,7 @@ MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname; ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName]; BaseName = [MetropolisFolder filesep ModelName];
load_last_mh_history_file(MetropolisFolder, ModelName); record=load_last_mh_history_file(MetropolisFolder, ModelName);
mh_number = record.LastFileNumber ; mh_number = record.LastFileNumber ;
bk_number = record.Nblck ; bk_number = record.Nblck ;

9
matlab/marginal_density.m
View File

@ -15,7 +15,7 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2005-2018 Dynare Team % Copyright © 2005-2023 Dynare Team
% %
% This file is part of Dynare. % 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 <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
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); MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname; ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName]; 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; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine; FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;

4
matlab/metropolis_draw.m
View File

@ -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 % 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. % This file is part of Dynare.
% %
@ -59,7 +59,7 @@ if init
FileName = M_.fname; FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName]; BaseName = [MetropolisFolder filesep FileName];
%load mh_history-file with info on what to load %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; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));

8
matlab/mh_autocorrelation_function.m
View File

@ -18,7 +18,7 @@ function mh_autocorrelation_function(options_,M_,estim_params_,type,blck,name1,n
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% Copyright © 2003-2017 Dynare Team % Copyright © 2003-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -48,17 +48,19 @@ end
% Get informations about the posterior draws: % Get informations about the posterior draws:
MetropolisFolder = CheckPath('metropolis',M_.dname); 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; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine; FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
nblck = size(record.LastParameters,1);
clear record; clear record;
% Get all the posterior draws: % 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: % Compute the autocorrelation function:
[autocov,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size); [autocov,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size);

4
matlab/posterior_sampler.m
View File

@ -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 % Then the comments write here can be used for all the other pairs of
% parallel functions and also for management functions. % parallel functions and also for management functions.
% Copyright © 2006-2017 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % 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); InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
% Load last mh history file % 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!!! % Only for test parallel results!!!

4
matlab/posterior_sampler_initialization.m
View File

@ -340,7 +340,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
elseif 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). % Here we consider previous mh files (previous mh did not crash).
disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...') 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 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.') error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.')
end end
@ -400,7 +400,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
elseif options_.mh_recover elseif options_.mh_recover
% The previous metropolis-hastings crashed before the end! I try to recover the saved draws... % The previous metropolis-hastings crashed before the end! I try to recover the saved draws...
disp('Estimation::mcmc: Recover mode!') 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. NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains.
options_.mh_nblck = NumberOfBlocks; options_.mh_nblck = NumberOfBlocks;

9
matlab/prior_posterior_statistics.m
View File

@ -19,7 +19,7 @@ function prior_posterior_statistics(type,dataset,dataset_info)
% See the comments in the posterior_sampler.m funtion. % 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. % This file is part of Dynare.
% %
@ -209,7 +209,8 @@ localVars.DirectoryName = DirectoryName;
if strcmpi(type,'posterior') if strcmpi(type,'posterior')
BaseName = [DirectoryName filesep M_.fname]; 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; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
@ -219,9 +220,9 @@ if strcmpi(type,'posterior')
mh_nblck = options_.mh_nblck; mh_nblck = options_.mh_nblck;
if B==NumberOfDraws*mh_nblck if B==NumberOfDraws*mh_nblck
% we load all retained MH runs ! % 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 for column=1:npar
x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
end end
else else
logpost=NaN(B,1); logpost=NaN(B,1);

2
matlab/selec_posterior_draws.m
View File

@ -67,7 +67,7 @@ ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName]; BaseName = [MetropolisFolder filesep ModelName];
% Get informations about the mcmc: % Get informations about the mcmc:
load_last_mh_history_file(MetropolisFolder, ModelName); record=load_last_mh_history_file(MetropolisFolder, ModelName);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));

5
matlab/trace_plot.m
View File

@ -53,16 +53,17 @@ end
% Get informations about the posterior draws: % Get informations about the posterior draws:
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
load_last_mh_history_file(MetropolisFolder, M_.fname); record=load_last_mh_history_file(MetropolisFolder, M_.fname);
FirstMhFile = 1; FirstMhFile = 1;
FirstLine = 1; FirstLine = 1;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
[mh_nblck] = size(record.LastParameters,2);
clear record; clear record;
% Get all the posterior draws: % 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: % Plot the posterior draws: