options crt;
name kalhp 'Kalman filter hyperparameter estimation with ML PROC';
?
?  This example shows how to use ML PROC to estimate the variance 
?  parameters of the Kalman Filter transition equation. The model 
?  is the standard Kalman Filter model (given in the Reference 
?  Manual) with n=T, the number of observations, and m=2 (2 beta
?  parameters). The matrix Q is assumed to be diagonal, but with 
?  unknown variances. 
?
?  Example by Clint Cummins (version 4.4, Jan. 1997)
?    
?-------------------------
? generate artificial data

const T,1000;
smpl 1,T;
random(seedin=974321) x eps eta1 eta2;
x = 2*x;  ? x ~ N(0,4)

? transition equation
eta1 = eta1*sqrt(2);   ? Q(1,1) = 2
eta2 = eta2*sqrt(5);   ? Q(2,2) = 5
? initial conditions
b1 = 1.5;
b2 = 2;
smpl 2,T;
b1 = b1(-1) + eta1;
b2 = b2(-1) + eta2;

? measurement equation
smpl 1,T;
y = b1 + b2*x + eps;

? end of data generation
?-----------------------

? The following code estimates the Kalman Filter model with an 
? unknown diagonal transition matrix, using the data generated 
? above. 

? First estimate the Kalman filter with Q assumed to be an identity 
? matrix (the default but not correct; note the poor log likelihood 
? compared to final estimate, for example).

kalman y c x;
rename @logl loglqid; ? Save log likelihood value for Q=identity matrix.

? Kalman filter with known true Q; this works for simulated data, but 
? of course not for real data, in general. 

read(nrow=2,type=diag) q;
2 5;
kalman(vtrans=q) y c x;
rename @logl loglqtru;  ? Save log likelihood value for true Q. 

? Estimate the diagonal of Q (using ML and the Proc kfq, defined below).

param q11,1 q22,4;
ml kfq q11 q22;  ? ML Proc feature in TSP (August 1996)

? Rerun the Kalman filter with the estimated Q; standard error 
? estimates for final Kalman filter will not be fully consistent, 
? given the two-stage estimation method.

set q(1,1) = q11;
set q(2,2) = q22;
kalman(vtrans=q) y c x;
rename @logl loglqest;
print loglqest loglqtru loglqid;

?--------------------------------------------------------------------
Proc kfq;  ? evaluate log likelihood as a function of q parameters. 

? check parameter constraints (variances>0).

  if (q11>0 & q22>0); then; do; 
    set q(1,1) = q11;
    set q(2,2) = q22;
    supres @coef;   ? Turn off all printing for logl evaluation.
    kalman(silent,vtrans=q) y c x;
    nosupres @coef;
  enddo ;

  else ; set @logl = @miss;

endproc kfq ;