diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 826fcc4d6..eab3545f6 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -264,13 +264,25 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
end
else
if kalman_algo == 1
- [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+ if missing_value
+ [alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother1(T,R,Q, ...
+ Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
+ else
+ [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q, ...
+ Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+ end
if all(alphahat(:)==0)
kalman_algo = 2;
end
end
if kalman_algo == 2
- [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+ if missing_value
+ [alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother3(T,R,Q, ...
+ Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
+ else
+ [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q, ...
+ Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+ end
end
if kalman_algo == 3
data1 = Y - trend;
diff --git a/matlab/describe_missing_data.m b/matlab/describe_missing_data.m
index b311acd80..97d16a9f5 100644
--- a/matlab/describe_missing_data.m
+++ b/matlab/describe_missing_data.m
@@ -17,6 +17,6 @@ number_of_observations = length(variable_index);
if ~missing_observations_counter
no_more_missing_observations = 0;
else
- tmp = find(missing_observations_counter>=(gend*nvarobs-number_of_observations))
+ tmp = find(missing_observations_counter>=(gend*nvarobs-number_of_observations));
no_more_missing_observations = tmp(1);
end
\ No newline at end of file
diff --git a/matlab/missing_DiffuseKalmanSmoother1.m b/matlab/missing_DiffuseKalmanSmoother1.m
new file mode 100644
index 000000000..d86830f52
--- /dev/null
+++ b/matlab/missing_DiffuseKalmanSmoother1.m
@@ -0,0 +1,207 @@
+function [alphahat,etahat,atilde, aK] = 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-2008 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+% 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+1);
+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) & td+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
new file mode 100644
index 000000000..60edcd8f0
--- /dev/null
+++ b/matlab/missing_DiffuseKalmanSmoother1_Z.m
@@ -0,0 +1,238 @@
+function [alphahat,etahat,atilde,P,aK,PK,d,decomp] = 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-2008 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+% 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) & td+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.m b/matlab/missing_DiffuseKalmanSmoother3.m
new file mode 100644
index 000000000..13efff232
--- /dev/null
+++ b/matlab/missing_DiffuseKalmanSmoother3.m
@@ -0,0 +1,284 @@
+function [alphahat,etahat,a, aK] = missing_DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf,data_index)
+% function [alphahat,etahat,a1, aK] = missing_DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf,data_index)
+% 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
+% 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})
+%
+% SPECIAL REQUIREMENTS
+% See "Filtering and Smoothing of State Vector for Diffuse State Space
+% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
+% Analysis, vol. 24(1), pp. 85-98).
+
+% Copyright (C) 2004-2008 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+% 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_
+
+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);
+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);
+ di = data_index{t}';
+ for i=di
+ 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));
+ % new terminiation criteria by M. Ratto
+ P0=Pinf(:,:,t);
+ % newRank = any(diag(P0(mf,mf))>crit);
+ % if newRank==0, id = i; end,
+ if ~isempty(options_.diffuse_d),
+ newRank = (icccrit)==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 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);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ 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 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);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ end
+ P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
+ notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t)))) crit
+ a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
+ end
+ end
+ a1(:,t+1) = T*a(:,t);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ 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)+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
+ 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
+ %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
+
+
diff --git a/matlab/missing_DiffuseKalmanSmoother3_Z.m b/matlab/missing_DiffuseKalmanSmoother3_Z.m
new file mode 100644
index 000000000..e7c74a7a4
--- /dev/null
+++ b/matlab/missing_DiffuseKalmanSmoother3_Z.m
@@ -0,0 +1,316 @@
+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-2008 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+% 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;
+aK = zeros(nk,mm,smpl+nk);
+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 = (icccrit)==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);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ 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 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);
+ for jnk=1:nk,
+ Pf = T*Pf*T' + QQ;
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ PK(jnk,:,:,t+jnk) = Pf;
+ end
+ P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
+ notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t)))) 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)+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
+ 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) + ...
+ 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);
+ di = flipud(data_index{t})';
+ for i = di
+ 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