From bf65cb4f10ae3d123d751adfce6300c2d1ee1f30 Mon Sep 17 00:00:00 2001
From: michel <michel@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Wed, 3 Dec 2008 07:30:26 +0000
Subject: [PATCH] various bug fixes for filter and smoother with missing
 observations

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2282 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/DiffuseKalmanSmoother3.m               |  1 +
 matlab/DiffuseKalmanSmoother3_Z.m             |  1 +
 ...ssing_observations_diffuse_kalman_filter.m |  2 +-
 matlab/missing_DiffuseKalmanSmoother3.m       | 55 ++++++++--------
 matlab/missing_DiffuseKalmanSmoother3_Z.m     | 65 ++++++++++---------
 5 files changed, 64 insertions(+), 60 deletions(-)

diff --git a/matlab/DiffuseKalmanSmoother3.m b/matlab/DiffuseKalmanSmoother3.m
index 2889c372e..f3b0bc57d 100644
--- a/matlab/DiffuseKalmanSmoother3.m
+++ b/matlab/DiffuseKalmanSmoother3.m
@@ -74,6 +74,7 @@ end
 
 Fstar   	= zeros(pp,smpl_diff);
 Finf		= zeros(pp,smpl_diff);
+Fi		= zeros(pp,smpl_diff);
 Ki       	= zeros(mm,pp,smpl);
 Li      	= zeros(mm,mm,pp,smpl);
 Linf    	= zeros(mm,mm,pp,smpl_diff);
diff --git a/matlab/DiffuseKalmanSmoother3_Z.m b/matlab/DiffuseKalmanSmoother3_Z.m
index c2f58b3bd..216f72636 100644
--- a/matlab/DiffuseKalmanSmoother3_Z.m
+++ b/matlab/DiffuseKalmanSmoother3_Z.m
@@ -83,6 +83,7 @@ end
 
 Fstar           = zeros(pp,smpl_diff);
 Finf            = zeros(pp,smpl_diff);
+Fi              = zeros(pp,smpl_diff);
 Ki              = zeros(mm,pp,smpl);
 Li              = zeros(mm,mm,pp,smpl);
 Linf            = zeros(mm,mm,pp,smpl_diff);
diff --git a/matlab/kalman/likelihood/missing_observations_diffuse_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_diffuse_kalman_filter.m
index 15108fdc2..a3d74c639 100644
--- a/matlab/kalman/likelihood/missing_observations_diffuse_kalman_filter.m
+++ b/matlab/kalman/likelihood/missing_observations_diffuse_kalman_filter.m
@@ -53,7 +53,7 @@ function [LIK, lik] = missing_observations_diffuse_kalman_filter(T,R,Q,H,Pinf,Ps
   oldK = 0;
   lik  = zeros(smpl+1,1);
   LIK  = Inf;
-  lik(smpl+1) = number_of_observations*pp*log(2*pi);
+  lik(smpl+1) = number_of_observations*log(2*pi);
   notsteady   = 1;
   reste       = 0;
   
diff --git a/matlab/missing_DiffuseKalmanSmoother3.m b/matlab/missing_DiffuseKalmanSmoother3.m
index 13efff232..140992e3b 100644
--- a/matlab/missing_DiffuseKalmanSmoother3.m
+++ b/matlab/missing_DiffuseKalmanSmoother3.m
@@ -75,6 +75,7 @@ end
 
 Fstar   	= zeros(pp,smpl_diff);
 Finf		= zeros(pp,smpl_diff);
+Fi		= zeros(pp,smpl_diff);
 Ki       	= zeros(mm,pp,smpl);
 Li      	= zeros(mm,mm,pp,smpl);
 Linf    	= zeros(mm,mm,pp,smpl_diff);
@@ -211,34 +212,34 @@ while notsteady & t<smpl
     aK(jnk,:,t+jnk) 	 	= T^jnk*a(:,t);
   end
   P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
-  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
-Fi_s = Fi(:,t);
-Ki_s = Ki(:,:,t);
-L_s  =Li(:,:,:,t);
-if t<smpl
-  t_steady = t+1;
-  P  = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
-  Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
-  Li  = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
-  Ki  = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
-end
-while t<smpl
-  t=t+1;
-  a(:,t) = a1(:,t);
-  di = data_index{t}';
-  for i=di
-    v(i,t)      = Y(i,t) - a(mf(i),t) - trend(i,t);
-    if Fi_s(i) > crit
-      a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
-    end
-  end
-  a1(:,t+1) = T*a(:,t);
-  for jnk=1:nk,
-    aK(jnk,:,t+jnk)	= T^jnk*a(:,t);
-  end
+%  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
 end
+% $$$ P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
+% $$$ Fi_s = Fi(:,t);
+% $$$ Ki_s = Ki(:,:,t);
+% $$$ L_s  =Li(:,:,:,t);
+% $$$ if t<smpl
+% $$$   t_steady = t+1;
+% $$$   P  = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
+% $$$   Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
+% $$$   Li  = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
+% $$$   Ki  = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
+% $$$ end
+% $$$ while t<smpl
+% $$$   t=t+1;
+% $$$   a(:,t) = a1(:,t);
+% $$$   di = data_index{t}';
+% $$$   for i=di
+% $$$     v(i,t)      = Y(i,t) - a(mf(i),t) - trend(i,t);
+% $$$     if Fi_s(i) > crit
+% $$$       a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
+% $$$     end
+% $$$   end
+% $$$   a1(:,t+1) = T*a(:,t);
+% $$$   for jnk=1:nk,
+% $$$     aK(jnk,:,t+jnk)	= T^jnk*a(:,t);
+% $$$   end
+% $$$ end
 ri=zeros(mm,1);
 t = smpl+1;
 while t>d+1
diff --git a/matlab/missing_DiffuseKalmanSmoother3_Z.m b/matlab/missing_DiffuseKalmanSmoother3_Z.m
index e7c74a7a4..c2da9ec43 100644
--- a/matlab/missing_DiffuseKalmanSmoother3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmoother3_Z.m
@@ -84,6 +84,7 @@ end
 
 Fstar           = zeros(pp,smpl_diff);
 Finf            = zeros(pp,smpl_diff);
+Fi              = zeros(pp,smpl);
 Ki              = zeros(mm,pp,smpl);
 Li              = zeros(mm,mm,pp,smpl);
 Linf            = zeros(mm,mm,pp,smpl_diff);
@@ -214,39 +215,39 @@ while notsteady & t<smpl
       PK(jnk,:,:,t+jnk) = Pf;
   end
   P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
-  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-end
-P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
-P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
-Fi_s = Fi(:,t);
-Ki_s = Ki(:,:,t);
-L_s  =Li(:,:,:,t);
-if t<smpl
-  P  = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
-  P1  = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
-  Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
-  Li  = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
-  Ki  = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
-end
-while t<smpl
-  t=t+1;
-  a(:,t) = a1(:,t);
-  di = data_index{t}';
-  for i=di
-    Zi = Z(i,:);
-    v(i,t)      = Y(i,t) - Zi*a(:,t);
-    if Fi_s(i) > crit
-      a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
-    end
-  end
-  a1(:,t+1) = T*a(:,t);
-  Pf          = P(:,:,t);
-  for jnk=1:nk,
-      Pf = T*Pf*T' + QQ;
-      aK(jnk,:,t+jnk) = T^jnk*a(:,t);
-      PK(jnk,:,:,t+jnk) = Pf;
-  end
+%  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
 end
+% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
+% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
+% $$$ Fi_s = Fi(:,t);
+% $$$ Ki_s = Ki(:,:,t);
+% $$$ L_s  =Li(:,:,:,t);
+% $$$ if t<smpl
+% $$$   P  = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
+% $$$   P1  = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
+% $$$   Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
+% $$$   Li  = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
+% $$$   Ki  = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
+% $$$ end
+% $$$ while t<smpl
+% $$$   t=t+1;
+% $$$   a(:,t) = a1(:,t);
+% $$$   di = data_index{t}';
+% $$$   for i=di
+% $$$     Zi = Z(i,:);
+% $$$     v(i,t)      = Y(i,t) - Zi*a(:,t);
+% $$$     if Fi_s(i) > crit
+% $$$       a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
+% $$$     end
+% $$$   end
+% $$$   a1(:,t+1) = T*a(:,t);
+% $$$   Pf          = P(:,:,t);
+% $$$   for jnk=1:nk,
+% $$$       Pf = T*Pf*T' + QQ;
+% $$$       aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+% $$$       PK(jnk,:,:,t+jnk) = Pf;
+% $$$   end
+% $$$ end
 ri=zeros(mm,1);
 t = smpl+1;
 while t > d+1