From c7fc3b3946a4d33799fb0b76b22bc327747c9816 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Argos=29?=
<stepan@adjemian.eu>
Date: Sun, 6 Oct 2024 22:33:11 +0200
Subject: [PATCH] WIP
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
[skip ci]
+ New routine get_user_defined_boundaries()
+ Add @dprior method for plotting prior densities.
- check_mode_file() signature change (do not pass options_ as an
argument, remove bayestopt_), bug fixes (wrong indexing of cell
arrays) and add unit tests.
+ set of private functions under the estimation subfolder (code factorization).
+ initialguess() returns the initial guess for the estimated parameters
+ @dprior/isempty() method, return true iff object is empty
+ derror() Throws an error without the stack trace under
MATLAB. Behaviour unchanged under Octave. Ideally,we should print
the stack trace when in debug mode. Can we test if we are in debug
mode without passing options_?
+ isbayesian() returns tre iff bayesian approach is considered.
- set_prior returns a dprior object.
- set_prior ne devrait en principe pas être appelé quand on ne
considère pas une approche bayésienne... C'est pourtant le cas car
cette fonction retourne aussi la condition initiale et les bornes de
l'optimisateur. Il faut casser cette fonction.
uniform_specification() was throwing an error if the uniform prior is
defined with the expectation and variance.
---
matlab/+identification/analysis.m | 7 +-
matlab/@dprior/bounds.m | 11 +-
matlab/@dprior/data_for_plotting_density.m | 106 ++++
matlab/@dprior/hyperparameters.m | 117 +++++
matlab/@dprior/isempty.m | 96 ++++
matlab/@dprior/subsref.m | 14 +-
matlab/default_option_values.m | 2 +-
matlab/derror.m | 38 ++
matlab/dynare.m | 12 +
.../estimation/PlotPosteriorDistributions.m | 25 +-
...check_bounds_and_definiteness_estimation.m | 27 +-
matlab/estimation/check_prior_bounds.m | 7 +-
matlab/estimation/dsge_likelihood.m | 81 ++-
matlab/estimation/dsge_var_likelihood.m | 72 ++-
matlab/estimation/dynare_estimation_1.m | 18 +-
matlab/estimation/dynare_estimation_init.m | 205 ++++----
.../estimation/get_user_defined_boundaries.m | 77 +++
matlab/estimation/initial_estimation_checks.m | 20 +-
matlab/estimation/initialguess.m | 80 +++
matlab/estimation/isbayesian.m | 150 ++++++
.../estimation/non_linear_dsge_likelihood.m | 8 +-
matlab/estimation/plot_priors.m | 31 +-
matlab/estimation/posterior_sampler.m | 4 +-
.../posterior_sampler_initialization.m | 4 +-
...ize_measurement_errors_covariance_matrix.m | 30 ++
.../set_corrn_observable_correspondence.m | 30 ++
matlab/estimation/private/set_jscale.m | 64 +++
.../set_number_of_estimated_parameters.m | 56 +++
.../set_nvn_observable_correspondence.m | 31 ++
.../private/set_optimization_boundaries.m | 66 +++
.../estimation/private/set_parameter_names.m | 68 +++
matlab/estimation/set_boundaries.m | 68 +++
matlab/estimation/set_prior.m | 419 +++++++--------
matlab/estimation/uniform_specification.m | 45 +-
matlab/isdebug.m | 39 ++
.../select_qz_criterium_value.m | 2 +-
matlab/utilities/estimation/check_mode_file.m | 475 +++++++++++++++---
.../estimation/check_prior_stderr_corr.m | 13 +-
.../check_varobs_are_endo_and_declared_once.m | 18 +-
tests/kalman_filter_fast/swk_0.mod | 2 +-
40 files changed, 2048 insertions(+), 590 deletions(-)
create mode 100644 matlab/@dprior/data_for_plotting_density.m
create mode 100644 matlab/@dprior/hyperparameters.m
create mode 100644 matlab/@dprior/isempty.m
create mode 100644 matlab/derror.m
create mode 100644 matlab/estimation/get_user_defined_boundaries.m
create mode 100644 matlab/estimation/initialguess.m
create mode 100644 matlab/estimation/isbayesian.m
create mode 100644 matlab/estimation/private/initialize_measurement_errors_covariance_matrix.m
create mode 100644 matlab/estimation/private/set_corrn_observable_correspondence.m
create mode 100644 matlab/estimation/private/set_jscale.m
create mode 100644 matlab/estimation/private/set_number_of_estimated_parameters.m
create mode 100644 matlab/estimation/private/set_nvn_observable_correspondence.m
create mode 100644 matlab/estimation/private/set_optimization_boundaries.m
create mode 100644 matlab/estimation/private/set_parameter_names.m
create mode 100644 matlab/estimation/set_boundaries.m
create mode 100644 matlab/isdebug.m
diff --git a/matlab/+identification/analysis.m b/matlab/+identification/analysis.m
index 7e7933b862..be686f4ee6 100644
--- a/matlab/+identification/analysis.m
+++ b/matlab/+identification/analysis.m
@@ -297,11 +297,10 @@ if info(1) == 0 %no errors in solution
options_.analytic_derivation = -2; %this sets asy_Hess=1 in dsge_likelihood.m
[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, options_.varobs);
dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments
+ Prior = dprior(bayestopt_, options_.prior_trunc, false);
derivatives_info.no_DLIK = 1;
- bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
- %note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.
- [~, info, ~, ~, AHess, ~, ~, M_, options_, ~, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state, derivatives_info); %non-used output variables need to be set for octave for some reason
- %note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
+ [~, info, ~, ~, AHess, ~, ~, M_, options_, ~, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, Prior, estim_params_, bayestopt_, Prior.bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state, derivatives_info); %non-used output variables need to be set for octave for some reason
+ %note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
options_.analytic_derivation = analytic_derivation; %reset option
AHess = -AHess; %take negative of hessian
if min(eig(AHess))<-tol_rank
diff --git a/matlab/@dprior/bounds.m b/matlab/@dprior/bounds.m
index d54cc7821d..ab8ca30e3c 100644
--- a/matlab/@dprior/bounds.m
+++ b/matlab/@dprior/bounds.m
@@ -41,8 +41,7 @@ if inplace
end
else
if ~isempty(o.lb) && (nargin<2 || isempty(truncation) || abs(truncation-o.trunc)<eps(1))
- lub = [o.lb, o.ub];
- o = lub;
+ o = struct('lb', o.lb, 'ub', o.ub);
return
end
end
@@ -125,7 +124,7 @@ if inplace
o.lb = lub(:,1);
o.ub = lub(:,2);
else
- o = lub;
+ o = struct('lb', lub(:,1), 'ub', lub(:,2));
end
return % --*-- Unit tests --*--
@@ -209,7 +208,7 @@ try
lb = o.lb;
ub = o.ub;
% Call the tested method
- lub = o.bounds(.1, false);
+ boundaries = o.bounds(.1, false);
t(1) = true;
catch
t(1) = false;
@@ -218,8 +217,8 @@ end
if t(1)
t(2) = all(abs(o.lb-lb)<2*eps(1));
t(3) = all(abs(o.ub-ub)<2*eps(1));
- t(4) = all(lub(:,1)>=lb);
- t(5) = all(lub(:,2)<=ub);
+ t(4) = all(boundaries.lb>=lb);
+ t(5) = all(boundaries.ub<=ub);
end
T = all(t);
diff --git a/matlab/@dprior/data_for_plotting_density.m b/matlab/@dprior/data_for_plotting_density.m
new file mode 100644
index 0000000000..35ddf3dac9
--- /dev/null
+++ b/matlab/@dprior/data_for_plotting_density.m
@@ -0,0 +1,106 @@
+function [x, f, xlb, xub] = data_for_plotting_density(o, idx, trunc, steps)
+
+% Evaluate the logged prior density at specified points.
+%
+% INPUTS
+% - o [dprior] object containing prior distribution information.
+% - idx [integer] scalar, index of the parameter (based on o.name order).
+% - trunc [double] scalar, the probability mass to trim from distribution tails (default: 1e-3).
+% - steps [integer] scalar, the number of evaluation points for the density (default: 200).
+%
+% OUTPUTS
+% - x [double] steps×1 vector representing the grid for density evaluation.
+% - f [double] steps×1 vector, density.
+% - xlb [double] scalar, lower bound of x.
+% - xub [double] scalar, upper bound of x.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ if nargin<4 || isempty(steps)
+ steps = 200;
+ end
+
+ if nargin<3 || isempty(trunc)
+ trunc = 1e-3;
+ end
+
+ if ismember(idx, o.idbeta)
+ density = @(x, a, b, aa, bb) betapdf((x-aa)/(bb-aa), a, b)/(bb-aa);
+ infbound = betainv(trunc, o.p6(idx), o.p7(idx))*(o.p4(idx)-o.p3(idx))+o.p3(idx);
+ supbound = betainv(1-trunc, o.p6(idx), o.p7(idx))*(o.p4(idx)-o.p3(idx))+o.p3(idx);
+ nargs = 4;
+ elseif ismember(idx, o.idgaussian)
+ density = @(x, a, b) normpdf(x, a, b);
+ infbound = max(o.p3(idx), norminv(trunc, o.p6(idx), o.p7(idx)));
+ supbound = min(o.p4(idx), norminv(1-trunc, o.p6(idx), o.p7(idx)));
+ nargs = 2;
+ elseif ismember(idx, o.idgamma)
+ density = @(x,a,b,c) gampdf(x-c, a, b);
+ infbound = gaminv(trunc, o.p6(idx), o.p7(idx))+o.p3(idx);
+ supbound = gaminv(1-trunc, o.p6(idx), o.p7(idx))+o.p3(idx);
+ nargs = 3;
+ elseif ismember(idx, o.idinvgamma1)
+ density = @(x, a, b, c) exp(lpdfig1(x-c, a, b));
+ % TODO: Check out why we increase the value of trunc here...
+ infbound = 1.0/sqrt(gaminv(1.0-10.0*trunc, o.p7(idx)/2.0, 2.0/o.p6(idx)))+o.p3(idx);
+ supbound = 1.0/sqrt(gaminv(10.0*trunc, o.p7(idx)/2.0, 2.0/o.p6(idx)))+o.p3(idx);
+ nargs = 3;
+ elseif ismember(idx, o.idinvgamma2)
+ density = @(x, a, b, c) xp(lpdfig2(x-c, a, b));
+ % TODO: Check out why we increase the value of trunc here...
+ infbound = 1.0/(gaminv(1.0-10.0*trunc, o.p7(idx)/2.0, 2.0/o.p6(idx)))+p3(idx);
+ supbound = 1.0/(gaminv(10.0*trunc, o.p7(idx)/2.0, 2.0/o.p6(idx)))+o.p3(idx);
+ nargs = 3;
+ elseif ismember(idx, o.iduniform)
+ density = @(x, a, b) x/(b-a);
+ infbound = o.p6(idx);
+ supbound = o.p7(idx);
+ nargs = 2;
+ elseif ismember(idx, o.idweibull)
+ density = @(x, a, b, c) exp(lpdfgweibull(x, a, b, c));
+ infbound = o.p3(idx)+wblinv(trunc, o.p6(idx), o.p7(idx));
+ supbound = o.p3(idx)+wblinv(1-trunc, o.p6(idx), o.p7(idx));
+ nargs = 3;
+ else
+ derror('Estimation:Prior:Bug', 'Unknown prior shape.')
+ end
+
+ x = linspace(infbound, supbound, steps);
+
+ xlb = x(1);
+ xub = x(end);
+
+ switch nargs
+ case 2
+ f = density(x, o.p6(idx), o.p7(idx));
+ case 3
+ f = density(x, o.p6(idx), o.p7(idx), o.p3(idx));
+ case 4
+ f = density(x, o.p6(idx), o.p7(idx), o.p3(idx), o.p4(idx));
+ otherwise
+ end
+
+ if ~ismember(idx, o.iduniform)
+ % Cut densities with mode at the lower or upper bound
+ [~, k] = max(f);
+ if k==1 || k==length(x)
+ f(f>10) = NaN;
+ end
+ end
+
+end
diff --git a/matlab/@dprior/hyperparameters.m b/matlab/@dprior/hyperparameters.m
new file mode 100644
index 0000000000..72b80985ea
--- /dev/null
+++ b/matlab/@dprior/hyperparameters.m
@@ -0,0 +1,117 @@
+function o = hyperparameters(o, varargin)
+
+% Compute hyperparameters of the prior distributions
+%
+% INPUTS
+% - o [dprior] Distribution specification for a n×1 vector of independent continuous random variables
+% - truncation [double] scalar, prior mass excluded in left and right parts of the prior distribution. Defaut: 0 if inplace==true.
+% - inplace [logical] scalar, if true modifies o in place.
+%
+% OUTPUTS
+% - lub [double] n×2 matrix, lower and upper bounds of the prior support
+%
+% REMARKS
+% o.p3 and o.p4 hold the lower and upper natural boundaries of the prior supports, o.lb and o.ub hold the lower and upper bounds
+% of the truncated prior distributions (excluding probability mass in the tails).
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ if mod(number,2) ~= 0
+ derror('Prior:hyperparameters', 'Input arguments must come by key/value pairs.')
+ end
+
+ if o.isbeta
+ if truncation<eps(1)
+ lub(o.idbeta,1) = o.p3(o.idbeta);
+ lub(o.idbeta,2) = o.p4(o.idbeta);
+ else
+ lub(o.idbeta,1) = betainv(truncation, o.p6(o.idbeta), o.p7(o.idbeta)).*(o.p4(o.idbeta)-o.p3(o.idbeta))+o.p3(o.idbeta);
+ lub(o.idbeta,2) = betainv(1.0-truncation, o.p6(o.idbeta), o.p7(o.idbeta)).*(o.p4(o.idbeta)-o.p3(o.idbeta))+o.p3(o.idbeta);
+ end
+ end
+
+ if o.isgamma
+ if truncation<eps(1)
+ lub(o.idgamma,1) = o.p3(o.idgamma);
+ lub(o.idgamma,2) = Inf;
+ else
+ lub(o.idgamma,1) = gaminv(truncation, o.p6(o.idgamma), o.p7(o.idgamma))+o.p3(o.idgamma);
+ lub(o.idgamma,2) = gaminv(1.0-truncation, o.p6(o.idgamma), o.p7(o.idgamma))+o.p3(o.idgamma);
+ end
+ end
+
+ if o.isgaussian
+ if truncation<eps(1)
+ lub(o.idgaussian,1) = max(-Inf, o.p3(o.idgaussian));
+ lub(o.idgaussian,2) = min(Inf, o.p4(o.idgaussian));
+ else
+ lub(o.idgaussian,1) = max(norminv(truncation, o.p6(o.idgaussian), o.p7(o.idgaussian)), o.p3(o.idgaussian));
+ lub(o.idgaussian,2) = min(norminv(1-truncation, o.p6(o.idgaussian), o.p7(o.idgaussian)), o.p4(o.idgaussian));
+ end
+ end
+
+ if o.isinvgamma1
+ if truncation<eps(1)
+ lub(o.idinvgamma1,1) = o.p3(o.idinvgamma1);
+ lub(o.idinvgamma1,2) = Inf;
+ else
+ lub(o.idinvgamma1,1) = 1.0./sqrt(gaminv(1.0-truncation, o.p7(o.idinvgamma1)/2.0, 2.0./o.p6(o.idinvgamma1)))+o.p3(o.idinvgamma1);
+ lub(o.idinvgamma1,2) = 1.0./sqrt(gaminv(truncation, o.p7(o.idinvgamma1)/2.0, 2.0./o.p6(o.idinvgamma1)))+o.p3(o.idinvgamma1);
+ end
+ end
+
+ if o.isuniform
+ if truncation<eps(1)
+ lub(o.iduniform,1) = o.p6(o.iduniform);
+ lub(o.iduniform,2) = o.p7(o.iduniform);
+ else
+ lub(o.iduniform,1) = o.p6(o.iduniform)+(o.p7(o.iduniform)-o.p6(o.iduniform))*truncation;
+ lub(o.iduniform,2) = o.p7(o.iduniform)-(o.p7(o.iduniform)-o.p6(o.iduniform))*truncation;
+ end
+ end
+
+ if o.isinvgamma2
+ if truncation<eps(1)
+ lub(o.idinvgamma2,1) = o.p3(o.idinvgamma2);
+ lub(o.idinvgamma2,2) = Inf;
+ else
+ lub(o.idinvgamma2,1) = 1.0./gaminv(1.0-truncation, o.p7(o.idinvgamma2)/2.0, 2.0./o.p6(o.idinvgamma2))+o.p3(o.idinvgamma2);
+ lub(o.idinvgamma2,2) = 1.0./gaminv(truncation, o.p7(o.idinvgamma2)/2.0, 2.0./o.p6(o.idinvgamma2))+o.p3(o.idinvgamma2);
+ end
+ end
+
+ if o.isweibull
+ if truncation<eps(1)
+ lub(o.idweibull,1) = o.p3(o.idweibull);
+ lub(o.idweibull,2) = Inf;
+ else
+ lub(o.idweibull,1) = o.p3(o.idweibull)+wblinv(truncation, o.p6(o.idweibull), o.p7(o.idweibull));
+ lub(o.idweibull,2) = o.p3(o.idweibull)+wblinv(1.0-truncation, o.p6(o.idweibull), o.p7(o.idweibull));
+ end
+ end
+
+ if inplace
+ o.lb = lub(:,1);
+ o.ub = lub(:,2);
+ else
+ o = lub;
+ end
+
+ return % --*-- Unit tests --*--
+
+end
diff --git a/matlab/@dprior/isempty.m b/matlab/@dprior/isempty.m
new file mode 100644
index 0000000000..e2a631db54
--- /dev/null
+++ b/matlab/@dprior/isempty.m
@@ -0,0 +1,96 @@
+function b = isempty(o)
+
+% Return true iff dprior object is empty.
+%
+% INPUTS
+% - o [dprior]
+%
+% OUTPUTS
+% - b [logical] scalar, true if o is empty
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ b = ~o.isbeta && ~o.isgaussian && ~o.isuniform && ~o.isgamma && ~o.isinvgamma1 && ~o.isinvgamma2 && ~o.isweibull;
+
+ return % --*-- Unit tests --*--
+
+ %@test:1
+ o = dprior();
+
+ try
+ b = o.isempty();
+ t(1) = true;
+ catch
+ t(1) = false;
+ end
+
+ if t(1)
+ t(2) = b;
+ end
+
+ T = all(t);
+ %@eof:1
+
+
+ %@test:2
+ prior_trunc = 1e-10;
+ p0 = [1; 2]; % Prior shape
+ p1 = .4*ones(2,1); % Prior mean
+ p2 = .2*ones(2,1); % Prior std.
+ p3 = NaN(2,1);
+ p4 = NaN(2,1);
+ p5 = NaN(2,1);
+ p6 = NaN(2,1);
+ p7 = NaN(2,1);
+
+ p3(1) = 0;
+ p4(1) = 1;
+ [p6(1), p7(1)] = beta_specification(p1(1), p2(1)^2, p3(1), p4(1));
+ p5(1) = compute_prior_mode([p6(1) p7(1)], 1);
+
+ p3(2) = 0;
+ p4(2) = Inf;
+ [p6(2), p7(2)] = gamma_specification(p1(2), p2(2)^2, p3(2), p4(2));
+ p5(2) = compute_prior_mode([p6(2) p7(2)], 2);
+
+ BayesInfo.pshape = p0;
+ BayesInfo.p1 = p1;
+ BayesInfo.p2 = p2;
+ BayesInfo.p3 = p3;
+ BayesInfo.p4 = p4;
+ BayesInfo.p5 = p5;
+ BayesInfo.p6 = p6;
+ BayesInfo.p7 = p7;
+
+ o = dprior(BayesInfo, prior_trunc);
+
+ try
+ b = o.isempty();
+ t(1) = true;
+ catch
+ t(2) = false;
+ end
+
+ if t(1)
+ t(2) = ~b;
+ end
+
+ T = all(t);
+ %@eof:2
+
+end
diff --git a/matlab/@dprior/subsref.m b/matlab/@dprior/subsref.m
index be49924f8f..5bab04aa3c 100644
--- a/matlab/@dprior/subsref.m
+++ b/matlab/@dprior/subsref.m
@@ -1,4 +1,4 @@
-function p = subsref(o, S)
+function varargout = subsref(o, S)
% Overload subsref method.
@@ -56,16 +56,16 @@ switch S(1).type
derror('dprior:wrongInputArguments', 'First input argument to ishape must be the name of a distribution (row char array).')
end
elseif isequal(S(1).subs, 'bounds') && length(S)==1
- p = bounds(o, [], false);
+ varargout{1} = bounds(o, [], false);
elseif ismember(S(1).subs, {'deviate','length','isempty'})
- p = feval(S(1).subs, o);
- elseif ismember(S(1).subs, {'sample', 'density', 'densities', 'moments', 'admissible', 'bounds'})
- p = feval(S(1).subs, o , S(2).subs{:});
+ varargout{1} = feval(S(1).subs, o);
+ elseif ismember(S(1).subs, {'sample', 'density', 'densities', 'moments', 'admissible', 'bounds', 'data_for_plotting_density'})
+ [varargout{1:nargout}] = feval(S(1).subs, o , S(2).subs{:});
elseif ismember(S(1).subs, {'mean', 'median', 'variance', 'mode'})
if (length(S)==2 && isempty(S(2).subs)) || length(S)==1
- p = feval(S(1).subs, o);
+ varargout{1} = feval(S(1).subs, o);
else
- p = feval(S(1).subs, o , S(2).subs{:});
+ varargout{1} = feval(S(1).subs, o , S(2).subs{:});
end
else
error('dprior::subsref: unknown method (%s).', S(1).subs)
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index c991a120f4..d006d5fe9f 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -439,7 +439,7 @@ options_.occbin.smoother.status=false;
options_.silent_optimizer = false;
% Prior restrictions
-options_.prior_restrictions.status = 0;
+options_.prior_restrictions.status = false;
options_.prior_restrictions.routine = [];
options_.mode_compute = 5;
diff --git a/matlab/derror.m b/matlab/derror.m
new file mode 100644
index 0000000000..2c5b50b890
--- /dev/null
+++ b/matlab/derror.m
@@ -0,0 +1,38 @@
+function derror(varargin)
+
+% Throws an error message, without stacktrace under MATLAB.
+%
+% INPUTS
+% - id [char] 1×m array, error id (of the form 'xyz:abcdef')
+% - msg [char] error message
+%
+% OUTPUTS
+% None.
+%
+% REMARKS
+% - First argument is optional.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ if isoctave() || isdebug()
+ error(varargin{:})
+ else
+ throwAsCaller(MException(varargin{:}))
+ end
+
+end
diff --git a/matlab/dynare.m b/matlab/dynare.m
index ce973cd45e..cffe619875 100644
--- a/matlab/dynare.m
+++ b/matlab/dynare.m
@@ -185,6 +185,7 @@ nolog = ismember('nolog', varargin) || ismember('nolog', file_opts);
onlymacro = ismember('onlymacro', varargin) || ismember('onlymacro', file_opts);
onlyjson = ismember('onlyjson', varargin) || ismember('onlyjson', file_opts);
fast = ismember('fast', varargin) || ismember('fast', file_opts);
+verbose = ismember('verbose', varargin) || ismember('verbose', file_opts); % TODO: Find another name for the option.
% Start journal
diary off
@@ -196,6 +197,17 @@ if ~nolog
diary(logfile)
end
+if verbose
+ setenv('DYNARE_DEBUG_MODE', 'true')
+else
+ warning('off', 'backtrace');
+ if isoctave() || ~matlab_ver_less_than('9.13')
+ unsetenv('DYNARE_DEBUG_MODE');
+ else
+ setenv('DYNARE_DEBUG_MODE', 'false');
+ end
+end
+
if preprocessoroutput
fprintf(['Starting Dynare (version ' dynare_version() ').\n']);
fprintf('Calling Dynare with arguments: ');
diff --git a/matlab/estimation/PlotPosteriorDistributions.m b/matlab/estimation/PlotPosteriorDistributions.m
index 633ab1cf56..e60925155a 100644
--- a/matlab/estimation/PlotPosteriorDistributions.m
+++ b/matlab/estimation/PlotPosteriorDistributions.m
@@ -1,21 +1,18 @@
-function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_)
-% oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_)
-% plots posterior distributions
+function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, Prior, oo_)
+
+% Plot posterior densities
%
% INPUTS
-% estim_params_ [structure] information on estimated paramters
-% M_ [structure] information on model
-% options_ [structure] information on options
-% bayestopt_ [structure] information on priors
-% oo_ [structure] results
+% - estim_params_ [struct] description of the estimated paramters.
+% - M_ [struct] description of the model.
+% - options_ [struct] Dynare's options.
+% - Prior [dprior] object containing information about the prior distribution.
+% - oo_ [struct] results.
%
% OUTPUTS
-% oo_ [structure] updated results
-%
-% SPECIAL REQUIREMENTS
-% none
+% - oo_ [struct] updated results.
-% Copyright © 2005-2023 Dynare Team
+% Copyright © 2005-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -65,7 +62,7 @@ for i=1:npar
hh_fig=dyn_figure(options_.nodisplay, 'Name', figurename);
end
[nam,texnam] = get_the_name(i, TeX, M_, estim_params_, options_.varobs);
- [x2, f2, ~, ~, binf2, bsup2] = draw_prior_density(i, bayestopt_);
+ [x2, f2, ~, ~, binf2, bsup2] = Prior.data_for_plotting_density(i);
top2 = max(f2);
if i <= nvx
name = M_.exo_names{estim_params_.var_exo(i,1)};
diff --git a/matlab/estimation/check_bounds_and_definiteness_estimation.m b/matlab/estimation/check_bounds_and_definiteness_estimation.m
index a25478026c..fd3195d9a2 100644
--- a/matlab/estimation/check_bounds_and_definiteness_estimation.m
+++ b/matlab/estimation/check_bounds_and_definiteness_estimation.m
@@ -33,16 +33,17 @@ function [fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xpar
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-fval = [];
-exit_flag = 1;
-info = zeros(4,1);
-Q=[];
-H=[];
+fval = [];
+exit_flag = true;
+info = zeros(4,1);
+Q = [];
+H = [];
+
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if any(xparam1<bounds.lb)
k = find(xparam1(:) < bounds.lb);
fval = Inf;
- exit_flag = 0;
+ exit_flag = false;
info(1) = 41;
info(4) = sum((bounds.lb(k)-xparam1(k)).^2);
return
@@ -52,7 +53,7 @@ end
if any(xparam1>bounds.ub)
k = find(xparam1(:)>bounds.ub);
fval = Inf;
- exit_flag = 0;
+ exit_flag = false;
info(1) = 42;
info(4) = sum((xparam1(k)-bounds.ub(k)).^2);
return
@@ -62,10 +63,10 @@ Q = M_.Sigma_e;
H = M_.H;
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_,'calibrated_covariances')
- [Q_is_positive_definite, penalty] = ispd(Q(estim_params_.Sigma_e_entries_to_check_for_positive_definiteness,estim_params_.Sigma_e_entries_to_check_for_positive_definiteness));
+ [Q_is_positive_definite, penalty] = ispd(Q(estim_params_.Sigma_e_entries_to_check_for_positive_definiteness, estim_params_.Sigma_e_entries_to_check_for_positive_definiteness));
if ~Q_is_positive_definite
fval = Inf;
- exit_flag = 0;
+ exit_flag = false;
info(1) = 43;
info(4) = penalty;
return
@@ -75,7 +76,7 @@ if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_,'calibrated_covari
if ~correct_flag
penalty = sum(Q(estim_params_.calibrated_covariances.position).^2);
fval = Inf;
- exit_flag = 0;
+ exit_flag = false;
info(1) = 71;
info(4) = penalty;
return
@@ -85,10 +86,10 @@ if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_,'calibrated_covari
end
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_,'calibrated_covariances_ME')
- [H_is_positive_definite, penalty] = ispd(H(estim_params_.H_entries_to_check_for_positive_definiteness,estim_params_.H_entries_to_check_for_positive_definiteness));
+ [H_is_positive_definite, penalty] = ispd(H(estim_params_.H_entries_to_check_for_positive_definiteness, estim_params_.H_entries_to_check_for_positive_definiteness));
if ~H_is_positive_definite
fval = Inf;
- exit_flag = 0;
+ exit_flag = false;
info(1) = 44;
info(4) = penalty;
return
@@ -98,7 +99,7 @@ if ~issquare(H) || estim_params_.ncn || isfield(estim_params_,'calibrated_covari
if ~correct_flag
penalty = sum(H(estim_params_.calibrated_covariances_ME.position).^2);
fval = Inf;
- exit_flag = 0;
+ exit_flag = false;
info(1) = 72;
info(4) = penalty;
return
diff --git a/matlab/estimation/check_prior_bounds.m b/matlab/estimation/check_prior_bounds.m
index bf77952d63..d677f161a5 100644
--- a/matlab/estimation/check_prior_bounds.m
+++ b/matlab/estimation/check_prior_bounds.m
@@ -37,7 +37,10 @@ if ~isempty(outside_bound_pars)
end
error(['Initial value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0. If you used the mode_file-option, check whether your mode-file is consistent with the priors.'])
end
-inadmissible_inverse_gamma_values=find(bayestopt_.pshape==4 & xparam1 == 0);
+
+% TODO: Check the test below. If the inverse gamma is shifted, it probably doesn't make much sense... Instead, we should test whether the initial guess is equal to the natural lower bound of the inverse gamma (in the case where the prior mode is at the lower bound).
+% TODO: Should we not check the same for densities with an asymptote at the lower or upper natural boundary (eg for the beta distribution)?
+inadmissible_inverse_gamma_values=find(bayestopt_.pshape & xparam1 == 0);
if ~isempty(inadmissible_inverse_gamma_values)
for ii=1:length(inadmissible_inverse_gamma_values)
inadmissible_inverse_gamma_par_names{ii,1}=get_the_name(inadmissible_inverse_gamma_values(ii),0,M_,estim_params_,options_.varobs);
@@ -47,4 +50,4 @@ if ~isempty(inadmissible_inverse_gamma_values)
disp_string=[disp_string,', ',inadmissible_inverse_gamma_par_names{ii,:}];
end
error(['Initial value(s) of ', disp_string ,' is zero. This is not allowed when using an inverse gamma prior.\n'])
-end
\ No newline at end of file
+end
diff --git a/matlab/estimation/dsge_likelihood.m b/matlab/estimation/dsge_likelihood.m
index 265373e1a9..fb65dfbfd2 100644
--- a/matlab/estimation/dsge_likelihood.m
+++ b/matlab/estimation/dsge_likelihood.m
@@ -1,24 +1,23 @@
-function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,derivatives_info)
-% [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,oo_] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,oo_,derivatives_info)
-% Evaluates the negative of the posterior kernel of a DSGE model using the specified
-% kalman_algo; the resulting posterior includes the 2*pi constant of the
+function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,derivatives_info)
+
+% Evaluates the negative of the posterior kernel of a DSGE model using the specified kalman_algo; the resulting posterior includes the 2*pi constant of the
% likelihood function
%
% INPUTS
-% - xparam1 [double] current values for the estimated parameters.
-% - dataset_ [structure] dataset after transformations
-% - dataset_info [structure] storing informations about the
-% sample; not used but required for interface
-% - options_ [structure] Matlab's structure describing the current options
-% - M_ [structure] Matlab's structure describing the model
-% - estim_params_ [structure] characterizing parameters to be estimated
-% - bayestopt_ [structure] describing the priors
-% - BoundsInfo [structure] containing prior bounds
-% - dr [structure] Reduced form model.
-% - endo_steady_state [vector] steady state value for endogenous variables
-% - exo_steady_state [vector] steady state value for exogenous variables
-% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
-% - derivatives_info [structure] derivative info for identification
+% - xparam1 [double] current values for the estimated parameters.
+% - dataset_ [struct] dataset after transformations
+% - dataset_info [struct] storing informations about the sample; not used but required for interface
+% - options_ [struct] Matlab's structure describing the current options
+% - M_ [struct] Matlab's structure describing the model
+% - Prior [dprior] Definition of the prior density
+% - estim_params_ [struct] characterizing parameters to be estimated
+% - bayestopt_ [struct] describing the priors
+% - BoundsInfo [struct] containing prior bounds
+% - dr [struct] Reduced form model.
+% - endo_steady_state [double] steady state value for endogenous variables
+% - exo_steady_state [double] steady state value for exogenous variables
+% - exo_det_steady_state [double] steady state value for exogenous deterministic variables
+% - derivatives_info [struct] derivative info for identification
%
% OUTPUTS
% - fval [double] scalar, value of minus the likelihood or posterior kernel.
@@ -35,7 +34,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,baye
%
% This function is called by: dynare_estimation_1, mode_check
% This function calls: dynare_resolve, lyapunov_symm, lyapunov_solver, compute_Pinf_Pstar, kalman_filter_d, missing_observations_kalman_filter_d,
-% univariate_kalman_filter_d, kalman_steady_state, get_perturbation_params_deriv, kalman_filter, missing_observations_kalman_filter, univariate_kalman_filter, priordens
+% univariate_kalman_filter_d, kalman_steady_state, get_perturbation_params_deriv, kalman_filter, missing_observations_kalman_filter, univariate_kalman_filter
% Copyright © 2004-2024 Dynare Team
%
@@ -57,17 +56,17 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,baye
% Initial author: stephane DOT adjemian AT univ DASH lemans DOT FR
% Initialization of the returned variables and others...
-fval = [];
+fval = [];
SteadyState = [];
trend_coeff = [];
-exit_flag = 1;
-info = zeros(4,1);
+exit_flag = false;
+info = zeros(4,1);
if options_.analytic_derivation
- DLIK = NaN(1,length(xparam1));
+ DLIK = NaN(1,length(xparam1));
else
- DLIK = [];
+ DLIK = [];
end
-Hess = [];
+Hess = [];
% Ensure that xparam1 is a column vector.
% (Don't do the transformation if xparam1 is empty, otherwise it would become a
@@ -90,19 +89,19 @@ if analytic_derivation
end
if nargout==1
- analytic_derivation=0;
+ analytic_derivation = 0;
end
if analytic_derivation
- kron_flag=options_.analytic_derivation_mode;
+ kron_flag = options_.analytic_derivation_mode;
end
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
-M_ = set_all_parameters(xparam1,estim_params_,M_);
+M_ = set_all_parameters(xparam1, estim_params_, M_);
-[fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, BoundsInfo);
+[fval, info, exit_flag, Q, H] = check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, BoundsInfo);
if info(1)
return
end
@@ -135,7 +134,7 @@ if options_.occbin.likelihood.status
occbin_.info= {options_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, M_, occbin_options, TTx, RRx, CCx,T0,R0};
else
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
- [T,R,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
+ [T, R, SteadyState, info, dr, M_.params] = dynare_resolve(M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, 'restrict');
occbin_.status = false;
end
@@ -153,30 +152,30 @@ if info(1)
else
info(4) = info(2);
end
- exit_flag = 0;
+ exit_flag = false;
if analytic_derivation
- DLIK=ones(length(xparam1),1);
+ DLIK = ones(length(xparam1), 1);
end
return
else
fval = Inf;
info(4) = 0.1;
- exit_flag = 0;
+ exit_flag = false;
if analytic_derivation
- DLIK=ones(length(xparam1),1);
+ DLIK = ones(length(xparam1), 1);
end
return
end
end
% check endogenous prior restrictions
-info=endogenous_prior_restrictions(T,R,M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
+info = endogenous_prior_restrictions(T, R, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if info(1)
fval = Inf;
- info(4)=info(2);
- exit_flag = 0;
+ info(4) = info(2);
+ exit_flag = false;
if analytic_derivation
- DLIK=ones(length(xparam1),1);
+ DLIK = ones(length(xparam1), 1);
end
return
end
@@ -804,10 +803,10 @@ likelihood = LIK;
% ------------------------------------------------------------------------------
if analytic_derivation
if full_Hess
- [lnprior, dlnprior, d2lnprior] = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
+ [lnprior, dlnprior, d2lnprior] = Prior.density(xparam1);
Hess = Hess - d2lnprior;
else
- [lnprior, dlnprior] = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
+ [lnprior, dlnprior] = Prior.density(xparam1);
end
if no_DLIK==0
DLIK = DLIK - dlnprior';
@@ -818,7 +817,7 @@ if analytic_derivation
Hess = dlik'*dlik;
end
else
- lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
+ lnprior = Prior.density(xparam1);
end
if options_.endogenous_prior==1
diff --git a/matlab/estimation/dsge_var_likelihood.m b/matlab/estimation/dsge_var_likelihood.m
index 39eace525f..aa74e6cd48 100644
--- a/matlab/estimation/dsge_var_likelihood.m
+++ b/matlab/estimation/dsge_var_likelihood.m
@@ -1,45 +1,40 @@
function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_u_tilde,iXX,prior] = dsge_var_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
-% [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_u_tilde,iXX,prior] = dsge_var_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+
% Evaluates the posterior kernel of the BVAR-DSGE model.
%
% INPUTS
-% o xparam1 [double] Vector of model's parameters.
-% o gend [integer] Number of observations (without conditionning observations for the lags).
-% o dataset_ [dseries] object storing the dataset
-% o dataset_info [structure] storing informations about the sample.
-% o options_ [structure] describing the options
-% o M_ [structure] decribing the model
-% o estim_params_ [structure] characterizing parameters to be estimated
-% o bayestopt_ [structure] describing the priors
-% o BoundsInfo [structure] containing prior bounds
-% o dr [structure] Reduced form model.
-% o endo_steady_state [vector] steady state value for endogenous variables
-% o exo_steady_state [vector] steady state value for exogenous variables
-% o exo_det_steady_state [vector] steady state value for exogenous deterministic variables
+% - xparam1 [double] Vector of model's parameters.
+% - gend [integer] Number of observations (without conditionning observations for the lags).
+% - dataset_ [dseries] object storing the dataset
+% - dataset_info [struct] storing informations about the sample.
+% - options_ [struct] describing the options
+% - M_ [struct] decribing the model
+% - Prior [dprior] definition of the prior density
+% - estim_params_ [struct] characterizing parameters to be estimated
+% - bayestopt_ [struct] describing the priors
+% - BoundsInfo [struct] containing prior bounds
+% - dr [struct] Reduced form model.
+% - endo_steady_state [double] steady state value for endogenous variables
+% - exo_steady_state [double] steady state value for exogenous variables
+% - exo_det_steady_state [double] steady state value for exogenous deterministic variables
%
% OUTPUTS
-% o fval [double] Value of the posterior kernel at xparam1.
-% o info [integer] Vector of informations about the penalty.
-% o exit_flag [integer] Zero if the function returns a penalty, one otherwise.
-% o grad [double] place holder for gradient of the likelihood
-% currently not supported by dsge_var
-% o hess [double] place holder for hessian matrix of the likelihood
-% currently not supported by dsge_var
-% o SteadyState [double] Steady state vector possibly recomputed
-% by call to dynare_resolve()
-% o trend_coeff [double] place holder for trend coefficients,
-% currently not supported by dsge_var
-% o PHI_tilde [double] Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1);
-% formula (28), DS (2004)
-% o SIGMA_u_tilde [double] Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1),
-% formula (29), DS (2004)
-% o iXX [double] inv(lambda*T*Gamma_XX^*+ X'*X)
-% o prior [double] a matlab structure describing the dsge-var prior
-% - SIGMA_u_star: prior covariance matrix, formula (23), DS (2004)
-% - PHI_star: prior autoregressive matrices, formula (22), DS (2004)
-% - ArtificialSampleSize: number of artificial observations from the prior (T^* in DS (2004))
-% - DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
-% - iGXX_star: theoretical covariance of regressors ({\Gamma_{XX}^*}^{-1} in DS (2004))
+% - fval [double] Value of the posterior kernel at xparam1.
+% - info [integer] Vector of informations about the penalty.
+% - exit_flag [integer] Zero if the function returns a penalty, one otherwise.
+% - grad [double] place holder for gradient of the likelihood, currently not supported by dsge_var
+% - hess [double] place holder for hessian matrix of the likelihood, currently not supported by dsge_var
+% - SteadyState [double] Steady state vector possibly recomputed by call to dynare_resolve()
+% - trend_coeff [double] place holder for trend coefficients, currently not supported by dsge_var
+% - PHI_tilde [double] Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1); formula (28), DS (2004)
+% - SIGMA_u_tilde [double] Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1); formula (29), DS (2004)
+% - iXX [double] inv(lambda*T*Gamma_XX^*+ X'*X)
+% - prior [double] a matlab structure describing the dsge-var prior
+% - SIGMA_u_star: prior covariance matrix, formula (23), DS (2004)
+% - PHI_star: prior autoregressive matrices, formula (22), DS (2004)
+% - ArtificialSampleSize: number of artificial observations from the prior (T^* in DS (2004))
+% - DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
+% - iGXX_star: theoretical covariance of regressors ({\Gamma_{XX}^*}^{-1} in DS (2004))
%
% ALGORITHMS
% Follows the computations outlined in Del Negro/Schorfheide (2004):
@@ -49,7 +44,7 @@ function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_
% SPECIAL REQUIREMENTS
% None.
-% Copyright © 2006-2023Dynare Team
+% Copyright © 2006-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -284,8 +279,7 @@ if imag(lik)~=0
end
% Add the (logged) prior density for the dsge-parameters.
-lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
-fval = (lik-lnprior);
+fval = lik-Prior.density(xparam1);
if isnan(fval)
fval = Inf;
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 0962b4c884..650e1eef89 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -120,7 +120,7 @@ else
objective_function = str2func('dsge_var_likelihood');
end
-[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = ...
+[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, Prior, estim_params_, bayestopt_, bounds] = ...
dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
if options_.dsge_var
@@ -166,7 +166,7 @@ end
%
try
- oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);
+ initial_estimation_checks(objective_function, xparam1, dataset_, dataset_info, M_, Prior, estim_params_, options_, bayestopt_, bounds, oo_);
catch e % if check fails, provide info on using calibration if present
if estim_params_.full_calibration_detected %calibrated model present and no explicit starting values
skipline(1);
@@ -237,7 +237,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
for optim_iter = 1:length(optimizer_vec)
current_optimizer = optimizer_vec{optim_iter};
- [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
if isnumeric(current_optimizer)
@@ -266,8 +266,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
- hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
- end
+ hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, Prior, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ end
hh = reshape(hh, nx, nx);
elseif isnumeric(current_optimizer) && isequal(current_optimizer,5)
% other numerical hessian options available with optimizer
@@ -322,7 +322,7 @@ if options_.mode_check.status && ~options_.mh_posterior_mode_estimation && ~issm
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
mode_check(objective_function,xparam1,hh,options_,M_,estim_params_,bayestopt_,bounds,false,...
- dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ dataset_, dataset_info, options_, M_, Prior, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_.analytic_derivation = ana_deriv_old;
end
@@ -369,7 +369,7 @@ if ~issmc(options_) && any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode
estim_params_nbr = size(xparam1,1);
if ispd(invhess)
log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh));
- likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
else
oo_.MarginalDensity.LaplaceApproximation = NaN;
@@ -439,7 +439,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
if options_.mh_replic
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
- posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString);
+ posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_, Prior, estim_params_,bayestopt_,oo_,dispString);
options_.analytic_derivation = ana_deriv_old;
end
end
@@ -478,7 +478,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
% Store posterior statistics by parameter name
oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, prior_dist_names);
if ~options_.nograph
- oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
+ oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, Prior, oo_);
end
% Store posterior mean in a vector and posterior variance in
% a matrix
diff --git a/matlab/estimation/dynare_estimation_init.m b/matlab/estimation/dynare_estimation_init.m
index 0f15716239..6229e95215 100644
--- a/matlab/estimation/dynare_estimation_init.m
+++ b/matlab/estimation/dynare_estimation_init.m
@@ -1,4 +1,4 @@
-function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, bounds] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
+function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, Prior, estim_params_,bayestopt_, bounds] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
% Performs initialization tasks before estimation or global sensitivity analysis
%
@@ -91,31 +91,30 @@ end
% Check the perturbation order for pruning (k order perturbation based nonlinear filters are not yet implemented for k>3).
if options_.order>3 && options_.particle.pruning
- error('Higher order nonlinear filters are not compatible with pruning option.')
+ derror('Estimation:nonlinear', 'Higher order nonlinear filters are not compatible with pruning option.')
end
% analytical derivation is not yet available for kalman_filter_fast
if options_.analytic_derivation && options_.fast_kalman_filter
- error(['estimation option conflict: analytic_derivation isn''t available ' ...
- 'for fast_kalman_filter'])
+ derror('Estimation:kalman', 'estimation option conflict: analytic_derivation isn''t available for fast_kalman_filter.')
end
% fast kalman filter is only available with kalman_algo == 0,1,3
if options_.fast_kalman_filter && ~ismember(options_.kalman_algo, [0,1,3])
- error(['estimation option conflict: fast_kalman_filter is only available ' ...
- 'with kalman_algo = 0, 1 or 3'])
+ derror('Estimation:kalman', 'estimation option conflict: fast_kalman_filter is only available with kalman_algo = 0, 1 or 3.')
end
% Set options_.lik_init equal to 3 if diffuse filter is used or kalman_algo refers to a diffuse filter algorithm.
if isequal(options_.diffuse_filter,1) || (options_.kalman_algo>2)
if isequal(options_.lik_init,2)
- error('options diffuse_filter, lik_init and/or kalman_algo have contradictory settings')
+ derror('Estimation:kalman', 'options diffuse_filter, lik_init and/or kalman_algo have contradictory settings.')
else
options_.lik_init = 3;
end
end
-options_=select_qz_criterium_value(options_);
+% Set the value of qz_criterium
+options_ = select_qz_criterium_value(options_);
% Set options related to filtered variables.
if isequal(options_.filtered_vars,0) && ~isempty(options_.filter_step_ahead) %require filtered var because filter_step_ahead was set
@@ -135,23 +134,40 @@ else
M_.dname = dname;
end
-% Set priors over the estimated parameters.
-if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
- [xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
-end
+estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+estim_params_ = set_nvn_observable_correspondence(estim_params_, M_.endo_names, options_.varobs);
+
+estim_params_ = set_corrn_observable_correspondence(estim_params_, M_.endo_names, options_.varobs);
+
+M_ = initialize_measurement_errors_covariance_matrix(estim_params_, M_, options_);
-if ~isempty(bayestopt_) && any(bayestopt_.pshape==0) && any(bayestopt_.pshape~=0)
- error('Estimation must be either fully ML or fully Bayesian. Maybe you forgot to specify a prior distribution.')
+bayesianestimation = isbayesian(estim_params_);
+
+if bayesianestimation
+ Prior = set_prior(estim_params_, M_, options_);
+ estim_params_ = set_jscale(estim_params_, true);
+ estim_params_ = set_parameter_names(estim_params_, M_.param_names, M_.exo_names, M_.endo_names, options_.varobs);
+else
+ estim_params_ = set_parameter_names(estim_params_, M_.param_names, M_.exo_names, M_.endo_names, options_.varobs);
+ Prior = dprior();
end
+
+xparam1 = initialguess(estim_params_, options_.initialize_estimated_parameters_with_the_prior_mode, Prior);
+
+[lb, ub] = get_user_defined_boundaries(estim_params_, M_, options_);
+bounds.lb = lb;
+bounds.ub = ub;
+
% Check if a _prior_restrictions.m file exists
if isfile(sprintf('%s_prior_restrictions.m', M_.fname))
- options_.prior_restrictions.status = 1;
+ options_.prior_restrictions.status = true;
options_.prior_restrictions.routine = str2func([M_.fname '_prior_restrictions']);
end
% Check that the provided mode_file is compatible with the current estimation settings.
if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0) && ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
- [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_);
+ [xparam1, hh] = check_mode_file(xparam1, hh, options_.mode_compute>0, options_.mode_file, estim_params_.names);
end
%check for calibrated covariances before updating parameters
@@ -174,26 +190,21 @@ if options_.use_calibration_initialization %set calibration as starting values
end
end
-if ~isempty(bayestopt_) && all(bayestopt_.pshape==0) && any(isnan(xparam1))
+if estim_params_.estimation && ~bayesianestimation && any(isnan(xparam1))
error('ML estimation requires all estimated parameters to be initialized, either in an estimated_params or estimated_params_init-block ')
end
if ~isempty(estim_params_) && ~(all(strcmp(fieldnames(estim_params_),'full_calibration_detected')) || (isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0))
- if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
+ if bayesianestimation
% Plot prior densities.
if ~options_.nograph && options_.plot_priors
- plot_priors(bayestopt_,M_,estim_params_,options_)
+ plot_priors(Prior, M_, estim_params_, options_)
end
- % Set prior bounds
- bounds = prior_bounds(bayestopt_, options_.prior_trunc);
- bounds.lb = max(bounds.lb,lb);
- bounds.ub = min(bounds.ub,ub);
- else % estimated parameters but no declared priors
- % No priors are declared so Dynare will estimate the model by
- % maximum likelihood with inequality constraints for the parameters.
+ else
+ % estimated parameters but no declared priors
+ % No priors are declared so Dynare will estimate the model by
+ % maximum likelihood
options_.mh_replic = 0;% No metropolis.
- bounds.lb = lb;
- bounds.ub = ub;
end
% Test if initial values of the estimated parameters are all between the prior lower and upper bounds.
if options_.use_calibration_initialization
@@ -473,7 +484,7 @@ end
estim_params_= get_matrix_entries_for_psd_check(M_,estim_params_);
if options_.load_results_after_load_mh
- if ~isfile(sprintf('%s%sOutput%s%s_results.mat', M_.dname, filesep, filesep, M°.fname))
+ if ~isfile(sprintf('%s%sOutput%s%s_results.mat', M_.dname, filesep, filesep, M_.fname))
fprintf('\ndynare_estimation_init:: You specified the load_results_after_load_mh, but no _results.mat-file\n')
fprintf('dynare_estimation_init:: was found. Results will be recomputed.\n')
options_.load_results_after_load_mh=0;
@@ -488,73 +499,73 @@ end
%% set heteroskedastic shocks
if options_.heteroskedastic_filter
if options_.fast_kalman_filter
- error('estimation option conflict: "heteroskedastic_filter" incompatible with "fast_kalman_filter"')
- end
- if options_.analytic_derivation
- error(['estimation option conflict: analytic_derivation isn''t available ' ...
- 'for heteroskedastic_filter'])
- end
-
- M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs+1);
- M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs+1);
-
- for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig)
- v = M_.heteroskedastic_shocks.Qvalue_orig(k);
- temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
- temp_periods=temp_periods(temp_periods>=options_.first_obs);
- M_.heteroskedastic_shocks.Qvalue(v.exo_id, temp_periods-(options_.first_obs-1)) = v.value^2;
- end
- for k=1:length(M_.heteroskedastic_shocks.Qscale_orig)
- v = M_.heteroskedastic_shocks.Qscale_orig(k);
- temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
- temp_periods=temp_periods(temp_periods>=options_.first_obs);
- M_.heteroskedastic_shocks.Qscale(v.exo_id, temp_periods-(options_.first_obs-1)) = v.scale^2;
- end
-
- if any(any(~isnan(M_.heteroskedastic_shocks.Qvalue) & ~isnan(M_.heteroskedastic_shocks.Qscale)))
- fprintf('\ndynare_estimation_init: With the option "heteroskedastic_shocks" you cannot define\n')
- fprintf('dynare_estimation_init: the scale and the value for the same shock \n')
- fprintf('dynare_estimation_init: in the same period!\n')
- error('Scale and value defined for the same shock in the same period with "heteroskedastic_shocks".')
- end
-end
-
-if (options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter) || (options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
- if isempty(options_.occbin.likelihood.IVF_shock_observable_mapping)
- options_.occbin.likelihood.IVF_shock_observable_mapping=find(diag(M_local.Sigma_e)~=0);
- else
- zero_var_shocks=find(diag(M_local.Sigma_e)==0);
- if any(ismember(options_.occbin.likelihood.IVF_shock_observable_mapping,zero_var_shocks))
- error('IVF-filter: an observable is mapped to a zero variance shock.')
- end
- end
- if ~isequal(M_.H,0)
+ error('estimation option conflict: "heteroskedastic_filter" incompatible with "fast_kalman_filter"')
+ end
+ if options_.analytic_derivation
+ error(['estimation option conflict: analytic_derivation isn''t available ' ...
+ 'for heteroskedastic_filter'])
+ end
+
+ M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs+1);
+ M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs+1);
+
+ for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig)
+ v = M_.heteroskedastic_shocks.Qvalue_orig(k);
+ temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
+ temp_periods=temp_periods(temp_periods>=options_.first_obs);
+ M_.heteroskedastic_shocks.Qvalue(v.exo_id, temp_periods-(options_.first_obs-1)) = v.value^2;
+ end
+ for k=1:length(M_.heteroskedastic_shocks.Qscale_orig)
+ v = M_.heteroskedastic_shocks.Qscale_orig(k);
+ temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
+ temp_periods=temp_periods(temp_periods>=options_.first_obs);
+ M_.heteroskedastic_shocks.Qscale(v.exo_id, temp_periods-(options_.first_obs-1)) = v.scale^2;
+ end
+
+ if any(any(~isnan(M_.heteroskedastic_shocks.Qvalue) & ~isnan(M_.heteroskedastic_shocks.Qscale)))
+ fprintf('\ndynare_estimation_init: With the option "heteroskedastic_shocks" you cannot define\n')
+ fprintf('dynare_estimation_init: the scale and the value for the same shock \n')
+ fprintf('dynare_estimation_init: in the same period!\n')
+ error('Scale and value defined for the same shock in the same period with "heteroskedastic_shocks".')
+ end
+ end
+
+ if (options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter) || (options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
+ if isempty(options_.occbin.likelihood.IVF_shock_observable_mapping)
+ options_.occbin.likelihood.IVF_shock_observable_mapping=find(diag(M_local.Sigma_e)~=0);
+ else
+ zero_var_shocks=find(diag(M_local.Sigma_e)==0);
+ if any(ismember(options_.occbin.likelihood.IVF_shock_observable_mapping,zero_var_shocks))
+ error('IVF-filter: an observable is mapped to a zero variance shock.')
+ end
+ end
+ if ~isequal(M_.H,0)
error('IVF-filter: Measurement erros are not allowed with the inversion filter.')
end
-end
-
-if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
- if ~isempty(options_.nk)
- fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead and filtered_variables. Disabling the option.\n')
- options_.nk=[];
- options_.filter_step_ahead=[];
- options_.filtered_vars=0;
- end
- if options_.filter_covariance
- fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')
- options_.filter_covariance=false;
- end
- if options_.smoothed_state_uncertainty
- fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n')
- options_.smoothed_state_uncertainty=false;
- end
-end
-
-if options_.occbin.smoother.status || options_.occbin.likelihood.status
- if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
- fprintf('dynare_estimation_init: OccBin smoother/filter: previous shocks(surprise) block detected. Deleting its content.\n')
- options_.occbin.simul.SHOCKS=zeros(1,M_.exo_nbr);
- options_.occbin.simul.exo_pos=1:M_.exo_nbr;
- M_.surprise_shocks=[];
- end
-end
+ end
+
+ if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
+ if ~isempty(options_.nk)
+ fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead and filtered variables. Disabling the option.\n')
+ options_.nk=[];
+ options_.filter_step_ahead=[];
+ options_.filtered_vars=0;
+ end
+ if options_.filter_covariance
+ fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')
+ options_.filter_covariance=false;
+ end
+ if options_.smoothed_state_uncertainty
+ fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n')
+ options_.smoothed_state_uncertainty=false;
+ end
+ end
+
+ if options_.occbin.smoother.status || options_.occbin.likelihood.status
+ if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
+ fprintf('dynare_estimation_init: OccBin smoother/filter: previous shocks(surprise) block detected. Deleting its content.\n')
+ options_.occbin.simul.SHOCKS=zeros(1,M_.exo_nbr);
+ options_.occbin.simul.exo_pos=1:M_.exo_nbr;
+ M_.surprise_shocks=[];
+ end
+ end
diff --git a/matlab/estimation/get_user_defined_boundaries.m b/matlab/estimation/get_user_defined_boundaries.m
new file mode 100644
index 0000000000..1d2bb517fc
--- /dev/null
+++ b/matlab/estimation/get_user_defined_boundaries.m
@@ -0,0 +1,77 @@
+function [lb, ub] = get_user_defined_boundaries(estim_params_, M_, options_)
+
+% Get optimization boundaries on the estimated parameters.
+%
+% INPUTS
+% - estim_params_ [struct] characterizing parameters to be estimated.
+% - M_ [struct] characterizing the model.
+% - options_ [struct] characterizing the options.
+%
+% OUTPUTS
+% - lb [double] n×1 vector, lower bounds.
+% - ub [double] n×1 vector, upper bounds.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+ nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+ mm = estim_params_.mm; % Total number of estimated parameters
+
+ lb = -Inf(mm, 1);
+ ub = Inf(mm, 1);
+
+ if estim_params_.nvx
+ lb(1:nn(1)) = max(estim_params_.var_exo(:,3), 0);
+ ub(1:nn(1)) = estim_params_.var_exo(:,4);
+ check_interval(lb(1:nn(1)), ub(1:nn(1)), 1:nn(1), estim_params_.names)
+ end
+
+ if estim_params_.nvn
+ lb(nn(1)+1:nn(2)) = max(estim_params_.var_endo(:,3), 0);
+ ub(nn(1)+1:nn(2)) = estim_params_.var_endo(:,4);
+ check_interval(lb(nn(1)+1:nn(2)), ub(nn(1)+1:nn(2)), nn(1)+1:nn(2), estim_params_.names)
+ end
+
+ if estim_params_.ncx
+ lb(nn(2)+1:nn(3)) = min(max(estim_params_.corrx(:,4), -1), 1);
+ ub(nn(2)+1:nn(3)) = max(min(estim_params_.corrx(:,5), 1), -1);
+ check_interval(lb(nn(2)+1:nn(3)), ub(nn(2)+1:nn(3)), nn(2)+1:nn(3), estim_params_.names)
+ end
+
+ if estim_params_.ncn
+ lb(nn(3)+1:nn(4)) = min(max(estim_params_.corrn(:,4), -1), 1);
+ ub(nn(3)+1:nn(4)) = max(min(estim_params_.corrn(:,5), 1), -1);
+ check_interval(lb(nn(3)+1:nn(4)), ub(nn(3)+1:nn(4)), nn(3)+1:nn(4), estim_params_.names)
+ end
+
+ if estim_params_.np
+ lb(nn(4)+1:mm) = estim_params_.param_vals(:,3);
+ ub(nn(4)+1:mm) = estim_params_.param_vals(:,4);
+ check_interval(lb(nn(4)+1:mm), ub(nn(4)+1:mm), nn(4)+1:mm, estim_params_.names)
+ end
+
+ function check_interval(a, b, idx, names)
+ problem = a>b;
+ if any(problem)
+ names = names(idx(problem));
+ derror('Estimation:wrongBoundaries', sprintf('Check the declared bounds on (%s), the lower bound must be smaller than the upper bound.', regexprep(sprintf('%s, ', names{:}), ', $', '')))
+ end
+ end
+
+end
diff --git a/matlab/estimation/initial_estimation_checks.m b/matlab/estimation/initial_estimation_checks.m
index bab50ffecd..1dcf31cfac 100644
--- a/matlab/estimation/initial_estimation_checks.m
+++ b/matlab/estimation/initial_estimation_checks.m
@@ -166,8 +166,8 @@ end
[oo_.steady_state] = check_steady_state_changes_parameters(M_,estim_params_,oo_,options_, [options_.diffuse_filter==0 options_.diffuse_filter==0] );
% check and display warning if negative values of stderr or corr params are outside unit circle for Bayesian estimation
-if any(bayestopt_.pshape)
- check_prior_stderr_corr(estim_params_,bayestopt_);
+if ~isempty(Prior)
+ check_prior_stderr_corr(estim_params_, Prior);
end
% display warning if some parameters are still NaN
@@ -199,13 +199,11 @@ end
% Evaluate the likelihood.
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation=0;
-if ~isequal(options_.mode_compute,11) || ...
- (isequal(options_.mode_compute,11) && isequal(options_.order,1))
- %shut off potentially automatic switch to diffuse filter for the
- %purpose of checking stochastic singularity
- use_univariate_filters_if_singularity_is_detected_old=options_.use_univariate_filters_if_singularity_is_detected;
- options_.use_univariate_filters_if_singularity_is_detected=0;
- [fval,info] = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+if ~isequal(options_.mode_compute,11) || (isequal(options_.mode_compute,11) && isequal(options_.order,1))
+ % Shut off potentially automatic switch to diffuse filter for the purpose of checking stochastic singularity
+ use_univariate_filters_if_singularity_is_detected_old = options_.use_univariate_filters_if_singularity_is_detected;
+ options_.use_univariate_filters_if_singularity_is_detected = false;
+ [fval, info] = feval(objective_function, xparam1, dataset_, dataset_info, options_, M_, Prior, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if info(1)==50
fprintf('\ninitial_estimation_checks:: The forecast error variance in the multivariate Kalman filter became singular.\n')
fprintf('initial_estimation_checks:: This is often a sign of stochastic singularity, but can also sometimes happen by chance\n')
@@ -217,8 +215,8 @@ if ~isequal(options_.mode_compute,11) || ...
message=get_error_message(info,options_);
error('initial_estimation_checks:: %s.',message)
end
- %reset options
- options_.use_univariate_filters_if_singularity_is_detected=use_univariate_filters_if_singularity_is_detected_old;
+ % Reset options
+ options_.use_univariate_filters_if_singularity_is_detected = use_univariate_filters_if_singularity_is_detected_old;
else
info=0;
fval = 0;
diff --git a/matlab/estimation/initialguess.m b/matlab/estimation/initialguess.m
new file mode 100644
index 0000000000..020518b22b
--- /dev/null
+++ b/matlab/estimation/initialguess.m
@@ -0,0 +1,80 @@
+function [xparam1, estim_params_] = initialguess(estim_params_, usepriormode, Prior)
+
+% Return the initial guess for the estimated parameters
+%
+% INPUTS
+% - estim_params_ [struct] characterizing parameters to be estimated
+% - Prior [dprior] definition of the prior (if any)
+%
+% OUTPUTS
+% - xparam1 [double] mm×1 vector, initial guess
+% - estim_params_ [struct] characterizing parameters to be estimated (updated)
+%
+% REMARKS
+% - This routine should be called after the instanciation of the dprior object
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+ nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+ mm = estim_params_.mm; % Total number of estimated parameters
+
+ xparam1 = NaN(mm, 1);
+
+ if estim_params_.nvx
+ xparam1(1:nn(1)) = estim_params_.var_exo(:,2);
+ end
+
+ if estim_params_.nvn
+ xparam1(nn(1)+1:nn(2)) = estim_params_.var_endo(:,2);
+ end
+
+ if estim_params_.ncx
+ xparam1(nn(2)+1:nn(3)) = estim_params_.corrx(:,3);
+ end
+
+ if estim_params_.ncn
+ xparam1(nn(3)+1:nn(4)) = estim_params_.corrn(:,3);
+ end
+
+ if estim_params_.np
+ xparam1(nn(4)+1:mm) = estim_params_.param_vals(:,2);
+ end
+
+ if ~Prior.isempty()
+ declared_initial_guess = ~isnan(xparam1);
+ if ~usepriormode
+ % Initialize the estimated parameters with the prior mean (unless an initial guess has been explicitly defined).
+ xparam1(~declared_initial_guess) = Prior.p1(~declared_initial_guess);
+ else
+ % Initialize the estimated parameters with the prior mode (unless an initial guess has been explicitly defined).
+ xparam1(~declared_initial_guess) = Prior.p5(~declared_initial_guess);
+ no_prior_mode = isnan(xparam1);
+ if any(no_prior_mode)
+ % If the prior mode does not exist, use the prior mean as an initial guess.
+ xparam1(no_prior_mode) = Prior.p1(no_prior_mode);
+ end
+ end
+ end
+
+ if any(isnan(xparam1))
+ derror('Estimation:initialguess', 'Initial guess for all the estimated parameters must be provided.');
+ end
+
+end
diff --git a/matlab/estimation/isbayesian.m b/matlab/estimation/isbayesian.m
new file mode 100644
index 0000000000..46804e5e30
--- /dev/null
+++ b/matlab/estimation/isbayesian.m
@@ -0,0 +1,150 @@
+function c = isbayesian(estim_params_)
+
+% Return true iff bayesian approach is used.
+%
+% INPUTS
+% - estim_params_ [struct] Output of the preprocessor related to the estimation.
+%
+% OUTPUTS
+% - c [logical] scalar.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ c = false;
+
+ b = false(5, 1);
+ i = true(5, 1);
+
+ [b, i] = ispriorinblock(1, estim_params_.var_exo(:,5), b, i);
+ [b, i] = ispriorinblock(2, estim_params_.var_endo(:,5), b, i);
+ [b, i] = ispriorinblock(3, estim_params_.corrx(:,6), b, i);
+ [b, i] = ispriorinblock(4, estim_params_.corrn(:,6), b, i);
+ [b, i] = ispriorinblock(5, estim_params_.param_vals(:,5), b, i);
+
+ b = b(i);
+
+ if all(b)
+ c = true;
+ elseif any(b)
+ throwAsCaller(MException('Prior:specificationError', 'All estimated parameters must have a prior or none must have a prior (ML).'));
+ end
+
+ return % --*-- Unit tests --*--
+
+
+ %@test:1
+ estimparams.var_exo = zeros(7, 10);
+ estimparams.var_exo(:,5) = [1;2;3;4;5;6;8];
+ estimparams.var_endo = zeros(7, 10);
+ estimparams.var_endo(:,5) = [1;2;3;4;5;6;8];
+ estimparams.corrx = zeros(7, 10);
+ estimparams.corrx(:,6) = [1;2;3;4;5;6;8];
+ estimparams.corrn = zeros(7, 10);
+ estimparams.corrn(:,6) = [1;2;3;4;5;6;8];
+ estimparams.param_vals = zeros(7, 10);
+ estimparams.param_vals(:,5) = [1;2;3;4;5;6;8];
+
+ try
+ b = isbayesian(estimparams);
+ t(1) = true;
+ catch
+ t(1) = false;
+ end
+
+ if t(1)
+ t(2) = b;
+ end
+
+ T = all(t);
+ %@eof:1
+
+
+ %@test:2
+ estimparams.var_exo = zeros(7, 10);
+ estimparams.var_endo = zeros(7, 10);
+ estimparams.var_endo(:,5) = [1;2;3;4;5;6;8];
+ estimparams.corrx = zeros(7, 10);
+ estimparams.corrx(:,6) = [1;2;3;4;5;6;8];
+ estimparams.corrn = zeros(7, 10);
+ estimparams.corrn(:,6) = [1;2;3;4;5;6;8];
+ estimparams.param_vals = zeros(7, 10);
+ estimparams.param_vals(:,5) = [1;2;3;4;5;6;8];
+
+ try
+ b = isbayesian(estimparams);
+ t(1) = false;
+ catch
+ t(1) = true;
+ end
+
+ T = all(t);
+ %@eof:2
+
+
+ %@test:3
+ estimparams.var_exo = zeros(7, 10);
+ estimparams.var_exo(:,5) = [1;2;3;4;5;6;8];
+ estimparams.param_vals = zeros(7, 10);
+ estimparams.param_vals(:,5) = [1;2;3;4;5;6;8];
+
+ try
+ b = isbayesian(estimparams);
+ t(1) = true;
+ catch
+ t(1) = false;
+ end
+
+ if t(1)
+ t(2) = b;
+ end
+
+ T = all(t);
+ %@eof:3
+
+
+ %@test:4
+ estimparams.var_exo = zeros(7, 10);
+ estimparams.var_exo(:,5) = [1;2;3;4;5;6;8];
+ estimparams.param_vals = zeros(7, 10);
+ estimparams.param_vals(:,5) = [1;2;3;4;5;6;8];
+ estimparams.param_vals(2:3,5) = 0;
+
+ try
+ b = isbayesian(estimparams);
+ t(1) = false;
+ catch
+ t(1) = true;
+ end
+
+ T = all(t);
+ %@eof:4
+
+end
+
+
+function [b, i] = ispriorinblock(type, data, b, i)
+ if isempty(data)
+ i(type) = false;
+ else
+ if all(data)
+ b(type) = true;
+ elseif any(data)
+ throwAsCaller(MException('Prior:specificationError', 'All estimated parameters must have a prior or none must have a prior (ML).'));
+ end
+ end
+end
diff --git a/matlab/estimation/non_linear_dsge_likelihood.m b/matlab/estimation/non_linear_dsge_likelihood.m
index a7b13c1798..5e7dcd3bca 100644
--- a/matlab/estimation/non_linear_dsge_likelihood.m
+++ b/matlab/estimation/non_linear_dsge_likelihood.m
@@ -1,4 +1,4 @@
-function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr] = non_linear_dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr] = non_linear_dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr] = non_linear_dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Evaluates the posterior kernel of a dsge model using a non linear filter.
%
@@ -8,6 +8,7 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr
% - dataset_info [struct] Matlab's structure describing the dataset
% - options_ [struct] Matlab's structure describing the options
% - M_ [struct] Matlab's structure describing the M_
+% - Prior [dprior] Definition of the prior density.
% - estim_params_ [struct] Matlab's structure describing the estimated_parameters
% - bayestopt_ [struct] Matlab's structure describing the priors
% - BoundsInfo [struct] Matlab's structure specifying the bounds on the paramater values
@@ -29,7 +30,7 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr
% - bayestopt_ [struct] See INPUTS section.
% - dr [struct] decision rule structure described in INPUTS section.
-% Copyright © 2010-2023 Dynare Team
+% Copyright © 2010-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -223,8 +224,7 @@ options_.warning_for_steadystate = 1;
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
-lnprior = priordens(xparam1(:),bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
-fval = (likelihood-lnprior);
+fval = likelihood-Prior.density(xparam1);
if isnan(fval)
fval = Inf;
diff --git a/matlab/estimation/plot_priors.m b/matlab/estimation/plot_priors.m
index 63b811fd8b..f2eea81e98 100644
--- a/matlab/estimation/plot_priors.m
+++ b/matlab/estimation/plot_priors.m
@@ -1,21 +1,18 @@
-function plot_priors(bayestopt_,M_,estim_params_,options_,optional_title)
-% plot_priors(bayestopt_,M_,estim_params_,options_,optional_title)
-% plots prior density
+function plot_priors(Prior, M_, estim_params_, options_, optional_title)
+
+% Plot prior densities
%
% INPUTS
-% o bayestopt_ [structure]
-% o M_ [structure]
-% o estim_params_ [structure]
-% o options_ [structure]
-% o optional_title [string]
-
-% OUTPUTS
-% None
+% - Prior [dprior] object containing prior distribution information.
+% - M_ [struct] description of the model.
+% - estim_params_ [struct] description of the estimated parameters.
+% - options_ [struct] Dynare's options.
+% - optional_title [char] 1×n array, title.
%
-% SPECIAL REQUIREMENTS
-% None
+% OUTPUTS
+% None
-% Copyright © 2004-2023 Dynare Team
+% Copyright © 2004-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -40,7 +37,7 @@ if nargin<5
else
figurename = optional_title;
end
-npar = length(bayestopt_.p1);
+npar = Prior.length();
[nbplt,nr,nc,~,~,nstar] = pltorg(npar);
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
@@ -59,7 +56,7 @@ for plt = 1:nbplt
for index=1:nstar0
names = [];
i = (plt-1)*nstar + index;
- [x,f] = draw_prior_density(i,bayestopt_);
+ [x, f] = Prior.data_for_plotting_density(i);
[nam,texnam] = get_the_name(i,TeX,M_,estim_params_,options_.varobs);
subplot(nr,nc,index)
hh_plt = plot(x,f,'-k','linewidth',2);
@@ -86,4 +83,4 @@ if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,' \n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
-end
\ No newline at end of file
+end
diff --git a/matlab/estimation/posterior_sampler.m b/matlab/estimation/posterior_sampler.m
index 2668ee1a58..b4bc463d3f 100644
--- a/matlab/estimation/posterior_sampler.m
+++ b/matlab/estimation/posterior_sampler.m
@@ -1,4 +1,4 @@
-function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString)
+function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,oo_,dispString)
% function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString)
% Random Walk Metropolis-Hastings algorithm.
%
@@ -57,7 +57,7 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun
vv = sampler_options.invhess;
% Initialization of the sampler
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
- posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_, dispString);
+ posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,oo_, dispString);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
diff --git a/matlab/estimation/posterior_sampler_initialization.m b/matlab/estimation/posterior_sampler_initialization.m
index 354b00916f..9526e60463 100644
--- a/matlab/estimation/posterior_sampler_initialization.m
+++ b/matlab/estimation/posterior_sampler_initialization.m
@@ -1,5 +1,5 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
- posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_, dispString)
+ posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_, dataset_info, options_, M_, Prior, estim_params_, bayestopt_, oo_, dispString)
% Posterior sampler initialization.
%
@@ -207,7 +207,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
end
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2(j,new_estimated_parameters) = candidate(new_estimated_parameters);
- ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,Prior,estim_params_,bayestopt_,mh_bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if isfinite(ilogpo2(j)) % log density is finite
fprintf(fidlog,[' Blck ' int2str(j) ':\n']);
for i=1:length(ix2(1,:))
diff --git a/matlab/estimation/private/initialize_measurement_errors_covariance_matrix.m b/matlab/estimation/private/initialize_measurement_errors_covariance_matrix.m
new file mode 100644
index 0000000000..f2bc931f8b
--- /dev/null
+++ b/matlab/estimation/private/initialize_measurement_errors_covariance_matrix.m
@@ -0,0 +1,30 @@
+function M_ = initialize_measurement_errors_covariance_matrix(estim_params_, M_, options_)
+
+% If no previously set measurement error, initialize the measurement error covariance matrix H.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+ if estim_params_.nvn && isequal(M_.H, 0)
+ n = length(options_.varobs);
+ M_.H = zeros(n, n);
+ M_.Correlation_matrix_ME = eye(n);
+ end
+
+end
diff --git a/matlab/estimation/private/set_corrn_observable_correspondence.m b/matlab/estimation/private/set_corrn_observable_correspondence.m
new file mode 100644
index 0000000000..2b3f7d8426
--- /dev/null
+++ b/matlab/estimation/private/set_corrn_observable_correspondence.m
@@ -0,0 +1,30 @@
+function estim_params_ = set_corrn_observable_correspondence(estim_params_, endo_names, varobs)
+
+% Stores indices, in options_.varobs, of the corresponding observed variables
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ if ~isfield(estim_params_, 'corrn_observable_correspondence')
+ estim_params_.corrn_observable_correspondence = NaN(estim_params_.ncn, 2);
+ for i=1:estim_params_.ncn
+ estim_params_.corrn_observable_correspondence(i,1) = find(strcmp(endo_names{estim_params_.corrn(i,1)}, varobs));
+ estim_params_.corrn_observable_correspondence(i,2) = find(strcmp(endo_names{estim_params_.corrn(i,2)}, varobs));
+ end
+ end
+
+end
diff --git a/matlab/estimation/private/set_jscale.m b/matlab/estimation/private/set_jscale.m
new file mode 100644
index 0000000000..b31c2fd799
--- /dev/null
+++ b/matlab/estimation/private/set_jscale.m
@@ -0,0 +1,64 @@
+function estim_params_ = set_jscale(estim_params_, overwrite)
+
+% Set the estimated parameter specific MH scaling parameter.
+%
+% INPUTS
+% - estim_params_ [struct] characterizing parameters to be estimated
+% - overwrite [logical] scalar, overwrite jscale field if true (default is false)
+%
+% OUTPUTS
+% - estim_params_ [struct] Updated characterization of the parameters to be estimated.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ if nargin<2
+ overwrite = false;
+ end
+
+ if ~isfield(estim_params_, 'jscale') || overwrite
+
+ estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+ nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+ mm = estim_params_.mm; % Total number of estimated parameters
+
+ estim_params_.jscale = ones(mm, 1);
+
+ if estim_params_.nvx
+ estim_params_.jscale(1:nn(1)) = estim_params_.var_exo(:,10);
+ end
+
+ if estim_params_.nvn
+ estim_params_.jscale(nn(1)+1:nn(2)) = estim_params_.var_endo(:,10);
+ end
+
+ if estim_params_.ncx
+ estim_params_.jscale(nn(2)+1:nn(3)) = estim_params_.corrx(:,11);
+ end
+
+ if estim_params_.ncn
+ estim_params_.jscale((nn(3)+1:nn(4))) = estim_params_.corrn(:,11);
+ end
+
+ if estim_params_.np
+ estim_params_.jscale(nn(4)+1:mm) = estim_params_.param_vals(:,10);
+ end
+
+ end
+
+end
diff --git a/matlab/estimation/private/set_number_of_estimated_parameters.m b/matlab/estimation/private/set_number_of_estimated_parameters.m
new file mode 100644
index 0000000000..342c17f444
--- /dev/null
+++ b/matlab/estimation/private/set_number_of_estimated_parameters.m
@@ -0,0 +1,56 @@
+function estim_params_ = set_number_of_estimated_parameters(estim_params_)
+
+% Add fields to estim_params_ structure (number of estimated parameters).
+%
+% INPUTS
+% - estim_params_ [structure] characterizing parameters to be estimated.
+%
+% OUTPUTS
+% - estim_params_ [structure] characterizing parameters to be estimated
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ if isempty(estim_params_)
+ estim_params_.estimation = false;
+ return
+ end
+
+ if ~isfield(estim_params_, 'nvx')
+ nvx = size(estim_params_.var_exo,1);
+ nvn = size(estim_params_.var_endo,1);
+ ncx = size(estim_params_.corrx,1);
+ ncn = size(estim_params_.corrn,1);
+ np = size(estim_params_.param_vals,1);
+ nn = cumsum([nvx; nvn; ncx; ncn; np]);
+ mm = nn(end);
+ estim_params_.nvx = nvx; % exogenous shock variances
+ estim_params_.nvn = nvn; % endogenous variances, i.e. measurement error
+ estim_params_.ncx = ncx; % exogenous shock correlations
+ estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
+ estim_params_.np = np; % other parameters of the model
+ estim_params_.nn = nn; % cumulated sum of parameter # by type
+ estim_params_.mm = mm; % total number of estimated parameters
+ end
+
+ if isfield(estim_params_, 'mm')
+ estim_params_.estimation = estim_params_.mm>0;
+ else
+ estim_params_.estimation = false;
+ end
+
+end
diff --git a/matlab/estimation/private/set_nvn_observable_correspondence.m b/matlab/estimation/private/set_nvn_observable_correspondence.m
new file mode 100644
index 0000000000..5adab22e51
--- /dev/null
+++ b/matlab/estimation/private/set_nvn_observable_correspondence.m
@@ -0,0 +1,31 @@
+function estim_params_ = set_nvn_observable_correspondence(estim_params_, endo_names, varobs)
+
+% Stores indices, in options_.varobs, of the corresponding observed variables
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ estim_params_.nvn_observable_correspondence = NaN(estim_params_.nvn, 1);
+ for i=1:estim_params_.nvn
+ obs = strcmp(endo_names{estim_params_.var_endo(i,1)}, varobs);
+ if ~any(obs)
+ derror('Estimation:wrongMeasurementEquation', sprintf('The variable %s has to be declared as observable since you assume a measurement error on it.', M_.endo_names{estim_params_.var_endo(i,1)}))
+ end
+ estim_params_.nvn_observable_correspondence(i) = find(obs);
+ end
+
+end
diff --git a/matlab/estimation/private/set_optimization_boundaries.m b/matlab/estimation/private/set_optimization_boundaries.m
new file mode 100644
index 0000000000..486a1f636f
--- /dev/null
+++ b/matlab/estimation/private/set_optimization_boundaries.m
@@ -0,0 +1,66 @@
+function [lb, ub] = set_optimization_boundaries(estim_params_, Prior)
+
+% Set optimization boundaries on the estimated parameters
+%
+% INPUTS:
+% - estim_params_ [struct] characterizing parameters to be estimated.
+% - Prior [dprior] Prior distributon (if any)
+%
+% OUTPUTS:
+% - lb [double] mm×1 vector, lower bounds
+% - ub [double] mm×1 vector, lower bounds
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+ mm = estim_params_.mm; % Total number of estimated parameters
+
+ ub = Inf(mm, 1); % Upper bound imposed for optimization.
+ lb = -Inf(mm, 1); % Lower bound imposed for optimization.
+
+ if estim_params_.nvx
+ ub(1:nn(1)) = estim_params_.var_exo(:,4);
+ lb(1:nn(1)) = estim_params_.var_exo(:,3);
+ end
+
+ if estim_params_.nvn
+ ub(nn(1)+1:nn(2)) = estim_params_.var_endo(:,4);
+ lb(nn(1)+1:nn(2)) = estim_params_.var_endo(:,3);
+ end
+
+ if estim_params_.ncx
+ ub(nn(2)+1:nn(3)) = max(min(estim_params_.corrx(:,5),1),-1);
+ lb(nn(2)+1:nn(3)) = min(max(estim_params_.corrx(:,4),-1),1);
+ end
+
+ if estim_params_.ncn
+ ub(nn(3)+1:nn(4)) = max(min(estim_params_.corrn(:,5),1),-1);
+ lb(nn(3)+1:nn(4)) = min(max(estim_params_.corrn(:,4),-1),1);
+ end
+
+ if estim_params_.np
+ ub(nn(4)+1:mm) = estim_params_.param_vals(:,4);
+ lb(nn(4)+1:mm) = estim_params_.param_vals(:,3);
+ end
+
+ if ~isempty(Prior)
+ ub = min(Prior.ub, ub);
+ lb = max(Prior.lb, lb);
+ end
+
+end
diff --git a/matlab/estimation/private/set_parameter_names.m b/matlab/estimation/private/set_parameter_names.m
new file mode 100644
index 0000000000..52a299d81b
--- /dev/null
+++ b/matlab/estimation/private/set_parameter_names.m
@@ -0,0 +1,68 @@
+function estim_params_ = set_parameter_names(estim_params_, paramnames, exonames, endonames, varobs)
+
+% Set names of the estimated parameters under a new field of estim_params_.
+%
+% INPUTS
+% - estim_params_ [struct]
+% - paramnames [cell] parameter names
+% - exonames [cell] names of the exogenous variables
+% - endonames [cell] names of the endogenous variables
+% - varobs [cell] names of the observed variables
+%
+% OUTPUTS
+% - estim_params_ [struct]
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+ if ~isfield('names', estim_params_)
+
+ nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+ mm = estim_params_.mm; % Total number of estimated parameters
+
+ estim_params_.names = cell(mm, 1); % names of the estimated parameters
+
+ if estim_params_.nvx
+ estim_params_.names(1:nn(1)) = exonames(estim_params_.var_exo(:,1));
+ end
+
+ if estim_params_.nvn
+ estim_params_.names(nn(1)+1:nn(2)) = varobs(estim_params_.nvn_observable_correspondence);
+ end
+
+ if estim_params_.ncx
+ for i = 1:estim_params_.ncx
+ estim_params_.names(nn(2)+i) = {sprintf('corr %s, %s', exonames{estim_params_.corrx(i,1)}, exonames{estim_params_.corrx(i,2)})};
+ end
+ end
+
+ if estim_params_.ncn
+ for i=1:estim_params_.ncn
+ id = estim_params_.corrn_observable_correspondence(i,:);
+ estim_params_.names(nn(3)+i) = {sprintf('corr %s, %s', varobs{id(1)}, varobs{id(2)})};
+ end
+ end
+
+ if estim_params_.np
+ estim_params_.names(nn(4)+1:mm) = paramnames(estim_params_.param_vals(:,1));
+ end
+
+ end
+
+end
diff --git a/matlab/estimation/set_boundaries.m b/matlab/estimation/set_boundaries.m
new file mode 100644
index 0000000000..a1a1e3da8c
--- /dev/null
+++ b/matlab/estimation/set_boundaries.m
@@ -0,0 +1,68 @@
+function [lb, ub] = set_boundaries(estim_params_, M_, options_)
+
+
+% Set optimization boundaries on the estimated parameters.
+%
+% INPUTS
+% - estim_params_ [struct] characterizing parameters to be estimated.
+% - M_ [struct] characterizing the model.
+% - options_ [struct] characterizing the options.
+%
+% OUTPUTS
+% - lb [double] n×1 vector, lower bounds.
+% - ub [double] n×1 vector, upper bounds.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ estim_params_ = set_number_of_estimated_parameters(estim_params_);
+
+ nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+ mm = estim_params_.mm; % Total number of estimated parameters
+
+ lb = -Inf(mm, 1);
+ ub = Inf(mm, 1);
+
+ if estim_params_.nvx
+ lb(1:nn(1)) = max(estim_params_.var_exo(:,3), 0);
+ ub(1:nn(1)) = estim_params_.var_exo(:,4);
+ problem = lb(1:nn(1))>ub(1:nn(1));
+ if any(problem)
+ derror('Estimation:wrongBoundaries', sprintf('Check the declared bounds on (%s), the lower bound must be smaller than the upper bound.', regexprep(sprintf('%s, ', p{[1, 3]}), ', $', '')))
+ end
+ end
+
+ if estim_params_.nvn
+ lb(nn(1)+1:nn(2)) = max(estim_params_.var_endo(:,3), 0);
+ ub(nn(1)+1:nn(2)) = estim_params_.var_endo(:,4);
+ end
+
+ if estim_params_.ncx
+ lb(nn(2)+1:nn(3)) = min(max(estim_params_.corrx(:,4), -1), 1);
+ ub(nn(2)+1:nn(3)) = max(min(estim_params_.corrx(:,5), 1), -1);
+ end
+
+ if estim_params_.ncn
+ lb(nn(3)+1:nn(4)) = min(max(estim_params_.corrn(:,4), -1), 1);
+ ub(nn(3)+1:nn(4)) = max(min(estim_params_.corrn(:,5), 1), -1);
+ end
+
+ if estim_params_.np
+ lb(nn(4)+1:mm) = estim_params_.param_vals(:,3);
+ ub(nn(4)+1:mm) = estim_params_.param_vals(:,4);
+ end
+end
diff --git a/matlab/estimation/set_prior.m b/matlab/estimation/set_prior.m
index 238aca5cd3..198c96496c 100644
--- a/matlab/estimation/set_prior.m
+++ b/matlab/estimation/set_prior.m
@@ -1,24 +1,16 @@
-function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_)
-% function [xparam1,estim_params_,bayestopt_,lb,ub, M_]=set_prior(estim_params_, M_, options_)
+function Prior = set_prior(estim_params_, M_, options_)
+
% sets prior distributions
%
% INPUTS
-% o estim_params_ [structure] characterizing parameters to be estimated.
-% o M_ [structure] characterizing the model.
-% o options_ [structure] characterizing the options.
+% - estim_params_ [struct] characterizing parameters to be estimated.
+% - M_ [struct] characterizing the model.
+% - options_ [struct] characterizing the options.
%
% OUTPUTS
-% o xparam1 [double] vector of parameters to be estimated (initial values)
-% o estim_params_ [structure] characterizing parameters to be estimated
-% o bayestopt_ [structure] characterizing priors
-% o lb [double] vector of lower bounds for the estimated parameters.
-% o ub [double] vector of upper bounds for the estimated parameters.
-% o M_ [structure] characterizing the model.
-%
-% SPECIAL REQUIREMENTS
-% None
+% - Prior [dprior]
-% Copyright © 2003-2023 Dynare Team
+% Copyright © 2003-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -35,254 +27,233 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-nvx = size(estim_params_.var_exo,1);
-nvn = size(estim_params_.var_endo,1);
-ncx = size(estim_params_.corrx,1);
-ncn = size(estim_params_.corrn,1);
-np = size(estim_params_.param_vals,1);
+estim_params_ = set_number_of_estimated_parameters(estim_params_);
-estim_params_.nvx = nvx; %exogenous shock variances
-estim_params_.nvn = nvn; %endogenous variances, i.e. measurement error
-estim_params_.ncx = ncx; %exogenous shock correlations
-estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
-estim_params_.np = np; % other parameters of the model
+nn = estim_params_.nn; % Cumulated sum of estimated parameters by type (structural innovations std., measurement errors std., structural innovations correlations, measurement errors correlations, and deep parameters)
+mm = estim_params_.mm; % Total number of estimated parameters
-xparam1 = [];
-ub = []; % Upper bound imposed for optimization.
-lb = []; % Lower bound imposed for optimization.
-bayestopt_.pshape = [];
-bayestopt_.p1 = []; % prior mean
-bayestopt_.p2 = []; % prior standard deviation
-bayestopt_.p3 = []; % lower bound of the distribution, only considering whether a generalized distribution is used, not when the prior is truncated
-bayestopt_.p4 = []; % upper bound of the distribution, only considering whether a generalized distribution is used, not when the prior is truncated
-bayestopt_.p5 = zeros(nvx+nvn+ncx+ncn+np,1); % prior mode
-bayestopt_.p6 = []; % first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
-bayestopt_.p7 = []; % second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
+prior.pshape = NaN(mm, 1); % Index for the prior distribution name:
+ %
+ % 1 → (Generalized) Beta distribution,
+ % 2 → Gamma distribution (possibly shifted)
+ % 3 → Gaussian distribution
+ % 4 → Inverse gamma, type I distribution (possibly shifted)
+ % 5 → Uniform distribution
+ % 6 → Inverse gamma, type II distribution (possibly shifted)
+ % 8 → Weibull distribution (possibly shifted)
+ %
+prior.p1 = NaN(mm, 1); % prior mean
+prior.p2 = NaN(mm, 1); % prior standard deviation
+prior.p3 = NaN(mm, 1); % lower bound of the distribution, only considering whether a generalized distribution is used, not when the prior is truncated
+prior.p4 = NaN(mm, 1); % upper bound of the distribution, only considering whether a generalized distribution is used, not when the prior is truncated
+prior.p5 = NaN(mm, 1); % prior mode
+prior.p6 = NaN(mm, 1); % first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
+prior.p7 = NaN(mm, 1); % second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
+prior.name = cell(mm, 1); % names of the estimated parameters
-bayestopt_.jscale = []; %jscale might subsequently be overwritten by mode_compute=6 or check_posterior_sampler_options
-bayestopt_.name = {};
-if nvx
- xparam1 = estim_params_.var_exo(:,2);
- ub = estim_params_.var_exo(:,4);
- lb = estim_params_.var_exo(:,3);
- bayestopt_.pshape = estim_params_.var_exo(:,5);
- bayestopt_.p1 = estim_params_.var_exo(:,6);
- bayestopt_.p2 = estim_params_.var_exo(:,7);
- bayestopt_.p3 = estim_params_.var_exo(:,8); %take generalized distribution into account
- bayestopt_.p4 = estim_params_.var_exo(:,9); %take generalized distribution into account
- bayestopt_.jscale = estim_params_.var_exo(:,10);
- bayestopt_.name = M_.exo_names(estim_params_.var_exo(:,1));
+if estim_params_.nvx
+ prior.pshape(1:nn(1)) = estim_params_.var_exo(:,5);
+ prior.p1(1:nn(1)) = estim_params_.var_exo(:,6);
+ prior.p2(1:nn(1)) = estim_params_.var_exo(:,7);
+ prior.p3(1:nn(1)) = estim_params_.var_exo(:,8);
+ prior.p4(1:nn(1)) = estim_params_.var_exo(:,9);
+ prior.name(1:nn(1)) = M_.exo_names(estim_params_.var_exo(:,1));
end
-if nvn
- estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable
- if isequal(M_.H,0) %if no previously set measurement error, initialize H
- nvarobs = length(options_.varobs);
- M_.H = zeros(nvarobs,nvarobs);
- M_.Correlation_matrix_ME = eye(nvarobs);
- end
- for i=1:nvn
- obsi_ = strmatch(M_.endo_names{estim_params_.var_endo(i,1)}, options_.varobs, 'exact');
- if isempty(obsi_)
- error(['The variable ' M_.endo_names{estim_params_.var_endo(i,1)} ' has to be declared as observable since you assume a measurement error on it.'])
- end
- estim_params_.nvn_observable_correspondence(i,1)=obsi_;
- end
- xparam1 = [xparam1; estim_params_.var_endo(:,2)];
- ub = [ub; estim_params_.var_endo(:,4)];
- lb = [lb; estim_params_.var_endo(:,3)];
- bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.var_endo(:,5)];
- bayestopt_.p1 = [ bayestopt_.p1; estim_params_.var_endo(:,6)];
- bayestopt_.p2 = [ bayestopt_.p2; estim_params_.var_endo(:,7)];
- bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)]; %take generalized distribution into account
- bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)]; %take generalized distribution into account
- bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
- bayestopt_.name = [ bayestopt_.name; options_.varobs(estim_params_.nvn_observable_correspondence)];
+
+if estim_params_.nvn
+ prior.pshape(nn(1)+1:nn(2)) = estim_params_.var_endo(:,5);
+ prior.p1(nn(1)+1:nn(2)) = estim_params_.var_endo(:,6);
+ prior.p2(nn(1)+1:nn(2)) = estim_params_.var_endo(:,7);
+ prior.p3(nn(1)+1:nn(2)) = estim_params_.var_endo(:,8);
+ prior.p4(nn(1)+1:nn(2)) = estim_params_.var_endo(:,9);
+ prior.name(nn(1)+1:nn(2)) = options_.varobs(estim_params_.nvn_observable_correspondence);
end
-if ncx
- xparam1 = [xparam1; estim_params_.corrx(:,3)];
- ub = [ub; max(min(estim_params_.corrx(:,5),1),-1)];
- lb = [lb; min(max(estim_params_.corrx(:,4),-1),1)];
- bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrx(:,6)];
- bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrx(:,7)];
- bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrx(:,8)];
- bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrx(:,9)]; %take generalized distribution into account
- bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrx(:,10)]; %take generalized distribution into account
- bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrx(:,11)];
- baseid = length(bayestopt_.name);
- bayestopt_.name = [bayestopt_.name; cell(ncx, 1)];
- for i = 1:ncx
- bayestopt_.name(baseid+i) = {sprintf('corr %s, %s', ...
- M_.exo_names{estim_params_.corrx(i,1)}, ...
- M_.exo_names{estim_params_.corrx(i,2)})};
+
+if estim_params_.ncx
+ prior.pshape(nn(2)+1:nn(3)) = estim_params_.corrx(:,6);
+ prior.p1(nn(2)+1:nn(3)) = estim_params_.corrx(:,7);
+ prior.p2(nn(2)+1:nn(3)) = estim_params_.corrx(:,8);
+ prior.p3(nn(2)+1:nn(3)) = estim_params_.corrx(:,9);
+ prior.p4(nn(2)+1:nn(3)) = estim_params_.corrx(:,10);
+ for i = 1:estim_params_.ncx
+ prior.name(nn(2)+i) = {sprintf('corr %s, %s', M_.exo_names{estim_params_.corrx(i,1)}, M_.exo_names{estim_params_.corrx(i,2)})};
end
end
-if ncn
- estim_params_.corrn_observable_correspondence=NaN(ncn,2);
- if isequal(M_.H,0)
- nvarobs = length(options_.varobs);
- M_.H = zeros(nvarobs,nvarobs);
- M_.Correlation_matrix_ME = eye(nvarobs);
- end
- xparam1 = [xparam1; estim_params_.corrn(:,3)];
- ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
- lb = [lb; min(max(estim_params_.corrn(:,4),-1),1)];
- bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrn(:,6)];
- bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrn(:,7)];
- bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrn(:,8)];
- bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrn(:,9)]; %take generalized distribution into account
- bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)]; %take generalized distribution into account
- bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
- baseid = length(bayestopt_.name);
- bayestopt_.name = [bayestopt_.name; cell(ncn, 1)];
- for i=1:ncn
- k1 = estim_params_.corrn(i,1);
- k2 = estim_params_.corrn(i,2);
- bayestopt_.name(baseid+i) = {sprintf('corr %s, %s', M_.endo_names{k1}, M_.endo_names{k2})};
- % find correspondence to varobs to construct H in set_all_paramters
- obsi1 = strmatch(M_.endo_names{k1}, options_.varobs, 'exact');
- obsi2 = strmatch(M_.endo_names{k2}, options_.varobs, 'exact');
- % save correspondence
- estim_params_.corrn_observable_correspondence(i,:)=[obsi1, obsi2];
+
+if estim_params_.ncn
+ prior.pshape(nn(3)+1:nn(4)) = estim_params_.corrn(:,6);
+ prior.p1(nn(3)+1:nn(4)) = estim_params_.corrn(:,7);
+ prior.p2(nn(3)+1:nn(4)) = estim_params_.corrn(:,8);
+ prior.p3(nn(3)+1:nn(4)) = estim_params_.corrn(:,9);
+ prior.p4(nn(3)+1:nn(4)) = estim_params_.corrn(:,10);
+ for i=1:estim_params_.ncn
+ id = estim_params_.corrn_observable_correspondence(i,:);
+ prior.name(nn(3)+i) = {sprintf('corr %s, %s', options_.varobs{id(1)}, options_.varobs{id(2)})};
end
end
-if np
- xparam1 = [xparam1; estim_params_.param_vals(:,2)];
- ub = [ub; estim_params_.param_vals(:,4)];
- lb = [lb; estim_params_.param_vals(:,3)];
- bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.param_vals(:,5)];
- bayestopt_.p1 = [ bayestopt_.p1; estim_params_.param_vals(:,6)];
- bayestopt_.p2 = [ bayestopt_.p2; estim_params_.param_vals(:,7)];
- bayestopt_.p3 = [ bayestopt_.p3; estim_params_.param_vals(:,8)]; %take generalized distribution into account
- bayestopt_.p4 = [ bayestopt_.p4; estim_params_.param_vals(:,9)]; %take generalized distribution into account
- bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.param_vals(:,10)];
- bayestopt_.name = [bayestopt_.name; M_.param_names(estim_params_.param_vals(:,1))];
+
+if estim_params_.np
+ prior.pshape(nn(4)+1:mm) = estim_params_.param_vals(:,5);
+ prior.p1(nn(4)+1:mm) = estim_params_.param_vals(:,6);
+ prior.p2(nn(4)+1:mm) = estim_params_.param_vals(:,7);
+ prior.p3(nn(4)+1:mm) = estim_params_.param_vals(:,8); %take generalized distribution into account
+ prior.p4(nn(4)+1:mm) = estim_params_.param_vals(:,9); %take generalized distribution into account
+ prior.name(nn(4)+1:mm) = M_.param_names(estim_params_.param_vals(:,1));
end
-bayestopt_.p6 = NaN(size(bayestopt_.p1)) ;
-bayestopt_.p7 = bayestopt_.p6 ;
+%
+% Throw an error if a prior standard deviation is negative
+%
-%% check for point priors and disallow them as they do not work with MCMC
-if any(bayestopt_.p2 ==0)
- error('Error in prior for %s: you cannot use a point prior in estimation. Either increase the prior standard deviation',...
- ' or fix the parameter completely.', bayestopt_.name{bayestopt_.p2 ==0})
+negative_std = prior.p2<0;
+
+if any(negative_std)
+ for i=1:mm
+ dprintf('Prior std. on %s is negative.', prior.name{i})
+ end
+ skipline()
+ derror('Prior:wrongDefinition', 'Check your priors. Negative prior standard deviation is not allowed.')
end
-% generalized location parameters by default for beta distribution
-k = find(bayestopt_.pshape == 1);
-k1 = find(isnan(bayestopt_.p3(k)));
-bayestopt_.p3(k(k1)) = zeros(length(k1),1);
-k1 = find(isnan(bayestopt_.p4(k)));
-bayestopt_.p4(k(k1)) = ones(length(k1),1);
+%
+% Throw an error if a prior has zero variance (degenerate priors will not work with MCMC)
+%
+
+zero_std = prior.p2<eps(1);
+
+if any(zero_std)
+ for i=1:mm
+ dprintf('Prior std. on %s is zero.', prior.name{i})
+ end
+ skipline()
+ derror('Prior:wrongDefinition', 'Check your priors. Degenerate priors are not allowed (increase the prior variance or calibrate the parameter).')
+end
+
+%
+% Set hyperparameters
+%
+% TODO: This should be done in the constructor of the dprior class.
+
+% [1] Beta distrubution
+k = find(prior.pshape == 1);
+% Set default boundaries
+k1 = find(isnan(prior.p3(k)));
+prior.p3(k(k1)) = zeros(length(k1),1);
+k1 = find(isnan(prior.p4(k)));
+prior.p4(k(k1)) = ones(length(k1),1);
+% Compute the mode and set hyperparameters
for i=1:length(k)
- [bayestopt_.p6(k(i)), bayestopt_.p7(k(i))] = beta_specification(bayestopt_.p1(k(i)), bayestopt_.p2(k(i))^2, bayestopt_.p3(k(i)), bayestopt_.p4(k(i)), bayestopt_.name{k(i)});
- if bayestopt_.p6(k(i))<1 || bayestopt_.p7(k(i))<1
- fprintf('Prior distribution for parameter %s has unbounded density!\n',bayestopt_.name{k(i)})
+ [prior.p6(k(i)), prior.p7(k(i))] = beta_specification(prior.p1(k(i)), prior.p2(k(i))^2, prior.p3(k(i)), prior.p4(k(i)), prior.name{k(i)});
+ if prior.p6(k(i))<1 || prior.p7(k(i))<1
+ dprintf('Prior distribution for parameter %s has unbounded density.', prior.name{k(i)})
end
- m = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) , bayestopt_.p4(k(i)) ],1);
- if length(m)==1
- bayestopt_.p5(k(i)) = m;
+ m = compute_prior_mode([prior.p6(k(i)), prior.p7(k(i)), prior.p3(k(i)) , prior.p4(k(i))], 1);
+ if isscalar(m)
+ prior.p5(k(i)) = m;
else
- disp(['Prior distribution for parameter ' bayestopt_.name{k(i)} ' has two modes!'])
- bayestopt_.p5(k(i)) = m(1);
+ dprintf('Prior distribution for parameter %s has two modes.', prior.name{k(i)})
+ prior.p5(k(i)) = m(1);
end
end
-% generalized location parameter by default for gamma distribution
-k = find(bayestopt_.pshape == 2);
-k1 = find(isnan(bayestopt_.p3(k)));
-k2 = find(isnan(bayestopt_.p4(k)));
-bayestopt_.p3(k(k1)) = zeros(length(k1),1);
-bayestopt_.p4(k(k2)) = Inf(length(k2),1);
+% [2] Gamma distribution
+k = find(prior.pshape == 2);
+% Set default boundaries
+k1 = find(isnan(prior.p3(k)));
+k2 = find(isnan(prior.p4(k)));
+prior.p3(k(k1)) = zeros(length(k1), 1);
+prior.p4(k(k2)) = Inf(length(k2), 1);
+% Compute the mode and set hyperparameters
for i=1:length(k)
- [bayestopt_.p6(k(i)), bayestopt_.p7(k(i))] = gamma_specification(bayestopt_.p1(k(i)), bayestopt_.p2(k(i))^2, bayestopt_.p3(k(i)), bayestopt_.name{k(i)});
- if bayestopt_.p6(k(i))<1
- fprintf('Prior distribution for parameter %s has unbounded density!\n',bayestopt_.name{k(i)})
+ [prior.p6(k(i)), prior.p7(k(i))] = gamma_specification(prior.p1(k(i)), prior.p2(k(i))^2, prior.p3(k(i)), prior.name{k(i)});
+ if prior.p6(k(i))<1
+ dprintf('Prior distribution for parameter %s has unbounded density.', prior.name{k(i)})
end
- bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 2) ;
+ prior.p5(k(i)) = compute_prior_mode([prior.p6(k(i)), prior.p7(k(i)), prior.p3(k(i))], 2);
end
-% truncation parameters by default for normal distribution
-k = find(bayestopt_.pshape == 3);
-k1 = find(isnan(bayestopt_.p3(k)));
-k2 = find(isnan(bayestopt_.p4(k)));
-bayestopt_.p3(k(k1)) = -Inf*ones(length(k1),1);
-bayestopt_.p4(k(k2)) = Inf*ones(length(k2),1);
-bayestopt_.p6(k) = bayestopt_.p1(k);
-bayestopt_.p7(k) = bayestopt_.p2(k);
-bayestopt_.p5(k) = bayestopt_.p1(k);
+% [3] Gaussian distribution
+k = find(prior.pshape == 3);
+k1 = find(isnan(prior.p3(k)));
+k2 = find(isnan(prior.p4(k)));
+prior.p3(k(k1)) = -Inf*ones(length(k1),1);
+prior.p4(k(k2)) = Inf*ones(length(k2),1);
+prior.p6(k) = prior.p1(k);
+prior.p7(k) = prior.p2(k);
+prior.p5(k) = prior.p1(k);
-% inverse gamma distribution (type 1)
-k = find(bayestopt_.pshape == 4);
-k1 = find(isnan(bayestopt_.p3(k)));
-k2 = find(isnan(bayestopt_.p4(k)));
-bayestopt_.p3(k(k1)) = zeros(length(k1),1);
-bayestopt_.p4(k(k2)) = Inf(length(k2),1);
+% [4] Inverse gamma distribution (type I)
+k = find(prior.pshape == 4);
+% Set default boundaries
+k1 = find(isnan(prior.p3(k)));
+k2 = find(isnan(prior.p4(k)));
+prior.p3(k(k1)) = zeros(length(k1),1);
+prior.p4(k(k2)) = Inf(length(k2),1);
+% Compute the mode and set hyperparameters
for i=1:length(k)
- [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = inverse_gamma_specification(bayestopt_.p1(k(i)), bayestopt_.p2(k(i))^2, bayestopt_.p3(k(i)), 1, false, bayestopt_.name{k(i)});
- bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4);
+ [prior.p6(k(i)), prior.p7(k(i))] = inverse_gamma_specification(prior.p1(k(i)), prior.p2(k(i))^2, prior.p3(k(i)), 1, false, prior.name{k(i)});
+ prior.p5(k(i)) = compute_prior_mode([prior.p6(k(i)), prior.p7(k(i)), prior.p3(k(i))], 4);
end
-% uniform distribution
-k = find(bayestopt_.pshape == 5);
-problem_parameters='';
+% [5] Uniform distribution
+k = find(prior.pshape == 5);
+problem_prior_specification = false(length(k), 1);
for i=1:length(k)
- [bayestopt_.p1(k(i)),bayestopt_.p2(k(i)),bayestopt_.p6(k(i)),bayestopt_.p7(k(i)),error_indicator] = ...
- uniform_specification(bayestopt_.p1(k(i)),bayestopt_.p2(k(i)),bayestopt_.p3(k(i)),bayestopt_.p4(k(i)));
- if error_indicator
- if isempty(problem_parameters)
- problem_parameters=[bayestopt_.name{k(i)}];
- else
- problem_parameters=[problem_parameters ', ' bayestopt_.name{k(i)}];
- end
- end
- bayestopt_.p3(k(i)) = bayestopt_.p6(k(i)) ;
- bayestopt_.p4(k(i)) = bayestopt_.p7(k(i)) ;
- bayestopt_.p5(k(i)) = NaN ;
+ [prior.p1(k(i)), prior.p2(k(i)), prior.p6(k(i)), prior.p7(k(i)), errorflag] = ...
+ uniform_specification(prior.p1(k(i)), prior.p2(k(i)), prior.p3(k(i)), prior.p4(k(i)), prior.name{k(i)});
+ problem_prior_specification(i) = errorflag;
+ prior.p3(k(i)) = prior.p6(k(i)) ;
+ prior.p4(k(i)) = prior.p7(k(i)) ;
+ prior.p5(k(i)) = NaN;
end
-if ~isempty(problem_parameters)
- error(['uniform_specification: You defined lower and upper bounds for parameter ', problem_parameters, '. In this case, you need to leave mean and standard deviation empty.'])
+if any(problem_prior_specification)
+ names = prior.name{k(problem_prior_specification)};
+ if isscalar(k)
+ errmsg = sprintf('uniform_specification: You defined lower and upper bounds for parameter %s. In this case, you need to leave mean and standard deviation empty.', names{1})
+ else
+ errmsg = sprintf('uniform_specification: You defined lower and upper bounds for parameters %s. In this case, you need to leave mean and standard deviation empty.', strrep(strtrim(sprintf('%s ', names{:})), ' ', ', '))
+ end
+ derror('Prior:wrongDefinition', errmsg)
end
-% inverse gamma distribution (type 2)
-k = find(bayestopt_.pshape == 6);
-k1 = find(isnan(bayestopt_.p3(k)));
-k2 = find(isnan(bayestopt_.p4(k)));
-bayestopt_.p3(k(k1)) = zeros(length(k1),1);
-bayestopt_.p4(k(k2)) = Inf(length(k2),1);
+% [6] Inverse gamma distribution (type II)
+k = find(prior.pshape == 6);
+k1 = find(isnan(prior.p3(k)));
+k2 = find(isnan(prior.p4(k)));
+prior.p3(k(k1)) = zeros(length(k1),1);
+prior.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
- [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
- inverse_gamma_specification(bayestopt_.p1(k(i)), bayestopt_.p2(k(i))^2, bayestopt_.p3(k(i)), 2, false, bayestopt_.name{k(i)});
- bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ;
+ [prior.p6(k(i)), prior.p7(k(i))] = ...
+ inverse_gamma_specification(prior.p1(k(i)), prior.p2(k(i))^2, prior.p3(k(i)), 2, false, prior.name{k(i)});
+ prior.p5(k(i)) = compute_prior_mode([ prior.p6(k(i)) , prior.p7(k(i)) , prior.p3(k(i)) ], 6) ;
end
-% Weibull distribution
-k = find(bayestopt_.pshape == 8);
-k1 = find(isnan(bayestopt_.p3(k)));
-k2 = find(isnan(bayestopt_.p4(k)));
-bayestopt_.p3(k(k1)) = zeros(length(k1),1);
-bayestopt_.p4(k(k2)) = Inf(length(k2),1);
+% [8] Weibull distribution
+k = find(prior.pshape == 8);
+k1 = find(isnan(prior.p3(k)));
+k2 = find(isnan(prior.p4(k)));
+prior.p3(k(k1)) = zeros(length(k1),1);
+prior.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
- [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = weibull_specification(bayestopt_.p1(k(i)), bayestopt_.p2(k(i))^2, bayestopt_.p3(k(i)), bayestopt_.name{k(i)});
- bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 8) ;
+ [prior.p6(k(i)),prior.p7(k(i))] = weibull_specification(prior.p1(k(i)), prior.p2(k(i))^2, prior.p3(k(i)), prior.name{k(i)});
+ prior.p5(k(i)) = compute_prior_mode([ prior.p6(k(i)) , prior.p7(k(i)) , prior.p3(k(i)) ], 8) ;
end
-k = find(isnan(xparam1));
-if ~isempty(k)
- xparam1(k) = bayestopt_.p1(k);
-end
+%
+% Return dprior object
+%
-if options_.initialize_estimated_parameters_with_the_prior_mode
- xparam1 = bayestopt_.p5;
- k = find(isnan(xparam1));% Because the uniform density do not have a mode!
- if ~isempty(k)
- xparam1(k) = bayestopt_.p1(k);
- end
-end
+Prior = dprior(prior, options_.prior_trunc, false);
+%
+% Save prior specification on disk
+%
+% TODO: Check if we really need this anymore.
+return
% I create subfolder M_.dname/prior if needed.
-CheckPath('prior',M_.dname);
+CheckPath('prior', M_.dname);
% I save the prior definition if the prior has changed.
priorfile = sprintf('%s/prior/definition.mat', M_.dname);
@@ -306,12 +277,12 @@ if isfile(priorfile)
prior_has_changed = true;
end
else
- prior_has_changed = 1;
+ prior_has_changed = true;
end
if prior_has_changed
- delete([M_.dname '/prior/definition.mat']);
- save([M_.dname '/prior/definition.mat'],'bayestopt_');
+ delete(priorfile);
+ save(priorfile, 'prior');
end
else
- save([M_.dname '/prior/definition.mat'],'bayestopt_');
+ save(priorfile, 'prior');
end
diff --git a/matlab/estimation/uniform_specification.m b/matlab/estimation/uniform_specification.m
index 0689b3028d..038871deec 100644
--- a/matlab/estimation/uniform_specification.m
+++ b/matlab/estimation/uniform_specification.m
@@ -1,23 +1,25 @@
-function [m,s,p6,p7,error_indicator] = uniform_specification(m,s,p3,p4)
+function [m, s, p6, p7, errorflag] = uniform_specification(m, s, p3, p4, name)
+
% Specification of the uniform density function parameters
%
% INPUTS
-% m: mean
-% s: standard deviation
-% p3: lower bound
-% p4: upper bound
-
+% - m [double] scalar, mean of the uniform random variable
+% - s [double] scalar, standard deviation of the uniform random variable
+% - p3 [double] scalar, lower bound of the unidorm random variable domain
+% - p4 [double] scalar, upper bound of the unidorm random variable domain
+% - name [char] 1×n array, name of the parameter
+%
% OUTPUTS
-% m: mean
-% s: standard deviation
-% p6: lower bound
-% p7: upper bound
-% error_indicator: whether inconsistent prior specification was used
+% - m [double] scalar, mean
+% - s [double] scalar, standard deviation
+% - p6 [double] scalar, lower bound
+% - p7 [double] scalar, upper bound
+% - errorflag [logical] scalar, true if inconsistent prior specification was used
%
% SPECIAL REQUIREMENTS
-% none
+% none
-% Copyright © 2004-2017 Dynare Team
+% Copyright © 2004-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -33,17 +35,22 @@ function [m,s,p6,p7,error_indicator] = uniform_specification(m,s,p3,p4)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-error_indicator=0;
+
+errorflag = false;
+
if ~(isnan(p3) || isnan(p4))
p6 = p3;
p7 = p4;
if ~isnan(m) || ~isnan(s)
- error_indicator=1;
+ errorflag = true;
end
- m = (p3+p4)/2;
- s = (p4-p3)/(sqrt(12));
+ % m and s are overwritten if (wrongly) already defined.
+ m = (p3+p4)/2;
+ s = (p4-p3)/sqrt(12);
else
+ if (isnan(m) || isnan(s))
+ error('Cannot compute the hyperparameters of the uniform distribution (missing prior moments for %s)', name)
+ end
p6 = m-s*sqrt(3);
p7 = m+s*sqrt(3);
-
-end
\ No newline at end of file
+end
diff --git a/matlab/isdebug.m b/matlab/isdebug.m
new file mode 100644
index 0000000000..ac597c1f98
--- /dev/null
+++ b/matlab/isdebug.m
@@ -0,0 +1,39 @@
+function b = isdebug()
+
+% Return true iff we are in debug mode.
+%
+% INPUTS
+% None
+%
+% OUTPUTS
+% - b [logical] scalar.
+%
+% REMARKS
+% - First argument is optional.
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+ debug = getenv('DYNARE_DEBUG_MODE');
+
+ if isempty(debug)
+ b = false;
+ else
+ b = isequal(debug, 'true');
+ end
+
+end
diff --git a/matlab/stochastic_solver/select_qz_criterium_value.m b/matlab/stochastic_solver/select_qz_criterium_value.m
index 36ce832da7..bb4d3bb48d 100644
--- a/matlab/stochastic_solver/select_qz_criterium_value.m
+++ b/matlab/stochastic_solver/select_qz_criterium_value.m
@@ -72,4 +72,4 @@ else
end
end
end
-end
\ No newline at end of file
+end
diff --git a/matlab/utilities/estimation/check_mode_file.m b/matlab/utilities/estimation/check_mode_file.m
index 040c74a80e..84704042b3 100644
--- a/matlab/utilities/estimation/check_mode_file.m
+++ b/matlab/utilities/estimation/check_mode_file.m
@@ -1,23 +1,19 @@
-function [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_)
-% function [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_)
-% -------------------------------------------------------------------------
+function [xparam1, hh] = check_mode_file(xparam1, hh, modecompute, modefile, enames)
+
% Check that the provided mode_file is compatible with the current estimation settings.
-% -------------------------------------------------------------------------
+%
% INPUTS
-% o xparam1: [vector] current vector of parameter values at the mode
-% o hh: [matrix] current hessian matrix at the mode
-% o options_: [structure] information about options
-% o bayestopt_: [structure] information about priors
-% -------------------------------------------------------------------------
+% - xparam1 [double] current vector of parameter values at the mode
+% - hh [double] current hessian matrix at the mode
+% - modecompute [logical] information about options
+% - modefile [char] row char array, name of the mode file
+% - enames [cell] list of estimated parameters (row char arrays)
+%
% OUTPUTS
-% o xparam1: [vector] updated vector of parameter values at the mode
-% o hh: [matrix] updated hessian matrix at the mode
-% -------------------------------------------------------------------------
-% This function is called by
-% o dynare_estimation_init.m
-% -------------------------------------------------------------------------
-
-% Copyright © 2023 Dynare Team
+% - xparam1: [vector] updated vector of parameter values at the mode
+% - hh: [matrix] updated hessian matrix at the mode
+
+% Copyright © 2023-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -35,61 +31,83 @@ function [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
number_of_estimated_parameters = length(xparam1);
-mode_file = load(options_.mode_file);
+
+if ~isfile(modefile)
+ derror('Estimation:modeFileMissing', sprintf('Cannot find %s (mode file).', modefile))
+else
+ mode_file = load(modefile);
+end
+
+if ~(issquare(hh) && rows(hh)==number_of_estimated_parameters && length(enames)==number_of_estimated_parameters)
+ derror('Estimation:WrongInputs', 'Wrong call to check_mode_file(). This is a bug.')
+end
+
+if isfield(mode_file, 'parameter_names') && length(mode_file.parameter_names)~=length(mode_file.xparam1)
+ derror('Estimation:modeFileCorrupted', sprintf('The ''mode_file'' (%s) is not valid (number of elements in xparam1 must match the number of declared parameters)', modefile))
+end
+
+if isfield(mode_file, 'hh') && ~issquare(mode_file.hh)
+ derror('Estimation:modeFileCorrupted', sprintf('The ''mode_file'' (%s) is not valid (hh must be a square matrix)', modefile))
+end
+
+if isfield(mode_file, 'hh') && rows(mode_file.hh)~=length(mode_file.xparam1)
+ derror('Estimation:modeFileCorrupted', sprintf('The ''mode_file'' (%s) is not valid (number of elements in xparam1 must match the number of rows in hh)', modefile))
+end
+
if number_of_estimated_parameters>length(mode_file.xparam1)
% More estimated parameters than parameters in the mode_file.
skipline()
- disp(['The ''mode_file'' ' options_.mode_file ' has been generated using another specification of the model or another model!'])
- disp(['Your file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate ' int2str(number_of_estimated_parameters) ' parameters:'])
+ dprintf('The ''mode_file'' %s has been generated using another specification of the model or another model.', modefile)
+ dprintf('Your file contains estimates for %u parameters, while you are attempting to estimate %u parameters:', length(mode_file.xparam1), number_of_estimated_parameters)
md = []; xd = [];
for i=1:number_of_estimated_parameters
- id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
- if isempty(id)
- disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded ''mode_file'' (prior mean or initialized values will be used, if possible).'])
+ id = strcmp(enames{i}, mode_file.parameter_names);
+ if ~any(id)
+ dprintf('--> Estimated parameter %s is not present in the loaded ''mode_file'' (prior mean or initialized values will be used, if possible).', enames{i})
else
xd = [xd; i];
- md = [md; id];
+ md = [md; find(id)];
end
end
for i=1:length(mode_file.xparam1)
- id = strmatch(mode_file.parameter_names{i},bayestopt_.name,'exact');
- if isempty(id)
- disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
+ id = strcmp(mode_file.parameter_names{i}, enames);
+ if ~any(id)
+ dprintf('--> Parameter %s is not estimated according to the current mod file.', mode_file.parameter_names{i})
end
end
- if ~options_.mode_compute
+ if ~modecompute
% The mode is not estimated.
- error('Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
+ derror('Estimation:modeFileIssue', 'Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
else
% The mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean or initialized values.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
- error('Please remove the ''mode_file'' option.')
+ derror('Estimation:modeFileIssue', 'Please remove the ''mode_file'' option.')
end
end
elseif number_of_estimated_parameters<length(mode_file.xparam1)
% Less estimated parameters than parameters in the mode_file.
skipline()
- disp(['The ''mode_file'' ' options_.mode_file ' has been generated using another specification of the model or another model!'])
- disp(['Your file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate only ' int2str(number_of_estimated_parameters) ' parameters:'])
+ dprintf('The ''mode_file'' %s has been generated using another specification of the model or another model.', modefile)
+ dprintf('Your file contains estimates for %u parameters, while you are attempting to estimate only %u parameters.', length(mode_file.xparam1), number_of_estimated_parameters)
md = []; xd = [];
for i=1:number_of_estimated_parameters
- id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
- if isempty(id)
- disp(['--> Estimated parameter ' deblank(bayestopt_.name(i,:)) ' is not present in the loaded ''mode_file'' (prior mean or initialized values will be used, if possible).'])
+ id = strcmp(enames{i}, mode_file.parameter_names);
+ if ~any(id)
+ dprintf('--> Estimated parameter %s is not present in the loaded ''mode_file'' (prior mean or initialized values will be used, if possible).', enames{i})
else
xd = [xd; i];
- md = [md; id];
+ md = [md; find(id)];
end
end
for i=1:length(mode_file.xparam1)
- id = strmatch(mode_file.parameter_names{i},bayestopt_.name,'exact');
- if isempty(id)
- disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
+ id = strcmp(mode_file.parameter_names{i}, enames);
+ if ~any(id)
+ dprintf('--> Parameter %s is not estimated according to the current mod file.', mode_file.parameter_names{i})
end
end
- if ~options_.mode_compute
+ if ~modecompute
% The posterior mode is not estimated. If possible, fix the mode_file.
if isequal(length(xd),number_of_estimated_parameters)
disp('==> Fix mode_file (remove unused parameters).')
@@ -98,7 +116,7 @@ elseif number_of_estimated_parameters<length(mode_file.xparam1)
hh = mode_file.hh(md,md);
end
else
- error('Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
+ derror('Estimation:modeFileIssue', 'Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
end
else
% The mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode_file using the prior mean or initialized values.
@@ -106,46 +124,47 @@ elseif number_of_estimated_parameters<length(mode_file.xparam1)
xparam1(xd) = mode_file.xparam1(md);
else
% None of the estimated parameters are present in the mode_file.
- error('Please remove the ''mode_file'' option.')
+ derror('Estimation:modeFileIssue', 'Please remove the ''mode_file'' option.')
end
end
else
% The number of declared estimated parameters match the number of parameters in the mode file.
% Check that the parameters in the mode file and according to the current mod file are identical.
if ~isfield(mode_file,'parameter_names')
- disp(['The ''mode_file'' ' options_.mode_file ' has been generated using an older version of Dynare. It cannot be verified if it matches the present model. Proceed at your own risk.'])
- mode_file.parameter_names=deblank(bayestopt_.name); %set names
+ dprintf('The ''mode_file'' %s has been generated using an older version of Dynare. It cannot be verified if it matches the present model. Proceed at your own risk.', modefile)
+ mode_file.parameter_names = enames; %set names
end
- if isequal(mode_file.parameter_names, bayestopt_.name)
+ if isequal(mode_file.parameter_names, enames)
xparam1 = mode_file.xparam1;
if isfield(mode_file,'hh')
hh = mode_file.hh;
end
+ elseif isempty(setdiff(mode_file.parameter_names, enames))
+ % Estimated parameters are the same but ordering is different
+ if isdebug()
+ dprintf('The parameters in the provided ''mode file'' (%s) are not in the same order as the declared estimated parameters.', modefile)
+ dprintf('Reorder parameters from %s are reordered...', modefile)
+ end
+ [~, id] = ismember(mode_file.parameter_names, enames);
+ xparam1 = mode_file.xparam1(id);
+ hh = mode_file.hh(id,id);
else
+ % Parameters are not the same
skipline()
- disp(['The ''mode_file'' ' options_.mode_file ' has been generated using another specification of the model or another model!'])
+ dprintf('The ''mode_file'' %s has been generated using another specification of the model or another model.', modefile)
% Check if this is only an ordering issue or if the missing parameters can be initialized with the prior mean.
md = []; xd = [];
for i=1:number_of_estimated_parameters
- id = strmatch(deblank(bayestopt_.name(i,:)), mode_file.parameter_names,'exact');
- if isempty(id)
- disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded ''mode_file''.'])
+ id = strcmp(enames{i}, mode_file.parameter_names);
+ if ~any(id)
+ dprintf('--> Estimated parameter %s is not present in the loaded ''mode_file''.', enames{i})
else
xd = [xd; i];
- md = [md; id];
+ md = [md; find(id)];
end
end
- if ~options_.mode_compute
- % The mode is not estimated
- if isequal(length(xd), number_of_estimated_parameters)
- % This is an ordering issue.
- xparam1 = mode_file.xparam1(md);
- if isfield(mode_file,'hh')
- hh = mode_file.hh(md,md);
- end
- else
- error('Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
- end
+ if ~modecompute
+ derror('Estimation:modeFileIssue', 'Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
else
% The mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode_file using the prior mean or initialized values.
if ~isempty(xd)
@@ -155,9 +174,339 @@ else
end
else
% None of the estimated parameters are present in the mode_file.
- error('Please remove the ''mode_file'' option.')
+ derror('Estimation:modeFileIssue', 'Please remove the ''mode_file'' option.')
end
end
end
end
-skipline()
\ No newline at end of file
+skipline()
+
+return % --*-- Unit tests --*--
+
+%@test:1
+xparam1 = ones(3,1);
+hh = randn(3,3);
+parameter_names = {'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+try
+ [xparam1, hh] = check_mode_file(rand(3, 1), randn(3, 3), false, 'modefile1.mat', {'a1'; 'a2'; 'a3'});
+ t(1) = false;
+catch ME
+ t(1) = true;
+end
+
+if t(1)
+ t(2) = strcmp(ME.identifier, 'Estimation:modeFileCorrupted');
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:1
+
+%@test:2
+xparam1 = ones(3,1);
+hh = randn(3,3);
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+try
+ [xparam1, hh] = check_mode_file(rand(3, 1), randn(3, 3), false, 'modefile2.mat', {'a1'; 'a2'; 'a3'});
+ t(1) = false;
+catch ME
+ t(1) = true;
+end
+
+if t(1)
+ t(2) = strcmp(ME.identifier, 'Estimation:modeFileMissing');
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:2
+
+%@test:3
+xparam1 = ones(3,1);
+hh = randn(2,3);
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+try
+ [xparam1, hh] = check_mode_file(rand(3, 1), randn(3, 3), false, 'modefile1.mat', {'a1'; 'a2'; 'a3'});
+ t(1) = false;
+catch ME
+ t(1) = true;
+end
+
+if t(1)
+ t(2) = strcmp(ME.identifier, 'Estimation:modeFileCorrupted');
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:3
+
+%@test:4
+xparam1 = ones(3,1);
+hh = randn(2,2);
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+try
+ [xparam1, hh] = check_mode_file(rand(3, 1), randn(3, 3), false, 'modefile1.mat', {'a1'; 'a2'; 'a3'});
+ t(1) = false;
+catch ME
+ t(1) = true;
+end
+
+if t(1)
+ t(2) = strcmp(ME.identifier, 'Estimation:modeFileCorrupted');
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:4
+
+
+%@test:5
+xparam1 = ones(3,1);
+hh = randn(3,3);
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = 2*ones(4,1);
+hh = 2*ones(4,4);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, false, 'modefile1.mat', {'a1'; 'a0'; 'a2'; 'a3'});
+ t(1) = false;
+catch ME
+ t(1) = true;
+end
+
+if t(1)
+ t(2) = strcmp(ME.identifier, 'Estimation:modeFileIssue');
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:5
+
+
+%@test:6
+xparam1 = ones(3,1);
+hh = randn(3,3);
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = 2*ones(4,1);
+hh = 2*ones(4,4);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a1'; 'a0'; 'a2'; 'a3'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [1; 2; 1; 1]);
+ t(3) = isequal(hh, 2*ones(4, 4));
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:6
+
+
+%@test:7
+xparam1 = ones(3,1);
+hh = randn(3,3);
+parameter_names = {'a1'; 'a2'; 'b3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = 2*ones(4,1);
+hh = 2*ones(4,4);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a1'; 'a0'; 'a2'; 'a3'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [1; 2; 1; 2]);
+ t(3) = isequal(hh, 2*ones(4, 4));
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:7
+
+
+%@test:8
+xparam1 = [1; 2; 3];
+hh = [1 2 3; 1 2 3; 1 2 3];
+parameter_names = {'a1'; 'a2'; 'b3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = [-1; -2];
+hh = 2*ones(2, 2);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a1'; 'a2'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [1; 2]);
+ t(3) = isequal(hh, 2*ones(2,2));
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:8
+
+
+%@test:9
+xparam1 = [1; 2; 3];
+hh = [1 2 3; 1 2 3; 1 2 3];
+parameter_names = {'a1'; 'a2'; 'b3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = [-1; -2];
+hh = 2*ones(2, 2);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a2'; 'a1'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [2; 1]);
+ t(3) = isequal(hh, 2*ones(2,2));
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:9
+
+
+%@test:10
+xparam1 = [1; 2; 3];
+hh = [1 2 3; 1 2 3; 1 2 3];
+parameter_names = {'a1'; 'a2'; 'b3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = [-1; -2];
+hh = 2*ones(2, 2);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a2'; 'b1'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [2; -2]);
+ t(3) = isequal(hh, 2*ones(2,2));
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:10
+
+
+%@test:11
+xparam1 = [1; 2; 3];
+hh = [1 2 3; 1 2 3; 1 2 3];
+parameter_names = {'a1'; 'a2'; 'b3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = [-1; -2];
+hh = 2*ones(2, 2);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, false, 'modefile1.mat', {'a2'; 'b1'});
+ t(1) = false;
+catch ME
+ t(1) = true;
+end
+
+if t(1)
+ t(2) = strcmp(ME.identifier, 'Estimation:modeFileIssue');
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:11
+
+
+%@test:12
+xparam1 = [1; 2; 3];
+hh = [1 2 3; 1 2 3; 1 2 3];
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = [-1; -2; -3];
+hh = 2*ones(3, 3);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a3'; 'a2'; 'a1'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [3; 2; 1]);
+ t(3) = isequal(hh, [3 2 1; 3 2 1; 3 2 1]);
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:12
+
+%@test:13
+xparam1 = [1; 2; 3];
+hh = [1 2 3; 1 2 3; 1 2 3];
+parameter_names = {'a1'; 'a2'; 'a3'};
+save('modefile1.mat', 'xparam1', 'hh', 'parameter_names')
+clear('xparam1', 'hh', 'parameter_names')
+
+xparam1 = [-1; -2; -3];
+hh = 2*ones(3, 3);
+
+try
+ [xparam1, hh] = check_mode_file(xparam1, hh, true, 'modefile1.mat', {'a3'; 'a2'; 'b1'});
+ t(1) = true;
+catch ME
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isequal(xparam1, [3; 2; -3]);
+ t(3) = isequal(hh, [3 2 2; 3 2 2; 2 2 2]);
+ delete modefile1.mat
+end
+
+T = all(t);
+%@eof:13
diff --git a/matlab/utilities/estimation/check_prior_stderr_corr.m b/matlab/utilities/estimation/check_prior_stderr_corr.m
index c724390cf2..b5a7485cd0 100644
--- a/matlab/utilities/estimation/check_prior_stderr_corr.m
+++ b/matlab/utilities/estimation/check_prior_stderr_corr.m
@@ -27,21 +27,24 @@ function check_prior_stderr_corr(estim_params_, Prior)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
nvx = estim_params_.nvx; % number of stderr parameters for structural shocks
-if nvx && any(bayestopt_.p3(1:nvx)<0)
+if nvx && any(Prior.p3(1:nvx)<0)
warning('Your prior allows for negative standard deviations for structural shocks. Due to working with variances, Dynare will be able to continue, but it is recommended to change your prior.')
end
offset = nvx;
+
nvn = estim_params_.nvn; % number of stderr parameters for measurement errors
-if nvn && any(bayestopt_.p3(1+offset:offset+nvn)<0)
+if nvn && any(Prior.p3(1+offset:offset+nvn)<0)
warning('Your prior allows for negative standard deviations for measurement error. Due to working with variances, Dynare will be able to continue, but it is recommended to change your prior.')
end
offset = nvx+nvn;
+
ncx = estim_params_.ncx; % number of corr parameters for structural shocks
-if ncx && (any(bayestopt_.p3(1+offset:offset+ncx)<-1) || any(bayestopt_.p4(1+offset:offset+ncx)>1))
+if ncx && (any(Prior.p3(1+offset:offset+ncx)<-1) || any(Prior.p4(1+offset:offset+ncx)>1))
warning('Your prior allows for correlations between structural shocks larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
end
offset = nvx+nvn+ncx;
+
ncn = estim_params_.ncn; % number of corr parameters for measurement errors
-if ncn && (any(bayestopt_.p3(1+offset:offset+ncn)<-1) || any(bayestopt_.p4(1+offset:offset+ncn)>1))
+if ncn && (any(Prior.p3(1+offset:offset+ncn)<-1) || any(Prior.p4(1+offset:offset+ncn)>1))
warning('Your prior allows for correlations between measurement errors larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
-end
\ No newline at end of file
+end
diff --git a/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m b/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m
index 1cb92635da..b9b0e2ceac 100644
--- a/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m
+++ b/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m
@@ -16,7 +16,7 @@ function check_varobs_are_endo_and_declared_once(varobs,endo_names)
% o dynare_estimation_init.m
% -------------------------------------------------------------------------
-% Copyright © 2023 Dynare Team
+% Copyright © 2023-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -33,18 +33,20 @@ function check_varobs_are_endo_and_declared_once(varobs,endo_names)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-number_of_observed_variables = length(varobs);
-for i = 1:number_of_observed_variables
- if ~any(strcmp(varobs{i},endo_names))
- error(['VAROBS: unknown variable (' varobs{i} ')!'])
- end
+list_of_missing_endogenous_variables = setdiff(varobs, endo_names);
+
+if ~isempty(list_of_missing_endogenous_variables)
+ throwAsCaller(MException('Estimation:varobsError', sprintf('Observed variables have to be declared as endogenous variables: %s.', strrep(strtrim(sprintf('%s ', list_of_missing_endogenous_variables{:})), ' ', ', '))));
end
+
+
% Check that a variable is not declared as observed more than once.
+number_of_observed_variables = length(varobs);
if length(unique(varobs))<length(varobs)
for i = 1:number_of_observed_variables
- if sum(strcmp(varobs{i},varobs)) > 1
+ if sum(strcmp(varobs{i},varobs)) > 1
error(['VAROBS: a variable cannot be declared as observed more than once (' varobs{i} ')!'])
end
end
-end
\ No newline at end of file
+end
diff --git a/tests/kalman_filter_fast/swk_0.mod b/tests/kalman_filter_fast/swk_0.mod
index e1f1c35419..700cf83124 100644
--- a/tests/kalman_filter_fast/swk_0.mod
+++ b/tests/kalman_filter_fast/swk_0.mod
@@ -133,7 +133,7 @@ options_.riccati_tol = 0;
// actual computation for timing
tic;
for i=1:100;
- fval = dsge_likelihood(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_,[]);
+ fval = dsge_likelihood(xparam1,dataset_,options_,M_,Prior,estim_params_,bayestopt_,oo_,[]);
end;
toc;
--
GitLab