Evaluates the posterior kernel of a dsge model using a non linear filter.
0001 function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = non_linear_dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) 0002 % Evaluates the posterior kernel of a dsge model using a non linear filter. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults}] =} non_linear_dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults}) 0006 %! @anchor{dsge_likelihood} 0007 %! @sp 1 0008 %! Evaluates the posterior kernel of a dsge model using a non linear filter. 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item xparam1 0014 %! Vector of doubles, current values for the estimated parameters. 0015 %! @item DynareDataset 0016 %! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}). 0017 %! @item DynareOptions 0018 %! Matlab's structure describing the options (initialized by dynare, see @ref{options_}). 0019 %! @item Model 0020 %! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}). 0021 %! @item EstimatedParamemeters 0022 %! Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}). 0023 %! @item BayesInfo 0024 %! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}). 0025 %! @item DynareResults 0026 %! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}). 0027 %! @end table 0028 %! @sp 2 0029 %! @strong{Outputs} 0030 %! @sp 1 0031 %! @table @ @var 0032 %! @item fval 0033 %! Double scalar, value of (minus) the likelihood. 0034 %! @item exit_flag 0035 %! Integer scalar, equal to zero if the routine return with a penalty (one otherwise). 0036 %! @item ys 0037 %! Vector of doubles, steady state level for the endogenous variables. 0038 %! @item trend_coeffs 0039 %! Matrix of doubles, coefficients of the deterministic trend in the measurement equation. 0040 %! @item info 0041 %! Integer scalar, error code. 0042 %! @table @ @code 0043 %! @item info==0 0044 %! No error. 0045 %! @item info==1 0046 %! The model doesn't determine the current variables uniquely. 0047 %! @item info==2 0048 %! MJDGGES returned an error code. 0049 %! @item info==3 0050 %! Blanchard & Kahn conditions are not satisfied: no stable equilibrium. 0051 %! @item info==4 0052 %! Blanchard & Kahn conditions are not satisfied: indeterminacy. 0053 %! @item info==5 0054 %! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure. 0055 %! @item info==6 0056 %! The jacobian evaluated at the deterministic steady state is complex. 0057 %! @item info==19 0058 %! The steadystate routine thrown an exception (inconsistent deep parameters). 0059 %! @item info==20 0060 %! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations). 0061 %! @item info==21 0062 %! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state. 0063 %! @item info==22 0064 %! The steady has NaNs. 0065 %! @item info==23 0066 %! M_.params has been updated in the steadystate routine and has complex valued scalars. 0067 %! @item info==24 0068 %! M_.params has been updated in the steadystate routine and has some NaNs. 0069 %! @item info==30 0070 %! Ergodic variance can't be computed. 0071 %! @item info==41 0072 %! At least one parameter is violating a lower bound condition. 0073 %! @item info==42 0074 %! At least one parameter is violating an upper bound condition. 0075 %! @item info==43 0076 %! The covariance matrix of the structural innovations is not positive definite. 0077 %! @item info==44 0078 %! The covariance matrix of the measurement errors is not positive definite. 0079 %! @item info==45 0080 %! Likelihood is not a number (NaN). 0081 %! @item info==45 0082 %! Likelihood is a complex valued number. 0083 %! @end table 0084 %! @item Model 0085 %! Matlab's structure describing the model (initialized by dynare, see @ref{M_}). 0086 %! @item DynareOptions 0087 %! Matlab's structure describing the options (initialized by dynare, see @ref{options_}). 0088 %! @item BayesInfo 0089 %! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}). 0090 %! @item DynareResults 0091 %! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}). 0092 %! @end table 0093 %! @sp 2 0094 %! @strong{This function is called by:} 0095 %! @sp 1 0096 %! @ref{dynare_estimation_1}, @ref{mode_check} 0097 %! @sp 2 0098 %! @strong{This function calls:} 0099 %! @sp 1 0100 %! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{priordens} 0101 %! @end deftypefn 0102 %@eod: 0103 0104 % Copyright (C) 2010, 2011, 2012 Dynare Team 0105 % 0106 % This file is part of Dynare. 0107 % 0108 % Dynare is free software: you can redistribute it and/or modify 0109 % it under the terms of the GNU General Public License as published by 0110 % the Free Software Foundation, either version 3 of the License, or 0111 % (at your option) any later version. 0112 % 0113 % Dynare is distributed in the hope that it will be useful, 0114 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0115 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0116 % GNU General Public License for more details. 0117 % 0118 % You should have received a copy of the GNU General Public License 0119 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0120 0121 % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr 0122 % frederic DOT karame AT univ DASH lemans DOT fr 0123 0124 % Declaration of the penalty as a persistent variable. 0125 persistent penalty 0126 persistent init_flag 0127 persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1 0128 persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations 0129 0130 % Initialization of the persistent variable. 0131 if ~nargin || isempty(penalty) 0132 penalty = 1e8; 0133 if ~nargin, return, end 0134 end 0135 if nargin==1 0136 penalty = xparam1; 0137 return 0138 end 0139 0140 % Initialization of the returned arguments. 0141 fval = []; 0142 ys = []; 0143 trend_coeff = []; 0144 exit_flag = 1; 0145 0146 % Set the number of observed variables 0147 nvobs = DynareDataset.info.nvobs; 0148 0149 %------------------------------------------------------------------------------ 0150 % 1. Get the structural parameters & define penalties 0151 %------------------------------------------------------------------------------ 0152 0153 % Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain. 0154 if (DynareOptions.mode_compute~=1) & any(xparam1<BayesInfo.lb) 0155 k = find(xparam1 < BayesInfo.lb); 0156 fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2); 0157 exit_flag = 0; 0158 info = 41; 0159 return 0160 end 0161 0162 % Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain. 0163 if (DynareOptions.mode_compute~=1) & any(xparam1>BayesInfo.ub) 0164 k = find(xparam1>BayesInfo.ub); 0165 fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2); 0166 exit_flag = 0; 0167 info = 42; 0168 return 0169 end 0170 0171 % Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H). 0172 Q = Model.Sigma_e; 0173 H = Model.H; 0174 for i=1:EstimatedParameters.nvx 0175 k =EstimatedParameters.var_exo(i,1); 0176 Q(k,k) = xparam1(i)*xparam1(i); 0177 end 0178 offset = EstimatedParameters.nvx; 0179 if EstimatedParameters.nvn 0180 for i=1:EstimatedParameters.nvn 0181 k = EstimatedParameters.var_endo(i,1); 0182 H(k,k) = xparam1(i+offset)*xparam1(i+offset); 0183 end 0184 offset = offset+EstimatedParameters.nvn; 0185 else 0186 H = zeros(nvobs); 0187 end 0188 0189 % Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite. 0190 if EstimatedParameters.ncx 0191 for i=1:EstimatedParameters.ncx 0192 k1 =EstimatedParameters.corrx(i,1); 0193 k2 =EstimatedParameters.corrx(i,2); 0194 Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2)); 0195 Q(k2,k1) = Q(k1,k2); 0196 end 0197 % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite) 0198 [CholQ,testQ] = chol(Q); 0199 if testQ 0200 % 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. 0201 a = diag(eig(Q)); 0202 k = find(a < 0); 0203 if k > 0 0204 fval = penalty+sum(-a(k)); 0205 exit_flag = 0; 0206 info = 43; 0207 return 0208 end 0209 end 0210 offset = offset+EstimatedParameters.ncx; 0211 end 0212 0213 % Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite. 0214 if EstimatedParameters.ncn 0215 for i=1:EstimatedParameters.ncn 0216 k1 = DynareOptions.lgyidx2varobs(EstimatedParameters.corrn(i,1)); 0217 k2 = DynareOptions.lgyidx2varobs(EstimatedParameters.corrn(i,2)); 0218 H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2)); 0219 H(k2,k1) = H(k1,k2); 0220 end 0221 % Try to compute the cholesky decomposition of H (possible iff H is positive definite) 0222 [CholH,testH] = chol(H); 0223 if testH 0224 % 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. 0225 a = diag(eig(H)); 0226 k = find(a < 0); 0227 if k > 0 0228 fval = penalty+sum(-a(k)); 0229 exit_flag = 0; 0230 info = 44; 0231 return 0232 end 0233 end 0234 offset = offset+EstimatedParameters.ncn; 0235 end 0236 0237 % Update estimated structural parameters in Mode.params. 0238 if EstimatedParameters.np > 0 0239 Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end); 0240 end 0241 0242 % Update Model.Sigma_e and Model.H. 0243 Model.Sigma_e = Q; 0244 Model.H = H; 0245 0246 %------------------------------------------------------------------------------ 0247 % 2. call model setup & reduction program 0248 %------------------------------------------------------------------------------ 0249 0250 % Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R). 0251 [T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict'); 0252 0253 if info(1) == 1 || info(1) == 2 || info(1) == 5 0254 fval = penalty+1; 0255 exit_flag = 0; 0256 return 0257 elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21 0258 fval = penalty+info(2); 0259 exit_flag = 0; 0260 return 0261 end 0262 0263 % Define a vector of indices for the observed variables. Is this really usefull?... 0264 BayesInfo.mf = BayesInfo.mf1; 0265 0266 % Define the deterministic linear trend of the measurement equation. 0267 if DynareOptions.noconstant 0268 constant = zeros(nvobs,1); 0269 else 0270 if DynareOptions.loglinear 0271 constant = log(SteadyState(BayesInfo.mfys)); 0272 else 0273 constant = SteadyState(BayesInfo.mfys); 0274 end 0275 end 0276 0277 % Define the deterministic linear trend of the measurement equation. 0278 if BayesInfo.with_trend 0279 trend_coeff = zeros(DynareDataset.info.nvobs,1); 0280 t = DynareOptions.trend_coeffs; 0281 for i=1:length(t) 0282 if ~isempty(t{i}) 0283 trend_coeff(i) = evalin('base',t{i}); 0284 end 0285 end 0286 trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs]; 0287 else 0288 trend = repmat(constant,1,DynareDataset.info.ntobs); 0289 end 0290 0291 % Get needed informations for kalman filter routines. 0292 start = DynareOptions.presample+1; 0293 np = size(T,1); 0294 mf = BayesInfo.mf; 0295 Y = transpose(DynareDataset.rawdata); 0296 0297 %------------------------------------------------------------------------------ 0298 % 3. Initial condition of the Kalman filter 0299 %------------------------------------------------------------------------------ 0300 0301 % Get decision rules and transition equations. 0302 dr = DynareResults.dr; 0303 0304 % Set persistent variables (first call). 0305 if isempty(init_flag) 0306 mf0 = BayesInfo.mf0; 0307 mf1 = BayesInfo.mf1; 0308 restrict_variables_idx = BayesInfo.restrict_var_list; 0309 observed_variables_idx = restrict_variables_idx(mf1); 0310 state_variables_idx = restrict_variables_idx(mf0); 0311 sample_size = size(Y,2); 0312 number_of_state_variables = length(mf0); 0313 number_of_observed_variables = length(mf1); 0314 number_of_structural_innovations = length(Q); 0315 init_flag = 1; 0316 end 0317 0318 ReducedForm.ghx = dr.ghx(restrict_variables_idx,:); 0319 ReducedForm.ghu = dr.ghu(restrict_variables_idx,:); 0320 ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:); 0321 ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); 0322 ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); 0323 ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx)); 0324 ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx); 0325 ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx)); 0326 ReducedForm.Q = Q; 0327 ReducedForm.H = H; 0328 ReducedForm.mf0 = mf0; 0329 ReducedForm.mf1 = mf1; 0330 0331 % Set initial condition. 0332 switch DynareOptions.particle.initialization 0333 case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model. 0334 StateVectorMean = ReducedForm.constant(mf0); 0335 StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',1e-12,1e-12); 0336 case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model). 0337 StateVectorMean = ReducedForm.constant(mf0); 0338 old_DynareOptionsperiods = DynareOptions.periods; 0339 DynareOptions.periods = 5000; 0340 y_ = simult(oo_.steady_state, dr); 0341 y_ = y_(state_variables_idx,2001:5000); 0342 StateVectorVariance = cov(y_'); 0343 DynareOptions.periods = old_DynareOptionsperiods; 0344 clear('old_DynareOptionsperiods','y_'); 0345 case 3% Initial state vector covariance is a diagonal matrix. 0346 StateVectorMean = ReducedForm.constant(mf0); 0347 StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables); 0348 otherwise 0349 error('Unknown initialization option!') 0350 end 0351 ReducedForm.StateVectorMean = StateVectorMean; 0352 ReducedForm.StateVectorVariance = StateVectorVariance; 0353 0354 %------------------------------------------------------------------------------ 0355 % 4. Likelihood evaluation 0356 %------------------------------------------------------------------------------ 0357 DynareOptions.warning_for_steadystate = 0; 0358 LIK = feval(DynareOptions.particle.algorithm,ReducedForm,Y,[],DynareOptions); 0359 if imag(LIK) 0360 likelihood = penalty; 0361 exit_flag = 0; 0362 elseif isnan(LIK) 0363 likelihood = penalty; 0364 exit_flag = 0; 0365 else 0366 likelihood = LIK; 0367 end 0368 DynareOptions.warning_for_steadystate = 1; 0369 % ------------------------------------------------------------------------------ 0370 % Adds prior if necessary 0371 % ------------------------------------------------------------------------------ 0372 lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); 0373 fval = (likelihood-lnprior);