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;