diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 6e1739698cc9af716094b9d94818fc08b6ac4443..09842c44b54492b115747a1819c756479f22d68f 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -238,7 +238,8 @@ if kalman_algo == 1 || kalman_algo == 3
fprintf('\nDsgeSmoother: Switching to univariate filter. This may be a sign of stochastic singularity.\n')
kalman_algo = 2;
elseif kalman_algo == 3
- fprintf('\nDsgeSmoother: Switching to univariate filter. This may be a sign of stochastic singularity.\n')
+ fprintf('\nDsgeSmoother: Switching to univariate filter. This is usually due to co-integration in diffuse filter,\n')
+ fprintf(' otherwise it may be a sign of stochastic singularity.\n')
kalman_algo = 4;
else
error('This case shouldn''t happen')
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index 6be075ad0eb5ea72492364e514dd448371f35ee0..d2cf00a695a794b38c2600d713d6065a41deab74 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -96,7 +96,7 @@ if np
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
- if options_.mh_replic
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
name = bayestopt_.name{ip};
@@ -139,7 +139,7 @@ if nvx
disp('standard deviation of shocks')
disp(tit2)
for i=1:nvx
- if options_.mh_replic
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
k = estim_params_.var_exo(i,1);
@@ -183,7 +183,7 @@ if nvn
disp(tit2)
ip = nvx+1;
for i=1:nvn
- if options_.mh_replic
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
@@ -222,7 +222,7 @@ if ncx
disp(tit2)
ip = nvx+nvn+1;
for i=1:ncx
- if options_.mh_replic
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig);
k1 = estim_params_.corrx(i,1);
@@ -272,7 +272,7 @@ if ncn
disp(tit2)
ip = nvx+nvn+ncx+1;
for i=1:ncn
- if options_.mh_replic
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
k1 = estim_params_.corrn(i,1);
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index cda0ecb02694b30bb7b6385b62746f61a12d671f..3cd3d5864cd521429454b541b59b758a0bf1a70a 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -195,38 +195,13 @@ end
%% Estimation of the posterior mode or likelihood mode
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
- %prepare settings for newrat
- if options_.mode_compute==5
- %get whether outer product Hessian is requested
- newratflag=[];
- if ~isempty(options_.optim_opt)
- options_list = read_key_value_string(options_.optim_opt);
- for i=1:rows(options_list)
- if strcmp(options_list{i,1},'Hessian')
- newratflag=options_list{i,2};
- end
- end
- end
- if options_.analytic_derivation
- options_analytic_derivation_old = options_.analytic_derivation;
- options_.analytic_derivation = -1;
- if ~isempty(newratflag) && newratflag~=0 %numerical hessian explicitly specified
- error('newrat: analytic_derivation is incompatible with numerical Hessian.')
- else %use default
- newratflag=0; %exclude DYNARE numerical hessian
- end
- elseif ~options_.analytic_derivation
- if isempty(newratflag)
- newratflag=options_.newrat.hess; %use default numerical dynare hessian
- end
- end
- end
[xparam1, fval, exitflag, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,options_.mode_compute,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
- if isnumeric(options_.mode_compute) && options_.mode_compute==5 && options_.analytic_derivation==-1 %reset options changed by newrat
- options_.analytic_derivation = options_analytic_derivation_old; %reset
+ if isnumeric(options_.mode_compute) && options_.mode_compute==5
+ newratflag = new_rat_hess_info.newratflag;
+ new_rat_hess_info = new_rat_hess_info.new_rat_hess_info;
elseif isnumeric(options_.mode_compute) && options_.mode_compute==6 %save scaling factor
save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
options_.mh_jscale = Scale;
@@ -241,7 +216,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
options_.analytic_derivation = ana_deriv_old;
elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood'))
- % with flag==0, we force to use the hessian from outer product gradient of optimizer 5
+ % with flag==0 or 2, we force to use the hessian from outer product gradient of optimizer 5
if options_.hessian.use_penalized_objective
penalized_objective_function = str2func('penalty_objective_function');
hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
@@ -260,7 +235,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
% with diagonal elements computed with numerical second order derivatives
%
% uses univariate filters, so to get max # of available
- % densitities for outer product gradient
+ % densities for outer product gradient
kalman_algo0 = options_.kalman_algo;
compute_hessian = 1;
if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4))
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 956072b8582deabdd3b749edc9b984541ddb4eb8..2b9c50d90bd4d34b2733e45acc1e6b098c01c468 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -274,7 +274,9 @@ switch minimizer_algorithm
if isempty(prior_information) %mr_hessian requires it, but can be NaN
prior_information.p2=NaN(n_params,1);
end
- if options_.analytic_derivation==-1 %set outside as code for use of analytic derivation
+ if options_.analytic_derivation
+ old_analytic_derivation = options_.analytic_derivation;
+ options_.analytic_derivation=-1; %force analytic outer product gradient hessian for each iteration
analytic_grad=1;
crit = options_.newrat.tolerance.f_analytic;
newratflag = 0; %analytical Hessian
@@ -317,7 +319,14 @@ switch minimizer_algorithm
hess_info.gstep=options_.gstep;
hess_info.htol = 1.e-4;
hess_info.h1=options_.gradient_epsilon*ones(n_params,1);
+ % here we force 7th input argument (flagg) to be 0, since outer product
+ % gradient Hessian is handled in dynare_estimation_1
[opt_par_values,hessian_mat,gg,fval,invhess,new_rat_hess_info] = newrat(objective_function,start_par_value,bounds,analytic_grad,crit,nit,0,Verbose,Save_files,hess_info,prior_information.p2,options_.gradient_epsilon,parameter_names,varargin{:}); %hessian_mat is the plain outer product gradient Hessian
+ new_rat_hess_info.new_rat_hess_info = new_rat_hess_info;
+ new_rat_hess_info.newratflag = newratflag;
+ if options_.analytic_derivation
+ options_.analytic_derivation = old_analytic_derivation;
+ end
case 6
if isempty(prior_information) %Inf will be reset
prior_information.p2=Inf(n_params,1);
diff --git a/matlab/utilities/dataset/lagged.m b/matlab/utilities/dataset/lagged.m
index 4f17d7eeb302d84c536050e83bcd7570d913f75e..276f6b785ec49761fc319f04e60be63d0aa0b67f 100644
--- a/matlab/utilities/dataset/lagged.m
+++ b/matlab/utilities/dataset/lagged.m
@@ -9,7 +9,7 @@ function xlag = lagged(x, n)
% OUTPUT
% xlag = backward shifted series
-% Copyright (C) 2017 Dynare Team
+% Copyright (C) 2017-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -31,4 +31,8 @@ if nargin==1
end
x=x(:);
-xlag=[NaN(n,1); x(1:end-n)];
\ No newline at end of file
+if n>0
+ xlag=[NaN(n,1); x(1:end-n)];
+else
+ xlag=[x(abs(n)+1:end); NaN(abs(n),1)];
+end
\ No newline at end of file
diff --git a/matlab/utilities/dataset/nanautocovariance.m b/matlab/utilities/dataset/nanautocovariance.m
index 715ae957f389f66d387053b014a6cd1bda43003b..eb7a81ab2b54e7bf13672f51db0a37b565bb4441 100644
--- a/matlab/utilities/dataset/nanautocovariance.m
+++ b/matlab/utilities/dataset/nanautocovariance.m
@@ -55,8 +55,9 @@ function autocov = nanautocovariance(data,order)
n = size(data,2);
missing = isanynan(data);
-
-autocov = zeros(n, n, order);
+autocov = nan(n, n, order);
+order = min(size(data,1)-2,order);
+autocov(:, :, 1:order)=0;
for lag=1:order
if missing