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

Added check for BK and steady state in the optimization of the prior density.

time-shift
Stéphane Adjemian (Charybdis) 2012-09-21 22:59:03 +02:00
parent af7ac86996
commit da159038e6
3 changed files with 103 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
matlab/get_prior_info.m
View File

@ -130,6 +130,7 @@ end
if info==2% Prior optimization.
% Initialize to the prior mode if possible
oo_.dr=set_state_space(oo_.dr,M_,options_);
k = find(~isnan(bayestopt_.p5));
xparam1(k) = bayestopt_.p5(k);
% Pertubation of the initial condition.
@ -139,7 +140,11 @@ if info==2% Prior optimization.
while look_for_admissible_initial_condition
xinit = xparam1+scale*randn(size(xparam1));
if all(xinit(:)>bayestopt_.p3) && all(xinit(:)<bayestopt_.p4)
look_for_admissible_initial_condition = 0;
M_ = set_all_parameters(xinit,estim_params_,M_);
[dr,INFO,M_,options_,oo_] = resol(0,M_,options_,oo_);
if ~INFO(1)
look_for_admissible_initial_condition = 0;
end
else
if iter == 2000;
scale = scale/1.1;
@ -155,7 +160,7 @@ if info==2% Prior optimization.
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4);
bayestopt_.p4,options_,M_,estim_params_,oo_);
% Display the results.
disp(' ')
disp(' ')

6
matlab/maximize_prior_density.m
View File

@ -1,5 +1,5 @@
function [xparams,lpd,hessian] = ...
maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound)
maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,DynareOptions,DynareModel,EstimatedParams,DynareResults)
% Maximizes the logged prior density using Chris Sims' optimization routine.
%
% INPUTS
@ -41,7 +41,7 @@ gradient_method = 2;
gradient_epsilon = 1e-6;
[lpd,xparams,grad,hessian,itct,fcount,retcodehat] = ...
csminwel1('minus_logged_prior_density',iparams,H0,[],crit,nit,gradient_method, gradient_epsilon, ...
prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound);
csminwel1('minus_logged_prior_density',iparams,H0,[],crit,nit,gradient_method, gradient_epsilon, ...
prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,DynareOptions,DynareModel,EstimatedParams,DynareResults);
lpd = -lpd;

96
matlab/minus_logged_prior_density.m
View File

@ -1,4 +1,4 @@
function [f,fake_1, fake_2, fake_3] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4)
function [fval,fake_1, fake_2, exit_flag ] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,DynareOptions,DynareModel,EstimatedParams,DynareResults)
% Evaluates minus the logged prior density.
%
% INPUTS
@ -28,8 +28,98 @@ function [f,fake_1, fake_2, fake_3] = minus_logged_prior_density(xparams,pshape,
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global objective_function_penalty_base
fake_1 = 1;
fake_2 = 1;
fake_3 = 1;
f = - priordens(xparams,pshape,p6,p7,p3,p4);
exit_flag = 1;
info = 0;
%------------------------------------------------------------------------------
% 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 ~isequal(DynareOptions.mode_compute,1) && any(xparams<p3)
k = find(xparams<p3);
fval = objective_function_penalty_base+sum((p3(k)-xparams(k)).^2);
exit_flag = 0;
info = 41;
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);
fval = objective_function_penalty_base+sum((xparams(k)-p4(k)).^2);
exit_flag = 0;
info = 42;
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);
Q = DynareModel.Sigma_e;
H = DynareModel.H;
% Test if Q is positive definite.
if EstimatedParams.ncx
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
[CholQ,testQ] = chol(Q);
if testQ
% 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.
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = objective_function_penalty_base+sum(-a(k));
exit_flag = 0;
info = 43;
return
end
end
end
% Test if H is positive definite.
if EstimatedParams.ncn
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
[CholH,testH] = chol(H);
if testH
% 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.
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = objective_function_penalty_base+sum(-a(k));
exit_flag = 0;
info = 44;
return
end
end
end
%-----------------------------
% 2. Check BK and steady state
%-----------------------------
M_ = set_all_parameters(xparams,EstimatedParams,DynareModel);
[dr,info,DynareModel,DynareOptions,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) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
== 8 || info(1) == 22 || info(1) == 24 || info(1) == 19
fval = objective_function_penalty_base+1;
info = info(1);
exit_flag = 0;
return
elseif info(1) == 3 || info(1) == 4 || info(1)==6 || info(1) == 20 || info(1) == 21 || info(1) == 23
fval = objective_function_penalty_base+info(2);
info = info(1);
exit_flag = 0;
return
end
fval = - priordens(xparams,pshape,p6,p7,p3,p4);