Arbitrary-Order Finite-Time Corrections for the Kramers–Moyal Operator

With the aim of improving the reconstruction of stochastic evolution equations from empirical time-series data, we derive a full representation of the generator of the Kramers–Moyal operator via a power-series expansion of the exponential operator. This expansion is necessary for deriving the different terms in a stochastic differential equation. With the full representation of this operator, we are able to separate finite-time corrections of the power-series expansion of arbitrary order into terms with and without derivatives of the Kramers–Moyal coefficients. We arrive at a closed-form solution expressed through conditional moments, which can be extracted directly from time-series data with a finite sampling intervals. We provide all finite-time correction terms for parametric and non-parametric estimation of the Kramers–Moyal coefficients for discontinuous processes which can be easily implemented—employing Bell polynomials—in time-series analyses of stochastic processes. With exemplary cases of insufficiently sampled diffusion and jump-diffusion processes, we demonstrate the advantages of our arbitrary-order finite-time corrections and their impact in distinguishing diffusion and jump-diffusion processes strictly from time-series data.


Introduction
The reconstruction of stochastic evolution equations from time-series data in terms of the Langevin equation and the corresponding Fokker-Planck equation is often challenged by the inevitably finite temporal sampling of time-series data. Moreover, the Fokker-Planck equation is restricted to continuous stochastic processes, i.e., diffusion, and thus cannot adequately describe discontinuous transitions in time-series data. A more general description of continuous and discontinuous stochastic processes can be constructed using the Kramers-Moyal equation [1,2] given by of which the Fokker-Planck equation is a particular case (D n ≡ 0 for n > 2). The Kramers-Moyal equation serves as a stepping stone to adequately describe time-series data with both diffusive and discontinuous characteristics, but it is nevertheless challenged by finitetime sampling in real-world data. Recent applications of the Kramers-Moyal equation include brain [3,4] and heart dynamics [5], stochastic harmonic oscillators [6], renewableenergy generation [7], solar irradiance [8], turbulence [9], nano-scale friction [10], and X-ray imaging [11,12]. Previous work has demonstrated that a finite sampling interval ∆t not only influences the first-and second-order Kramers-Moyal (KM) coefficients [13] but also causes nonvanishing, higher-order (>2) coefficients [4,[14][15][16][17][18][19][20]. A more recent example is the jumpdiffusion process discussed in reference [3] dX t = a(X t , t)dt + b(X t , t)dW(t) + ξ(X t , t)dJ(t), where a(X t , t) is a drift function, b(X t , t) is the diffusion associated with an uncorrelated Brownian motion or Wiener process W(t), J(t) is a Poisson process with jump rate λ, independent of W(t), and ξ(X t , t) is Gaussian-distributed N (0, s) with zero mean and variance s. For such jump-diffusion processes, additional influences of the finite temporal sampling need to be taken into account. As shown in reference [8], jump events produce terms of order O(∆t) in the KM coefficients of even orders and the jump rate and amplitude induce terms of order O(∆t 2 ) in all coefficients. These influences are heightened for bivariate jump-diffusion processes [21], since terms of order O(∆t i ), i ≥ 3 impact higher-order (≥4) coefficients [22]. Although most of the aforementioned studies reported on finite-time corrections for KM coefficients and/or conditional moments of various orders, we still lack an explicit arbitrary-order correction or a closed-form solution in which the conditional moments are represented as functions of the KM coefficients and vice versa. In this article, we derive a full expansion of the generator of the Kramers-Moyal operator in exponential form for onedimensional Markovian processes. This is equivalent to van Kampen's system-size expansion, which is taken over a finite time interval τ [23,24]. The derivations presented henceforth are generally applicable to Markovian diffusion as well as jump-diffusion processes.
On a more general level, our solution is an explicit approximate solution of the Kramers-Moyal equation [1,2], which generalises the Fokker-Planck equation for discontinuous processes [13,23,25]. Our approximation of the Kramers-Moyal operator can be taken as an arbitrary order. In particular, we focus on the solution of this partial differential equation by representing the Kramers-Moyal operator in an exponential form and equating the conditional moments with the KM coefficients after representing the exponential operator as a power series. This representation of the exponential operator can similarly be used in other problems with an equivalent formulation [26][27][28] or similar discontinuous stochastic processes with different jump distributions, e.g., the Gamma distribution [29,30].

Mathematical Background
The Fokker-Planck(-Kolmogorov) equation (Kolmogorov forward equation or Smoluchowski equation) for the conditional probability density p(x, t + τ|x , t), that is wellknown within the fields of physics and mathematics, yields the propagation in time and space of any diffusion (thus continuous) process, is given by [31] We restrict our investigation to stationary processes, hence D n (x, t) = D n (x). Equation (2) describes the evolution of, for instance, a Brownian particle (for the case D 1 (x) = 0), which results in the known heat equation, or more complicated Markovian motions with drift.
Here, one recognises the function D 1 (x), the first KM coefficient, commonly denoted as drift, and the function D 2 (x), the second KM coefficient, commonly denoted as diffusion or volatility. The Fokker-Planck equation is, nevertheless, only valid for continuous motions and thus cannot describe jump-diffusion processes as in the case in Equation (1) or other stochastic motions with discontinuous paths.
A more general equation-the so-called Kramers-Moyal equation-takes higher-order KM coefficients D n (x), n ∈ N into account where L KM denotes the Kramers-Moyal operator defined as the power series [1,2] which we will subsequently solve for τ and an appropriate starting condition by exponentiating the Kramers-Moyal operator L KM . When examining a stochastic process in terms of time-series data, there is no direct access to the KM coefficients D n (x) but rather to the conditional moments of the data. The nth-order conditional moment M n (x, τ) is given by The KM coefficients D n (x) can be retrieved from the conditional moments M n (x, τ) via When dealing with real-world data, we do not have access to infinite temporal resolution, meaning that the above limit τ → 0 is not possible. A best-case scenario is to analyse the smallest possible temporal differences. If the data are sampled at ∆t time steps, take In order to non-parametrically retrieve the conditional moments M n (x, τ) from data, a set of histogram or Nadaraya-Watson estimators can be utilised (see Refs. [29,32] for details). Here, we will focus not on how to estimate the conditional moments but rather on how to derive a set of finite-time corrections to estimate the KM coefficients from conditional moments. These can be retrieved from data with software packages like kramersmoyal [33] or JumpDiff [34] in Python or Langevin [35] in R.

The Formal Solution of the Kramers-Moyal Equation and Its Approximations
First, we explicitly derive the corrective terms and subsequently link these to the results in reference [36], connecting them to the relation between statistical cumulants and moments [13].
Let us assume a well-defined initial state of the Kramers-Moyal equation be given by δ(x − x ). The formal solution of the time-dependent Kramers-Moyal equation (3) is given by where p(x, t + τ|x , t) is a normalisable function, such that ∞ −∞ p(x, t + τ|x , t)dx = 1, ∀(t, τ). We will now proceed to show the first-, second-, third-order, and arbitrary-order approximation to the solution of this partial differential equation with this particular initial condition.

The First-and Second-Order Approximations
The first-order approximation of the formal solution of Equation (5) is given by where the large square brackets indicate that the derivation operation is limited to the terms within the brackets. The superscript [1] indicates the order of approximation. The second-order approximation is obtained in a similar fashion, now including the quadratic term from the exponential representation Equation (5), i.e., To alleviate the notation, we refer to the KM coefficient without explicit state dependencies, i.e., D n . The second-order approximation M [2] n (x , τ) of the n-th conditional moment in Equation (4) reads The first integral is only non-vanishing if n − p − m = 0 and the second integral is only non-vanishing if n − p − s = 0, with s < n. Hence, Separating the terms between those with explicit derivatives of the KM coefficients and those without, it is immediately clear that the second-order approximation follows a structure given by the partial ordinary Bell polynomialsB n,m [37] where the summation is taken over The first-and second-order approximations can be written with the help of the partial ordinary Bell polynomials with m = 1 and m = 2, respectively, where Φ [2] n incorporates all derivatives of the KM coefficients from the 2nd-order corrections, and is given by To simplify the description, we introduce a short-hand notation and take the superscript . These results are in line with those reported for diffusion-type processes [16][17][18][19]38], where the Kramers-Moyal operator L KM = L FP reduces to the Fokker-Planck operator and we are solely left with the first two KM coefficients, as in Equation (2). In particular, applying the second-order approximation in Equation (8) to the two first KM coefficients results in and truncating the sums at second order yields the expressions in reference [16].

The Third-Order Approximation
Before we introduce the general formalism for the arbitrary-order approximation, we explicitly derive the third-order approximation which leads to M [3] n (x , τ) =M [1] n (x , τ) + M [2] n (x , τ) [1] n (x , τ) + M [2] n (x , τ) n (x , τ) Notice that the first integral is only non-vanishing for the combination q + p + m = n, which can again be expressed via the partial ordinary Bell polynomialB n,m , where m = 3, for the third-order approximation. The second expression requires q + s + k = n as well as p + r = s ∧ m − r = k. Separating these again into two expressions, one with and another without derivatives, we can express the third-order approximation as M [3] n (x , τ) = (n!) τB n,1 (D 1 , . . . , where Φ [3] n incorporates all derivatives of the KM coefficients from the third-order corrections Here, we compare our derivation to the derivation of third-order approximation in Gottschall and Peinke [16]. We note that our derivation takes the general form of the Kramers-Moyal operator, to which the Fokker-Planck operator is circumscribed. From Equation (9), we derive an identical expression for the Fokker-Planck operator reported in reference [16]. Since the Fokker-Planck operator is limited to second-order terms, i.e., D n ≡ 0 for n ≥ 3, the sum in Equation (10) can be express in full. For the first conditional moment M [3] 1 , we obtain the corrective terms Φ which is identical to Equation (A1) in the Appendix of reference [16]. Similarly, for the second conditional moment M [3] 2 , we obtain the corrective terms Φ which is in agreement with Equation (A2) in the Appendix of reference [16]. A similar derivation can be found in Appendix B of reference [8], which also yields congruent findings for the first two conditional moments of jump-diffusion processes. However, no explicit expression for all terms is given in either publication. As a simple rule of thumb, one can confer if the result is correct, as follows: the sum of the order of the KM coefficients subtracted by the derivation operation must equal n, the order of the conditional moment being calculated. In the notation used in this work, the sum of subscripts minus the sum of superscripts must equal the order n of the coefficient under investigation.

Arbitrary-Order Approximation
We now derive the arbitrary-order corrections of the Kramers-Moyal operator. This is done by induction from the previous derivations, whilst disregarding any emerging terms with derivatives of the KM coefficients with σ(k, n) a partition of a set of k ∈ N obeying Equation (7). This, in turn, is the same as a collection of partial Bell polynomials, namely where we combine terms with derivatives in Φ [k] n . If we disregard the derivative terms, the summation has an upper bound, namely m ≤ n. This is directly seen as the Bell polynomials are similarly bounded, and thus we arrive at neglecting the derivative terms Φ.
From the perspective of estimation, the aim is to determine the KM coefficients D n (x ), however what we have expressed here is the relation of the conditional moments M n (x , τ).
As we now have an explicit relation in terms of partial Bell polynomials, we will invert the relation and express the KM coefficients D n (x ) as functions of the conditional moments M 1 (x , τ), . . . , M n (x , τ).
Note that the first conditional moment M 1 (x , τ) is solely a function of the first KM coefficient D 1 (x ). The second conditional moment M 2 (x , τ) is a function of the second KM coefficient D 2 (x ), and by substitution, a function of the first conditional moment M 1 (x , τ), given by Equation (11). Subsequently M 3 (x , τ) is a function of D 3 (x ), M 2 (x , τ), and M 1 (x , τ). Thus, by recursively substituting the n − 1 KM coefficients by their expressions via the conditional moments, we obtain a relation of D n (x ) as a function of the M n (x , τ), M n−1 (x , τ), . . . , M 1 (x , τ) conditional moments.
We can then utilise the reciprocal relations of the partial exponential Bell polynomials: for a set of variables y 1 , . . . , y n , defined as functions of n other variables x 1 , . . . , x n given by the inverse relation holds With this, we can finally express any KM coefficients D n (x ) from the nth-order power series expansion, neglecting the derivative terms Φ, as We note here that these relations are equivalent to the relation between cumulants and (non-central) moments of a probability distribution [13,36]. Let M(y) be the momentgenerating function, such that with µ n the (non-central) moments and K(x) the cumulant-generating function. For n < 4, the cumulants κ n and the central moments are the same (e.g., the mean and variance). This is not the same for higher cumulants and moments. The relation between the cumulants κ n and the (non-central) moments µ n is given by the reciprocal relation of the Bell polynomials, as in Equations (12) and (13). This is in line with our exponential representation of the Kramers-Moyal operator. Here, the KM coefficients are the cumulants (with the exception of the τ term).

Exemplary Cases with Constant Diffusion and Constant Jumps
Here, we present two illustrative examples: first, a constant diffusion process, the Ornstein-Uhlenbeck process; secondly, we augment this process with jumps to obtain a jump-diffusion process. We implement the corrective terms derived thus far to show the impact of the finite-time corrections. This choice of parameters, i.e., constant diffusion and constant jumps, considerably simplifies Equation (1) to where −aX t is the state-dependent linear drift function, with a > 0, also denoted meanreverting strength, b > 0 a constant diffusion, W(t) a Brownian motion or Wiener process, ξ a state-independent and normally distributed jump amplitude with zero mean and variance s, and J(t) a Poisson process with jump rate λ. Note that the conventional Ornstein-Uhlenbeck process is recovered if we omit the jump process. We have derived an expression for the conditional moments M n (x, τ) as a function of the KM coefficients D n (x), given by Equation (11), which is valid for any Markovian diffusion or jump-diffision process. For our particular application to the Poissonian jumpdiffusion process in Equation (1) we require at least the first six KM coefficients/first six moments. These are given by We invert this expression explicitly using Equation (14) and report on the KM coefficients as functions of the conditional moments, which are given by We again note that these expressions are valid for any case of diffusion and jumpdiffusion processes. In the first case, where there are no jump terms in Equation (15), i.e., the Ornstein-Uhlenbeck process, we know that all KM coefficients D n (x) with n ≥ 3 are zero. However, this is not the case when estimating the coefficients from time-series data, i.e., from one realisation of the stochastic process sampled at finite resolution. It is common to find that these terms do not vanish due to finite-time effects. In our second case with a jump-diffusion process, the KM coefficients D n (x) with n ≥ 3 can be related directly to the jump parameters. These relations were derived in reference [3], and are given by where ξ 2n = (2n!) 2 n (n!) ξ 2 n = (2n!) 2 n (n!) s n , for Gaussian distributions with zero mean and variance s.
We will now compare the derived theoretical corrections to KM coefficients estimated from numerically generated time-series data. In Figures 1 and 2, we display the second-, fourth-, and sixth-order KM coefficients D 2 (x), D 4 (x), and D 6 (x) estimated with the firstorder, second-order, and full-order approximations given by Equation (16) (or in general Equation (14)). The full-order approximations have the same order as the KM coefficients, i.e, second-, fourth-, and sixth-order approximation for D 2 (x), D 4 (x), and D 6 (x), respectively. For the data shown in Figure 1, we use a Euler-Maruyama scheme to numerically integrate an Ornstein-Uhlenbeck process Equation (15) (without the jump terms) with parameters: drift a = 1.0 and diffusion b = 0.5 (λ = s = 0.0). We numerically integrate this process with a coarse time-step ∆t = 0.1 to deliberately emphasise the finite-time effects on the aforementioned KM coefficients. For example, the second-order KM coefficients D 2 (x) takes a quadratic form, despite the fact that the diffusion term is constant. The KM coefficients D 4 (x) and D 6 (x) are not truly zero, as would be expected for purely diffusive processes [39,40], due to the finite-time effects, but the full-order finite-time correction approximates the theoretical values with far greater detail. , and D 6 (x). There are no corrections to the drift term D 1 (x). Note that for D 2 (x) the 2nd-order and the full-order corrections are identical. The impact of solely using the second-order approximation for D 6 (x) is also evident. The coefficients D 4 (x) and D 6 (x) are theoretically zero. In all cases, the improvements to the estimation of the respective KM coefficients D n (x) are clear. The grey lines indicate the theoretical values. The numerical integration has a total time of 5 × 10 5 and a time step ∆t = 0.1 (5 × 10 6 datapoints) with a Euler-Maruyama scheme [33]. Consistent results are obtained when considering a much finer numerical time step ∆t of integration and subsequently down-sampling the data.
For the data shown in Figure 2 we follow a similar approach, now augmenting the Ornstein-Uhlenbeck process with Poissonian jumps, i.e., as given in Equation (15). The parameters are as follows: drift a = 0.5, diffusion b = 0.5, jump amplitude with a Gaussian distribution with variance s = 0.75 and zero mean, a Poissonian jump rate λ = 0.6, and a time step ∆t = 0.05. For this process, we know the higher-order KM coefficients D 4 (x) and D 6 (x) reflect the presence of discontinuous paths, which, for our particular case of the Poissonian jump-Ornstein-Uhlenbeck process, we know the explicit inversion in Equation (17) (cf. reference [3]). For the chosen coarse time step, we notice that the estimations do not correspond exactly with the theoretical values, regardless of the order of finite-time correction chosen. This can likely be traced back to the limitations of the Kramers-Moyal equation to fully capture discontinuous stochastic processes (cf. reference [41]). Nevertheless, the higher-order finite-time corrections approximate the theoretical values with greater accuracy. The abscissas are re-scaled by the standard deviation σ x of x. From left to right are shown the second-, fourth-, and sixth-order KM coefficients D 2 (x), D 4 (x), and D 6 (x). There are no corrections to the drift term D 1 (x). Note that for D 2 (x) the 2nd-order and the full-order corrections are identical. The impact of solely using the second-order approximation for D 6 (x) is also evident. The coefficients D 4 (x) and D 6 (x) are theoretically zero. In all cases, the improvements to the estimation of the respective KM coefficients D n (x) are clear. The grey lines indicate the theoretical values. The numerical integration has a total time of 5 × 10 5 and a time step ∆t = 0.1 (5 × 10 6 datapoints) with a Euler-Maruyama scheme [33].
We note here that the parameter estimation from data heavily depends on the number of data points and the sampling rate of numerically simulated or real-world time-series data. Real-world time-series data can often be sampled at higher sampling rates, but not always in such a large number of datapoints. A closer inspection of the limitations of both the sampling rate and the number of data points in parameter estimation is necessary, but falls outside the scope of this publication. Moreover, it should be emphasised that, prior to any examination of time-series data within the purview of either the Fokker-Planck or the Kramers-Moyal equation, the Markov property of the data must be account for, i.e., a vanishing memory of the increments of the data. This can be examined, for example, via the Chapman-Kolmogorov equality [13].
Summarising our findings, we conclude that our proposed arbitrary-order finite-time corrections considerably help in differentiating one-dimensional purely diffusive processes and jump-diffusion processes, as these accurately show that higher-order KM coefficients D n (x), n ≥ 2 vanish for purely diffusive processes. These arbitrary-order finite-time corrections should now also be considered for N-dimensional stochastic processes. A first examination of the second-order finite-time corrections for two-dimensional processes was recently addressed in reference [22]. Note that the one-dimensional second-order finite-time correction for these KM coefficients was recently addressed in another publication [34]. Here, it is extended to arbitrary order.

Implementation: Symbolic Calculations in Python
In this section, we implement the main results from above to compute the moments into available software packages, e.g., kramersmoyal [33] and JumpDiff [34] in Python or Langevin [35] in R, or any self-made parametric or non-parametric estimator. In order to facilitate numerical implementations of the higher-order corrections, we include a short Python script to obtain the non-derivative corrections to any desired order, with the desired truncation of the power-series expansion.

Conclusions
We have presented a set of arbitrary-order finite-time corrections to the Kramers-Moyal operator, solved by exponentiating the Kramers-Moyal operator, equivalent to van Kampen's system-size expansion. We expressed the exponential operator as a power series and worked out each element of the series, ultimately combining it in a series representation via the partial Bell polynomials. We obtained a closed form for the set of arbitrary-order finite-time corrections relating the conditional moments to the Kramers-Moyal coefficients. Moreover, by representing the arbitrary-order finite-time corrections with partial Bell polynomials, we derived a reciprocal relation for the conditional moments and the Kramers-Moyal coefficients. This provided a closed-form representation of the Kramers-Moyal coefficients via conditional moments, which is crucial for timeseries data estimation. We included two illustrative cases of poorly sampled diffusion and jump-diffusion processes with constant diffusion and constant jumps, demonstrating the suitability of our corrections for a non-parametric estimation of higher-order Kramers-Moyal coefficients. Our corrections approximated the theoretical values with a high degree of accuracy and help to distinguish processes with and without jumps. We are confident that our arbitrary-order finite-time corrections contribute to an improved reconstruction of stochastic evolution equations from empirical time-series data.