# Maximum Likelihood Estimation for State Space Models using BFGS

I have written some code that can do Kalman filtering (using a number of different Kalman-type filters [Information Filter et al.]) for Linear Gaussian State Space Analysis for an n-dimensional state vector.

Let us consider the linear Gaussian State Space model

$$y_t = \mathbf{Z}_{t}\alpha_{t} + \epsilon_{t},$$ $$\alpha_{t + 1} = \mathbf{T}_{t}\alpha_{t} + \mathbf{R}_{t}\eta_{t},$$

where $y_{t}$ is our observation vector, $\alpha_{t}$ our state vector at time step $t$. The quantities in bold are the transformation matrices of the state space model which are set according to the characteristics of the system under consideration. We also have

$$\epsilon_{t} \sim NID(0, \mathbf{H}_{t}),$$ $$\eta_{t} \sim NID(0, \mathbf{Q}_{t}),$$ $$\alpha_{1} \sim NID(a_{1}, \mathbf{P}_{1}).$$

where $t = 1,\ldots, n$. Now, I have derived and implemented the recursion for the Kalman Filter for this generic state space model by guessing the initial parameters and variance matrices $\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$.

where the points are the Nile River water levels for Jan over 100 years, the line is the Kalamn Estimated state, and the dashed lines are the 90% confidence levels.

Now, in order to estimate the initial parameters of my model ($\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$) I have adopted the Expectation Maximisation method to maximise the loglikelihood function

$$\log L(Y_{n}|\psi) = -\frac{np}{2}\log(2\pi) - \frac{1}{2}\sum^{n}_{t = 1}(log|\mathbf{F}_{t}| + v^{\mathsf{T}}_{t}\mathbf{F}_{t}^{-1}v_{t}), \;\;\;\; (1)$$

with respect to parameter vector $\psi$. Where $v_{t}$ is the state error and $\mathbf{F}_{t}$ is the state error variance. The vector $\psi$ contains the unknown parameters (i.e., either the measurement and state disturbance variances themselves, or a re-parametrisation of these same variances, see below)

A crucial role in the maximisation of (1) is played by the gradient or score vector

$$\partial_{1}(\psi) = \partial \log L(Y_{n}|\psi)/\partial\psi,\;\;\;\; (2)$$

At least two algorithms can be used to maximise (1). The first is the so-called EM (Expectation-Maximisation) algorithm, and the second is the BFGS (Broyden-Fletcher-Goldfarb-Shannon) algorithm. I have successfully implemented the EM algorithm and this is providing correct results but is very slow; in light of this I now wish to implement the BFGS algorithm as this has better performance.

So, from (2) using Taylor's theorem we can write

$$\partial_{1}(\psi) \simeq \tilde{\partial}_{1}(\psi) + \tilde{\partial}_{2}(\psi)(\psi - \tilde{\psi}),$$

where $\partial_{2}(\psi)$ is the Hessian matrix or second derivative of the loglikelihood function wrt $\psi$ and $\tilde{\partial}_{N}(\psi) = \partial_{N}(\psi)|_{\psi = \tilde{\psi}}$. So now we can write our updated parameter vector as

$$\bar{\psi} = \tilde{\psi} - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi) = \tilde{\psi} - s\tilde{\pi(\psi)}.\;\;\;\; (3)$$

where here $\tilde{\pi(\psi)} = - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi)$ and $s$ is our step size. So the basic outline of the BFGS algorithm to get our parameter vector and maximize (1) using BFGS is:

1. Initialise parameter vector $\psi = \psi^{*}$.

2. Then apply the Kalman and disturbance smoothing filters and thus for score vector (2) at $\psi = \psi^{*}$. When running the Kalman Filter evaluate the log-likelihood function.

3. Use (3) to obtain new values for $\psi$ given by $\psi^{+}$ and go back to step 2 until the value of loglikelihood function (1) no longer improves.

This seems fine, but my problem is how to calculate the step size $s$ in (3) above. I know a line search is required, but remember $\psi$ is a vector and it is not clear what function I can use in a traditional line search technique in order to evaluate $s$.

Question: what line search method should I use in order to evaluate $s$?

Ps. Assume that I would like to use the zoom line search algorithm (algorithm 3.5 of Numerical Optimization by Nocedal & Wright).

Further to the very helpful answer by @LKlevin, I want to classify my problems more explicitly...

$^{*}$Most (L-)BFGS methods accept vector functions, that is a function $f$ that depends on a n-dimensional (in this case let us assume two-dimensions) parameter vector $\mathbf{x} = (x_{1}, x_{2})^{\mathsf{T}}$, so if we had $F(x_{1}, x_{2}) = 2x_{2}^{2} + 4x_{2}$ you would pass into the BFGS algo. $F$ and $\nabla{F}$. This is fine and the text I reference as a PDF above tells me all I need to know here in terms of using BFGS to minimise this function.

However, in the case I have above, from (3) we have

$$\bar{\psi} = \tilde{\psi} - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi) = \tilde{\psi} - s\tilde{\pi(\psi)},$$

which can be written as

$$\left[\begin{array}{c} \psi_{1}\\ \psi_{2}\\ \end{array}\right] = \left[\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2}\\ \end{array}\right] - \mathcal{H}\,\partial_{1}(\psi)\left[\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2}\\ \end{array}\right] = \left[\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2}\\ \end{array}\right] - s\left[\begin{array}{c} \tilde{\pi}_{1}(\psi_{1})\\ \tilde{\pi}_{2}(\psi_{2})\\ \end{array}\right],$$

where $\mathcal{H}$ is the inverse Hessian matrix.

So the problem seems to consist of $n$ linear optimizations, not a single equation depending on a vector of parameters. My problem here is the use of the vectors $\psi$ in this case above.

So, can I optimise this without splitting the n-dimensional (2 dimensional in the example above) vector into $n$ equations which I can optimise using functions as described in the paragraph $^{*}$ above?

Another issue is that the loglikelihood function is not explicitly dependent on $\psi$, so how can i use this in the line search method?

Ok, so Kalman filters are not my forte, but non-linear optimization is, so let's give it a go.

You want to obtain $\psi$, not by solving your equation (3), but by maximizing (1), so

$\psi^+ = arg\max_\psi (\log L(Y_n | \psi))$.

Your stepsize $s$ is in theory found by

$s^+ = arg\max_s (\log L(Y_n | \psi+s\pi))$, however in practice you don't need to find the exact minimum, it simply needs to fulfil the Wolfe conditions, as described in Nocedal & Wright.

BFGS algorithms have the advantage that the step size is typically around 1, and the simple Armijo linesearch works quite well. You will probably get faster convergence with the zoom linesearch.

This article has a nice description of a Kalman filter using the LBFGS algorithm.

• Thanks very much for your reply. My problem is not with the method so much. It is more that I the function (1) is not an explicit function of the vector $\psi$ AND $\psi$ is a vector. All BFGS methods i have implemented before have accepted 'vector's but these were parameter vectors, that is a function $F(x_{1}, x_{2}) = 2 x_{1}^{2} + 5x_{2}$, here we actually have a vector in the true sense; (3) is a vector equation, of single variables... I will try and classify the problem further in an extension to the question; give me a minute. – MoonKnight Jun 24 '14 at 11:02
• That $\psi$ is a vector should not matter; as you note in your question, BFGS is a nonlinear programming algorithm that can be applied to scalar objective functions with vector-valued inputs. If you can calculate (1) at all from $\psi$ as an input, and you can calculate reasonable gradient information (needed to approximate the Hessian matrix), you should have all the information you need to use the LBFGS algorithm. – Geoff Oxberry Jun 24 '14 at 21:36
• I take the point that I am minimizing the loglikelihood function and this is the same as what I have outlined in the comment above. However, the likelihood function (1) is not explicitly dependent on $\psi$. In my case $\psi_{1}$ to $\psi_{n}$ are the diagonal components of my variance matrix $\mathbf{H}$ and $\mathbf{F}(\mathbf{H})$. So in order to satisfy the Wolfe conditions I will have to check whether the loglikelihood increases from the previous iteration of the algorithm; but what do I use to update $log(L)$? Can I construct and new $\mathbf{H}$ from $\psi_{1-n}$ and use that? – MoonKnight Jun 26 '14 at 8:27
• You would construct a new $\bf H$ from $\psi_{n-1}+s \pi$ and use that to evaluate (1). Consider if you can reformulate the problem to calculate $\log L_n(s)$ fast. For the BFGS algorithm to be truely efficient, evaluating (1) should be much faster than calculating the gradient (2), as is typically the case. – LKlevin Jun 26 '14 at 12:17