# Kalman filter, Multivariate gaussian log-likelihood, Newton finite difference

For parameter estimation using Kalman filter technique I have obtained the negative Log-likelihood of mutivariate gaussian. My problem is how to obtain the gradient and hessian from this log-likelihood by newton finite difference method for updating the parameters e.g H*deltatheta = G, where H is hessian and G is gradient of log-likelihood corresponding to the parameters.I am using Matlab. If possible please give some example in matlab code.

## 2 Answers

In addition to the EM algorithm already mentioned in other answer, a quasi-Newton algorithm can be used to numerically optimize your likelihood function. The BFGS algorithm for unconstrained optimization and the L-BFGS-B algorithm for bound constrained optimization are widely used to this end.

The functions fminunc and fmincon do this job in Matlab. Searching for the keywords Matlab and BFGS or L-BFGS-B you may probably find other implementations or interfaces and examples.

As proposed by @Arthur B., the EM algorithm is another alternative that you could use. It has some theoretical advantages: e.g., the likelihood always increase at each step of the algorithm and the derivation of the algorithm leads to closed-form expressions for the updating equations. It is also relatively robust to poorly chosen starting values. Some disadvantages are observed in practice, nonetheless. In the context of structural time series models, it converges slowly, especially as the algorithm approaches to the local optimum. The same may happen in other contexts but I can tell it by experience only in the context of structural time series models. This is already noticed in this paper, which is the basis of the reference given by Arthur:

Shumway, R.H. and Stoffer, D.S. (1982). An approach to time series smoothing and forecasting using the EM algorithm. Journal of Time Series Analysis, 3, 253-264.

In fact, the authors of this paper suggest using the EM algorithm for the first iterations of procedure and switch to another optimization algorithm after some iterations. Good candidates for these other algorithms are the BFGS or the L-BFGS-B algorithms mentioned above.

The paper mentioned by Arthur is a recent and interesting reference but I think that the issue of convergence (in terms of number of iterations required until convergence rather than computational time) remains still as a practical drawback. But I didn't read the paper thoroughly and didn't try the proposed implementation. It is worth a try to get a hands-on feel for the performance.

Edit

You may be interested in this paper:

Koopman, S. J. and Shephard, N. (1992). "Exact score for time series models in state space form". Biometrika (1992) 79 (4) pp.823-826. http://doi.org/10.1093/biomet/79.4.823.
Link to paper from author's website.

In this work, the authors describe a method to compute the gradient of the likelihood function for Gaussian state space models.

• (+1) one thing though, EM is not (to my mind) strictly an "alternative" to [L-]BFGS in this context [State Space Methods]. The implementation of the EM algorithm to start convergence before switching to BFGS is not only an option, it is required for arbitrary initial conditions (in my experience of implementing both optimizations for n-dimensional State Space models). If there is another way of consistently initializing the BFGS algo. using derived initial conditions [not using the EM for N iterations] I would like to hear it! :] – MoonKnight Jan 14 '15 at 10:55
• I wouldn't say it is required but I agree that using the EM algorithm can be a helpful way to get starting values to be passed to the BFGS optimization algorithm. In the context of the basic structural model (level, plus trend plus seasonal component) the EM algorithm is rarely employed. However, sometimes it is hard to achieve convergence in the optimization algorithm, in these cases, having an alternative technique such as the EM algorithm can be of great help. – javlacalle Jan 14 '15 at 14:26

You should use the EM algorithm to fit the Kalman filter, rather than gradient descent. Check this paper for ideas on implementation details.

• 1) The updating equation sketched by the OP includes information from the gradient and also from the Hessian as a way to obtain the direction vector. Thus, it seems that the OP has in mind some kind of quasi-Newton algorithm rather than the gradient descent algorithm. 2) The OP seems more interested in a ready-to-use package or interface than on implementation details. Is there an implementation of this paper available in Matlab? – javlacalle Nov 6 '14 at 20:11
• I am trying to implement this paper ieeexplore.ieee.org/xpl/… . I am struck in finding the hessian and gradient of log-likelihood function by finite difference. – Mithun Mondal Nov 9 '14 at 14:06
• If you're using finite differences, you're doing it horribly wrong. This is all analytical. – Arthur B. Nov 9 '14 at 16:06
• Could you please give some idea/paper for finding the gradient and hessian of the log-likelihood function of multivariate gaussian. I am totally confused – Mithun Mondal Nov 10 '14 at 13:58
• Why do you need the gradient? Just do EM. – Arthur B. Nov 10 '14 at 14:00