diff --git a/doc/dynare.texi b/doc/dynare.texi
index ba6d7bcc8..4a312ee6e 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -4494,12 +4494,20 @@ algorithm. For the time being, @code{mh_replic} should be larger than
@code{1200}. Default: @code{20000}
@item sub_draws = @var{INTEGER}
-@anchor{sub_draws} number of draws from the Metropolis iterations that
-are used to compute posterior distribution of various objects (smoothed
-variable, smoothed shocks, forecast, moments, IRF). @code{sub_draws} should be smaller than
-the total number of Metropolis draws available. Default:
-@code{min(1200,0.25*Total number of draws)}
+@anchor{sub_draws} number of draws from the MCMC that are used to
+compute posterior distribution of various objects (smoothed variable,
+smoothed shocks, forecast, moments, IRF). The draws used to compute
+these posterior moments are sampled uniformly in the estimated empirical
+posterior distribution (@i{ie} draws of the MCMC). @code{sub_draws}
+should be smaller than the total number of MCMC draws available.
+Default: @code{min(posterior_max_subsample_draws,0.25*Total number of
+draws)}
+@item posterior_max_subsample_draws = @var{INTEGER}
+@anchor{posterior_max_subsample_draws} maximum number of draws from the
+MCMC used to compute posterior distribution of various objects (smoothed
+variable, smoothed shocks, forecast, moments, IRF), if not overriden by
+option @ref{sub_draws}.
@item mh_nblocks = @var{INTEGER}
@anchor{mh_nblocks} Number of parallel chains for Metropolis-Hastings algorithm. Default:
diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index 018c6ad3a..cf9737f0f 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -95,9 +95,8 @@ end
delete([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']);
delete([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']);
if strcmpi(type,'posterior')
- load_last_mh_history_file(MhDirectoryName, M_.fname);
- TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
- NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+ B = options_.sub_draws;
+ options_.B = B;
elseif strcmpi(type,'gsa')
RootDirectoryName = CheckPath('gsa',M_.dname);
if options_.opt_gsa.pprior
@@ -107,13 +106,10 @@ elseif strcmpi(type,'gsa')
end
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat istable
- NumberOfDraws=size(x,1);
- B=NumberOfDraws; options_.B = B;
+ B=size(x,1); options_.B = B;
else% type = 'prior'
- NumberOfDraws = options_.prior_draws;
-end
-if ~strcmpi(type,'gsa')
- B = min([round(.5*NumberOfDraws),500]); options_.B = B;
+ B = options_.prior_draws;
+ options_.B = B;
end
try
delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 0578791d7..d925acc58 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -827,14 +827,23 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
else
load([M_.fname '_results'],'oo_');
end
- metropolis_draw(1);
+ error_flag = metropolis_draw(1);
if options_.bayesian_irf
+ if error_flag
+ error('Estimation::mcmc: I cannot compute the posterior IRFs!')
+ end
PosteriorIRF('posterior');
end
if options_.moments_varendo
+ if error_flag
+ error('Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!')
+ end
oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
end
if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
+ if error_flag
+ error('Estimation::mcmc: I cannot compute the posterior statistics!')
+ end
prior_posterior_statistics('posterior',dataset_);
end
xparam = get_posterior_parameters('mean');
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index de4741c13..04b2242a2 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -355,7 +355,6 @@ options_.dataset.xls_range = NaN;
options_.Harvey_scale_factor = 10;
options_.MaxNumberOfBytes = 1e6;
options_.MaximumNumberOfMegaBytes = 111;
-options_.PosteriorSampleSize = 1000;
options_.analytic_derivation = 0;
options_.analytic_derivation_mode = 0;
options_.bayesian_irf = 0;
@@ -399,6 +398,7 @@ options_.presample = 0;
options_.prior_trunc = 1e-10;
options_.smoother = 0;
options_.student_degrees_of_freedom = 3;
+options_.posterior_max_subsample_draws = 1200;
options_.sub_draws = [];
options_.use_mh_covariance_matrix = 0;
options_.gradient_method = 2;
diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m
index 8e1f10453..771ba0a43 100644
--- a/matlab/metropolis_draw.m
+++ b/matlab/metropolis_draw.m
@@ -32,6 +32,9 @@ function [xparams, logpost]=metropolis_draw(init)
global options_ estim_params_ M_
persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns
+xparams = 0;
+logpost = 0;
+
if init
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
@@ -39,7 +42,6 @@ if init
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
-
MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName];
@@ -52,6 +54,18 @@ if init
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
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, round(.25*NumberOfDraws));
+ else
+ if options_.sub_draws>NumberOfDraws
+ skipline()
+ disp(['Estimation::mcmc: The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws) ')!'])
+ disp('Estimation::mcmc: 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
end
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index df978c647..0382f923f 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -17,7 +17,7 @@ function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_)
% along with Dynare. If not, see .
info = check_posterior_analysis_data(type,M_);
-SampleSize = options_.PosteriorSampleSize;
+SampleSize = options_.sub_draws;
switch info
case 0
disp('check_posterior_analysis_data:: Can''t find any mcmc file!')
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index fa4ff29ad..f9192164e 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -69,21 +69,7 @@ maxlag = M_.maximum_endo_lag;
if strcmpi(type,'posterior')
DirectoryName = CheckPath('metropolis',M_.dname);
- load_last_mh_history_file(DirectoryName, M_.fname);
- 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);
- clear record;
- if ~isempty(options_.sub_draws)
- B = options_.sub_draws;
- if B > NumberOfDraws
- B = NumberOfDraws;
- end
- else
- B = min(1200, round(0.25*NumberOfDraws));
- end
+ B = options_.sub_draws;
elseif strcmpi(type,'gsa')
RootDirectoryName = CheckPath('gsa',M_.dname);
if options_.opt_gsa.pprior
@@ -99,8 +85,7 @@ elseif strcmpi(type,'gsa')
x=lpmat(istable,:);
end
clear lpmat lpmat0 istable
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
+ B = size(x,1);
elseif strcmpi(type,'prior')
DirectoryName = CheckPath('prior',M_.dname);
B = options_.prior_draws;
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index ea7d82e0d..3e95184c9 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -110,7 +110,7 @@ class ParsingDriver;
%token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS TOLF
%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LYAPUNOV
%token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT
-%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS
+%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS
%token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE
%token NAME
@@ -1633,6 +1633,7 @@ estimation_options : o_datafile
| o_geweke_interval
| o_mcmc_jumping_covariance
| o_irf_plot_threshold
+ | o_posterior_max_subsample_draws
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -2423,6 +2424,7 @@ o_subsample_name : symbol EQUAL date_expr ':' date_expr
;
o_conf_sig : CONF_SIG EQUAL non_negative_number { driver.option_num("conf_sig", $3); };
o_mh_replic : MH_REPLIC EQUAL INT_NUMBER { driver.option_num("mh_replic", $3); };
+o_posterior_max_subsample_draws : POSTERIOR_MAX_SUBSAMPLE_DRAWS EQUAL INT_NUMBER { driver.option_num("posterior_max_subsample_draws", $3); };
o_mh_drop : MH_DROP EQUAL non_negative_number { driver.option_num("mh_drop", $3); };
o_mh_jscale : MH_JSCALE EQUAL non_negative_number { driver.option_num("mh_jscale", $3); };
o_optim : OPTIM EQUAL '(' optim_options ')';
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 7343a5ad1..9b6b5ed46 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -284,6 +284,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
dsge_var {return token::DSGE_VAR;}
dsge_varlag {return token::DSGE_VARLAG;}
moments_varendo {return token::MOMENTS_VARENDO;}
+posterior_max_subsample_draws {return token::POSTERIOR_MAX_SUBSAMPLE_DRAWS;}
filtered_vars {return token::FILTERED_VARS;}
filter_step_ahead {return token::FILTER_STEP_AHEAD;}
relative_irf {return token::RELATIVE_IRF;}