From 6ea053e47d4542b02bc75e14c1d008c639d2340b Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@ens.fr>
Date: Thu, 13 Jan 2011 11:38:12 +0100
Subject: [PATCH] correcting bugs in diffuse Kalman smoother

---
 matlab/DiffuseKalmanSmoother1_Z.m  |  5 ++-
 matlab/DiffuseKalmanSmoother3.m    | 56 +++++++++-------------------
 matlab/DiffuseKalmanSmootherH1_Z.m | 59 +++++++++++++++++++++---------
 3 files changed, 61 insertions(+), 59 deletions(-)

diff --git a/matlab/DiffuseKalmanSmoother1_Z.m b/matlab/DiffuseKalmanSmoother1_Z.m
index 854d9b446..35fc69b30 100644
--- a/matlab/DiffuseKalmanSmoother1_Z.m
+++ b/matlab/DiffuseKalmanSmoother1_Z.m
@@ -101,7 +101,7 @@ while rank(Pinf(:,:,t+1),kalman_tol) && (t<smpl)
             % The univariate diffuse kalman filter should be used.
             return
         else
-            Fstar(:,:,t)  = Z*Pstar(:,:,t)*Z' + H;
+            Fstar(:,:,t)  = Z*Pstar(:,:,t)*Z';
             if rcond(Fstar(:,:,t)) < kalman_tol
                 if ~all(abs(Fstar(:,:,t))<kalman_tol)
                     % The univariate diffuse kalman filter should be used.
@@ -133,7 +133,8 @@ while rank(Pinf(:,:,t+1),kalman_tol) && (t<smpl)
         Linf(:,:,t)         = T - Kinf(:,:,t)*Z;
         Fstar(:,:,t)        = Z*Pstar(:,:,t)*Z';
         Kstar(:,:,t)        = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
-        Pstar(:,:,t+1)      = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*Fstar(:,:,t)*Kstar(:,:,t)' + QQ;
+        %        Pstar(:,:,t+1)      = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*Fstar(:,:,t)*Kstar(:,:,t)' + QQ;
+        Pstar(:,:,t+1)      = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-T*Pinf(:,:,t)*Z'*Kstar(:,:,t)' + QQ;
         Pinf(:,:,t+1)       = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
     end
 
diff --git a/matlab/DiffuseKalmanSmoother3.m b/matlab/DiffuseKalmanSmoother3.m
index 8675118c2..03d64d676 100644
--- a/matlab/DiffuseKalmanSmoother3.m
+++ b/matlab/DiffuseKalmanSmoother3.m
@@ -68,25 +68,18 @@ a               = zeros(mm,smpl);
 a1              = zeros(mm,smpl+1);
 aK              = zeros(nk,mm,smpl+nk);
 PK              = zeros(nk,mm,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);
-Fi              = zeros(pp,smpl_diff);
+Fstar           = zeros(pp,smpl);
+Finf            = zeros(pp,smpl);
+Fi              = zeros(pp,smpl);
 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);
+Linf            = zeros(mm,mm,pp,smpl);
+L0              = zeros(mm,mm,pp,smpl);
+Kstar           = zeros(mm,pp,smpl);
 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;
+Pstar           = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
+Pinf            = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
 Pstar1          = Pstar;
 Pinf1           = Pinf;
 crit            = options_.kalman_tol;
@@ -132,29 +125,7 @@ while newRank & t < smpl
             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 = (icc<options_.diffuse_d);  
-                %if newRank & any(diag(P0(mf,mf))>crit)==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
@@ -170,13 +141,20 @@ while newRank & t < smpl
     for jnk=2:nk
         aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
     end
+    if newRank
+        oldRank = rank(P0,crit);
+    else
+        oldRank = 0;
+    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);
+    if newRank
         newRank   = rank(P0,crit1);
     end
+    if oldRank ~= newRank
+        disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!')   
+    end
 end
 
 
diff --git a/matlab/DiffuseKalmanSmootherH1_Z.m b/matlab/DiffuseKalmanSmootherH1_Z.m
index 63c4f6c8e..7f60c7e4c 100644
--- a/matlab/DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/DiffuseKalmanSmootherH1_Z.m
@@ -97,25 +97,48 @@ t = 0;
 while rank(Pinf(:,:,t+1),crit1) & t<smpl
     t = t+1;
     v(:,t)= Y(:,t) - Z*a(:,t);
-    F = Z*Pinf(:,:,t)*Z';
-    if rcond(F) < crit
-        return          
+    Finf = Z*Pinf(:,:,t)*Z';
+    if rcond(Finf) < kalman_tol
+        if ~all(abs(Finf(:)) < kalman_tol)
+            % The univariate diffuse kalman filter should be used.
+            return
+        else
+            Fstar(:,:,t)  = Z*Pstar(:,:,t)*Z' + H;
+            if rcond(Fstar(:,:,t)) < kalman_tol
+                if ~all(abs(Fstar(:,:,t))<kalman_tol)
+                    % The univariate diffuse kalman filter should be used.
+                    return
+                else
+                    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)*Z'*iFstar(:,:,t);
+                Pinf(:,:,t+1)   = T*Pinf(:,:,t)*transpose(T);
+                Pstar(:,:,t+1)  = T*(Pstar(:,:,t)-Pstar(:,:,t)*Z'*Kstar(:,:,t)')*T'+QQ;
+                a(:,:,t+1)        = T*(a(:,:,t)+Kstar(:,:,t)*v(:,t));
+            end
+        end
+    else
+
+        iFinf(:,:,t)        = inv(F);
+        PZI                 = Pinf(:,:,t)*Z'*iFinf(:,:,t);
+        atilde(:,t)         = a(:,t) + PZI*v(:,t);
+        Kinf(:,:,t)         = T*PZI;
+        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
+            aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
+        end
+        Linf(:,:,t)         = T - Kinf(:,:,t)*Z;
+        Fstar(:,:,t)        = Z*Pstar(:,:,t)*Z' + H;
+        Kstar(:,:,t)        = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
+        Pstar(:,:,t+1)      = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-T*Pinf(:,:,t)*Z'*Kstar(:,:,t)' + QQ;
+        Pinf(:,:,t+1)       = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
     end
-    iFinf(:,:,t)        = inv(F);
-    PZI                 = Pinf(:,:,t)*Z'*iFinf(:,:,t);
-    atilde(:,t)         = a(:,t) + PZI*v(:,t);
-    Kinf(:,:,t)         = T*PZI;
-    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
-        aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
-    end
-    Linf(:,:,t)         = T - Kinf(:,:,t)*Z;
-    Fstar(:,:,t)        = Z*Pstar(:,:,t)*Z' + H;
-    Kstar(:,:,t)        = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
-    Pstar(:,:,t+1)      = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*F*Kstar(:,:,t)' + QQ;
-    Pinf(:,:,t+1)       = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
 end
 d = t;
 P(:,:,d+1) = Pstar(:,:,d+1);