diff --git a/matlab/DiffuseLikelihood1.m b/matlab/DiffuseLikelihood1.m deleted file mode 100644 index f39132ab4d3932073a867d291c3bcc52f81d79e6..0000000000000000000000000000000000000000 --- 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 2edfc729eae9a44d427ba65a0120b236d9a88bcd..0000000000000000000000000000000000000000 --- 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 65b7591cb575877b4d6c04d4ba09996e75a3204d..0000000000000000000000000000000000000000 --- 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 f75377d300ef4a3bb98e41c390baaaf9342bbb14..0000000000000000000000000000000000000000 --- 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 88c4e619c68ae7bbbf7373e983de14320ca373db..0000000000000000000000000000000000000000 --- 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 6765a5be63e57e0e69b077199ce1e7b750d37836..0000000000000000000000000000000000000000 --- 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 f549b164da2ce142769d25ee02c2e95776d8f82b..0000000000000000000000000000000000000000 --- 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 c802cbf2f6e86e83a9d1b0f5789c137ff3f84a53..0000000000000000000000000000000000000000 --- 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 5d5a1c502b1cf3a95cb1445b283827633b34308a..0000000000000000000000000000000000000000 --- 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.