2012-04-27 15:57:58 +02:00
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)
2011-09-19 16:38:38 +02:00
% Evaluates the posterior kernel of a dsge model.
%@info:
2011-12-26 15:42:35 +01:00
%! @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})
2011-12-26 15:39:27 +01:00
%! @anchor{dsge_likelihood}
2011-09-19 16:38:38 +02:00
%! @sp 1
%! Evaluates the posterior kernel of a dsge model.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item xparam1
%! Vector of doubles, current values for the estimated parameters.
%! @item DynareDataset
%! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item Model
%! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
%! @item EstimatedParamemeters
%! Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @item derivates_flag
%! Integer scalar, flag for analytical derivatives of the likelihood.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item fval
%! Double scalar, value of (minus) the likelihood.
%! @item exit_flag
%! Integer scalar, equal to zero if the routine return with a penalty (one otherwise).
%! @item ys
%! Vector of doubles, steady state level for the endogenous variables.
%! @item trend_coeffs
%! Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
%! @item info
%! Integer scalar, error code.
%! @table @ @code
%! @item info==0
%! No error.
%! @item info==1
%! The model doesn't determine the current variables uniquely.
%! @item info==2
%! MJDGGES returned an error code.
%! @item info==3
%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
%! @item info==4
%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==6
%! The jacobian evaluated at the deterministic steady state is complex.
%! @item info==19
%! The steadystate routine thrown an exception (inconsistent deep parameters).
%! @item info==20
%! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations).
%! @item info==21
%! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.
%! @item info==22
%! The steady has NaNs.
%! @item info==23
%! M_.params has been updated in the steadystate routine and has complex valued scalars.
%! @item info==24
%! M_.params has been updated in the steadystate routine and has some NaNs.
%! @item info==30
%! Ergodic variance can't be computed.
%! @item info==41
%! At least one parameter is violating a lower bound condition.
%! @item info==42
%! At least one parameter is violating an upper bound condition.
%! @item info==43
%! The covariance matrix of the structural innovations is not positive definite.
%! @item info==44
%! The covariance matrix of the measurement errors is not positive definite.
%! @item info==45
%! Likelihood is not a number (NaN).
2012-06-14 15:11:36 +02:00
%! @item info==46
2011-09-19 16:38:38 +02:00
%! Likelihood is a complex valued number.
2012-06-14 15:11:36 +02:00
%! @item info==47
%! Posterior kernel is not a number (logged prior density is NaN)
%! @item info==48
%! Posterior kernel is a complex valued number (logged prior density is complex).
2011-09-19 16:38:38 +02:00
%! @end table
%! @item Model
%! Matlab's structure describing the model (initialized by dynare, see @ref{M_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @item DLIK
%! Vector of doubles, score of the likelihood.
%! @item AHess
%! Matrix of doubles, asymptotic hessian matrix.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{dynare_estimation_1}, @ref{mode_check}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
2011-09-19 17:01:24 +02:00
%! @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}
2011-09-19 16:38:38 +02:00
%! @end deftypefn
%@eod:
2008-08-01 14:40:33 +02:00
2012-06-08 18:22:34 +02:00
% Copyright (C) 2004-2012 Dynare Team
2008-08-01 14:40:33 +02:00
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
2011-09-19 16:38:38 +02:00
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
2012-08-28 12:17:07 +02:00
global objective_function_penalty_base
2011-09-19 16:38:38 +02:00
2011-09-21 08:30:26 +02:00
% Initialization of the returned variables and others...
fval = [ ] ;
ys = [ ] ;
trend_coeff = [ ] ;
exit_flag = 1 ;
info = 0 ;
2011-11-02 11:10:58 +01:00
singularity_flag = 0 ;
2012-01-08 21:55:02 +01:00
DLIK = [ ] ;
2012-04-27 15:57:58 +02:00
Hess = [ ] ;
2012-01-08 21:55:02 +01:00
if DynareOptions . estimation_dll
[ fval , exit_flag , ys , trend_coeff , info , params , H , Q ] ...
= logposterior ( xparam1 , DynareDataset , DynareOptions , Model , ...
EstimatedParameters , BayesInfo , DynareResults ) ;
2013-04-18 18:27:54 +02:00
mexErrCheck ( ' logposterior' , exit_flag ) ;
2012-01-08 21:55:02 +01:00
Model . params = params ;
if ~ isequal ( Model . H , 0 )
Model . H = H ;
end
Model . Sigma_e = Q ;
DynareResults . dr . ys = ys ;
return
end
2011-09-21 08:30:26 +02:00
2011-09-19 16:38:38 +02:00
% Set flag related to analytical derivatives.
2012-04-27 15:57:58 +02:00
analytic_derivation = DynareOptions . analytic_derivation ;
2012-06-08 11:33:33 +02:00
if analytic_derivation && DynareOptions . loglinear
error ( ' The analytic_derivation and loglinear options are not compatible' )
end
2012-04-27 15:57:58 +02:00
if nargout == 1 ,
2011-03-18 11:05:40 +01:00
analytic_derivation = 0 ;
2011-09-19 16:38:38 +02:00
end
2012-06-14 15:11:36 +02:00
2012-07-05 10:14:10 +02:00
if analytic_derivation ,
kron_flag = DynareOptions . analytic_derivation_mode ;
end
2009-12-16 18:17:34 +01:00
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
2011-09-19 16:38:38 +02:00
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if ~ isequal ( DynareOptions . mode_compute , 1 ) && any ( xparam1 < BayesInfo . lb )
k = find ( xparam1 < BayesInfo . lb ) ;
2012-08-28 12:17:07 +02:00
fval = objective_function_penalty_base + sum ( ( BayesInfo . lb ( k ) - xparam1 ( k ) ) .^ 2 ) ;
2011-09-19 16:38:38 +02:00
exit_flag = 0 ;
2013-05-02 14:31:29 +02:00
info = 41 ;
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
DLIK = ones ( length ( xparam1 ) , 1 ) ;
end
2011-09-19 16:38:38 +02:00
return
2009-12-16 18:17:34 +01:00
end
2011-09-19 16:38:38 +02:00
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if ~ isequal ( DynareOptions . mode_compute , 1 ) && any ( xparam1 > BayesInfo . ub )
k = find ( xparam1 > BayesInfo . ub ) ;
2012-08-28 12:17:07 +02:00
fval = objective_function_penalty_base + sum ( ( xparam1 ( k ) - BayesInfo . ub ( k ) ) .^ 2 ) ;
2011-09-19 16:38:38 +02:00
exit_flag = 0 ;
2013-05-02 14:31:29 +02:00
info = 42 ;
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
DLIK = ones ( length ( xparam1 ) , 1 ) ;
end
2011-09-19 16:38:38 +02:00
return
2009-12-16 18:17:34 +01:00
end
2011-09-19 16:38:38 +02:00
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
2012-06-07 15:12:10 +02:00
Model = set_all_parameters ( xparam1 , EstimatedParameters , Model ) ;
2011-09-19 16:38:38 +02:00
Q = Model . Sigma_e ;
H = Model . H ;
2012-06-07 15:12:10 +02:00
% Test if Q is positive definite.
2011-09-19 16:38:38 +02:00
if EstimatedParameters . ncx
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
2005-09-11 11:04:41 +02:00
[ CholQ , testQ ] = chol ( Q ) ;
2011-09-19 16:38:38 +02:00
if testQ
% 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.
2009-12-16 18:17:34 +01:00
a = diag ( eig ( Q ) ) ;
k = find ( a < 0 ) ;
if k > 0
2012-08-28 12:17:07 +02:00
fval = objective_function_penalty_base + sum ( - a ( k ) ) ;
2011-09-19 16:38:38 +02:00
exit_flag = 0 ;
2009-12-16 18:17:34 +01:00
info = 43 ;
return
end
2005-09-11 11:04:41 +02:00
end
2009-12-16 18:17:34 +01:00
end
2011-09-19 16:38:38 +02:00
2012-06-07 15:12:10 +02:00
% Test if H is positive definite.
2011-09-19 16:38:38 +02:00
if EstimatedParameters . ncn
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
2005-09-11 11:04:41 +02:00
[ CholH , testH ] = chol ( H ) ;
if testH
2011-12-26 17:44:04 +01:00
% 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.
2009-12-16 18:17:34 +01:00
a = diag ( eig ( H ) ) ;
k = find ( a < 0 ) ;
if k > 0
2012-08-28 12:17:07 +02:00
fval = objective_function_penalty_base + sum ( - a ( k ) ) ;
2011-09-19 16:38:38 +02:00
exit_flag = 0 ;
2009-12-16 18:17:34 +01:00
info = 44 ;
return
end
2005-09-11 11:04:41 +02:00
end
2009-12-16 18:17:34 +01:00
end
2011-09-19 16:38:38 +02:00
2009-12-16 18:17:34 +01:00
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
2011-01-26 13:50:11 +01:00
2011-09-19 16:38:38 +02:00
% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
2011-09-22 11:17:31 +02:00
[ T , R , SteadyState , info , Model , DynareOptions , DynareResults ] = dynare_resolve ( Model , DynareOptions , DynareResults , ' restrict' ) ;
2011-09-19 16:38:38 +02:00
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
2012-08-06 12:00:03 +02:00
if info ( 1 ) == 1 || info ( 1 ) == 2 || info ( 1 ) == 5 || info ( 1 ) == 7 || info ( 1 ) ...
2012-09-14 17:01:05 +02:00
== 8 || info ( 1 ) == 22 || info ( 1 ) == 24 || info ( 1 ) == 19
2012-08-28 12:17:07 +02:00
fval = objective_function_penalty_base + 1 ;
2011-09-19 16:38:38 +02:00
info = info ( 1 ) ;
exit_flag = 0 ;
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
DLIK = ones ( length ( xparam1 ) , 1 ) ;
end
2009-12-16 18:17:34 +01:00
return
2012-09-14 17:01:05 +02:00
elseif info ( 1 ) == 3 || info ( 1 ) == 4 || info ( 1 ) == 6 || info ( 1 ) == 20 || info ( 1 ) == 21 || info ( 1 ) == 23
2012-08-28 12:17:07 +02:00
fval = objective_function_penalty_base + info ( 2 ) ;
2011-09-19 16:38:38 +02:00
info = info ( 1 ) ;
exit_flag = 0 ;
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
DLIK = ones ( length ( xparam1 ) , 1 ) ;
end
2009-12-16 18:17:34 +01:00
return
end
2011-09-19 16:38:38 +02:00
% Define a vector of indices for the observed variables. Is this really usefull?...
BayesInfo . mf = BayesInfo . mf1 ;
% Define the constant vector of the measurement equation.
if DynareOptions . noconstant
2011-09-19 17:50:23 +02:00
constant = zeros ( DynareDataset . info . nvobs , 1 ) ;
2011-09-19 16:38:38 +02:00
else
if DynareOptions . loglinear
constant = log ( SteadyState ( BayesInfo . mfys ) ) ;
2009-12-16 18:17:34 +01:00
else
2011-09-19 16:38:38 +02:00
constant = SteadyState ( BayesInfo . mfys ) ;
2009-12-16 18:17:34 +01:00
end
end
2011-09-19 16:38:38 +02:00
% Define the deterministic linear trend of the measurement equation.
if BayesInfo . with_trend
2011-09-19 17:50:23 +02:00
trend_coeff = zeros ( DynareDataset . info . nvobs , 1 ) ;
2011-09-19 16:38:38 +02:00
t = DynareOptions . trend_coeffs ;
2007-10-01 11:39:32 +02:00
for i = 1 : length ( t )
2009-12-16 18:17:34 +01:00
if ~ isempty ( t { i } )
trend_coeff ( i ) = evalin ( ' base' , t { i } ) ;
end
2005-09-11 11:04:41 +02:00
end
2011-09-19 16:38:38 +02:00
trend = repmat ( constant , 1 , DynareDataset . info . ntobs ) + trend_coeff * [ 1 : DynareDataset . info . ntobs ] ;
2009-12-16 18:17:34 +01:00
else
2011-09-19 16:38:38 +02:00
trend = repmat ( constant , 1 , DynareDataset . info . ntobs ) ;
2009-12-16 18:17:34 +01:00
end
2011-09-19 16:38:38 +02:00
% Get needed informations for kalman filter routines.
start = DynareOptions . presample + 1 ;
Z = BayesInfo . mf ; % old mf
no_missing_data_flag = ~ DynareDataset . missing . state ;
mm = length ( T ) ; % old np
pp = DynareDataset . info . nvobs ;
rr = length ( Q ) ;
kalman_tol = DynareOptions . kalman_tol ;
riccati_tol = DynareOptions . riccati_tol ;
Y = DynareDataset . data - trend ;
2009-12-16 18:17:34 +01:00
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
2011-09-19 16:38:38 +02:00
kalman_algo = DynareOptions . kalman_algo ;
2011-10-22 15:26:07 +02:00
% resetting measurement errors covariance matrix for univariate filters
if ( kalman_algo == 2 ) || ( kalman_algo == 4 )
if isequal ( H , 0 )
2012-04-27 15:02:43 +02:00
H = zeros ( pp , 1 ) ;
2011-11-02 11:10:58 +01:00
mmm = mm ;
2011-10-22 15:26:07 +02:00
else
if all ( all ( abs ( H - diag ( diag ( H ) ) ) < 1e-14 ) ) % ie, the covariance matrix is diagonal...
H = diag ( H ) ;
2012-06-14 15:11:36 +02:00
mmm = mm ;
2011-10-22 15:26:07 +02:00
else
2011-11-02 11:10:58 +01:00
Z = [ Z , eye ( pp ) ] ;
T = blkdiag ( T , zeros ( pp ) ) ;
Q = blkdiag ( Q , H ) ;
R = blkdiag ( R , eye ( pp ) ) ;
Pstar = blkdiag ( Pstar , H ) ;
Pinf = blckdiag ( Pinf , zeros ( pp ) ) ;
2012-04-27 15:02:43 +02:00
H = zeros ( pp , 1 ) ;
2011-11-02 11:10:58 +01:00
mmm = mm + pp ;
2011-10-22 15:26:07 +02:00
end
end
end
2011-09-20 10:13:01 +02:00
diffuse_periods = 0 ;
2012-01-22 18:37:29 +01:00
correlated_errors_have_been_checked = 0 ;
2012-01-22 18:59:19 +01:00
singular_diffuse_filter = 0 ;
2011-09-19 16:38:38 +02:00
switch DynareOptions . lik_init
case 1 % Standard initialization with the steady state of the state equation.
if kalman_algo ~= 2
% Use standard kalman filter except if the univariate filter is explicitely choosen.
2009-12-16 18:17:34 +01:00
kalman_algo = 1 ;
end
2012-03-06 12:03:23 +01:00
if DynareOptions . lyapunov_fp == 1
2012-04-20 19:23:00 +02:00
Pstar = lyapunov_symm ( T , Q , DynareOptions . lyapunov_fixed_point_tol , DynareOptions . lyapunov_complex_threshold , 3 , R ) ;
elseif DynareOptions . lyapunov_db == 1
Pstar = disclyap_fast ( T , R * Q * R ' , DynareOptions . lyapunov_doubling_tol ) ;
elseif DynareOptions . lyapunov_srs == 1
Pstar = lyapunov_symm ( T , Q , DynareOptions . lyapunov_fixed_point_tol , DynareOptions . lyapunov_complex_threshold , 4 , R ) ;
2012-03-06 12:03:23 +01:00
else
Pstar = lyapunov_symm ( T , R * Q * R ' , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
end ;
2011-09-19 16:38:38 +02:00
Pinf = [ ] ;
a = zeros ( mm , 1 ) ;
Zflag = 0 ;
case 2 % Initialization with large numbers on the diagonal of the covariance matrix if the states (for non stationary models).
2009-12-16 18:17:34 +01:00
if kalman_algo ~= 2
2011-09-19 16:38:38 +02:00
% Use standard kalman filter except if the univariate filter is explicitely choosen.
2009-12-16 18:17:34 +01:00
kalman_algo = 1 ;
end
2011-09-19 16:38:38 +02:00
Pstar = DynareOptions . Harvey_scale_factor * eye ( mm ) ;
Pinf = [ ] ;
a = zeros ( mm , 1 ) ;
Zflag = 0 ;
case 3 % Diffuse Kalman filter (Durbin and Koopman)
% Use standard kalman filter except if the univariate filter is explicitely choosen.
2012-01-22 18:37:29 +01:00
if kalman_algo == 0
2009-12-16 18:17:34 +01:00
kalman_algo = 3 ;
2012-06-14 15:11:36 +02:00
elseif ~ ( ( kalman_algo == 3 ) || ( kalman_algo == 4 ) )
2012-01-22 18:37:29 +01:00
error ( [ ' diffuse filter: options_.kalman_algo can only be equal ' ...
' to 0 (default), 3 or 4' ] )
2009-12-16 18:17:34 +01:00
end
2012-01-22 18:37:29 +01:00
2011-09-19 16:38:38 +02:00
[ Z , T , R , QT , Pstar , Pinf ] = schur_statespace_transformation ( Z , T , R , Q , DynareOptions . qz_criterium ) ;
Zflag = 1 ;
% Run diffuse kalman filter on first periods.
if ( kalman_algo == 3 )
% Multivariate Diffuse Kalman Filter
if no_missing_data_flag
2012-04-27 15:02:43 +02:00
[ dLIK , dlik , a , Pstar ] = kalman_filter_d ( Y , 1 , size ( Y , 2 ) , ...
2011-09-19 16:38:38 +02:00
zeros ( mm , 1 ) , Pinf , Pstar , ...
kalman_tol , riccati_tol , DynareOptions . presample , ...
T , R , Q , H , Z , mm , pp , rr ) ;
else
2012-04-27 15:02:43 +02:00
[ dLIK , dlik , a , Pstar ] = missing_observations_kalman_filter_d ( DynareDataset . missing . aindex , DynareDataset . missing . number_of_observations , DynareDataset . missing . no_more_missing_observations , ...
2011-09-19 16:38:38 +02:00
Y , 1 , size ( Y , 2 ) , ...
zeros ( mm , 1 ) , Pinf , Pstar , ...
kalman_tol , riccati_tol , DynareOptions . presample , ...
T , R , Q , H , Z , mm , pp , rr ) ;
end
2012-04-27 15:02:43 +02:00
diffuse_periods = length ( dlik ) ;
2011-09-19 16:38:38 +02:00
if isinf ( dLIK )
% Go to univariate diffuse filter if singularity problem.
2012-01-22 18:59:19 +01:00
singular_diffuse_filter = 1 ;
2011-09-19 16:38:38 +02:00
end
end
2012-01-22 18:59:19 +01:00
if singular_diffuse_filter || ( kalman_algo == 4 )
2011-09-19 16:38:38 +02:00
% Univariate Diffuse Kalman Filter
2012-01-22 18:37:29 +01:00
if isequal ( H , 0 )
2012-04-27 15:02:43 +02:00
H1 = zeros ( pp , 1 ) ;
2012-01-22 18:37:29 +01:00
mmm = mm ;
else
if all ( all ( abs ( H - diag ( diag ( H ) ) ) < 1e-14 ) ) % ie, the covariance matrix is diagonal...
2012-01-22 18:59:19 +01:00
H1 = diag ( H ) ;
2012-06-14 15:11:36 +02:00
mmm = mm ;
2011-11-02 11:10:58 +01:00
else
2012-01-22 18:37:29 +01:00
Z = [ Z , eye ( pp ) ] ;
T = blkdiag ( T , zeros ( pp ) ) ;
Q = blkdiag ( Q , H ) ;
R = blkdiag ( R , eye ( pp ) ) ;
Pstar = blkdiag ( Pstar , H ) ;
Pinf = blckdiag ( Pinf , zeros ( pp ) ) ;
2012-04-27 15:02:43 +02:00
H1 = zeros ( pp , 1 ) ;
2012-01-22 18:37:29 +01:00
mmm = mm + pp ;
2011-11-02 11:10:58 +01:00
end
end
2012-01-22 18:37:29 +01:00
% no need to test again for correlation elements
correlated_errors_have_been_checked = 1 ;
2012-04-27 15:02:43 +02:00
[ dLIK , dlik , a , Pstar ] = univariate_kalman_filter_d ( DynareDataset . missing . aindex , ...
2012-01-22 18:37:29 +01:00
DynareDataset . missing . number_of_observations , ...
DynareDataset . missing . no_more_missing_observations , ...
Y , 1 , size ( Y , 2 ) , ...
zeros ( mmm , 1 ) , Pinf , Pstar , ...
kalman_tol , riccati_tol , DynareOptions . presample , ...
2012-01-22 18:59:19 +01:00
T , R , Q , H1 , Z , mmm , pp , rr ) ;
2012-04-27 15:02:43 +02:00
diffuse_periods = length ( dlik ) ;
2011-09-19 16:38:38 +02:00
end
2012-09-14 17:05:35 +02:00
if isnan ( dLIK ) ,
info = 45 ;
2012-10-01 14:23:21 +02:00
fval = objective_function_penalty_base + 100 ;
2012-09-14 17:05:35 +02:00
exit_flag = 0 ;
return
end
2011-09-19 16:38:38 +02:00
case 4 % Start from the solution of the Riccati equation.
2011-11-02 11:10:58 +01:00
if kalman_algo ~= 2
2011-06-23 23:39:15 +02:00
kalman_algo = 1 ;
end
2011-06-28 14:46:43 +02:00
if isequal ( H , 0 )
2012-05-30 17:59:05 +02:00
[ err , Pstar ] = kalman_steady_state ( transpose ( T ) , R * Q * transpose ( R ) , transpose ( build_selection_matrix ( Z , mm , length ( Z ) ) ) ) ;
2011-06-28 14:46:43 +02:00
else
2012-05-30 17:59:05 +02:00
[ err , Pstar ] = kalman_steady_state ( transpose ( T ) , R * Q * transpose ( R ) , transpose ( build_selection_matrix ( Z , mm , length ( Z ) ) ) , H ) ;
2011-06-28 14:46:43 +02:00
end
if err
2011-12-26 15:39:27 +01:00
disp ( [ ' dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!' ] ) ;
2011-09-19 16:38:38 +02:00
DynareOptions . lik_init = 1 ;
Pstar = lyapunov_symm ( T , R * Q * R ' , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
2011-06-28 14:46:43 +02:00
end
2012-06-08 14:23:18 +02:00
Pinf = [ ] ;
2012-05-30 17:59:05 +02:00
a = zeros ( mm , 1 ) ;
Zflag = 0 ;
2011-09-19 16:38:38 +02:00
otherwise
2011-12-26 15:39:27 +01:00
error ( ' dsge_likelihood:: Unknown initialization approach for the Kalman filter!' )
2009-12-16 18:17:34 +01:00
end
2011-06-28 14:46:43 +02:00
2012-06-08 08:50:07 +02:00
if analytic_derivation ,
offset = EstimatedParameters . nvx ;
offset = offset + EstimatedParameters . nvn ;
offset = offset + EstimatedParameters . ncx ;
offset = offset + EstimatedParameters . ncn ;
2011-04-15 15:26:33 +02:00
no_DLIK = 0 ;
2012-04-27 15:57:58 +02:00
full_Hess = analytic_derivation == 2 ;
asy_Hess = analytic_derivation == - 2 ;
2012-06-08 14:23:18 +02:00
outer_product_gradient = analytic_derivation == - 1 ;
2012-04-27 15:57:58 +02:00
if asy_Hess ,
analytic_derivation = 1 ;
end
2012-06-08 14:23:18 +02:00
if outer_product_gradient ,
analytic_derivation = 1 ;
end
2011-04-15 15:26:33 +02:00
DLIK = [ ] ;
AHess = [ ] ;
2012-08-21 16:00:55 +02:00
iv = DynareResults . dr . restrict_var_list ;
2011-11-07 09:19:36 +01:00
if nargin < 8 || isempty ( derivatives_info )
2011-09-19 16:38:38 +02:00
[ A , B , nou , nou , Model , DynareOptions , DynareResults ] = dynare_resolve ( Model , DynareOptions , DynareResults ) ;
if ~ isempty ( EstimatedParameters . var_exo )
indexo = EstimatedParameters . var_exo ( : , 1 ) ;
2011-04-12 18:18:18 +02:00
else
indexo = [ ] ;
end
2011-09-19 16:38:38 +02:00
if ~ isempty ( EstimatedParameters . param_vals )
indparam = EstimatedParameters . param_vals ( : , 1 ) ;
2011-04-12 18:18:18 +02:00
else
indparam = [ ] ;
end
2011-11-05 11:08:00 +01:00
if full_Hess ,
2012-08-21 16:00:55 +02:00
[ dum , DT , DOm , DYss , dum2 , D2T , D2Om , D2Yss ] = getH ( A , B , Model , DynareResults , DynareOptions , kron_flag , indparam , indexo , iv ) ;
clear dum dum2 ;
2011-11-05 11:08:00 +01:00
else
2012-08-21 16:00:55 +02:00
[ dum , DT , DOm , DYss ] = getH ( A , B , Model , DynareResults , DynareOptions , kron_flag , indparam , indexo , iv ) ;
2011-11-05 11:08:00 +01:00
end
2011-03-18 11:05:40 +01:00
else
2012-08-21 16:00:55 +02:00
DT = derivatives_info . DT ( iv , iv , : ) ;
DOm = derivatives_info . DOm ( iv , iv , : ) ;
DYss = derivatives_info . DYss ( iv , : ) ;
2011-11-05 11:08:00 +01:00
if isfield ( derivatives_info , ' full_Hess' ) ,
full_Hess = derivatives_info . full_Hess ;
end
if full_Hess ,
D2T = derivatives_info . D2T ;
D2Om = derivatives_info . D2Om ;
D2Yss = derivatives_info . D2Yss ;
end
if isfield ( derivatives_info , ' no_DLIK' ) ,
2011-04-15 15:26:33 +02:00
no_DLIK = derivatives_info . no_DLIK ;
end
2011-09-19 16:38:38 +02:00
clear ( ' derivatives_info' ) ;
2011-03-18 11:05:40 +01:00
end
2011-04-12 18:18:18 +02:00
DYss = [ zeros ( size ( DYss , 1 ) , offset ) DYss ] ;
2012-06-08 15:26:14 +02:00
DH = zeros ( [ length ( H ) , length ( H ) , length ( xparam1 ) ] ) ;
2011-03-18 11:05:40 +01:00
DQ = zeros ( [ size ( Q ) , length ( xparam1 ) ] ) ;
DP = zeros ( [ size ( T ) , length ( xparam1 ) ] ) ;
2011-11-05 11:08:00 +01:00
if full_Hess ,
for j = 1 : size ( D2Yss , 1 ) ,
tmp ( j , : , : ) = blkdiag ( zeros ( offset , offset ) , squeeze ( D2Yss ( j , : , : ) ) ) ;
end
D2Yss = tmp ;
2012-08-21 16:00:55 +02:00
D2H = sparse ( size ( D2Om , 1 ) , size ( D2Om , 2 ) ) ; %zeros([size(H),length(xparam1),length(xparam1)]);
D2P = sparse ( size ( D2Om , 1 ) , size ( D2Om , 2 ) ) ; %zeros([size(T),length(xparam1),length(xparam1)]);
jcount = 0 ;
2011-11-05 11:08:00 +01:00
end
2012-09-14 17:07:38 +02:00
if DynareOptions . lik_init == 1 ,
2011-09-19 16:38:38 +02:00
for i = 1 : EstimatedParameters . nvx
k = EstimatedParameters . var_exo ( i , 1 ) ;
2011-03-18 11:05:40 +01:00
DQ ( k , k , i ) = 2 * sqrt ( Q ( k , k ) ) ;
2011-09-19 16:38:38 +02:00
dum = lyapunov_symm ( T , DOm ( : , : , i ) , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
2012-08-21 16:00:55 +02:00
% kk = find(abs(dum) < 1e-12);
% dum(kk) = 0;
2011-03-18 11:05:40 +01:00
DP ( : , : , i ) = dum ;
2011-11-05 11:08:00 +01:00
if full_Hess
for j = 1 : i ,
2012-08-21 16:00:55 +02:00
jcount = jcount + 1 ;
dum = lyapunov_symm ( T , dyn_unvech ( D2Om ( : , jcount ) ) , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
% kk = (abs(dum) < 1e-12);
% dum(kk) = 0;
D2P ( : , jcount ) = dyn_vech ( dum ) ;
% D2P(:,:,j,i)=dum;
2011-11-05 11:08:00 +01:00
end
end
2011-03-18 11:05:40 +01:00
end
2012-09-14 17:07:38 +02:00
end
2011-09-19 16:38:38 +02:00
offset = EstimatedParameters . nvx ;
for i = 1 : EstimatedParameters . nvn
k = EstimatedParameters . var_endo ( i , 1 ) ;
2011-03-18 11:05:40 +01:00
DH ( k , k , i + offset ) = 2 * sqrt ( H ( k , k ) ) ;
2011-11-05 11:08:00 +01:00
if full_Hess
D2H ( k , k , i + offset , i + offset ) = 2 ;
end
2011-03-18 11:05:40 +01:00
end
2011-09-19 16:38:38 +02:00
offset = offset + EstimatedParameters . nvn ;
2012-09-14 17:07:38 +02:00
if DynareOptions . lik_init == 1 ,
2011-09-19 16:38:38 +02:00
for j = 1 : EstimatedParameters . np
dum = lyapunov_symm ( T , DT ( : , : , j + offset ) * Pstar * T ' + T * Pstar * DT ( : , : , j + offset ) ' + DOm ( : , : , j + offset ) , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
2012-08-21 16:00:55 +02:00
% kk = find(abs(dum) < 1e-12);
% dum(kk) = 0;
2011-03-18 11:05:40 +01:00
DP ( : , : , j + offset ) = dum ;
2011-11-05 11:08:00 +01:00
if full_Hess
DTj = DT ( : , : , j + offset ) ;
DPj = dum ;
2012-08-02 14:57:20 +02:00
for i = 1 : j + offset ,
2012-08-21 16:00:55 +02:00
jcount = jcount + 1 ;
2012-08-02 14:57:20 +02:00
DTi = DT ( : , : , i ) ;
DPi = DP ( : , : , i ) ;
2012-08-21 16:00:55 +02:00
D2Tij = reshape ( D2T ( : , jcount ) , size ( T ) ) ;
D2Omij = dyn_unvech ( D2Om ( : , jcount ) ) ;
2011-11-05 11:08:00 +01:00
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 ;
dum = lyapunov_symm ( T , tmp , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
2012-08-21 16:00:55 +02:00
% dum(abs(dum)<1.e-12) = 0;
D2P ( : , jcount ) = dyn_vech ( dum ) ;
% D2P(:,:,j+offset,i) = dum;
2011-11-05 11:08:00 +01:00
end
end
2011-03-18 11:05:40 +01:00
end
2012-09-14 17:07:38 +02:00
end
2012-04-27 15:57:58 +02:00
if analytic_derivation == 1 ,
2012-06-08 14:23:18 +02:00
analytic_deriv_info = { analytic_derivation , DT , DYss , DOm , DH , DP , asy_Hess } ;
2012-04-27 15:57:58 +02:00
else
analytic_deriv_info = { analytic_derivation , DT , DYss , DOm , DH , DP , D2T , D2Yss , D2Om , D2H , D2P } ;
2012-08-21 16:00:55 +02:00
clear DT DYss DOm DH DP D2T D2Yss D2Om D2H D2P ,
2012-04-27 15:57:58 +02:00
end
else
analytic_deriv_info = { 0 } ;
2011-03-18 11:05:40 +01:00
end
2009-12-16 18:17:34 +01:00
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
2011-09-19 16:38:38 +02:00
if ( ( kalman_algo == 1 ) || ( kalman_algo == 3 ) ) % Multivariate Kalman Filter
2009-12-16 18:17:34 +01:00
if no_missing_data_flag
2012-06-11 10:39:25 +02:00
if DynareOptions . block
2011-10-28 22:36:40 +02:00
[ 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 ) ;
2011-09-23 18:21:04 +02:00
mexErrCheck ( ' block_kalman_filter' , err ) ;
2011-09-20 14:18:31 +02:00
else
2012-04-27 15:02:43 +02:00
[ LIK , lik ] = kalman_filter ( Y , diffuse_periods + 1 , size ( Y , 2 ) , ...
2011-09-23 18:09:06 +02:00
a , Pstar , ...
kalman_tol , riccati_tol , ...
DynareOptions . presample , ...
2012-04-27 15:57:58 +02:00
T , Q , R , H , Z , mm , pp , rr , Zflag , diffuse_periods , ...
2012-06-14 15:11:36 +02:00
analytic_deriv_info { : } ) ;
2011-03-18 11:05:40 +01:00
end
2009-12-16 18:17:34 +01:00
else
2012-06-11 10:39:25 +02:00
if 0 %DynareOptions.block
[ err , LIK , lik ] = block_kalman_filter ( DynareDataset . missing . aindex , DynareDataset . missing . number_of_observations , DynareDataset . missing . no_more_missing_observations , ...
T , R , Q , H , Pstar , Y , start , Z , kalman_tol , riccati_tol , Model . nz_state_var , Model . n_diag , Model . nobs_non_statevar ) ;
else
[ 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 ) , ...
2011-09-19 16:38:38 +02:00
a , Pstar , ...
kalman_tol , DynareOptions . riccati_tol , ...
DynareOptions . presample , ...
T , Q , R , H , Z , mm , pp , rr , Zflag , diffuse_periods ) ;
2012-09-27 14:47:21 +02:00
end
2009-12-16 18:17:34 +01:00
end
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
LIK1 = LIK ;
LIK = LIK1 { 1 } ;
2012-06-08 14:23:18 +02:00
lik1 = lik ;
lik = lik1 { 1 } ;
2012-04-27 15:57:58 +02:00
end
2009-12-16 18:17:34 +01:00
if isinf ( LIK )
2012-09-27 14:47:21 +02:00
if DynareOptions . use_univariate_filters_if_singularity_is_detected
if kalman_algo == 1
kalman_algo = 2 ;
else
kalman_algo = 4 ;
end
2011-11-02 14:02:12 +01:00
else
2012-09-27 14:47:21 +02:00
if isinf ( LIK )
info = 66 ;
fval = objective_function_penalty_base + 1 ;
exit_flag = 0 ;
return
end
2011-11-02 14:02:12 +01:00
end
2009-12-16 18:17:34 +01:00
else
2011-09-19 16:38:38 +02:00
if DynareOptions . lik_init == 3
LIK = LIK + dLIK ;
2012-04-27 15:02:43 +02:00
if analytic_derivation == 0 && nargout == 2 ,
lik = [ dlik ; lik ] ;
end
2009-12-16 18:17:34 +01:00
end
end
end
2011-09-19 16:38:38 +02:00
2012-01-22 18:37:29 +01:00
if ( kalman_algo == 2 ) || ( kalman_algo == 4 )
2011-10-22 15:26:07 +02:00
% Univariate Kalman Filter
2012-06-14 15:11:36 +02:00
% resetting measurement error covariance matrix when necessary %
2012-01-22 18:37:29 +01:00
if ~ correlated_errors_have_been_checked
2011-10-22 15:26:07 +02:00
if isequal ( H , 0 )
2013-05-17 23:54:17 +02:00
H1 = zeros ( pp , 1 ) ;
2011-11-02 11:10:58 +01:00
mmm = mm ;
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
DH = zeros ( pp , length ( xparam1 ) ) ;
end
2011-10-22 15:26:07 +02:00
else
if all ( all ( abs ( H - diag ( diag ( H ) ) ) < 1e-14 ) ) % ie, the covariance matrix is diagonal...
2013-05-17 23:54:17 +02:00
H1 = diag ( H ) ;
2011-11-02 11:10:58 +01:00
mmm = mm ;
2012-06-08 15:26:14 +02:00
clear tmp
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
for j = 1 : pp ,
tmp ( j , : ) = DH ( j , j , : ) ;
end
DH = tmp ;
end
2011-10-22 15:26:07 +02:00
else
2011-11-02 11:10:58 +01:00
Z = [ Z , eye ( pp ) ] ;
T = blkdiag ( T , zeros ( pp ) ) ;
Q = blkdiag ( Q , H ) ;
R = blkdiag ( R , eye ( pp ) ) ;
Pstar = blkdiag ( Pstar , H ) ;
Pinf = blckdiag ( Pinf , zeros ( pp ) ) ;
2013-05-17 23:54:17 +02:00
H1 = zeros ( pp , 1 ) ;
2011-11-02 11:10:58 +01:00
mmm = mm + pp ;
2011-10-22 15:26:07 +02:00
end
end
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
analytic_deriv_info { 5 } = DH ;
end
2009-12-16 18:17:34 +01:00
end
2011-10-22 15:26:07 +02:00
2012-04-27 15:02:43 +02:00
[ 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 ) , ...
2011-09-19 16:38:38 +02:00
a , Pstar , ...
DynareOptions . kalman_tol , ...
DynareOptions . riccati_tol , ...
DynareOptions . presample , ...
2013-05-17 23:54:17 +02:00
T , Q , R , H1 , Z , mmm , pp , rr , Zflag , diffuse_periods , analytic_deriv_info { : } ) ;
2012-04-27 15:57:58 +02:00
if analytic_derivation ,
LIK1 = LIK ;
LIK = LIK1 { 1 } ;
2012-06-08 14:23:18 +02:00
lik1 = lik ;
lik = lik1 { 1 } ;
2012-04-27 15:57:58 +02:00
end
2011-09-19 16:38:38 +02:00
if DynareOptions . lik_init == 3
LIK = LIK + dLIK ;
2012-04-27 15:02:43 +02:00
if analytic_derivation == 0 && nargout == 2 ,
lik = [ dlik ; lik ] ;
end
2009-12-16 18:17:34 +01:00
end
end
2011-09-19 16:38:38 +02:00
2012-04-27 15:57:58 +02:00
if analytic_derivation
if no_DLIK == 0
DLIK = LIK1 { 2 } ;
% [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
end
2012-06-08 14:23:18 +02:00
if full_Hess ,
2012-04-27 15:57:58 +02:00
Hess = - LIK1 { 3 } ;
% [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);
% Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss);
end
if asy_Hess ,
2012-06-08 15:26:14 +02:00
% if ~((kalman_algo==2) || (kalman_algo==4)),
% [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
% else
2012-06-08 14:23:18 +02:00
Hess = LIK1 { 3 } ;
2012-06-08 15:26:14 +02:00
% end
2012-04-27 15:57:58 +02:00
end
end
2010-02-05 23:18:08 +01:00
if isnan ( LIK )
2011-09-19 16:38:38 +02:00
info = 45 ;
2012-10-24 10:02:17 +02:00
fval = objective_function_penalty_base + 100 ;
2011-09-19 16:38:38 +02:00
exit_flag = 0 ;
2010-02-05 23:18:08 +01:00
return
end
2012-06-14 15:11:36 +02:00
2010-02-05 23:18:08 +01:00
if imag ( LIK ) ~= 0
2012-06-14 15:11:36 +02:00
info = 46 ;
2012-10-24 10:02:17 +02:00
fval = objective_function_penalty_base + 100 ;
2012-06-14 15:11:36 +02:00
exit_flag = 0 ;
return
2009-12-16 18:17:34 +01:00
end
2011-09-19 16:38:38 +02:00
2012-06-14 15:11:36 +02:00
likelihood = LIK ;
2009-12-16 18:17:34 +01:00
% ------------------------------------------------------------------------------
2011-09-19 16:38:38 +02:00
% 5. Adds prior if necessary
2009-12-16 18:17:34 +01:00
% ------------------------------------------------------------------------------
2011-11-07 09:19:36 +01:00
if analytic_derivation
if full_Hess ,
[ lnprior , dlnprior , d2lnprior ] = priordens ( xparam1 , BayesInfo . pshape , BayesInfo . p6 , BayesInfo . p7 , BayesInfo . p3 , BayesInfo . p4 ) ;
2012-04-27 15:57:58 +02:00
Hess = Hess - d2lnprior ;
2011-11-07 09:19:36 +01:00
else
[ lnprior , dlnprior ] = priordens ( xparam1 , BayesInfo . pshape , BayesInfo . p6 , BayesInfo . p7 , BayesInfo . p3 , BayesInfo . p4 ) ;
end
if no_DLIK == 0
DLIK = DLIK - dlnprior ' ;
end
2012-06-08 14:23:18 +02:00
if outer_product_gradient ,
dlik = lik1 { 2 } ;
dlik = [ - dlnprior ; dlik ( start : end , : ) ] ;
Hess = dlik ' * dlik ;
end
2011-11-07 09:19:36 +01:00
else
lnprior = priordens ( xparam1 , BayesInfo . pshape , BayesInfo . p6 , BayesInfo . p7 , BayesInfo . p3 , BayesInfo . p4 ) ;
end
2013-03-17 22:51:23 +01:00
if DynareOptions . endogenous_prior == 1
2013-03-18 23:46:15 +01:00
if DynareOptions . lik_init == 2 || DynareOptions . lik_init == 3
error ( ' Endogenous prior not supported with non-stationary models' )
else
[ lnpriormom ] = endogenous_prior ( Y , Pstar , BayesInfo , H ) ;
fval = ( likelihood - lnprior - lnpriormom ) ;
end
2013-03-17 22:51:23 +01:00
else
fval = ( likelihood - lnprior ) ;
end
2011-09-19 16:38:38 +02:00
2012-06-14 15:11:36 +02:00
if isnan ( fval )
info = 47 ;
2012-10-24 10:02:17 +02:00
fval = objective_function_penalty_base + 100 ;
2012-06-14 15:11:36 +02:00
exit_flag = 0 ;
return
end
if imag ( fval ) ~= 0
info = 48 ;
2012-10-24 10:02:17 +02:00
fval = objective_function_penalty_base + 100 ;
2012-06-14 15:11:36 +02:00
exit_flag = 0 ;
return
end
2011-09-19 16:38:38 +02:00
% Update DynareOptions.kalman_algo.
DynareOptions . kalman_algo = kalman_algo ;
2012-04-27 15:02:43 +02:00
if analytic_derivation == 0 && nargout == 2 ,
lik = lik ( start : end , : ) ;
DLIK = [ - lnprior ; lik ( : ) ] ;
2012-04-27 15:57:58 +02:00
end