2012-01-04 14:25:21 +01:00
function [fval,llik,cost_flag,ys,trend_coeff,info] = dsge_likelihood_hh ( xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
2011-11-05 10:46:16 +01:00
% function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
2011-10-03 12:19:41 +02:00
% Evaluates the posterior kernel of a dsge model.
%
% INPUTS
2008-12-12 08:39:03 +01:00
% xparam1 [double] vector of model parameters.
% gend [integer] scalar specifying the number of observations.
% data [double] matrix of data
% data_index [cell] cell of column vectors
% number_of_observations [integer]
2011-10-03 12:19:41 +02:00
% no_more_missing_observations [integer]
% OUTPUTS
2008-12-12 08:39:03 +01:00
% fval : value of the posterior kernel at xparam1.
2011-11-05 10:46:16 +01:00
% llik : probabilities at each time point
2008-12-12 08:39:03 +01:00
% cost_flag : zero if the function returns a penalty, one otherwise.
2008-07-02 13:39:58 +02:00
% ys : steady state of original endogenous variables
% trend_coeff :
2008-12-12 08:39:03 +01:00
% info : vector of informations about the penalty:
% 41: one (many) parameter(s) do(es) not satisfied the lower bound
% 42: one (many) parameter(s) do(es) not satisfied the upper bound
2011-10-03 12:19:41 +02:00
%
2008-07-02 13:39:58 +02:00
% SPECIAL REQUIREMENTS
2008-12-12 08:39:03 +01:00
%
2008-07-02 13:39:58 +02:00
2011-02-04 17:27:33 +01:00
% Copyright (C) 2004-2011 Dynare Team
2008-12-12 08:39:03 +01: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/>.
2008-07-02 13:39:58 +02:00
2011-10-03 12:19:41 +02:00
% Declaration of the penalty as a persistent variable.
2012-01-09 21:23:17 +01:00
% Persistent variable 'penalty' is used to compute an endogenous penalty to
% the value 'fval' when various conditions are encountered. These conditions
% set also 'exit_flag' equal to 0 instead of 1. It is only when
% dsge_likelihood_hh() is called by an newrat() called by
% dynare_estimation_1() that 'exit_flag' is ignored and penalized 'fval' is
% actually used.
% In that case, 'penalty' is properly initialized, at the very end of the
% present function, by a call to dsge_likelihood_hh() made in
% initial_estimation_checks(). If a condition triggers exit_flag ==
% 0, initial_estimation_checks() triggers an error.
% In summary, an initial call to the present function, without triggering
% any condition, guarantees that 'penalty' is properly initialized when needed.
persistent penalty
2011-10-03 12:19:41 +02:00
2011-11-05 10:46:16 +01:00
% Initialization of the returned variables and others...
fval = [ ] ;
ys = [ ] ;
trend_coeff = [ ] ;
cost_flag = 1 ;
2011-02-04 17:17:48 +01:00
llik = NaN ;
2011-11-05 10:46:16 +01:00
info = 0 ;
singularity_flag = 0 ;
2011-10-03 12:19:41 +02:00
if DynareOptions . block == 1
error ( ' DsgeLikelihood_hh:: This routine (called if mode_compute==5) is not compatible with the block option!' )
end
2011-02-04 17:17:48 +01:00
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
2011-10-03 12:19:41 +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 ) ;
fval = penalty + sum ( ( BayesInfo . lb ( k ) - xparam1 ( k ) ) .^ 2 ) ;
exit_flag = 0 ;
2011-02-04 17:17:48 +01:00
info = 41 ;
2011-10-03 12:19:41 +02:00
return
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +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 ) ;
fval = penalty + sum ( ( xparam1 ( k ) - BayesInfo . ub ( k ) ) .^ 2 ) ;
exit_flag = 0 ;
2011-02-04 17:17:48 +01:00
info = 42 ;
2011-10-03 12:19:41 +02:00
return
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Q = Model . Sigma_e ;
H = Model . H ;
for i = 1 : EstimatedParameters . nvx
k = EstimatedParameters . var_exo ( i , 1 ) ;
2008-07-02 13:39:58 +02:00
Q ( k , k ) = xparam1 ( i ) * xparam1 ( i ) ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
offset = EstimatedParameters . nvx ;
if EstimatedParameters . nvn
for i = 1 : EstimatedParameters . nvn
k = EstimatedParameters . var_endo ( i , 1 ) ;
2011-02-04 17:17:48 +01:00
H ( k , k ) = xparam1 ( i + offset ) * xparam1 ( i + offset ) ;
2008-07-02 13:39:58 +02:00
end
2011-10-03 12:19:41 +02:00
offset = offset + EstimatedParameters . nvn ;
else
H = zeros ( DynareDataset . info . nvobs ) ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
if EstimatedParameters . ncx
for i = 1 : EstimatedParameters . ncx
k1 = EstimatedParameters . corrx ( i , 1 ) ;
k2 = EstimatedParameters . corrx ( i , 2 ) ;
2011-02-04 17:17:48 +01:00
Q ( k1 , k2 ) = xparam1 ( i + offset ) * sqrt ( Q ( k1 , k1 ) * Q ( k2 , k2 ) ) ;
Q ( k2 , k1 ) = Q ( k1 , k2 ) ;
2008-07-02 13:39:58 +02:00
end
2011-10-03 12:19:41 +02:00
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
2008-07-02 13:39:58 +02:00
[ CholQ , testQ ] = chol ( Q ) ;
2011-10-03 12:19:41 +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.
2011-02-04 17:17:48 +01:00
a = diag ( eig ( Q ) ) ;
k = find ( a < 0 ) ;
if k > 0
2012-01-08 17:58:22 +01:00
fval = penalty + sum ( - a ( k ) ) ;
2011-10-03 12:19:41 +02:00
exit_flag = 0 ;
2011-02-04 17:17:48 +01:00
info = 43 ;
return
end
2008-07-02 13:39:58 +02:00
end
2011-10-03 12:19:41 +02:00
offset = offset + EstimatedParameters . ncx ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
if EstimatedParameters . ncn
for i = 1 : EstimatedParameters . ncn
k1 = DynareOptions . lgyidx2varobs ( EstimatedParameters . corrn ( i , 1 ) ) ;
k2 = DynareOptions . lgyidx2varobs ( EstimatedParameters . corrn ( i , 2 ) ) ;
2011-02-04 17:17:48 +01:00
H ( k1 , k2 ) = xparam1 ( i + offset ) * sqrt ( H ( k1 , k1 ) * H ( k2 , k2 ) ) ;
H ( k2 , k1 ) = H ( k1 , k2 ) ;
2008-07-02 13:39:58 +02:00
end
2011-10-03 12:19:41 +02:00
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
2008-07-02 13:39:58 +02:00
[ CholH , testH ] = chol ( H ) ;
if testH
2011-10-03 12:19:41 +02:00
% 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.
2011-02-04 17:17:48 +01:00
a = diag ( eig ( H ) ) ;
k = find ( a < 0 ) ;
if k > 0
2012-01-08 17:58:22 +01:00
fval = penalty + sum ( - a ( k ) ) ;
2011-10-03 12:19:41 +02:00
exit_flag = 0 ;
2011-02-04 17:17:48 +01:00
info = 44 ;
return
end
2008-07-02 13:39:58 +02:00
end
2011-10-03 12:19:41 +02:00
offset = offset + EstimatedParameters . ncn ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
% Update estimated structural parameters in Mode.params.
if EstimatedParameters . np > 0
Model . params ( EstimatedParameters . param_vals ( : , 1 ) ) = xparam1 ( offset + 1 : end ) ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
% Update Model.Sigma_e and Model.H.
Model . Sigma_e = Q ;
Model . H = H ;
2011-02-04 17:17:48 +01:00
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
2011-10-03 12:19:41 +02:00
% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
[ T , R , SteadyState , info , Model , DynareOptions , DynareResults ] = dynare_resolve ( Model , DynareOptions , DynareResults , ' restrict' ) ;
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info ( 1 ) == 1 || info ( 1 ) == 2 || info ( 1 ) == 5 || info ( 1 ) == 22 || info ( 1 ) == 24
fval = penalty + 1 ;
info = info ( 1 ) ;
exit_flag = 0 ;
2011-02-04 17:17:48 +01:00
return
2011-10-03 12:19:41 +02:00
elseif info ( 1 ) == 3 || info ( 1 ) == 4 || info ( 1 ) == 6 || info ( 1 ) == 19 || info ( 1 ) == 20 || info ( 1 ) == 21 || info ( 1 ) == 23
fval = penalty + info ( 2 ) ;
info = info ( 1 ) ;
exit_flag = 0 ;
2011-02-04 17:17:48 +01:00
return
end
2011-10-03 12:19:41 +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
constant = zeros ( DynareDataset . info . nvobs , 1 ) ;
else
if DynareOptions . loglinear
constant = log ( SteadyState ( BayesInfo . mfys ) ) ;
2011-02-04 17:17:48 +01:00
else
2011-10-03 12:19:41 +02:00
constant = SteadyState ( BayesInfo . mfys ) ;
2011-02-04 17:17:48 +01:00
end
end
2011-10-03 12:19:41 +02:00
% Define the deterministic linear trend of the measurement equation.
if BayesInfo . with_trend
trend_coeff = zeros ( DynareDataset . info . nvobs , 1 ) ;
t = DynareOptions . trend_coeffs ;
2008-07-02 13:39:58 +02:00
for i = 1 : length ( t )
2011-02-04 17:17:48 +01:00
if ~ isempty ( t { i } )
trend_coeff ( i ) = evalin ( ' base' , t { i } ) ;
end
2008-07-02 13:39:58 +02:00
end
2011-10-03 12:19:41 +02:00
trend = repmat ( constant , 1 , DynareDataset . info . ntobs ) + trend_coeff * [ 1 : DynareDataset . info . ntobs ] ;
2011-02-04 17:17:48 +01:00
else
2011-10-03 12:19:41 +02:00
trend = repmat ( constant , 1 , DynareDataset . info . ntobs ) ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +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 ;
2011-02-04 17:17:48 +01:00
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
2011-10-03 12:19:41 +02:00
kalman_algo = DynareOptions . kalman_algo ;
2011-11-05 10:46:16 +01:00
% resetting measurement errors covariance matrix for univariate filters
if ( kalman_algo == 2 ) || ( kalman_algo == 4 )
if isequal ( H , 0 )
H = zeros ( nobs , 1 ) ;
mmm = mm ;
else
if all ( all ( abs ( H - diag ( diag ( H ) ) ) < 1e-14 ) ) % ie, the covariance matrix is diagonal...
H = diag ( H ) ;
mmm = mm ;
else
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 ) ) ;
H = zeros ( nobs , 1 ) ;
mmm = mm + pp ;
end
end
end
2011-10-03 12:19:41 +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-10-03 12:19:41 +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.
2011-02-04 17:17:48 +01:00
kalman_algo = 1 ;
end
2011-10-03 12:19:41 +02:00
Pstar = lyapunov_symm ( T , R * Q * R ' , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
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).
2011-02-04 17:17:48 +01:00
if kalman_algo ~= 2
2011-10-03 12:19:41 +02:00
% Use standard kalman filter except if the univariate filter is explicitely choosen.
2011-02-04 17:17:48 +01:00
kalman_algo = 1 ;
end
2011-10-03 12:19:41 +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
2011-02-04 17:17:48 +01:00
kalman_algo = 3 ;
2012-01-22 18:37:29 +01:00
elseif ~ ( ( kalman_algo == 3 ) || ( kalman_algo == 4 ) )
error ( [ ' diffuse filter: options_.kalman_algo can only be equal ' ...
' to 0 (default), 3 or 4' ] )
2011-02-04 17:17:48 +01:00
end
2012-01-22 18:37:29 +01:00
2011-10-03 12:19:41 +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
[ dLIK , dlik , a , Pstar ] = kalman_filter_d ( Y , 1 , size ( Y , 2 ) , ...
2012-01-22 18:37:29 +01:00
zeros ( mm , 1 ) , Pinf , Pstar , ...
kalman_tol , riccati_tol , DynareOptions . presample , ...
T , R , Q , H , Z , mm , pp , rr ) ;
2011-10-03 12:19:41 +02:00
else
[ dLIK , dlik , a , Pstar ] = missing_observations_kalman_filter_d ( DynareDataset . missing . aindex , DynareDataset . missing . number_of_observations , DynareDataset . missing . no_more_missing_observations , ...
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
diffuse_periods = length ( dlik ) ;
if isinf ( dLIK )
% Go to univariate diffuse filter if singularity problem.
2012-01-22 18:59:19 +01:00
singular_diffuse_filter
2011-10-03 12:19:41 +02:00
end
2011-02-04 17:17:48 +01:00
end
2012-01-22 18:59:19 +01:00
if singular_diffuse_filter || ( kalman_algo == 4 )
2011-10-03 12:19:41 +02:00
% Univariate Diffuse Kalman Filter
2012-01-22 18:37:29 +01:00
if isequal ( H , 0 )
2012-01-22 18:59:19 +01:00
H1 = zeros ( nobs , 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-01-22 18:37:29 +01:00
mmm = mm ;
2011-11-05 10:46:16 +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-01-22 18:59:19 +01:00
H1 = zeros ( nobs , 1 ) ;
2012-01-22 18:37:29 +01:00
mmm = mm + pp ;
2011-11-05 10:46:16 +01:00
end
2011-02-04 17:17:48 +01:00
end
2012-01-22 18:37:29 +01:00
% no need to test again for correlation elements
correlated_errors_have_been_checked = 1 ;
2012-01-22 18:59:19 +01:00
[ dLIK , dlik , a , Pstar ] = univariate_kalman_filter_d ( DynareDataset . missing . aindex , ...
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 , ...
T , R , Q , H1 , Z , mmm , pp , rr ) ;
2011-10-03 12:19:41 +02:00
diffuse_periods = length ( dlik ) ;
end
case 4 % Start from the solution of the Riccati equation.
2011-11-05 10:46:16 +01:00
if kalman_algo ~= 2
2011-10-03 12:19:41 +02:00
kalman_algo = 1 ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
if isequal ( H , 0 )
2011-11-05 10:46:16 +01:00
[ err , Pstar ] = kalman_steady_state ( transpose ( T ) , R * Q * transpose ( R ) , transpose ( build_selection_matrix ( Z , np , length ( Z ) ) ) ) ;
2011-02-04 17:17:48 +01:00
else
2011-11-05 10:46:16 +01:00
[ err , Pstar ] = kalman_steady_state ( transpose ( T ) , R * Q * transpose ( R ) , transpose ( build_selection_matrix ( Z , np , length ( Z ) ) ) , H ) ;
2011-10-03 12:19:41 +02:00
end
if err
disp ( [ ' DsgeLikelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!' ] ) ;
DynareOptions . lik_init = 1 ;
Pstar = lyapunov_symm ( T , R * Q * R ' , DynareOptions . qz_criterium , DynareOptions . lyapunov_complex_threshold ) ;
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
Pinf = [ ] ;
otherwise
error ( ' DsgeLikelihood:: Unknown initialization approach for the Kalman filter!' )
2011-02-04 17:17:48 +01:00
end
2011-10-03 12:19:41 +02:00
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if ( ( kalman_algo == 1 ) || ( kalman_algo == 3 ) ) % Multivariate Kalman Filter
2011-02-04 17:17:48 +01:00
if no_missing_data_flag
2011-10-03 12:19:41 +02:00
[ LIK , lik ] = kalman_filter ( Y , diffuse_periods + 1 , size ( Y , 2 ) , ...
a , Pstar , ...
kalman_tol , riccati_tol , ...
DynareOptions . presample , ...
T , Q , R , H , Z , mm , pp , rr , Zflag , diffuse_periods ) ;
2011-02-04 17:17:48 +01:00
else
2011-10-03 12:19:41 +02:00
[ 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 ) , ...
a , Pstar , ...
kalman_tol , DynareOptions . riccati_tol , ...
DynareOptions . presample , ...
T , Q , R , H , Z , mm , pp , rr , Zflag , diffuse_periods ) ;
2011-02-04 17:17:48 +01:00
end
if isinf ( LIK )
2011-11-05 10:46:16 +01:00
if kalman_algo == 1
kalman_algo = 2 ;
else
kalman_algo = 4 ;
end
2011-10-03 12:19:41 +02:00
singularity_flag = 1 ;
else
if DynareOptions . lik_init == 3
LIK = LIK + dLIK ;
lik = [ dlik ; lik ] ;
end
2011-02-04 17:17:48 +01:00
end
end
2011-10-03 12:19:41 +02:00
2012-01-22 18:37:29 +01:00
if ( kalman_algo == 2 ) || ( kalman_algo == 4 )
2011-11-05 10:46:16 +01:00
% Univariate Kalman Filter
% resetting measurement error covariance matrix when necessary %
2012-01-22 18:37:29 +01:00
if ~ correlated_errors_have_been_checked
2011-11-05 10:46:16 +01:00
if isequal ( H , 0 )
H = zeros ( nobs , 1 ) ;
2011-10-03 12:19:41 +02:00
mmm = mm ;
2011-02-04 17:17:48 +01:00
else
2011-11-05 10:46:16 +01:00
if all ( all ( abs ( H - diag ( diag ( H ) ) ) < 1e-14 ) ) % ie, the covariance matrix is diagonal...
H = diag ( H ) ;
mmm = mm ;
else
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 ) ) ;
H = zeros ( nobs , 1 ) ;
mmm = mm + pp ;
end
2011-02-04 17:17:48 +01:00
end
end
2011-11-05 10:46:16 +01:00
2011-10-03 12:19:41 +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 ) , ...
a , Pstar , ...
DynareOptions . kalman_tol , ...
DynareOptions . riccati_tol , ...
DynareOptions . presample , ...
T , Q , R , H , Z , mmm , pp , rr , diffuse_periods ) ;
if DynareOptions . lik_init == 3
LIK = LIK + dLIK ;
lik = [ dlik ; lik ] ;
2011-02-04 17:17:48 +01:00
end
end
2011-10-03 12:19:41 +02:00
if isnan ( LIK )
info = 45 ;
exit_flag = 0 ;
return
end
if imag ( LIK ) ~= 0
likelihood = penalty ;
2011-02-04 17:17:48 +01:00
else
likelihood = LIK ;
end
2011-10-03 12:19:41 +02:00
2011-02-04 17:17:48 +01:00
% ------------------------------------------------------------------------------
2011-10-03 12:19:41 +02:00
% 5. Adds prior if necessary
2011-02-04 17:17:48 +01:00
% ------------------------------------------------------------------------------
2011-10-03 12:19:41 +02:00
lnprior = priordens ( xparam1 , BayesInfo . pshape , BayesInfo . p6 , BayesInfo . p7 , BayesInfo . p3 , BayesInfo . p4 ) ;
2011-02-04 17:17:48 +01:00
fval = ( likelihood - lnprior ) ;
2011-10-03 12:19:41 +02:00
% Update DynareOptions.kalman_algo.
DynareOptions . kalman_algo = kalman_algo ;
% Update the penalty.
penalty = fval ;
% Add the prior density at the top of the vector for the density of each observation.
2011-02-14 11:48:29 +01:00
lik = lik ( start : end , : ) ;
2011-10-03 12:19:41 +02:00
llik = [ - lnprior ; lik ( : ) ] ;