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