% $Id: kalman.tex 207 2011-12-18 18:11:16Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%


\section{The Kalman filter}

The Kalman filter is a device for estimating parameters in
a class of ``time-series''-like models that are put into state-space
form. We have used the 
notation from \cite{harvey1990}, Chapter~3.
\X{state-space form}
\X{Kalman filter}
The general state-space form is an  multivariate time series
$$y_t=Z_t\alpha_t+d_t+\epsilon_t$$ 
where $Z_t$ is an $N\times m$ matrix, $d_t$ is an $N$-dimensional
vector, $y_t$ is an $N$-dimensional vector, and $\epsilon_t$ is
a set of serially uncorrelated $N$-dimensional random vectors
with mean $0$ and correlation $H_t$.
The elements of $\alpha_t$ are not observable, but are
assumed to be generated by a first-order Markov process
$$\alpha_t=T_t\alpha_{t-1}+c_t+R_t\eta_t$$
where $T_t$ is an $m\times m$ matrix, $c_t$ is an $m\times 1$
vector, $R_t$ is an $m\times g$ matrix and $\eta_t$
is a $g\times 1$ vector of serially uncorrelated 
random vectors with mean $0$ and covariance matrix $H_t$.
The specification of the state-space system is completed by
two further assumptions:
\begin{enumerate}
  \item The initial state vector $\alpha_0$ has a mean of $a_0$ and a 
variance of~$P_0$.
  \item The random vectors $\epsilon_t$ and $\eta_t$ are uncorrelated
with each other, and uncorrelated with the initial state.
\end{enumerate}

In applications of the model, many of the
parameters $Z_t$, $d_t$, $H_t$, $T_t$, $c_t$, $R_t$, and $Q_t$
may be independent of $t$, in which case we will write them
without the subscript. Also, $R$ may be the identity matrix,
in which case we will omit it.

As a simple example of such a model, consider the (2-dimensional) random walk observed with error: %fixed.
 \begin{align}
  \nonumber \alpha_t =&\alpha_{t-1}+\eta_t\\
   y_t =&\alpha_t+\epsilon_t
 \end{align}
For this model, the following parameters are fixed
$$ T=\begin{pmatrix}1 & 0\\
                    0 & 1 \end{pmatrix}
  \quad
   Z=\begin{pmatrix}1 & 0\cr
                    0 & 1 \end{pmatrix}
  \quad
   R=\begin{pmatrix}1 & 0\cr
                    0 & 1\end{pmatrix}
  \quad
  d=(0,0)
  \quad
  c=(0,0)
$$
while the covariance matrices $Q$ and $H$ are estimated. The true values of $Q$ and $H$ used in the simulation were %fixed
$$
   Q=\begin{pmatrix}1 & 0.8\cr
                   0.8 & 1 \end{pmatrix}
  \quad
   H=\begin{pmatrix}3 & -2.5\cr
                   -2.5 & 3 \end{pmatrix}
$$
and the initial value of $a$ is $(0,0)$.%fixed
 
\section{Equations for the Kalman filter}

For a moment, go back to the general state-space model.
Given $a_0$ and $P_0$, we recursively calculate the a number of quantities
via the relationships
 \begin{align}
   \nonumber  \aone &= T_t\, a_{t-1}+c_t\\[3pt]
    \nonumber   \Pone &=T_t\, P_{t-1}\, T^\prime_t +R_t\, Q_t\, R^\prime_t\\[3pt]
   \nonumber    v_t &=y_t-Z_t\, \aone-d_t\\[3pt]
    \nonumber   F_t &=Z_t\, \Pone\, Z_t^\prime+H_t\\[3pt]
    \nonumber   a_t &=\aone-\Pone\, Z_t^\prime\, F_t^{-1}\, v_t\\[3pt]
      P_t &= \Pone-\Pone\, Z_t^\prime\, F_t^{-1}\, Z_t\, \Pone
 \end{align}
The log-likelihood function for the models parameters is
given by:
$$\log L = -\tfrac{NT}{ 2}\log 2\pi -0.5 \sum_{t=1}^T \log|F_t|
        -0.5 \sum_{t=1}^T v_t F_t^{-1} v_t
$$

The \textsc{tpl} file for the random walk Kalman filter code
follows:
\begin{lstlisting}
DATA_SECTION
  init_int nt
  init_int N
  init_int m
  int m1
  init_matrix Y(1,nt,1,N)
  matrix P0(1,m,1,m)
 !! P0.initialize(); 
 !! m1=m*(m+1)/2; 
PARAMETER_SECTION
  init_bounded_vector Qcoff(1,m1,-10.,10.1)
  init_bounded_vector Hcoff(1,m1,-10.,10.1)
  init_vector a0(1,m)
  matrix T(1,m,1,m)
  matrix TT(1,m,1,m)
  vector d(1,N)
  vector c(1,m)
  matrix chQ(1,m,1,m)
  sdreport_matrix Q(1,m,1,m)
  matrix chH(1,N,1,N)
  sdreport_matrix H(1,N,1,N)
  matrix Z(1,N,1,m)
  matrix TZ(1,m,1,N)
  objective_function_value f
 LOCAL_CALCS
   d.initialize();
   c.initialize();
   Z.initialize();
   Z(1,1)=1; Z(2,2)=1;
   T.initialize();
   T(1,1)=1; T(2,2)=1;
   TZ=trans(Z);
   TT=trans(T);
PROCEDURE_SECTION
  setup_Q();
  setup_H();
  f+=kalman_filter();
  cout << " f = " << f << endl;

FUNCTION setup_Q
  chQ.initialize();
  int ii=1;
  for (int i=1;i<=m;i++)
    for (int j=1;j<=i;j++)
      chQ(i,j)=Qcoff(ii++);
  Q=chQ*trans(chQ);    
FUNCTION setup_H
  chH.initialize();
  int ii=1;
  for (int i=1;i<=N;i++)
    for (int j=1;j<=i;j++)
      chH(i,j)=Hcoff(ii++);
  H=chH*trans(chH);    

FUNCTION dvariable kalman_filter(void)  
  dvar3_array P(0,nt,1,m,1,m);
  dvar3_array P1(1,nt,1,m,1,m);
  dvar3_array F(1,nt,1,N,1,N);
  dvar3_array Finv(1,nt,1,N,1,N);
  dvar_matrix Ptemp(1,m,1,m);
  dvar_matrix a(0,nt,1,m);
  dvar_matrix a1(1,nt,1,m);
  dvar_matrix v(1,nt,1,N);
  a(0)=a0;
  P(0)=P0;
  // This is the Kalman filter recursion. The objects tmp1
  // and tmp2 hold common calculations to optimize a bit
  int t;
  for (t=1;t<=nt;t++)
  {
    a1(t)=T*a(t-1)+c;
    P1(t)=T*P(t-1)*TT+Q;
    dvar_vector pred_y=Z*a1(t)+d;
    v(t)=Y(t)-pred_y;
    dvar_matrix tmp1=P1(t)*TZ;
    F(t)=Z*tmp1+H;
    Finv(t)=inv(F(t));
    dvar_matrix tmp2= tmp1*Finv(t);
    P(t)=P1(t)-tmp2*Z*P1(t);
    a(t)=a1(t)+tmp2*v(t);
  }
  int sgn=0;
  dvariable f=0.0;
  for (t=1;t<=nt;t++)
    f+=ln_det(F(t),sgn)+v(t)*Finv(t)*v(t);
  return f;  
TOP_OF_MAIN_SECTION
  arrmblsize=20000000;
  gradient_structure::set_CMPDIF_BUFFER_SIZE(3000000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1000000);
\end{lstlisting} 
This example was deliberately not optimized as much as it could be,
in order to retain the flavor of the more general state-space
problem. For example, since $T$ is the identity matrix and $c$ is the
zero vector, the line of code 
\begin{lstlisting}
  a1(t)=T*a(t-1)+c;
\end{lstlisting} 
reduces to
\begin{lstlisting}
  a1(t)=a(t-1);
\end{lstlisting} 
The parameters being estimated are  $a_0$, $Q$, and~$H$. 

To parameterize the covariance matrices, the Choleski decomposition
parameterization was used. This ensures that the covariance matrices are
positive (semi-) definite. The technique can be seen in the function
\texttt{setup\_Q}. 
The lower triangular matrix \texttt{ch\_Q} is filled with parameters
from a bounded vector:
\begin{lstlisting}
FUNCTION setup_Q
  chQ.initialize();
  int ii=1;
  for (int i=1;i<=m;i++)
    for (int j=1;j<=i;j++)
      chQ(i,j)=Qcoff(ii++);
  Q=chQ*trans(chQ);   // chQ is the choleski decomposition of Q 
\end{lstlisting}
Notice that the bounded vector \texttt{Qcoff} has slightly asymmetric
bounds. This is a simple way to ensure that its initial value is not
identically zero, which would lead to a singular covariance matrix.
\begin{lstlisting}
  init_bounded_vector Qcoff(1,m1,-10.,10.1)
\end{lstlisting}

The model parameters, standard deviations, and 
correlations are reproduced from the
standard \textsc{admb} report.

\begin{smallcode}
 index   name    value      std.dev    7      8      9     10     11     12     13     14     15      16   
    7   a0    -1.1682e+00 9.0191e-01  1.000
    8   a0     1.2218e+00 8.6442e-01  0.352  1.000
    9   Q      9.9468e-01 1.0862e-01  0.059 -0.006  1.000
   10   Q      7.8808e-01 7.8737e-02  0.038  0.028  0.683  1.000
   11   Q      7.8808e-01 7.8737e-02  0.038  0.028  0.683  1.000  1.000
   12   Q      8.7279e-01 9.6118e-02 -0.018  0.069  0.185  0.721  0.721  1.000
   13   H      3.1352e+00 1.8123e-01 -0.015 -0.007 -0.305 -0.136 -0.136 -0.018  1.000
   14   H     -2.7119e+00 1.4922e-01 -0.021  0.001 -0.102 -0.238 -0.238 -0.139 -0.692  1.000
   15   H     -2.7119e+00 1.4922e-01 -0.021  0.001 -0.102 -0.238 -0.238 -0.139 -0.692  1.000  1.000
   16   H      3.2264e+00 1.7936e-01  0.015 -0.029 -0.031 -0.121 -0.121 -0.249  0.370 -0.698 -0.698  1.000
\end{smallcode}


\section{Parameterizing the covariance\br matrix parameterizations}

The Choleski decomposition parameterization merely ensures that the 
matrix is positive semi-definite. By adding a small positive number to the
diagonal elements, one can ensure that the covariance matrix is
positive definite, and can speed up and improve the stability of the
estimation. Of course, what is meant by ``small'' will depend on the
particular problem being considered. 
A modified form of the routine \texttt{setup\_Q}
follows:
\begin{lstlisting}
FUNCTION setup_Q
  int i;
  chQ.initialize();
  int ii=1;
  for (i=1;i<=m;i++)
    for (int j=1;j<=i;j++)
      chQ(i,j)=Qcoff(ii++);
  Q=chQ*trans(chQ);   // chQ is the choleski decomposition of Q 
  for (i=1;i<=m;i++)
    Q(i,i)+=0.1;   // make Q positive definite
\end{lstlisting}
Performing this modification for the present model for both $Q$ and
$H$ causes the program to converge about twice as fast.

