dsge_likelihood

 Evaluates the posterior kernel of a dsge model.

0001 function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
0002 % Evaluates the posterior kernel of a dsge model.
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},@var{DLIK},@var{AHess}] =} dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults},@var{derivatives_flag})
0006 %! @anchor{dsge_likelihood}
0007 %! @sp 1
0008 %! Evaluates the posterior kernel of a dsge model.
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 %! @item derivates_flag
0028 %! Integer scalar, flag for analytical derivatives of the likelihood.
0029 %! @end table
0030 %! @sp 2
0031 %! @strong{Outputs}
0032 %! @sp 1
0033 %! @table @ @var
0034 %! @item fval
0035 %! Double scalar, value of (minus) the likelihood.
0036 %! @item exit_flag
0037 %! Integer scalar, equal to zero if the routine return with a penalty (one otherwise).
0038 %! @item ys
0039 %! Vector of doubles, steady state level for the endogenous variables.
0040 %! @item trend_coeffs
0041 %! Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
0042 %! @item info
0043 %! Integer scalar, error code.
0044 %! @table @ @code
0045 %! @item info==0
0046 %! No error.
0047 %! @item info==1
0048 %! The model doesn't determine the current variables uniquely.
0049 %! @item info==2
0050 %! MJDGGES returned an error code.
0051 %! @item info==3
0052 %! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
0053 %! @item info==4
0054 %! Blanchard & Kahn conditions are not satisfied: indeterminacy.
0055 %! @item info==5
0056 %! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
0057 %! @item info==6
0058 %! The jacobian evaluated at the deterministic steady state is complex.
0059 %! @item info==19
0060 %! The steadystate routine thrown an exception (inconsistent deep parameters).
0061 %! @item info==20
0062 %! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations).
0063 %! @item info==21
0064 %! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.
0065 %! @item info==22
0066 %! The steady has NaNs.
0067 %! @item info==23
0068 %! M_.params has been updated in the steadystate routine and has complex valued scalars.
0069 %! @item info==24
0070 %! M_.params has been updated in the steadystate routine and has some NaNs.
0071 %! @item info==30
0072 %! Ergodic variance can't be computed.
0073 %! @item info==41
0074 %! At least one parameter is violating a lower bound condition.
0075 %! @item info==42
0076 %! At least one parameter is violating an upper bound condition.
0077 %! @item info==43
0078 %! The covariance matrix of the structural innovations is not positive definite.
0079 %! @item info==44
0080 %! The covariance matrix of the measurement errors is not positive definite.
0081 %! @item info==45
0082 %! Likelihood is not a number (NaN).
0083 %! @item info==45
0084 %! Likelihood is a complex valued number.
0085 %! @end table
0086 %! @item Model
0087 %! Matlab's structure describing the model (initialized by dynare, see @ref{M_}).
0088 %! @item DynareOptions
0089 %! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
0090 %! @item BayesInfo
0091 %! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
0092 %! @item DynareResults
0093 %! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
0094 %! @item DLIK
0095 %! Vector of doubles, score of the likelihood.
0096 %! @item AHess
0097 %! Matrix of doubles, asymptotic hessian matrix.
0098 %! @end table
0099 %! @sp 2
0100 %! @strong{This function is called by:}
0101 %! @sp 1
0102 %! @ref{dynare_estimation_1}, @ref{mode_check}
0103 %! @sp 2
0104 %! @strong{This function calls:}
0105 %! @sp 1
0106 %! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{schur_statespace_transformation}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
0107 %! @end deftypefn
0108 %@eod:
0109 
0110 % Copyright (C) 2004-2011 Dynare Team
0111 %
0112 % This file is part of Dynare.
0113 %
0114 % Dynare is free software: you can redistribute it and/or modify
0115 % it under the terms of the GNU General Public License as published by
0116 % the Free Software Foundation, either version 3 of the License, or
0117 % (at your option) any later version.
0118 %
0119 % Dynare is distributed in the hope that it will be useful,
0120 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0121 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0122 % GNU General Public License for more details.
0123 %
0124 % You should have received a copy of the GNU General Public License
0125 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0126 
0127 % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
0128 
0129 % Declaration of the penalty as a persistent variable.
0130 
0131 % Persistent variable 'penalty' is used to compute an endogenous penalty to
0132 % the value 'fval' when various conditions are encountered. These conditions
0133 % set also 'exit_flag' equal to 0 instead of 1.  It is only when
0134 % dsge_likelihood() is called by an optimizer called by
0135 % dynare_estimation_1() that 'exit_flag' is ignored and penalized 'fval' is
0136 % actually used.
0137 % In that case, 'penalty' is properly initialized, at the very end of the
0138 % present function, by a call to dsge_likelihood() made in
0139 % initial_estimation_checks(). If a condition triggers exit_flag ==
0140 % 0, initial_estimation_checks() triggers an error.
0141 % In summary, an initial call to the present function, without triggering
0142 % any condition, guarantees that 'penalty' is properly initialized when needed.
0143 
0144 persistent penalty
0145 
0146 % Initialization of the returned variables and others...
0147 fval        = [];
0148 ys          = [];
0149 trend_coeff = [];
0150 exit_flag   = 1;
0151 info        = 0;
0152 singularity_flag = 0;
0153 DLIK        = [];
0154 Hess       = [];
0155 
0156 if DynareOptions.estimation_dll
0157     [fval,exit_flag,ys,trend_coeff,info,params,H,Q] ...
0158         = logposterior(xparam1,DynareDataset, DynareOptions,Model, ...
0159                           EstimatedParameters,BayesInfo,DynareResults);
0160     Model.params = params;
0161     if ~isequal(Model.H,0)
0162         Model.H = H;
0163     end
0164     Model.Sigma_e = Q;
0165     DynareResults.dr.ys = ys;
0166     return
0167 end
0168 
0169 % Set flag related to analytical derivatives.
0170 analytic_derivation = DynareOptions.analytic_derivation;
0171 if nargout==1,
0172     analytic_derivation=0;
0173 end
0174     
0175 %------------------------------------------------------------------------------
0176 % 1. Get the structural parameters & define penalties
0177 %------------------------------------------------------------------------------
0178 
0179 % Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
0180 if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BayesInfo.lb)
0181     k = find(xparam1<BayesInfo.lb);
0182     fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
0183     exit_flag = 0;
0184     info = 41;
0185     if analytic_derivation,
0186         DLIK=ones(length(xparam1),1);
0187     end
0188     return
0189 end
0190 
0191 % Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
0192 if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BayesInfo.ub)
0193     k = find(xparam1>BayesInfo.ub);
0194     fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
0195     exit_flag = 0;
0196     info = 42;
0197     if analytic_derivation,
0198         DLIK=ones(length(xparam1),1);
0199     end
0200     return
0201 end
0202 
0203 % Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
0204 Q = Model.Sigma_e;
0205 H = Model.H;
0206 for i=1:EstimatedParameters.nvx
0207     k =EstimatedParameters.var_exo(i,1);
0208     Q(k,k) = xparam1(i)*xparam1(i);
0209 end
0210 offset = EstimatedParameters.nvx;
0211 if EstimatedParameters.nvn
0212     for i=1:EstimatedParameters.nvn
0213         k = EstimatedParameters.var_endo(i,1);
0214         H(k,k) = xparam1(i+offset)*xparam1(i+offset);
0215     end
0216     offset = offset+EstimatedParameters.nvn;
0217 else
0218     H = zeros(DynareDataset.info.nvobs);
0219 end
0220 
0221 % Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
0222 if EstimatedParameters.ncx
0223     for i=1:EstimatedParameters.ncx
0224         k1 =EstimatedParameters.corrx(i,1);
0225         k2 =EstimatedParameters.corrx(i,2);
0226         Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
0227         Q(k2,k1) = Q(k1,k2);
0228     end
0229     % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
0230     [CholQ,testQ] = chol(Q);
0231     if testQ
0232         % 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.
0233         a = diag(eig(Q));
0234         k = find(a < 0);
0235         if k > 0
0236             fval = penalty+sum(-a(k));
0237             exit_flag = 0;
0238             info = 43;
0239             return
0240         end
0241     end
0242     offset = offset+EstimatedParameters.ncx;
0243 end
0244 
0245 % Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
0246 if EstimatedParameters.ncn
0247     for i=1:EstimatedParameters.ncn
0248         k1 = DynareOptions.lgyidx2varobs(EstimatedParameters.corrn(i,1));
0249         k2 = DynareOptions.lgyidx2varobs(EstimatedParameters.corrn(i,2));
0250         H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
0251         H(k2,k1) = H(k1,k2);
0252     end
0253     % Try to compute the cholesky decomposition of H (possible iff H is positive definite)
0254     [CholH,testH] = chol(H);
0255     if testH
0256         % 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.
0257         a = diag(eig(H));
0258         k = find(a < 0);
0259         if k > 0
0260             fval = penalty+sum(-a(k));
0261             exit_flag = 0;
0262             info = 44;
0263             return
0264         end
0265     end
0266     offset = offset+EstimatedParameters.ncn;
0267 end
0268 
0269 % Update estimated structural parameters in Mode.params.
0270 if EstimatedParameters.np > 0
0271     Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end);
0272 end
0273 
0274 % Update Model.Sigma_e and Model.H.
0275 Model.Sigma_e = Q;
0276 Model.H = H;
0277 
0278 %------------------------------------------------------------------------------
0279 % 2. call model setup & reduction program
0280 %------------------------------------------------------------------------------
0281 
0282 % Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
0283 [T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
0284 
0285 % Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
0286 if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 22 || info(1) == 24
0287     fval = penalty+1;
0288     info = info(1);
0289     exit_flag = 0;
0290     if analytic_derivation,
0291         DLIK=ones(length(xparam1),1);
0292     end
0293     return
0294 elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21  || info(1) == 23
0295     fval = penalty+info(2);
0296     info = info(1);
0297     exit_flag = 0;
0298     if analytic_derivation,
0299         DLIK=ones(length(xparam1),1);
0300     end
0301     return
0302 end
0303 
0304 % Define a vector of indices for the observed variables. Is this really usefull?...
0305 BayesInfo.mf = BayesInfo.mf1;
0306 
0307 % Define the constant vector of the measurement equation.
0308 if DynareOptions.noconstant
0309     constant = zeros(DynareDataset.info.nvobs,1);
0310 else
0311     if DynareOptions.loglinear
0312         constant = log(SteadyState(BayesInfo.mfys));
0313     else
0314         constant = SteadyState(BayesInfo.mfys);
0315     end
0316 end
0317 
0318 % Define the deterministic linear trend of the measurement equation.
0319 if BayesInfo.with_trend
0320     trend_coeff = zeros(DynareDataset.info.nvobs,1);
0321     t = DynareOptions.trend_coeffs;
0322     for i=1:length(t)
0323         if ~isempty(t{i})
0324             trend_coeff(i) = evalin('base',t{i});
0325         end
0326     end
0327     trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs];
0328 else
0329     trend = repmat(constant,1,DynareDataset.info.ntobs);
0330 end
0331 
0332 % Get needed informations for kalman filter routines.
0333 start = DynareOptions.presample+1;
0334 Z = BayesInfo.mf; % old mf
0335 no_missing_data_flag = ~DynareDataset.missing.state;
0336 mm = length(T); % old np
0337 pp = DynareDataset.info.nvobs;
0338 rr = length(Q);
0339 kalman_tol = DynareOptions.kalman_tol;
0340 riccati_tol = DynareOptions.riccati_tol;
0341 Y   = DynareDataset.data-trend;
0342 
0343 %------------------------------------------------------------------------------
0344 % 3. Initial condition of the Kalman filter
0345 %------------------------------------------------------------------------------
0346 kalman_algo = DynareOptions.kalman_algo;
0347 
0348 % resetting measurement errors covariance matrix for univariate filters
0349 if (kalman_algo == 2) || (kalman_algo == 4)
0350     if isequal(H,0)
0351         H = zeros(pp,1);
0352         mmm = mm;
0353     else
0354         if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
0355             H = diag(H);
0356             mmm = mm; 
0357         else
0358             Z = [Z, eye(pp)];
0359             T = blkdiag(T,zeros(pp));
0360             Q = blkdiag(Q,H);
0361             R = blkdiag(R,eye(pp));
0362             Pstar = blkdiag(Pstar,H);
0363             Pinf  = blckdiag(Pinf,zeros(pp));
0364             H = zeros(pp,1);
0365             mmm   = mm+pp;
0366         end
0367     end
0368 end
0369 
0370 
0371 diffuse_periods = 0;
0372 correlated_errors_have_been_checked = 0;
0373 singular_diffuse_filter = 0;
0374 switch DynareOptions.lik_init
0375   case 1% Standard initialization with the steady state of the state equation.
0376     if kalman_algo~=2
0377         % Use standard kalman filter except if the univariate filter is explicitely choosen.
0378         kalman_algo = 1;
0379     end
0380     if DynareOptions.lyapunov_fp == 1
0381         Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 3, R);
0382     elseif DynareOptions.lyapunov_db == 1
0383         Pstar = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol);
0384     elseif DynareOptions.lyapunov_srs == 1
0385         Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 4, R);
0386     else
0387         Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
0388     end;
0389     Pinf  = [];
0390     a     = zeros(mm,1);
0391     Zflag = 0;
0392   case 2% Initialization with large numbers on the diagonal of the covariance matrix if the states (for non stationary models).
0393     if kalman_algo ~= 2
0394         % Use standard kalman filter except if the univariate filter is explicitely choosen.
0395         kalman_algo = 1;
0396     end
0397     Pstar = DynareOptions.Harvey_scale_factor*eye(mm);
0398     Pinf  = [];
0399     a     = zeros(mm,1);
0400     Zflag = 0;
0401   case 3% Diffuse Kalman filter (Durbin and Koopman)
0402         % Use standard kalman filter except if the univariate filter is explicitely choosen.
0403     if kalman_algo == 0
0404         kalman_algo = 3;
0405     elseif ~((kalman_algo == 3) || (kalman_algo == 4)) 
0406             error(['diffuse filter: options_.kalman_algo can only be equal ' ...
0407                    'to 0 (default), 3 or 4'])
0408     end
0409 
0410     [Z,T,R,QT,Pstar,Pinf] = schur_statespace_transformation(Z,T,R,Q,DynareOptions.qz_criterium);
0411     Zflag = 1;
0412     % Run diffuse kalman filter on first periods.
0413     if (kalman_algo==3)
0414         % Multivariate Diffuse Kalman Filter
0415         if no_missing_data_flag
0416             [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
0417                                                        zeros(mm,1), Pinf, Pstar, ...
0418                                                        kalman_tol, riccati_tol, DynareOptions.presample, ...
0419                                                        T,R,Q,H,Z,mm,pp,rr);
0420         else
0421             [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
0422                                                               Y, 1, size(Y,2), ...
0423                                                               zeros(mm,1), Pinf, Pstar, ...
0424                                                               kalman_tol, riccati_tol, DynareOptions.presample, ...
0425                                                               T,R,Q,H,Z,mm,pp,rr);
0426         end
0427         diffuse_periods = length(dlik);
0428         if isinf(dLIK)
0429             % Go to univariate diffuse filter if singularity problem.
0430             singular_diffuse_filter = 1;
0431         end
0432     end
0433     if singular_diffuse_filter || (kalman_algo==4)
0434         % Univariate Diffuse Kalman Filter
0435         if isequal(H,0)
0436             H1 = zeros(pp,1);
0437             mmm = mm;
0438         else
0439             if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
0440                 H1 = diag(H);
0441                 mmm = mm; 
0442             else
0443                 Z = [Z, eye(pp)];
0444                 T = blkdiag(T,zeros(pp));
0445                 Q = blkdiag(Q,H);
0446                 R = blkdiag(R,eye(pp));
0447                 Pstar = blkdiag(Pstar,H);
0448                 Pinf  = blckdiag(Pinf,zeros(pp));
0449                 H1 = zeros(pp,1);
0450                 mmm   = mm+pp;
0451             end
0452         end
0453         % no need to test again for correlation elements
0454         correlated_errors_have_been_checked = 1;
0455 
0456         [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,...
0457                                                         DynareDataset.missing.number_of_observations,...
0458                                                         DynareDataset.missing.no_more_missing_observations, ...
0459                                                         Y, 1, size(Y,2), ...
0460                                                         zeros(mmm,1), Pinf, Pstar, ...
0461                                                         kalman_tol, riccati_tol, DynareOptions.presample, ...
0462                                                         T,R,Q,H1,Z,mmm,pp,rr);
0463         diffuse_periods = length(dlik);
0464     end
0465   case 4% Start from the solution of the Riccati equation.
0466     if kalman_algo ~= 2
0467         kalman_algo = 1;
0468     end
0469     if isequal(H,0)
0470         [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))));
0471     else
0472         [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))),H);
0473     end
0474     if err
0475         disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
0476         DynareOptions.lik_init = 1;
0477         Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
0478     end
0479     Pinf  = [];
0480   otherwise
0481     error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!')
0482 end
0483 
0484 if analytic_derivation
0485     no_DLIK = 0;
0486     full_Hess = analytic_derivation==2;
0487     asy_Hess = analytic_derivation==-2;
0488     if asy_Hess,
0489         analytic_derivation=1;
0490     end
0491     DLIK = [];
0492     AHess = [];
0493     if nargin<8 || isempty(derivatives_info)
0494         [A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
0495         if ~isempty(EstimatedParameters.var_exo)
0496             indexo=EstimatedParameters.var_exo(:,1);
0497         else
0498             indexo=[];
0499         end
0500         if ~isempty(EstimatedParameters.param_vals)
0501             indparam=EstimatedParameters.param_vals(:,1);
0502         else
0503             indparam=[];
0504         end
0505 
0506         if full_Hess,
0507             [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
0508         else
0509             [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
0510         end
0511     else
0512         DT = derivatives_info.DT;
0513         DOm = derivatives_info.DOm;
0514         DYss = derivatives_info.DYss;
0515         if isfield(derivatives_info,'full_Hess'),
0516             full_Hess = derivatives_info.full_Hess;
0517         end
0518         if full_Hess,
0519         D2T = derivatives_info.D2T;
0520         D2Om = derivatives_info.D2Om;
0521         D2Yss = derivatives_info.D2Yss;
0522         end
0523         if isfield(derivatives_info,'no_DLIK'),
0524             no_DLIK = derivatives_info.no_DLIK;
0525         end
0526         clear('derivatives_info');
0527     end
0528     iv = DynareResults.dr.restrict_var_list;
0529     DYss = [zeros(size(DYss,1),offset) DYss];
0530     DT = DT(iv,iv,:);
0531     DOm = DOm(iv,iv,:);
0532     DYss = DYss(iv,:);
0533     DH=zeros([size(H),length(xparam1)]);
0534     DQ=zeros([size(Q),length(xparam1)]);
0535     DP=zeros([size(T),length(xparam1)]);
0536     if full_Hess,
0537         for j=1:size(D2Yss,1),
0538         tmp(j,:,:) = blkdiag(zeros(offset,offset), squeeze(D2Yss(j,:,:)));
0539         end
0540         D2Yss = tmp;
0541         D2T = D2T(iv,iv,:,:);
0542         D2Om = D2Om(iv,iv,:,:);
0543         D2Yss = D2Yss(iv,:,:);
0544         D2H=zeros([size(H),length(xparam1),length(xparam1)]);
0545         D2P=zeros([size(T),length(xparam1),length(xparam1)]);
0546     end
0547     for i=1:EstimatedParameters.nvx
0548         k =EstimatedParameters.var_exo(i,1);
0549         DQ(k,k,i) = 2*sqrt(Q(k,k));
0550         dum =  lyapunov_symm(T,DOm(:,:,i),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
0551         kk = find(abs(dum) < 1e-12);
0552         dum(kk) = 0;
0553         DP(:,:,i)=dum;
0554         if full_Hess
0555         for j=1:i,
0556             dum =  lyapunov_symm(T,D2Om(:,:,i,j),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
0557             kk = (abs(dum) < 1e-12);
0558             dum(kk) = 0;
0559             D2P(:,:,i,j)=dum;
0560             D2P(:,:,j,i)=dum;
0561         end
0562         end
0563     end
0564     offset = EstimatedParameters.nvx;
0565     for i=1:EstimatedParameters.nvn
0566         k = EstimatedParameters.var_endo(i,1);
0567         DH(k,k,i+offset) = 2*sqrt(H(k,k));
0568         if full_Hess
0569         D2H(k,k,i+offset,i+offset) = 2;
0570         end
0571     end
0572     offset = offset + EstimatedParameters.nvn;
0573     for j=1:EstimatedParameters.np
0574         dum =  lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
0575         kk = find(abs(dum) < 1e-12);
0576         dum(kk) = 0;
0577         DP(:,:,j+offset)=dum;
0578         if full_Hess
0579         DTj = DT(:,:,j+offset);
0580         DPj = dum;
0581         for i=1:j,
0582             DTi = DT(:,:,i+offset);
0583             DPi = DP(:,:,i+offset);
0584             D2Tij = D2T(:,:,i,j);
0585             D2Omij = D2Om(:,:,i,j);
0586             tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij;
0587             dum = lyapunov_symm(T,tmp,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
0588             dum(abs(dum)<1.e-12) = 0;
0589             D2P(:,:,i+offset,j+offset) = dum;
0590             D2P(:,:,j+offset,i+offset) = dum;
0591         end
0592         end
0593     end
0594     if analytic_derivation==1,
0595         analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP};
0596     else
0597         analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P};
0598     end
0599 else
0600     analytic_deriv_info={0};
0601 end
0602 
0603 %------------------------------------------------------------------------------
0604 % 4. Likelihood evaluation
0605 %------------------------------------------------------------------------------
0606 
0607 if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
0608     if no_missing_data_flag
0609         if DynareOptions.block == 1
0610             [err, LIK] = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
0611             mexErrCheck('block_kalman_filter', err);
0612         else
0613             [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
0614                                 a,Pstar, ...
0615                                 kalman_tol, riccati_tol, ...
0616                                 DynareOptions.presample, ...
0617                                 T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
0618                                 analytic_deriv_info{:}); 
0619                             
0620         end
0621     else
0622         [LIK,lik] = missing_observations_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
0623                                                a, Pstar, ...
0624                                                kalman_tol, DynareOptions.riccati_tol, ...
0625                                                DynareOptions.presample, ...
0626                                                T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods);
0627     end
0628     if analytic_derivation,
0629         LIK1=LIK;
0630         LIK=LIK1{1};
0631     end
0632     if isinf(LIK)
0633         if kalman_algo == 1
0634             kalman_algo = 2;
0635         else
0636             kalman_algo = 4;
0637         end
0638     else
0639         if DynareOptions.lik_init==3
0640             LIK = LIK + dLIK;
0641             if analytic_derivation==0 && nargout==2,
0642                 lik = [dlik; lik];
0643             end
0644         end
0645     end
0646 end
0647 
0648 if (kalman_algo==2) || (kalman_algo==4)
0649     % Univariate Kalman Filter
0650     % resetting measurement error covariance matrix when necessary                                                           %
0651     if ~correlated_errors_have_been_checked
0652         if isequal(H,0)
0653             H = zeros(pp,1);
0654             mmm = mm;
0655             if analytic_derivation,
0656                 DH = zeros(pp,length(xparam1));
0657             end
0658         else
0659             if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
0660                 H = diag(H);
0661                 mmm = mm;
0662                 if analytic_derivation,
0663                     for j=1:pp,
0664                         tmp(j,:)=DH(j,j,:);
0665                     end
0666                     DH=tmp;
0667                 end
0668             else
0669                 Z = [Z, eye(pp)];
0670                 T = blkdiag(T,zeros(pp));
0671                 Q = blkdiag(Q,H);
0672                 R = blkdiag(R,eye(pp));
0673                 Pstar = blkdiag(Pstar,H);
0674                 Pinf  = blckdiag(Pinf,zeros(pp));
0675                 H = zeros(pp,1);
0676                 mmm   = mm+pp;
0677             end
0678         end
0679         if analytic_derivation,
0680             analytic_deriv_info{5}=DH;
0681         end
0682     end
0683 
0684     [LIK, lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
0685                                        a,Pstar, ...
0686                                        DynareOptions.kalman_tol, ...
0687                                        DynareOptions.riccati_tol, ...
0688                                        DynareOptions.presample, ...
0689                                        T,Q,R,H,Z,mmm,pp,rr,Zflag,diffuse_periods,analytic_deriv_info{:});
0690     if analytic_derivation,
0691         LIK1=LIK;
0692         LIK=LIK1{1};
0693     end
0694     if DynareOptions.lik_init==3
0695         LIK = LIK+dLIK;
0696         if analytic_derivation==0 && nargout==2,
0697             lik = [dlik; lik];
0698         end
0699     end
0700 end
0701 
0702 if analytic_derivation
0703     if no_DLIK==0
0704         DLIK = LIK1{2};
0705         %                 [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
0706     end
0707     if full_Hess,
0708         Hess = -LIK1{3};
0709         %                     [Hess, DLL] = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol);
0710         %                     Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss);
0711     end
0712     if asy_Hess,
0713         [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
0714     end
0715 end
0716 
0717 
0718 if isnan(LIK)
0719     info = 45;
0720     exit_flag = 0;
0721     return
0722 end
0723 if imag(LIK)~=0
0724     likelihood = penalty;
0725 else
0726     likelihood = LIK;
0727 end
0728 
0729 % ------------------------------------------------------------------------------
0730 % 5. Adds prior if necessary
0731 % ------------------------------------------------------------------------------
0732 if analytic_derivation
0733     if full_Hess,
0734         [lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
0735         Hess = Hess - d2lnprior;
0736     else
0737         [lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
0738     end
0739     if no_DLIK==0
0740         DLIK = DLIK - dlnprior';
0741     end
0742 else
0743     lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
0744 end
0745 fval    = (likelihood-lnprior);
0746 
0747 % Update DynareOptions.kalman_algo.
0748 DynareOptions.kalman_algo = kalman_algo;
0749 
0750 % Update the penalty.
0751 penalty = fval;
0752 
0753 if analytic_derivation==0 && nargout==2,
0754     lik=lik(start:end,:);
0755     DLIK=[-lnprior; lik(:)];
0756 end