From 301a45ef5959ecff80c237536d26cb72e333efdc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?= Date: Wed, 10 May 2023 08:56:35 +0200 Subject: [PATCH] Add members to @dprior class. --- matlab/@dprior/dprior.m | 105 ++++++++++++++++------------ matlab/maximize_prior_density.m | 33 ++++++--- matlab/minus_logged_prior_density.m | 84 +++++++++++----------- matlab/optimize_prior.m | 39 +++++------ 4 files changed, 145 insertions(+), 116 deletions(-) diff --git a/matlab/@dprior/dprior.m b/matlab/@dprior/dprior.m index b237fd62d..781aac885 100644 --- a/matlab/@dprior/dprior.m +++ b/matlab/@dprior/dprior.m @@ -1,10 +1,13 @@ classdef dprior properties - p6 = []; % Prior first hyperparameter. - p7 = []; % Prior second hyperparameter. + p1 = []; % Prior mean. + p2 = []; % Prior stddev. p3 = []; % Lower bound of the prior support. p4 = []; % Upper bound of the prior support. + p5 = []; % Prior mode. + p6 = []; % Prior first hyperparameter. + p7 = []; % Prior second hyperparameter. lb = []; % Truncated prior lower bound. ub = []; % Truncated prior upper bound. iduniform = []; % Index for the uniform priors. @@ -38,31 +41,34 @@ classdef dprior % % REQUIREMENTS % None. - o.p6 = BayesInfo.p6; - o.p7 = BayesInfo.p7; - o.p3 = BayesInfo.p3; - o.p4 = BayesInfo.p4; - bounds = prior_bounds(BayesInfo, PriorTrunc); - o.lb = bounds.lb; - o.ub = bounds.ub; - if nargin>2 && Uniform - prior_shape = repmat(5, length(o.p6), 1); - else - prior_shape = BayesInfo.pshape; - end - o.idbeta = find(prior_shape==1); - if ~isempty(o.idbeta) - o.isbeta = true; - end - o.idgamma = find(prior_shape==2); - if ~isempty(o.idgamma) - o.isgamma = true; - end - o.idgaussian = find(prior_shape==3); - if ~isempty(o.idgaussian) - o.isgaussian = true; - end - o.idinvgamma1 = find(prior_shape==4); + if isfield(BayesInfo, 'p1'), o.p1 = BayesInfo.p1; end + if isfield(BayesInfo, 'p2'), o.p2 = BayesInfo.p2; end + if isfield(BayesInfo, 'p3'), o.p3 = BayesInfo.p3; end + if isfield(BayesInfo, 'p4'), o.p4 = BayesInfo.p4; end + if isfield(BayesInfo, 'p5'), o.p5 = BayesInfo.p5; end + if isfield(BayesInfo, 'p6'), o.p6 = BayesInfo.p6; end + if isfield(BayesInfo, 'p7'), o.p7 = BayesInfo.p7; end + bounds = prior_bounds(BayesInfo, PriorTrunc); + o.lb = bounds.lb; + o.ub = bounds.ub; + if nargin>2 && Uniform + prior_shape = repmat(5, length(o.p6), 1); + else + prior_shape = BayesInfo.pshape; + end + o.idbeta = find(prior_shape==1); + if ~isempty(o.idbeta) + o.isbeta = true; + end + o.idgamma = find(prior_shape==2); + if ~isempty(o.idgamma) + o.isgamma = true; + end + o.idgaussian = find(prior_shape==3); + if ~isempty(o.idgaussian) + o.isgaussian = true; + end + o.idinvgamma1 = find(prior_shape==4); if ~isempty(o.idinvgamma1) o.isinvgamma1 = true; end @@ -80,6 +86,19 @@ classdef dprior end end % dprior (constructor) + function p = subsref(o, S) + switch S(1).type + case '.' + if ismember(S(1).subs, {'p1','p2','p3','p4','p5','p6','p7','lb','ub'}) + r = builtin('subsref', o, S(1)); + else + error('dprior::subsref: unknown method.') + end + otherwise + error('dprior::subsref: case not implemented.') + end + end + function p = draw(o) % Return a random draw from the prior distribution. % @@ -196,7 +215,7 @@ classdef dprior % REMARKS % Second order derivatives holder, d2lpd, has the same rank and shape than dlpd because the priors are % independent (we would have to use a matrix if non orthogonal priors were allowed in Dynare). - % + % % EXAMPLE % % >> Prior = dprior(bayestopt_, options_.prior_trunc); @@ -288,21 +307,21 @@ classdef dprior case 1 lpd = lpd + sum(lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1))); if isinf(lpd), return, end - case 2 - [tmp, dlpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1)); - lpd = lpd + sum(tmp); - if isinf(lpd), return, end + case 2 + [tmp, dlpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1)); + lpd = lpd + sum(tmp); + if isinf(lpd), return, end case 3 [tmp, dlpd(o.idinvgamma1), d2lpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1)); - lpd = lpd + sum(tmp); - if isinf(lpd), return, end + lpd = lpd + sum(tmp); + if isinf(lpd), return, end case 4 [tmp, dlpd(o.idinvgamma1), d2lpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1)); - lpd = lpd + sum(tmp); - if isinf(lpd) - info = o.idinvgamma1(isinf(tmp)); - return - end + lpd = lpd + sum(tmp); + if isinf(lpd) + info = o.idinvgamma1(isinf(tmp)); + return + end end end if o.isinvgamma2 @@ -316,14 +335,14 @@ classdef dprior if isinf(lpd), return, end case 3 [tmp, dlpd(o.idinvgamma2), d2lpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2)); - lpd = lpd + sum(tmp); - if isinf(lpd), return, end + lpd = lpd + sum(tmp); + if isinf(lpd), return, end case 4 [tmp, dlpd(o.idinvgamma2), d2lpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2)); lpd = lpd + sum(tmp); if isinf(lpd) info = o.idinvgamma2(isinf(tmp)); - return + return end end end @@ -345,7 +364,7 @@ classdef dprior lpd = lpd + sum(tmp); if isinf(lpd) info = o.idweibull(isinf(tmp)); - return + return end end end diff --git a/matlab/maximize_prior_density.m b/matlab/maximize_prior_density.m index 4e7543410..06d879f7d 100644 --- a/matlab/maximize_prior_density.m +++ b/matlab/maximize_prior_density.m @@ -1,9 +1,18 @@ -function [xparams,lpd,hessian_mat] = ... - maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,DynareOptions,DynareModel,BayesInfo,EstimatedParams,DynareResults) +function [xparams, lpd, hessian_mat] = ... + maximize_prior_density(iparams, names, DynareOptions, DynareModel, Prior, EstimatedParams, DynareResults) + % Maximizes the logged prior density using Chris Sims' optimization routine. % % INPUTS -% iparams [double] vector of initial parameters. +% - iparams [double] vector of initial parameters. +% - Prior [dprior] vector specifying prior densities shapes. +% - DynareOptions [struct] Options, AKA options_ +% - DynareModel [struct] Model description, AKA M_ +% - EstimatedParams [struct] Info about estimated parameters, AKA estimated_params_ +% - DynareResults [struct] Results, AKA oo_ + +% +% % prior_shape [integer] vector specifying prior densities shapes. % prior_hyperparameter_1 [double] vector, first hyperparameter. % prior_hyperparameter_2 [double] vector, second hyperparameter. @@ -32,10 +41,18 @@ function [xparams,lpd,hessian_mat] = ... % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[xparams, lpd, exitflag, hessian_mat]=dynare_minimize_objective('minus_logged_prior_density', ... - iparams, DynareOptions.mode_compute, DynareOptions, [prior_inf_bound, prior_sup_bound], ... - BayesInfo.name, BayesInfo, [], ... - prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound, ... - DynareOptions,DynareModel,EstimatedParams,DynareResults); +[xparams, lpd, exitflag, hessian_mat] = dynare_minimize_objective('minus_logged_prior_density', ... + iparams, ... + DynareOptions.mode_compute, ... + DynareOptions, ... + [Prior.p3, Prior.p4], ... + BayesInfo.name, ... + [], ... + [], ... + Prior, + DynareOptions, ... + DynareModel, ... + EstimatedParams, ... + DynareResults); lpd = -lpd; diff --git a/matlab/minus_logged_prior_density.m b/matlab/minus_logged_prior_density.m index 19203297a..41e5a5f54 100644 --- a/matlab/minus_logged_prior_density.m +++ b/matlab/minus_logged_prior_density.m @@ -1,19 +1,20 @@ -function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,DynareOptions,DynareModel,EstimatedParams,DynareResults) +function [fval, info, exitflag, ~, ~] = minus_logged_prior_density(xparams, Prior, DynareOptions, DynareModel, EstimatedParams, DynareResults) + % Evaluates minus the logged prior density. % % INPUTS -% xparams [double] vector of parameters. -% pshape [integer] vector specifying prior densities shapes. -% p6 [double] vector, first hyperparameter. -% p7 [double] vector, second hyperparameter. -% p3 [double] vector, prior's lower bound. -% p4 [double] vector, prior's upper bound. +% - xparams [double] vector of parameters. +% - Prior [dprior] vector specifying prior densities shapes. +% - DynareOptions [struct] Options, AKA options_ +% - DynareModel [struct] Model description, AKA M_ +% - EstimatedParams [struct] Info about estimated parameters, AKA estimated_params_ +% - DynareResults [struct] Results, AKA oo_ % % OUTPUTS -% f [double] value of minus the logged prior density. -% info [double] vector: second entry stores penalty, first entry the error code. -% -% Copyright © 2009-2017 Dynare Team +% - fval [double] value of minus the logged prior density. +% - info [double] 4×1 vector, second entry stores penalty, first entry the error code, last entry a penalty (used for optimization). + +% Copyright © 2009-2023 Dynare Team % % This file is part of Dynare. % @@ -30,10 +31,7 @@ function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparam % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -fake_1 = 1; -fake_2 = 1; - -exit_flag = 1; +exitflag = true; info = zeros(4,1); %------------------------------------------------------------------------------ @@ -41,74 +39,75 @@ info = zeros(4,1); %------------------------------------------------------------------------------ % Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain. -if ~isequal(DynareOptions.mode_compute,1) && any(xparamsp4) - k = find(xparams>p4); +if ~isequal(DynareOptions.mode_compute,1) && any(xparams>Prior.p4) + k = find(xparams>Prior.p4); fval = Inf; - exit_flag = 0; + exitflag = false; info(1) = 42; - info(4) = sum((xparams(k)-p4(k)).^2); + info(4) = sum((xparams(k)-Prior.p4(k)).^2); return end % Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H). -DynareModel = set_all_parameters(xparams,EstimatedParams,DynareModel); +DynareModel = set_all_parameters(xparams, EstimatedParams, DynareModel); Q = DynareModel.Sigma_e; H = DynareModel.H; % Test if Q is positive definite. -if ~issquare(Q) || EstimatedParams.ncx || isfield(EstimatedParams,'calibrated_covariances') +if ~issquare(Q) || EstimatedParams.ncx || isfield(EstimatedParams, 'calibrated_covariances') % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite) [Q_is_positive_definite, penalty] = ispd(Q); if ~Q_is_positive_definite - % The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. + % The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the + % eigenvalues of this matrix in order to build the endogenous penalty. fval = Inf; - exit_flag = 0; + exitflag = false; info(1) = 43; info(4) = penalty; return end - if isfield(EstimatedParams,'calibrated_covariances') - correct_flag=check_consistency_covariances(Q); + if isfield(EstimatedParams, 'calibrated_covariances') + correct_flag = check_consistency_covariances(Q); if ~correct_flag penalty = sum(Q(EstimatedParams.calibrated_covariances.position).^2); fval = Inf; - exit_flag = 0; + exitflag = false; info(1) = 71; info(4) = penalty; return end end - end % Test if H is positive definite. -if ~issquare(H) || EstimatedParams.ncn || isfield(EstimatedParams,'calibrated_covariances_ME') +if ~issquare(H) || EstimatedParams.ncn || isfield(EstimatedParams, 'calibrated_covariances_ME') [H_is_positive_definite, penalty] = ispd(H); if ~H_is_positive_definite - % The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. + % The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues + % of this matrix in order to build the endogenous penalty. fval = Inf; - exit_flag = 0; + exitflag = false; info(1) = 44; info(4) = penalty; return end - if isfield(EstimatedParams,'calibrated_covariances_ME') - correct_flag=check_consistency_covariances(H); + if isfield(EstimatedParams, 'calibrated_covariances_ME') + correct_flag = check_consistency_covariances(H); if ~correct_flag penalty = sum(H(EstimatedParams.calibrated_covariances_ME.position).^2); fval = Inf; - exit_flag = 0; + exitflag = false; info(1) = 72; info(4) = penalty; return @@ -121,8 +120,7 @@ end % 2. Check BK and steady state %----------------------------- -M_ = set_all_parameters(xparams,EstimatedParams,DynareModel); -[dr,info,DynareModel,DynareResults] = resol(0,DynareModel,DynareOptions,DynareResults); +[dr, info, DynareModel, DynareResults] = resol(0, DynareModel, DynareOptions, DynareResults); % Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). if info(1) @@ -132,16 +130,14 @@ if info(1) %meaningful second entry of output that can be used fval = Inf; info(4) = info(2); - exit_flag = 0; + exitflag = false; return else fval = Inf; info(4) = 0.1; - exit_flag = 0; + exitflag = false; return end end - - -fval = - priordens(xparams,pshape,p6,p7,p3,p4); \ No newline at end of file +fval = - Prior.density(xparams); diff --git a/matlab/optimize_prior.m b/matlab/optimize_prior.m index cfad0c090..12e89c892 100644 --- a/matlab/optimize_prior.m +++ b/matlab/optimize_prior.m @@ -1,4 +1,4 @@ -function optimize_prior(DynareOptions, ModelInfo, DynareResults, BayesInfo, EstimationInfo) +function optimize_prior(DynareOptions, ModelInfo, DynareResults, BayesInfo, EstimationInfo, pnames) % This routine computes the mode of the prior density using an optimization algorithm. @@ -19,24 +19,25 @@ function optimize_prior(DynareOptions, ModelInfo, DynareResults, BayesInfo, Esti % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +DynareResults.dr = set_state_space(DynareResults.dr, ModelInfo, DynareOptions); + % Initialize to the prior mean -DynareResults.dr = set_state_space(DynareResults.dr,ModelInfo,DynareOptions); -xparam1 = BayesInfo.p1; +xparam1 = Prior.p1; % Pertubation of the initial condition. -look_for_admissible_initial_condition = 1; scale = 1.0; iter = 0; +look_for_admissible_initial_condition = true; scale = 1.0; iter = 0; while look_for_admissible_initial_condition xinit = xparam1+scale*randn(size(xparam1)); - if all(xinit(:)>BayesInfo.p3) && all(xinit(:)Prior.p3) && all(xinit