diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 1e20e1214f49705b77d46e3e8324505dc3a62f54..5ff8ab2978804a80db003fc7ef08619687235712 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -189,6 +189,7 @@ elseif options_.lik_init == 2           % Old Diffuse Kalman filter
     Pstar = options_.Harvey_scale_factor*eye(np);
     Pinf        = [];
 elseif options_.lik_init == 3           % Diffuse Kalman filter
+    my_mf = mf;
     if kalman_algo ~= 4
         kalman_algo = 3;
     else
@@ -201,9 +202,10 @@ elseif options_.lik_init == 3           % Diffuse Kalman filter
             R  = blkdiag(R,eye(vobs));
             H   = zeros(vobs,vobs);
             Z   = [Z, eye(vobs)];
+            my_mf = find(any(Z))';
         end
     end
-    [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium);
+    [Pstar,Pinf] = compute_Pinf_Pstar(my_mf,T,R,Q,options_.qz_criterium);
 elseif options_.lik_init == 4           % Start from the solution of the Riccati equation.
     Pstar = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
     Pinf  = [];
diff --git a/matlab/compute_Pinf_Pstar.m b/matlab/compute_Pinf_Pstar.m
index 90da303d6442638cacaeff884538705a4d4dc221..de026577cf443f653c317a7068aec81fdf68a036 100644
--- a/matlab/compute_Pinf_Pstar.m
+++ b/matlab/compute_Pinf_Pstar.m
@@ -169,6 +169,7 @@ if np0
     QTinf(1:np0,np0+1:np0+nk) = STinf0;
     QTinf([indx0(:); indx(:)],:) = QTinf;
     STinf1 = [zeros(np0+np,np0) [STinf0; eye(nk); zeros(np-nk,nk)] zeros(np0+np,np-nk)];
+    mf = ismember([indx0(:); indx(:)],mf);
     for k = 1:nk
         if norm(QTinf(mf,:)*ST([indx0(:); indx(:)],k+np0)) < 1e-8
             Pinf(k+np0,k+np0) = 0;