From dcc46959b4805c093693b9351e187c0208c0364d Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@ens.fr>
Date: Thu, 23 Dec 2010 08:34:13 +0100
Subject: [PATCH] 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                 |  4 ++++
 matlab/dynare_estimation_1.m   | 24 ++++++++++++++++++++----
 matlab/global_initialization.m |  2 +-
 matlab/osr.m                   |  4 ++++
 matlab/stoch_simul.m           |  4 ++++
 5 files changed, 33 insertions(+), 5 deletions(-)

diff --git a/matlab/check.m b/matlab/check.m
index 6268f23db..950652e25 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 e0ec06c86..8b0d238ab 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 68b0291d2..5d5799f1c 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 671ceebc8..256a0c50f 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 615b7cd25..8091f31dd 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