Parameter estimation for stochastic differential equation from discrete observations

Suppose we have a time-series $$x(t_i)$$ at discrete times $$t_i$$ and we want to estimate the parameters of an underlying SDE corresponding to this time-series:

$$dx_t = f(x_t,\theta)dt + \sigma(x_t,\theta)dB_t$$

1. What is the best method (according to the type of data/model) ?

2. Is there a Python or Matlab toolbox doing the job ? I have had a look at SDE Toolbox for Matlab but it is not clear to me how to use it, so it would be appreciated to provide an example.

I think it can be useful to have an overview of this topic since it may be interesting for many people. Thank you.

You question is extremely broad. You're asking about the class of all SDEs. I'm assuming that you're interested in numerically integrating some form of SDE in order to estimate parameters.

A comprehensive introduction to the subject of SDE integration is Kloeden & Platen's book Numerical Solution of Stochastic Differential Equations, which, style- and code-wise is a bit dated now, but is still good. For the Matlab user, another fine (and shorter) introduction is this paper:

Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302

The demonstration Matlab files listed in the paper can be found here now.

Best method? Like with ODEs, there is no best method. It depends on the system. Adaptive step-size Runge-Kutta methods work for a huge class of ODE problems, but SDEs, with their noise/diffusion term, are more complicated. Without knowing anything about your system (the diffusion function $\sigma(x_t,\theta,t)$ in particular) or what stochastic formulation you're using, I can't say much. To start, use the Euler-Maruyama method if you have an Itô SDE or additive noise (i.e., the diffusion does not depend on the state variable, $\sigma(x_t,\theta,t) = \sigma(\theta,t)$) and the Euler-Heun method if you have a Stratonovich-formulated SDE with non-additive noise. These are the workhorses. Higher-order schemes are trickier to implement and are usually designed for particular types of systems/noise.

Finally, if you want a Matlab Toolbox that is still under active development, is much faster, and will be easy to use if you've ever used Matlab's ODE suite (e.g., ode45), try my SDETools Matlab Toolbox for the Numerical Solution of Stochastic Differential Equations hosted on GitHub.

• Thank you, but your answer concerns the simulation of SDEs rather than the parameter estimation. Oct 14 '13 at 17:33

Part 1: Methods

There is no best method, but most methods end up in quite a similar area. Here's an old review which is still quite good, though some of the language and specific methods need to be updated. The idea is either to develop some loss function against your data or take a Bayesian route. Let's look at the two separately.

Loss functions are an optimization-based approach. You essentially choose properties about the data you want to match. For example, method of moments and generalized method of moments are simply the process:

1. Solve N times
2. Compare the difference of some average quantity of your Monte Carlo solution against the data. Make the difference your scalar loss function $$L(y)$$ where $$y$$ is your numerical solution.
3. Use whatever optimization algorithm to find the $$y$$ that minimizes the function $$L(y)$$

Then you mix and match. Use whatever SDE discretization method you like with your chosen loss function and local/global optimization algorithm, and that's a potential parameter estimation method. The simplest is Euler-Maruyama with $$L_2$$ loss on the expected value:

$$L(y) = \sum_i \Vert E[y(t_i)] - E[d_i] \Vert$$

for discrete data points $$d_i$$ at time $$t_i$$, and then throw that into an optimization package. Inside of that is a parameter $$N$$ for how many times you need to solve the SDE, and the higher $$N$$ is the more accurate your expected value is. Doing this on the means is generally referred to as the method of moments approach, where you could also do some $$E[h(y(t_i))]$$ which is then generalized method of moments.

Very much related to method of moments approaches is likelihood-based approaches. In this case, you assume that your data goes according to some distribution, solve your SDE $$N$$ times, compare the two distributions, and then use the distance of the distributions (KL-divergence, etc.) or the likelihood to calculate the fit. In this case you get another $$L(y)$$ which you put into an optimizer. To do this, you have to have some pre-specified likelihood as your input data, and usually this is done by assuming that the data has a distribution and performing a maximum-likelihood fit. For example, you can assume that your data is normally distributed at each $$t_i$$, then find the mean and variance (those are the sufficient statistics for MLE of the normal distribution), and so then you assume that at time $$t_i$$ you must have $$N(\mu_i,\sigma_i)$$, and calculate the loss against this. This approach can make it much easier to incorporate distributional assumptions about the errors, and allows correlations to be taken into account by using multivariate distributions. In some sense, it's a super-set of the previous approach. For ODEs, minimizing the $$L_2$$-loss is equivalent to maximizing the log-liklihood under the assumption of constant $$\sigma$$. While I don't know of a result for SDEs, I assume there's something similar here.

Let me make a quick remark that you don't have to sample $$N$$ solutions. Instead, you can sample from the solution's distribution directly using "exact" methods (see this link or by solving PDEs for the distribution directly, though these methods can have trouble scaling to high-dimensional problems. But yes, mix and match these techniques with choices for loss functions and optimizers and it covers a large amount of what people have done.

The other approach is Bayesian. In this case, you again need to assume some likelihood on your data, but now you assume some prior distribution on your parameters, and then use an MCMC technique where you

1. Sample parameters from the prior
2. Calculate the likelihood for those parameters (solving the equations $$N$$ times, etc)
3. Use the MCMC rules to choose whether to reject the sample
4. Choose new parameters to step to, and repeat.

Just like with ODEs, this process is much more expensive but instead of getting a point estimate out you get posterior distributions for your possible parameters.

While there are quite a few published methods, when you dig into the details it's usually some variant of how to make the choices for the steps in these overarching ideas.

Part 2: Toolboxes

I do not know of a toolbox in MATLAB or Python for this. Using what I described above you can put together some SDE solver toolbox, such as @horchler's or JiTCSDE, with an optimization routine. However, there are much more developed tools in Julia (including a more substantial selection of SDE solvers, and the built-in parameter estimation tooling from DifferentialEquations.jl currently handles the optimization-based approaches (here's an example), with the Bayesian approaches slated as coming soon.