Kalman Filter: Diffuse Kalman Filtering
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: KALEX33 */
/* TITLE: Kalman Filter: Diffuse Kalman Filtering */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: */
/* DATA: */
/* */
/* SUPPORT: GJW UPDATE: APRIL 1996 */
/* REF: */
/* MISC: */
/****************************************************************/
title 'Diffuse Kalman Filtering';
proc iml;
z_1 = 0; z_2 = 0;
do i = 1 to 30;
z = 1.5*z_1 - .5*z_2 + rannor(1234567);
z_2 = z_1;
z_1 = z;
x = z + .8*rannor(1234578);
if ( i > 10 ) then y = y // x;
end;
run;
t = nrow(y);
h = { 1 0 };
f = { 1.5 -.5, 1 0 };
rt = .64;
vt = diag({1 0});
ny = nrow(h);
nz = ncol(h);
nb = nz;
nd = nz;
a = j(nz,1,0);
b = j(ny,1,0);
int = j(ny+nz,nb,0);
coef = f // h;
var = ( vt || j(nz,ny,0) ) //
( j(ny,nz,0) || rt );
intd = j(nz+nb,1,0);
coefd = i(nz) // j(nb,nd,0);
at = j(t*nz,nd+1,0);
mt = j(t*nz,nz,0);
qt = j(t*(nd+1),nd+1,0);
n0 = -1;
call kaldff(kaldff_p,dvpred,initial,s2,y,0,int,
coef,var,intd,coefd,n0,at,mt,qt);
call kalcvf(kalcvf_p,vpred,filt,vfilt,y,0,a,f,b,h,var);
print kalcvf_p kaldff_p;
d = 0;
do i = 1 to t;
dt = h*mt[(i-1)*nz+1:i*nz,]*h` + rt;
d = d + log(det(dt));
end;
s = qt[(t-1)*(nd+1)+1:t*(nd+1)-1,1:nd];
likl = -(t*log(s2) + d)/2;
likl_dff = likl - log(det(s))/2;
print 'Log L ' likl,
'Diffuse Log L' likl_dff;
quit;
title; title2;