From 37f14e9bc9aa295723fb3533b608d929b918f09d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
<stephane.adjemian@ens.fr>
Date: Fri, 12 Nov 2010 17:20:02 +0100
Subject: [PATCH] Removed unused routines for (diffuse) kalman filter
evaluations.
---
matlab/DiffuseLikelihood1.m | 130 ---------------------
matlab/DiffuseLikelihood1_Z.m | 130 ---------------------
matlab/DiffuseLikelihood3.m | 195 -------------------------------
matlab/DiffuseLikelihood3_Z.m | 180 ----------------------------
matlab/DiffuseLikelihoodH1.m | 131 ---------------------
matlab/DiffuseLikelihoodH1_Z.m | 132 ---------------------
matlab/DiffuseLikelihoodH3.m | 178 ----------------------------
matlab/DiffuseLikelihoodH3_Z.m | 182 -----------------------------
matlab/DiffuseLikelihoodH3corr.m | 115 ------------------
9 files changed, 1373 deletions(-)
delete mode 100644 matlab/DiffuseLikelihood1.m
delete mode 100644 matlab/DiffuseLikelihood1_Z.m
delete mode 100644 matlab/DiffuseLikelihood3.m
delete mode 100644 matlab/DiffuseLikelihood3_Z.m
delete mode 100644 matlab/DiffuseLikelihoodH1.m
delete mode 100644 matlab/DiffuseLikelihoodH1_Z.m
delete mode 100644 matlab/DiffuseLikelihoodH3.m
delete mode 100644 matlab/DiffuseLikelihoodH3_Z.m
delete mode 100644 matlab/DiffuseLikelihoodH3corr.m
diff --git a/matlab/DiffuseLikelihood1.m b/matlab/DiffuseLikelihood1.m
deleted file mode 100644
index f39132ab4..000000000
--- a/matlab/DiffuseLikelihood1.m
+++ /dev/null
@@ -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.
\ No newline at end of file
diff --git a/matlab/DiffuseLikelihood1_Z.m b/matlab/DiffuseLikelihood1_Z.m
deleted file mode 100644
index 2edfc729e..000000000
--- a/matlab/DiffuseLikelihood1_Z.m
+++ /dev/null
@@ -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.
\ No newline at end of file
diff --git a/matlab/DiffuseLikelihood3.m b/matlab/DiffuseLikelihood3.m
deleted file mode 100644
index 65b7591cb..000000000
--- a/matlab/DiffuseLikelihood3.m
+++ /dev/null
@@ -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.
diff --git a/matlab/DiffuseLikelihood3_Z.m b/matlab/DiffuseLikelihood3_Z.m
deleted file mode 100644
index f75377d30..000000000
--- a/matlab/DiffuseLikelihood3_Z.m
+++ /dev/null
@@ -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.
diff --git a/matlab/DiffuseLikelihoodH1.m b/matlab/DiffuseLikelihoodH1.m
deleted file mode 100644
index 88c4e619c..000000000
--- a/matlab/DiffuseLikelihoodH1.m
+++ /dev/null
@@ -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.
diff --git a/matlab/DiffuseLikelihoodH1_Z.m b/matlab/DiffuseLikelihoodH1_Z.m
deleted file mode 100644
index 6765a5be6..000000000
--- a/matlab/DiffuseLikelihoodH1_Z.m
+++ /dev/null
@@ -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.
\ No newline at end of file
diff --git a/matlab/DiffuseLikelihoodH3.m b/matlab/DiffuseLikelihoodH3.m
deleted file mode 100644
index f549b164d..000000000
--- a/matlab/DiffuseLikelihoodH3.m
+++ /dev/null
@@ -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.
diff --git a/matlab/DiffuseLikelihoodH3_Z.m b/matlab/DiffuseLikelihoodH3_Z.m
deleted file mode 100644
index c802cbf2f..000000000
--- a/matlab/DiffuseLikelihoodH3_Z.m
+++ /dev/null
@@ -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.
-
diff --git a/matlab/DiffuseLikelihoodH3corr.m b/matlab/DiffuseLikelihoodH3corr.m
deleted file mode 100644
index 5d5a1c502..000000000
--- a/matlab/DiffuseLikelihoodH3corr.m
+++ /dev/null
@@ -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.
--
GitLab