\documentclass[a4paper,11pt]{article}

\usepackage{verbatim}
\begin{document}

\section*{RFC: Implementation in \texttt{gretl} of the Kalman filter
  at the user level}

The Kalman filter is one of those functionalities that may come in
handy in a number of cases: of course we already use internally it for
ARMA estimation, but for example, it'd be nice to have a tool like
that for people who do DSGE models or maybe for people who would like
to implement ARMA estimation with series with missing values and so
on.

\subsection*{Notation}

Notation \emph{is} a problem. It seems that in econometrics everyone
is happy with $y = X \beta + u$, but we can't, as a community, make up
our minds on a standard for state-space models. Harvey (1991),
Hamilton (1994) and Harvey \& Proietti (2005) all use different
conventions.

After much hesitation, I'll settle on the one DSG Pollock uses in his
1999 book, with slight variations. A state-space model can be written
down as
\begin{eqnarray}
  \label{eq:obseq}
  \mathbf{y}_t & = & \mathbf{H}_t \xi_t + \mathbf{d}_t + \eta_t \\
  \label{eq:straneq}
  \xi_{t} & = & \mathbf{\Phi}_{t} \xi_{t-1} + \mathbf{c}_t + \nu_t
\end{eqnarray}
The symbol $T$ (not in bold) indicates the sample size. The
observables $\mathbf{y}_t$ are an $n$-vector (in many cases, $n=1$,
but not necessarily), while the unobservable vector of states
$\xi_t$ has $m$ elements.%; usually, if not always, $n < m$.
The covariance matrices for $\eta_t$ and $\nu_t$ are $\mathbf{\Omega}_t$ and
$\mathbf{\Psi}_t$, respectively.

In the special case when $\mathbf{\Phi}_t = \mathbf{\Phi}$, $\mathbf{H}_t =
\mathbf{H}$ and so on, the model will be called \emph{time-invariant}.

We have two options: we may go forward in time and do the
\emph{filtering} proper, or after that we go back and do the
\emph{smoothing}. 

To quote Pollock (p. 241), 
\begin{quote}
[\ldots] it is assumed that prior information on the previous state vector
$\xi_0$ is available in the form of an unbiased estimate $\mathbf{x}_0$ which has been drawn from
a distribution with a mean of $\xi_0$ and a dispersion matrix of $\mathbf{P}_0$.
\end{quote}
It may happen that $\mathbf{P}_{0} \to \infty$ for diffuse priors. In
this case, we need to check Koopman (1997); in the standard case,
however, the prediction equations are:
\begin{eqnarray}
  \label{eq:filtstates}
  E_{t-1}(\xi_t) = \mathbf{x}_{t|t-1} & = & 
  \mathbf{\Phi}_t \mathbf{x}_{t-1} + \mathbf{c}_t  \label{eq:statepred}\\
  V_{t-1}(\xi_t) = \mathbf{P}_{t|t-1} & = &
  \mathbf{\Phi}_t \mathbf{P}_{t-1} \mathbf{\Phi}'_t + \mathbf{\Psi}_t \label{eq:Ppred}\\
  E_{t-1}(\mathbf{y}_t) = \hat{\mathbf{y}}_{t|t-1} & = & \mathbf{H}_t \mathbf{x}_{t|t-1} + \mathbf{d}_t  \\
  \mathbf{e}_t & = & \mathbf{y}_t - \hat{\mathbf{y}}_{t|t-1} \\
  V_{t-1}(\mathbf{e}_t) & = & \mathbf{F}_t = \mathbf{H}_t \mathbf{P}_{t|t-1} \mathbf{H}_t' + \mathbf{\Omega}_t
 \end{eqnarray}
so the updating equations are:
\begin{eqnarray}
  \mathbf{x}_t & = & 
  \mathbf{x}_{t|t-1} + 
  \mathbf{K}_{t} \mathbf{e}_t \label{eq:stateupd}\\
  \mathbf{P}_{t} &  = & 
  \mathbf{P}_{t|t-1} - 
  \mathbf{P}_{t|t-1} \mathbf{H}_t' \mathbf{F}^{-1}_t \mathbf{H}_t \mathbf{P}_{t|t-1} \label{eq:Pupd}
\end{eqnarray}
where $\mathbf{K}_{t} = \mathbf{P}_{t|t-1} \mathbf{H}_t'
\mathbf{\Omega}_t^{-1}$ is the Kalman gain.  Equations
(\ref{eq:statepred})--(\ref{eq:stateupd}) and
(\ref{eq:Ppred})--(\ref{eq:Pupd}) can be written as a single set of
recursions going directly from $\mathbf{x}_{t-1}$ to $\mathbf{x}_t$
and from $\mathbf{P}_{t-1}$ to $\mathbf{P}_{t}$.

Note that in practice we may
want to implement something computationally smarter, like the
information filter (some bibliographic search is needed here, I can't
remember exactly what it looks like).

I can't be bothered right now to write the smoothing recursions, but
it's not terribly important.

\subsection*{User syntax}

The idea is to define a command block like we do with \texttt{mle} or
\texttt{gmm}:
\begin{verbatim}
list y = x1 x2
# or possibly
matrix Y = { x1 x2 }

kalman [ --smooth ]
  obsy y
  obsmat H
  obsex d
  obsvar Omega
  strmat Phi
  strex c
  strvar Psi
  inistate x0
  iniprec P0
end kalman
\end{verbatim}
where \texttt{obsy} stands for ``observable $\mathbf{y}$'',
\texttt{obsmat} stands for ``matrix in the observation equation''
(\ref{eq:obseq}), \texttt{obsex} for ``exogenous term in the
observation equation'' and \texttt{obsvar} for ``variance of the
observation equation''. A similar conventions would hold for the state
transition equation (\ref{eq:straneq}), with \texttt{str} replacing
\texttt{obs}. The strings \texttt{inistate} and \texttt{iniprec}
should be self-explanatory.

The elements \texttt{obsex}, \texttt{strex} and \texttt{inistate} need
not be compulsory. If absent, they are understood to be 0. More or
less the same may go for \texttt{iniprec}, if the initialisation is
done via a diffuse prior (the simplest thing would be to use somehing
like $P_0 = 10^{30}\cdot I$, but smarter methods are available --
again, Koopman (1997) should come to the rescue).

For time-invariant models, the arguments to \texttt{obsmat, strmat}
etcetera must be matrices of the appropriate size. In the case of
time-varying models, we have two options:
\begin{enumerate}
\item pass matrices with $T$ rows and an appropriate number of
  columns, where each row is the vectorization of the corresponding
  matrix at time $t$; for example, if $\mathbf{\Phi}_t$ is time-varying,
  then the argument to \texttt{tmat} will have to be a $T \times m^2$
  matrix. It'll be gretl's job to reshape the pertinent row of such a matrix.
\item pass matrices of the ``correct'' size, but allow them to be the
  result of user-defined functions. This is clearly more general, as
  it's easy to encapsulate the previous option into this one. However,
  several interesting possibilities are opened here, such as allowing,
  say, the elements of $\mathbf{H}_t$ to depend on $\mathbf{e}_{t-1},
  \mathbf{y}_{t-1}$, structural parameters (as can be the case with
  DSGE models) and so on. In this case, a mechanism would be needed
  for keeping track of the time index (maybe giving the \texttt{t}
  variable a special meaning in this case? defining a new dollar
  accessor \texttt{\$t} only valid inside Kalman blocks?).
\end{enumerate}

Of course, one also needs to retrieve the interesting quantities: this
could be done via via the existing accessors \texttt{\$uhat} for
$\mathbf{e}_t$ and \texttt{\$sigma} (other candidates may be
\texttt{\$vcv} and \texttt{\$h} --- but I personally like
\texttt{\$sigma} better) for $P_t$. When $n = 1$, these would be
series; otherwise \texttt{\$uhat} would be a $T \times n$ matrix (a
bit like we do for single equations/multiple equation models), while
\texttt{\$sigma} should be a $T \times \frac{n (n+1)}{2}$ matrix, with
the $t$-th row holding $\mathrm{vech}(\mathbf{F}_t)'$. New accessors
should be devoted to returning the states and their variances:
suitable candidates could be and \texttt{\$states} for $\xi_t$ and
\texttt{\$stvars} for $\mathrm{vech}(\mathbf{P}_t)'$.


\subsection*{References}

Hamilton, James D. (1994), \emph{Time Series Analysis}, Princeton
University Press.
\smallskip

\noindent
Harvey, A.C. (1991), \emph{Forecasting, Structural Time Series Models
  and the Kalman Filter}, Cambridge University Press.
\smallskip

\noindent
Harvey, A.C. and T. Proietti (2005), \emph{Readings in Unobserved
  Component Models}, Oxford University Press.
\smallskip

\noindent
Koopman, S. J. (1997), ``Exact Initial Kalman Filtering and Smoothing
for Nonstationary Time Series Models'', Journal of the American
Statistical Association, Vol. 92, No. 440, pp. 1630--1638 \smallskip

\noindent
Pollock, D.S.G. (1999), \emph{A Handbook of Time-Series Analysis,
  Signal Processing and Dynamics}, Academic Press.


\end{document}


Andrew C. Harvey 

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: