Merge branch 'factorize_checks' into 'master'

Factorize prior bound and definiteness checks

See merge request Dynare/dynare!1736
time-shift
Sébastien Villemot 2020-06-23 16:17:55 +00:00
commit bcc213644f
5 changed files with 115 additions and 199 deletions

View File

@ -0,0 +1,109 @@
function [fval,info,exit_flag,M_,Q,H]=check_bounds_and_definiteness_estimation(xparam1, M_, options_, estim_params_, bayestopt_)
% function [fval,info,exit_flag]=check_bounds_and_definiteness_estimation(xparam1, M_, options_, estim_params_, bayestopt_)
% Checks whether parameter vector satisfies
%
% INPUTS
% - xparam1 [double] n by 1 vector, estimated parameters.
% - M_ [struct] Matlab's structure describing the Model.
% - options_ [struct] Matlab's structure describing the options.
% - estim_params_ [struct] Matlab's structure describing the estimated_parameters.
% - bayestopt_ [struct] Matlab's structure specifying the bounds on the paramater values (initialized by dynare,aka bayesopt_).
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
% - info [integer] 4 by 1 vector, informations resolution of the model and evaluation of the likelihood.
% - exit_flag [integer] scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood).
% - Q [matrix] Covariance matrix of structural shocks
% - H [matrix] Covariance matrix of measurement errors
% Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
fval = [];
exit_flag = 1;
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<bayestopt_.lb)
k = find(xparam1(:) < bayestopt_.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4) = sum((bayestopt_.lb(k)-xparam1(k)).^2);
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if any(xparam1>bayestopt_.ub)
k = find(xparam1(:)>bayestopt_.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4) = sum((xparam1(k)-bayestopt_.ub(k)).^2);
return
end
M_ = set_all_parameters(xparam1,estim_params_,M_);
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));
if ~Q_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 43;
info(4) = penalty;
return
end
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;
info(1) = 71;
info(4) = penalty;
return
end
end
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));
if ~H_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 44;
info(4) = penalty;
return
end
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;
info(1) = 72;
info(4) = penalty;
return
end
end
end

View File

@ -181,85 +181,11 @@ end
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if isestimation(DynareOptions) && ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BoundsInfo.lb)
k = find(xparam1<BoundsInfo.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4)= sum((BoundsInfo.lb(k)-xparam1(k)).^2);
if analytic_derivation
DLIK=ones(length(xparam1),1);
end
[fval,info,exit_flag,Model,Q,H]=check_bounds_and_definiteness_estimation(xparam1, Model, DynareOptions, EstimatedParameters, BoundsInfo);
if info(1)
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if isestimation(DynareOptions) && ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BoundsInfo.ub)
k = find(xparam1>BoundsInfo.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4)= sum((xparam1(k)-BoundsInfo.ub(k)).^2);
if analytic_derivation
DLIK=ones(length(xparam1),1);
end
return
end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
Q = Model.Sigma_e;
H = Model.H;
% Test if Q is positive definite.
if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Q(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness));
if ~Q_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 43;
info(4) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances')
correct_flag=check_consistency_covariances(Q);
if ~correct_flag
penalty = sum(Q(EstimatedParameters.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 71;
info(4) = penalty;
return
end
end
end
% Test if H is positive definite.
if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
[H_is_positive_definite, penalty] = ispd(H(EstimatedParameters.H_entries_to_check_for_positive_definiteness,EstimatedParameters.H_entries_to_check_for_positive_definiteness));
if ~H_is_positive_definite
fval = Inf;
info(1) = 44;
info(4) = penalty;
exit_flag = 0;
return
end
if isfield(EstimatedParameters,'calibrated_covariances_ME')
correct_flag=check_consistency_covariances(H);
if ~correct_flag
penalty = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 72;
info(4) = penalty;
return
end
end
end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------

View File

@ -101,38 +101,11 @@ mYX = evalin('base', 'mYX');
mXY = evalin('base', 'mXY');
mXX = evalin('base', 'mXX');
% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain.
if isestimation(DynareOptions) && DynareOptions.mode_compute ~= 1 && any(xparam1 < BoundsInfo.lb)
k = find(xparam1 < BoundsInfo.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4)= sum((BoundsInfo.lb(k)-xparam1(k)).^2);
[fval,info,exit_flag,Model,Q]=check_bounds_and_definiteness_estimation(xparam1, Model, DynareOptions, EstimatedParameters, BoundsInfo);
if info(1)
return
end
% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain.
if isestimation(DynareOptions) && DynareOptions.mode_compute ~= 1 && any(xparam1 > BoundsInfo.ub)
k = find(xparam1 > BoundsInfo.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4) = sum((xparam1(k)-BoundsInfo.ub(k)).^2);
return
end
% Get the variance of each structural innovation.
Q = Model.Sigma_e;
for i=1:EstimatedParameters.nvx
k = EstimatedParameters.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i);
end
offset = EstimatedParameters.nvx;
% Update Model.params and Model.Sigma_e.
Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end);
Model.Sigma_e = Q;
% Get the weight of the dsge prior.
dsge_prior_weight = Model.params(dsge_prior_weight_idx);

View File

@ -63,77 +63,11 @@ end
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if isestimation(DynareOptions) && (DynareOptions.mode_compute~=1) && any(xparam1<BoundsInfo.lb)
k = find(xparam1(:) < BoundsInfo.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4) = sum((BoundsInfo.lb(k)-xparam1(k)).^2);
[fval,info,exit_flag,Model,Q,H]=check_bounds_and_definiteness_estimation(xparam1, Model, DynareOptions, EstimatedParameters, BoundsInfo);
if info(1)
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if isestimation(DynareOptions) && (DynareOptions.mode_compute~=1) && any(xparam1>BoundsInfo.ub)
k = find(xparam1(:)>BoundsInfo.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4) = sum((xparam1(k)-BoundsInfo.ub(k)).^2);
return
end
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
Q = Model.Sigma_e;
H = Model.H;
if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Q(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness));
if ~Q_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 43;
info(4) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances')
correct_flag=check_consistency_covariances(Q);
if ~correct_flag
penalty = sum(Q(EstimatedParameters.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 71;
info(4) = penalty;
return
end
end
end
if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
[H_is_positive_definite, penalty] = ispd(H(EstimatedParameters.H_entries_to_check_for_positive_definiteness,EstimatedParameters.H_entries_to_check_for_positive_definiteness));
if ~H_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 44;
info(4) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances_ME')
correct_flag=check_consistency_covariances(H);
if ~correct_flag
penalty = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 72;
info(4) = penalty;
return
end
end
end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------

View File

@ -1,26 +0,0 @@
function b = isestimation(option)
% Returns true if we are currently estimating a model.
% Copyright (C) 2017 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 <http://www.gnu.org/licenses/>.
b = false;
if ischar(option.mode_compute) || option.mode_compute || option.mh_replic
b = true;
end