diff --git a/doc/dynare.texi b/doc/dynare.texi
index 6a814468930106ae7fe116fbccb9b082b269081b..82b0dfa3bf4bcb14760cffcd9123da712df178ac 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5301,10 +5301,30 @@ variables. This is equivalent to manually listing all the observed
 variables after the @code{estimation} command.
 
 @item number_of_particles = @var{INTEGER}
+@anchor{number_of_particles}
 Number of particles used when evaluating the likelihood of a non linear state space model. Default: @code{1000}.
 
+@item resampling = @var{OPTION}
+@anchor{resampling}
+Determines if resampling of the particles is done. Possible values for @var{OPTION} are:
+
+@table @code
+
+@item none
+No resampling.
+
+@item systematic
+Resampling at each iteration, this is the default value.
+
+@item generic
+Resampling if and only if the effective sample size is below a certain level defined by @code{resampling_threshold*number_of_particles}.
+
 @end table
 
+
+@end table
+
+
 @customhead{Note}
 
 If no @code{mh_jscale} parameter is used in estimated_params, the
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 1a1840c915ef7c0730b6c0cd2539e6341ad404b4..f7afd262a4e98e359ec90607aa80644dffcbf705 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -236,7 +236,10 @@ particle.unscented.alpha = .01;
 particle.unscented.beta = 2;
 particle.unscented.kappa = 0;
 % Configuration of resampling in case of particles
-particle.resampling.status = 'systematic'; % 'none', 'generic', 'smoothed'
+particle.resampling.status.systematic = 1;
+particle.resampling.status.none = 0;
+particle.resampling.status.generic = 0;
+
 particle.resampling.neff_threshold = .5;
 % Choice of the resampling method
 particle.resampling.method1 = 'traditional' ;
diff --git a/matlab/particle/conditional_particle_filter.m b/matlab/particle/conditional_particle_filter.m
index 941e01331ecd2c3f625607a509a8b2ed6d1af9b6..36244809b0d3765fb70d236a4ff047dd237e20d1 100644
--- a/matlab/particle/conditional_particle_filter.m
+++ b/matlab/particle/conditional_particle_filter.m
@@ -113,8 +113,7 @@ for t=1:sample_size
     SumSampleWeights = sum(SampleWeights) ;
     lik(t) = log(SumSampleWeights) ; 
     SampleWeights = SampleWeights./SumSampleWeights ;		
-    if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
-        strcmp(DynareOptions.particle.resampling.status,'systematic')
+    if (DynareOptions.particle.resampling.status.generic && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size) || DynareOptions.particle.resampling.status.systematic
         ks = ks + 1 ;
         StateParticles = resample(StateParticles',SampleWeights',DynareOptions)';
         SampleWeights = ones(1,number_of_particles)/number_of_particles ;
diff --git a/matlab/particle/gaussian_filter.m b/matlab/particle/gaussian_filter.m
index ebe2b1d31aa6bb333dae703490c43fa1324a7078..04edbcc5d018163f7755d6a3f76d0dc5e87c1b9d 100644
--- a/matlab/particle/gaussian_filter.m
+++ b/matlab/particle/gaussian_filter.m
@@ -146,14 +146,13 @@ for t=1:sample_size
         lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ; 
         SampleWeights = SampleWeights./SumSampleWeights ;		
         Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
-        if (Neff<.5*sample_size && strcmpi(DynareOptions.particle.resampling.status,'generic')) || ... 
-                                   strcmpi(DynareOptions.particle.resampling.status,'systematic')
+        if (Neff<.5*sample_size && DynareOptions.particle.resampling.status.generic) || DynareOptions.particle.resampling.status.systematic
             ks = ks + 1 ;
             StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ; 
             StateVectorMean = mean(StateParticles,2) ;
             StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
             SampleWeights = 1/number_of_particles ;
-        elseif strcmpi(DynareOptions.particle.resampling.status,'none')
+        elseif DynareOptions.particle.resampling.status.none
             StateVectorMean = (sampleWeights*StateParticles)' ;
             temp = sqrt(SampleWeights').*StateParticles ;
             StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )';
diff --git a/matlab/particle/gaussian_mixture_filter.m b/matlab/particle/gaussian_mixture_filter.m
index 8f05ea390c6f1d1c32a427b46cf1ff643172e5fe..0506b3c507ad5f8ac58dc55bcc10d9e70fca0940 100644
--- a/matlab/particle/gaussian_mixture_filter.m
+++ b/matlab/particle/gaussian_mixture_filter.m
@@ -202,13 +202,13 @@ for t=1:sample_size
         %estimate(t,:,1) = SampleWeights*StateParticles' ;
         % Resampling if needed of required 
         Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
-        if (Neff<.5*sample_size && strcmpi(DynareOptions.particle.resampling.status,'generic')) || strcmpi(DynareOptions.particle.resampling.status,'systematic')
+        if (DynareOptions.particle.resampling.status.generic && Neff<.5*sample_size) || DynareOptions.particle.resampling.status.systematic
             ks = ks + 1 ;
             StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ; 
             StateVectorMean = mean(StateParticles,2) ;
             StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )';
             SampleWeights = 1/number_of_particles ;
-        elseif strcmpi(DynareOptions.particle.resampling.status,'none')
+        elseif DynareOptions.particle.resampling.status.none
             StateVectorMean = StateParticles*sampleWeights ;
             temp = sqrt(SampleWeights').*StateParticles ;
             StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';
diff --git a/matlab/particle/sequential_importance_particle_filter.m b/matlab/particle/sequential_importance_particle_filter.m
index 6474f6cdab7ce2d8215c81eb97e7305db8f9a91d..3c86f174625bdc9a76e666c5af58204e149eda9d 100644
--- a/matlab/particle/sequential_importance_particle_filter.m
+++ b/matlab/particle/sequential_importance_particle_filter.m
@@ -158,8 +158,7 @@ for t=1:sample_size
     wtilde = weights.*exp(lnw-dfac);
     lik(t) = log(sum(wtilde))+dfac;
     weights = wtilde/sum(wtilde);
-    if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
-        strcmp(DynareOptions.particle.resampling.status,'systematic')
+    if (DynareOptions.particle.resampling.status.generic && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size) || DynareOptions.particle.resampling.status.systematic
         if pruning
             temp = resample([tmp(mf0,:)' tmp_(mf0,:)'],weights',DynareOptions);
             StateVectors = temp(:,1:number_of_state_variables)';
@@ -168,7 +167,7 @@ for t=1:sample_size
             StateVectors = resample(tmp(mf0,:)',weights',DynareOptions)';
         end
         weights = ones(1,number_of_particles)/number_of_particles;
-    elseif strcmp(DynareOptions.particle.resampling.status,'none')
+    elseif DynareOptions.particle.resampling.status.none
         StateVectors = tmp(mf0,:);
         if pruning
             StateVectors_ = tmp_(mf0,:);
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 44b42e37aeb94c3253126256a9fb49ca3e3f1eee..2362e0f1ec877d7bbf0b691e3d3248022a07c174 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -107,7 +107,7 @@ class ParsingDriver;
 %token MFS MH_CONF_SIG 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 MOMENT_CALIBRATION
-%token NUMBER_OF_PARTICLES
+%token NUMBER_OF_PARTICLES RESAMPLING SYSTEMATIC GENERIC
 %token <string_val> NAME
 %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
 %token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS
@@ -1678,6 +1678,7 @@ estimation_options : o_datafile
                    | o_consider_all_endogenous
                    | o_consider_only_observed
 		   | o_number_of_particles
+		   | o_resampling
                    ;
 
 list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -2641,6 +2642,9 @@ o_bvar_prior_train : BVAR_PRIOR_TRAIN EQUAL INT_NUMBER { driver.option_num("bvar
 o_bvar_replic : BVAR_REPLIC EQUAL INT_NUMBER { driver.option_num("bvar_replic", $3); };
 
 o_number_of_particles : NUMBER_OF_PARTICLES EQUAL INT_NUMBER { driver.option_num("particle.number_of_particles", $3); };
+o_resampling : RESAMPLING EQUAL SYSTEMATIC
+              | RESAMPLING EQUAL NONE {driver.option_num("particle.resampling.status.systematic", "0"); driver.option_num("particle.resampling.status.none", "1"); }
+              | RESAMPLING EQUAL GENERIC {driver.option_num("particle.resampling.status.systematic", "0"); driver.option_num("particle.resampling.status.generic", "1"); };
 
 o_gsa_identification : IDENTIFICATION EQUAL INT_NUMBER { driver.option_num("identification", $3); }; /*not in doc */
 o_gsa_morris : MORRIS EQUAL INT_NUMBER { driver.option_num("morris", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 37386530a0ac57fbbda0b7fb7709a0e59b5f58fe..f1fab51561fc8a4cb39338014c84eacc0a003348 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -352,7 +352,13 @@ 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>use_univariate_filters_if_singularity_is_detected {return token::USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED;}
 <DYNARE_STATEMENT>hybrid {return token::HYBRID;}
 <DYNARE_STATEMENT>default {return token::DEFAULT;}
+
 <DYNARE_STATEMENT>number_of_particles {return token::NUMBER_OF_PARTICLES;}
+<DYNARE_STATEMENT>resampling {return token::RESAMPLING;}
+<DYNARE_STATEMENT>systematic {return token::SYSTEMATIC;}
+<DYNARE_STATEMENT>none {return token::NONE;}
+<DYNARE_STATEMENT>generic {return token::GENERIC;}
+
 <DYNARE_STATEMENT>alpha {
   yylval->string_val = new string(yytext);
   return token::ALPHA;