Merge pull request #1217 from JohannesPfeifer/qz_estimation

Make sure diffuse_filter triggers qz_criterium>1

doc/dynare.texi

@@ -7015,6 +7015,12 @@ more details.
 @item filter_decomposition
 @xref{filter_decomposition}.
 
+@item diffuse_filter = @var{INTEGER}
+@xref{diffuse_filter}.
+
+@item diffuse_kalman_tol = @var{DOUBLE}
+@xref{diffuse_kalman_tol}.
+
 @end table
 
 @end deffn

matlab/dynare_estimation_1.m

@@ -31,6 +31,9 @@ function dynare_estimation_1(var_list_,dname)
 
 global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
 
+%store qz_criterium
+qz_criterium_old=options_.qz_criterium;
+
 % Set particle filter flag.
 if options_.order > 1
     if options_.particle.status && options_.order==2  
@@ -163,6 +166,8 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m
         [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value);
         [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
     end
+    %reset qz_criterium
+    options_.qz_criterium=qz_criterium_old;
     return
 end
 
@@ -421,6 +426,8 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
     %% Here I discard first mh_drop percent of the draws:
     CutSample(M_, options_, estim_params_);
     if options_.mh_posterior_mode_estimation
+        %reset qz_criterium
+        options_.qz_criterium=qz_criterium_old;
         return
     else
         if ~options_.nodiagnostic && options_.mh_replic>0
@@ -467,6 +474,8 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
 end
 
 if options_.particle.status
+    %reset qz_criterium
+    options_.qz_criterium=qz_criterium_old;
     return
 end
 
@@ -713,3 +722,5 @@ if np > 0
     save([M_.fname '_pindx.mat'] ,'pindx');
 end
 
+%reset qz_criterium
+options_.qz_criterium=qz_criterium_old;

matlab/dynare_estimation_init.m

@@ -123,22 +123,7 @@ if isequal(options_.diffuse_filter,1) || (options_.kalman_algo>2)
     end
 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 isequal(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
+options_=select_qz_criterium_value(options_);
 
 % Set options related to filtered variables.
 if ~isequal(options_.filtered_vars,0) && isempty(options_.filter_step_ahead)

matlab/evaluate_smoother.m

@@ -51,6 +51,9 @@ global options_ M_ bayestopt_ oo_ estim_params_   % estim_params_ may be emty
 
 persistent dataset_ dataset_info
 
+%store qz_criterium
+qz_criterium_old=options_.qz_criterium;
+
 if ischar(parameters) && strcmp(parameters,'calibration')
     options_.smoother=1;
 end
@@ -100,4 +103,7 @@ end
 
 if nargout==2
    Smoothed_variables_declaration_order_deviation_form=atT(oo_.dr.inv_order_var(bayestopt_.smoother_var_list),:);
-end
+end
+
+%reset qz_criterium
+options_.qz_criterium=qz_criterium_old;

matlab/imcforecast.m

@@ -119,6 +119,10 @@ if estimated_model
     data_index = dataset_info.missing.aindex;
     gend = dataset_.nobs;
     missing_value = dataset_info.missing.state;
+    
+    %store qz_criterium
+    qz_criterium_old=options_.qz_criterium;
+    options_=select_qz_criterium_value(options_);
     [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam,gend,data,data_index,missing_value);
     %get constant part
     if options_.noconstant
@@ -146,16 +150,16 @@ if estimated_model
     trend = constant(oo_.dr.order_var,:);
     InitState(:,1) = atT(:,end);
 else
+    qz_criterium_old=options_.qz_criterium;
+    if isempty(options_.qz_criterium)
+        options_.qz_criterium = 1+1e-6;
+    end
     graph_title='Calibration';
     if ~isfield(oo_.dr,'kstate')
         error('You need to call stoch_simul before conditional_forecast')
     end
 end
 
-if isempty(options_.qz_criterium)
-    options_.qz_criterium = 1+1e-6;
-end
-
 [T,R,ys,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
 
 if ~isdiagonal(M_.Sigma_e)
@@ -273,4 +277,7 @@ end
 forecasts.graph.title=graph_title;
 forecasts.graph.fname=M_.fname;
 
+%reset qz_criterium
+options_.qz_criterium=qz_criterium_old;
+
 save('conditional_forecasts.mat','forecasts');

matlab/select_qz_criterium_value.m

@@ -0,0 +1,56 @@
+function options_=select_qz_criterium_value(options_)
+% function options_=select_qz_criterium_value(options_)
+% set the value of options_.qz_criterium depending on the Kalman filter used
+%
+% INPUTS
+%   options_:   Dynare options structure
+%
+% OUTPUTS
+%   options_:   Dynare options structure
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2016 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+
+% 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
+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'])
+    end
+else
+    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

preprocessor/DynareBison.yy

@@ -2579,6 +2579,8 @@ calib_smoother_option : o_filtered_vars
                       | o_loglinear
                       | o_first_obs
                       | o_filter_decomposition
+                      | o_diffuse_kalman_tol
+                      | o_diffuse_filter
                       ;
 
 extended_path : EXTENDED_PATH ';'

tests/fs2000/fs2000_calib.mod

@@ -70,3 +70,6 @@ varobs gp_obs gy_obs;
 
 calib_smoother(datafile=fsdat_simul, filtered_vars, filter_step_ahead = [3:4]) m P c e W R k d n l y dA;
 
+%run calibrated diffuse filter
+bayestopt_=[];
+calib_smoother(diffuse_filter,datafile=fsdat_simul, filtered_vars, filter_step_ahead = [3:4]) m P c e W R k d n l y dA;