diff --git a/matlab/DiffuseKalmanSmoother1.m b/matlab/DiffuseKalmanSmoother1.m
deleted file mode 100644
index e947bbe854ad06276cfd646c8d98115c970d3cd5..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmoother1.m
+++ /dev/null
@@ -1,177 +0,0 @@
-function [alphahat,etahat,atilde, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-
-% function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% trend
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% mf: observed variables index in the state vector
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% modified by M. Ratto:
-% new output argument aK (1-step to k-step predictions)
-% new options_.nk: the max step ahed prediction in aK (default is 4)
-% new crit1 value for rank of Pinf
-% it is assured that P is symmetric
-
-
-global options_
-
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-aK = zeros(nk,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-Fstar = zeros(pp,pp,smpl);
-iFinf = zeros(pp,pp,smpl);
-K = zeros(mm,pp,smpl);
-L = zeros(mm,mm,smpl);
-Linf = zeros(mm,mm,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-crit = options_.kalman_tol;
-crit1 = 1.e-8;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl+1);
-
-Z = zeros(pp,mm);
-for i=1:pp;
- Z(i,mf(i)) = 1;
-end
-
-t = 0;
-while rank(Pinf(:,:,t+1),crit1) && t<smpl
- t = t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- if rcond(Pinf(mf,mf,t)) < crit
- return
- end
- iFinf(:,:,t) = inv(Pinf(mf,mf,t));
- PZI = Pinf(:,mf,t)*iFinf(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- Kinf(:,:,t) = T*PZI;
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Linf(:,:,t) = T - Kinf(:,:,t)*Z;
- Fstar(:,:,t) = Pstar(mf,mf,t);
- Kstar(:,:,t) = (T*Pstar(:,mf,t)-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf,t)*transpose(Kinf(:,:,t))-Kinf(:,:,t)*Pinf(mf,mf,t)*transpose(Kstar(:,:,t)) + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf,t)*transpose(Kinf(:,:,t));
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-iFinf = iFinf(:,:,1:d);
-Linf = Linf(:,:,1:d);
-Fstar = Fstar(:,:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- if rcond(P(mf,mf,t)) < crit
- return
- end
- iF(:,:,t) = inv(P(mf,mf,t));
- PZI = P(:,mf,t)*iF(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- K(:,:,t) = T*PZI;
- L(:,:,t) = T-K(:,:,t)*Z;
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-if t<smpl
- PZI_s = PZI;
- K_s = K(:,:,t);
- iF_s = iF(:,:,t);
- P_s = P(:,:,t+1);
- t_steady = t+1;
- P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
- iF = cat(3,iF(:,:,1:t),repmat(inv(P_s(mf,mf)),[1 1 smpl-t_steady+1]));
- L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t_steady+1]));
- K = cat(3,K(:,:,1:t),repmat(T*P_s(:,mf)*iF_s,[1 1 smpl-t_steady+1]));
-end
-while t<smpl
- t=t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- a(:,t+1) = T*a(:,t) + K_s*v(:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
-end
-t = smpl+1;
-while t>d+1
- t = t-1;
- r(:,t) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t+1);
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-if d
- r0 = zeros(mm,d+1);
- r0(:,d+1) = r(:,d+1);
- r1 = zeros(mm,d+1);
- for t = d:-1:1
- r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
- r1(:,t) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t+1)) + Linf(:,:,t)'*r1(:,t+1);
- alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
- etahat(:,t) = QRt*r0(:,t);
- end
-end
\ No newline at end of file
diff --git a/matlab/DiffuseKalmanSmoother1_Z.m b/matlab/DiffuseKalmanSmoother1_Z.m
deleted file mode 100644
index 5ed71db973a7fa1ffc99be9f50bf398ec9b9360b..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmoother1_Z.m
+++ /dev/null
@@ -1,239 +0,0 @@
-function [alphahat,etahat,atilde,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl)
-
-% function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-%
-% OUTPUTS
-% alphahat: smoothed variables (a_{t|T})
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t)
-% (meaningless for periods 1:d)
-% P: 3D array of one-step ahead forecast error variance
-% matrices
-% PK: 4D array of k-step ahead forecast error variance
-% matrices (meaningless for periods 1:d)
-% d: number of periods where filter remains in diffuse part
-% (should be equal to the order of integration of the model)
-% decomp: decomposition of the effect of shocks on filtered values
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% modified by M. Ratto:
-% new output argument aK (1-step to k-step predictions)
-% new options_.nk: the max step ahed prediction in aK (default is 4)
-% new crit1 value for rank of Pinf
-% it is assured that P is symmetric
-
-
-global options_
-
-d = 0;
-decomp = [];
-nk = options_.nk;
-kalman_tol = options_.kalman_tol;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-aK = zeros(nk,mm,smpl+nk);
-PK = zeros(nk,mm,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-Fstar = zeros(pp,pp,smpl);
-iFinf = zeros(pp,pp,smpl);
-K = zeros(mm,pp,smpl);
-L = zeros(mm,mm,smpl);
-Linf = zeros(mm,mm,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-crit = options_.kalman_tol;
-crit1 = 1.e-8;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl+1);
-
-t = 0;
-while rank(Pinf(:,:,t+1),kalman_tol) && (t<smpl)
- t = t+1;
- v(:,t)= Y(:,t) - Z*a(:,t);
- Finf = Z*Pinf(:,:,t)*Z' ;
- if rcond(Finf) < kalman_tol
- if ~all(abs(Finf(:)) < kalman_tol)
- % The univariate diffuse kalman filter should be used.
- return
- else
- Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
- if rcond(Fstar(:,:,t)) < kalman_tol
- if ~all(abs(Fstar(:,:,t))<kalman_tol)
- % The univariate diffuse kalman filter should be used.
- return
- else
- a(:,:,t+1) = T*a(:,:,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- end
- else
- iFstar = inv(Fstar(:,:,t));
- Kstar(:,:,t) = Pstar(:,:,t)*Z'*iFstar(:,:,t);
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*Z'*Kstar(:,:,t)')*T'+QQ;
- a(:,:,t+1) = T*(a(:,:,t)+Kstar(:,:,t)*v(:,t));
- end
- end
- else
- iFinf(:,:,t) = inv(Finf);
- PZI = Pinf(:,:,t)*Z'*iFinf(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- Kinf(:,:,t) = T*PZI;
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- % isn't a meaningless as long as we are in the diffuse part? MJ
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Linf(:,:,t) = T - Kinf(:,:,t)*Z;
- Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
- Kstar(:,:,t) = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
- % Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*Fstar(:,:,t)*Kstar(:,:,t)' + QQ;
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-T*Pinf(:,:,t)*Z'*Kstar(:,:,t)' + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
- end
-
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-iFinf = iFinf(:,:,1:d);
-Linf = Linf(:,:,1:d);
-Fstar = Fstar(:,:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- v(:,t) = Y(:,t) - Z*a(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- F = Z*P(:,:,t)*Z';
- if rcond(F) < crit
- return
- end
- iF(:,:,t) = inv(F);
- PZI = P(:,:,t)*Z'*iF(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- K(:,:,t) = T*PZI;
- L(:,:,t) = T-K(:,:,t)*Z;
- a(:,t+1) = T*atilde(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk,
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
- P(:,:,t+1) = T*P(:,:,t)*T'-T*P(:,:,t)*Z'*K(:,:,t)' + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-if t<smpl
- PZI_s = PZI;
- K_s = K(:,:,t);
- iF_s = iF(:,:,t);
- P_s = P(:,:,t+1);
- P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
- iF = cat(3,iF(:,:,1:t),repmat(iF_s,[1 1 smpl-t]));
- L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t]));
- K = cat(3,K(:,:,1:t),repmat(T*P_s*Z'*iF_s,[1 1 smpl-t]));
-end
-while t<smpl
- t=t+1;
- v(:,t) = Y(:,t) - Z*a(:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- a(:,t+1) = T*atilde(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
-end
-t = smpl+1;
-while t>d+1
- t = t-1;
- r(:,t) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t+1);
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-if d
- r0 = zeros(mm,d+1);
- r0(:,d+1) = r(:,d+1);
- r1 = zeros(mm,d+1);
- for t = d:-1:1
- r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
- r1(:,t) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t+1)) + Linf(:,:,t)'*r1(:,t+1);
- alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
- etahat(:,t) = QRt*r0(:,t);
- end
-end
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = Z'*iF(:,:,t)*v(:,t);
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- for h = 1:nk
- for j=1:rr
- eta=zeros(rr,1);
- eta(j) = eta_tm1t(j);
- decomp(h,:,j,t+h) = T^(h-1)*P(:,:,t)*Z'*ZRQinv*Z*R*eta;
- end
- end
- end
-end
diff --git a/matlab/DiffuseKalmanSmoother3.m b/matlab/DiffuseKalmanSmoother3.m
deleted file mode 100644
index f3917e1aa9eb6b232b7a368306bef77925f51af4..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmoother3.m
+++ /dev/null
@@ -1,297 +0,0 @@
-function [alphahat,etahat,a,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% function [alphahat,etahat,a1, aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% trend
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% mf: observed variables index in the state vector
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% a: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% Modified by M. Ratto
-% New output argument aK: 1-step to nk-stpe ahed predictions)
-% New input argument nk: max order of predictions in aK
-% New option options_.diffuse_d where the DKF stops (common with
-% diffuselikelihood3)
-% New icc variable to count number of iterations for Finf steps
-% Pstar % Pinf simmetric
-% New termination of DKF iterations based on options_.diffuse_d
-% Li also stored during DKF iterations !!
-% some bugs corrected in the DKF part of the smoother (Z matrix and
-% alphahat)
-
-global options_
-
-d=0;
-decomp=[];
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl);
-a1 = zeros(mm,smpl+1);
-aK = zeros(nk,mm,smpl+nk);
-PK = zeros(nk,mm,mm,smpl+nk);
-Fstar = zeros(pp,smpl);
-Finf = zeros(pp,smpl);
-Fi = zeros(pp,smpl);
-Ki = zeros(mm,pp,smpl);
-Li = zeros(mm,mm,pp,smpl);
-Linf = zeros(mm,mm,pp,smpl);
-L0 = zeros(mm,mm,pp,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-P1 = P;
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-Pstar1 = Pstar;
-Pinf1 = Pinf;
-crit = options_.kalman_tol;
-crit1 = 1.e-6;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl+1);
-
-Z = zeros(pp,mm);
-for i=1:pp;
- Z(i,mf(i)) = 1;
-end
-
-t = 0;
-icc=0;
-newRank = rank(Pinf(:,:,1),crit1);
-while newRank && t < smpl
- t = t+1;
- a(:,t) = a1(:,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
- Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1));
- Pstar1(:,:,t) = Pstar(:,:,t);
- Pinf1(:,:,t) = Pinf(:,:,t);
- for i=1:pp
- v(i,t) = Y(i,t)-a(mf(i),t)-trend(i,t);
- Fstar(i,t) = Pstar(mf(i),mf(i),t);
- Finf(i,t) = Pinf(mf(i),mf(i),t);
- Kstar(:,i,t) = Pstar(:,mf(i),t);
- if Finf(i,t) > crit && newRank, % original MJ: if Finf(i,t) > crit
- icc=icc+1;
- Kinf(:,i,t) = Pinf(:,mf(i),t);
- Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
- L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Z(i,:)/Finf(i,t);
- a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) + ...
- Kinf(:,i,t)*transpose(Kinf(:,i,t))*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
- (Kstar(:,i,t)*transpose(Kinf(:,i,t)) +...
- Kinf(:,i,t)*transpose(Kstar(:,i,t)))/Finf(i,t);
- Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*transpose(Kinf(:,i,t))/Finf(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
- Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1));
- P0=Pinf(:,:,t);
- elseif Fstar(i,t) > 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].
- Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother
- a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*transpose(Kstar(:,i,t))/Fstar(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- if newRank
- oldRank = rank(P0,crit);
- else
- oldRank = 0;
- end
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- P0 = Pinf(:,:,t+1);
- if newRank
- newRank = rank(P0,crit1);
- end
- if oldRank ~= newRank
- disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!')
- end
-end
-
-
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-Linf = Linf(:,:,:,1:d);
-L0 = L0(:,:,:,1:d);
-Fstar = Fstar(:,1:d);
-Finf = Finf(:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-Pstar1 = Pstar1(:,:,1:d);
-Pinf1 = Pinf1(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- a(:,t) = a1(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- P1(:,:,t) = P(:,:,t);
- for i=1:pp
- v(i,t) = Y(i,t) - a(mf(i),t) - trend(i,t);
- Fi(i,t) = P(mf(i),mf(i),t);
- Ki(:,i,t) = P(:,mf(i),t);
- if Fi(i,t) > crit
- Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
- a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
- P(:,:,t) = P(:,:,t) - Ki(:,i,t)*transpose(Ki(:,i,t))/Fi(i,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- end
- end
- a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk,
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
- P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
-Fi_s = Fi(:,t);
-Ki_s = Ki(:,:,t);
-L_s =Li(:,:,:,t);
-if t<smpl
- t_steady = t+1;
- P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
- Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
- Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
- Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
-end
-while t<smpl
- t=t+1;
- a(:,t) = a1(:,t);
- for i=1:pp
- v(i,t) = Y(i,t) - a(mf(i),t) - trend(i,t);
- if Fi_s(i) > crit
- a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
- end
- end
- a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
-end
-ri=zeros(mm,1);
-t = smpl+1;
-while t>d+1
- t = t-1;
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri;
- end
- end
- r(:,t) = ri;
- alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
- ri = T'*ri;
-end
-if d
- r0 = zeros(mm,d);
- r0(:,d) = ri;
- r1 = zeros(mm,d);
- for t = d:-1:1
- for i=pp:-1:1
- if Finf(i,t) > crit && ~(t==d && i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination
- %r1(:,t) = transpose(Z)*v(:,t)/Finf(i,t) + ... BUG HERE in transpose(Z)
- r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
- L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t);
- r0(:,t) = Linf(:,:,i,t)'*r0(:,t);
- elseif Fstar(i,t) > crit % step needed whe Finf == 0
- r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+Li(:,:,i,t)'*r0(:,t);
- end
- end
- alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
- r(:,t) = r0(:,t);
- etahat(:,t) = QRt*r(:,t);
- if t > 1
- r0(:,t-1) = T'*r0(:,t);
- r1(:,t-1) = T'*r1(:,t);
- end
- end
-end
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = zeros(mm,1);
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri_d;
- end
- end
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- for h = 1:nk
- for j=1:rr
- eta=zeros(rr,1);
- eta(j) = eta_tm1t(j);
- decomp(h,:,j,t+h) = Ttok*P1(:,:,t)*Z'*ZRQinv*Z*R*eta;
- end
- Ttok = T*Ttok;
- end
- end
-end
-
diff --git a/matlab/DiffuseKalmanSmoother3_Z.m b/matlab/DiffuseKalmanSmoother3_Z.m
deleted file mode 100644
index 6f42616c7543c6858586a6e3d6931347bded09dc..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmoother3_Z.m
+++ /dev/null
@@ -1,314 +0,0 @@
-function [alphahat,etahat,a,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl)
-% function [alphahat,etahat,a1,P,aK,PK,d,decomp_filt] = DiffuseKalmanSmoother3(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% a: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-% (meaningless for periods 1:d)
-% P: 3D array of one-step ahead forecast error variance
-% matrices
-% PK: 4D array of k-step ahead forecast error variance
-% matrices (meaningless for periods 1:d)
-% d: number of periods where filter remains in diffuse part
-% (should be equal to the order of integration of the model)
-% decomp: decomposition of the effect of shocks on filtered values
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% Modified by M. Ratto
-% New output argument aK: 1-step to nk-stpe ahed predictions)
-% New input argument nk: max order of predictions in aK
-% New option options_.diffuse_d where the DKF stops (common with
-% diffuselikelihood3)
-% New icc variable to count number of iterations for Finf steps
-% Pstar % Pinf simmetric
-% New termination of DKF iterations based on options_.diffuse_d
-% Li also stored during DKF iterations !!
-% some bugs corrected in the DKF part of the smoother (Z matrix and
-% alphahat)
-
-global options_
-
-d = 0;
-decomp = [];
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl);
-a1 = zeros(mm,smpl+1);
-aK = zeros(nk,mm,smpl+nk);
-
-if isempty(options_.diffuse_d),
- smpl_diff = 1;
-else
- smpl_diff=rank(Pinf1);
-end
-
-Fstar = zeros(pp,smpl_diff);
-Finf = zeros(pp,smpl_diff);
-Fi = zeros(pp,smpl);
-Ki = zeros(mm,pp,smpl);
-%Li = zeros(mm,mm,pp,smpl);
-%Linf = zeros(mm,mm,pp,smpl_diff);
-%L0 = zeros(mm,mm,pp,smpl_diff);
-Kstar = zeros(mm,pp,smpl_diff);
-P = zeros(mm,mm,smpl+1);
-P1 = P;
-PK = zeros(nk,mm,mm,smpl+nk);
-Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
-Pstar1 = Pstar;
-Pinf1 = Pinf;
-crit = options_.kalman_tol;
-crit1 = 1.e-6;
-steady = smpl;
-rr = size(Q,1); % number of structural shocks
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl);
-
-t = 0;
-icc=0;
-newRank = rank(Pinf(:,:,1),crit1);
-while newRank && t < smpl
- t = t+1;
- a(:,t) = a1(:,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
- Pstar1(:,:,t) = Pstar(:,:,t);
- Pinf1(:,:,t) = Pinf(:,:,t);
- for i=1:pp
- Zi = Z(i,:);
- v(i,t) = Y(i,t)-Zi*a(:,t);
- Fstar(i,t) = Zi*Pstar(:,:,t)*Zi';
- Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
- Kstar(:,i,t) = Pstar(:,:,t)*Zi';
- if Finf(i,t) > crit && newRank
- icc=icc+1;
- Kinf(:,i,t) = Pinf(:,:,t)*Zi';
- % Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
- % L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Zi/Finf(i,t);
- a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) + ...
- Kinf(:,i,t)*Kinf(:,i,t)'*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
- (Kstar(:,i,t)*Kinf(:,i,t)' +...
- Kinf(:,i,t)*Kstar(:,i,t)')/Finf(i,t);
- Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
- % new terminiation criteria by M. Ratto
- P0=Pinf(:,:,t);
- if ~isempty(options_.diffuse_d),
- newRank = (icc<options_.diffuse_d);
- if newRank && (any(diag(Z*P0*Z')>crit)==0 && rank(P0,crit1)==0);
- disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF')
- options_.diffuse_d = icc;
- newRank=0;
- end
- else
- newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1));
- if newRank==0,
- options_.diffuse_d = icc;
- end
- end,
- % end new terminiation criteria by M. Ratto
- elseif Fstar(i,t) > 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].
- % Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother
- a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
- P0=Pinf(:,:,t+1);
- if newRank,
- newRank = rank(P0,crit1);
- end
-end
-
-
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-%Linf = Linf(:,:,:,1:d);
-%L0 = L0(:,:,:,1:d);
-Fstar = Fstar(:,1:d);
-Finf = Finf(:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-Pstar1 = Pstar1(:,:,1:d);
-Pinf1 = Pinf1(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- a(:,t) = a1(:,t);
- P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
- P1(:,:,t) = P(:,:,t);
- for i=1:pp
- Zi = Z(i,:);
- v(i,t) = Y(i,t) - Zi*a(:,t);
- Fi(i,t) = Zi*P(:,:,t)*Zi';
- Ki(:,i,t) = P(:,:,t)*Zi';
- if Fi(i,t) > crit
- % Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
- a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
- P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t);
- P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
- end
- end
- a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk,
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
- P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
- % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
-% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
-% $$$ Fi_s = Fi(:,t);
-% $$$ Ki_s = Ki(:,:,t);
-% $$$ L_s =Li(:,:,:,t);
-% $$$ if t<smpl
-% $$$ P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
-% $$$ P1 = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
-% $$$ Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
-% $$$ Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
-% $$$ Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
-% $$$ end
-% $$$ while t<smpl
-% $$$ t=t+1;
-% $$$ a(:,t) = a1(:,t);
-% $$$ for i=1:pp
-% $$$ Zi = Z(i,:);
-% $$$ v(i,t) = Y(i,t) - Zi*a(:,t);
-% $$$ if Fi_s(i) > crit
-% $$$ a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
-% $$$ end
-% $$$ end
-% $$$ a1(:,t+1) = T*a(:,t);
-% $$$ Pf = P(:,:,t);
-% $$$ for jnk=1:nk,
-% $$$ Pf = T*Pf*T' + QQ;
-% $$$ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
-% $$$ PK(jnk,:,:,t+jnk) = Pf;
-% $$$ end
-% $$$ end
-ri=zeros(mm,1);
-t = smpl+1;
-while t > d+1
- t = t-1;
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri = Z(i,:)'/Fi(i,t)*v(i,t)+ri-Ki(:,i,t)'*ri/Fi(i,t)*Z(i,:)';
- end
- end
- r(:,t) = ri;
- alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
- ri = T'*ri;
-end
-if d
- r0 = zeros(mm,d);
- r0(:,d) = ri;
- r1 = zeros(mm,d);
- for t = d:-1:1
- for i=pp:-1:1
- % if Finf(i,t) > crit && ~(t==d && i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination
- if Finf(i,t) > crit
- r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
- (Kinf(:,i,t)'*Fstar(i,t)/Finf(i,t)-Kstar(:,i,t)')*r0(:,t)/Finf(i,t)*Z(i,:)' + ...
- r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)';
- r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)';
- elseif Fstar(i,t) > crit % step needed whe Finf == 0
- r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+r0(:,t)-(Kstar(:,i,t)'*r0(:,t))/Fstar(i,t)*Z(i,:)';
- end
- end
- alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
- r(:,t) = r0(:,t);
- etahat(:,t) = QRt*r(:,t);
- if t > 1
- r0(:,t-1) = T'*r0(:,t);
- r1(:,t-1) = T'*r1(:,t);
- end
- end
-end
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = zeros(mm,1);
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)';
- end
- end
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- AAA = P1(:,:,t)*Z'*ZRQinv*Z*R;
- for h = 1:nk
- BBB = Ttok*AAA;
- for j=1:rr
- decomp(h,:,j,t+h) = eta_tm1t(j)*BBB(:,j);
- end
- Ttok = T*Ttok;
- end
- end
-end
diff --git a/matlab/DiffuseKalmanSmootherH1.m b/matlab/DiffuseKalmanSmootherH1.m
deleted file mode 100644
index 738a5fa3b0ba87303cee6d894de67e79e332fbeb..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmootherH1.m
+++ /dev/null
@@ -1,179 +0,0 @@
-function [alphahat,epsilonhat,etahat,atilde, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% function [alphahat,epsilonhat,etahat,atilde, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% Computes the diffuse kalman smoother with 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 variance of measurement errors
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% trend
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% mf: observed variables index in the state vector
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% epsilonhat:smoothed measurement errors
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t)}
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% modified by M. Ratto:
-% new output argument aK (1-step to k-step predictions)
-% new options_.nk: the max step ahed prediction in aK (default is 4)
-% new crit1 value for rank of Pinf
-% it is assured that P is symmetric
-
-global options_
-
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-aK = zeros(nk,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-Fstar = zeros(pp,pp,smpl);
-iFinf = zeros(pp,pp,smpl);
-K = zeros(mm,pp,smpl);
-L = zeros(mm,mm,smpl);
-Linf = zeros(mm,mm,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-crit = options_.kalman_tol;
-crit1 = 1.e-8;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-epsilonhat = zeros(size(Y));
-r = zeros(mm,smpl+1);
-
-Z = zeros(pp,mm);
-for i=1:pp;
- Z(i,mf(i)) = 1;
-end
-
-t = 0;
-while rank(Pinf(:,:,t+1),crit1) && t<smpl
- t = t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- if rcond(Pinf(mf,mf,t)) < crit
- return
- end
- iFinf(:,:,t) = inv(Pinf(mf,mf,t));
- PZI = Pinf(:,mf,t)*iFinf(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- Kinf(:,:,t) = T*PZI;
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Linf(:,:,t) = T - Kinf(:,:,t)*Z;
- Fstar(:,:,t) = Pstar(mf,mf,t) + H;
- Kstar(:,:,t) = (T*Pstar(:,mf,t)-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf,t)*transpose(Kinf(:,:,t))-Kinf(:,:,t)*Pinf(mf,mf,t)*transpose(Kstar(:,:,t)) + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf,t)*transpose(Kinf(:,:,t));
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-iFinf = iFinf(:,:,1:d);
-Linf = Linf(:,:,1:d);
-Fstar = Fstar(:,:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- if rcond(P(mf,mf,t) + H) < crit
- return
- end
- iF(:,:,t) = inv(P(mf,mf,t) + H);
- PZI = P(:,mf,t)*iF(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- K(:,:,t) = T*PZI;
- L(:,:,t) = T-K(:,:,t)*Z;
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-if t<smpl
- PZI_s = PZI;
- K_s = K(:,:,t);
- iF_s = iF(:,:,t);
- P_s = P(:,:,t+1);
- t_steady = t+1;
- P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
- iF = cat(3,iF(:,:,1:t),repmat(inv(P_s(mf,mf)+H),[1 1 smpl-t_steady+1]));
- L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t_steady+1]));
- K = cat(3,K(:,:,1:t),repmat(T*P_s(:,mf)*iF_s,[1 1 smpl-t_steady+1]));
-end
-while t<smpl
- t=t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- atilde(:,t) = a(:,t) + PZI_s*v(:,t);
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
-end
-t = smpl+1;
-while t>d+1
- t = t-1;
- r(:,t) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t+1);
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-if d
- r0 = zeros(mm,d+1);
- r0(:,d+1) = r(:,d+1);
- r1 = zeros(mm,d+1);
- for t = d:-1:1
- r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
- r1(:,t) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t+1)) + Linf(:,:,t)'*r1(:,t+1);
- alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
- etahat(:,t) = QRt*r0(:,t);
- end
-end
-epsilonhat = Y-alphahat(mf,:)-trend;
diff --git a/matlab/DiffuseKalmanSmootherH1_Z.m b/matlab/DiffuseKalmanSmootherH1_Z.m
deleted file mode 100644
index 130e64d33a8055c95d66ae88e66b82cace1ca8d8..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmootherH1_Z.m
+++ /dev/null
@@ -1,241 +0,0 @@
-function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,d,decomp] = DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
-
-% function [alphahat,etahat,atilde, aK] = DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
-% Computes the diffuse kalman smoother 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
-% H: pp*pp matrix variance of measurement errors
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-%
-% OUTPUTS
-% alphahat: smoothed variables (a_{t|T})
-% epsilonhat:smoothed measurement errors
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t)
-% (meaningless for periods 1:d)
-% P: 3D array of one-step ahead forecast error variance
-% matrices
-% PK: 4D array of k-step ahead forecast error variance
-% matrices (meaningless for periods 1:d)
-% d: number of periods where filter remains in diffuse part
-% (should be equal to the order of integration of the model)
-% decomp: decomposition of the effect of shocks on filtered values
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% modified by M. Ratto:
-% new output argument aK (1-step to k-step predictions)
-% new options_.nk: the max step ahed prediction in aK (default is 4)
-% new crit1 value for rank of Pinf
-% it is assured that P is symmetric
-
-
-global options_
-
-d = 0;
-decomp = [];
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-aK = zeros(nk,mm,smpl+nk);
-PK = zeros(nk,mm,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-Fstar = zeros(pp,pp,smpl);
-iFinf = zeros(pp,pp,smpl);
-K = zeros(mm,pp,smpl);
-L = zeros(mm,mm,smpl);
-Linf = zeros(mm,mm,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-crit = options_.kalman_tol;
-crit1 = 1.e-8;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-epsilonhat = zeros(size(Y));
-r = zeros(mm,smpl+1);
-
-t = 0;
-while rank(Pinf(:,:,t+1),crit1) && t<smpl
- t = t+1;
- v(:,t)= Y(:,t) - Z*a(:,t);
- Finf = Z*Pinf(:,:,t)*Z';
- if rcond(Finf) < kalman_tol
- if ~all(abs(Finf(:)) < kalman_tol)
- % The univariate diffuse kalman filter should be used.
- return
- else
- Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H;
- if rcond(Fstar(:,:,t)) < kalman_tol
- if ~all(abs(Fstar(:,:,t))<kalman_tol)
- % The univariate diffuse kalman filter should be used.
- return
- else
- a(:,:,t+1) = T*a(:,:,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- end
- else
- iFstar = inv(Fstar(:,:,t));
- Kstar(:,:,t) = Pstar(:,:,t)*Z'*iFstar(:,:,t);
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*Z'*Kstar(:,:,t)')*T'+QQ;
- a(:,:,t+1) = T*(a(:,:,t)+Kstar(:,:,t)*v(:,t));
- end
- end
- else
-
- iFinf(:,:,t) = inv(F);
- PZI = Pinf(:,:,t)*Z'*iFinf(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- Kinf(:,:,t) = T*PZI;
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- % isn't a meaningless as long as we are in the diffuse part? MJ
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Linf(:,:,t) = T - Kinf(:,:,t)*Z;
- Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H;
- Kstar(:,:,t) = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-T*Pinf(:,:,t)*Z'*Kstar(:,:,t)' + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
- end
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-iFinf = iFinf(:,:,1:d);
-Linf = Linf(:,:,1:d);
-Fstar = Fstar(:,:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- v(:,t) = Y(:,t) - Z*a(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- F = Z*P(:,:,t)*Z' + H;
- if rcond(F) < crit
- return
- end
- iF(:,:,t) = inv(F);
- PZI = P(:,:,t)*Z'*iF(:,:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- K(:,:,t) = T*PZI;
- L(:,:,t) = T-K(:,:,t)*Z;
- a(:,t+1) = T*atilde(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
- P(:,:,t+1) = T*P(:,:,t)*T'-T*P(:,:,t)*Z'*K(:,:,t)' + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-if t<smpl
- PZI_s = PZI;
- K_s = K(:,:,t);
- iF_s = iF(:,:,t);
- P_s = P(:,:,t+1);
- P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
- iF = cat(3,iF(:,:,1:t),repmat(iF_s,[1 1 smpl-t]));
- L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t]));
- K = cat(3,K(:,:,1:t),repmat(T*P_s*Z'*iF_s,[1 1 smpl-t]));
-end
-while t<smpl
- t=t+1;
- v(:,t) = Y(:,t) - Z*a(:,t);
- atilde(:,t) = a(:,t) + PZI*v(:,t);
- a(:,t+1) = T*atilde(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
-end
-t = smpl+1;
-while t>d+1
- t = t-1;
- r(:,t) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t+1);
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-if d
- r0 = zeros(mm,d+1);
- r0(:,d+1) = r(:,d+1);
- r1 = zeros(mm,d+1);
- for t = d:-1:1
- r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
- r1(:,t) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t+1)) + Linf(:,:,t)'*r1(:,t+1);
- alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
- etahat(:,t) = QRt*r0(:,t);
- end
-end
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = Z'*iF(:,:,t)*v(:,t);
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- for h = 1:nk
- for j=1:rr
- eta=zeros(rr,1);
- eta(j) = eta_tm1t(j);
- decomp(h,:,j,t+h) = T^(h-1)*P(:,:,t)*Z'*ZRQinv*Z*R*eta;
- end
- end
- end
-end
-epsilonhat = Y-Z*alphahat;
diff --git a/matlab/DiffuseKalmanSmootherH3.m b/matlab/DiffuseKalmanSmootherH3.m
deleted file mode 100644
index ebc69d65166f13e53bafd9a2eea37d2c15b4078a..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmootherH3.m
+++ /dev/null
@@ -1,282 +0,0 @@
-function [alphahat,epsilonhat,etahat,a, aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% function [alphahat,epsilonhat,etahat,a1, aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% trend
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% mf: observed variables index in the state vector
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% epsilonhat:smoothed measurement error
-% etahat: smoothed shocks
-% a: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-%
-% 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-2010 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/>.
-
-% Modified by M. Ratto
-% New output argument aK: 1-step to nk-stpe ahed predictions)
-% New input argument nk: max order of predictions in aK
-% New global variable id_ where the DKF stops (common with
-% diffuselikelihood3)
-% New icc variable to count number of iterations for Finf steps
-% Pstar % Pinf simmetric
-% New termination of DKF iterations based on id_
-% Li also stored during DKF iterations !!
-% some bugs corrected in the DKF part of the smoother (Z matrix and
-% alphahat)
-
-global options_
-
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl);
-a1 = zeros(mm,smpl+1);
-aK = zeros(nk,mm,smpl+nk);
-
-if isempty(options_.diffuse_d),
- smpl_diff = 1;
-else
- smpl_diff=rank(Pinf1);
-end
-
-Fstar = zeros(pp,smpl_diff);
-Finf = zeros(pp,smpl_diff);
-Ki = zeros(mm,pp,smpl);
-Li = zeros(mm,mm,pp,smpl);
-Linf = zeros(mm,mm,pp,smpl_diff);
-L0 = zeros(mm,mm,pp,smpl_diff);
-Kstar = zeros(mm,pp,smpl_diff);
-P = zeros(mm,mm,smpl+1);
-P1 = P;
-Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
-Pstar1 = Pstar;
-Pinf1 = Pinf;
-crit = options_.kalman_tol;
-crit1 = 1.e-6;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-epsilonhat = zeros(size(Y));
-r = zeros(mm,smpl+1);
-
-Z = zeros(pp,mm);
-for i=1:pp;
- Z(i,mf(i)) = 1;
-end
-
-t = 0;
-icc=0;
-newRank = rank(Pinf(:,:,1),crit1);
-while newRank && t < smpl
- t = t+1;
- a(:,t) = a1(:,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
- Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1));
- Pstar1(:,:,t) = Pstar(:,:,t);
- Pinf1(:,:,t) = Pinf(:,:,t);
- for i=1:pp
- v(i,t) = Y(i,t)-a(mf(i),t)-trend(i,t);
- Fstar(i,t) = Pstar(mf(i),mf(i),t) + H(i,i);
- Finf(i,t) = Pinf(mf(i),mf(i),t);
- Kstar(:,i,t) = Pstar(:,mf(i),t);
- if Finf(i,t) > crit && newRank, % original MJ: if Finf(i,t) > crit
- icc=icc+1;
- Kinf(:,i,t) = Pinf(:,mf(i),t);
- Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
- L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Z(i,:)/Finf(i,t);
- a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) + ...
- Kinf(:,i,t)*transpose(Kinf(:,i,t))*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
- (Kstar(:,i,t)*transpose(Kinf(:,i,t)) +...
- Kinf(:,i,t)*transpose(Kstar(:,i,t)))/Finf(i,t);
- Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*transpose(Kinf(:,i,t))/Finf(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
- Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1));
- % new terminiation criteria by M. Ratto
- P0=Pinf(:,:,t);
- % newRank = any(diag(P0(mf,mf))>crit);
- % if newRank==0, options_.diffuse_d = i; end,
- if ~isempty(options_.diffuse_d),
- newRank = (icc<options_.diffuse_d);
- %if newRank && any(diag(P0(mf,mf))>crit)==0;
- if newRank && (any(diag(P0(mf,mf))>crit)==0 && rank(P0,crit1)==0);
- disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF')
- options_.diffuse_d = icc;
- newRank=0;
- end
- else
- %newRank = any(diag(P0(mf,mf))>crit);
- newRank = (any(diag(P0(mf,mf))>crit) | rank(P0,crit1));
- if newRank==0,
- options_.diffuse_d = icc;
- end
- end,
- % if newRank==0,
- % options_.diffuse_d=i; % this line is buggy!
- % end
- % end new terminiation criteria by M. Ratto
- elseif Fstar(i,t) > 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].
- Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother
- a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*transpose(Kstar(:,i,t))/Fstar(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1)
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- P0=Pinf(:,:,t+1);
- if newRank,
- %newRank = any(diag(P0(mf,mf))>crit);
- newRank = rank(P0,crit1);
- end
-end
-
-
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-Linf = Linf(:,:,:,1:d);
-L0 = L0(:,:,:,1:d);
-Fstar = Fstar(:,1:d);
-Finf = Finf(:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-Pstar1 = Pstar1(:,:,1:d);
-Pinf1 = Pinf1(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- a(:,t) = a1(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- P1(:,:,t) = P(:,:,t);
- for i=1:pp
- v(i,t) = Y(i,t) - a(mf(i),t) - trend(i,t);
- Fi(i,t) = P(mf(i),mf(i),t)+H(i,i);
- Ki(:,i,t) = P(:,mf(i),t);
- if Fi(i,t) > crit
- Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
- a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
- P(:,:,t) = P(:,:,t) - Ki(:,i,t)*transpose(Ki(:,i,t))/Fi(i,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
-Fi_s = Fi(:,t);
-Ki_s = Ki(:,:,t);
-L_s =Li(:,:,:,t);
-if t<smpl
- t_steady = t+1;
- P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
- Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
- Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
- Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
-end
-while t<smpl
- t=t+1;
- a(:,t) = a1(:,t);
- for i=1:pp
- v(i,t) = Y(i,t) - a(mf(i),t) - trend(i,t);
- if Fi_s(i) > crit
- a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
-end
-ri=zeros(mm,1);
-t = smpl+1;
-while t>d+1
- t = t-1;
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri=Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri;
- end
- end
- r(:,t) = ri;
- alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
- ri = T'*ri;
-end
-if d
- r0 = zeros(mm,d);
- r0(:,d) = ri;
- r1 = zeros(mm,d);
- for t = d:-1:1
- for i=pp:-1:1
- if Finf(i,t) > crit && ~(t==d && i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination
- r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
- L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t);
- r0(:,t) = Linf(:,:,i,t)'*r0(:,t);
- elseif Fstar(i,t) > crit % step needed whe Finf == 0
- r0(:,t)=Z(i,:)'/Fstar(i,t)*v(i,t)+Li(:,:,i,t)'*r0(:,t);
- end
- end
- alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
- r(:,t) = r0(:,t);
- etahat(:,t) = QRt*r(:,t);
- if t > 1
- r0(:,t-1) = transpose(T)*r0(:,t);
- r1(:,t-1) = transpose(T)*r1(:,t);
- end
- end
-end
-epsilonhat = Y-alphahat(mf,:)-trend;
-
-
diff --git a/matlab/DiffuseKalmanSmootherH3_Z.m b/matlab/DiffuseKalmanSmootherH3_Z.m
deleted file mode 100644
index cde9c56ea4a0455e7629240915b01003977fd6ed..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmootherH3_Z.m
+++ /dev/null
@@ -1,319 +0,0 @@
-function [alphahat,epsilonhat,etahat,a,P,aK,PK,d,decomp] = DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
-% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp_filt] = DiffuseKalmanSmootherH3(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% epsilonhat:smoothed measurement error
-% a: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-% (meaningless for periods 1:d)
-% P: 3D array of one-step ahead forecast error variance
-% matrices
-% PK: 4D array of k-step ahead forecast error variance
-% matrices (meaningless for periods 1:d)
-% d: number of periods where filter remains in diffuse part
-% (should be equal to the order of integration of the model)
-% decomp: decomposition of the effect of shocks on filtered values
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% Modified by M. Ratto
-% New output argument aK: 1-step to nk-stpe ahed predictions)
-% New input argument nk: max order of predictions in aK
-% New option options_.diffuse_d where the DKF stops (common with
-% diffuselikelihood3)
-% New icc variable to count number of iterations for Finf steps
-% Pstar % Pinf simmetric
-% New termination of DKF iterations based on options_.diffuse_d
-% Li also stored during DKF iterations !!
-% some bugs corrected in the DKF part of the smoother (Z matrix and
-% alphahat)
-
-global options_
-
-d = 0;
-decomp = [];
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl);
-a1 = zeros(mm,smpl+1);
-aK = zeros(nk,mm,smpl+nk);
-
-if isempty(options_.diffuse_d),
- smpl_diff = 1;
-else
- smpl_diff=rank(Pinf1);
-end
-
-Fstar = zeros(pp,smpl_diff);
-Finf = zeros(pp,smpl_diff);
-Ki = zeros(mm,pp,smpl);
-Li = zeros(mm,mm,pp,smpl);
-Linf = zeros(mm,mm,pp,smpl_diff);
-L0 = zeros(mm,mm,pp,smpl_diff);
-Kstar = zeros(mm,pp,smpl_diff);
-P = zeros(mm,mm,smpl+1);
-P1 = P;
-PK = zeros(nk,mm,mm,smpl+nk);
-Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
-Pstar1 = Pstar;
-Pinf1 = Pinf;
-crit = options_.kalman_tol;
-crit1 = 1.e-6;
-steady = smpl;
-rr = size(Q,1); % number of structural shocks
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-epsilonhat = zeros(size(Y));
-r = zeros(mm,smpl);
-
-t = 0;
-icc=0;
-newRank = rank(Pinf(:,:,1),crit1);
-while newRank && t < smpl
- t = t+1;
- a(:,t) = a1(:,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
- Pstar1(:,:,t) = Pstar(:,:,t);
- Pinf1(:,:,t) = Pinf(:,:,t);
- for i=1:pp
- Zi = Z(i,:);
- v(i,t) = Y(i,t)-Zi*a(:,t);
- Fstar(i,t) = Zi*Pstar(:,:,t)*Zi' +H(i,i);
- Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
- Kstar(:,i,t) = Pstar(:,:,t)*Zi';
- if Finf(i,t) > crit && newRank
- icc=icc+1;
- Kinf(:,i,t) = Pinf(:,:,t)*Zi';
- Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
- L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Zi/Finf(i,t);
- a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) + ...
- Kinf(:,i,t)*Kinf(:,i,t)'*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
- (Kstar(:,i,t)*Kinf(:,i,t)' +...
- Kinf(:,i,t)*Kstar(:,i,t)')/Finf(i,t);
- Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
- % new terminiation criteria by M. Ratto
- P0=Pinf(:,:,t);
- if ~isempty(options_.diffuse_d),
- newRank = (icc<options_.diffuse_d);
- if newRank && (any(diag(Z*P0*Z')>crit)==0 && rank(P0,crit1)==0);
- disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF')
- options_.diffuse_d = icc;
- newRank=0;
- end
- else
- newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1));
- if newRank==0,
- options_.diffuse_d = icc;
- end
- end,
- % end new terminiation criteria by M. Ratto
- elseif Fstar(i,t) > 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].
- Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother
- a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
- P0=Pinf(:,:,t+1);
- if newRank,
- newRank = rank(P0,crit1);
- end
-end
-
-
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-Linf = Linf(:,:,:,1:d);
-L0 = L0(:,:,:,1:d);
-Fstar = Fstar(:,1:d);
-Finf = Finf(:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-Pstar1 = Pstar1(:,:,1:d);
-Pinf1 = Pinf1(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- a(:,t) = a1(:,t);
- P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
- P1(:,:,t) = P(:,:,t);
- for i=1:pp
- Zi = Z(i,:);
- v(i,t) = Y(i,t) - Zi*a(:,t);
- Fi(i,t) = Zi*P(:,:,t)*Zi' + H(i,i);
- Ki(:,i,t) = P(:,:,t)*Zi';
- if Fi(i,t) > crit
- Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
- a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
- P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t);
- P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
- end
- end
- a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk,
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)) ;
- end
- end
- P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
-P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
-Fi_s = Fi(:,t);
-Ki_s = Ki(:,:,t);
-L_s =Li(:,:,:,t);
-if t<smpl
- P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
- P1 = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
- Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
- Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
- Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
-end
-while t<smpl
- t=t+1;
- a(:,t) = a1(:,t);
- for i=1:pp
- Zi = Z(i,:);
- v(i,t) = Y(i,t) - Zi*a(:,t);
- if Fi_s(i) > crit
- a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
- end
- end
- a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk,
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
-end
-ri=zeros(mm,1);
-t = smpl+1;
-while t>d+1
- t = t-1;
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri;
- end
- end
- r(:,t) = ri;
- alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
- ri = T'*ri;
-end
-if d
- r0 = zeros(mm,d);
- r0(:,d) = ri;
- r1 = zeros(mm,d);
- for t = d:-1:2
- for i=pp:-1:1
- % if Finf(i,t) > crit && ~(t==d && i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination
- if Finf(i,t) > crit
- r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
- L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t);
- r0(:,t) = Linf(:,:,i,t)'*r0(:,t);
- elseif Fstar(i,t) > crit % step needed whe Finf == 0
- r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+Li(:,:,i,t)'*r0(:,t);
- end
- end
- alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
- r(:,t) = r0(:,t);
- etahat(:,t) = QRt*r(:,t);
- if t > 1
- r0(:,t-1) = T'*r0(:,t);
- r1(:,t-1) = T'*r1(:,t);
- end
- end
-end
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = zeros(mm,1);
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri_d;
- end
- end
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- for h = 1:nk
- for j=1:rr
- eta=zeros(rr,1);
- eta(j) = eta_tm1t(j);
- decomp(h,:,j,t+h) = Ttok*P1(:,:,t)*Z'*ZRQinv*Z*R*eta;
- end
- Ttok = T*Ttok;
- end
- end
-end
-
-epsilonhat = Y-Z*alphahat;
diff --git a/matlab/DiffuseKalmanSmootherH3corr.m b/matlab/DiffuseKalmanSmootherH3corr.m
deleted file mode 100644
index f54ec2a0aeeb6344a429ce0db68e76ae40c3e68e..0000000000000000000000000000000000000000
--- a/matlab/DiffuseKalmanSmootherH3corr.m
+++ /dev/null
@@ -1,237 +0,0 @@
-function [alphahat,epsilonhat,etahat,a,aK] = DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-
-% function [alphahat,epsilonhat,etahat,a1] = DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% trend
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% mf: observed variables index in the state vector
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% epsilonhat:smoothed measurement error
-% etahat: smoothed shocks
-% a: matrix of updated variables (a_{t|t})
-% aK: matrix of one step ahead filtered state variables (a_{t+k|t})
-
-% SPECIAL REQUIREMENTS
-% See "Fast Filtering and Smoothing for Multivariate State Space
-% Models", S.J. Koopman and J. Durbin (2000, in Journal of Time Series
-% Analysis, vol. 21(3), pp. 281-296).
-
-% Copyright (C) 2004-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-global options_;
-
-nk = options_.nk;
-rr = size(Q,1);
-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(Pinf1,1) % Otherwise Pinf = 0 (no unit root)
- Pinf1 = cat(1,cat(2,Pinf1,zeros(mm,pp)),zeros(pp,mm+pp));
-end
-Pstar1 = cat(1,cat(2,Pstar1,zeros(mm,pp)),cat(2,zeros(pp,mm),H));
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-Pstar1 = Pstar;
-Pinf1 = Pinf;
-v = zeros(pp,smpl);
-a = zeros(mm+pp,smpl);
-a1 = zeros(mm+pp,smpl+1);
-aK = zeros(nk,mm,smpl+nk);
-Fstar = zeros(pp,smpl);
-Finf = zeros(pp,smpl);
-Fi = zeros(pp,smpl);
-Ki = zeros(mm+pp,pp,smpl);
-Li = zeros(mm+pp,mm+pp,pp,smpl);
-Linf = zeros(mm+pp,mm+pp,pp,smpl);
-L0 = zeros(mm+pp,mm+pp,pp,smpl);
-Kstar = zeros(mm+pp,pp,smpl);
-Kinf = zeros(mm+pp,pp,smpl);
-P = zeros(mm+pp,mm+pp,smpl+1);
-P1 = zeros(mm+pp,mm+pp,smpl+1);
-crit = options_.kalman_tol;
-steady = smpl;
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm+pp,smpl);
-etahat = zeros(rr,smpl);
-epsilonhat = zeros(pp,smpl);
-r = zeros(mm+pp,smpl+1);
-Z = zeros(pp,mm+pp);
-for i=1:pp;
- Z(i,mf(i)) = 1;
- Z(i,mm+i) = 1;
-end
-%% [1] Kalman filter...
-t = 0;
-newRank = rank(Pinf(:,:,1),crit);
-while newRank && t < smpl
- t = t+1;
- a(:,t) = a1(:,t);
- Pstar1(:,:,t) = Pstar(:,:,t);
- Pinf1(:,:,t) = Pinf(:,:,t);
- for i=1:pp
- v(i,t) = Y(i,t)-a(mf(i),t)-a(mm+i,t)-trend(i,t);
- Fstar(i,t) = Pstar(mf(i),mf(i),t)+Pstar(mm+i,mm+i,t);
- Finf(i,t) = Pinf(mf(i),mf(i),t);
- Kstar(:,i,t) = Pstar(:,mf(i),t)+Pstar(:,mm+i,t);
- if Finf(i,t) > crit
- Kinf(:,i,t) = Pinf(:,mf(i),t);
- Linf(:,:,i,t) = eye(mm+pp) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
- L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Z(i,:)/Finf(i,t);
- a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) + ...
- Kinf(:,i,t)*transpose(Kinf(:,i,t))*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
- (Kstar(:,i,t)*transpose(Kinf(:,i,t)) +...
- Kinf(:,i,t)*transpose(Kstar(:,i,t)))/Finf(i,t);
- Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*transpose(Kinf(:,i,t))/Finf(i,t);
- 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].
- a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*transpose(Kstar(:,i,t))/Fstar(i,t);
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1)
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
- P0=Pinf(:,:,t+1);
- newRank = ~all(abs(P0(:))<crit);
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-Linf = Linf(:,:,:,1:d);
-L0 = L0(:,:,:,1:d);
-Fstar = Fstar(:,1:d);
-Finf = Finf(:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-Pstar1 = Pstar1(:,:,1:d);
-Pinf1 = Pinf1(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- a(:,t) = a1(:,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- P1(:,:,t) = P(:,:,t);
- for i=1:pp
- v(i,t) = Y(i,t) - a(mf(i),t) - a(mm+i,t) - trend(i,t);
- Fi(i,t) = P(mf(i),mf(i),t)+P(mm+i,mm+i,t);
- Ki(:,i,t) = P(:,mf(i),t)+P(:,mm+i,t);
- if Fi(i,t) > crit
- Li(:,:,i,t) = eye(mm+pp)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
- a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
- P(:,:,t) = P(:,:,t) - Ki(:,i,t)*transpose(Ki(:,i,t))/Fi(i,t);
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
- notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
-Fi_s = Fi(:,t);
-Ki_s = Ki(:,:,t);
-L_s =Li(:,:,:,t);
-if t<smpl
- t_steady = t+1;
- P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
- Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
- Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
- Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
-end
-while t<smpl
- t=t+1;
- a(:,t) = a1(:,t);
- for i=1:pp
- v(i,t) = Y(i,t) - a(mf(i),t) - a(mm+i,t) - trend(i,t);
- if Fi_s(i) > crit
- a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
-end
-
-%% [2] Kalman smoother...
-ri=zeros(mm,1);
-t = smpl+1;
-while t>d+1
- t = t-1;
- for i=pp:-1:1
- if Fi(i,t) > crit
- ri=Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri;
- end
- end
- r(:,t) = ri(:,t);
- alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
- tmp = QRt*r(:,t);
- etahat(:,t) = tmp(1:rr);
- epsilonhat(:,t) = tmp(rr+(1:pp));
- ri = T'*ri;
-end
-if d
- r0 = zeros(mm+pp,d);
- r0(:,d) = ri;
- r1 = zeros(mm+pp,d);
- for t = d:-1:1
- for i=pp:-1:1
- if Finf(i,t) > crit
- r1(:,t) = transpose(Z)*v(:,t)/Finf(i,t) + ...
- L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t);
- r0(:,t) = transpose(Linf(:,:,i,t))*r0(:,t);
- end
- end
- alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
- r(:,t-1) = r0(:,t);
- tmp = QRt*r(:,t);
- etahat(:,t) = tmp(1:rr);
- epsilonhat(:,t) = tmp(rr+(1:pp));
- if t > 1
- r0(:,t-1) = T'*r0(:,t);
- r1(:,t-1) = T'*r1(:,t);
- end
- end
-end
-alphahat = alphahat(1:mm,:);
diff --git a/matlab/kalman/smoother/kalman_smoother.m b/matlab/kalman/smoother/kalman_smoother.m
deleted file mode 100644
index 1e9bbb2e015de7e2bab61b3658db248b74ba94b5..0000000000000000000000000000000000000000
--- a/matlab/kalman/smoother/kalman_smoother.m
+++ /dev/null
@@ -1,188 +0,0 @@
-function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = kalman_smoother(T,R,Q,H,P0,Y,start,mf,kalman_tol,riccati_tol)
-
-% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = kalman_smoother(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol)
-% Computes the kalman smoother of a stationary state space model.
-%
-% INPUTS
-% T [double] mm*mm transition matrix of the state equation.
-% R [double] mm*rr matrix, mapping structural innovations to state variables.
-% Q [double] rr*rr covariance matrix of the structural innovations.
-% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
-% P0 [double] mm*mm variance-covariance matrix with stationary variables
-% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
-% start [integer] scalar, likelihood evaluation starts at 'start'.
-% mf [integer] pp*1 vector of indices.
-% kalman_tol [double] scalar, tolerance parameter (rcond).
-% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
-%
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-global options_
-
-option_filter_covariance = options_.filter_covariance;
-option_filter_decomposition = options_.filter_decomposition;
-
-nk = options_.nk;
-smpl = size(Y,2); % Sample size.
-mm = size(T,2); % Number of state variables.
-pp = size(Y,1); % Maximum number of
- % observed variables.
-rr = size(Q,1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-K = zeros(mm,pp,smpl);
-aK = zeros(nk,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-QQ = R*Q*R';
-QRt = Q*R';
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-epsilonhat = zeros(rr,smpl);
-r = zeros(mm,smpl+1);
-oldK = 0;
-
-if option_filter_covariance
- PK = zeros(nk,mm,mm,smpl+nk);
-else
- PK = [];
-end
-
-if option_filter_decomposition
- decomp = zeros(nk,mm,rr,smpl+nk);
-else
- decomp = [];
-end
-
-P(:,:,1) = P0;
-
-t = 0;
-notsteady = 1;
-F_singular = 1;
-while notsteady && t<smpl
- t = t+1;
- v(:,t) = Y(:,t)-a(mf,t);
- F = P(mf,mf,t) + H;
- if rcond(F) < kalman_tol
- if ~all(abs(F(:))<kalman_tol)
- return
- else
- atilde(:,t) = a(:,t);
- a(:,t+1) = T*a(:,t);
- P(:,:,t+1) = T*P(:,:,t)*T'+QQ;
- end
- else
- F_singular = 0;
- iF(:,:,t) = inv(F);
- K1 = P(:,mf,t)*iF(:,:,t);
- atilde(:,t) = a(:,t) + K1*v(:,t);
- K(:,:,t) = T*K1;
- a(:,t+1) = T*atilde(:,t);
- P(:,:,t+1) = (T*P(:,:,t)-K(:,:,t)*P(mf,:,t))*T'+QQ;
- end
- aK(1,:,t+1) = a(:,t+1);
- if option_filter_covariance
- Pf = P(:,:,t);
- Pf = T*Pf*T' + QQ;
- PK(1,:,:,t+1) = Pf;
- end
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- if option_filter_covariance
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- end
- end
- notsteady = max(max(abs(K(:,:,t)-oldK))) > riccati_tol;
- oldK = K(:,:,t);
-end
-
-if F_singular
- error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-if t < smpl
- t0 = t;
- while t < smpl
- t = t+1;
- v(:,t) = Y(:,t)-a(mf,t);
- atilde(:,t) = a(:,t) + K1*v(:,t);
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- if option_filter_covariance
- Pf = P(:,:,t);
- Pf = T*Pf*T' + QQ;
- PK(1,:,:,t+1) = Pf;
- end
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- if option_filter_covariance
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- end
- end
- end
-
- K= cat(3,K(:,:,1:t0),repmat(K(:,:,t0),[1 1 smpl-t0+1]));
- P = cat(3,P(:,:,1:t0),repmat(P(:,:,t0),[1 1 smpl-t0+1]));
- iF = cat(3,iF(:,:,1:t0),repmat(iF(:,:,t0),[1 1 smpl-t0+1]));
-end
-
-t = smpl+1;
-while t>1
- t = t-1;
- r(:,t) = T'*r(:,t+1);
- r(mf,t) = r(mf,t)+iF(:,:,t)*v(:,t) - K(:,:,t)'*r(:,t+1);
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-epsilonhat = Y-alphahat(mf,:);
-
-if option_filter_decomposition
- ZRQinv = inv(QQ(mf,mf));
- for t = 1:smpl
- % calculate eta_tm1t
- eta = QRt(:,mf)*iF(:,:,t)*v(:,t);
- AAA = P(:,mf,t)*ZRQinv*bsxfun(@times,R(mf,:),eta');
- % calculate decomposition
- Ttok = eye(mm,mm);
- decomp(1,:,:,t+1) = AAA;
- for h = 2:nk
- AAA = T*AAA;
- decomp(h,:,:,t+h) = AAA;
- end
- end
-end
-
-if ~option_filter_covariance
- P = [];
-end
-
diff --git a/matlab/missing_DiffuseKalmanSmoother1.m b/matlab/missing_DiffuseKalmanSmoother1.m
deleted file mode 100644
index 52a9cf1ec593ccdc8b41d28d568be370e387221c..0000000000000000000000000000000000000000
--- a/matlab/missing_DiffuseKalmanSmoother1.m
+++ /dev/null
@@ -1,207 +0,0 @@
-function [alphahat,etahat,atilde, aK] = missing_DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf,data_index)
-
-% function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% trend
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% mf: observed variables index in the state vector
-% data_index [cell] 1*smpl cell of column vectors of indices.
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% modified by M. Ratto:
-% new output argument aK (1-step to k-step predictions)
-% new options_.nk: the max step ahed prediction in aK (default is 4)
-% new crit1 value for rank of Pinf
-% it is assured that P is symmetric
-
-
-global options_
-
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-aK = zeros(nk,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-Fstar = zeros(pp,pp,smpl);
-iFinf = zeros(pp,pp,smpl);
-K = zeros(mm,pp,smpl);
-L = zeros(mm,mm,smpl);
-Linf = zeros(mm,mm,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-crit = options_.kalman_tol;
-crit1 = 1.e-8;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl+1);
-
-Z = zeros(pp,mm);
-for i=1:pp;
- Z(i,mf(i)) = 1;
-end
-
-t = 0;
-while rank(Pinf(:,:,t+1),crit1) && t<smpl
- t = t+1;
- di = data_index{t};
- if isempty(di)
- atilde(:,t) = a(:,t);
- Linf(:,:,t) = T;
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T' + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
- else
- mf1 = mf(di);
- v(di,t)= Y(di,t) - a(mf1,t) - trend(di,t);
- if rcond(Pinf(mf1,mf1,t)) < crit
- return
- end
- iFinf(di,di,t) = inv(Pinf(mf1,mf1,t));
- PZI = Pinf(:,mf1,t)*iFinf(di,di,t);
- atilde(:,t) = a(:,t) + PZI*v(di,t);
- Kinf(:,di,t) = T*PZI;
- a(:,t+1) = T*atilde(:,t);
- Linf(:,:,t) = T - Kinf(:,di,t)*Z(di,:);
- Fstar(di,di,t) = Pstar(mf1,mf1,t);
- Kstar(:,di,t) = (T*Pstar(:,mf1,t)-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf1,t)*transpose(Kinf(:,di,t))-Kinf(:,di,t)*Pinf(mf1,mf1,t)*transpose(Kstar(:,di,t)) + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf1,t)* ...
- transpose(Kinf(:,di,t));
- end
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-iFinf = iFinf(:,:,1:d);
-Linf = Linf(:,:,1:d);
-Fstar = Fstar(:,:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- di = data_index{t};
- if isempty(di)
- atilde(:,t) = a(:,t);
- L(:,:,t) = T;
- P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
- else
- mf1 = mf(di);
- v(di,t) = Y(di,t) - a(mf1,t) - trend(di,t);
- if rcond(P(mf1,mf1,t)) < crit
- return
- end
- iF(di,di,t) = inv(P(mf1,mf1,t));
- PZI = P(:,mf1,t)*iF(di,di,t);
- atilde(:,t) = a(:,t) + PZI*v(di,t);
- K(:,di,t) = T*PZI;
- L(:,:,t) = T-K(:,di,t)*Z(di,:);
- a(:,t+1) = T*atilde(:,t);
- end
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ;
- % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-% $$$ if t<smpl
-% $$$ PZI_s = PZI;
-% $$$ K_s = K(:,:,t);
-% $$$ iF_s = iF(:,:,t);
-% $$$ P_s = P(:,:,t+1);
-% $$$ t_steady = t+1;
-% $$$ P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
-% $$$ iF = cat(3,iF(:,:,1:t),repmat(inv(P_s(mf,mf)),[1 1 smpl-t_steady+1]));
-% $$$ L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t_steady+1]));
-% $$$ K = cat(3,K(:,:,1:t),repmat(T*P_s(:,mf)*iF_s,[1 1 smpl-t_steady+1]));
-% $$$ end
-% $$$ while t<smpl
-% $$$ t=t+1;
-% $$$ v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
-% $$$ atilde(:,t) = a(:,t) + PZI*v(:,t);
-% $$$ a(:,t+1) = T*a(:,t) + K_s*v(:,t);
-% $$$ aK(1,:,t+1) = a(:,t+1);
-% $$$ for jnk=2:nk,
-% $$$ aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
-% $$$ end
-% $$$ end
-t = smpl+1;
-while t>d+1
- t = t-1;
- di = data_index{t};
- if isempty(di)
- r(:,t) = L(:,:,t)'*r(:,t+1);
- else
- r(:,t) = Z(di,:)'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1);
- end
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-if d
- r0 = zeros(mm,d+1);
- r0(:,d+1) = r(:,d+1);
- r1 = zeros(mm,d+1);
- for t = d:-1:1
- r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
- di = data_index{t};
- if isempty(di)
- r1(:,t) = Linf(:,:,t)'*r1(:,t+1);
- else
- r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ...
- + Linf(:,:,t)'*r1(:,t+1);
- end
- alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
- etahat(:,t) = QRt*r0(:,t);
- end
-end
diff --git a/matlab/missing_DiffuseKalmanSmoother1_Z.m b/matlab/missing_DiffuseKalmanSmoother1_Z.m
deleted file mode 100644
index 47b1102f079ab4ba6df33a179ef838597a2729f6..0000000000000000000000000000000000000000
--- a/matlab/missing_DiffuseKalmanSmoother1_Z.m
+++ /dev/null
@@ -1,241 +0,0 @@
-function [alphahat,etahat,atilde,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother1_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl,data_index)
-
-% function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% data_index [cell] 1*smpl cell of column vectors of indices.
-%
-% OUTPUTS
-% alphahat: smoothed variables (a_{t|T})
-% etahat: smoothed shocks
-% atilde: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t)
-% (meaningless for periods 1:d)
-% P: 3D array of one-step ahead forecast error variance
-% matrices
-% PK: 4D array of k-step ahead forecast error variance
-% matrices (meaningless for periods 1:d)
-% d: number of periods where filter remains in diffuse part
-% (should be equal to the order of integration of the model)
-% decomp: decomposition of the effect of shocks on filtered values
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% modified by M. Ratto:
-% new output argument aK (1-step to k-step predictions)
-% new options_.nk: the max step ahed prediction in aK (default is 4)
-% new crit1 value for rank of Pinf
-% it is assured that P is symmetric
-
-
-global options_
-
-d = 0;
-decomp = [];
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-atilde = zeros(mm,smpl);
-aK = zeros(nk,mm,smpl+nk);
-PK = zeros(nk,mm,mm,smpl+nk);
-iF = zeros(pp,pp,smpl);
-Fstar = zeros(pp,pp,smpl);
-iFinf = zeros(pp,pp,smpl);
-K = zeros(mm,pp,smpl);
-L = zeros(mm,mm,smpl);
-Linf = zeros(mm,mm,smpl);
-Kstar = zeros(mm,pp,smpl);
-P = zeros(mm,mm,smpl+1);
-Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
-crit = options_.kalman_tol;
-crit1 = 1.e-8;
-steady = smpl;
-rr = size(Q,1);
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl+1);
-
-t = 0;
-while rank(Pinf(:,:,t+1),crit1) && t<smpl
- t = t+1;
- di = data_index{t};
- if isempty(di)
- atilde(:,t) = a(:,t);
- Linf(:,:,t) = T;
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T' + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
- else
- ZZ = Z(di,:);
- v(di,t)= Y(di,t) - ZZ*a(:,t);
- F = ZZ*Pinf(:,:,t)*ZZ';
- if rcond(F) < crit
- return
- end
- iFinf(di,di,t) = inv(F);
- PZI = Pinf(:,:,t)*ZZ'*iFinf(di,di,t);
- atilde(:,t) = a(:,t) + PZI*v(di,t);
- Kinf(:,di,t) = T*PZI;
- Linf(:,:,t) = T - Kinf(:,di,t)*ZZ;
- Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ';
- Kstar(:,di,t) = (T*Pstar(:,:,t)*ZZ'-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t);
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*ZZ'*Kinf(:,di,t)'-Kinf(:,di,t)*F*Kstar(:,di,t)' + QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*ZZ'*Kinf(:,di,t)';
- end
- a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- % isn't a meaningless as long as we are in the diffuse part? MJ
- for jnk=2:nk,
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
-end
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-iFinf = iFinf(:,:,1:d);
-Linf = Linf(:,:,1:d);
-Fstar = Fstar(:,:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
- di = data_index{t};
- if isempty(di)
- atilde(:,t) = a(:,t);
- L(:,:,t) = T;
- P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
- else
- ZZ = Z(di,:);
- v(di,t) = Y(di,t) - ZZ*a(:,t);
- F = ZZ*P(:,:,t)*ZZ';
- if rcond(F) < crit
- return
- end
- iF(di,di,t) = inv(F);
- PZI = P(:,:,t)*ZZ'*iF(di,di,t);
- atilde(:,t) = a(:,t) + PZI*v(di,t);
- K(:,di,t) = T*PZI;
- L(:,:,t) = T-K(:,di,t)*ZZ;
- P(:,:,t+1) = T*P(:,:,t)*T'-T*P(:,:,t)*ZZ'*K(:,di,t)' + QQ;
- end
- a(:,t+1) = T*atilde(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
- % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-% $$$ if t<smpl
-% $$$ PZI_s = PZI;
-% $$$ K_s = K(:,:,t);
-% $$$ iF_s = iF(:,:,t);
-% $$$ P_s = P(:,:,t+1);
-% $$$ P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
-% $$$ iF = cat(3,iF(:,:,1:t),repmat(iF_s,[1 1 smpl-t]));
-% $$$ L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t]));
-% $$$ K = cat(3,K(:,:,1:t),repmat(T*P_s*Z'*iF_s,[1 1 smpl-t]));
-% $$$ end
-% $$$ while t<smpl
-% $$$ t=t+1;
-% $$$ v(:,t) = Y(:,t) - Z*a(:,t);
-% $$$ atilde(:,t) = a(:,t) + PZI*v(:,t);
-% $$$ a(:,t+1) = T*atilde(:,t);
-% $$$ Pf = P(:,:,t);
-% $$$ for jnk=1:nk,
-% $$$ Pf = T*Pf*T' + QQ;
-% $$$ aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
-% $$$ PK(jnk,:,:,t+jnk) = Pf;
-% $$$ end
-% $$$ end
-t = smpl+1;
-while t>d+1
- t = t-1;
- di = data_index{t};
- if isempty(di)
- r(:,t) = L(:,:,t)'*r(:,t+1);
- else
- ZZ = Z(di,:);
- r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1);
- end
- alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
-end
-if d
- r0 = zeros(mm,d+1);
- r0(:,d+1) = r(:,d+1);
- r1 = zeros(mm,d+1);
- for t = d:-1:1
- r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
- di = data_index{t};
- if isempty(di)
- r1(:,t) = Linf(:,:,t)'*r1(:,t+1);
- else
- r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ...
- + Linf(:,:,t)'*r1(:,t+1);
- end
- alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
- etahat(:,t) = QRt*r0(:,t);
- end
-end
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = Z'*iF(:,:,t)*v(:,t);
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- for h = 1:nk
- for j=1:rr
- eta=zeros(rr,1);
- eta(j) = eta_tm1t(j);
- decomp(h,:,j,t+h) = T^(h-1)*P(:,:,t)*Z'*ZRQinv*Z*R*eta;
- end
- end
- end
-end
diff --git a/matlab/missing_DiffuseKalmanSmoother3_Z.m b/matlab/missing_DiffuseKalmanSmoother3_Z.m
deleted file mode 100644
index 75ff3a7bdac14a7355d969a7537e0a0f3a2fd608..0000000000000000000000000000000000000000
--- a/matlab/missing_DiffuseKalmanSmoother3_Z.m
+++ /dev/null
@@ -1,323 +0,0 @@
-function [alphahat,etahat,a,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl,data_index)
-% function [alphahat,etahat,a1,P,aK,PK,d,decomp_filt] = missing_DiffuseKalmanSmoother3_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl)
-% Computes the diffuse kalman smoother 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
-% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
-% Pstar1: mm*mm variance-covariance matrix with stationary variables
-% Y: pp*1 vector
-% pp: number of observed variables
-% mm: number of state variables
-% smpl: sample size
-% data_index [cell] 1*smpl cell of column vectors of indices.
-%
-% OUTPUTS
-% alphahat: smoothed state variables (a_{t|T})
-% etahat: smoothed shocks
-% a: matrix of updated variables (a_{t|t})
-% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
-% (meaningless for periods 1:d)
-% P: 3D array of one-step ahead forecast error variance
-% matrices
-% PK: 4D array of k-step ahead forecast error variance
-% matrices (meaningless for periods 1:d)
-% d: number of periods where filter remains in diffuse part
-% (should be equal to the order of integration of the model)
-% decomp: decomposition of the effect of shocks on filtered values
-%
-% 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-2011 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-% Modified by M. Ratto
-% New output argument aK: 1-step to nk-stpe ahed predictions)
-% New input argument nk: max order of predictions in aK
-% New option options_.diffuse_d where the DKF stops (common with
-% diffuselikelihood3)
-% New icc variable to count number of iterations for Finf steps
-% Pstar % Pinf simmetric
-% New termination of DKF iterations based on options_.diffuse_d
-% Li also stored during DKF iterations !!
-% some bugs corrected in the DKF part of the smoother (Z matrix and
-% alphahat)
-
-global options_
-
-d = 0;
-decomp = [];
-nk = options_.nk;
-spinf = size(Pinf1);
-spstar = size(Pstar1);
-v = zeros(pp,smpl);
-a = zeros(mm,smpl);
-a1 = zeros(mm,smpl+1);
-aK = zeros(nk,mm,smpl+nk);
-
-if isempty(options_.diffuse_d),
- smpl_diff = 1;
-else
- smpl_diff=rank(Pinf1);
-end
-
-Fstar = zeros(pp,smpl_diff);
-Finf = zeros(pp,smpl_diff);
-Fi = zeros(pp,smpl);
-Ki = zeros(mm,pp,smpl);
-%Li = zeros(mm,mm,pp,smpl);
-%Linf = zeros(mm,mm,pp,smpl_diff);
-%L0 = zeros(mm,mm,pp,smpl_diff);
-Kstar = zeros(mm,pp,smpl_diff);
-P = zeros(mm,mm,smpl+1);
-P1 = P;
-PK = zeros(nk,mm,mm,smpl+nk);
-Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1;
-Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
-Pstar1 = Pstar;
-Pinf1 = Pinf;
-crit = options_.kalman_tol;
-crit1 = 1.e-6;
-steady = smpl;
-rr = size(Q,1); % number of structural shocks
-QQ = R*Q*transpose(R);
-QRt = Q*transpose(R);
-alphahat = zeros(mm,smpl);
-etahat = zeros(rr,smpl);
-r = zeros(mm,smpl);
-
-t = 0;
-icc=0;
-newRank = rank(Pinf(:,:,1),crit1);
-while newRank && t < smpl
- t = t+1;
- a(:,t) = a1(:,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
- Pstar1(:,:,t) = Pstar(:,:,t);
- Pinf1(:,:,t) = Pinf(:,:,t);
- di = data_index{t}';
- for i=di
- Zi = Z(i,:);
- v(i,t) = Y(i,t)-Zi*a(:,t);
- Fstar(i,t) = Zi*Pstar(:,:,t)*Zi';
- Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
- Kstar(:,i,t) = Pstar(:,:,t)*Zi';
- if Finf(i,t) > crit && newRank
- icc=icc+1;
- Kinf(:,i,t) = Pinf(:,:,t)*Zi';
- % Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
- % L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Zi/Finf(i,t);
- a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) + ...
- Kinf(:,i,t)*Kinf(:,i,t)'*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
- (Kstar(:,i,t)*Kinf(:,i,t)' +...
- Kinf(:,i,t)*Kstar(:,i,t)')/Finf(i,t);
- Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
- % new terminiation criteria by M. Ratto
- P0=Pinf(:,:,t);
- if ~isempty(options_.diffuse_d),
- newRank = (icc<options_.diffuse_d);
- if newRank && (any(diag(Z*P0*Z')>crit)==0 && rank(P0,crit1)==0);
- disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF')
- options_.diffuse_d = icc;
- newRank=0;
- end
- else
- newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1));
- if newRank==0,
- options_.diffuse_d = icc;
- end
- end,
- % end new terminiation criteria by M. Ratto
- elseif Fstar(i,t) > 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].
- %Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother
- a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
- Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t);
- Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
- end
- end
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
- Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
- P0=Pinf(:,:,t+1);
- if newRank,
- newRank = rank(P0,crit1);
- end
-end
-
-
-d = t;
-P(:,:,d+1) = Pstar(:,:,d+1);
-%Linf = Linf(:,:,:,1:d);
-%L0 = L0(:,:,:,1:d);
-Fstar = Fstar(:,1:d);
-Finf = Finf(:,1:d);
-Kstar = Kstar(:,:,1:d);
-Pstar = Pstar(:,:,1:d);
-Pinf = Pinf(:,:,1:d);
-Pstar1 = Pstar1(:,:,1:d);
-Pinf1 = Pinf1(:,:,1:d);
-notsteady = 1;
-while notsteady && t<smpl
- t = t+1;
- a(:,t) = a1(:,t);
- P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
- P1(:,:,t) = P(:,:,t);
- di = data_index{t}';
- for i=di
- Zi = Z(i,:);
- v(i,t) = Y(i,t) - Zi*a(:,t);
- Fi(i,t) = Zi*P(:,:,t)*Zi';
- Ki(:,i,t) = P(:,:,t)*Zi';
- if Fi(i,t) > crit
- %Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
- a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
- P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t);
- P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
- end
- end
- a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk
- Pf = T*Pf*T' + QQ;
- PK(jnk,:,:,t+jnk) = Pf;
- if jnk>1
- aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
- end
- end
- P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
- % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
-% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
-% $$$ Fi_s = Fi(:,t);
-% $$$ Ki_s = Ki(:,:,t);
-% $$$ L_s =Li(:,:,:,t);
-% $$$ if t<smpl
-% $$$ P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
-% $$$ P1 = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
-% $$$ Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
-% $$$ Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
-% $$$ Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
-% $$$ end
-% $$$ while t<smpl
-% $$$ t=t+1;
-% $$$ a(:,t) = a1(:,t);
-% $$$ di = data_index{t}';
-% $$$ for i=di
-% $$$ Zi = Z(i,:);
-% $$$ v(i,t) = Y(i,t) - Zi*a(:,t);
-% $$$ if Fi_s(i) > crit
-% $$$ a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
-% $$$ end
-% $$$ end
-% $$$ a1(:,t+1) = T*a(:,t);
-% $$$ Pf = P(:,:,t);
-% $$$ for jnk=1:nk,
-% $$$ Pf = T*Pf*T' + QQ;
-% $$$ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
-% $$$ PK(jnk,:,:,t+jnk) = Pf;
-% $$$ end
-% $$$ end
-ri=zeros(mm,1);
-t = smpl+1;
-while t > d+1
- t = t-1;
- di = flipud(data_index{t})';
- for i = di
- if Fi(i,t) > crit
- ri = Z(i,:)'/Fi(i,t)*v(i,t)+ri-Ki(:,i,t)'*ri/Fi(i,t)*Z(i,:)';
- end
- end
- r(:,t) = ri;
- alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
- etahat(:,t) = QRt*r(:,t);
- ri = T'*ri;
-end
-if d
- r0 = zeros(mm,d);
- r0(:,d) = ri;
- r1 = zeros(mm,d);
- for t = d:-1:1
- di = flipud(data_index{t})';
- for i = di
- % if Finf(i,t) > crit && ~(t==d && i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination
- if Finf(i,t) > crit
- r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
- (Kinf(:,i,t)'*Fstar(i,t)/Finf(i,t)-Kstar(:,i,t)')*r0(:,t)/Finf(i,t)*Z(i,:)' + ...
- r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)';
- r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)';
- elseif Fstar(i,t) > crit % step needed whe Finf == 0
- r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+r0(:,t)-(Kstar(:,i,t)'*r0(:,t))/Fstar(i,t)*Z(i,:)';
- end
- end
- alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
- r(:,t) = r0(:,t);
- etahat(:,t) = QRt*r(:,t);
- if t > 1
- r0(:,t-1) = T'*r0(:,t);
- r1(:,t-1) = T'*r1(:,t);
- end
- end
-end
-
-disp('smoother done');
-
-if nargout > 7
- decomp = zeros(nk,mm,rr,smpl+nk);
- ZRQinv = inv(Z*QQ*Z');
- for t = max(d,1):smpl
- ri_d = zeros(mm,1);
- di = flipud(data_index{t})';
- for i = di
- if Fi(i,t) > crit
- ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)';
- end
- end
-
- % calculate eta_tm1t
- eta_tm1t = QRt*ri_d;
- % calculate decomposition
- Ttok = eye(mm,mm);
- AAA = P1(:,:,t)*Z'*ZRQinv*Z*R;
- for h = 1:nk
- BBB = Ttok*AAA;
- for j=1:rr
- decomp(h,:,j,t+h) = eta_tm1t(j)*BBB(:,j);
- end
- Ttok = T*Ttok;
- end
- end
-end