/************     Kalman Filter (Nixon93)     ****************************
**                                                                      **
**          KALMAN 1.0 - Time Varying Regression Program                **
**              revised: 12/23/93                                       **
**                                                                      **
**                 David C. Nixon                                       **
**                 Department of Political Science                      **
**                 Washington University                                **
**                 St. Louis, MO 63130                                  **
**                 nixon@fcrc-next.ecs.wustl.edu                        **
**                                                                      **
** An excellent reference for the Kalman filter (and the basis          **
**  for this program both in terms of the algebra and the notation)     **
**  is:                                                                 **
**   Richard J. Meinhold and Nozer D. Singpurwalla. 1983.               **
**     "Understanding the Kalman Filter." American Statistician.        **
**     37:123-7.                                                        **
**                                                                      **
** The data in KAL.DAT is a time series of market returns with a        **
**  constant and one covariate - valuated market return.  More          **
**  complicated models (such as additional covariates or a multi-       **
**  variate model may require some program doctoring.                   **
**                                                                      **
** The program requires you to specify a starting value for             **
**  theta, the vector of parameters on the predictors                   **
**  v0, the prior for the variance matrix of the observation equation   **
**  w0, the prior for the variance matrix of the system equation        **
**                                                                      **
** The tradeoff in specifying priors is between prediction accuracy and **
**  estimation precision.  Thus, a very small prior for v0 and some-    **
**  what larger prior for w0 results in a tighter confidence interval   **
**  on theta, but a looser confidence interval for Y.  In this program  **
**  the effects of the priors are fairly constrained because the        **
**  variance matrices for v and w are estimated, rather than assuming   **
**  we know them, as Meinhold and Singpurwalla assume.                  **
**                                                                      **
** There are two output files                                           **
**     KAL.OUT is an ascii file of the parameters and standard          **
**       deviations, along with a counter                               **
**     KAL.PLT, which is generated on the screen, is a visual           **
**       for the data in kal.out.  It is an HPLG plotter file, which    **
**       can be loaded into Wordperfect (sorry, that's what I use)      **
**                                                                      **
*************************************************************************/

  dataset = "c:\\kalman\\kal.dat";
       open f1=^dataset;
          nc=colsf(f1);
          n = getnr(6,nc);
          z={};
          do until eof(f1);
             z=z|readr(f1,n);
          endo;
        closeall f1;

output file = kal.out reset;



let theta0[2,1] = 0 1;
let sigma0[2,2] = 1 0,
                  0 1;
v0 = .05;                        @ n x n @   @ specify priors here @
let w0[2,2] = .7 0               @ k x k @   @                     @
              0 .7;

T=rows(z);
w = zeros(rows(z),2);
v = zeros(rows(z),1);
theta = zeros(rows(z),2);

let wmean={0, 0};
vmean=0;
let swsq[2,2]=0;
svsq=0;

i=1;
do while i LE rows(z);
y = z[i,4];
x = z[i,3 2];
yearmo = z[i,1];

e = y - x*theta0;
r = sigma0 + w0;

theta1 = theta0 + r*x'*inv(v0 + x*r*x')*e;
sigma1 = r - r*x'*inv(v0 + x*r*x')*x*r;

/********** ESTIMATION OF VARIANCES ****************/
wt = theta1-theta0;
vt = y - x*theta1;
v[i,1] = vt;
w[i,1] = wt[1,1];
w[i,2] = wt[2,1];
wmean = wmean + (1/T)*wt;
vmean = vmean + (1/T)*vt;
j=1;
let swsq = 0 0;
svsq = 0;
do while j le i;
swsq = swsq + (w[j,.] - wmean')'*(w[j,.] - wmean');
svsq = svsq + (v[j,1] - vmean)*(v[j,1] - vmean)';
j=j+1;
endo;
wvar = swsq/i;
vvar = svsq/i;
w0 = wvar;
v0 = vvar;
/**************************************************/

sdev1 = sqrt(sigma1[1,1]);
sdev2 = sqrt(sigma1[2,2]);

format /rd 7,4;
print theta1' sdev1 sdev2 yearmo;

theta0=theta1;
sigma0=sigma1;

theta[i,.]=theta1';
i = i+1;
endo;



library pgraph;
graphset;
_pdate="";
fonts ("complex");
title("KALMAN Filter Estimates");
xlabel("Time");
ylabel("Theta");
_plctrl = 0;
_pltype = {6 6};

@ graphprt("-c=3 -cf=kal.plt -w=3");    @

xy(z[1:rows(z),1],theta);