jumpdiff : A Python Library for Statistical Inference of Jump-Diffusion Processes in Observational or Experimental Data Sets

We introduce a Python library, called jumpdiff , which includes all necessary functions to assess jump-diffusion processes. This library includes functions which compute a set of non-parametric estimators of all contributions composing a jump-diffusion process, namely the drift, the diffusion, and the stochastic jump strengths. Having a set of measurements from a jump-diffusion process, jumpdiff is able to retrieve the evolution equation producing data series statistically equivalent to the series of measurements. The back-end calculations are based on second-order corrections of the conditional moments expressed from the series of Kramers-Moyal coefficients. Additionally, the library is also able to test if stochastic jump contributions are present in the dynamics underlying a set of measurements. Finally, we introduce a simple iterative method for deriving second-order corrections of any Kramers-Moyal coefficient.


Introduction
Models of complex systems based on stochastic processes are ubiquitous across several research fields. However, the complexity of the stochastic contributions lays often beyond pure diffusive processes, such as Brownian motion and random walks, and might include contributions from discontinuous jumps (Pascucci 2011). Jump-diffusion processes incorporate both diffusive contributions as well as stochastic jumps, making them a natural extension of pure diffusive processes (Aït-Sahalia 2004). These processes are well suited to describe data which exhibits contributions from both continuous and discontinuous stochastic noise. However, there is a lack of reliable computational libraries oriented at estimating all necessary parameters that quantitatively describe jump-diffusion processes. The problem of deriving jump-diffusion equations has been scarcely addressed (Aït-Sahalia and Lo 1998;Friedrich, Peinke, Sahimi, and Tabar 2011;Anvari, Tabar, Peinke, and Lehnertz 2016b). Moreover, proper numerical routines which cover all different kinds of parameteric contributions are still lacking. The main aim of this paper is to introduce and describe jumpdiff, a Python (Van Rossum et al. 2011) library, which is able to derive the evolution equation governing specific series of measurements from a jump-diffusion process without any a priori knowledge of the process.
In recent years, with the goal of going beyond the narrow scope of Gaussian processes, extensions of classical stochastic processes have included jumps, both with the practical aim of better describing the natural world, or because of the riddling mathematical hardship and the interesting results obtained. As a direct extension of the classical Langevin equation, jumpdiffusion processes came into focus via the Kramers-Moyal equation (Anvari et al. 2016b;Rydin Gorjão, Heysel, Lehnertz, and Tabar 2019). This family of processes embodies both a conventional Gaussian (diffusion) contribution as well as Poissonian jump contributions. From a time series analysis point-of-view, jump-diffusion processes are particularly hard to analyze and have required a revised interpretation under the classic Kramers-Moyal expansion (Kramers 1940;Moyal 1949). This hardship is given by the simultaneous presence of jumps and diffusion and required a revised interpretation under the classic Kramers-Moyal expansion, particularly when the goal is to extract separately the continuous and the discontinuous contributions from a data set. Jump-diffusion processes have found application in mathematical finance, particularly option pricing (Duffie, Pan, and Singleton 2000;Andersen, Benzoni, and Lund 2002;Johannes 2004; Aït-Sahalia and Jacod 2009), electricity markets (Cartea and Figueroa 2005), early-warning signal identification (Dakos et al. 2012), soil moisture dynamics (Daly and Porporato 2006), solar radiation and EEG recordings (Anvari et al. 2016a,b;Lehnertz, Zabawa, and Tabar 2018), and neural activity (Giraudo and Sacerdote 1997), amongst others.
In this paper we introduce the Python library jumpdiff, oriented towards non-parametric estimation of jump-diffusion processes via the Kramers-Moyal equation, which we will discuss shortly. While being a new library, it incorporates some of the recent theoretical and numerical achievements in stochastic methods (Anvari et al. 2016b;Tabar 2019). For example, Langevin processes have been implemented in an R (R Core Team 2022) package and an Python library, cf. Rinn, Lind, Wächter, and Peinke (2016) and Rydin Gorjão and Meirinhos (2019), respectively, enabling to model time series with drift and diffusion strengths. These packages can similarly non-parametrically estimate the strength of drift, diffusion, and other elements in data, but do not close the gap of directly evaluating and extracting the parameters of the analyzed time series. Other packages, such as Julia's (Bezanson, Edelman, Karpinski, and Shah 2017) DiffEqJump.jl (Rackauckas and Nie 2017, see https://github.com/SciML/DiffEqJump.jl), MATLAB's (The MathWorks Inc. 2021) PROJ_Option_Pricing_Matlab (Kirkby, Nguyen, Cui, Zhang, Deng, and Aguilar 2021), or C++'s DerivativesPricing (Gosain 2020) focus primarily on generating time series. Library jumpdiff allows not only to generate a series of values from a synthetic time series but also to recover all terms composing the evolution equation governing empirical sets of measurements from real jump-diffusion processes.
For the implementation of the jumpdiff library, we address two aspects in this paper. First, at the theoretical level, we introduce a set of higher-order finite-time corrections to the Kramers-Moyal equation, which are particularly relevant for non-parametric estimation of jumps, and in the presence of measurement noise (Böttcher, Peinke, Kleinhans, Friedrich, Lind, and Haase 2006;Lehle 2011Lehle , 2013Scholz et al. 2017). Second, at the numerical level, we implement a kernel density Nadaraya-Watson estimator (Nadaraya 1964;Watson 1964) to calculate the Kramers-Moyal coefficients, employing the aforementioned developed finite-time higher-order corrections in the library, and design a simple method to extract the jump contributions. Library jumpdiff enables the extraction of the parameters of the system in a non-parametric manner, as well as the distinction between pure diffusions and jump-diffusions. The library is suited for the family of Poissonian jump processes and represents a step forward towards more complex stochastic processes, namely Lévy-like processes (Siegert and Friedrich 2001;Lubashevsky, Friedrich, and Heuer 2009;Zaburdaev, Denisov, and Klafter 2015;Zan, Xu, Kurths, Chechkin, and Metzler 2020).
The paper is organized as follows: In Section 2 we introduce the main theoretical background behind jump-diffusion processes and outline the numerical approximation of Kramers-Moyal coefficients included in the jumpdiff library. In Section 3 a full description of all functions composing the library is presented, as well as a simple how-to-use guide. Section 4 concludes the paper, putting some aspects in perspective. Readers familiar with the theory and interested in the implementation and usage of the library jumpdiff only may directly go to Section 3. A more descriptive and hands-on explanation can be found in the documentation of the library at https://jumpdiff.readthedocs.io/, and the repository of the library is hosted at https://github.com/LRydin/jumpdiff.

Evolution and inference of jump-diffusion processes
In the following we introduce the stochastic differential equation describing a jump-diffusion process. We describe it together with its counterpart, the Kramers-Moyal equation, a partial differential equation that generalizes Langevin processes and can include discontinuous elements. Then we briefly explain how the terms, so-called coefficients, in the Kramers-Moyal equation link to the parameters of the jump-diffusion process. This is central for two reasons: First, the non-parametric technique employed in jumpdiff is based on the Kramers-Moyal equation. Second, we present a set of corrective terms to the estimation of the parameters via the Kramers-Moyal equation, which we include in jumpdiff to improve the estimation quality of the parameters.

Theoretical aspects behind jump-diffusion processes
In this section we consider the stochastic evolution of a time-continuous Markov process, X(t) ∈ R, that is governed by three independent contributions: one drift strength, one diffusive strength, and one Poissonian (jump) strength. The evolution equation of such a variable reads: where a(x, t) is the drift strength, b(x, t) is the diffusion or volatility, W (t) is a Wiener process, and J(t) is a time-homogeneous Poisson jump process with rate λ(x, t) and an amplitude ξ, which is normally distributed as ξ ∼ N (0, σ 2 ξ ). Jump-diffusion processes governed by (1) are a generalization of diffusion processes, since one recovers the latter when σ ξ = 0 (similarly λ(x, t) = 0). Below, we make use of the Kramers-Moyal expansion to connect the orders of the expansion with each parameter of the jump-diffusion process, thus enabling their estimation.
In order to link (1) to the evolution of the probability p(x, t+τ |x ′ , t), we employ the Kramers- where L KM denotes the Kramers-Moyal operator defined as (Risken 1996) Functions D m (x), called Kramers-Moyal coefficients, relate to the m-order conditional moment M m (x, τ ), given by For simplicity we drop the t-dependency focusing on stationary processes. The Kramers-Moyal coefficients D m (x) are thus defined for any integer m as The jump-diffusion process defined in (1) is linked to a particular case of the Kramers-Moyal expansion defined in (2) and (3), namely With this, we can invert the problem and ask instead the question: If we have a single realization of a stochastic process X(t), can we use the estimation of the Kramers-Moyal coefficients to uncover the parameters of the analyzed realization? The answer for jumpdiffusion processes is straightforward, and the parameters defining (1) are given by: Notice the relation of the parameters of the jump elements in (1) to the Kramers-Moyal coefficients in (5) is given by higher-order terms D 2m = (2m!) σ m ξ λ, for m ≥ 3 (Anvari et al. 2016b).
With this at hand, we can equate that estimating the drift strength a(x), the diffusion strength b(x), and the jump element ξ, translates to estimating the Kramers-Moyal coefficients. The job now is to provide an accurate estimation and to perform the inversion from Kramers-Moyal coefficients to the aforementioned parameters, which is the task of library jumpdiff.
It is important to notice that for purely diffusive processes, (2) reduces to the Fokker-Planck-Kolmogorov equation (sometimes denoted solely Fokker-Planck or Smoluchowski equation; Risken 1996), and the estimates of the two first Kramers-Moyal coefficients are sufficient to fully describe a diffusive process. In the case of jump-diffusion processes, higher-order Kramers-Moyal coefficients need to be taken into account as they are non-vanishing. For instance, for the jump-diffusion processes in (1), one can invert the problem from Kramers-Moyal coefficients to the stochastic parameters knowing the first six Kramers-Moyal coefficients (Anvari et al. 2016b).
Moreover, as we shall detail in the subsequent section, up to now, higher-order Kramers-Moyal coefficients were estimated from first-order estimations of the conditional moments. Here, we derive the second-order approximations, which imply a more cumbersome analytical approach to the Kramers-Moyal equation 2. These second-order approximations are crucial to improve the estimation of jump-diffusion processes.

Numerical computation of the Kramers-Moyal coefficients
The numerical computation of the Kramers-Moyal coefficients is based on the numerical computation of conditional moments M m (x, t), which can be estimated directly from a set of measurements X(t). Indeed, the instantaneous time rate of the moment of order m for the process X(t), conditioned to a specific value x, is given by with ⟨X(t)⟩ denoting the average of X(t), for all measured t.
The conditional moments can be expressed as sums of products of Kramers-Moyal coefficients, derived from the formal solution of (2), namely Depending on the number of terms used from the sum in (9) one obtains different orders of approximation of the Kramers-Moyal operator. Here we will consider first-and second-order approximations.

First-order approximation
The first-order approximation is given by where expressing the Kramers-Moyal coefficients from (5), via the moments in (4), yields The full derivation is given in Appendix A. The derivation is not cumbersome and well known to the community 1 , thus of little interest to show here.

Second-order approximation
Higher-order approximations of the conditional moment are especially relevant when handling low-sampled data. The second-order approximation of the conditional moments takes in one additional term from the sum in (9), namely where naturally the first-order terms ∼ τ are present at the first-order of the approximation. The derivation, found in full detail in Appendix A, involves a set of relations between the moments M m (x, τ ) of a given order m and the Kramers-Moyal coefficients D n (x) of all orders 0 < n ≤ m. Inverting these relations with respect to the Kramers-Moyal coefficients yields where F m (x, τ ) denotes the second-order approximation, in comparison with the first-order approximation given above in (11). The second-order approximations F m (x, τ ) are given by (dependencies removed for clarity) Here naturally the first term on each right-hand side is the first-order approximation. This second-order approximation neglects terms including derivatives of the Kramers-Moyal coefficients, which enables one to express the nth Kramers-Moyal coefficient as a function of conditional moments up to order n − 1. In this way, we provide a general formula for improving the estimates of Kramers-Moyal coefficient, taken as linear approximations of the corresponding conditional moment. The full derivation is given in Appendix A.

Implementation of the functions in library jumpdiff
In this section we introduce jumpdiff and some of its main functions and utilities. The simplest way to install jumpdiff is via the Python Package Index (PyPI), simply using $ pip install jumpdiff or with the readers' favorite Python package installation program. The library has three main functionalities for addressing jump-diffusion processes, namely: • Generate sample trajectories of jump-diffusion processes. • Extract all parameters for a single trajectory of a jump-diffusion process.
• Evaluate if a trajectory is a purely diffusion process or a jump-diffusion process.
We note here that the user interested in evaluating their own data should start with the last mentioned step, i.e., to firstly evaluate if their data is described by a jump-diffusion process.
In this section, we describe in detail each one of these functionalities, explaining how we employ the most important functions in jumpdiff to estimate the parameters of a jumpdiffusion process. All functions included in library jumpdiff are listed in Table 1.

From equation to data: Generating sample trajectories
To start off, import the library into a working Python environment, using

>>> import jumpdiff as jd
Subsequently, all functions of jumpdiff can be access via jd. To generate a sample trajectory of a jump-diffusion process, call the function jd_process(). We will generate a single trajectory of the process, and subsequently employ the non-parametric estimators in jumpdiff to retrieve the parameters of the jump-diffusion process generated, described in the following subsections. Let us take a simple example and generate a trajectory of an Ornstein-Uhlenbeck process with Poissonian jumps, given by equation Naturally we also need to specify the values of θ, the drift strength, σ, the diffusion strength, and σ ξ and λ, the standard deviation of ξ and the jump rate of J(t). We can use jd_process() for integrating a jump-diffusion process, with a number of points N = 1×10 7 (t = 10 4 ) and a time-step of ∆t = 0.001. This is implemented in the following way: First, specify the integration time and time sampling: >>> time = 10000 >>> delta_t = 0.001

Then define the drift function a(x) and the diffusion function b(x):
>>> def a(x): ... return -0.5*x >>> def b(x): ... return 0.75 Define jump amplitude and rate >>> xi = 1.5 >>> lamb = 1.75 and generate the jump-diffusion process: >>> X = jd.jd_process(time, delta_t, a, b, xi, lamb) In Figure 1 we plot the the trajectory generated with this code alongside with a zoomed region where a jump can be distinctly seen.
While in this example an Ornstein-Uhlenbeck process, i.e., a linear drift strength and constant diffusion strength, is chosen with a Poissonian jump strength, other higher polynomial

From data to the jump-diffusion equation
Having described how to generate series of values from a jump-diffusion equation, using the function jd_process(), we now consider the inverse problem: Starting from that series of values, derive the parameters of the jump-diffusion equation. To that end, we extract the Kramers-Moyal coefficients D m (x) of this exact simulated process X, such that we can employ the non-parametric estimation method and compare if the parameters extracted match the parameters chosen.
We now employ the function moments() on the generated time series X to extract the Kramers-Moyal coefficients. This we achieve using: First, we extract the Kramers-Moyal coefficients and state-space without second-order corrections: >>> edge, simple_mom = jd.moments (timeseries=X, correction=False) and then with second-order corrections: >>> edge, mom = jd.moments(timeseries=X, correction=True) In the case above we include the estimation of the Kramers-Moyal coefficients both with and without the corrective terms we designed. Figure 2 shows the derived results for the inverse problem, using first-order (dashed lines) and second-order corrections (solid lines).
The theoretical values are indicated by the dotted lines.
In order to showcase the importance of the corrective terms (correction = True), we show, for the generated process X, the estimation of the Kramers-Moyal coefficients while downsampling data. We consider different sampling rates, namely s f = 1, 0.05, and 0.01, i.e., taking the full data, only every 20th datapoint, and only every 100th datapoint. As can be seen in Figure 2, the second-order corrective term improves the estimate of all Kramers-Moyal coefficients D m (x), m ≥ 2, especially in the cases with lowest sampling rates. We can see in Figure 2 that the function moments() has estimated up to six Kramers-Moyal coefficients D n (x). We will now turn -after a short paragraph -to two specific functions in jumpdiff that allow us to directly extract the jump amplitude σ 2 ξ and the jump rate λ, as these are the contributions from the discontinuous part to the process.
Before we turn to the estimation of the jump elements, note that technically, the function moments() performs a kernel-density estimation employing a Nadaraya-Watson estimator (Nadaraya 1964;Watson 1964) of the different orders of the moments. By default it employs an Epanechnikov kernel and uses a convolution method to perform the kernel-density estimation of the moments. Four other kernels are also available, namely Gaussian, uniform, triangular, and quartic kernels. Moreover, the function moments() includes the input parameter lag, which is a time-lag that enables one to estimate the conditional moments at different time-lags τ . If left unspecified, it assumes the shortest increment of the time series, i.e., the time series sampling rate τ = 1/s f . This is especially suited for evaluating the conditional moments in the limiting case of τ → 0, since numerical accuracy is bounded by the sampling rate s f . This evaluation is done by plotting a few time-lags, e.g., τ = 1/s f , . . . , 10/s f , and extrapolate the limit τ → 0 (Böttcher et al. 2006;Lind, Haase, Böttcher, Peinke, Kleinhans, and Friedrich 2010).
In Figure 2 we can see the estimation of the Kramers-Moyal coefficients, with first-and secondorder corrections. The results reproduce well the theoretical values (dotted lines), allowing us to extract the parameters of our process via (7). Retrieving the drift and diffusion functions, a(x, t) and b(x, t) respectively, are known problems, which imply studying the first and second Kramers-Moyal coefficients. From (7) one recovers also the jump amplitude σ 2 ξ and the jump rate λ of the time series, by considering higher-order conditional moments. In jumpdiff, functions jump_rate() and jump_amplitude() implement (7d) and (7c), respectively. To estimate the jump amplitude σ 2 ξ and the jump rate λ of the time series, use >>> xi_est = jd.jump_amplitude(moments=simple_mom) to estimate the jump amplitude and >>> lamb_est = jd.jump_rate(moments=simple_mom) to estimate the jump rate. This yields a numerical estimation of each of the parameters. The user can as well utilize mom instead of simple_mom.
Here we leave an important note for the user, regarding the impact of a low number of data points in a time series. The estimation of the jump amplitude σ 2 ξ and the jump rate λ is strongly dependent on the number of data points in each time series. In Figure 3 we numerically integrate (15), as before, and employ the estimators of the jump amplitude and jump rate, i.e., functions jump_rate() and jump_amplitude(), for a time series with increasing number of data points N . The estimators converge to the theoretical values (dashed line) with increasing accuracy with the average number on jumps ⟨λ⟩ in the time series. It is central to understand that short time series cause an under-estimation of σ 2 ξ and simultaneously an −0.5x 0.75 1.5 1.75 Estimated −0.496x 0.760 1.524 1.802 Table 2: Parameter estimation for the process (15) in Figure 1, generated with jd_process(), with indicated parameters. first-order second-order 10 3 10 4 10 5 10 3 10 4 10 5 Figure 3: Estimation of the jump amplitudeσ 2 ξ (left) and the jump rateλ (right) from 20 numerically simulated jump-diffusion processes with increasing time lengths N . The data was generated by integrating (15), with the σ 2 ξ = 1.5 and λ = 1.75 (dashed lines), and each term was estimated using the functions jump_amplitude() and jump_rate(). The estimates are shown as a function of the number of points N ∈ [1×10 5 , 5 × 10 7 ], always with a time-step ∆t = 0.001. Top axis denotes the average number of jumps ⟨λ⟩ in the respective numerically integrated time series. Each point is an average over 10 iterations. Standard deviations depicted in the shaded areas.
over-estimation of λ. This also springboard us to the next important question: How can we know our time series is a jump-diffusion?
Before we move on, we can compare the theoretical values used in (15) and estimated parameters via fitting the slopes of D 1 (x) and the constant values of D 2 (x) in Figure 2, as well as the used functions jump_rate() and jump_amplitude(). The values can be found in Table 2 and agree with our constructed model. For a general method to recover all Kramers-Moyal coefficients, for stochastic processes of any dimension, see Rydin Gorjão and Meirinhos (2019).

Distinguishing between purely diffusive and jump-diffusion processes
As mentioned at the start of this section, when analyzing real measurements or observational data, the following steps are crucial to understand whether the data at hand fall into the category of a jump-diffusion process or not.
The presence of a discontinuous term in (1) is the fundamental addition to a general diffusion process. However, assuming such Ansatz for purely diffusive process may lead to spurious jump terms. To avoid this, one fundamental question to ask is whether we are in the presence of a pure diffusion or a jump-diffusion process, in order to choose ab initio the proper Ansatz.
To differentiate diffusive from jump-diffusion processes, Lehnertz et al. (2018) introduced a simple criterion, which is based on the scaling of the ratio of the fourth-and sixth-order conditional moments with the scale τ , denoted the Q-ratio, given by If the process is purely diffusive Q(x, τ ) = τ (b(x)) 2 (linear function in τ ), whereas if the process has a jump term, Q(x, τ ) = σ 2 ξ (constant, independent of τ ). This criterion can be employed directly for any time series and is implemented as the function q_ratio() in jumpdiff.
To analyze the scaling of the Q-ratio, follow the recipe: First import numpy for generating logarithmic spaced arrays and matplotlib for plotting in double logarithmic scale >>> import numpy as np >>> import matplotlib.pyplot as plt Then take a sequence of integer lags >>> lag = np.logspace (0, 4, 25, dtype=int) and recover the Q-ratio of the time series X using >>> lag, Q = jd.q_ratio(lag, X) before plotting in a log-log scale >>> plt.loglog(lag, Q) Figure 4 showcases how the Q-ratio changes between purely diffusive to jump-diffusion processes. In this manner, we consider the same process we have introduced in (15), but we will vary the the amplitude of the jumps σ 2 ξ ∈ [0, 1]. Figure 4 illustrates the implementation of the function q_ratio(), plotting it for several time-lags τ in a log-log plot. If there is a linear dependence on τ the time series X(t) is a pure diffusive process, whereas if the plot is approximately flat, showing a constant Q-ratio the time series should have a jump term. Notice that increasing the jump rate from ξ = 0 (pure diffusion) to ξ = b(x) = 1 (same amplitude as diffusion strength), the Q-ratio changes from linearly depending on τ to a constant. The different relation between the Q-ratio and τ should help the user distinguish which type of data they are analyzing, and thus help assert whether their data follows a purely diffusive process or is better described by a jump-diffusion.

Discussion and conclusions
Jump-diffusion processes are stochastic models able to describe processes riddled with jumps. Alongside with regular diffusion processes, their elegance lies in the possibility of estimating, non-parametrically, the parameters via the Kramers-Moyal equation. Access to higher-order moments of the Kramers-Moyal equation is computationally feasible, and thus permits a oneto-one correspondence between model parameters and the estimators (cf. (7)). These, on the other hand, are hampered by the scarcity of jumps. Where the diffusive strengths are ubiquitous, the jumps take place sparsely in time. To this effect, the set of second-order corrections to the Kramers-Moyal operator implemented here are crucial in improving the estimation of the parameters of the model.
The use of jump-diffusion models as descriptives of processes beyond diffusions have seen application across different research fields. The elegance of these processes lies not only in their general applicability, but also in the ease of using non-parametric parameter estimation. In this sense, the presented second-order corrections are of substantial relevance, especially taking into account that jumps occur sparsely. Here we present a two-fold project, developing second-order corrections of the Kramers-Moyal expansion of jump-diffusion processes and designing an easy-to-employ Python library. jumpdiff requires minimal mathematical knowledge to be used, and thus can find application even in non-mathematically oriented research fields.
The limitations and finite scope of the library here described motivate further development, namely towards more general jump types, having amplitudes and rates which are time (and space) dependent. Moreover, beyond Poissonian jumps, new libraries can be developed, for instance to approach Lévy-like processes.

A. Second-order corrections of Kramers-Moyal coefficients
For the first-order approximation of conditional moments M n , we substitute in (4) the firstorder approximation of the exponential operator, namely The last integral is given by: This approximation is the one used in Anvari et al. (2016b).
For the second-order approximation of conditional moments M n , we substitute in (4) the second-order approximation of the exponential operator, namely Finally, inverting (25) yields (14).
The formulas for the relations (25) and the corrections F n given by (14) are given in symbolic Python by the functions m_formula and f_formula, respectively. This numerical procedure is generalized to any maximal order N needed for the data analysis. For jump-diffusion processes N = 6 is sufficient.

B. Higher-order corrections and Kramers-Moyal coefficients
A more accurate approximation is possible by combining (23) and (24). More precisely, the steps to solve numerically (23) are the following ones: • Solve (14) introducing the conditional moments, extracted directly from the data, up to order N , and derive the Kramers-Moyal coefficients D n (x) as in (13).
• Compute the derivatives up to order N d (see discussion below).
• Introduce the derivatives of the Kramers-Moyal coefficients and the empirical conditional moments in (23) and solve it with respect to the Kramers-Moyal coefficients.
• Repeat steps S1 and S2 until the Kramers-Moyal coefficients converge within a pre-given numerical accuracy.
Moreover, since the Kramers-Moyal coefficients are typically polynomials of lower order, not larger than five or six, the derivative order N d is a finite number, which leads to a simplification of (23). Namely, the derivative in the sum obeys 0 ≤ p + m − n ≤ N d . Thus, p 0 ≤ p ≤ N d + n − m. Since p 0 = max (1, n − m) one has N d + n − m ≥ 1 and therefore the sum over m is bounded by 1 ≤ m ≤ N d + n − 1. Since N d > 1, p 0 = 1 and therefore the sum over p is also bounded by 1 ≤ p ≤ N d + n − m.
Introducing these bounds in (23) and following steps S0-S3 above, yields the second-order approximation of the Kramers-Moyal coefficients.
In the case of (26) we have k = 2.
The ordinary Bell's polynomials fulfill a reciprocal relation, namely any sum of ordinary Bell's polynomials of the form y n = n k=1 B n,k (x 1 , . . . , x n−k+1 )