diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index de04e839850365b5f1b7b9df97d4f62672a35564..93f5e2335f432b92331d6772e5961719f49adaf0 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -733,7 +733,7 @@ if do_bayesian_estimation_mcmc
oo_.mom.posterior.metropolis = oo_load_mh.oo_.mom.posterior.metropolis;
end
end
- [error_flag,~,options_mom_]= metropolis_draw(1,options_mom_,estim_params_,M_);
+ [options_mom_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_mom_); %check whether number of sub_draws is feasible
if ~(~isempty(options_mom_.sub_draws) && options_mom_.sub_draws==0)
% THIS IS PROBABLY NOT USEFUL HERE AND CAN BE REMOVED (PREPROCESSOR: REMOVE bayesian_irf, moments_varendo)
%if options_mom_.bayesian_irf
diff --git a/matlab/estimation/CutSample.m b/matlab/estimation/CutSample.m
index 462aa070a4c3906dd6537c2063f0c3902796cb28..4cbbbc887712e44db34d183ce6e1240ccdeeb270 100644
--- a/matlab/estimation/CutSample.m
+++ b/matlab/estimation/CutSample.m
@@ -46,7 +46,7 @@ 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' ]);
-if ~length(mh_files)
+if isempty(mh_files)
error('%s: I can''t find MH file to load here!',dispString)
end
diff --git a/matlab/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index 9a4b67b3e44d31fb34931ce82c1ea4a5d84c9596..fabf390af4934d131ec41357988d4d042396a6cc 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -37,24 +37,28 @@ function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstMhFil
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if ishssmc(options_)
- % Load draws from the posterior distribution
- pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
- posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
- if column==0
- draws = posterior.tlogpostkernel;
- else
- draws = transpose(posterior.particles(column,:));
- end
-elseif isdime(options_)
- posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
- tune = posterior.tune;
- chains = posterior.chains(end-tune:end,:,:);
- if column>0
- chains = reshape(chains, [], size(chains, 3));
- draws = chains(:,column);
+if issmc(options_)
+ if ishssmc(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
+ posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
+ if column==0
+ draws = posterior.tlogpostkernel;
+ else
+ draws = transpose(posterior.particles(column,:));
+ end
+ elseif isdime(options_)
+ posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
+ tune = posterior.tune;
+ chains = posterior.chains(end-tune:end,:,:);
+ if column>0
+ chains = reshape(chains, [], size(chains, 3));
+ draws = chains(:,column);
+ else
+ draws = posterior.lprobs;
+ end
else
- draws = posterior.lprobs;
+ error('GetAllPosteriorDraws:: case should not happen. Please contact the developers')
end
else
iline = FirstLine;
diff --git a/matlab/estimation/GetOneDraw.m b/matlab/estimation/GetOneDraw.m
deleted file mode 100644
index 8c92a84ba70dac74984995a35d658ca067dd4c54..0000000000000000000000000000000000000000
--- a/matlab/estimation/GetOneDraw.m
+++ /dev/null
@@ -1,47 +0,0 @@
-function [xparams, logpost] = GetOneDraw(distribution, M_, estim_params_, oo_, options_, bayestopt_)
-
-% draws one parameter vector and its posterior from MCMC or the prior
-%
-% INPUTS
-% type: [string] 'posterior': draw from MCMC draws
-% 'prior': draw from prior
-% M_ [structure] Definition of the model
-% estim_params_ [structure] characterizing parameters to be estimated
-% oo_ [structure] Storage of results
-% options_ [structure] Options
-% bayestopt_ [structure] describing the priors
-%
-% OUTPUTS
-% xparams: vector of estimated parameters (drawn from posterior or prior distribution)
-% logpost: log of the posterior density of this parameter vector
-%
-% SPECIAL REQUIREMENTS
-% none
-
-% Copyright © 2005-2023 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-if isprior(distribution)
- xparams = distribution.draw();
- if nargout>1
- logpost = evaluate_posterior_kernel(xparams, M_, estim_params_, oo_, options_, bayestopt_);
- end
-elseif ischar(distribution) && strcmpi(distribution, 'posterior')
- [xparams, logpost] = metropolis_draw(0);
-else
- error('GetOneDraw:: Wrong inputs.')
-end
diff --git a/matlab/estimation/GetPosteriorMeanVariance.m b/matlab/estimation/GetPosteriorMeanVariance.m
index 6cb5a7f4ab5eaaec2ee0e70ba183e93b7b401115..9a1c1205129f6cc24b7a0745b15e92ec94db039b 100644
--- a/matlab/estimation/GetPosteriorMeanVariance.m
+++ b/matlab/estimation/GetPosteriorMeanVariance.m
@@ -1,7 +1,6 @@
function [mean, variance] = GetPosteriorMeanVariance(options_, M_)
% [mean,variance] = GetPosteriorMeanVariance(options_, M_)
% Computes the posterior mean and variance
-% (+updates of oo_ & TeX output).
%
% INPUTS
% - options_ [struct] Dynare's options.
@@ -28,14 +27,25 @@ function [mean, variance] = GetPosteriorMeanVariance(options_, M_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if ishssmc(options_)
- % Load draws from the posterior distribution
- pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
- posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
- % Compute the posterior mean
- mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
- % Compute the posterior covariance
- variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
+if issmc(options_)
+ if ishssmc(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+ posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ % Compute the posterior mean
+ mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
+ % Compute the posterior covariance
+ variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
+ elseif isdime(options_)
+ posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+ tune = posterior.tune;
+ chains = transpose(reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3)));
+ mean = sum(chains, 2)/size(chains,1);
+ % Compute the posterior covariance
+ variance = (chains-mean)*(chains-mean)'/size(chains,1);
+ else
+ error('GetPosteriorMeanVariance:: case should not happen. Please contact the developers')
+ end
else
MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname;
diff --git a/matlab/estimation/PosteriorIRF.m b/matlab/estimation/PosteriorIRF.m
index d53bc344d17795840f4a8fe34d2bf9d7dc5c3e0f..d6cf33ea5be1c63cf6d2b153a24958f7536058d9 100644
--- a/matlab/estimation/PosteriorIRF.m
+++ b/matlab/estimation/PosteriorIRF.m
@@ -156,10 +156,7 @@ localVars.npar = npar;
localVars.type=type;
if strcmpi(type,'posterior')
- while b<B
- b = b + 1;
- x(b,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
- end
+ x=get_posterior_subsample(M_,options_,B);
end
if ~strcmpi(type,'prior')
diff --git a/matlab/estimation/check_posterior_analysis_data.m b/matlab/estimation/check_posterior_analysis_data.m
index c8a8c71f2bb47d5c155a16305d5b435be6090a9d..90c836afc9b759f5311fcb7ccca5aa34a0f203ab 100644
--- a/matlab/estimation/check_posterior_analysis_data.m
+++ b/matlab/estimation/check_posterior_analysis_data.m
@@ -1,11 +1,12 @@
-function [info,description] = check_posterior_analysis_data(type,M_)
-% function [info,description] = check_posterior_analysis_data(type,M_)
+function [info,description] = check_posterior_analysis_data(type,M_,options_)
+% function [info,description] = check_posterior_analysis_data(type,M_,options_)
% Checks the status of posterior analysis and in particular if files need to be
% created or updated; called by posterior_analysis.m
%
% Inputs:
% type [string] name of the posterior moment considered
% M_ [structure] Dynare model structure
+% options_ [structure] Dynare options structure
%
% Outputs:
% info [scalar] return code
@@ -17,7 +18,7 @@ function [info,description] = check_posterior_analysis_data(type,M_)
% info = 6; % Ok (nothing to do ;-)
% description [string] Message corresponding to info
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -39,7 +40,17 @@ if nargout>1
description = '';
end
-[MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
+if ~issmc(options_)
+ [MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
+else
+ if ishssmc(options_)
+ [MetropolisFolder, info] = CheckPath('hssmc',M_.dname);
+ elseif isdime(options_)
+ [MetropolisFolder, info] = CheckPath('dime',M_.dname);
+ else
+ error('check_posterior_analysis_data:: case should not happen. Please contact the developers')
+ end
+end
% Get informations about mcmc files.
if info
@@ -47,9 +58,6 @@ if info
return
end
-mhname = get_name_of_the_last_mh_file(M_);
-mhdate = get_date_of_a_file([MetropolisFolder filesep mhname]);
-
% Get informations about _posterior_draws files.
drawsinfo = dir([ MetropolisFolder filesep M_.fname '_posterior_draws*.mat']);
if isempty(drawsinfo)
@@ -59,7 +67,19 @@ if isempty(drawsinfo)
end
return
else
- number_of_last_posterior_draws_file = length(drawsinfo);
+ mhname = get_name_of_the_last_mh_file(M_);
+ if ~issmc(options_)
+ mhdate = get_date_of_a_file([MetropolisFolder filesep mhname]);
+ number_of_last_posterior_draws_file = length(drawsinfo);
+ else
+ if ishssmc(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+ mhdate = get_date_of_a_file(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ elseif isdime(options_)
+ mhdate = get_date_of_a_file(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+ end
+ end
pddate = get_date_of_a_file([ MetropolisFolder filesep M_.fname '_posterior_draws' int2str(number_of_last_posterior_draws_file) '.mat']);
if pddate<mhdate
info = 2; % _posterior_draws files have to be updated.
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 8d1efefbc4b2c4cc48206352a3923d96b2522f4e..23436980c5d829d76763792e7aa57b3c66f1109a 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -483,9 +483,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
end
% Store posterior mean in a vector and posterior variance in
% a matrix
- if ~isdime(options_)
- [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] = GetPosteriorMeanVariance(options_, M_);
- end
+ [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] = GetPosteriorMeanVariance(options_, M_);
elseif options_.load_mh_file && options_.load_results_after_load_mh
% load fields from previous MCMC run stored in results-file
field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
@@ -506,28 +504,22 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
oo_.posterior.metropolis=oo_load_mh.oo_.posterior.metropolis;
end
end
- if ~issmc(options_)
- [error_flag, ~, options_]= metropolis_draw(1, options_, estim_params_, M_);
- else
- error_flag=false;
- end
+ [options_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_); %check whether number is feasible
if ~(~isempty(options_.sub_draws) && options_.sub_draws==0)
if options_.bayesian_irf
- if ~issmc(options_)
- if error_flag
- error('%s: I cannot compute the posterior IRFs!',dispString)
- end
- oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
- else
- fprintf('%s: SMC does not yet support the bayesian_irf option. Skipping computation.\n',dispString);
+ if error_flag
+ error('%s: I cannot compute the posterior IRFs!',dispString)
end
+ oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
end
if options_.moments_varendo
- if ~issmc(options_)
- if error_flag
- error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
- end
- if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results
+ if error_flag
+ error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
+ end
+ if options_.load_mh_file
+ if issmc(options_)
+ error('%s: SMC does not yet support the load_mh_file option.\n',dispString);
+ elseif options_.mh_replic==0 %user wants to recompute results for standard MCMC
[MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
if ~info
generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
@@ -546,10 +538,8 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
end
end
end
- oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
- else
- fprintf('%s: SMC does not yet support the moments_varendo option. Skipping computation.\n',dispString);
end
+ oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
end
if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
if ~ishssmc(options_) && ~isdime(options_)
diff --git a/matlab/estimation/execute_prior_posterior_function.m b/matlab/estimation/execute_prior_posterior_function.m
index f23d4664babdb2bb206837bc6671e1a6abeb5253..898cf63005952007353f40404ec49f83ee749613 100644
--- a/matlab/estimation/execute_prior_posterior_function.m
+++ b/matlab/estimation/execute_prior_posterior_function.m
@@ -17,7 +17,7 @@ function oo_=execute_prior_posterior_function(posterior_function_name,M_,options
% OUTPUTS
% oo_ [structure] Matlab/Octave structure gathering the results (initialized by dynare, see @ref{oo_}).
-% Copyright © 2013-2023 Dynare Team
+% Copyright © 2013-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -53,9 +53,8 @@ if strcmpi(type,'posterior')
% Get informations about the _posterior_draws files.
% discard first mh_drop percent of the draws:
CutSample(M_, options_, 'prior_posterior_function');
- % initialize metropolis draws
options_.sub_draws = n_draws; % set draws for sampling; changed value is not returned to base workspace
- [error_flag, ~, options_] = metropolis_draw(1, options_, estim_params_, M_);
+ [options_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_); %check whether number is feasible
if error_flag
error('EXECUTE_POSTERIOR_FUNCTION: The draws could not be initialized')
end
@@ -69,7 +68,7 @@ elseif strcmpi(type,'prior')
error('The prior distributions are not properly set up.')
end
end
- if exist([M_.fname '_prior_restrictions.m'])
+ if exist([M_.fname '_prior_restrictions.m'],"file")
warning('prior_function currently does not support endogenous prior restrictions. They will be ignored. Consider using a posterior_function with nobs=1.')
end
Prior = dprior(bayestopt_, options_.prior_trunc);
@@ -77,13 +76,11 @@ else
error('EXECUTE_POSTERIOR_FUNCTION: Unknown type!')
end
+
if strcmpi(type, 'prior')
parameter_mat = Prior.draws(n_draws);
else
- parameter_mat = NaN(length(bayestopt_.p6), n_draws);
- for i = 1:n_draws
- parameter_mat(:,i) = GetOneDraw(type, M_, estim_params_, oo_, options_, bayestopt_);
- end
+ parameter_mat=get_posterior_subsample(M_,options_,n_draws)';
end
% Get output size
diff --git a/matlab/estimation/get_posterior_subsample.m b/matlab/estimation/get_posterior_subsample.m
new file mode 100644
index 0000000000000000000000000000000000000000..5df136f5a2f3f9bc71891037b91f00cadd02b58f
--- /dev/null
+++ b/matlab/estimation/get_posterior_subsample.m
@@ -0,0 +1,105 @@
+function [pdraws, log_posterior]=get_posterior_subsample(M_,options_,SampleSize)
+% function [pdraws, log_posterior]=get_posterior_subsample(M_,options_,SampleSize)
+% Checks the status of posterior analysis and in particular if files need to be
+% created or updated; called by posterior_analysis.m
+%
+% Inputs:
+% M_ [structure] Dynare model structure
+% options_ [structure] Dynare options structure
+% SampleSize [integer] number of desired subdraws
+%
+% Outputs:
+% pdraws [double] SampleSize by npar matrix of parameter draws
+% log_posterior [double] SampleSize by 1 vector of log posterior
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+if ~issmc(options_)
+ % Get informations about the mcmc:
+ record=load_last_mh_history_file([M_.dname filesep 'metropolis'], M_.fname);
+ npar=size(record.InitialParameters,2);
+ FirstMhFile = record.KeepedDraws.FirstMhFile;
+ FirstLine = record.KeepedDraws.FirstLine;
+ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+ NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+ MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+
+ % Randomly select draws in the posterior distribution:
+ SampleAddress = zeros(SampleSize,4);
+ for iter = 1:SampleSize
+ ChainNumber = ceil(rand*options_.mh_nblck);
+ DrawNumber = ceil(rand*NumberOfDraws);
+ SampleAddress(iter,1) = DrawNumber;
+ SampleAddress(iter,2) = ChainNumber;
+ if DrawNumber <= MAX_nruns-FirstLine+1 % in first file where FirstLine may be bigger than 1
+ MhFileNumber = FirstMhFile;
+ MhLineNumber = FirstLine+DrawNumber-1;
+ else
+ DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
+ MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
+ MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
+ end
+ SampleAddress(iter,3) = MhFileNumber;
+ SampleAddress(iter,4) = MhLineNumber;
+ end
+ SampleAddress = sortrows(SampleAddress,[3 2]);
+
+ % Selected draws in the posterior distribution, and if drsize>0
+ % reduced form solutions, are saved on disk.
+ pdraws = NaN(SampleSize,npar);
+ if nargout>1
+ log_posterior = NaN(SampleSize,1);
+ end
+ old_mhfile = 0;
+ old_mhblck = 0;
+ for iter = 1:SampleSize
+ mhfile = SampleAddress(iter,3);
+ mhblck = SampleAddress(iter,2);
+ if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
+ temp=load([M_.dname filesep 'metropolis' filesep M_.fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2','logpo2');
+ end
+ pdraws(iter,:) = temp.x2(SampleAddress(iter,4),:);
+ if nargout>1
+ log_posterior(iter,1) = temp.logpo2(SampleAddress(iter,4));
+ end
+ old_mhfile = mhfile;
+ old_mhblck = mhblck;
+ end
+else
+ if ishssmc(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+ posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ NumberOfDraws=size(posterior.particles,2);
+ DrawNumbers = ceil(rand(SampleSize,1)*NumberOfDraws);
+ pdraws = transpose(posterior.particles(:,DrawNumbers));
+ if nargout>1
+ log_posterior = posterior.tlogpostkernel(DrawNumbers);
+ end
+ elseif isdime(options_)
+ posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+ tune = posterior.tune;
+ chains = reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3));
+ NumberOfDraws=size(chains,1);
+ DrawNumbers = ceil(rand(SampleSize,1)*NumberOfDraws);
+ pdraws = chains(DrawNumbers,:);
+ if nargout>1
+ log_posterior = posterior.lprobs(DrawNumbers);
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/estimation/metropolis_draw.m b/matlab/estimation/metropolis_draw.m
deleted file mode 100644
index 505d1dbeb6b5deaa7964a366b882f5fabf2cd2a7..0000000000000000000000000000000000000000
--- a/matlab/estimation/metropolis_draw.m
+++ /dev/null
@@ -1,99 +0,0 @@
-function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
-% function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
-% Builds draws from metropolis
-%
-% INPUTS:
-% init: scalar equal to
-% 1: first call to store the required information
-% on files, lines, and chains to be read
-% in persistent variables to make them available
-% for future calls
-% 0: load a parameter draw
-% Additional Inputs required for initialization
-% options_ [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
-% M_ [structure] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
-% estim_params_ [structure] Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
-%
-% OUTPUTS:
-% xparams: if init==0: vector of estimated parameters
-% if init==1: error flaog
-% logpost: log of posterior density
-% options_: [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
-%
-% SPECIAL REQUIREMENTS
-%
-% Requires CutSample to be run before in order to set up mh_history-file
-
-% Copyright © 2003-2023 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns
-
-xparams = 0;
-logpost = 0;
-
-if init
- %get number of parameters
- nvx = estim_params_.nvx;
- nvn = estim_params_.nvn;
- ncx = estim_params_.ncx;
- ncn = estim_params_.ncn;
- np = estim_params_.np ;
- npar = nvx+nvn+ncx+ncn+np;
- %get path of metropolis files
- MetropolisFolder = CheckPath('metropolis',M_.dname);
- FileName = M_.fname;
- BaseName = [MetropolisFolder filesep FileName];
- %load mh_history-file with info on what to load
- record=load_last_mh_history_file(MetropolisFolder, FileName);
- FirstMhFile = record.KeepedDraws.FirstMhFile;
- FirstLine = record.KeepedDraws.FirstLine;
- 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 ?
- mh_nblck = options_.mh_nblck;
- % set sub_draws option if empty
- if isempty(options_.sub_draws)
- options_.sub_draws = min(options_.posterior_max_subsample_draws, ceil(.25*NumberOfDraws));
- else
- if options_.sub_draws>NumberOfDraws*mh_nblck
- skipline()
- disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*mh_nblck) ')!'])
- disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
- skipline()
- xparams = 1; % xparams is interpreted as an error flag
- end
- end
- return
-else %not initialization, return one draw
- %get random draw from random chain
- ChainNumber = ceil(rand*mh_nblck);
- DrawNumber = ceil(rand*NumberOfDraws);
-
- if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line
- MhFilNumber = FirstMhFile;
- MhLine = FirstLine+DrawNumber-1;
- else %draw in other file
- DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
- MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
- MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
- end
- %load parameters and posterior
- load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
- xparams = x2(MhLine,:);
- logpost= logpo2(MhLine);
-end
\ No newline at end of file
diff --git a/matlab/estimation/posterior_analysis.m b/matlab/estimation/posterior_analysis.m
index 72f4713da3bdd09b24a987195124ba3f0566a261..a412e4277d661eae64b0aa254e433fb6d608b25d 100644
--- a/matlab/estimation/posterior_analysis.m
+++ b/matlab/estimation/posterior_analysis.m
@@ -17,7 +17,7 @@ function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_,estim_para
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-info = check_posterior_analysis_data(type,M_);
+info = check_posterior_analysis_data(type,M_,options_);
SampleSize = options_.sub_draws;
switch info
case 0
diff --git a/matlab/estimation/prior_posterior_statistics.m b/matlab/estimation/prior_posterior_statistics.m
index 6971f9465660c4f7e6b624b54786b382638b83bd..7ff142419183a00f86bae58a5fe04eed88c22427 100644
--- a/matlab/estimation/prior_posterior_statistics.m
+++ b/matlab/estimation/prior_posterior_statistics.m
@@ -217,10 +217,7 @@ if strcmpi(type,'posterior')
x(:,column) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
end
else
- logpost=NaN(B,1);
- for b=1:B
- [x(b,:), logpost(b)] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
- end
+ [x, logpost]=get_posterior_subsample(M_,options_,B);
end
localVars.logpost=logpost;
end
diff --git a/matlab/estimation/selec_posterior_draws.m b/matlab/estimation/selec_posterior_draws.m
index 624c6e54f6e648553ec6d882860faba041cc983e..e001e1bc1abf5a54fe2deb6ea0d04ecaeb1d037a 100644
--- a/matlab/estimation/selec_posterior_draws.m
+++ b/matlab/estimation/selec_posterior_draws.m
@@ -2,8 +2,7 @@ function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,
% Selects a sample of draws from the posterior distribution and if nargin>1
% saves the draws in _pdraws mat files (metropolis folder). If drsize>0
% the dr structure, associated to the parameters, is also saved in _pdraws.
-% This routine is more efficient than metropolis_draw.m because here an
-% _mh file cannot be opened twice.
+% This routine assures an _mh file cannot be opened twice.
%
% INPUTS
% o M_ [structure] Matlab's structure describing the model
@@ -55,7 +54,7 @@ switch nargin
case 8
info = 0;
case 9
- MAX_mega_bytes = 10;% Should be an option...
+ MAX_mega_bytes = 100;% Should be an option...
if drsize>0
info=2;
else
@@ -66,100 +65,55 @@ switch nargin
error('selec_posterior_draws:: Unexpected number of input arguments!')
end
-MetropolisFolder = CheckPath('metropolis',M_.dname);
-ModelName = M_.fname;
-BaseName = [MetropolisFolder filesep ModelName];
-
-% Get informations about the mcmc:
-record=load_last_mh_history_file(MetropolisFolder, ModelName);
-FirstMhFile = record.KeepedDraws.FirstMhFile;
-FirstLine = record.KeepedDraws.FirstLine;
-TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
-mh_nblck = options_.mh_nblck;
-
-% Randomly select draws in the posterior distribution:
-SampleAddress = zeros(SampleSize,4);
-for i = 1:SampleSize
- ChainNumber = ceil(rand*mh_nblck);
- DrawNumber = ceil(rand*NumberOfDraws);
- SampleAddress(i,1) = DrawNumber;
- SampleAddress(i,2) = ChainNumber;
- if DrawNumber <= MAX_nruns-FirstLine+1
- MhFileNumber = FirstMhFile;
- MhLineNumber = FirstLine+DrawNumber-1;
+if ~issmc(options_)
+ MetropolisFolder = CheckPath('metropolis',M_.dname);
+else
+ if ishssmc(options_)
+ [MetropolisFolder] = CheckPath('metropolis',M_.dname);
+ elseif isdime(options_)
+ [MetropolisFolder] = CheckPath('dime',M_.dname);
else
- DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
- MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
- MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
+ error('check_posterior_analysis_data:: case should not happen. Please contact the developers')
end
- SampleAddress(i,3) = MhFileNumber;
- SampleAddress(i,4) = MhLineNumber;
end
-SampleAddress = sortrows(SampleAddress,[3 2]);
+
+ModelName = M_.fname;
+BaseName = [MetropolisFolder filesep ModelName];
+
+parameter_draws=get_posterior_subsample(M_,options_,SampleSize);
% Selected draws in the posterior distribution, and if drsize>0
% reduced form solutions, are saved on disk.
if info
%delete old stale files before creating new ones
delete_stale_file([BaseName '_posterior_draws*.mat'])
- if SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file.
- pdraws = cell(SampleSize,info);
- old_mhfile = 0;
- old_mhblck = 0;
- for i = 1:SampleSize
- mhfile = SampleAddress(i,3);
- mhblck = SampleAddress(i,2);
- if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
- load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
- end
- pdraws(i,1) = {x2(SampleAddress(i,4),:)};
- if info==2
- M_ = set_parameters_locally(M_,pdraws{i,1});
- [dr,~,M_.params] =compute_decision_rules(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
- pdraws(i,2) = { dr };
- end
- old_mhfile = mhfile;
- old_mhblck = mhblck;
- end
- clear('x2')
- save([BaseName '_posterior_draws1.mat'],'pdraws','estim_params_')
- else% The posterior draws are saved in xx files.
- NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
- NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
- NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
- linee = 0;
- fnum = 1;
+ NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
+ NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
+ NumberOfLinesLastFile = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
+ linee = 0;
+ fnum = 1;
+ if fnum < NumberOfFiles
pdraws = cell(NumberOfDrawsPerFile,info);
- old_mhfile = 0;
- old_mhblck = 0;
- for i=1:SampleSize
- linee = linee+1;
- mhfile = SampleAddress(i,3);
- mhblck = SampleAddress(i,2);
- if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
- load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
- end
- pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
- if info==2
- M_ = set_parameters_locally(M_,pdraws{linee,1});
- [dr,~,M_.params] = compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
- pdraws(linee,2) = { dr };
- end
- old_mhfile = mhfile;
- old_mhblck = mhblck;
- if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
- linee = 0;
- save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
- fnum = fnum+1;
- if fnum < NumberOfFiles
- pdraws = cell(NumberOfDrawsPerFile,info);
- else
- pdraws = cell(NumberOfLines,info);
- end
+ else
+ pdraws = cell(NumberOfLinesLastFile,info);
+ end
+ for i=1:SampleSize
+ linee = linee+1;
+ pdraws(i,1) = {parameter_draws(i,:)};
+ if info==2
+ M_ = set_parameters_locally(M_,pdraws{linee,1});
+ [dr,~,M_.params] = compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
+ pdraws(linee,2) = { dr };
+ end
+ if (fnum < NumberOfFiles && linee == NumberOfDrawsPerFile) || (fnum <= NumberOfFiles && i==SampleSize)
+ linee = 0;
+ save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
+ fnum = fnum+1;
+ if fnum < NumberOfFiles
+ pdraws = cell(NumberOfDrawsPerFile,info);
+ else
+ pdraws = cell(NumberOfLinesLastFile,info);
end
end
- save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
end
end
\ No newline at end of file
diff --git a/matlab/estimation/set_number_of_subdraws.m b/matlab/estimation/set_number_of_subdraws.m
new file mode 100644
index 0000000000000000000000000000000000000000..027303d16738245daf60b51c9638b00eafaf04e3
--- /dev/null
+++ b/matlab/estimation/set_number_of_subdraws.m
@@ -0,0 +1,78 @@
+function [sub_draws, error_flag]=set_number_of_subdraws(M_,options_)
+% function [sub_draws, error_flag]=set_number_of_subdraws(M_,options_)
+% Set option field sub_draws based on size of available sample
+% Inputs:
+% M_ [structure] Dynare model structure
+% options_ [structure] Dynare options structure
+%
+% Outputs:
+% sub_draws [scalar] number of sub-draws
+% error_flag [boolean] error indicate of insufficient draws are
+% available
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+
+error_flag=false;
+if ~issmc(options_)
+ record=load_last_mh_history_file([M_.dname filesep 'metropolis'], M_.fname);
+ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+ NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+ if isempty(options_.sub_draws)
+ sub_draws = min(options_.posterior_max_subsample_draws, ceil(.25*NumberOfDraws));
+ else
+ if options_.sub_draws>NumberOfDraws*record.Nblck
+ skipline()
+ disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*options_.mh_nblck) ')!'])
+ disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
+ skipline()
+ error_flag = true;
+ sub_draws=options_.sub_draws; %pass through
+ else
+ sub_draws=options_.sub_draws; %pass through
+ end
+ end
+else
+ if ishssmc(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+ posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ NumberOfDraws = size(posterior.particles,2);
+ elseif isdime(options_)
+ posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+ tune = posterior.tune;
+ chains = reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3));
+ NumberOfDraws=size(chains,1);
+ else
+ error('set_number_of_subdraws:: case should not happen. Please contact the developers')
+ end
+ if isempty(options_.sub_draws)
+ sub_draws = min(options_.posterior_max_subsample_draws, ceil(.25*NumberOfDraws));
+ else
+ if options_.sub_draws>NumberOfDraws
+ skipline()
+ disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws) ')!'])
+ disp('You can either change the value of sub_draws or run another MCMC.')
+ skipline()
+ error_flag = true;
+ sub_draws=options_.sub_draws; %pass through
+ else
+ sub_draws=options_.sub_draws; %pass through
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/list_of_functions_to_be_cleared.m b/matlab/list_of_functions_to_be_cleared.m
index 90c7187b3c8c34dae908937f592ff520c492f735..dc1a5eadcf86bd551cfc23a038927d1b6f0056f3 100644
--- a/matlab/list_of_functions_to_be_cleared.m
+++ b/matlab/list_of_functions_to_be_cleared.m
@@ -1 +1 @@
-list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', '+gsa/prior_draw.m', '+identification/analysis.m', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', '+occbin/solver.m','+occbin/mkdatap_anticipated_dyn.m','+occbin/mkdatap_anticipated_2constraints_dyn.m','+occbin/match_function.m','+occbin/solve_one_constraint.m','+occbin/solve_two_constraint.m','+occbin/plot/shock_decomposition.m'};
+list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', '+gsa/prior_draw.m', '+identification/analysis.m', 'computeDLIK', 'univariate_computeDLIK', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', '+occbin/solver.m','+occbin/mkdatap_anticipated_dyn.m','+occbin/mkdatap_anticipated_2constraints_dyn.m','+occbin/match_function.m','+occbin/solve_one_constraint.m','+occbin/solve_two_constraint.m','+occbin/plot/shock_decomposition.m'};