Closes #567.

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
+@anchor{sub_draws}  number of  draws  from  the MCMC  that  are used  to
-are used to compute posterior distribution of various objects (smoothed
+compute posterior  distribution of  various objects  (smoothed variable,
-variable, smoothed shocks, forecast, moments, IRF). @code{sub_draws} should be smaller than
+smoothed shocks,  forecast, moments,  IRF).  The  draws used  to compute
-the total number of Metropolis draws available. Default:
+these posterior moments are sampled uniformly in the estimated empirical
-@code{min(1200,0.25*Total number of draws)}
+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:

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);
+    B = options_.sub_draws;
-    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    options_.B = B;
-    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
 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=size(x,1); options_.B = B;
-    B=NumberOfDraws; options_.B = B;
 else% type = 'prior'
-    NumberOfDraws = options_.prior_draws;
+    B = options_.prior_draws;
-end
+    options_.B = B;
-if ~strcmpi(type,'gsa')
-    B = min([round(.5*NumberOfDraws),500]); options_.B = B;
 end
 try 
     delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])

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');

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;

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
 

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 <http://www.gnu.org/licenses/>.
 
 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!')

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);
+    B = options_.sub_draws;
-    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
 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 = size(x,1);
-    B=NumberOfDraws; 
 elseif strcmpi(type,'prior')
     DirectoryName = CheckPath('prior',M_.dname);
     B = options_.prior_draws;

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 <string_val> 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 ')';

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
 <DYNARE_STATEMENT>dsge_var 	{return token::DSGE_VAR;}
 <DYNARE_STATEMENT>dsge_varlag 	{return token::DSGE_VARLAG;}
 <DYNARE_STATEMENT>moments_varendo {return token::MOMENTS_VARENDO;}
+<DYNARE_STATEMENT>posterior_max_subsample_draws	{return token::POSTERIOR_MAX_SUBSAMPLE_DRAWS;}
 <DYNARE_STATEMENT>filtered_vars	{return token::FILTERED_VARS;}
 <DYNARE_STATEMENT>filter_step_ahead	{return token::FILTER_STEP_AHEAD;}
 <DYNARE_STATEMENT>relative_irf 	{return token::RELATIVE_IRF;}