diff --git a/matlab/check.m b/matlab/check.m
index 6268f23dba12134b082cd39b159ee7a95bf07c1c..950652e2563824b897cbeaa5dba04c9cd754190a 100644
--- a/matlab/check.m
+++ b/matlab/check.m
@@ -39,6 +39,10 @@ end
 
 options_.order = 1;
 
+if isempty(options_.qz_criterium)
+    options_.qz_criterium = 1+1e-6;
+end
+
 [dr, info] = resol(oo_.steady_state,1);
 
 oo_.dr = dr;
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index e0ec06c86944cc078baf26546cac435eb29dacf2..8b0d238abc9be308d797a66389e4503ef16b2bb2 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -54,6 +54,22 @@ if (options_.diffuse_filter==1) && (options_.lik_init==1)
     options_.lik_init = 3;
 end
 
+%% If options_.lik_init == 1
+%%  set by default options_.qz_criterium to 1-1e-6 
+%%  and check options_.qz_criterium < 1-eps if options_.lik_init == 1
+%% Else set by default options_.qz_criterium to 1+1e-6
+if options_.lik_init == 1
+    if isempty(options_.qz_criterium)
+        options_.qz_criterium = 1-1e-6;
+    elseif options_.qz_criterium > 1-eps
+        error(['estimation: option qz_criterium is too large for estimating ' ...
+               'a stationary model. If your model contains unit roots, use ' ...
+               'option diffuse_filter'])
+    end
+elseif isempty(options_.qz_criterium)
+    options_.qz_criterium = 1+1e-6;
+end
+
 %% If the data are prefiltered then there must not be constants in the
 %% measurement equation of the DSGE model or in the DSGE-VAR model.
 if options_.prefilter == 1
@@ -587,7 +603,7 @@ if options_.mode_compute > 0 && ~options_.mh_posterior_mode_estimation
         neps=10;
         %  Set input parameters. 
         maxy=0;
-        eps=1.0e-9;
+        epsilon=1.0e-9;
         rt_=.10;
         t=15.0;
         ns=10;
@@ -602,7 +618,7 @@ if options_.mode_compute > 0 && ~options_.mh_posterior_mode_estimation
         
         vm=1*ones(npar,1);
         disp(['number of parameters= ' num2str(npar) 'max= '  num2str(maxy) 't=  ' num2str(t)]);
-        disp(['rt_=  '  num2str(rt_) 'eps=  '  num2str(eps) 'ns=  '  num2str(ns)]);
+        disp(['rt_=  '  num2str(rt_) 'eps=  '  num2str(epsilon) 'ns=  '  num2str(ns)]);
         disp(['nt=  '  num2str(nt) 'neps= '   num2str(neps) 'maxevl=  '  num2str(maxevl)]);
         %      disp(['iprint=   '   num2str(iprint) 'seed=   '   num2str(seed)]);
         disp '  ';
@@ -615,10 +631,10 @@ if options_.mode_compute > 0 && ~options_.mh_posterior_mode_estimation
         
         %  keyboard 
         if ~options_.dsge_var
-            [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(fh,xparam1,maxy,rt_,eps,ns,nt ...
+            [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(fh,xparam1,maxy,rt_,epsilon,ns,nt ...
                                                               ,neps,maxevl,LB,UB,c,idisp ,t,vm,gend,data,data_index,number_of_observations,no_more_missing_observations);
         else
-            [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(fh,xparam1,maxy,rt_,eps,ns,nt ...
+            [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(fh,xparam1,maxy,rt_,epsilon,ns,nt ...
                                                               ,neps,maxevl,LB,UB,c,idisp ,t,vm,gend);
         end
       otherwise
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 68b0291d2f9dae68ed838fbedd7b8268a745b2ad..5d5799f1c24df619afbbb3b36620f9c4b6f00660 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -46,7 +46,7 @@ options_.initval_file = 0;
 options_.Schur_vec_tol = 1e-11; % used to find nonstationary variables
                                 % in Schur decomposition of the
                                 % transition matrix
-options_.qz_criterium = 1.000001;
+options_.qz_criterium = [];
 options_.lyapunov_complex_threshold = 1e-15;
 options_.solve_tolf = eps^(1/3);
 options_.solve_tolx = eps^(2/3);
diff --git a/matlab/osr.m b/matlab/osr.m
index 671ceebc8dd7b75d8d194ec69a925f4e121bc755..256a0c50f33629c706ca800be5425b349ac2a9c2 100644
--- a/matlab/osr.m
+++ b/matlab/osr.m
@@ -32,6 +32,10 @@ options_ = set_default_option(options_,'hp_ngrid',512);
 options_ = set_default_option(options_,'simul',0);
 options_ = set_default_option(options_,'periods',1);
 
+if isempty(options_.qz_criterium)
+    options_.qz_criterium = 1+1e-6;
+end
+
 make_ex_;
 
 np = size(params,1);
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 615b7cd256274670887bae3ca0bcf2f4d4fe849a..8091f31ddb3a4a711639357f4865b9accca4a196 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -31,6 +31,10 @@ elseif options_.order == 3
     options_.k_order_solver = 1;
 end
 
+if isempty(options_.qz_criterium)
+    options_.qz_criterium = 1+1e-6;
+end
+
 if options_.partial_information == 1 || options_.ACES_solver == 1
     PI_PCL_solver = 1;
     if options_.order ~= 1