Estimation of drift and diffusion functions from time series data: A maximum likelihood framework

Complex systems are characterized by a huge number of degrees of freedom often interacting in a non-linear manner. In many cases macroscopic states, however, can be characterized by a small number of order parameters that obey stochastic dynamics in time. Recently techniques for the estimation of the corresponding stochastic differential equations from measured data have been introduced. This contribution develops a framework for the estimation of the functions and their respective (Bayesian posterior) confidence regions based on likelihood estimators. In succession approximations are introduced that significantly improve the efficiency of the estimation procedure. While being consistent with standard approaches to the problem this contribution solves important problems concerning the applicability and the accuracy of estimated parameters.

💡 Research Summary

The paper addresses the problem of inferring the drift D^{(1)}(x,t) and diffusion D^{(2)}(x,t) functions that govern the stochastic dynamics of macroscopic order parameters in complex systems. While the traditional “direct estimation” approach evaluates the Kramers‑Moyal coefficients by taking the limit τ→0 of conditional moments, it suffers from sensitivity to sampling frequency, kernel choice, finite‑sample effects, and measurement noise. The authors propose a maximum‑likelihood (ML) framework that works directly with the conditional transition probability p(x_i | x_{i‑1}; Ω, τ) of a Markovian process, where Ω denotes a set of parameters that define the drift and diffusion functions.

Using Bayes’ theorem with uniform priors, the posterior probability of Ω given the observed trajectory {x_i} is proportional to the joint likelihood, i.e. the product of the transition densities. The log‑likelihood L(Ω)=∑{i=1}^N log p(x_i | x{i‑1}; Ω, τ) is maximized to obtain the ML estimate Ω̂. The authors emphasize that, when the exact finite‑time propagator is known analytically (as for the Ornstein‑Uhlenbeck (OU) process), the maximization is trivial and yields closed‑form expressions for the optimal parameters (γ̂, Q̂). For more general drift and diffusion functions, the propagator must be approximated.

Two approximations are introduced to make the method computationally feasible for large data sets:

1. Small‑τ Gaussian Approximation – Assuming τ is sufficiently small, the transition density is approximated by a Gaussian whose mean and variance are given by the first two Kramers‑Moyal coefficients multiplied by τ. This yields an explicit form for p(x_i | x_{i‑1}) that is accurate for high‑frequency data and dramatically reduces computational cost compared with solving the full Fokker‑Planck equation.

2. Non‑Parametric Local Likelihood – The state space is discretized (e.g., via histograms or kernel bins). Within each bin the drift and diffusion are treated as locally constant, and the local log‑likelihood is maximized independently. This reproduces the results of direct estimation while providing a natural Bayesian framework for confidence assessment.

Confidence intervals are derived in a Bayesian sense (Bayesian posterior confidence intervals, BpCIs). The authors define a confidence region C as the set of parameters for which the log‑likelihood does not fall below L(Ω̂) − R, where R is a threshold. Using Wilks’ theorem, R can be linked to a desired confidence level ν via the χ² distribution (R_W(ν)=−½ F^{-1}_{χ²}(ν, 1)). For the OU example, BpCIs computed with R_W(0.95) contain the true parameters for all sample sizes examined, and the empirical coverage matches the nominal 95 % level for N≥100. For very small samples the authors recommend a direct numerical integration of the posterior to verify coverage.

The paper presents extensive numerical experiments. Table I shows that as the number of observations N grows from 10² to 10⁶, the ML estimates of γ and Q converge to their true values (γ=1, Q=1) and the BpCIs shrink accordingly. The OU case also admits analytical ML formulas (γ̂=τ⁻¹ log