From 1519e4673a19a16f8b165f3be55c73ed765bfef6 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 29 Oct 2020 17:12:11 +0100 Subject: [PATCH 1/9] trap case when order is too big wrt sample size --- matlab/utilities/dataset/nanautocovariance.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/utilities/dataset/nanautocovariance.m b/matlab/utilities/dataset/nanautocovariance.m index 715ae957f..bc1cf582d 100644 --- a/matlab/utilities/dataset/nanautocovariance.m +++ b/matlab/utilities/dataset/nanautocovariance.m @@ -55,7 +55,7 @@ function autocov = nanautocovariance(data,order) n = size(data,2); missing = isanynan(data); - +order = min(size(data,1)-2,order); autocov = zeros(n, n, order); for lag=1:order From ebec3ede8afaaffe8b90e682c48c3c69e5ad6a56 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 29 Oct 2020 17:13:00 +0100 Subject: [PATCH 2/9] allow negative lag (i.e. a lead) --- matlab/utilities/dataset/lagged.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/matlab/utilities/dataset/lagged.m b/matlab/utilities/dataset/lagged.m index 4f17d7eeb..18e3ca98b 100644 --- a/matlab/utilities/dataset/lagged.m +++ b/matlab/utilities/dataset/lagged.m @@ -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 From b9a6462462060ccc8dacd7c36d8836ac4814c006 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 29 Oct 2020 17:55:34 +0100 Subject: [PATCH 3/9] bug fix in trapping option newratflag --- matlab/dynare_estimation_1.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index cda0ecb02..ef9c2b3f2 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -240,7 +240,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation [~,~,~,~,hh] = feval(objective_function,xparam1, ... 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')) + 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 if options_.hessian.use_penalized_objective penalized_objective_function = str2func('penalty_objective_function'); From 22ef0b1413704cd2dba1f6796db7a055eccfd273 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 29 Oct 2020 17:57:17 +0100 Subject: [PATCH 4/9] get statistics also when load_mh_file --- matlab/GetPosteriorParametersStatistics.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index 6be075ad0..b32d7056e 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 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}; From 15ab85788c1a9e2cd77580dfb9ee11dea97ffdf1 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 29 Oct 2020 19:21:06 +0100 Subject: [PATCH 5/9] adjusted message in case of diffuse filter, where cointegration is most of the time the reason for such a warning --- matlab/DsgeSmoother.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 6e1739698..09842c44b 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') From d56bba724c52c8e58766799bfc8d046e226bb507 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 13 Feb 2021 11:27:18 +0100 Subject: [PATCH 6/9] autocov is initialized NaN with the original order and then only loop up to updted order compatible with sample size. --- matlab/utilities/dataset/nanautocovariance.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/utilities/dataset/nanautocovariance.m b/matlab/utilities/dataset/nanautocovariance.m index bc1cf582d..eb7a81ab2 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 = nan(n, n, order); order = min(size(data,1)-2,order); -autocov = zeros(n, n, order); +autocov(:, :, 1:order)=0; for lag=1:order if missing From 915502ce1309e92221cb6c05b5a481394793cfac Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 13 Feb 2021 11:29:40 +0100 Subject: [PATCH 7/9] copyright year and indentation --- matlab/utilities/dataset/lagged.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/matlab/utilities/dataset/lagged.m b/matlab/utilities/dataset/lagged.m index 18e3ca98b..276f6b785 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. % @@ -32,7 +32,7 @@ end x=x(:); if n>0 -xlag=[NaN(n,1); x(1:end-n)]; + xlag=[NaN(n,1); x(1:end-n)]; else -xlag=[x(abs(n)+1:end); NaN(abs(n),1)]; + xlag=[x(abs(n)+1:end); NaN(abs(n),1)]; end \ No newline at end of file From 5c3b0ebc9525d2847756cd4dd5c4862f1d02be92 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 13 Feb 2021 16:44:50 +0100 Subject: [PATCH 8/9] better handle newrat options, avoiding to treat newratflag in two places --- matlab/dynare_estimation_1.m | 37 +++---------------- .../optimization/dynare_minimize_objective.m | 11 +++++- 2 files changed, 16 insertions(+), 32 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index ef9c2b3f2..3cd3d5864 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; @@ -240,8 +215,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation [~,~,~,~,hh] = feval(objective_function,xparam1, ... 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 + elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood')) + % 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 956072b85..2b9c50d90 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); From faa203e0292ec31e53b5f42d33d51f3e2dfa48f4 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Sat, 13 Feb 2021 18:53:17 +0100 Subject: [PATCH 9/9] condition computation of posterior statistics also to (options_.load_mh_file && ~options_.load_results_after_load_mh) and apply it to all cases --- matlab/GetPosteriorParametersStatistics.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index b32d7056e..d2cf00a69 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 || options_.load_mh_file + 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);