diff --git a/doc/dynare.texi b/doc/dynare.texi
index bafb4438a0ac617da8448c4e335c9899b01f3396..9b85534ae783bb4d55074d23facd0ce34869c733 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3090,6 +3090,13 @@ problems. Default: @code{1.000001} (except when estimating with
 Number of simulated series used to compute the IRFs. Default: @code{1}
 if @code{order}=@code{1}, and @code{50} otherwise.
 
+@item simul_replic = @var{INTEGER}
+Number of series to simulate when empirical moments are requested
+(@i{i.e.} @code{periods} > 0). Note that if this option is greater
+than @code{1}, the additional series will not be used for computing
+the empirical moments but will simply be saved in binary form to the
+file @file{@var{FILENAME}_simul}. Default: @code{1}.
+
 @item solve_algo = @var{INTEGER}
 @xref{solve_algo}, for the possible values and their meaning.
 
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 0fa8a012ab182a32890371652f41e07e4b8e6f8e..1711db5b35a9ea5fec19617f1422d4ec97f74abb 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -238,7 +238,7 @@ options_.pruning = 0;
 options_.solve_algo = 2;
 options_.linear = 0;
 options_.replic = 50;
-options_.simul_replic = 0;
+options_.simul_replic = 1;
 options_.drop = 100;
 % if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path.
 % There exists now qz/mjdgges.m that contains the calls to the old Sims code
diff --git a/matlab/simult.m b/matlab/simult.m
index 827c3b7d4ad8915a900b346b837af3fe32a310f8..12196d4324961c96487f4afc6bddcc45421c302a 100644
--- a/matlab/simult.m
+++ b/matlab/simult.m
@@ -35,10 +35,6 @@ global M_ options_ oo_
 order = options_.order;
 replic = options_.simul_replic;
 
-if replic == 0
-    replic = 1;
-end
-
 if replic > 1
     fname = [M_.fname,'_simul'];
     fh = fopen(fname,'w+');
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index e1d2a6041279f75aaaafd131436b861f972aa28e..3c298b87c08f8420b201acc8a058157d32eb25cd 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -122,7 +122,7 @@ class ParsingDriver;
 %token PRINT PRIOR_MC PRIOR_TRUNC PRIOR_MODE PRIOR_MEAN POSTERIOR_MODE POSTERIOR_MEAN POSTERIOR_MEDIAN PRUNING
 %token <string_val> QUOTED_STRING
 %token QZ_CRITERIUM FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE
-%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY
+%token RELATIVE_IRF REPLIC SIMUL_REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY
 %token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO SOLVER_PERIODS
 %token STDERR STEADY STOCH_SIMUL SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME
 %token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY DISCRETIONARY_TOL
@@ -926,6 +926,7 @@ stoch_simul_options : o_dr_algo
                     | o_periods
                     | o_simul
                     | o_simul_seed
+                    | o_simul_replic
                     | o_qz_criterium
                     | o_print
                     | o_noprint
@@ -2191,6 +2192,7 @@ o_markowitz : MARKOWITZ EQUAL non_negative_number { driver.option_num("markowitz
 o_minimal_solving_periods : MINIMAL_SOLVING_PERIODS EQUAL non_negative_number { driver.option_num("minimal_solving_periods", $3); };
 o_mfs : MFS EQUAL INT_NUMBER { driver.mfs($3); };
 o_simul : SIMUL; // Do nothing, only here for backward compatibility
+o_simul_replic : SIMUL_REPLIC EQUAL INT_NUMBER { driver.option_num("replic_simul", $3); };
 o_simul_seed : SIMUL_SEED EQUAL INT_NUMBER { driver.error("'simul_seed' option is no longer supported; use 'set_dynare_seed' command instead"); } ;
 o_qz_criterium : QZ_CRITERIUM EQUAL non_negative_number { driver.option_num("qz_criterium", $3); };
 o_file : FILE EQUAL filename { driver.option_str("file", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 43d8d35da17b42ec37449b977a610c4b7bc2ffbb..4e965e71aa3ab02ad53307747d56a4bf312a0d31 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -517,6 +517,7 @@ string eofbuff;
 <DYNARE_STATEMENT>simul_seed {return token::SIMUL_SEED;}
 <DYNARE_STATEMENT>qz_criterium {return token::QZ_CRITERIUM;}
 <DYNARE_STATEMENT>simul {return token::SIMUL;}
+<DYNARE_STATEMENT>simul_replic {return token::SIMUL_REPLIC;}
 <DYNARE_STATEMENT>xls_sheet {return token::XLS_SHEET;}
 <DYNARE_STATEMENT>xls_range {return token::XLS_RANGE;}
 <DYNARE_STATEMENT>mh_recover {return token::MH_RECOVER;}