diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index 80b5c6efa..d5f5cde9f 100644
--- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
@@ -43,6 +43,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal
 %   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). 
+%   Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
+%   Second Edition, Ch. 5
 
 % Copyright (C) 2004-2016 Dynare Team
 %
@@ -83,16 +85,16 @@ K               = zeros(mm,pp,smpl);
 L               = zeros(mm,mm,smpl);
 Linf            = zeros(mm,mm,smpl);
 Kstar           = zeros(mm,pp,smpl);
+Kinf            = 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;
-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(rr,smpl);
+epsilonhat      = zeros(rr,smpl);
 r               = zeros(mm,smpl+1);
 
 t = 0;
@@ -100,54 +102,57 @@ while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
     t = t+1;
     di = data_index{t};
     if isempty(di)
-        atilde(:,t) = a(:,t);
+        %no observations, propagate forward without updating based on
+        %observations
+        atilde(:,t)     = a(:,t);
+        a(:,t+1)        = T*atilde(:,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);
-        Finf = ZZ*Pinf(:,:,t)*ZZ';
-        if rcond(Finf) < diffuse_kalman_tol
-            if ~all(abs(Finf(:)) < diffuse_kalman_tol)
+        ZZ = Z(di,:);                                                       %span selector matrix
+        v(di,t)= Y(di,t) - ZZ*a(:,t);                                       %get prediction error v^(0) in (5.13) DK (2012)
+        Finf = ZZ*Pinf(:,:,t)*ZZ';                                          % (5.7) in DK (2012)
+        if rcond(Finf) < diffuse_kalman_tol                                 %F_{\infty,t} = 0 
+            if ~all(abs(Finf(:)) < diffuse_kalman_tol)                      %rank-deficient but not rank 0
                 % The univariate diffuse kalman filter should be used.
                 alphahat = Inf;
                 return
-            else
-                Fstar(:,:,t)  = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);
-                if rcond(Fstar(:,:,t)) < kalman_tol
+            else                                                            %rank of F_{\infty,t} is 0
+                Fstar(:,:,t)  = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);             % (5.7) in DK (2012)
+                if rcond(Fstar(:,:,t)) < kalman_tol                         %F_{*} is singular
                     if ~all(abs(Fstar(:,:,t))<kalman_tol)
                         % The univariate diffuse kalman filter should be used.
                         alphahat = Inf;
                         return
-                    else
+                    else %rank 0
                         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)*ZZ'*iFstar(:,:,t);
-                    Pinf(:,:,t+1)   = T*Pinf(:,:,t)*transpose(T);
-                    Pstar(:,:,t+1)  = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kstar(:,:,t)')*T'+QQ;
-                    a(:,t+1)        = T*(a(:,t)+Kstar(:,:,t)*v(:,t));
+                    Kstar(:,:,t)  = Pstar(:,:,t)*ZZ'*iFstar(:,:,t);         %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0)
+                    Pinf(:,:,t+1)   = T*Pinf(:,:,t)*transpose(T);           % DK (2012), 5.16
+                    Pstar(:,:,t+1)  = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kstar(:,:,t)')*T'+QQ;    % (5.17) DK (2012) with L_0 plugged in
+                    a(:,t+1)        = T*(a(:,t)+Kstar(:,:,t)*v(:,t));       % (5.13) DK (2012)
                 end
             end
         else
+            %see notes in kalman_filter_d.m for details of computations
             iFinf(di,di,t)  = inv(Finf);
-            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' + H(di,di);
-            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)'-T*Pinf(:,:,t)*ZZ'*Kstar(:,di,t)' + QQ;
-            Pinf(:,:,t+1)   = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*ZZ'*Kinf(:,di,t)';
+            Kinf(:,di,t)    = Pinf(:,:,t)*ZZ'*iFinf(di,di,t);               %define Kinf=T^{-1}*K_0 with M_{\infty}=Pinf*Z'
+            atilde(:,t)     = a(:,t) + Kinf(:,di,t)*v(di,t);
+            Linf(:,:,t)     = T - T*Kinf(:,di,t)*ZZ;                        %L^(0) in DK (2012), eq. 5.12
+            Fstar(di,di,t)  = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);               %(5.7) DK(2012)
+            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)
+            Pstar(:,:,t+1)  = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kinf(:,di,t)'-Pinf(:,:,t)*ZZ'*Kstar(:,di,t)')*T' + QQ; %(5.14) DK(2012)
+            Pinf(:,:,t+1)   = T*(Pinf(:,:,t)-Pinf(:,:,t)*ZZ'*Kinf(:,di,t)')*T';     %(5.14) DK(2012)
         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,
+        for jnk=2:nk
             aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
         end
     end
@@ -163,12 +168,12 @@ Pinf  = Pinf(:,:,1:d);
 notsteady = 1;
 while notsteady && t<smpl
     t = t+1;
-    P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
+    P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));                   % make sure P is symmetric
     di = data_index{t};
     if isempty(di)
-        atilde(:,t) = a(:,t);
+        atilde(:,t)     = a(:,t);
         L(:,:,t)        = T;
-        P(:,:,t+1)      = T*P(:,:,t)*T' + QQ;
+        P(:,:,t+1)      = T*P(:,:,t)*T' + QQ;                               %p. 111, DK(2012)
     else
         ZZ = Z(di,:);
         v(di,t)      = Y(di,t) - ZZ*a(:,t);
@@ -220,34 +225,37 @@ end
 % $$$   PK(jnk,:,:,t+jnk) = Pf;
 % $$$     end
 % $$$ end
+%% backward pass
 t = smpl+1;
 while t>d+1
     t = t-1;
     di = data_index{t};
     if isempty(di)
+        % in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains
         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);
+        r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1);              %DK (2012), eq. 4.38
     end
-    alphahat(:,t)       = a(:,t) + P(:,:,t)*r(:,t);
-    etahat(:,t) = QRt*r(:,t);
+    alphahat(:,t)       = a(:,t) + P(:,:,t)*r(:,t);                         %DK (2012), eq. 4.35
+    etahat(:,t) = QRt*r(:,t);                                               %DK (2012), eq. 4.63
 end
-if d
+if d %diffuse periods
+    % initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23
     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);
+        r0(:,t) = Linf(:,:,t)'*r0(:,t+1);                                   % DK (2012), eq. 5.21 where L^(0) is named Linf 
         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);
+            r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ...
+                      + 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)
         end
-        alphahat(:,t)   = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
-        etahat(:,t)             = QRt*r0(:,t);
+        alphahat(:,t)   = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); % DK (2012), eq. 5.23
+        etahat(:,t)     = QRt*r0(:,t); % DK (2012), p. 135
     end
 end
 
@@ -260,7 +268,6 @@ if decomp_flag
         eta_tm1t = QRt*Z(di,:)'*iF(di,di,t)*v(di,t);
         AAA = P(:,:,t)*Z(di,:)'*ZRQinv(di,di)*bsxfun(@times,Z(di,:)*R,eta_tm1t');
         % calculate decomposition
-        Ttok = eye(mm,mm); 
         decomp(1,:,:,t+1) = AAA;
         for h = 2:nk
             AAA = T*AAA;