stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

adding functions and fixing bugs for diffuse filter 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1924 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2008-07-02 12:53:24 +00:00
parent 5ba344d75c
commit fe54c7ca07
5 changed files with 292 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
matlab/DiffuseLikelihoodH1_Z.m Normal file
View File

@ -0,0 +1,119 @@
function [LIK, lik] = DiffuseLikelihood1_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
% function [LIK, lik] = DiffuseLikelihood1_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).
%
% part of DYNARE, copyright Dynare Team (2004-2008)
% Gnu Public License.
% 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,1);
LIK = Inf;
lik(smpl+1) = smpl*pp*log(2*pi);
notsteady = 1;
crit = options_.kalman_tol;
reste = 0;
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
reste = smpl-t;
while t < smpl
t = t+1;
v = Y(:,t)-Z*a;
a = T*(a+K*v);
lik(t) = v'*iF*v;
end
lik(t) = lik(t) + reste*log(dF);
LIK = .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl);% Minus the log-likelihood.

169
matlab/DiffuseLikelihoodH3_Z.m Normal file
View File

@ -0,0 +1,169 @@
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).
%
% part of DYNARE, copyright Dynare Team (2004-2008)
% Gnu Public License.
% 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,1);
lik(smpl+1) = smpl*pp*log(2*pi); %% the constant of minus two times the log-likelihood
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
LIK = .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl);

6
matlab/DsgeLikelihood.m
View File

@ -161,9 +161,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
Pstar = 10*eye(np);
Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter
if kalman_algo == 1
kalman_algo == 3
end
kalman_algo = 3;
[QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
[QT,ST] = ordschur(QT,ST,e1);
@ -280,7 +278,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
end
elseif kalman_algo == 2
LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start);
elseif options_.kalman_algo == 3
elseif kalman_algo == 3
data1 = data - trend;
LIK = DiffuseLikelihood1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
if isinf(LIK)

2
matlab/DsgeSmoother.m
View File

@ -196,7 +196,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
nobs,np,smpl);
end
end
elseif options_.kalman_algo == 4
elseif kalman_algo == 4
data1 = Y - trend;
if ~estim_params.ncn
[alphahat,epsilonhat,etahat,ahat,aK] = ...

2
matlab/initial_estimation_checks.m
View File

@ -55,7 +55,7 @@ function initial_estimation_checks(xparam1,gend,data)
% Check if the steady state obtained from the _steadystate file is a
% steady state.
check1 = 0;
if isfield(options_,'unit_root_vars')
if isfield(options_,'unit_root_vars') & options_.diffuse_filter == 0
if isempty(options_.unit_root_vars)
check1 = max(abs(feval([M_.fname '_static'],...
oo_.steady_state,...