diff --git a/doc/dynare.texi b/doc/dynare.texi
index 3449c209c37def776cccf40342c3b15cb66d7c85..9e980d681e594bfe660aada6366ecfd7c401a3ef 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3480,6 +3480,7 @@ end;
 
 @deffn Command check ;
 @deffnx Command check (solve_algo = @var{INTEGER}) ;
+@anchor{check}
 
 @descriptionhead
 
@@ -6345,6 +6346,22 @@ Sets the method for approximating the particle distribution. Possible values for
 @item cpf_weights = @var{OPTION}
 @anchor{cpf_weights} Controls the method used to update the weights in conditional particle filter, possible values are @code{amisanotristani} (@cite{Amisano et al (2010)}) or @code{murrayjonesparslow} (@cite{Murray et al. (2013)}). Default value is @code{amisanotristani}.
 
+@item nonlinear_filter_initialization = @var{INTEGER}
+@anchor{nonlinear_filter_initialization} Sets the initial condition of the
+nonlinear filters. By default the nonlinear filters are initialized with the
+unconditional covariance matrix of the state variables, computed with the
+reduced form solution of the first order approximation of the model. If
+@code{nonlinear_filter_initialization=2}, the nonlinear filter is instead
+initialized with a covariance matrix estimated with a stochastic simulation of
+the reduced form solution of the second order approximation of the model. Both
+these initializations assume that the model is stationary, and cannot be used
+if the model has unit roots (which can be seen with the @ref{check} command
+prior to estimation). If the model has stochastic trends, user must use
+@code{nonlinear_filter_initialization=3}, the filters are then initialized with
+an identity matrix for the covariance matrix of the state variables. Default
+value is @code{nonlinear_filter_initialization=1} (initialization based on the
+first order approximation of the model).
+
 @end table
 
 
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index 3debf8b7fccc4b7670ab635d1fbdf966e68613b3..57e02a0b7a0a6c0968ebcfbe592c95911fdfd869 100644
--- a/matlab/non_linear_dsge_likelihood.m
+++ b/matlab/non_linear_dsge_likelihood.m
@@ -316,7 +316,8 @@ switch DynareOptions.particle.initialization
     StateVectorVariance = cov(y_');
     DynareOptions.periods = old_DynareOptionsperiods;
     clear('old_DynareOptionsperiods','y_');
-  case 3% Initial state vector covariance is a diagonal matrix.
+  case 3% Initial state vector covariance is a diagonal matrix (to be used
+        % if model has stochastic trends).
     StateVectorMean = ReducedForm.constant(mf0);
     StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables);
   otherwise
diff --git a/matlab/select_qz_criterium_value.m b/matlab/select_qz_criterium_value.m
index 09282747ea32c93cb7d3777d3ea7b7c5b8eb2ac8..b2e4225ea1aff9c9798a15b5991ad4d666301daf 100644
--- a/matlab/select_qz_criterium_value.m
+++ b/matlab/select_qz_criterium_value.m
@@ -36,21 +36,36 @@ function options_=select_qz_criterium_value(options_)
 %     set by default options_.qz_criterium to 1+1e-6
 stack = dbstack;
 
-if isequal(options_.lik_init,1)
-    if isempty(options_.qz_criterium)
-        options_.qz_criterium = 1-1e-6;
-    elseif options_.qz_criterium > 1-eps
-        error([stack(2).file ': option qz_criterium is too large for estimating/smoothing ' ...
-               'a stationary model. If your model contains unit roots, use ' ...
-               'option diffuse_filter'])
+if options_.particle.status
+    % Non linear filter
+    if isequal(options_.particle.initialization, 3)
+        if isempty(options_.qz_criterium)
+            options_.qz_criterium = 1+1e-6;
+        else
+            if options_.qz_criterium <= 1
+                fprintf('\n%s:: You set nonlinear_filter_initialization equal to 3, it is assumed that you try to estimate a non stationary model. Resetting it to 1+1e-6.\n', stack(2).file)
+                options_.qz_criterium = 1+1e-6;
+            end
+        end
     end
 else
-    if isempty(options_.qz_criterium)
-        options_.qz_criterium = 1+1e-6;
+    % Linear filter
+    if isequal(options_.lik_init,1)
+        if isempty(options_.qz_criterium)
+            options_.qz_criterium = 1-1e-6;
+        elseif options_.qz_criterium > 1-eps
+            error([stack(2).file ': option qz_criterium is too large for estimating/smoothing ' ...
+                   'a stationary model. If your model contains unit roots, use ' ...
+                   'option diffuse_filter'])
+        end
     else
-        if options_.qz_criterium <= 1
-            fprintf('\n%s:: diffuse filter is incompatible with a qz_criterium<=1. Resetting it to 1+1e-6.\n',stack(2).file)
+        if isempty(options_.qz_criterium)
             options_.qz_criterium = 1+1e-6;
+        else
+            if options_.qz_criterium <= 1
+                fprintf('\n%s:: diffuse filter is incompatible with a qz_criterium<=1. Resetting it to 1+1e-6.\n',stack(2).file)
+                options_.qz_criterium = 1+1e-6;
+            end
         end
     end
-end
+end
\ No newline at end of file
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index e5f71d6633d99aa324bd866eb7ad39c3bd691978..c495801a03ecbbfb772cddd5ac706ea2636f9993 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -111,7 +111,7 @@ class ParsingDriver;
 %token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO CONTEMPORANEOUS_CORRELATION DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL RAFTERY_LEWIS_QRS RAFTERY_LEWIS_DIAGNOSTICS MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
 %token NUMBER_OF_PARTICLES RESAMPLING SYSTEMATIC GENERIC RESAMPLING_THRESHOLD RESAMPLING_METHOD KITAGAWA STRATIFIED SMOOTH
 %token CPF_WEIGHTS AMISANOTRISTANI MURRAYJONESPARSLOW
-%token FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
+%token NONLINEAR_FILTER_INITIALIZATION FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
 %token <string_val> NAME
 %token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN INIT_STATE
 %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
@@ -1818,6 +1818,7 @@ estimation_options : o_datafile
 		   | o_resampling_threshold
 		   | o_resampling_method
 		   | o_filter_algorithm
+                   | o_nonlinear_filter_initialization
                    | o_cpf_weights
 		   | o_proposal_approximation
 		   | o_distribution_approximation
@@ -3085,6 +3086,7 @@ o_resampling_method : RESAMPLING_METHOD EQUAL KITAGAWA {driver.option_num("parti
 o_cpf_weights : CPF_WEIGHTS EQUAL AMISANOTRISTANI {driver.option_num("particle.cpf_weights_method.amisanotristani", "1"); driver.option_num("particle.cpf_weights_method.murrayjonesparslow", "0"); }
               | CPF_WEIGHTS EQUAL MURRAYJONESPARSLOW {driver.option_num("particle.cpf_weights_method.amisanotristani", "0"); driver.option_num("particle.cpf_weights_method.murrayjonesparslow", "1"); };
 o_filter_algorithm : FILTER_ALGORITHM EQUAL symbol { driver.option_str("particle.filter_algorithm", $3); };
+o_nonlinear_filter_initialization : NONLINEAR_FILTER_INITIALIZATION EQUAL INT_NUMBER { driver.option_num("particle.initialization", $3); };
 o_proposal_approximation : PROPOSAL_APPROXIMATION EQUAL CUBATURE {driver.option_num("particle.proposal_approximation.cubature", "1"); driver.option_num("particle.proposal_approximation.unscented", "0"); driver.option_num("particle.proposal_approximation.montecarlo", "0");}
 		| PROPOSAL_APPROXIMATION EQUAL UNSCENTED {driver.option_num("particle.proposal_approximation.cubature", "0"); driver.option_num("particle.proposal_approximation.unscented", "1"); driver.option_num("particle.proposal_approximation.montecarlo", "0");}
 		| PROPOSAL_APPROXIMATION EQUAL MONTECARLO {driver.option_num("particle.proposal_approximation.cubature", "0"); driver.option_num("particle.proposal_approximation.unscented", "0"); driver.option_num("particle.proposal_approximation.montecarlo", "1");} ;
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 299d31eddaccf8af091a53cde112e86973cd8aef..0a1588c5ae79d33eeb578b6501dcdb9f51dd0c81 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -398,6 +398,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>amisanotristani {return token::AMISANOTRISTANI;}
 <DYNARE_STATEMENT>murrayjonesparslow {return token::MURRAYJONESPARSLOW;}
 <DYNARE_STATEMENT>filter_algorithm {return token::FILTER_ALGORITHM;}
+<DYNARE_STATEMENT>nonlinear_filter_initialization {return token::NONLINEAR_FILTER_INITIALIZATION;}
 <DYNARE_STATEMENT>proposal_approximation {return token::PROPOSAL_APPROXIMATION;}
 <DYNARE_STATEMENT>cubature {return token::CUBATURE;}
 <DYNARE_STATEMENT>unscented {return token::UNSCENTED;}