Skip to content
Snippets Groups Projects
Commit cfe687d4 authored by MichelJuillard's avatar MichelJuillard
Browse files

adapted dsge-var code for new handling of penalties

parent 91f3cc5c
Branches
No related tags found
No related merge requests found
...@@ -194,7 +194,7 @@ while fpar<B ...@@ -194,7 +194,7 @@ while fpar<B
end end
if MAX_nirfs_dsgevar if MAX_nirfs_dsgevar
IRUN = IRUN+1; IRUN = IRUN+1;
[fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); [fval,junk1,junk2,cost_flag,SteadyState,junk3,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight)); DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1))); SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1)));
......
function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults) function [fval,grad,hess,exit_flag,SteadyState,trend_coeff,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
% Evaluates the posterior kernel of the bvar-dsge model. % Evaluates the posterior kernel of the bvar-dsge model.
% %
% INPUTS % INPUTS
...@@ -8,6 +8,10 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelih ...@@ -8,6 +8,10 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelih
% OUTPUTS % OUTPUTS
% o fval [double] Value of the posterior kernel at xparam1. % o fval [double] Value of the posterior kernel at xparam1.
% o cost_flag [integer] Zero if the function returns a penalty, one otherwise. % o cost_flag [integer] Zero if the function returns a penalty, one otherwise.
% o SteadyState [double] Steady state vector possibly recomputed
% by call to dynare_results()
% o trend_coeff [double] place holder for trend coefficients,
% currently not supported by dsge_var
% o info [integer] Vector of informations about the penalty. % o info [integer] Vector of informations about the penalty.
% o PHI [double] Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1). % o PHI [double] Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1).
% o SIGMAu [double] Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1). % o SIGMAu [double] Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1).
...@@ -34,17 +38,18 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelih ...@@ -34,17 +38,18 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelih
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global objective_function_penalty_base
persistent dsge_prior_weight_idx persistent dsge_prior_weight_idx
grad=[]; grad=[];
hess=[]; hess=[];
exit_flag = []; exit_flag = [];
info = []; info = 0;
PHI = []; PHI = [];
SIGMAu = []; SIGMAu = [];
iXX = []; iXX = [];
prior = []; prior = [];
SteadyState = [];
trend_coeff = [];
% Initialization of of the index for parameter dsge_prior_weight in Model.params. % Initialization of of the index for parameter dsge_prior_weight in Model.params.
if isempty(dsge_prior_weight_idx) if isempty(dsge_prior_weight_idx)
...@@ -82,19 +87,21 @@ exit_flag = 1; ...@@ -82,19 +87,21 @@ exit_flag = 1;
% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain. % Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain.
if DynareOptions.mode_compute ~= 1 && any(xparam1 < BoundsInfo.lb) if DynareOptions.mode_compute ~= 1 && any(xparam1 < BoundsInfo.lb)
k = find(xparam1 < BoundsInfo.lb); fval = Inf;
fval = objective_function_penalty_base+sum((BoundsInfo.lb(k)-xparam1(k)).^2);
exit_flag = 0; exit_flag = 0;
info = 41; info(1) = 41;
k = find(xparam1 < BoundsInfo.lb);
info(2) = sum((BoundsInfo.lb(k)-xparam1(k)).^2);
return; return;
end end
% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain. % Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain.
if DynareOptions.mode_compute ~= 1 && any(xparam1 > BoundsInfo.ub) if DynareOptions.mode_compute ~= 1 && any(xparam1 > BoundsInfo.ub)
k = find(xparam1 > BoundsInfo.ub); fval = Inf;
fval = objective_function_penalty_base+sum((xparam1(k)-BoundsInfo.ub(k)).^2);
exit_flag = 0; exit_flag = 0;
info = 42; info(1) = 42;
k = find(xparam1 > BoundsInfo.ub);
info(2) = sum((xparam1(k)-BoundsInfo.ub(k)).^2);
return; return;
end end
...@@ -114,12 +121,14 @@ Model.Sigma_e = Q; ...@@ -114,12 +121,14 @@ Model.Sigma_e = Q;
dsge_prior_weight = Model.params(dsge_prior_weight_idx); dsge_prior_weight = Model.params(dsge_prior_weight_idx);
% Is the dsge prior proper? % Is the dsge prior proper?
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/NumberOfObservations; if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/ ...
fval = objective_function_penalty_base+abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables)); NumberOfObservations;
fval = Inf;
exit_flag = 0; exit_flag = 0;
info = 51; info(1) = 51;
info(2)=dsge_prior_weight; info(2) = abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs; % info(2)=dsge_prior_weight;
% info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs;
return return
end end
...@@ -134,13 +143,13 @@ end ...@@ -134,13 +143,13 @@ end
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). % 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 || ... if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 8 || ...
info(1) == 22 || info(1) == 24 || info(1) == 25 || info(1) == 10 info(1) == 22 || info(1) == 24 || info(1) == 25 || info(1) == 10
fval = objective_function_penalty_base+1; fval = Inf;
info = info(1); info(2) = 0.1;
exit_flag = 0; exit_flag = 0;
return return
elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21 elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21
fval = objective_function_penalty_base+info(2); fval = Inf;
info = info(1); info(2) = 0.1;
exit_flag = 0; exit_flag = 0;
return return
end end
...@@ -209,8 +218,9 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model ...@@ -209,8 +218,9 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model
SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0'); SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0');
[SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu); [SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu);
if ~SIGMAu_is_positive_definite if ~SIGMAu_is_positive_definite
fval = objective_function_penalty_base + penalty; fval = Inf;
info = 52; info(1) = 52;
info(2) = penalty;
exit_flag = 0; exit_flag = 0;
return; return;
end end
...@@ -240,15 +250,17 @@ else% Evaluation of the likelihood of the dsge-var model when the dsge prior wei ...@@ -240,15 +250,17 @@ else% Evaluation of the likelihood of the dsge-var model when the dsge prior wei
end end
if isnan(lik) if isnan(lik)
info = 45; info(1) = 45;
fval = objective_function_penalty_base + 100; info(2) = 0.1;
fval = Inf;
exit_flag = 0; exit_flag = 0;
return return
end end
if imag(lik)~=0 if imag(lik)~=0
info = 46; info(1) = 46;
fval = objective_function_penalty_base + 100; info(2) = 0.1;
fval = Inf;
exit_flag = 0; exit_flag = 0;
return return
end end
...@@ -258,20 +270,22 @@ lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo ...@@ -258,20 +270,22 @@ lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo
fval = (lik-lnprior); fval = (lik-lnprior);
if isnan(fval) if isnan(fval)
info = 47; info(1) = 47;
fval = objective_function_penalty_base + 100; info(2) = 0.1;
fval = Inf;
exit_flag = 0; exit_flag = 0;
return return
end end
if imag(fval)~=0 if imag(fval)~=0
info = 48; info(1) = 48;
fval = objective_function_penalty_base + 100; info(2) = 0.1;
fval = Inf;
exit_flag = 0; exit_flag = 0;
return return
end end
if (nargout == 8) if (nargout == 10)
if isinf(dsge_prior_weight) if isinf(dsge_prior_weight)
iXX = iGXX; iXX = iGXX;
else else
...@@ -279,7 +293,7 @@ if (nargout == 8) ...@@ -279,7 +293,7 @@ if (nargout == 8)
end end
end end
if (nargout==9) if (nargout==11)
if isinf(dsge_prior_weight) if isinf(dsge_prior_weight)
iXX = iGXX; iXX = iGXX;
else else
...@@ -292,3 +306,7 @@ if (nargout==9) ...@@ -292,3 +306,7 @@ if (nargout==9)
prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables; prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
prior.iGXX = iGXX; prior.iGXX = iGXX;
end end
if fval == Inf
pause
end
\ No newline at end of file
...@@ -47,7 +47,7 @@ if ~options_.noconstant ...@@ -47,7 +47,7 @@ if ~options_.noconstant
bvar.NumberOfVariables; bvar.NumberOfVariables;
end end
[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults); [fval,grad,hess,cost_flag,SteadyState,trend_coeff,info,Steady,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
% Conditionnal posterior density of the lagged matrices (given Sigma) -> % Conditionnal posterior density of the lagged matrices (given Sigma) ->
% Matric-variate normal distribution. % Matric-variate normal distribution.
......
...@@ -38,9 +38,9 @@ badg=0; ...@@ -38,9 +38,9 @@ badg=0;
for i=1:n for i=1:n
xiold = x(i); xiold = x(i);
x(i) = xiold+h; x(i) = xiold+h;
f1 = penalty_objective_function(x, fcn, penatly, varargin{:}); f1 = penalty_objective_function(x, fcn, penalty, varargin{:});
x(i) = xiold-h; x(i) = xiold-h;
f2 = penalty_objective_function(x, fcn, penatly, varargin{:}); f2 = penalty_objective_function(x, fcn, penalty, varargin{:});
g0 = (f1-f2)/H; g0 = (f1-f2)/H;
if abs(g0)< 1e15 if abs(g0)< 1e15
g(i)=g0; g(i)=g0;
......
function [fval,DLIK,Hess,exit_flag] = objective_function_penalty(x0,fcn,penalty,varargin) function [fval,DLIK,Hess,exit_flag] = objective_function_penalty(x0,fcn,penalty,varargin)
[fval,DLIK,Hess,exit_flag,SteadyState,trend_coeff,info] = fcn(x0,varargin{:}); [fval,DLIK,Hess,exit_flag,SteadyState,trend_coeff,info] = fcn(x0,varargin{:});
if info(1) ~= 0 if info(1) ~= 0
fval = penalty + info(2); fval = penalty + info(2);
end end
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment