stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add members to @dprior class.

remove-priordens
Stéphane Adjemian (Ryûk) 2023-05-10 08:56:35 +02:00
parent 394c8e21b8
commit b448780c34
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
4 changed files with 143 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
matlab/@dprior/dprior.m
View File

@ -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,10 +41,13 @@ classdef dprior
%
% REQUIREMENTS
% None.
o.p6 = bayestopt_.p6;
o.p7 = bayestopt_.p7;
o.p3 = bayestopt_.p3;
o.p4 = bayestopt_.p4;
if isfield(bayestopt_, 'p1'), o.p1 = bayestopt_.p1; end
if isfield(bayestopt_, 'p2'), o.p2 = bayestopt_.p2; end
if isfield(bayestopt_, 'p3'), o.p3 = bayestopt_.p3; end
if isfield(bayestopt_, 'p4'), o.p4 = bayestopt_.p4; end
if isfield(bayestopt_, 'p5'), o.p5 = bayestopt_.p5; end
if isfield(bayestopt_, 'p6'), o.p6 = bayestopt_.p6; end
if isfield(bayestopt_, 'p7'), o.p7 = bayestopt_.p7; end
bounds = prior_bounds(bayestopt_, PriorTrunc);
o.lb = bounds.lb;
o.ub = bounds.ub;
@ -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
@ -623,18 +642,18 @@ end % classdef --*-- Unit tests --*--
%$ end
%$ end
%$
%$ BayesInfo.pshape = p0;
%$ BayesInfo.p1 = p1;
%$ BayesInfo.p2 = p2;
%$ BayesInfo.p3 = p3;
%$ BayesInfo.p4 = p4;
%$ BayesInfo.p5 = p5;
%$ BayesInfo.p6 = p6;
%$ BayesInfo.p7 = p7;
%$ bayestopt_.pshape = p0;
%$ bayestopt_.p1 = p1;
%$ bayestopt_.p2 = p2;
%$ bayestopt_.p3 = p3;
%$ bayestopt_.p4 = p4;
%$ bayestopt_.p5 = p5;
%$ bayestopt_.p6 = p6;
%$ bayestopt_.p7 = p7;
%$
%$ % Call the tested routine
%$ try
%$ Prior = dprior(BayesInfo, prior_trunc, false);
%$ Prior = dprior(bayestopt_, prior_trunc, false);
%$
%$ % Compute density at the prior mode
%$ lpdstar = Prior.density(p5);
@ -718,18 +737,18 @@ end % classdef --*-- Unit tests --*--
%$ end
%$ end
%$
%$ BayesInfo.pshape = p0;
%$ BayesInfo.p1 = p1;
%$ BayesInfo.p2 = p2;
%$ BayesInfo.p3 = p3;
%$ BayesInfo.p4 = p4;
%$ BayesInfo.p5 = p5;
%$ BayesInfo.p6 = p6;
%$ BayesInfo.p7 = p7;
%$ bayestopt_.pshape = p0;
%$ bayestopt_.p1 = p1;
%$ bayestopt_.p2 = p2;
%$ bayestopt_.p3 = p3;
%$ bayestopt_.p4 = p4;
%$ bayestopt_.p5 = p5;
%$ bayestopt_.p6 = p6;
%$ bayestopt_.p7 = p7;
%$
%$ % Call the tested routine
%$ try
%$ Prior = dprior(BayesInfo, prior_trunc, false);
%$ Prior = dprior(bayestopt_, prior_trunc, false);
%$ mu = NaN(14,1);
%$
%$ for i=1:14

33
matlab/maximize_prior_density.m
View File

@ -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,options_,M_,bayestopt_,estim_params_,oo_)
function [xparams, lpd, hessian_mat] = ...
maximize_prior_density(iparams, names, options_, M_, Prior, estim_params_, oo_)
% 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.
@ -37,10 +46,18 @@ function [xparams,lpd,hessian_mat] = ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[xparams, lpd, exitflag, hessian_mat]=dynare_minimize_objective('minus_logged_prior_density', ...
iparams, options_.mode_compute, options_, [prior_inf_bound, prior_sup_bound], ...
bayestopt_.name, bayestopt_, [], ...
prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound, ...
options_,M_,estim_params_,oo_);
[xparams, lpd, exitflag, hessian_mat] = dynare_minimize_objective('minus_logged_prior_density', ...
iparams, ...
options_.mode_compute, ...
options_, ...
[Prior.p3, Prior.p4], ...
names, ...
[], ...
[], ...
Prior,
options_, ...
M_, ...
estim_params_, ...
oo_);
lpd = -lpd;

87
matlab/minus_logged_prior_density.m
View File

@ -1,24 +1,19 @@
function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,options_,M_,estim_params_,oo_)
% [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,options_,M_,estim_params_,oo_)
function [fval, info, exitflag, ~, ~] = minus_logged_prior_density(xparams, Prior, options_, M_, estim_params_, oo_)
% 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.
% prior_sup_bound [double] vector, prior's upper bound.
% options_ [structure] describing the options
% M_ [structure] describing the model
% estim_params_ [structure] characterizing parameters to be estimated
% oo_ [structure] storing the results
% - 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.
%
% - 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.
@ -36,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 <https://www.gnu.org/licenses/>.
fake_1 = [];
fake_2 = [];
exit_flag = 1;
exitflag = true;
info = zeros(4,1);
%------------------------------------------------------------------------------
@ -47,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(options_.mode_compute,1) && any(xparams<p3)
k = find(xparams<p3);
if ~isequal(options_.mode_compute, 1) && any(xparams<Prior.p3)
k = find(xparams<Prior.p3);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 41;
info(4) = sum((p3(k)-xparams(k)).^2);
info(4) = sum((Prior.p3(k)-xparams(k)).^2);
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if ~isequal(options_.mode_compute,1) && any(xparams>p4)
k = find(xparams>p4);
if ~isequal(options_.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).
M_ = set_all_parameters(xparams,estim_params_,M_);
M_ = set_all_parameters(xparams, estim_params_, M_);
Q = M_.Sigma_e;
H = M_.H;
% Test if Q is positive definite.
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_,'calibrated_covariances')
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_, '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(estim_params_,'calibrated_covariances')
correct_flag=check_consistency_covariances(Q);
if isfield(estim_params_, 'calibrated_covariances')
correct_flag = check_consistency_covariances(Q);
if ~correct_flag
penalty = sum(Q(estim_params_.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 71;
info(4) = penalty;
return4
end
end
end
% Test if H is positive definite.
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_,'calibrated_covariances_ME')
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_, '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(estim_params_,'calibrated_covariances_ME')
correct_flag=check_consistency_covariances(H);
if isfield(estim_params_, 'calibrated_covariances_ME')
correct_flag = check_consistency_covariances(H);
if ~correct_flag
penalty = sum(H(estim_params_.calibrated_covariances_ME.position).^2);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 72;
info(4) = penalty;
return
@ -127,7 +120,7 @@ end
% 2. Check BK and steady state
%-----------------------------
[~,info] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
[~, info] = resol(0, M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
@ -137,14 +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);
fval = - Prior.density(xparams);

41
matlab/optimize_prior.m
View File

@ -1,5 +1,5 @@
function optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
% optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
function optimize_prior(options_, M_, oo_, Prior, estim_params_, pnames)
% This routine computes the mode of the prior density using an optimization algorithm.
%
% INPUTS
@ -26,24 +26,25 @@ function optimize_prior(options_, M_, oo_, bayestopt_, 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/>.
oo_.dr = set_state_space(oo_.dr, M_, options_);
% Initialize to the prior mean
oo_.dr = set_state_space(oo_.dr,M_);
xparam1 = bayestopt_.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(:)>bayestopt_.p3) && all(xinit(:)<bayestopt_.p4)
M_ = set_all_parameters(xinit,estim_params_,M_);
[oo_.dr,INFO,M_.params] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if all(xinit>Prior.p3) && all(xinit<Prior.p4)
M_ = set_all_parameters(xinit, estim_params_, M_);
[dr, INFO, M_, oo_] = resol(0, M_, options_, oo_);
if ~INFO(1)
look_for_admissible_initial_condition = 0;
look_for_admissible_initial_condition = false;
end
else
if iter == 2000
if iter==2000
scale = scale/1.1;
iter = 0;
iter = 0;
else
iter = iter+1;
end
@ -52,23 +53,19 @@ end
% Maximization of the prior density
[xparams, lpd, hessian_mat] = ...
maximize_prior_density(xinit, bayestopt_.pshape, ...
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4,options_,M_,bayestopt_,estim_params_,oo_);
maximize_prior_density(xinit, pnames, options_, M_, Prior, estim_params_, oo_);
% Display the results.
% Display results.
skipline(2)
disp('------------------')
disp('PRIOR OPTIMIZATION')
disp('------------------')
skipline()
for i = 1:length(xparams)
disp(['deep parameter ' int2str(i) ': ' get_the_name(i,0,M_,estim_params_,options_) '.'])
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
disp([' Prior mode .............. ' num2str(bayestopt_.p5(i)) '.'])
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
dprintf('deep parameter %u: %s.', i, get_the_name(i, 0, M_, estim_params_, options_))
dprintf(' Initial condition ........ %s.', num2str(xinit(i)))
dprintf(' Prior mode ............... %s.', num2str(Prior.p5(i)))
dprintf(' Optimized prior mode ..... %s.', num2str(xparams(i)))
skipline()
end
skipline()
skipline()