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

Add members to @dprior class.

remove-particles-submodule
Stéphane Adjemian (Ryûk) 2023-05-10 08:56:35 +02:00
parent 5375070fa3
commit 301a45ef59
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
4 changed files with 145 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
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,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

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,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 <https://www.gnu.org/licenses/>.
[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;

84
matlab/minus_logged_prior_density.m
View File

@ -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 <https://www.gnu.org/licenses/>.
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(xparams<p3)
k = find(xparams<p3);
if ~isequal(DynareOptions.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(DynareOptions.mode_compute,1) && any(xparams>p4)
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);
fval = - Prior.density(xparams);

39
matlab/optimize_prior.m
View File

@ -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 <https://www.gnu.org/licenses/>.
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(:)<BayesInfo.p4)
ModelInfo = set_all_parameters(xinit,EstimationInfo,ModelInfo);
[dr,INFO,ModelInfo,DynareResults] = resol(0,ModelInfo,DynareOptions,DynareResults);
if all(xinit>Prior.p3) && all(xinit<Prior.p4)
ModelInfo = set_all_parameters(xinit, EstimationInfo, ModelInfo);
[dr, INFO, ModelInfo, DynareResults] = resol(0, ModelInfo, DynareOptions, DynareResults);
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
@ -45,23 +46,19 @@ end
% Maximization of the prior density
[xparams, lpd, hessian_mat] = ...
maximize_prior_density(xinit, BayesInfo.pshape, ...
BayesInfo.p6, ...
BayesInfo.p7, ...
BayesInfo.p3, ...
BayesInfo.p4,DynareOptions,ModelInfo,BayesInfo,EstimationInfo,DynareResults);
maximize_prior_density(xinit, pnames, DynareOptions, ModelInfo, Prior, EstimationInfo, DynareResults);
% 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,ModelInfo,EstimationInfo,DynareOptions) '.'])
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
disp([' Prior mode .............. ' num2str(BayesInfo.p5(i)) '.'])
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
dprintf('deep parameter %u: %s.', i, get_the_name(i, 0, ModelInfo, EstimationInfo, DynareOptions))
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()