From 8dc96cafa45289d1d1971ef551bc87f699371b33 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Sun, 19 Jun 2016 17:13:04 +0200
Subject: [PATCH] If F is identically 0 in Kalman filter, discard parameter
 draw instead of treating current observation as unobserved

See discussion on mailing list 18/06/2016
---
 matlab/kalman/likelihood/kalman_filter.m                     | 5 +++--
 matlab/kalman/likelihood/kalman_filter_d.m                   | 5 ++---
 matlab/kalman/likelihood/kalman_filter_fast.m                | 5 +++--
 .../kalman/likelihood/missing_observations_kalman_filter.m   | 5 +++--
 .../kalman/likelihood/missing_observations_kalman_filter_d.m | 5 ++---
 5 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/matlab/kalman/likelihood/kalman_filter.m b/matlab/kalman/likelihood/kalman_filter.m
index 14208f471..5c375a241 100644
--- a/matlab/kalman/likelihood/kalman_filter.m
+++ b/matlab/kalman/likelihood/kalman_filter.m
@@ -172,10 +172,11 @@ while notsteady && t<=last
     sig=sqrt(diag(F));
     if any(diag(F)<kalman_tol) || rcond(F./(sig*sig')) < kalman_tol
         if ~all(abs(F(:))<kalman_tol)
+            %use diffuse filter
             return
         else
-            a = T*a;
-            P = T*P*transpose(T)+QQ;
+            %pathological case, discard draw
+            return
         end
     else
         F_singular = 0;
diff --git a/matlab/kalman/likelihood/kalman_filter_d.m b/matlab/kalman/likelihood/kalman_filter_d.m
index 4a2c94cf1..a1e3281d7 100644
--- a/matlab/kalman/likelihood/kalman_filter_d.m
+++ b/matlab/kalman/likelihood/kalman_filter_d.m
@@ -81,9 +81,8 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
                     % The univariate diffuse kalman filter should be used.
                     return
                 else                                                    %rank 0
-                    a = T*a;
-                    Pstar = T*Pstar*transpose(T)+QQ;
-                    Pinf  = T*Pinf*transpose(T);                        % (5.16) DK (2012)
+                   %pathological case, discard draw
+                   return
                 end
             else
                 iFstar = inv(Fstar);
diff --git a/matlab/kalman/likelihood/kalman_filter_fast.m b/matlab/kalman/likelihood/kalman_filter_fast.m
index 6eeb90f76..a60011ade 100644
--- a/matlab/kalman/likelihood/kalman_filter_fast.m
+++ b/matlab/kalman/likelihood/kalman_filter_fast.m
@@ -152,10 +152,11 @@ while notsteady && t<=last
     end
     if rcond(F) < kalman_tol
         if ~all(abs(F(:))<kalman_tol)
+            % The univariate diffuse kalman filter should be used.            
             return
         else
-            a = T*a;
-            P = T*P*transpose(T)+QQ;
+            %pathological case, discard draw
+            return
         end
     else
         F_singular = 0;
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
index fab3bfb89..6b70a59af 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
@@ -105,10 +105,11 @@ while notsteady && t<=last
         sig=sqrt(diag(F));
         if any(diag(F)<kalman_tol) || rcond(F./(sig*sig')) < kalman_tol
             if ~all(abs(F(:))<kalman_tol)
+                %use diffuse filter
                 return
             else
-                a = T*a;
-                P = T*P*transpose(T)+QQ;
+                %pathological case, discard draw
+                return
             end
         else
             F_singular = 0;
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
index a61c3e15c..acfa5563c 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
@@ -100,9 +100,8 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
                         % The univariate diffuse kalman filter should be used.
                         return
                     else %rank 0
-                        a = T*a;
-                        Pstar = T*Pstar*transpose(T)+QQ;
-                        Pinf  = T*Pinf*transpose(T);                        % (5.16) DK (2012)
+                        %pathological case, discard draw
+                        return
                     end
                 else
                     iFstar = inv(Fstar);