Learning Stochastic Dynamics from Data

Stochastic Differential Equation (SDE) is a fundamental modeling tool in various science and engineering fields Evans (2013); Särkkä” & Solin (2019). Compared with traditional deterministic models which often fall short in capturing the stochasticity in the nature, the significance of SDE lies in the ability to model complex systems influenced by random perturbations. Hence, SDE provides insights into the behavior of such systems under uncertainty. By incorporating a random component, typically through a Brownian motion, SDE provides a more realistic and flexible framework for simulating and predicting the behavior of these complex and dynamic systems.

We consider the models following the SDE of the following form, $d{\mathbf{x}}_{t}={\bm{f}}({\mathbf{x}}_{t})dt+d{\mathbf{w}}_{t}$, where the state vector ${\mathbf{x}}_{t}\in\mathbb{R}^{d}$, the drift term ${\bm{f}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, and the stochastic noise ${\mathbf{w}}_{t}$ is a Brownian motion with covariance matrix ${\bm{D}}({\bm{x}})$ which also depends on the state. In the field of mathematical finance, there are several SDE models that are widely used, for example, Black-Scholes Model Black & Scholes (1973) in option pricing, Vasicek Model Vasicek (1977) for analyzing interest rate dynamic and Heston Model Heston (2015) for volatility studies. Apart from the field of finance, SDE has its application in physics Coffey et al. (2004); Ebeling et al. (2008) and biology Székely & Burrage (2014); Dingli & Pacheco (2011) where SDE is used to study evolution of particles subjected to random forces and modeling physical and biological systems.

The calibration of these models is crucial for their effective application in areas like derivative pricing and asset allocation in finance, as well as in the analysis of particle systems in quantum mechanics and the study of biological system behaviors in biology. This requires the use of diverse statistical and mathematical techniques to ensure the models’ outputs align with empirical data, thereby enhancing their predictive and explanatory power. Since these SDE models have explicit function form of both drift and diffusion terms, the calibration or estimation of parameters is mostly done by minimizing the least square error between observation and model prediction Mrázek & Pospíšil (2017); Abu-Mostafa (2001). SDE has also been studied for more general cases where the drift term does not have an explicit form where it satisfies ${\bm{f}}={\bm{f}}({\mathbf{x}}_{t},\theta)$ with $\theta$ being a vector a unknown parameters. These models can be estimated with the maximum likelihood estimator approach presented in Phillips (1972). The other approach is to obtain the likelihood function by Radon Nikodym derivative over the whole observation $\{{\mathbf{x}}_{t}\}_{t=0}^{T}$. Inspired by theorem $7.4$ in Särkkä” & Solin (2019) and the Girsanov theorem, we derive a likelihood function that incorporates state-dependent correlated noise. By utilizing such likelihood function, we are able to capture the essential structure of ${\bm{f}}$ from data with complex noise structure.

The remainder of the paper is structured as follows, section 2 outlines the framework we use to learn the drift term and the noise level (if not state dependent), we demonstrates the effectiveness of our learning by testing it on various cases summarized in section 4 and additional examples presented in Appendix B, we conclude our paper in section 5 with a few pointers for ongoing and future developments.

We consider the following SDE

where ${\bm{f}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a drift term, and ${\mathbf{w}}$ represents the Brownian noise with a symmetric positive definite covariance matrix ${\bm{D}}={\bm{D}}({\bm{x}}):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$.

We consider the scenario when we are given continuous observation data in the form of $\{{\mathbf{x}}_{t},d{\mathbf{x}}_{t}\}_{t\in[0,T]}$ for ${\mathbf{x}}_{0}\sim\mu_{0}$. We will find the minimizer to the following loss function

for $\tilde{\bm{f}}\in{\mathcal{H}}$; the function space ${\mathcal{H}}$ is designed to be convex and compact and it is also data-driven. This loss functional is derived from the Girsanov theorem as well as inspiration from Theorem $7.4$ from Särkkä” & Solin (2019).

In the case of uncorrelated noise, i.e. ${\bm{D}}({\bm{x}})=\sigma^{2}({\bm{x}}){\mathbf{I}}$, where ${\mathbf{I}}$ is the $d\times d$ identity matrix and $\sigma:\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$ is a scalar function depends on the state, representing the noise level, 2 can be simplified to

Moreover, we provide three different performance measures of our estimated drift. First, if we have access to original drift function ${\bm{f}}$, then we will use the following error to compute the difference between $\hat{\bm{f}}$ (our estimator) to ${\bm{f}}$ with the following norm

Here the weighted measure $\rho$ is defined on $\Omega$, where it defines the region of ${\mathbf{x}}$ explored due to the dynamics defined by equation 1; therefore $\rho$ is given as follows

The norm given by equation 4 is useful only from the theoretical perspective, i.e. showing convergence. In real life situation, ${\bm{f}}$ is most likely non-accessible. Hence we will look at a performance measure that compare the difference between ${\mathbf{X}}({\bm{f}},{\mathbf{x}}_{0},T)=\{{\mathbf{x}}_{t}\}_{t\in[0,T]}$ (the observed trajectory evolves from ${\mathbf{x}}_{0}\sim\mu_{0}$ with the unknown ${\bm{f}}$) and $\hat{\mathbf{X}}(\hat{\bm{f}},{\mathbf{x}}_{0},T)=\{\hat{\mathbf{x}}_{t}\}_{t\in[0,T]}$ (the estimated trajectory evolves from the same ${\mathbf{x}}_{0}$ with the learned $\hat{\bm{f}}$ and the same random noise as used by the original dynamics). Then, the difference between the two trajectories is measured as follows

However, comparing two sets of trajectories (even with the same initial condition) on the same random noise is not realistic. We compare the distribution of the trajectories over different initial conditions and all possible noise at some chosen time snapshots using the Wasserstein distance.

When the noise level becomes a constant, i.e. $\sigma({\mathbf{x}})=\sigma>0$, we end up a much simpler loss

which has been investigated in Lu et al. (2022). System identification of the drift term has been studied in many different scenarios, e.g. identification by enforcing sparsity such as SINDy Brunton et al. (2016), neural network based methods such as NeuralODE Chen et al. (2018) and PINN Raissi et al. (2019), regression based Cucker & Smale (2002), and high-dimensional reduction framework Lu et al. (2019). The uniqueness of our method is that we incorporate the covariance matrix into the learning and hence improving the estimation especially when the noise is correlated.

In this section, we demonstrate the application of our trajectory-based method for estimating drift functions, showcasing a variety of examples. We explore drift functions ranging from polynomials in one and two dimensions to combinations of polynomials and trigonometric functions, as well as a deep learning approach in one dimension. The observations, serving as the input dataset for testing our method, are generated by the Euler-Maruyama scheme, utilizing the drift functions as we just mentioned. The basis space ${\mathcal{H}}$ is constructed employing either B-spline or piecewise polynomial methods for maximum degree p-max equals $2$. For higher order dimensions where $d\geq 2$, each basis function is derived through a tensor grid product, utilizing one-dimensional basis defined by knots that segment the domain in each dimension. The parameters for the following examples are listed in table 1. Due to space constraints, additional examples are provided in appendix B.

The estimation results are evaluated using several different metrics. We record the noise terms, $d{\mathbf{w}}_{t}$, from the trajectory generation process and compare the trajectories produced by the estimated drift functions, $\hat{{\bm{f}}}$, under identical noise conditions. We examine trajectory-wise errors using equation 6 with relative trajectory error and plot both ${\bm{f}}$ and $\hat{{\bm{f}}}$ to calculate the relative $L^{2}$ error using 4, where $\rho$ is derived by 5. Additionally, we assess the distribution-wise discrepancies between observed and estimated results, computing the Wasserstein distance at various time steps.

We have presented a novel approach of learning the drift and diffusion from a noise-guide likelihood function. Such learning can handle general noise structure as long as the covariance information of the noise is known. Unknown noise structure when the covariance is a simple scalar is still learnable through quadratic variation. However learning state-dependent covariance noise is ongoing. Testing on high-dimensional drift term is ongoing. Our loss only incorporates likelihood on the observation data. It is possible to add various penalty terms to enhance the model performance, such as the sparsity assumption enforced in SINDy and other related deep learning methods such as NeuralODE and PINN.