Semiparametric maximum likelihood reconstruction of stochastic differential equations driven by white and correlated noise Open Access

The oldest and most established scientific modeling strategy is the physics-based or forward or direct modeling approach. Dynamical models are built from the known governing equations of motion of the system derived from first principles. A second and complementary approach is offered by data-driven or inverse modeling; here, dynamical models are extracted from datasets of the underlying system using statistical techniques. Data-driven modeling is useful if the equations of motion of the system are not readily known and, in particular, if some sort of coarse-graining is involved, that is, if the goal is to derive reduced or effective equations modeling particular macroscopic variables or modes of the system. Such equations will then naturally have a stochastic nature and incorporate memory effects.

The reconstruction of the drift and diffusion functions in Eq. (1) from data is a non-trivial estimation problem. Possible approaches include the direct estimation method,^4,5 maximum likelihood estimation,^6,7 Kalman filtering^8,9 and Markov chain Monte Carlo methods¹⁰ in a Bayesian framework. All these techniques have advantages and disadvantages regarding issues such as finite sampling interval, partial observation of the system, observational noise, and computation time.

Estimating the drift and diffusion functions in a non-Markovian Langevin equation from time series data of only $x$ is an even more challenging task than reconstructing a Markovian Langevin equation. An extension of the direct estimation method was proposed;¹³ it is based on a perturbation approach, admits a non-zero correlation time of the driving noise, but still needs it to be much smaller than the correlation time of $x$⁠. A direct estimation method without this restriction was introduced recently.¹⁴ However, the proposed optimization procedure does not appear to be straightforward; a complicated regularization guided by trial and error for each example system is necessary. Moreover, very large datasets are required for the method to be applicable.

This study introduces and explores a methodology for reconstructing drift and diffusion functions in Markovian and non-Markovian SDEs based on the unscented Kalman filter. The estimation is semiparametric as the drift and diffusion functions are expanded into nested sets of basis functions of increasing complexity. Unlike previous work using the state augmentation technique for parameter estimation,^8,9,15,16 the Kalman filter is here only used for state estimation and propagation of the probability density of the state. The model parameters are determined by maximizing the likelihood of the observations based on the predictive probability density generated by the Kalman filter.

The remainder of the paper is organized as follows: The method is described in detail in Sec. II. In Sec. III, the example systems used to validate the method and their properties are illustrated. Then, the results are given and discussed in Sec. IV. Section V concludes.

Depending on the domain of the system, the method could be extended by using also other types of basis functions such as spline, trigonometric, logarithmic, or rational functions which would then also be included in the orthonormalization, but this will not be discussed here. A similar strategy of combining different types of basis functions was recently applied in the context of maximum likelihood probability density estimation.¹⁷ 

In the non-Markovian case, the correlation time parameter $γ$ is also to be estimated. It is represented as $γ= e η$ to guarantee positivity. The total parameter vector of the Langevin model is $ω= ( A 1 , … , A I , B 1 , … , B J ) T$ in the white noise case and $ω= ( A 1 , … , A I , B 1 , … , B J , η ) T$ in the correlated noise case. There are no constraints on any of the parameters. When quoting values of expansion coefficients, we always refer to the parameters ${ α i } i ∈ I$ and ${ β j } j ∈ J$ of the monomial expansions in Eqs. (11) and (12), which are linked by a linear transformation to the parameters ${ A i } i = 1 I$ and ${ B j } j = 1 J$⁠.

The unscented Kalman filter^18,19 (UKF) is used for state estimation and propagation of state probability densities in the Langevin models. The UKF allows for recursive estimation of unobserved states in stochastic nonlinear systems from incomplete observations. It keeps the full nonlinear system dynamics rather than linearizing it but truncates the filter probability density to a Gaussian in each step by only propagating means and covariances.

An ensemble of size $M$ of parameter vectors consistent with the parameter uncertainty is generated as ${ ω ^ + L v i } i = 1 M$⁠, where $v i$ are column vectors of length $p$ of independent standard normal random variables and $L$ is obtained from the Cholesky decomposition $C= L L T$⁠. This ensemble converts to an ensemble of drift and diffusion functions from which confidence intervals can be constructed. In this study, ensembles of size $M=1000$ are used to calculate 90% confidence intervals, delimited by the 5% and 95% quantiles.

A possible ill-conditioning or ill-posedness of the parameter estimation problem can be monitored with the (spectral) condition number of $J$ given as $κ= λ max/ λ min$⁠, where $λ max$ and $λ min$ are the largest and smallest eigenvalues of the symmetric and positive definite matrix $J$⁠.

The estimation methodology is explored and validated on a couple of simulated datasets from prototype example systems with known dynamics of increasing complexity. Similar systems were discussed in Ref. 14.

For each system, first the model selection is performed by increasing the number of drift and diffusion basis functions in steps of one and finding the model with the smallest BIC. Then, the reconstructed drift and diffusion functions for the chosen model are examined. Additionally, we look at some important summary statistics in long-term integrations of the true and the reconstructed model: the (stationary) probability density function of the process $x(t)$⁠, the probability density function of the increment process $x(t+δt)−x(t)$⁠, and the autocorrelation function of $x(t)$⁠. These summary statistics are calculated from very long time series of $4× 10 6$ data points at a sampling interval of 0.025. The probability densities are estimated using a recently proposed semiparametric maximum likelihood estimator¹⁷ unless they are available analytically. The correlation times of the systems are estimated by evaluating the integral of Eq. (4) from the sample autocorrelation function using the (composite) trapezoidal quadrature rule. In practice, a finite upper limit of the integral has to be adopted. As the sample autocorrelation function is rather noisy at large time lags, even if calculated from a very long time series, we monitor the correlation time as a function of the upper limit and look for a plateau.

We first consider a linear model with strongly multiplicative noise.

The parameter values $a=1$⁠, $x 0=1$⁠, $b=1$⁠, $c=1$⁠, and $x ∗=2$ are chosen here. The system exhibits a markedly positively skewed probability density; the mode of the distribution is below the mean value, and there is a long right tail in the region of high noise level. The correlation time of the system is $T x=1$⁠. The reconstruction is based on a time series of length $N=5000$ with sampling interval $δt=0.1$⁠.

Next, we examine a nonlinear bistable system with parabolic diffusion function. The drift function has odd and the diffusion function has even symmetry.

The parameter setting here is $a=1$⁠, $b=−1$⁠, $c=0.8$⁠, and $d=0.25$⁠. The dynamics are bistable with stable fixed points at $x 0 ( 1 , 2 )=±1$⁠; the probability density is bimodal with the modes at $x m ( 1 , 2 )≈±0.87$⁠. The correlation time is $T x=2.3$⁠. A time series of length $N=10000$ at a sampling interval of $δt=0.2$ is used for the reconstruction.

Finally, we consider the Rayleigh process, here extended to a setting with multiplicative and correlated noise. This system has a non-polynomial drift function.

The parameter setting $a=1$⁠, $b=1$⁠, and $c=1$ is adopted here. The system has a stable fixed point at $x 0=1$⁠; the mode of the probability density is at $x m≈0.78$⁠. The correlation time of the system is $T x=0.6$⁠. A dataset of length $N=10000$ with sampling interval $δt=0.05$ is used for the estimation.

The model selection is illustrated in Tables I–VIII. The log-likelihood function at the maximum likelihood solution $ω ^$ and the BIC are listed when successively increasing the number of drift and diffusion basis functions. The values are rounded to integers as more precision is not needed here. For all example systems, the sampling interval $δt$ is divided into $L=10$ subintervals for propagating the probability densities of the states with the UKF (see Sec. II B).

For all example systems with a polynomial drift function (A.1, A.2, B.1, and B.2) the correct drift model is clearly chosen. The rational drift function of the Rayleigh process is approximated with a third-order polynomial in the white noise case and a fourth-order polynomial in the correlated noise case.

For all example systems, the true diffusion function lies outside the model classes covered by semiparametric modeling. For systems A.1, A.2, B.1, and B.2, the logarithm of the diffusion function can be well approximated on the data range by a quadratic polynomial; for system C.1 a fourth-order, and for system C.2, a fifth-order polynomial is chosen.

In the examples B.1 and B.2, the pattern of symmetry in both drift and diffusion is clearly visible (Tables III and V). A significant model improvement is only achieved when adding a basis function with the correct symmetry. This suggests rerunning the model selection with only odd drift and even diffusion basis functions, which leads to more parsimonious models (Tables IV and VI).

A comprehensive investigation of the model selection would involve running ensembles of realizations of the example systems and monitoring how often which model is chosen in the reconstruction as a function of the length of the data set $N$ and the sampling interval $δt$⁠. We leave this to future work.

The reconstructed drift and diffusion functions of the various example systems are displayed in Figs. 1 and 2 together with pointwise 90% confidence intervals. A canonical choice for the block length in the uncertainty estimates (see Sec. II D) is $l b=[ T x/δt]$⁠, the largest integer smaller or equal to $T x/δt$⁠, that is, the number of data points per correlation time. Outside the range of the learning dataset, a linear extrapolation is applied to the drift functions and a constant extrapolation to the diffusion functions. For the systems driven by correlated noise (A.2, B.2 and C.2), the estimates of the correlation time parameter $γ$ with their uncertainties are listed in Table IX.

Within the range of the learning data sets both the drift and diffusion functions are very accurately identified. In the long tails of the linear models and the Rayleigh processes the errors and the uncertainties are somewhat larger. The correlation time parameter $γ$ is also reasonably well estimated in all cases. The uncertainties of the diffusion function estimates are clearly smaller for the white noise models than for the correlated noise models, whereas the drift functions are slightly more accurately estimated in the correlated noise models. The reason for this is that in the white noise case process increments are in lowest order determined by the noise contribution and the drift enters only in the next order, whereas in the correlated noise case, drift and diffusion contribute in the same order.

Due to the symmetries of the drift and diffusion functions, the probability density in the systems B.1 and B.2 is symmetric about zero. The same holds for the reconstructed system preserving the symmetries of the drift and diffusion functions. The symmetry of the probability density is prescribed in the semiparametric estimator¹⁷ as there may be an imbalance between the two states due to their metastability even in very long time series.

For example systems C.1 and C.2, we apply a reflecting boundary condition at zero when integrating the reconstructed models. Otherwise, there would be a finite probability of reaching negative values of $x$ as the singularity of the drift function at zero cannot be captured in the reconstruction. This probability is very small, though, and the boundary at zero is actually never hit in our integrations. We here also apply a reflecting boundary condition at the largest data point in the learning dataset. The estimate of the drift function in the right tail of the data set is very uncertain due to the high noise level there and the extrapolation may lead to a positive drift for very large $x$ and, thus, the escape of trajectories to infinity.

The summary statistics as estimated from a long-term integration of the true and reconstructed models are displayed in Figs. 3 and 4. We observe excellent agreement in all these statistics between the true and the reconstructed systems with only very small errors. Even details of heavy-tailed increment probability densities are very well captured.

The time series here extend over between 185 and 870 times the correlation time of $x$ and between 10 and 27 samples per correlation time are used. The method is fairly robust with respect to these choices. It still provides virtually unbiased estimates of the drift and diffusion functions from shorter time series; obviously, the uncertainty increases with decreasing time series length $N$⁠. Practically useful estimates can be obtained down to a length of about 50 times the correlation time. If the sampling interval $δt$ becomes too large, the model parameter estimates are biased as the Gaussian approximation of the UKF is no longer sufficiently valid. When and how strongly this effect shows up depends on the system under consideration, in particular, on the strength of the nonlinearity and the state-dependence of the driving noise.

A methodology for reconstructing Markovian and non-Markovian SDEs from time series data is introduced and discussed. It makes use of the unscented Kalman filter to estimate states and propagate their probability densities and then employs the filter predictions to estimate model parameters, leading to a principled maximum likelihood technique that naturally comes with uncertainty estimates. The drift and diffusion functions are represented in a semiparametric manner.

The technique is validated on a couple of synthetic time series datasets with known underlying dynamics ranging from linear systems, to nonlinear bistable systems, to intermittent systems. The maximization of the likelihood function is robust and fast. Drift and diffusion functions are accurately reconstructed, even if they lie outside of the prescribed polynomial model class.

The computational effort of the present method is considerably higher than that of direct estimation methods,^13,14 but it is still not too onerous and much lower than the computational cost of Markov chain Monte Carlo techniques.¹⁰ The direct estimation method of Ref. 13 has the restriction that the correlation time of the driving noise is required to be much shorter than that of the variable $x$⁠. The system reconstruction with the present method is very accurate for datasets of moderate length, whereas direct estimation methods^13,14 require extremely long datasets to produce any reasonable results. Furthermore, the present method provides a straightforward optimization procedure without the need to tune hyperparameters by hand for each example system and it comes naturally with uncertainty estimates for the drift and diffusion functions, both of which is not the case with direct estimation methods.^13,14 In summary, the Kalman filter maximum likelihood approach to identification of nonlinear stochastic systems appears to offer a very good compromise in terms of computation time, accuracy, and data requirements.

A detailed study of the model selection and reconstruction accuracy depending on the length of the time series and the sampling interval using multiple realizations of example systems is left for future work.

The authors have no conflicts to disclose.

F. Kwasniok: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review and editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.