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..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,10 +31,10 @@ 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];
-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..b7713b66d 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -1,22 +1,23 @@
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.
-% Copyright © 2006-2018 Dynare Team
+% Copyright © 2006-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -38,26 +39,24 @@ 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;
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_);
@@ -69,14 +68,12 @@ 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*options_.mh_nblck;
+num_draws=NumberOfDraws*mh_nblck;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
if hpd_draws<2
@@ -97,7 +94,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 +103,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 +137,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 +149,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 +181,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 +190,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 +220,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 +237,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 +270,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 +285,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..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]
@@ -16,7 +17,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 +39,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
@@ -71,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));
@@ -81,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_);
@@ -91,7 +79,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
@@ -120,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
@@ -230,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.')
@@ -269,7 +253,6 @@ if isnumeric(options_.parallel)
clear fout
% Parallel execution!
else
- ModelName = ModelName;
if ~isempty(M_.bvar)
ModelName = [ModelName '_bvar'];
end
@@ -404,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/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..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
@@ -11,7 +12,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.
%
@@ -28,14 +29,13 @@ 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);
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..48014ba61 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -1,21 +1,21 @@
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
-% Copyright © 2005-2018 Dynare Team
+% Copyright © 2005-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -33,20 +33,18 @@ 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;
+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 533d79ad3..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:
@@ -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,11 +59,9 @@ 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));
- 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/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..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.
@@ -37,7 +36,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 +61,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..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);
@@ -340,7 +333,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 +393,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..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
@@ -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.
%
@@ -208,20 +208,19 @@ localVars.ifil=ifil;
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));
- LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
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);
@@ -317,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
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: