The default for options_.qz_criterium is now different for different contexts:

qz_criterium = 1+1e-6
- stoch_simul
- osr
- check
- dynare_estimation_1 with lik_init =2 or 3
qz_criterium = 1+1e-6 for dynare_estimation_1 with lik_init = 1

This correct a bug in estimation and computation of the posterior distribution of moments of endogenous variables when the priors permit roots too close to 1.

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;

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

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);

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);

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