1function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults] = non_linear_dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults) 2 3% Evaluates the posterior kernel of a dsge model using a non linear filter. 4% 5% INPUTS 6% - xparam1 [double] n×1 vector, estimated parameters. 7% - DynareDataset [struct] Matlab's structure containing the dataset (initialized by dynare, aka dataset_). 8% - DatasetInfo [struct] Matlab's structure describing the dataset (initialized by dynare, aka dataset_info). 9% - DynareOptions [struct] Matlab's structure describing the options (initialized by dynare, aka options_). 10% - Model [struct] Matlab's structure describing the Model (initialized by dynare, aka M_). 11% - EstimatedParameters [struct] Matlab's structure describing the estimated_parameters (initialized by dynare, aka estim_params_). 12% - BayesInfo [struct] Matlab's structure describing the priors (initialized by dynare,aka bayesopt_). 13% - BoundsInfo [struct] Matlab's structure specifying the bounds on the paramater values (initialized by dynare,aka bayesopt_). 14% - DynareResults [struct] Matlab's structure gathering the results (initialized by dynare,aka oo_). 15% 16% OUTPUTS 17% - fval [double] scalar, value of the likelihood or posterior kernel. 18% - info [integer] 4×1 vector, informations resolution of the model and evaluation of the likelihood. 19% - exit_flag [integer] scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood). 20% - DLIK [double] Empty array. 21% - Hess [double] Empty array. 22% - ys [double] Empty array. 23% - trend_coeff [double] Empty array. 24% - Model [struct] Updated Model structure described in INPUTS section. 25% - DynareOptions [struct] Updated DynareOptions structure described in INPUTS section. 26% - BayesInfo [struct] See INPUTS section. 27% - DynareResults [struct] Updated DynareResults structure described in INPUTS section. 28 29% Copyright (C) 2010-2019 Dynare Team 30% 31% This file is part of Dynare. 32% 33% Dynare is free software: you can redistribute it and/or modify 34% it under the terms of the GNU General Public License as published by 35% the Free Software Foundation, either version 3 of the License, or 36% (at your option) any later version. 37% 38% Dynare is distributed in the hope that it will be useful, 39% but WITHOUT ANY WARRANTY; without even the implied warranty of 40% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 41% GNU General Public License for more details. 42% 43% You should have received a copy of the GNU General Public License 44% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 45 46% Initialization of the returned arguments. 47fval = []; 48ys = []; 49trend_coeff = []; 50exit_flag = 1; 51DLIK = []; 52Hess = []; 53 54% Ensure that xparam1 is a column vector. 55xparam1 = xparam1(:); 56 57% Issue an error if loglinear option is used. 58if DynareOptions.loglinear 59 error('non_linear_dsge_likelihood: It is not possible to use a non linear filter with the option loglinear!') 60end 61 62%------------------------------------------------------------------------------ 63% 1. Get the structural parameters & define penalties 64%------------------------------------------------------------------------------ 65 66% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain. 67if isestimation(DynareOptions) && (DynareOptions.mode_compute~=1) && any(xparam1<BoundsInfo.lb) 68 k = find(xparam1(:) < BoundsInfo.lb); 69 fval = Inf; 70 exit_flag = 0; 71 info(1) = 41; 72 info(4) = sum((BoundsInfo.lb(k)-xparam1(k)).^2); 73 return 74end 75 76% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain. 77if isestimation(DynareOptions) && (DynareOptions.mode_compute~=1) && any(xparam1>BoundsInfo.ub) 78 k = find(xparam1(:)>BoundsInfo.ub); 79 fval = Inf; 80 exit_flag = 0; 81 info(1) = 42; 82 info(4) = sum((xparam1(k)-BoundsInfo.ub(k)).^2); 83 return 84end 85 86Model = set_all_parameters(xparam1,EstimatedParameters,Model); 87 88Q = Model.Sigma_e; 89H = Model.H; 90 91if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances') 92 [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)); 93 if ~Q_is_positive_definite 94 fval = Inf; 95 exit_flag = 0; 96 info(1) = 43; 97 info(4) = penalty; 98 return 99 end 100 if isfield(EstimatedParameters,'calibrated_covariances') 101 correct_flag=check_consistency_covariances(Q); 102 if ~correct_flag 103 penalty = sum(Q(EstimatedParameters.calibrated_covariances.position).^2); 104 fval = Inf; 105 exit_flag = 0; 106 info(1) = 71; 107 info(4) = penalty; 108 return 109 end 110 end 111 112end 113 114if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME') 115 [H_is_positive_definite, penalty] = ispd(H(EstimatedParameters.H_entries_to_check_for_positive_definiteness,EstimatedParameters.H_entries_to_check_for_positive_definiteness)); 116 if ~H_is_positive_definite 117 fval = Inf; 118 exit_flag = 0; 119 info(1) = 44; 120 info(4) = penalty; 121 return 122 end 123 if isfield(EstimatedParameters,'calibrated_covariances_ME') 124 correct_flag=check_consistency_covariances(H); 125 if ~correct_flag 126 penalty = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2); 127 fval = Inf; 128 exit_flag = 0; 129 info(1) = 72; 130 info(4) = penalty; 131 return 132 end 133 end 134 135end 136 137%------------------------------------------------------------------------------ 138% 2. call model setup & reduction program 139%------------------------------------------------------------------------------ 140 141% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R). 142[dr, info, Model, DynareOptions, DynareResults] = resol(0, Model, DynareOptions, DynareResults); 143 144if info(1) 145 if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 || ... 146 info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ... 147 info(1) == 81 || info(1) == 84 || info(1) == 85 148 %meaningful second entry of output that can be used 149 fval = Inf; 150 info(4) = info(2); 151 exit_flag = 0; 152 return 153 else 154 fval = Inf; 155 info(4) = 0.1; 156 exit_flag = 0; 157 return 158 end 159end 160 161% Define a vector of indices for the observed variables. Is this really usefull?... 162BayesInfo.mf = BayesInfo.mf1; 163 164% Get needed informations for kalman filter routines. 165start = DynareOptions.presample+1; 166Y = transpose(DynareDataset.data); 167 168%------------------------------------------------------------------------------ 169% 3. Initial condition of the Kalman filter 170%------------------------------------------------------------------------------ 171 172mf0 = BayesInfo.mf0; 173mf1 = BayesInfo.mf1; 174restrict_variables_idx = dr.restrict_var_list; 175state_variables_idx = restrict_variables_idx(mf0); 176number_of_state_variables = length(mf0); 177 178ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx)); 179ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx); 180ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx)); 181ReducedForm.Q = Q; 182ReducedForm.H = H; 183ReducedForm.mf0 = mf0; 184ReducedForm.mf1 = mf1; 185 186if DynareOptions.k_order_solver && ~(DynareOptions.particle.pruning && DynareOptions.order==2) 187 ReducedForm.use_k_order_solver = true; 188 ReducedForm.dr = dr; 189else 190 ReducedForm.use_k_order_solver = false; 191 ReducedForm.ghx = dr.ghx(restrict_variables_idx,:); 192 ReducedForm.ghu = dr.ghu(restrict_variables_idx,:); 193 ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:); 194 ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); 195 ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); 196end 197 198% Set initial condition. 199switch DynareOptions.particle.initialization 200 case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model. 201 StateVectorMean = ReducedForm.constant(mf0); 202 StateVectorVariance = lyapunov_symm(dr.ghx(mf0,:), dr.ghu(mf0,:)*Q*dr.ghu(mf0,:)', DynareOptions.lyapunov_fixed_point_tol, ... 203 DynareOptions.qz_criterium, DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); 204 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). 205 StateVectorMean = ReducedForm.constant(mf0); 206 old_DynareOptionsperiods = DynareOptions.periods; 207 DynareOptions.periods = 5000; 208 old_DynareOptionspruning = DynareOptions.pruning; 209 DynareOptions.pruning = DynareOptions.particle.pruning; 210 y_ = simult(DynareResults.steady_state, dr,Model,DynareOptions,DynareResults); 211 y_ = y_(state_variables_idx,2001:5000); 212 StateVectorVariance = cov(y_'); 213 DynareOptions.periods = old_DynareOptionsperiods; 214 DynareOptions.pruning = old_DynareOptionspruning; 215 clear('old_DynareOptionsperiods','y_'); 216 case 3% Initial state vector covariance is a diagonal matrix (to be used 217 % if model has stochastic trends). 218 StateVectorMean = ReducedForm.constant(mf0); 219 StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables); 220 otherwise 221 error('Unknown initialization option!') 222end 223ReducedForm.StateVectorMean = StateVectorMean; 224ReducedForm.StateVectorVariance = StateVectorVariance; 225 226%------------------------------------------------------------------------------ 227% 4. Likelihood evaluation 228%------------------------------------------------------------------------------ 229DynareOptions.warning_for_steadystate = 0; 230[s1,s2] = get_dynare_random_generator_state(); 231LIK = feval(DynareOptions.particle.algorithm, ReducedForm, Y, start, DynareOptions.particle, DynareOptions.threads, DynareOptions, Model); 232set_dynare_random_generator_state(s1,s2); 233if imag(LIK) 234 likelihood = Inf; 235 info(1) = 46; 236 info(4) = 0.1; 237 exit_flag = 0; 238elseif isnan(LIK) 239 likelihood = Inf; 240 info(1) = 45; 241 info(4) = 0.1; 242 exit_flag = 0; 243else 244 likelihood = LIK; 245end 246DynareOptions.warning_for_steadystate = 1; 247% ------------------------------------------------------------------------------ 248% Adds prior if necessary 249% ------------------------------------------------------------------------------ 250lnprior = priordens(xparam1(:),BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); 251fval = (likelihood-lnprior); 252 253if isnan(fval) 254 fval = Inf; 255 info(1) = 47; 256 info(4) = 0.1; 257 exit_flag = 0; 258 return 259end 260 261if ~isreal(fval) 262 fval = Inf; 263 info(1) = 48; 264 info(4) = 0.1; 265 exit_flag = 0; 266 return 267end 268 269if isinf(LIK) 270 fval = Inf; 271 info(1) = 50; 272 info(4) = 0.1; 273 exit_flag = 0; 274 return 275end 276