Variational inference of the drift function for stochastic differential equations driven by Lévy processes¹¹footnotemark: 1

Stochastic dynamical systems are often used to describe random phenomena in real life. Usually, stochastic differential equation[1, 2] is used to model the stochastic process, but one of the important issues is the fitting of the model and the observed data, that is, the system identification in the data. The stochastic differential equation consists of a drift function and a diffusion part. For the thermodynamic equilibrium model, the diffusion function is proportional to the identity matrix, and the drift function is the gradient of potential energy. A simple and known method for estimating the drift function from the data is obtained from[3], which takes advantage of the fact that the potential function is proportional to the logarithm of the stationary density. Therefore, we can get an explicit representation of the drift function by using the kernel density (KDE) estimator. However, this estimator is only based on the distribution of the data, completely ignoring the time sequence of the observations and the time lag between the data. Kernel density estimator is non-parameterized, it does not require the drift function satisfies the specific parameterized form. For non-equilibrium model, in general, the explicit representation between the drift function and the density is unknown. At the same time, for the higher dimensional case, the convergence of the kernel density to the real density will be slow with the increase of sample data[4].

There are many methods for estimating the parameters of the drift function[5, 6]. For these parametric and non-parametric methods, we have to deal with the problem of low data sampling rate[7]. For example, the non-parametric method of conditional expectation is used[8], which requires numerical solution of the Kolmogorov backward equation at a given time interval, and also requires a large amount of data to calculate the conditional expectation. The Bayesian estimator using the Gaussian process prior gives a good estimate of the complete path under dense observations[9]. However, generally speaking, this method requires that unobserved diffusion paths be classified as hidden random variables between adjacent observations, which may lead to time-consuming calculations or require further approximation. Papaspiliopoulos et al.[9] introduced the Monte Carlo Gibbs sampling technique, which switched back and forth between hidden path sampling and drift function sampling of the process. For the processing of hidden paths, the expectation maximization method[10] can be used to approximate the hidden process with linear stochastic differential equations. However, these methods only consider the stochastic differential equations driven by Brownian motion. In nature, there exists a class of stochastic differential equations driven by non-Gaussian Lévy motion. The Lévy motion[11, 12], like the Brownian motion, is a type of Markov process with independent and stationary increments, which has a wide applications in the fields of physics, biology, earth science, etc. Therefore, it is of great significance to study the data recognition problem of this kind of the stochastic differential equation. For the identification problem of stochastic differential equations driven by Lévy motion, Li and Duan[13] studied this equation by using parameterized identification method based on Kramers-Moyal coefficient. However, due to the need to solve the conditional expectation value, it requires a high amount of data and increases the calculation cost.

In this paper, the main purpose is to construct an effective nonparametric inference of the drift function of a class of stochastic differential equations driven by $\alpha$-stable Lévy motion, which can reduce the requirement for data sample size to a certain extent. In our method, we construct the Kullback-Leibler divergence[14] of the transition density through the diffusion process with two different drift functions, and use the variational formula based on the stationary Fokker-Planck equation[15] to derive the explicit form of the drift function. Combined with data information, we use empirical distribution to replace stationary density, and further minimize empirical functional to obtain target results. However, if the empirical functional has a penalty function term based on the kernel, then the method can be extended to nonparametric inference. In addition, the properties of $\alpha$-stable Lévy motion will result in a large variation due to the different $\alpha$ values, we also carry out numerical experiments for different parameters $\alpha$.

This paper is arranged as follows. In section 2, We derive the variational formula of the Fokker-Planck equation corresponding to the stochastic differential equation driven by an $\alpha$-stable Lévy motion, and we use empirical distribution to replace stationary density and make full use of data information to give the expression of the drift function. In section 3, We show the feasibility of this variational method through numerical application examples. We end the paper some conclusions and discussions in section 4.

In this work, we consider the stochastic differential equation for the dynamics of a $d$-dimensional diffusion processes $X_{t}\in\mathbb{R}^{d}$ given by

where $g(\cdot)\in\mathbb{R}^{d}$ is the drift function, the diffusion function $\sigma(\cdot)$ is the $d\times{k}$ dimensional matrix, $B_{t}$ is the Brownian motion in $\mathbb{R}^{k}$ and $L_{t}^{\alpha}$ is a symmetric $\alpha$-stable L$\mathrm{\acute{e}}$vy motion in $\mathbb{R}^{d}$ with the generating triplet $(0,0,\nu_{\alpha})$. The jump measure[16]

with $c(n,\alpha)=\frac{\alpha\Gamma((n+\alpha)/2)}{2^{1-\alpha}\pi^{n/2}\Gamma(1-\alpha/2)}.$ The processes $B_{t}$ and $L_{t}^{\alpha}$ are taken to be independent.

We will give a simple numerical example in this section to illustrate the results through the use of different amounts of data and the comparison of different $\alpha$ values.

From Figure (1) and (2), we can see that the size of the data and the different $\alpha$ values ??have an impact on the inference of the drift function. In this chapter, our method only considers the data obtained from one observation of the process. When the value of $\alpha$ is the same, the estimator of the drift function becomes better as the amount of data increases. And for the same amount of data, we’re going to consider the selection of $\alpha$ values of 1, 1.25, 1.5 and 1.75, respectively. As the value of $\alpha$ increases, the sample trajectory of the process is smoother, making the inference result of the drift function better. Of course, we also know from the numerical results that the estimation of the drift function needs to be improved, such as when $\alpha$ is close to 1, the estimation effect of the drift function is poor, so we can improve the numerical results by increasing the amount of sample data. In addition, the method to infer the drift function only. For diffusion function, we are still in the further study. On the one hand, we hope to obtain the diffusion function inference in the case of additive and multiplicative noise when the drift function is known. On the other hand, we also hope to obtain an optimization method to infer the drift function and the diffusion function at the same time.

In this paper, we have proposed a way to obtain the drift function estimator for a class of stochastic differential equations driven by $\alpha$-stable Lévy motion from data. We construct a variational formula based on the stationary Fokker-Planck equation, replace the stationary density with an empirical distribution, and derive the drift function expression. Another contribution is that our method only needs to observe the dense data of the process once, which reduces the cost of data processing, and is more in line with the actual observation situation, so it is more conducive to processing the actual data.

In addition, there are some interesting extensions for this study. We only studied the nonparametric estimation of the drift function of the stochastic dynamical system driven by Lévy motion. It is worth thinking about the inference of the diffusion function for additive and multiplicative cases. On the other hand, we can also explore the method of extracting the dynamics of this type of the stochastic system from the data, so as to gain insight into the stochastic phenomenon in related fields.

We would like to thank Dr. Xiaoli Chen, Dr. Wei Wei, Dr. Cheng Fang for helpful discussions. This work was supported by the National Natural Science Foundation of China (NSFC) grants 12001213, 11801192, 11771449 and National Science Foundation (NSF) grant 1620449.

For every $t\in[0,T]$, assuming that the dense path observation $X_{0:T}$ of equation (1) is obtained, we can derive the likelihood function. We discretize the time interval $[0,T]$, take the time step as $\Delta{t}$ and use the fact that $\Delta{t}\rightarrow{0}$, then the transition density of the equation (1) satisfies the Gaussian distribution. We find that the negative log-likelihood part that depends on the drift function $g$ can be approximated by the following equation

Where the inner product $\langle{u,v}\rangle\doteq{u\cdot{D^{-1}}v}$, and the corresponding square norm is $\|u\|^{2}\doteq{u\cdot{D^{-1}}u}$.

For a sufficiently small $\Delta{t}$, using the Gaussian form of the transition density, we can give the Kullback-Leibler divergence between the path probabilities of the two diffusion processes with drift functions of $g(x)$ and $r(x)$ respectively, where $g(x)=f(x)+r(x)$, which can be expressed as

Assume that the two diffusion processes have the same diffusion function $D(x)$ and non-stochastic initial state [19]. The density $p_{t}(x)$ is the probability density of the diffusion process with drift function $g$ at time $t$. Furthermore, we assume that the process reaches a stationary state and has a stationary density $p(x)$, thus the relative entropy rate [20] can be obtained

The comparison with equation (3) shows that when $A=D$, the minimization of equation (5) will result in that the stationary density process with drift function $g$ is close enough to the process with drift function $r$ in the case of relative entropy. Therefore, we can understand this as the generalized minimum relative entropy (Kullback Leibler divergence) under the constraints given by the stationary density.

We will show how to obtain an explicit representation of the estimator from equations (13) and (15). In order to simplify the notation, we make some notation simplification for the expressions covered in the equation. Let the vectors of $\nabla\psi(x_{i})$ and drift function $r(x_{i})$ at all data points be

The kernel function $B$ and the block-diagonal weight matrix $A$ are defined as

and define

then we have

We can rewrite (15) as

Where $I$ represents the identity matrix. Equation (13) can be similarly written as

and $k(x)=(\nabla_{x^{\prime}}K(x,x^{\prime})|_{x^{\prime}=x_{1}},\nabla_{x^{\prime}}K(x,x^{\prime})|_{x^{\prime}=x_{2}},\cdots,\nabla_{x^{\prime}}K(x,x^{\prime})|_{x^{\prime}=x_{n}}).$ Further, we can obtain the following result via solving the function $b$ of equation (23),

Combined with equation (24), the explicit form of the estimator $\psi$ can be solved