Skip to content

GitLab

D
dynare
Commit ef1cf5f0 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files
Options

Merge remote-tracking branch 'marco/master'

parents a9839333 9d595340
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 848 additions and 165 deletions
matlab/AHessian.m
View file @ ef1cf5f0
function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
% function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
% function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
%
% computes the asymptotic hessian matrix of the log-likelihood function of
% a state space model (notation as in kalman_filter.m in DYNARE
......@@ -28,6 +28,7 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
k = size(DT,3); % number of structural parameters
smpl = size(Y,2); % Sample size.
pp = size(Y,1); % Maximum number of observed variables.
mm = size(T,2); % Number of state variables.
a = zeros(mm,1); % State vector.
Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
......@@ -36,6 +37,11 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
notsteady = 1; % Steady state flag.
F_singular = 1;
lik = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
if nargout > 1,
DLIK = zeros(k,1); % Initialization of the score.
end
AHess = zeros(k,k); % Initialization of the Hessian
Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k);
......@@ -59,17 +65,21 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
F_singular = 0;
iF = inv(F);
K = P(:,mf)*iF;
lik(t) = log(det(F))+transpose(v)*iF*v;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start && nargout > 1
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
vecDPmf = reshape(DP(mf,mf,:),[],k);
iPmf = inv(P(mf,mf));
% iPmf = inv(P(mf,mf));
if t>=start
AHess = AHess + Dv'*iPmf*Dv + .5*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
end
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
......@@ -92,13 +102,23 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
for ii = 1:k
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start && nargout >1
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
if t>=start
AHess = AHess + Dv'*iPmf*Dv;
AHess = AHess + Dv'*iF*Dv;
end
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
if nargout > 1
for ii = 1:k
% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
end
end
AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
lik(t0:smpl) = lik(t0:smpl) + log(det(F));
% for ii = 1:k;
% for jj = 1:ii
% H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
......@@ -107,6 +127,13 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
end
AHess = -AHess;
if nargout > 1,
DLIK = DLIK/2;
end
% adding log-likelihhod constants
lik = (lik + pp*log(2*pi))/2;
LIK = sum(lik(start:end)); % Minus the log-likelihood.
% end of main function
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
......
matlab/DsgeLikelihood.m
View file @ ef1cf5f0
......@@ -439,9 +439,10 @@ end
if analytic_derivation
no_DLIK = 0;
full_Hess = 0;
DLIK = [];
AHess = [];
if nargin<7 || isempty(derivatives_info)
if nargin<8 || isempty(derivatives_info)
[A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
if ~isempty(EstimatedParameters.var_exo)
indexo=EstimatedParameters.var_exo(:,1);
......@@ -453,12 +454,25 @@ if analytic_derivation
else
indparam=[];
end
[dum, DT, DOm, DYss] = getH(A,B,Model,DynareResults,0,indparam,indexo);
if full_Hess,
[dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
else
[dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
end
else
DT = derivatives_info.DT;
DOm = derivatives_info.DOm;
DYss = derivatives_info.DYss;
if isfield(derivatives_info,'no_DLIK')
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'),
no_DLIK = derivatives_info.no_DLIK;
end
clear('derivatives_info');
......@@ -471,6 +485,17 @@ if analytic_derivation
DH=zeros([size(H),length(xparam1)]);
DQ=zeros([size(Q),length(xparam1)]);
DP=zeros([size(T),length(xparam1)]);
if full_Hess,
for j=1:size(D2Yss,1),
tmp(j,:,:) = blkdiag(zeros(offset,offset), squeeze(D2Yss(j,:,:)));
end
D2Yss = tmp;
D2T = D2T(iv,iv,:,:);
D2Om = D2Om(iv,iv,:,:);
D2Yss = D2Yss(iv,:,:);
D2H=zeros([size(H),length(xparam1),length(xparam1)]);
D2P=zeros([size(T),length(xparam1),length(xparam1)]);
end
for i=1:EstimatedParameters.nvx
k =EstimatedParameters.var_exo(i,1);
DQ(k,k,i) = 2*sqrt(Q(k,k));
......@@ -478,11 +503,23 @@ if analytic_derivation
kk = find(abs(dum) < 1e-12);
dum(kk) = 0;
DP(:,:,i)=dum;
if full_Hess
for j=1:i,
dum = lyapunov_symm(T,D2Om(:,:,i,j),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
kk = (abs(dum) < 1e-12);
dum(kk) = 0;
D2P(:,:,i,j)=dum;
D2P(:,:,j,i)=dum;
end
end
end
offset = EstimatedParameters.nvx;
for i=1:EstimatedParameters.nvn
k = EstimatedParameters.var_endo(i,1);
DH(k,k,i+offset) = 2*sqrt(H(k,k));
if full_Hess
D2H(k,k,i+offset,i+offset) = 2;
end
end
offset = offset + EstimatedParameters.nvn;
for j=1:EstimatedParameters.np
......@@ -490,6 +527,21 @@ if analytic_derivation
kk = find(abs(dum) < 1e-12);
dum(kk) = 0;
DP(:,:,j+offset)=dum;
if full_Hess
DTj = DT(:,:,j+offset);
DPj = dum;
for i=1:j,
DTi = DT(:,:,i+offset);
DPi = DP(:,:,i+offset);
D2Tij = D2T(:,:,i,j);
D2Omij = D2Om(:,:,i,j);
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);
dum(abs(dum)<1.e-12) = 0;
D2P(:,:,i+offset,j+offset) = dum;
D2P(:,:,j+offset,i+offset) = dum;
end
end
end
end
......@@ -515,6 +567,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
if nargout==11
[AHess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
if full_Hess,
Hess = 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
end
end
else
......@@ -587,7 +643,19 @@ end
% ------------------------------------------------------------------------------
% 5. Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
if analytic_derivation
if full_Hess,
[lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
AHess = Hess + d2lnprior;
else
[lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
end
if no_DLIK==0
DLIK = DLIK - dlnprior';
end
else
lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
end
fval = (likelihood-lnprior);
% Update DynareOptions.kalman_algo.
......
matlab/DsgeLikelihood_hh.m
View file @ ef1cf5f0
function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Evaluates the posterior kernel of a dsge model.
%
% INPUTS
......@@ -11,6 +11,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,D
% no_more_missing_observations [integer]
% OUTPUTS
% fval : value of the posterior kernel at xparam1.
% llik : probabilities at each time point
% cost_flag : zero if the function returns a penalty, one otherwise.
% ys : steady state of original endogenous variables
% trend_coeff :
......@@ -52,11 +53,14 @@ if nargin==1
return
end
fval = [];
ys = [];
trend_coeff = [];
cost_flag = 1;
% Initialization of the returned variables and others...
fval = [];
ys = [];
trend_coeff = [];
cost_flag = 1;
llik=NaN;
info = 0;
singularity_flag = 0;
if DynareOptions.block == 1
error('DsgeLikelihood_hh:: This routine (called if mode_compute==5) is not compatible with the block option!')
......@@ -159,8 +163,6 @@ end
Model.Sigma_e = Q;
Model.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
......@@ -223,8 +225,31 @@ Y = DynareDataset.data-trend;
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
kalman_algo = DynareOptions.kalman_algo;
% 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
diffuse_periods = 0;
switch DynareOptions.lik_init
case 1% Standard initialization with the steady state of the state equation.
......@@ -271,20 +296,32 @@ switch DynareOptions.lik_init
if isinf(dLIK)
% Go to univariate diffuse filter if singularity problem.
kalman_algo = 4;
singularity_flag = 1;
end
end
if (kalman_algo==4)
% Univariate Diffuse Kalman Filter
if no_correlation_flag
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));
mmm = mm+pp;
if singularity_flag
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
% no need to test again for correlation elements
singularity_flag = 0;
end
[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), ...
......@@ -294,13 +331,13 @@ switch DynareOptions.lik_init
diffuse_periods = length(dlik);
end
case 4% Start from the solution of the Riccati equation.
if kalman_algo ~= 2
if kalman_algo ~= 2
kalman_algo = 1;
end
if isequal(H,0)
[err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,length(mf))));
[err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))));
else
[err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,length(mf))),H);
[err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))),H);
end
if err
disp(['DsgeLikelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
......@@ -316,8 +353,6 @@ end
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
singularity_flag = 0;
if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
if no_missing_data_flag
[LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
......@@ -333,6 +368,11 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods);
end
if isinf(LIK)
if kalman_algo == 1
kalman_algo = 2;
else
kalman_algo = 4;
end
singularity_flag = 1;
else
if DynareOptions.lik_init==3
......@@ -342,21 +382,30 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
end
if ( (singularity_flag) || (kalman_algo==2) || (kalman_algo==4) )% Univariate Kalman Filter
if ( singularity_flag || (kalman_algo==2) || (kalman_algo==4) )
% Univariate Kalman Filter
% resetting measurement error covariance matrix when necessary %
if singularity_flag
if no_correlation
if isequal(H,0)
H = zeros(nobs,1);
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));
mmm = mm+pp;
a = [a; zeros(pp,1)];
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
[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, ...
......
matlab/disp_identification.m
View file @ ef1cf5f0
......@@ -188,7 +188,7 @@ if any(idemoments.ino),
iweak = length(find(idemoments.jweak_pair(:,j)));
if iweak,
[jx,jy]=find(jmap_pair==j);
jstore=[jstore jx(1) jy(1)]';
jstore=[jstore' jx(1) jy(1)]';
if SampleSize > 1
disp([' [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str((iweak)/SampleSize*100),'% of MC runs!' ])
else
......
matlab/dynare_estimation_1.m
View file @ ef1cf5f0
......@@ -353,7 +353,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
end
if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
mode_check('objective_function',xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
end
if ~options_.mh_posterior_mode_estimation
......
matlab/evaluate_steady_state_file.m
View file @ ef1cf5f0
......@@ -41,6 +41,8 @@ function [ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,f
if steadystate_flag == 1
% old format
assignin('base','tmp_00_',params);
evalin('base','M_.params=tmp_00_; clear(''tmp_00_'')');
h_steadystate = str2func([fname '_steadystate']);
[ys,check] = h_steadystate(ys_init, exo_ss);
params1 = evalin('base','M_.params');
......
matlab/getH.m
View file @ ef1cf5f0
This diff is collapsed.
matlab/get_Hessian.m 0 → 100644
View file @ ef1cf5f0
function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
% function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
%
% computes the hessian matrix of the log-likelihood function of
% a state space model (notation as in kalman_filter.m in DYNARE
% Thanks to Nikolai Iskrev
%
% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
% dimension equal to the number of structural parameters
% NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last
% two dimensions equal to the number of structural parameters
% Copyright (C) 2011 Dynare Team
%
% 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/>.
k = size(DT,3); % number of structural parameters
smpl = size(Y,2); % Sample size.
pp = size(Y,1); % Maximum number of observed variables.
mm = size(T,2); % Number of state variables.
a = zeros(mm,1); % State vector.
Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index.
oldK = 0;
notsteady = 1; % Steady state flag.
F_singular = 1;
Hess = zeros(k,k); % Initialization of the Hessian
Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k);
D2a = zeros(mm,k,k); % State vector.
D2v = zeros(length(mf),k,k);
C = zeros(length(mf),mm);
for ii=1:length(mf); C(ii,mf(ii))=1;end % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
dC = zeros(length(mf),mm,k);
d2C = zeros(length(mf),mm,k,k);
s = zeros(pp,1); % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
ds = zeros(pp,1,k);
d2s = zeros(pp,1,k,k);
% for ii = 1:k
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end
while notsteady & t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+Om;
end
else
F_singular = 0;
iF = inv(F);
K = P(:,mf)*iF;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
[D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,mf,iF,K,DK);
tmp = (a+K*v);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
% dai = da(:,:,ii);
dKi = DK(:,:,ii);
diFi = -iF*DF(:,:,ii)*iF;
dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
for jj = 1:ii
dFj = DF(:,:,jj);
diFj = -iF*DF(:,:,jj)*iF;
dKj = DK(:,:,jj);
d2Kij = D2K(:,:,jj,ii);
d2Fij = D2F(:,:,jj,ii);
d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
d2vij = -D2Yss(mf,jj,ii) - D2a(mf,jj,ii);
d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
end
Da(:,ii) = DT(:,:,ii)*tmp + T*dtmpi;
end
% vecDPmf = reshape(DP(mf,mf,:),[],k);
% iPmf = inv(P(mf,mf));
if t>=start
Hess = Hess + Hesst;
end
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
D2P = D2P1;
end
notsteady = max(max(abs(K-oldK))) > riccati_tol;
oldK = K;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
if t < smpl
t0 = t+1;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
tmp = (a+K*v);
for ii = 1:k,
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
dKi = DK(:,:,ii);
diFi = -iF*DF(:,:,ii)*iF;
dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
for jj = 1:ii,
dFj = DF(:,:,jj);
diFj = -iF*DF(:,:,jj)*iF;
dKj = DK(:,:,jj);
d2Kij = D2K(:,:,jj,ii);
d2Fij = D2F(:,:,jj,ii);
d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
d2vij = -D2Yss(mf,jj,ii) - D2a(mf,jj,ii);
d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
end
Da(:,ii) = DT(:,:,ii)*tmp + T*dtmpi;
end
if t>=start
Hess = Hess + Hesst;
end
a = T*(a+K*v);
end
% Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
% for ii = 1:k;
% for jj = 1:ii
% H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
% end
% end
end
Hess = Hess + tril(Hess,-1)';
Hess = -Hess/2;
% end of main function
function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
% computes (i,j) term in the Hessian
Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
% end of getHesst_ij
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
k = size(DT,3);
tmp = P-K*P(mf,:);
for ii = 1:k
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
end
% end of computeDKalman
function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0);
% computes the second derivatives of the Kalman matrices
% note: A=T in main func.
k = size(dA,3);
tmp = P0-K0*P0(mf,:);
[ns,no] = size(K0);
% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);