I write the kalman's flater for this problem myself in the second part. It is much faster (3~10 time) that using SAS kalman's call routine.


data simudata;
   n=500;
   seed=1230;
   e10=rannor(seed);
   e20=rannor(seed);
   e30=rannor(seed);
   e40=rannor(seed);
   x0=0;
   do i=-200 to n;
      x=0.9*x0+2*rannor(seed);

      e1=0.95*e10+rannor(seed);
      y1=2.5+0.95*x+e1;

      e2=0.9*e20+rannor(seed);
      y2=2.0+1.1*x+e2;

       e3=0.85*e30+rannor(seed);
      y3=1.5+1.25*x+e3;

       e4=0.8*e40+rannor(seed);
      y4=2.5+0.7*x+e4;

      if i>0 then output;
      x0=x;  e10=e1; e20=e2; e30=e3; e40=e4;
   end;
   keep y1-y4 x;
run;

proc corr data=simudata;
var y1-y4 x;
run;

proc iml;

start lik(pr) global(y,z0,vz0);

   KSI0=z0`;
   P0=vz0;
   F=diag(pr[1:5]);
   A = pr[10:13]; ;
   B = pr[6:9]||I(4);
   Q = diag((pr[14:18])##2);
   R = J(4,4,0);

   NRY = nrow(y);
   sum = 0;

   DO i=1 to NRY;
      ERR= y[i, ]` - A - B*KSI0;
      OMEGA=B*P0*B` + R;
      OMEGA1=inv(OMEGA);
      K=P0*B`* OMEGA1 ;
      KSI=KSI0 + K*ERR;
      P = P0 - K*B*P0;
      KSI0=F*KSI;
      P0 = F*P*F` + Q;
      sum=sum + log(det(OMEGA)) + ERR`*OMEGA1*ERR;
    end;
    return(-1*sum);
finish;

     use simudata var{y1 y2 y3 y4};
       read all into y;

pr={ 0.9
      0.8    0.8    0.8    0.8
      0.95    0.9    0.9    0.9
      2 2 2 2
      1.5
      1    1    1    1
      };

   z0 = j(1,5,0);
   vz0= I(5);
;

   logl = lik(pr);
   print logl;
    gamma=pr[1];
    phi=pr[2:5];
    beta=pr[6:9];
    alpha=pr[10:13];
    omega=pr[14];
    sigma=pr[15:18];


    optn = {1 2 };
  
     con={ . . . . . 0.95 . . . 2.5 . . . . . . . .,
           . . . . . 0.95 . . . 2.5 . . . . . . . .
          };

  call nlpdd(rc, xr,"lik", pr , optn, con)  ;

    gamma_e=xr[1];
    phi_e=xr[2:5];
    beta_e=xr[6:9];
    alpha_e=xr[10:13];
    omega_e=xr[14];
    sigma_e=xr[15:18];

    gamma_true={0.9};
    phi_true={0.95,0.90,0.85,0.80};
    beta_true={0.95, 1.1, 1.25,0.7};
    alpha_true={2.5,2,1.5,2.5};
    omega_true={2};
    sigma_true={1,1,1,1};


    print gamma_true gamma gamma_e;
    print phi_true phi phi_e;
    print beta_true beta beta_e;
    print alpha_true alpha alpha_e;
    print omega_true omega omega_e;
    print sigma_true sigma sigma_e;
    quit;

proc iml;
start lik(pr) global(y,z0,vz0);
   
   f=diag(pr[1:5]);
   a = j(nrow(f),1,0);
   h=pr[6:9]||I(4);
   b =pr[10:13] ;
   var1=(pr[14:18]// {0.000,0.000,0.000,0.000} )##2;
   var= diag(var1);
   
    nz = nrow(f);
    n = nrow(y);
    k = ncol(y);

    *print f a h b var z0 vz0;

    const = k*log(8*atan(1));

       call kalcvf(pred,vpred,filt,vfilt,y,0,a,f,b,h,var,z0,vz0);
    et = y - pred*h` - repeat(b`,n);
    *et = y - pred*h`;
    sum1 = 0;
    sum2 = 0;
    do i = 1 to n;
       vpred_i = vpred[(i-1)*nz+1:i*nz,];
       et_i = et[i,];
       ft = h*vpred_i*h` + var[nz+1:nz+k,nz+1:nz+k];
       sum1 = sum1 + log(det(ft));
       sum2 = sum2 + et_i*inv(ft)*et_i`;
    end;
    *return(-.5*const-.5*(sum1+sum2)/n);
    return(-1*(sum1+sum2));
finish;


       use simudata var{y1 y2 y3 y4};
       read all into y;

pr={ 0.9
      0.8    0.8    0.8    0.8
      0.95    0.9    0.9    0.9
      2 2 2 2
      1.5
      1    1    1    1
      };

   z0 = j(1,5,0);
   vz0= I(5);

   logl = lik(pr);

    print logl;

      optn = {1 2 };
  
     con={ . . . . . 0.95 . . . 2.5 . . . . . . . .,
           . . . . . 0.95 . . . 2.5 . . . . . . . .
          };

  call nlpdd(rc, xr,"lik", pr , optn, con)  ;

   logl = lik(xr);
   print logl;

    quit;