Removed unused routines for (diffuse) kalman filter evaluations.
Stéphane Adjemian (Charybdis)
parent
8ebc36df00
commit
37f14e9bc9
|
|
@ -1,130 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,Y,trend,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,Y,trend,start)
|
||||
% Computes the diffuse likelihood without measurement error, in the case of a non-singular var-cov matrix
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% trend
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2004-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
mf = bayestopt_.mf;
|
||||
smpl = size(Y,2);
|
||||
mm = size(T,2);
|
||||
pp = size(Y,1);
|
||||
a = zeros(mm,1);
|
||||
dF = 1;
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
LIK = Inf;
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
while rank(Pinf,crit) & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-a(mf)-trend(:,t);
|
||||
Finf = Pinf(mf,mf);
|
||||
if rcond(Finf) < crit
|
||||
if ~all(abs(Finf(:)) < crit)
|
||||
return
|
||||
else
|
||||
iFstar = inv(Pstar(mf,mf));
|
||||
dFstar = det(Pstar(mf,mf));
|
||||
Kstar = Pstar(:,mf)*iFstar;
|
||||
lik(t) = log(dFstar) + transpose(v)*iFstar*v;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
Pstar = T*(Pstar-Pstar(:,mf)*transpose(Kstar))*transpose(T)+QQ;
|
||||
a = T*(a+Kstar*v);
|
||||
end
|
||||
else
|
||||
lik(t) = log(det(Finf));
|
||||
iFinf = inv(Finf);
|
||||
Kinf = Pinf(:,mf)*iFinf; %% premultiplication by the transition matrix T is removed (stephane)
|
||||
Fstar = Pstar(mf,mf);
|
||||
Kstar = (Pstar(:,mf)-Kinf*Fstar)*iFinf; %% premultiplication by the transition matrix T is removed (stephane)
|
||||
Pstar = T*(Pstar-Pstar(:,mf)*transpose(Kinf)-Pinf(:,mf)*transpose(Kstar))*transpose(T)+QQ;
|
||||
Pinf = T*(Pinf-Pinf(:,mf)*transpose(Kinf))*transpose(T);
|
||||
a = T*(a+Kinf*v);
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
F_singular = 1;
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-a(mf)-trend(:,t);
|
||||
F = Pstar(mf,mf);
|
||||
oldPstar = Pstar;
|
||||
dF = det(F);
|
||||
if rcond(F) < crit
|
||||
if ~all(abs(F(:))<crit)
|
||||
return
|
||||
else
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T)+QQ;
|
||||
end
|
||||
else
|
||||
F_singular = 0;
|
||||
iF = inv(F);
|
||||
lik(t) = log(dF)+transpose(v)*iF*v;
|
||||
K = Pstar(:,mf)*iF; %% premultiplication by the transition matrix T is removed (stephane)
|
||||
a = T*(a+K*v); %% --> factorization of the transition matrix...
|
||||
Pstar = T*(Pstar-K*Pstar(mf,:))*transpose(T)+QQ; %% ... idem (stephane)
|
||||
end
|
||||
notsteady = ~(max(max(abs(Pstar-oldPstar)))<crit);
|
||||
end
|
||||
if F_singular == 1
|
||||
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)-trend(:,t);
|
||||
a = T*(a+K*v);
|
||||
lik(t) = transpose(v)*iF*v;
|
||||
end
|
||||
lik(t0:smpl) = lik(t0:smpl) + log(dF);
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,130 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihood1_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihood1_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
|
||||
% Computes the diffuse likelihood without measurement error, in the case of a non-singular var-cov matrix
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% Z: pp,mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2004-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
smpl = size(Y,2);
|
||||
mm = size(T,2);
|
||||
pp = size(Y,1);
|
||||
a = zeros(mm,1);
|
||||
dF = 1;
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
LIK = Inf;
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
while rank(Pinf,crit) & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-Z*a;
|
||||
Finf = Z*Pinf*Z';
|
||||
if rcond(Finf) < crit
|
||||
if ~all(abs(Finf(:)) < crit)
|
||||
return
|
||||
else
|
||||
Fstar = Z*Pstar*Z';
|
||||
iFstar = inv(Fstar);
|
||||
dFstar = det(Fstar);
|
||||
Kstar = Pstar*Z'*iFstar;
|
||||
lik(t) = log(dFstar) + v'*iFstar*v;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ;
|
||||
a = T*(a+Kstar*v);
|
||||
end
|
||||
else
|
||||
lik(t) = log(det(Finf));
|
||||
iFinf = inv(Finf);
|
||||
Kinf = Pinf*Z'*iFinf;
|
||||
Fstar = Z*Pstar*Z';
|
||||
Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf;
|
||||
Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ;
|
||||
Pinf = T*(Pinf-Pinf*Z'*Kinf')*T';
|
||||
a = T*(a+Kinf*v);
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
F_singular = 1;
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-Z*a;
|
||||
F = Z*Pstar*Z';
|
||||
oldPstar = Pstar;
|
||||
dF = det(F);
|
||||
if rcond(F) < crit
|
||||
if ~all(abs(F(:))<crit)
|
||||
return
|
||||
else
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*T'+QQ;
|
||||
end
|
||||
else
|
||||
F_singular = 0;
|
||||
iF = inv(F);
|
||||
lik(t) = log(dF)+v'*iF*v;
|
||||
K = Pstar*Z'*iF;
|
||||
a = T*(a+K*v);
|
||||
Pstar = T*(Pstar-K*Z*Pstar)*T'+QQ;
|
||||
end
|
||||
notsteady = ~(max(max(abs(Pstar-oldPstar)))<crit);
|
||||
end
|
||||
if F_singular == 1
|
||||
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)-Z*a;
|
||||
a = T*(a+K*v);
|
||||
lik(t) = v'*iF*v;
|
||||
end
|
||||
lik(t0:smpl) = lik(t0:smpl) + log(dF);
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,195 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,Y,trend,start)%//Z,T,R,Q,Pinf,Pstar,Y)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,Y,trend,start)
|
||||
% Computes the diffuse likelihood without measurement error, in the case of
|
||||
% a singular var-cov matrix.
|
||||
% Univariate treatment of multivariate time series.
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% trend
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2004-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output [October 2005]
|
||||
% changes by M. Ratto [April 2005]
|
||||
% introduced new options options_.diffuse_d for termination of DKF
|
||||
% new icc counter for Finf steps in DKF
|
||||
% new termination for DKF
|
||||
% likelihood terms for Fstar must be cumulated in DKF also when Pinf is non
|
||||
% zero.
|
||||
% [4/5/2005] correctyed bug in the modified verson of Ratto for rank of Pinf
|
||||
% introduced a specific crit1 for the DKF termination
|
||||
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
mf = bayestopt_.mf;
|
||||
pp = size(Y,1);
|
||||
mm = size(T,1);
|
||||
smpl = size(Y,2);
|
||||
a = zeros(mm,1);
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
crit1 = 1.e-6;
|
||||
newRank = rank(Pinf,crit1);
|
||||
icc=0;
|
||||
while newRank & t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t)-a(mf(i))-trend(i,t);
|
||||
Fstar = Pstar(mf(i),mf(i));
|
||||
Finf = Pinf(mf(i),mf(i));
|
||||
Kstar = Pstar(:,mf(i));
|
||||
if Finf > crit & newRank, %added newRank criterion
|
||||
icc=icc+1;
|
||||
Kinf = Pinf(:,mf(i));
|
||||
a = a + Kinf*v(i)/Finf;
|
||||
Pstar = Pstar + Kinf*transpose(Kinf)*Fstar/(Finf*Finf) - ...
|
||||
(Kstar*transpose(Kinf)+Kinf*transpose(Kstar))/Finf;
|
||||
Pinf = Pinf - Kinf*transpose(Kinf)/Finf;
|
||||
lik(t) = lik(t) + log(Finf);
|
||||
% start new termination criterion for DKF
|
||||
if ~isempty(options_.diffuse_d),
|
||||
newRank = (icc<options_.diffuse_d);
|
||||
%if newRank & any(diag(Pinf(mf,mf))>crit)==0; % M. Ratto this line is BUGGY
|
||||
if newRank & (any(diag(Pinf(mf,mf))>crit)==0 & rank(Pinf,crit1)==0);
|
||||
options_.diffuse_d = icc;
|
||||
newRank=0;
|
||||
disp('WARNING: Change in OPTIONS_.DIFFUSE_D in univariate DKF')
|
||||
disp(['new OPTIONS_.DIFFUSE_D = ',int2str(icc)])
|
||||
disp('You may have to reset the optimisation')
|
||||
end
|
||||
else
|
||||
%newRank = any(diag(Pinf(mf,mf))>crit); % M. Ratto this line is BUGGY
|
||||
newRank = (any(diag(Pinf(mf,mf))>crit) | rank(Pinf,crit1));
|
||||
if newRank==0,
|
||||
P0= T*Pinf*transpose(T);
|
||||
%newRank = any(diag(P0(mf,mf))>crit); % M. Ratto this line is BUGGY
|
||||
newRank = (any(diag(P0(mf,mf))>crit) | rank(P0,crit1)); % M. Ratto 11/10/2005
|
||||
if newRank==0,
|
||||
options_.diffuse_d = icc;
|
||||
end
|
||||
end
|
||||
end,
|
||||
% end new termination and checks for DKF
|
||||
elseif Fstar > crit
|
||||
%% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
|
||||
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
|
||||
%% rank(Pinf)=0. [stéphane,11-03-2004].
|
||||
%if rank(Pinf,crit) == 0
|
||||
% the likelihood terms should alwasy be cumulated, not only
|
||||
% when Pinf=0, otherwise the lik would depend on the ordering
|
||||
% of observed variables
|
||||
% presample options can be used to ignore initial time points
|
||||
lik(t) = lik(t) + log(Fstar) + v(i)*v(i)/Fstar;
|
||||
%end
|
||||
a = a + Kstar*v(i)/Fstar;
|
||||
Pstar = Pstar - Kstar*transpose(Kstar)/Fstar;
|
||||
else
|
||||
%disp(['zero F term in DKF for observed ',int2str(i),' ',num2str(Fstar)])
|
||||
end
|
||||
end
|
||||
% if all(abs(Pinf(:))<crit),
|
||||
% oldRank = 0;
|
||||
% else
|
||||
% oldRank = rank(Pinf,crit);
|
||||
% end
|
||||
if newRank,
|
||||
oldRank = rank(Pinf,crit1);
|
||||
else
|
||||
oldRank = 0;
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T)+QQ;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
% if all(abs(Pinf(:))<crit),
|
||||
% newRank = 0;
|
||||
% else
|
||||
% newRank = rank(Pinf,crit);
|
||||
% end
|
||||
if newRank,
|
||||
newRank = rank(Pinf,crit1); % new crit1 is used
|
||||
end
|
||||
if oldRank ~= newRank
|
||||
disp('DiffuseLiklihood3 :: T does influence the rank of Pinf!')
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
oldP = Pstar;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t) - a(mf(i)) - trend(i,t);
|
||||
Fi = Pstar(mf(i),mf(i));
|
||||
if Fi > crit
|
||||
Ki = Pstar(:,mf(i));
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*transpose(Ki)/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
else
|
||||
%disp(['zero F term for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T) + QQ;
|
||||
notsteady = ~(max(max(abs(Pstar-oldP)))<crit);
|
||||
end
|
||||
while t < smpl
|
||||
t = t+1;
|
||||
Pstar = oldP;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t) - a(mf(i)) - trend(i,t);
|
||||
Fi = Pstar(mf(i),mf(i));
|
||||
if Fi > crit
|
||||
Ki = Pstar(:,mf(i));
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*transpose(Ki)/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
else
|
||||
%disp(['zero F term for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
end
|
||||
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,180 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihood3_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihood3_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
|
||||
% Computes the diffuse likelihood without measurement error, in the case of
|
||||
% a singular var-cov matrix.
|
||||
% Univariate treatment of multivariate time series.
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% Z: pp*mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2004-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output [October 2005]
|
||||
% changes by M. Ratto [April 2005]
|
||||
% introduced new options options_.diffuse_d for termination of DKF
|
||||
% new icc counter for Finf steps in DKF
|
||||
% new termination for DKF
|
||||
% likelihood terms for Fstar must be cumulated in DKF also when Pinf is non
|
||||
% zero.
|
||||
% [4/5/2005] correctyed bug in the modified verson of Ratto for rank of Pinf
|
||||
% introduced a specific crit1 for the DKF termination
|
||||
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
pp = size(Y,1);
|
||||
mm = size(T,1);
|
||||
smpl = size(Y,2);
|
||||
a = zeros(mm,1);
|
||||
QQ = R*Q*R';
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
crit1 = 1.e-6;
|
||||
newRank = rank(Pinf,crit1);
|
||||
icc=0;
|
||||
while newRank & t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
Zi = Z(i,:);
|
||||
v(i) = Y(i,t)-Zi*a;
|
||||
Fstar = Zi*Pstar*Zi';
|
||||
Finf = Zi*Pinf*Zi';
|
||||
Kstar = Pstar*Zi';
|
||||
if Finf > crit & newRank
|
||||
icc=icc+1;
|
||||
Kinf = Pinf*Zi';
|
||||
a = a + Kinf*v(i)/Finf;
|
||||
Pstar = Pstar + Kinf*Kinf'*Fstar/(Finf*Finf) - ...
|
||||
(Kstar*Kinf'+Kinf*Kstar')/Finf;
|
||||
Pinf = Pinf - Kinf*Kinf'/Finf;
|
||||
lik(t) = lik(t) + log(Finf);
|
||||
if ~isempty(options_.diffuse_d),
|
||||
newRank = (icc<options_.diffuse_d);
|
||||
if newRank & (any(diag(Z*Pinf*Z')>crit)==0 & rank(Pinf,crit1)==0);
|
||||
options_.diffuse_d = icc;
|
||||
newRank=0;
|
||||
disp('WARNING: Change in OPTIONS_.DIFFUSE_D in univariate DKF')
|
||||
disp(['new OPTIONS_.DIFFUSE_D = ',int2str(icc)])
|
||||
disp('You may have to reset the optimisation')
|
||||
end
|
||||
else
|
||||
newRank = (any(diag(Z*Pinf*Z')>crit) | rank(Pinf,crit1));
|
||||
if newRank==0,
|
||||
P0= T*Pinf*T';
|
||||
newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1)); % M. Ratto 11/10/2005
|
||||
if newRank==0,
|
||||
options_.diffuse_d = icc;
|
||||
end
|
||||
end
|
||||
end,
|
||||
elseif Fstar > crit
|
||||
%% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
|
||||
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
|
||||
%% rank(Pinf)=0. [stéphane,11-03-2004].
|
||||
%if rank(Pinf,crit) == 0
|
||||
% the likelihood terms should alwasy be cumulated, not only
|
||||
% when Pinf=0, otherwise the lik would depend on the ordering
|
||||
% of observed variables
|
||||
% presample options can be used to ignore initial time points
|
||||
lik(t) = lik(t) + log(Fstar) + v(i)*v(i)/Fstar;
|
||||
a = a + Kstar*v(i)/Fstar;
|
||||
Pstar = Pstar - Kstar*(Kstar'/Fstar);
|
||||
else
|
||||
%disp(['zero F term in DKF for observed ',int2str(i),' ',num2str(Fstar)])
|
||||
end
|
||||
end
|
||||
if newRank,
|
||||
oldRank = rank(Pinf,crit1);
|
||||
else
|
||||
oldRank = 0;
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*T'+QQ;
|
||||
Pinf = T*Pinf*T';
|
||||
if newRank,
|
||||
newRank = rank(Pinf,crit1);
|
||||
end
|
||||
if oldRank ~= newRank
|
||||
disp('DiffuseLiklihood3 :: T does influence the rank of Pinf!')
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
oldP = Pstar;
|
||||
for i=1:pp
|
||||
Zi = Z(i,:);
|
||||
v(i) = Y(i,t) - Zi*a;
|
||||
Fi = Zi*Pstar*Zi';
|
||||
if Fi > crit
|
||||
Ki = Pstar*Zi';
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*(Ki'/Fi);
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
else
|
||||
%disp(['zero F term for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*T' + QQ;
|
||||
notsteady = ~(max(max(abs(Pstar-oldP)))<crit);
|
||||
end
|
||||
while t < smpl
|
||||
t = t+1;
|
||||
Pstar = oldP;
|
||||
for i=1:pp
|
||||
Zi = Z(i,:);
|
||||
v(i) = Y(i,t) - Zi*a;
|
||||
Fi = Zi*Pstar*Zi';
|
||||
if Fi > crit
|
||||
Ki = Pstar*Zi';
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*Ki'/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
else
|
||||
%disp(['zero F term for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,131 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,Y,trend,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,Y,trend,start)
|
||||
% Computes the diffuse likelihood (H=measurement error) in the case of a non-singular var-cov matrix
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% H: pp*pp matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% trend
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2005-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
mf = bayestopt_.mf;
|
||||
smpl = size(Y,2);
|
||||
mm = size(T,2);
|
||||
pp = size(Y,1);
|
||||
a = zeros(mm,1);
|
||||
dF = 1;
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
LIK = Inf;
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
while rank(Pinf,crit) & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-a(mf)-trend(:,t);
|
||||
Finf = Pinf(mf,mf);
|
||||
if rcond(Finf) < crit
|
||||
if ~all(abs(Finf(:))<crit)
|
||||
return
|
||||
else
|
||||
iFstar = inv(Pstar(mf,mf)+H);
|
||||
dFstar = det(Pstar(mf,mf)+H);
|
||||
Kstar = Pstar(:,mf)*iFstar;
|
||||
lik(t) = log(dFstar) + transpose(v)*iFstar*v;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
Pstar = T*(Pstar-Pstar(:,mf)*transpose(Kstar))*transpose(T)+QQ;
|
||||
a = T*(a+Kstar*v);
|
||||
end
|
||||
else
|
||||
lik(t) = log(det(Finf));
|
||||
iFinf = inv(Finf);
|
||||
Kinf = Pinf(:,mf)*iFinf; %% premultiplication by the transition matrix T is removed (stephane)
|
||||
Fstar = Pstar(mf,mf)+H;
|
||||
Kstar = (Pstar(:,mf)-Kinf*Fstar)*iFinf; %% premultiplication by the transition matrix T is removed (stephane)
|
||||
Pstar = T*(Pstar-Pstar(:,mf)*transpose(Kinf)-Pinf(:,mf)*transpose(Kstar))*transpose(T)+QQ;
|
||||
Pinf = T*(Pinf-Pinf(:,mf)*transpose(Kinf))*transpose(T);
|
||||
a = T*(a+Kinf*v);
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
F_singular = 1;
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-a(mf)-trend(:,t);
|
||||
F = Pstar(mf,mf)+H;
|
||||
oldPstar = Pstar;
|
||||
dF = det(F);
|
||||
if rcond(F) < crit
|
||||
if ~all(abs(F(:))<crit)
|
||||
return
|
||||
else
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T)+QQ;
|
||||
end
|
||||
else
|
||||
F_singular = 0;
|
||||
iF = inv(F);
|
||||
lik(t) = log(dF)+transpose(v)*iF*v;
|
||||
K = Pstar(:,mf)*iF; %% premultiplication by the transition matrix T is removed (stephane)
|
||||
a = T*(a+K*v); %% --> factorization of the transition matrix...
|
||||
Pstar = T*(Pstar-K*Pstar(mf,:))*transpose(T)+QQ; %% ... idem (stephane)
|
||||
end
|
||||
notsteady = ~(max(max(abs(Pstar-oldPstar)))<crit);
|
||||
end
|
||||
if F_singular == 1
|
||||
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)-trend(:,t);
|
||||
a = T*(a+K*v);
|
||||
lik(t) = transpose(v)*iF*v;
|
||||
end
|
||||
lik(t0:smpl) = lik(t0:smpl) + log(dF);
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,132 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihoodH1_Z(T,Z,R,Q,H,Pinf,Pstar,Y,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihoodH1_Z(T,Z,R,Q,H,Pinf,Pstar,Y,start)
|
||||
% Computes the diffuse likelihood (H=measurement error) in the case of a non-singular var-cov matrix
|
||||
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% Z: pp,mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% H: pp*pp matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2004-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
smpl = size(Y,2);
|
||||
mm = size(T,2);
|
||||
pp = size(Y,1);
|
||||
a = zeros(mm,1);
|
||||
dF = 1;
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
LIK = Inf;
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
while rank(Pinf,crit) & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-Z*a;
|
||||
Finf = Z*Pinf*Z';
|
||||
if rcond(Finf) < crit
|
||||
if ~all(abs(Finf(:)) < crit)
|
||||
return
|
||||
else
|
||||
Fstar = Z*Pstar*Z'+H;
|
||||
iFstar = inv(Fstar);
|
||||
dFstar = det(Fstar);
|
||||
Kstar = Pstar*Z'*iFstar;
|
||||
lik(t) = log(dFstar) + v'*iFstar*v;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ;
|
||||
a = T*(a+Kstar*v);
|
||||
end
|
||||
else
|
||||
lik(t) = log(det(Finf));
|
||||
iFinf = inv(Finf);
|
||||
Kinf = Pinf*Z'*iFinf;
|
||||
Fstar = Z*Pstar*Z'+H;
|
||||
Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf;
|
||||
Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ;
|
||||
Pinf = T*(Pinf-Pinf*Z'*Kinf')*T';
|
||||
a = T*(a+Kinf*v);
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
F_singular = 1;
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
v = Y(:,t)-Z*a;
|
||||
F = Z*Pstar*Z'+H;
|
||||
oldPstar = Pstar;
|
||||
dF = det(F);
|
||||
if rcond(F) < crit
|
||||
if ~all(abs(F(:))<crit)
|
||||
return
|
||||
else
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*T'+QQ;
|
||||
end
|
||||
else
|
||||
F_singular = 0;
|
||||
iF = inv(F);
|
||||
lik(t) = log(dF)+v'*iF*v;
|
||||
K = Pstar*Z'*iF;
|
||||
a = T*(a+K*v);
|
||||
Pstar = T*(Pstar-K*Z*Pstar)*T'+QQ;
|
||||
end
|
||||
notsteady = ~(max(max(abs(Pstar-oldPstar)))<crit);
|
||||
end
|
||||
if F_singular == 1
|
||||
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)-Z*a;
|
||||
a = T*(a+K*v);
|
||||
lik(t) = v'*iF*v;
|
||||
end
|
||||
lik(t0:smpl) = lik(t0:smpl) + log(dF);
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,178 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,Y,trend,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,Y,trend,start)
|
||||
% Computes the diffuse likelihood with measurement error, in the case of
|
||||
% a singular var-cov matrix.
|
||||
% Univariate treatment of multivariate time series.
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% H: pp*pp matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% trend
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2005-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output [October 2005]
|
||||
% changes by M. Ratto
|
||||
% introduced new global variable id_ for termination of DKF
|
||||
% introduced a persistent fmax, in order to keep track the max order of
|
||||
% magnitude of the 'zero' values in Pinf at DKF termination
|
||||
% new icc counter for Finf steps in DKF
|
||||
% new termination for DKF
|
||||
% likelihood terms for Fstar must be cumulated in DKF also when Pinf is non
|
||||
% zero. this bug is fixed.
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
mf = bayestopt_.mf;
|
||||
pp = size(Y,1);
|
||||
mm = size(T,1);
|
||||
smpl = size(Y,2);
|
||||
a = zeros(mm,1);
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
crit1 = 1.e-6;
|
||||
newRank = rank(Pinf,crit1);
|
||||
icc = 0;
|
||||
while newRank & t < smpl %% Matrix Finf is assumed to be zero
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t)-a(mf(i))-trend(i,t);
|
||||
Fstar = Pstar(mf(i),mf(i))+H(i,i);
|
||||
Finf = Pinf(mf(i),mf(i));
|
||||
Kstar = Pstar(:,mf(i));
|
||||
if Finf > crit & newRank
|
||||
icc = icc + 1;
|
||||
Kinf = Pinf(:,mf(i));
|
||||
a = a + Kinf*v(i)/Finf;
|
||||
Pstar = Pstar + Kinf*transpose(Kinf)*Fstar/(Finf*Finf) - ...
|
||||
(Kstar*transpose(Kinf)+Kinf*transpose(Kstar))/Finf;
|
||||
Pinf = Pinf - Kinf*transpose(Kinf)/Finf;
|
||||
lik(t) = lik(t) + log(Finf);
|
||||
% start new termination criterion for DKF
|
||||
if ~isempty(options_.diffuse_d),
|
||||
newRank = (icc<options_.diffuse_d);
|
||||
%if newRank & any(diag(Pinf(mf,mf))>crit)==0; % M. Ratto this line is BUGGY
|
||||
if newRank & (any(diag(Pinf(mf,mf))>crit)==0 & rank(Pinf,crit1)==0);
|
||||
options_.diffuse_d = icc;
|
||||
newRank=0;
|
||||
disp('WARNING: Change in OPTIONS_.DIFFUSE_D in univariate DKF')
|
||||
disp(['new OPTIONS_.DIFFUSE_D = ',int2str(icc)])
|
||||
disp('You may have to reset the optimisation')
|
||||
end
|
||||
else
|
||||
%newRank = any(diag(Pinf(mf,mf))>crit); % M. Ratto this line is BUGGY
|
||||
newRank = (any(diag(Pinf(mf,mf))>crit) | rank(Pinf,crit1));
|
||||
if newRank==0,
|
||||
P0= T*Pinf*transpose(T);
|
||||
%newRank = any(diag(P0(mf,mf))>crit); % M. Ratto this line is BUGGY
|
||||
newRank = (any(diag(P0(mf,mf))>crit) | rank(P0,crit1)); % M. Ratto 10 Oct 2005
|
||||
if newRank==0,
|
||||
options_.diffuse_d = icc;
|
||||
end
|
||||
end
|
||||
end,
|
||||
% end new termination and checks for DKF and fmax
|
||||
elseif Finf > crit
|
||||
%% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
|
||||
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
|
||||
%% rank(Pinf)=0. [stéphane,11-03-2004].
|
||||
%if rank(Pinf) == 0
|
||||
% the likelihood terms should alwasy be cumulated, not only
|
||||
% when Pinf=0, otherwise the lik would depend on the ordering
|
||||
% of observed variables
|
||||
lik(t) = lik(t) + log(Fstar) + v(i)*v(i)/Fstar;
|
||||
%end
|
||||
a = a + Kstar*v(i)/Fstar;
|
||||
Pstar = Pstar - Kstar*transpose(Kstar)/Fstar;
|
||||
else
|
||||
% disp(['zero F term in DKF for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
if newRank
|
||||
oldRank = rank(Pinf,crit1);
|
||||
else
|
||||
oldRank = 0;
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T)+QQ;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
if newRank
|
||||
newRank = rank(Pinf,crit1);
|
||||
end
|
||||
if oldRank ~= newRank
|
||||
disp('DiffuseLiklihoodH3 :: T does influence the rank of Pinf!')
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t) - a(mf(i)) - trend(i,t);
|
||||
Fi = Pstar(mf(i),mf(i))+H(i,i);
|
||||
if Fi > crit
|
||||
Ki = Pstar(:,mf(i));
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*transpose(Ki)/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
end
|
||||
end
|
||||
oldP = Pstar;
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T) + QQ;
|
||||
notsteady = ~(max(max(abs(Pstar-oldP)))<crit);
|
||||
end
|
||||
while t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t) - a(mf(i)) - trend(i,t);
|
||||
Fi = Pstar(mf(i),mf(i))+H(i,i);
|
||||
if Fi > crit
|
||||
Ki = Pstar(:,mf(i));
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*transpose(Ki)/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
|
@ -1,182 +0,0 @@
|
|||
function [LIK, lik] = DiffuseLikelihoodH3_Z(T,Z,R,Q,H,Pinf,Pstar,Y,start)
|
||||
|
||||
% function [LIK, lik] = DiffuseLikelihoodH3_A(T,R,Q,H,Pinf,Pstar,Y,start)
|
||||
% Computes the diffuse likelihood without measurement error, in the case of
|
||||
% a singular var-cov matrix.
|
||||
% Univariate treatment of multivariate time series.
|
||||
%
|
||||
% INPUTS
|
||||
% T: mm*mm matrix
|
||||
% Z: pp*mm matrix
|
||||
% R: mm*rr matrix
|
||||
% Q: rr*rr matrix
|
||||
% H: pp*pp matrix
|
||||
% Pinf: mm*mm diagonal matrix with with q ones and m-q zeros
|
||||
% Pstar: mm*mm variance-covariance matrix with stationary variables
|
||||
% Y: pp*1 vector
|
||||
% start: likelihood evaluation at 'start'
|
||||
%
|
||||
% OUTPUTS
|
||||
% LIK: likelihood
|
||||
% lik: density vector in each period
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% See "Filtering and Smoothing of State Vector for Diffuse State Space
|
||||
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
|
||||
% Analysis, vol. 24(1), pp. 85-98).
|
||||
|
||||
% Copyright (C) 2004-2008 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/>.
|
||||
|
||||
% M. Ratto added lik in output [October 2005]
|
||||
% changes by M. Ratto [April 2005]
|
||||
% introduced new options options_.diffuse_d for termination of DKF
|
||||
% new icc counter for Finf steps in DKF
|
||||
% new termination for DKF
|
||||
% likelihood terms for Fstar must be cumulated in DKF also when Pinf is non
|
||||
% zero.
|
||||
% [4/5/2005] correctyed bug in the modified verson of Ratto for rank of Pinf
|
||||
% introduced a specific crit1 for the DKF termination
|
||||
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
pp = size(Y,1);
|
||||
mm = size(T,1);
|
||||
smpl = size(Y,2);
|
||||
a = zeros(mm,1);
|
||||
QQ = R*Q*R';
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
crit1 = 1.e-6;
|
||||
newRank = rank(Pinf,crit1);
|
||||
icc=0;
|
||||
while newRank & t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
Zi = Z(i,:);
|
||||
v(i) = Y(i,t)-Zi*a;
|
||||
Fstar = Zi*Pstar*Zi'+H(i,i);
|
||||
Finf = Zi*Pinf*Zi';
|
||||
Kstar = Pstar*Zi';
|
||||
if Finf > crit & newRank
|
||||
icc=icc+1;
|
||||
Kinf = Pinf*Zi';
|
||||
a = a + Kinf*v(i)/Finf;
|
||||
Pstar = Pstar + Kinf*Kinf'*Fstar/(Finf*Finf) - ...
|
||||
(Kstar*Kinf'+Kinf*Kstar')/Finf;
|
||||
Pinf = Pinf - Kinf*Kinf'/Finf;
|
||||
lik(t) = lik(t) + log(Finf);
|
||||
if ~isempty(options_.diffuse_d),
|
||||
newRank = (icc<options_.diffuse_d);
|
||||
if newRank & (any(diag(Z*Pinf*Z')>crit)==0 & rank(Pinf,crit1)==0);
|
||||
options_.diffuse_d = icc;
|
||||
newRank=0;
|
||||
disp('WARNING: Change in OPTIONS_.DIFFUSE_D in univariate DKF')
|
||||
disp(['new OPTIONS_.DIFFUSE_D = ',int2str(icc)])
|
||||
disp('You may have to reset the optimisation')
|
||||
end
|
||||
else
|
||||
newRank = (any(diag(Z*Pinf*Z')>crit) | rank(Pinf,crit1));
|
||||
if newRank==0,
|
||||
P0= T*Pinf*T';
|
||||
newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1)); % M. Ratto 11/10/2005
|
||||
if newRank==0,
|
||||
options_.diffuse_d = icc;
|
||||
end
|
||||
end
|
||||
end,
|
||||
elseif Fstar > crit
|
||||
%% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
|
||||
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
|
||||
%% rank(Pinf)=0. [stéphane,11-03-2004].
|
||||
%if rank(Pinf,crit) == 0
|
||||
% the likelihood terms should alwasy be cumulated, not only
|
||||
% when Pinf=0, otherwise the lik would depend on the ordering
|
||||
% of observed variables
|
||||
% presample options can be used to ignore initial time points
|
||||
lik(t) = lik(t) + log(Fstar) + v(i)*v(i)/Fstar;
|
||||
a = a + Kstar*v(i)/Fstar;
|
||||
Pstar = Pstar - Kstar*Kstar'/Fstar;
|
||||
else
|
||||
%disp(['zero F term in DKF for observed ',int2str(i),' ',num2str(Fstar)])
|
||||
end
|
||||
end
|
||||
if newRank,
|
||||
oldRank = rank(Pinf,crit1);
|
||||
else
|
||||
oldRank = 0;
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*T'+QQ;
|
||||
Pinf = T*Pinf*T';
|
||||
if newRank,
|
||||
newRank = rank(Pinf,crit1);
|
||||
end
|
||||
if oldRank ~= newRank
|
||||
disp('DiffuseLiklihoodH3 :: T does influence the rank of Pinf!')
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
oldP = Pstar;
|
||||
for i=1:pp
|
||||
Zi = Z(i,:);
|
||||
v(i) = Y(i,t) - Zi*a;
|
||||
Fi = Zi*Pstar*Zi'+H(i,i);
|
||||
if Fi > crit
|
||||
Ki = Pstar*Zi';
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*Ki'/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
else
|
||||
%disp(['zero F term for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*T' + QQ;
|
||||
notsteady = ~(max(max(abs(Pstar-oldP)))<crit);
|
||||
end
|
||||
while t < smpl
|
||||
t = t+1;
|
||||
Pstar = oldP;
|
||||
for i=1:pp
|
||||
Zi = Z(i,:);
|
||||
v(i) = Y(i,t) - Zi*a;
|
||||
Fi = Zi*Pstar*Zi'+H(i,i);
|
||||
if Fi > crit
|
||||
Ki = Pstar*Zi';
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*Ki'/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
else
|
||||
%disp(['zero F term for observed ',int2str(i),' ',num2str(Fi)])
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
|
||||
|
|
@ -1,115 +0,0 @@
|
|||
function [LIK lik] = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,Y,trend,start)
|
||||
% Same as DiffuseLikelihoodH3 but allows correlation between the measurement
|
||||
% errors (this is not a problem with the multivariate approach).
|
||||
|
||||
% Copyright (C) 2004 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/>.
|
||||
|
||||
global bayestopt_ options_
|
||||
|
||||
mf = bayestopt_.mf;
|
||||
pp = size(Y,1);
|
||||
mm = size(T,1);
|
||||
rr = size(Q,1);
|
||||
smpl = size(Y,2);
|
||||
T = cat(1,cat(2,T,zeros(mm,pp)),zeros(pp,mm+pp));
|
||||
R = cat(1,cat(2,R,zeros(mm,pp)),cat(2,zeros(pp,rr),eye(pp)));
|
||||
Q = cat(1,cat(2,Q,zeros(rr,pp)),cat(2,zeros(pp,rr),H));
|
||||
if size(Pinf,1) % Otherwise Pinf = 0 (no unit root)
|
||||
Pinf = cat(1,cat(2,Pinf,zeros(mm,pp)),zeros(pp,mm+pp));
|
||||
end
|
||||
Pstar = cat(1,cat(2,Pstar,zeros(mm,pp)),cat(2,zeros(pp,mm),H));
|
||||
a = zeros(mm+pp,1);
|
||||
QQ = R*Q*transpose(R);
|
||||
t = 0;
|
||||
lik = zeros(smpl,1);
|
||||
notsteady = 1;
|
||||
crit = options_.kalman_tol;
|
||||
newRank = rank(Pinf,crit);
|
||||
|
||||
while rank(Pinf,crit) & t < smpl %% Matrix Finf is assumed to be zero
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t)-a(mf(i))-a(mm+i)-trend(i,t);
|
||||
Fstar = Pstar(mf(i),mf(i))+Pstar(mm+i,mm+i);
|
||||
Finf = Pinf(mf(i),mf(i));
|
||||
Kstar = Pstar(:,mf(i))+Pstar(:,mm+i);
|
||||
if Finf > crit
|
||||
Kinf = Pinf(:,mf(i));
|
||||
a = a + Kinf*v(i)/Finf;
|
||||
Pstar = Pstar + Kinf*transpose(Kinf)*Fstar/(Finf*Finf) - ...
|
||||
(Kstar*transpose(Kinf)+Kinf*transpose(Kstar))/Finf;
|
||||
Pinf = Pinf - Kinf*transpose(Kinf)/Finf;
|
||||
lik(t) = lik(t) + log(Finf);
|
||||
else %% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
|
||||
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
|
||||
%% rank(Pinf)=0. [stéphane,11-03-2004].
|
||||
if rank(Pinf) == 0
|
||||
lik(t) = lik(t) + log(Fstar) + v(i)*v(i)/Fstar;
|
||||
end
|
||||
a = a + Kstar*v(i)/Fstar;
|
||||
Pstar = Pstar - Kstar*transpose(Kstar)/Fstar;
|
||||
end
|
||||
oldRank = rank(Pinf,crit);
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T)+QQ;
|
||||
Pinf = T*Pinf*transpose(T);
|
||||
newRank = rank(Pinf,crit);
|
||||
if oldRank ~= newRank
|
||||
disp('DiffuseLiklihoodH3 :: T does influence the rank of Pinf!')
|
||||
end
|
||||
end
|
||||
end
|
||||
if t == smpl
|
||||
error(['There isn''t enough information to estimate the initial' ...
|
||||
' conditions of the nonstationary variables']);
|
||||
end
|
||||
while notsteady & t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t) - a(mf(i)) - trend(i,t) -a(mm+i);
|
||||
Fi = Pstar(mf(i),mf(i))+Pstar(mm+i,mm+i);
|
||||
if Fi > crit
|
||||
Ki = Pstar(:,mf(i))+Pstar(:,mm+i);
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*transpose(Ki)/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
end
|
||||
end
|
||||
oldP = Pstar;
|
||||
a = T*a;
|
||||
Pstar = T*Pstar*transpose(T) + QQ;
|
||||
notsteady = ~(max(max(abs(Pstar-oldP)))<crit);
|
||||
end
|
||||
while t < smpl
|
||||
t = t+1;
|
||||
for i=1:pp
|
||||
v(i) = Y(i,t) - a(mf(i)) - trend(i,t) - a(mm+i);
|
||||
Fi = Pstar(mf(i),mf(i))+Pstar(mm+i,mm+i);
|
||||
if Fi > crit
|
||||
Ki = Pstar(:,mf(i))+Pstar(:,mm+i);
|
||||
a = a + Ki*v(i)/Fi;
|
||||
Pstar = Pstar - Ki*transpose(Ki)/Fi;
|
||||
lik(t) = lik(t) + log(Fi) + v(i)*v(i)/Fi;
|
||||
end
|
||||
end
|
||||
a = T*a;
|
||||
end
|
||||
% adding log-likelihhod constants
|
||||
lik = (lik + pp*log(2*pi))/2;
|
||||
|
||||
LIK = sum(lik(start:end)); % Minus the log-likelihood.
|
||||
Loading…
Reference in New Issue