diff --git a/doc/dynare.texi b/doc/dynare.texi
index ea05ab2c7f763906d2f521d212e0871acc00e0d7..2aa309a6268147acca104a99253a3ccb1587d124 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3801,6 +3801,14 @@ For internal use and testing only.
 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 distirbution 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)}
+
+
 @item mh_nblocks = @var{INTEGER}
 Number of parallel chains for Metropolis-Hastings algorithm. Default:
 @code{2}
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index a9f9f97259209a2acf8d455a233e694b1f140ffc..1df1e8c8045cc8a775985dd0b810fa800c9703a8 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -219,7 +219,7 @@ options_.presample = 0;
 options_.prior_trunc = 1e-10;
 options_.smoother = 0;
 options_.student_degrees_of_freedom = 3;
-options_.subdraws = [];
+options_.sub_draws = [];
 options_.use_mh_covariance_matrix = 0;
 options_.gradient_method = 2;
 options_.gradient_epsilon = 1e-6;
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 87f0986617dd20d7df9596e49f9df9e05a03b27f..eded9eedbf6a2fdc9ffdd2b5578a9349447d476f 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -77,8 +77,8 @@ if strcmpi(type,'posterior')
     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
     clear record;
-    if ~isempty(options_.subdraws)
-        B = options_.subdraws;
+    if ~isempty(options_.sub_draws)
+        B = options_.sub_draws;
         if B > NumberOfDraws
             B = NumberOfDraws;
         end
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 6dc6706cfe4b258437d704726c9e7079eeed1c36..0503ea970ac53ad0b335892384b84f6f2c28471b 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -109,7 +109,7 @@ class ParsingDriver;
 %token 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 MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
-%token MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER
+%token MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS
 %token <string_val> NAME
 %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NOCORR NODIAGNOSTIC NOFUNCTIONS
 %token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF
@@ -1219,6 +1219,7 @@ estimation_options : o_datafile
                    | o_conditional_variance_decomposition
                    | o_cova_compute
                    | o_irf_shocks
+                   | o_sub_draws
                    ;
 
 list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -1922,7 +1923,7 @@ o_diffuse_filter: DIFFUSE_FILTER {driver.option_num("diffuse_filter", "1"); };
 o_plot_priors: PLOT_PRIORS EQUAL INT_NUMBER {driver.option_num("plot_priors", $3); };
 o_aim_solver: AIM_SOLVER {driver.option_num("aim_solver", "1"); };
 o_partial_information : PARTIAL_INFORMATION {driver.option_num("partial_information", "1"); };
-
+o_sub_draws: SUB_DRAWS EQUAL INT_NUMBER {driver.option_num("sub_draws",$3);};
 o_planner_discount : PLANNER_DISCOUNT EQUAL expression { driver.declare_optimal_policy_discount_factor_parameter($3); };
 
 o_bvar_prior_tau : BVAR_PRIOR_TAU EQUAL signed_number { driver.option_num("bvar_prior_tau", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 0ec27c8197839b230b3da24671f2b9bab50c2879..e9230dd4ee0aa7760052cacba479d84add899301 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -230,6 +230,7 @@ string eofbuff;
 <DYNARE_STATEMENT>nocorr	{return token::NOCORR;}
 <DYNARE_STATEMENT>optim		{return token::OPTIM;}
 <DYNARE_STATEMENT>periods	{return token::PERIODS;}
+<DYNARE_STATEMENT>sub_draws	{return token::SUB_DRAWS;}
 <DYNARE_STATEMENT>minimal_solving_periods {return token::MINIMAL_SOLVING_PERIODS;}
 <DYNARE_STATEMENT>markowitz	{return token::MARKOWITZ;}
 <DYNARE_STATEMENT>marginal_density {return token::MARGINAL_DENSITY;}