missing_DiffuseKalmanSmootherH3_Z

 function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag)
 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
    H:        pp*1  vector of 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
    data_index                   [cell]      1*smpl cell of column vectors of indices.
    nk        number of forecasting periods
    kalman_tol   tolerance for zero divider 
    decomp_flag  if true, compute filter decomposition

 OUTPUTS
    alphahat: smoothed state variables (a_{t|T})
    epsilonhat: measurement errors
    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)
    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).

0001 function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag)
0002 % function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag)
0003 % Computes the diffuse kalman smoother without measurement error, in the case of a singular var-cov matrix.
0004 % Univariate treatment of multivariate time series.
0005 %
0006 % INPUTS
0007 %    T:        mm*mm matrix
0008 %    Z:        pp*mm matrix
0009 %    R:        mm*rr matrix
0010 %    Q:        rr*rr matrix
0011 %    H:        pp*1  vector of variance of measurement errors
0012 %    Pinf1:    mm*mm diagonal matrix with with q ones and m-q zeros
0013 %    Pstar1:   mm*mm variance-covariance matrix with stationary variables
0014 %    Y:        pp*1 vector
0015 %    pp:       number of observed variables
0016 %    mm:       number of state variables
0017 %    smpl:     sample size
0018 %    data_index                   [cell]      1*smpl cell of column vectors of indices.
0019 %    nk        number of forecasting periods
0020 %    kalman_tol   tolerance for zero divider
0021 %    decomp_flag  if true, compute filter decomposition
0022 %
0023 % OUTPUTS
0024 %    alphahat: smoothed state variables (a_{t|T})
0025 %    epsilonhat: measurement errors
0026 %    etahat:   smoothed shocks
0027 %    a:        matrix of updated variables (a_{t|t})
0028 %    aK:       3D array of k step ahead filtered state variables (a_{t+k|t})
0029 %              (meaningless for periods 1:d)
0030 %    P:        3D array of one-step ahead forecast error variance
0031 %              matrices
0032 %    PK:       4D array of k-step ahead forecast error variance
0033 %              matrices (meaningless for periods 1:d)
0034 %    decomp:   decomposition of the effect of shocks on filtered values
0035 %
0036 % SPECIAL REQUIREMENTS
0037 %   See "Filtering and Smoothing of State Vector for Diffuse State Space
0038 %   Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
0039 %   Analysis, vol. 24(1), pp. 85-98).
0040 
0041 % Copyright (C) 2004-2011 Dynare Team
0042 %
0043 % This file is part of Dynare.
0044 %
0045 % Dynare is free software: you can redistribute it and/or modify
0046 % it under the terms of the GNU General Public License as published by
0047 % the Free Software Foundation, either version 3 of the License, or
0048 % (at your option) any later version.
0049 %
0050 % Dynare is distributed in the hope that it will be useful,
0051 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0052 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0053 % GNU General Public License for more details.
0054 %
0055 % You should have received a copy of the GNU General Public License
0056 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0057 
0058 % Modified by M. Ratto
0059 % New output argument aK: 1-step to nk-stpe ahed predictions)
0060 % New input argument nk: max order of predictions in aK
0061 
0062 d = 0;
0063 decomp = [];
0064 spinf           = size(Pinf1);
0065 spstar          = size(Pstar1);
0066 v               = zeros(pp,smpl);
0067 a               = zeros(mm,smpl);
0068 a1              = zeros(mm,smpl+1);
0069 aK          = zeros(nk,mm,smpl+nk);
0070 
0071 Fstar           = zeros(pp,smpl);
0072 Finf            = zeros(pp,smpl);
0073 Fi              = zeros(pp,smpl);
0074 Ki              = zeros(mm,pp,smpl);
0075 Kstar           = zeros(mm,pp,smpl);
0076 P               = zeros(mm,mm,smpl+1);
0077 P1              = P;
0078 PK              = zeros(nk,mm,mm,smpl+nk);
0079 Pstar           = zeros(spstar(1),spstar(2),smpl); Pstar(:,:,1) = Pstar1;
0080 Pinf            = zeros(spinf(1),spinf(2),smpl); Pinf(:,:,1) = Pinf1;
0081 Pstar1          = Pstar;
0082 Pinf1           = Pinf;
0083 crit1       = 1.e-6;
0084 steady          = smpl;
0085 rr              = size(Q,1); % number of structural shocks
0086 QQ              = R*Q*transpose(R);
0087 QRt                     = Q*transpose(R);
0088 alphahat        = zeros(mm,smpl);
0089 etahat          = zeros(rr,smpl);
0090 epsilonhat      = zeros(rr,smpl);
0091 r               = zeros(mm,smpl);
0092 
0093 t = 0;
0094 icc=0;
0095 newRank   = rank(Pinf(:,:,1),crit1);
0096 while newRank && t < smpl
0097     t = t+1;
0098     a(:,t) = a1(:,t);
0099     Pstar1(:,:,t) = Pstar(:,:,t);
0100     Pinf1(:,:,t) = Pinf(:,:,t);
0101     di = data_index{t}';
0102     for i=di
0103         Zi = Z(i,:);
0104         v(i,t)      = Y(i,t)-Zi*a(:,t);
0105         Fstar(i,t)  = Zi*Pstar(:,:,t)*Zi' +H(i);
0106         Finf(i,t)   = Zi*Pinf(:,:,t)*Zi';
0107         Kstar(:,i,t) = Pstar(:,:,t)*Zi';
0108         if Finf(i,t) > kalman_tol && newRank
0109             icc=icc+1;
0110             Kinf(:,i,t)       = Pinf(:,:,t)*Zi';
0111             Kinf_Finf         = Kinf(:,i,t)/Finf(i,t);
0112             a(:,t)            = a(:,t) + Kinf_Finf*v(i,t);
0113             Pstar(:,:,t)      = Pstar(:,:,t) + ...
0114                 Kinf(:,i,t)*Kinf_Finf'*(Fstar(i,t)/Finf(i,t)) - ...
0115                 Kstar(:,i,t)*Kinf_Finf' - ...
0116                 Kinf_Finf*Kstar(:,i,t)';
0117             Pinf(:,:,t)       = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t);
0118         elseif Fstar(i,t) > kalman_tol 
0119             a(:,t)            = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
0120             Pstar(:,:,t)      = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t);
0121         end
0122     end
0123     if newRank
0124         oldRank = rank(Pinf(:,:,t),crit1);
0125     else
0126         oldRank = 0;
0127     end
0128     a1(:,t+1) = T*a(:,t);
0129     aK(1,:,t+1) = a1(:,t+1); 
0130     for jnk=2:nk
0131         aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
0132     end
0133     Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
0134     Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
0135     P0=Pinf(:,:,t+1);
0136     if newRank,
0137         newRank       = rank(Pinf(:,:,t+1),crit1);
0138     end
0139     if oldRank ~= newRank
0140         disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!')   
0141     end
0142 end
0143 
0144 
0145 d = t;
0146 P(:,:,d+1) = Pstar(:,:,d+1);
0147 Fstar = Fstar(:,1:d);
0148 Finf = Finf(:,1:d);
0149 Kstar = Kstar(:,:,1:d);
0150 Pstar = Pstar(:,:,1:d);
0151 Pinf  = Pinf(:,:,1:d);
0152 Pstar1 = Pstar1(:,:,1:d);
0153 Pinf1  = Pinf1(:,:,1:d);
0154 notsteady = 1;
0155 while notsteady && t<smpl
0156     t = t+1;
0157     a(:,t) = a1(:,t);
0158     P1(:,:,t) = P(:,:,t);
0159     di = data_index{t}';
0160     for i=di
0161         Zi = Z(i,:);
0162         v(i,t)  = Y(i,t) - Zi*a(:,t);
0163         Fi(i,t) = Zi*P(:,:,t)*Zi' + H(i);
0164         Ki(:,i,t) = P(:,:,t)*Zi';
0165         if Fi(i,t) > kalman_tol
0166             a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
0167             P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t);
0168         end
0169     end
0170     a1(:,t+1) = T*a(:,t);
0171     Pf          = P(:,:,t);
0172     aK(1,:,t+1) = a1(:,t+1); 
0173     for jnk=1:nk
0174         Pf = T*Pf*T' + QQ;
0175         PK(jnk,:,:,t+jnk) = Pf;
0176         if jnk>1
0177             aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));    
0178         end
0179     end
0180     P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
0181     %  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol);
0182 end
0183 % $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
0184 % $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
0185 % $$$ Fi_s = Fi(:,t);
0186 % $$$ Ki_s = Ki(:,:,t);
0187 % $$$ L_s  =Li(:,:,:,t);
0188 % $$$ if t<smpl
0189 % $$$   P  = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
0190 % $$$   P1  = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
0191 % $$$   Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
0192 % $$$   Li  = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
0193 % $$$   Ki  = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
0194 % $$$ end
0195 % $$$ while t<smpl
0196 % $$$   t=t+1;
0197 % $$$   a(:,t) = a1(:,t);
0198 % $$$   di = data_index{t}';
0199 % $$$   for i=di
0200 % $$$     Zi = Z(i,:);
0201 % $$$     v(i,t)      = Y(i,t) - Zi*a(:,t);
0202 % $$$     if Fi_s(i) > kalman_tol
0203 % $$$       a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
0204 % $$$     end
0205 % $$$   end
0206 % $$$   a1(:,t+1) = T*a(:,t);
0207 % $$$   Pf          = P(:,:,t);
0208 % $$$   for jnk=1:nk,
0209 % $$$       Pf = T*Pf*T' + QQ;
0210 % $$$       aK(jnk,:,t+jnk) = T^jnk*a(:,t);
0211 % $$$       PK(jnk,:,:,t+jnk) = Pf;
0212 % $$$   end
0213 % $$$ end
0214 ri=zeros(mm,1);
0215 t = smpl+1;
0216 while t > d+1
0217     t = t-1;
0218     di = flipud(data_index{t})';
0219     for i = di
0220         if Fi(i,t) > kalman_tol
0221             ri = Z(i,:)'/Fi(i,t)*v(i,t)+ri-Ki(:,i,t)'*ri/Fi(i,t)*Z(i,:)';
0222         end
0223     end
0224     r(:,t) = ri;
0225     alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
0226     etahat(:,t) = QRt*r(:,t);
0227     ri = T'*ri;
0228 end
0229 if d
0230     r0 = zeros(mm,d); 
0231     r0(:,d) = ri;
0232     r1 = zeros(mm,d);
0233     for t = d:-1:1
0234         di = flipud(data_index{t})';
0235         for i = di
0236             if Finf(i,t) > kalman_tol 
0237                 r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
0238                           (Kinf(:,i,t)'*Fstar(i,t)/Finf(i,t)-Kstar(:,i,t)')*r0(:,t)/Finf(i,t)*Z(i,:)' + ...
0239                           r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)';
0240                 r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)';
0241             elseif Fstar(i,t) > kalman_tol % step needed whe Finf == 0
0242                 r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+r0(:,t)-(Kstar(:,i,t)'*r0(:,t))/Fstar(i,t)*Z(i,:)';
0243             end
0244         end
0245         alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
0246         r(:,t)        = r0(:,t);
0247         etahat(:,t)   = QRt*r(:,t);
0248         if t > 1
0249             r0(:,t-1) = T'*r0(:,t);
0250             r1(:,t-1) = T'*r1(:,t);
0251         end
0252     end
0253 end
0254 
0255 if decomp_flag
0256     decomp = zeros(nk,mm,rr,smpl+nk);
0257     ZRQinv = inv(Z*QQ*Z');
0258     for t = max(d,1):smpl
0259         ri_d = zeros(mm,1);
0260         di = flipud(data_index{t})';
0261         for i = di
0262             if Fi(i,t) > kalman_tol
0263                 ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)';
0264             end
0265         end
0266         
0267         % calculate eta_tm1t
0268         eta_tm1t = QRt*ri_d;
0269         % calculate decomposition
0270         Ttok = eye(mm,mm); 
0271         AAA = P1(:,:,t)*Z'*ZRQinv*Z*R;
0272         for h = 1:nk
0273             BBB = Ttok*AAA;
0274             for j=1:rr
0275                 decomp(h,:,j,t+h) = eta_tm1t(j)*BBB(:,j);
0276             end
0277             Ttok = T*Ttok;
0278         end
0279     end
0280 end
0281 
0282 epsilonhat = Y - Z*alphahat;