1function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag) 2 3% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag) 4% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix. 5% 6% INPUTS 7% T: mm*mm matrix 8% Z: pp*mm matrix 9% R: mm*rr matrix 10% Q: rr*rr matrix 11% H: pp*pp matrix variance of measurement errors 12% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros 13% Pstar1: mm*mm variance-covariance matrix with stationary variables 14% Y: pp*1 vector 15% pp: number of observed variables 16% mm: number of state variables 17% smpl: sample size 18% data_index [cell] 1*smpl cell of column vectors of indices. 19% nk number of forecasting periods 20% kalman_tol tolerance for reciprocal condition number 21% diffuse_kalman_tol tolerance for reciprocal condition number (for Finf) and the rank of Pinf 22% decomp_flag if true, compute filter decomposition 23% state_uncertainty_flag if true, compute uncertainty about smoothed 24% state estimate 25% 26% OUTPUTS 27% alphahat: smoothed variables (a_{t|T}) 28% epsilonhat:smoothed measurement errors 29% etahat: smoothed shocks 30% atilde: matrix of updated variables (a_{t|t}) 31% aK: 3D array of k step ahead filtered state variables (a_{t+k|t) 32% (meaningless for periods 1:d) 33% P: 3D array of one-step ahead forecast error variance 34% matrices 35% PK: 4D array of k-step ahead forecast error variance 36% matrices (meaningless for periods 1:d) 37% decomp: decomposition of the effect of shocks on filtered values 38% V: 3D array of state uncertainty matrices 39% 40% Notes: 41% Outputs are stored in decision-rule order, i.e. to get variables in order of declaration 42% as in M_.endo_names, ones needs code along the lines of: 43% variables_declaration_order(dr.order_var,:) = alphahat 44% 45% SPECIAL REQUIREMENTS 46% See "Filtering and Smoothing of State Vector for Diffuse State Space 47% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series 48% Analysis, vol. 24(1), pp. 85-98). 49% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press, 50% Second Edition, Ch. 5 51 52% Copyright (C) 2004-2018 Dynare Team 53% 54% This file is part of Dynare. 55% 56% Dynare is free software: you can redistribute it and/or modify 57% it under the terms of the GNU General Public License as published by 58% the Free Software Foundation, either version 3 of the License, or 59% (at your option) any later version. 60% 61% Dynare is distributed in the hope that it will be useful, 62% but WITHOUT ANY WARRANTY; without even the implied warranty of 63% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 64% GNU General Public License for more details. 65% 66% You should have received a copy of the GNU General Public License 67% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 68 69% modified by M. Ratto: 70% new output argument aK (1-step to k-step predictions) 71% new options_.nk: the max step ahed prediction in aK (default is 4) 72% new crit1 value for rank of Pinf 73% it is assured that P is symmetric 74 75d = 0; 76decomp = []; 77spinf = size(Pinf1); 78spstar = size(Pstar1); 79v = zeros(pp,smpl); 80a = zeros(mm,smpl+1); 81atilde = zeros(mm,smpl); 82aK = zeros(nk,mm,smpl+nk); 83PK = zeros(nk,mm,mm,smpl+nk); 84iF = zeros(pp,pp,smpl); 85Fstar = zeros(pp,pp,smpl); 86iFstar = zeros(pp,pp,smpl); 87iFinf = zeros(pp,pp,smpl); 88K = zeros(mm,pp,smpl); 89L = zeros(mm,mm,smpl); 90Linf = zeros(mm,mm,smpl); 91Lstar = zeros(mm,mm,smpl); 92Kstar = zeros(mm,pp,smpl); 93Kinf = zeros(mm,pp,smpl); 94P = zeros(mm,mm,smpl+1); 95Pstar = zeros(spstar(1),spstar(2),smpl+1); 96Pstar(:,:,1) = Pstar1; 97Pinf = zeros(spinf(1),spinf(2),smpl+1); 98Pinf(:,:,1) = Pinf1; 99rr = size(Q,1); 100QQ = R*Q*transpose(R); 101QRt = Q*transpose(R); 102alphahat = zeros(mm,smpl); 103etahat = zeros(rr,smpl); 104epsilonhat = zeros(rr,smpl); 105r = zeros(mm,smpl+1); 106Finf_singular = zeros(1,smpl); 107if state_uncertainty_flag 108 V = zeros(mm,mm,smpl); 109 N = zeros(mm,mm,smpl+1); 110else 111 V=[]; 112end 113 114t = 0; 115while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl 116 t = t+1; 117 di = data_index{t}; 118 if isempty(di) 119 %no observations, propagate forward without updating based on 120 %observations 121 atilde(:,t) = a(:,t); 122 a(:,t+1) = T*atilde(:,t); 123 Linf(:,:,t) = T; 124 Pstar(:,:,t+1) = T*Pstar(:,:,t)*T' + QQ; 125 Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'; 126 else 127 ZZ = Z(di,:); %span selector matrix 128 v(di,t)= Y(di,t) - ZZ*a(:,t); %get prediction error v^(0) in (5.13) DK (2012) 129 Finf = ZZ*Pinf(:,:,t)*ZZ'; % (5.7) in DK (2012) 130 if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0 131 if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0 132 % The univariate diffuse kalman filter should be used. 133 alphahat = Inf; 134 return 135 else %rank of F_{\infty,t} is 0 136 Finf_singular(1,t) = 1; 137 Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di); % (5.7) in DK (2012) 138 if rcond(Fstar(di,di,t)) < kalman_tol %F_{*} is singular 139 if ~all(all(abs(Fstar(di,di,t))<kalman_tol)) 140 % The univariate diffuse kalman filter should be used. 141 alphahat = Inf; 142 return 143 else %rank 0 144 a(:,t+1) = T*a(:,t); 145 Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+QQ; 146 Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); 147 end 148 else 149 iFstar(di,di,t) = inv(Fstar(di,di,t)); 150 Kstar(:,di,t) = Pstar(:,:,t)*ZZ'*iFstar(di,di,t); %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0) 151 Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); % DK (2012), 5.16 152 Lstar(:,:,t) = T - T*Kstar(:,di,t)*ZZ; %L^(0) in DK (2012), eq. 5.12 153 Pstar(:,:,t+1) = T*Pstar(:,:,t)*Lstar(:,:,t)'+QQ; % (5.17) DK (2012) 154 a(:,t+1) = T*(a(:,t)+Kstar(:,di,t)*v(di,t)); % (5.13) DK (2012) 155 end 156 end 157 else 158 %see notes in kalman_filter_d.m for details of computations 159 iFinf(di,di,t) = inv(Finf); 160 Kinf(:,di,t) = Pinf(:,:,t)*ZZ'*iFinf(di,di,t); %define Kinf=T^{-1}*K_0 with M_{\infty}=Pinf*Z' 161 atilde(:,t) = a(:,t) + Kinf(:,di,t)*v(di,t); 162 Linf(:,:,t) = T - T*Kinf(:,di,t)*ZZ; %L^(0) in DK (2012), eq. 5.12 163 Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di); %(5.7) DK(2012) 164 Kstar(:,di,t) = (Pstar(:,:,t)*ZZ'-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t); %(5.12) DK(2012) with Kstar=T^{-1}*K^(1); note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf", but it is correct in their appendix 165 Pstar(:,:,t+1) = T*Pstar(:,:,t)*Linf(:,:,t)'-T*Kinf(:,di,t)*Finf*Kstar(:,di,t)'*T' + QQ; %(5.14) DK(2012) 166 Pinf(:,:,t+1) = T*Pinf(:,:,t)*Linf(:,:,t)'; %(5.14) DK(2012) 167 end 168 a(:,t+1) = T*atilde(:,t); 169 aK(1,:,t+1) = a(:,t+1); 170 % isn't a meaningless as long as we are in the diffuse part? MJ 171 for jnk=2:nk 172 aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1)); 173 end 174 end 175end 176d = t; 177P(:,:,d+1) = Pstar(:,:,d+1); 178iFinf = iFinf(:,:,1:d); 179iFstar= iFstar(:,:,1:d); 180Linf = Linf(:,:,1:d); 181Lstar = Lstar(:,:,1:d); 182Kstar = Kstar(:,:,1:d); 183Pstar = Pstar(:,:,1:d); 184Pinf = Pinf(:,:,1:d); 185notsteady = 1; 186while notsteady && t<smpl 187 t = t+1; 188 P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); % make sure P is symmetric 189 di = data_index{t}; 190 if isempty(di) 191 atilde(:,t) = a(:,t); 192 L(:,:,t) = T; 193 P(:,:,t+1) = T*P(:,:,t)*T' + QQ; %p. 111, DK(2012) 194 else 195 ZZ = Z(di,:); 196 v(di,t) = Y(di,t) - ZZ*a(:,t); 197 F = ZZ*P(:,:,t)*ZZ' + H(di,di); 198 sig=sqrt(diag(F)); 199 200 if any(diag(F)<kalman_tol) || rcond(F./(sig*sig')) < kalman_tol 201 alphahat = Inf; 202 return 203 end 204 iF(di,di,t) = inv(F./(sig*sig'))./(sig*sig'); 205 PZI = P(:,:,t)*ZZ'*iF(di,di,t); 206 atilde(:,t) = a(:,t) + PZI*v(di,t); 207 K(:,di,t) = T*PZI; 208 L(:,:,t) = T-K(:,di,t)*ZZ; 209 P(:,:,t+1) = T*P(:,:,t)*L(:,:,t)' + QQ; 210 end 211 a(:,t+1) = T*atilde(:,t); 212 Pf = P(:,:,t); 213 aK(1,:,t+1) = a(:,t+1); 214 for jnk=1:nk 215 Pf = T*Pf*T' + QQ; 216 PK(jnk,:,:,t+jnk) = Pf; 217 if jnk>1 218 aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1)); 219 end 220 end 221 % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol); 222end 223% $$$ if t<smpl 224% $$$ PZI_s = PZI; 225% $$$ K_s = K(:,:,t); 226% $$$ iF_s = iF(:,:,t); 227% $$$ P_s = P(:,:,t+1); 228% $$$ P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t])); 229% $$$ iF = cat(3,iF(:,:,1:t),repmat(iF_s,[1 1 smpl-t])); 230% $$$ L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t])); 231% $$$ K = cat(3,K(:,:,1:t),repmat(T*P_s*Z'*iF_s,[1 1 smpl-t])); 232% $$$ end 233% $$$ while t<smpl 234% $$$ t=t+1; 235% $$$ v(:,t) = Y(:,t) - Z*a(:,t); 236% $$$ atilde(:,t) = a(:,t) + PZI*v(:,t); 237% $$$ a(:,t+1) = T*atilde(:,t); 238% $$$ Pf = P(:,:,t); 239% $$$ for jnk=1:nk, 240% $$$ Pf = T*Pf*T' + QQ; 241% $$$ aK(jnk,:,t+jnk) = T^jnk*atilde(:,t); 242% $$$ PK(jnk,:,:,t+jnk) = Pf; 243% $$$ end 244% $$$ end 245%% backward pass; r_T and N_T, stored in entry (smpl+1) were initialized at 0 246t = smpl+1; 247while t>d+1 248 t = t-1; 249 di = data_index{t}; 250 if isempty(di) 251 % in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains 252 r(:,t) = L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0 253 if state_uncertainty_flag 254 N(:,:,t)=L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 with Z=0 255 end 256 else 257 ZZ = Z(di,:); 258 r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 259 if state_uncertainty_flag 260 N(:,:,t)=ZZ'*iF(di,di,t)*ZZ+L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 261 end 262 end 263 alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); %DK (2012), eq. 4.35 264 etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63 265 if state_uncertainty_flag 266 V(:,:,t) = P(:,:,t)-P(:,:,t)*N(:,:,t)*P(:,:,t); %DK (2012), eq. 4.43 267 end 268end 269if d %diffuse periods 270 % initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23 271 r0 = zeros(mm,d+1); 272 r0(:,d+1) = r(:,d+1); %set r0_{d}, i.e. shifted by one period 273 r1 = zeros(mm,d+1); %set r1_{d}, i.e. shifted by one period 274 if state_uncertainty_flag 275 %N_0 at (d+1) is N(d+1), so we can use N for continuing and storing N_0-recursion 276 N_1=zeros(mm,mm,d+1); %set N_1_{d}=0, i.e. shifted by one period, below DK (2012), eq. 5.26 277 N_2=zeros(mm,mm,d+1); %set N_2_{d}=0, i.e. shifted by one period, below DK (2012), eq. 5.26 278 end 279 for t = d:-1:1 280 di = data_index{t}; 281 if isempty(di) 282 r1(:,t) = Linf(:,:,t)'*r1(:,t+1); 283 else 284 if ~Finf_singular(1,t) 285 r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf 286 r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ... 287 + Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1) 288 if state_uncertainty_flag 289 L_1=(-T*Kstar(:,di,t)*Z(di,:)); % noting that Kstar=T^{-1}*K^(1) 290 N(:,:,t)=Linf(:,:,t)'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.19, noting that L^(0) is named Linf 291 N_1(:,:,t)=Z(di,:)'*iFinf(di,di,t)*Z(di,:)+Linf(:,:,t)'*N_1(:,:,t+1)*Linf(:,:,t)... 292 +L_1'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.29; note that, compared to DK (2003) this drops the term (L_1'*N(:,:,t+1)*Linf(:,:,t))' in the recursion due to it entering premultiplied by Pinf when computing V, and Pinf*Linf'*N=0 293 N_2(:,:,t)=Z(di,:)'*(-iFinf(di,di,t)*Fstar(di,di,t)*iFinf(di,di,t))*Z(di,:) ... 294 + Linf(:,:,t)'*N_2(:,:,t+1)*Linf(:,:,t)... 295 + Linf(:,:,t)'*N_1(:,:,t+1)*L_1... 296 + L_1'*N_1(:,:,t+1)'*Linf(:,:,t)... 297 + L_1'*N(:,:,t+1)*L_1; % DK (2012), eq. 5.29 298 end 299 else 300 r0(:,t) = Z(di,:)'*iFstar(di,di,t)*v(di,t)-Lstar(:,:,t)'*r0(:,t+1); % DK (2003), eq. (14) 301 r1(:,t) = T'*r1(:,t+1); % DK (2003), eq. (14) 302 if state_uncertainty_flag 303 N(:,:,t)=Z(di,:)'*iFstar(di,di,t)*Z(di,:)... 304 +Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) 305 N_1(:,:,t)=T'*N_1(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) 306 N_2(:,:,t)=T'*N_2(:,:,t+1)*T'; % DK (2003), eq. (14) 307 end 308 end 309 end 310 alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); % DK (2012), eq. 5.23 311 etahat(:,t) = QRt*r0(:,t); % DK (2012), p. 135 312 if state_uncertainty_flag 313 V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)... 314 -(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'... 315 - Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)... 316 - Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30 317 end 318 end 319end 320 321if decomp_flag 322 decomp = zeros(nk,mm,rr,smpl+nk); 323 ZRQinv = inv(Z*QQ*Z'); 324 for t = max(d,1):smpl 325 di = data_index{t}; 326 % calculate eta_tm1t 327 eta_tm1t = QRt*Z(di,:)'*iF(di,di,t)*v(di,t); 328 AAA = P(:,:,t)*Z(di,:)'*ZRQinv(di,di)*bsxfun(@times,Z(di,:)*R,eta_tm1t'); 329 % calculate decomposition 330 decomp(1,:,:,t+1) = AAA; 331 for h = 2:nk 332 AAA = T*AAA; 333 decomp(h,:,:,t+h) = AAA; 334 end 335 end 336end 337 338epsilonhat = Y-Z*alphahat; 339