Analysis and data-driven reconstruction of bivariate jump-diffusion processes
Leonardo Rydin Gorjão, Jan Heysel, Klaus Lehnertz, M. Reza Rahimi Tabar
AAnalysis and data-driven reconstruction of bivariate jump-diffusion processes
Leonardo Rydin Gorj˜ao,
1, 2, 3, 4, ∗ Jan Heysel,
1, 2, † Klaus Lehnertz,
1, 2, 5, ‡ and M. Reza Rahimi Tabar
6, 7, § Department of Epileptology, University of Bonn, Venusberg Campus 1, 53127 Bonn, Germany Helmholtz Institute for Radiation and Nuclear Physics, University of Bonn, Nussallee 14–16, 53115 Bonn, Germany Forschungszentrum J¨ulich, Institute for Energy and Climate Research - SystemsAnalysis and Technology Evaluation (IEK-STE), 52428 J¨ulich, Germany Institute for Theoretical Physics, University of Cologne, 50937 K¨oln, Germany Interdisciplinary Centre for Complex Systems, University of Bonn, Br¨uhler Straße 7, 53175 Bonn, Germany Institute of Physics and ForWind, Carl von Ossietzky University of Oldenburg,Carl-von-Ossietzky-Straße 9–11, 26111 Oldenburg, Germany Department of Physics, Sharif University of Technology, 11365-9161, Tehran, Iran
We introduce the bivariate jump-diffusion process, comprising two-dimensional diffusion and two-dimensional jumps, that can be coupled to one another. We present a data-driven, non-parametricestimation procedure of higher-order (up to 8) Kramers–Moyal coefficients that allows one to re-construct relevant aspects of the underlying jump-diffusion processes and to recover the underlyingparameters. The procedure is validated with numerically integrated data using synthetic bivariatetime series from continuous and discontinuous processes. We further evaluate the possibility of es-timating the parameters of the jump-diffusion model via data-driven analyses of the higher-orderKramers–Moyal coefficients, and the limitations arising from the scarcity of points in the data ordisproportionate parameters in the system.
I. INTRODUCTION
Research over the last two decades has demonstratedthe high suitability of the network paradigm in advanc-ing our understanding of natural and man-made complexdynamical systems [1–7]. With this paradigm, a systemcomponent is represented by a vertex and interactionsbetween components are conveyed by edges connectingvertices, and graph theory provides a large repertoire ofmethods to characterize networks on various scales.Characterizing properties of interactions using theknowledge of the dynamics of each of the componentsis key to understanding real-world systems. To achievethis goal, a large number of time-series-analysis meth-ods has been developed that originate from synchroniza-tion theory, nonlinear dynamics, information theory, andfrom statistical physics (for an overview, see Refs. [8–15]).Some of these methods make rather strict assumptionsabout the dynamics of network components generatingthe time series and many approaches preferentially focuson the low-dimensional deterministic part of the dynam-ics.Real-world systems, however, are typically influencedby random forcing and interactions between constituentsare highly non-linear, which results in very complex,stochastic, and non-stationary system behavior that ex-hibits both deterministic and stochastic features. Aimingat determining characteristics and strength of fluctuat-ing forces as well as at assessing properties of non-linearinteractions, the analysis of such systems is associated ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] with the problem of retrieving a stochastic dynamicalsystem from measured time series. There is a substan-tial existing literature [16–19] for the modeling of com-plex dynamical systems which employs the conventionalLangevin equation that is based on the first- and second-order Kramers–Moyal (KM) coefficients, known as driftand diffusion terms. All functions and parameters of thismodeling can be found directly from the measured timeseries employing a widely used non-parametric approach.There are by now only few studies that make use ofthis ansatz to characterize interactions between stochas-tic processes [20–24].Despite its successful application in diverse scien-tific fields, growing evidence indicates that the contin-uous stochastic modeling of time series of complex sys-tems (white-noise-driven Langevin equation) should ac-count for the presence of discontinuous jump compo-nents [19, 25–33]. In this context, the jump-diffusionmodel [34–37] was shown to provide a theoretical toolto study processes of known and unknown nature thatexhibit jumps. It allows one to separate the determin-istic drift term as well as different stochastic behaviors,namely diffusive and jumpy behavior [19, 32, 33]. More-over, all of the unknown functions and coefficients ofa dynamical stochastic equation that describe a jump-diffusion process can be derived directly from measuredtime series. This approach involves estimating higher-order ( ≥
3) KM coefficients and it provides an intuitivephysical meaning of these coefficients.The focus of this paper is to introduce a method toinvestigate bivariate time-series with discontinuous jumpcomponents. We begin with an overview of bivariate dif-fusion processes that exhibit the known relation betweenthe parameters and functions in stochastic modeling andthe KM coefficients. Exemplary processes are portrayed,and we propose a measure to judge the quality of our re- a r X i v : . [ n li n . AO ] S e p construction procedure. We then present bivariate jump-diffusion processes alongside with the associated KM ex-pansion [32], and investigate the suitability of our recon-struction procedure using various examples. We concludethis paper by summarizing our findings. II. STATE OF THE ART: BIVARIATEJUMP-DIFFUSION MODEL
A bivariate jump-diffusion process can be modeled via[19, 32] y (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) d y ( t )d y ( t ) (cid:19) = N (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) N N (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) drift d t + g (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) g , g , g , g , (cid:19) d w (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) d w d w (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) diffusion + ξ (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) ξ , ξ , ξ , ξ , (cid:19) d J (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) d J d J (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) Poissonian jumps , (1)where all the elements of vectors N , d J , and d w aswell as of matrices g and ξ may, in general, be state-and time-dependent (dependencies not shown for con-venience of notation). The drift coefficient is a two-dimensional vector N = ( N , N ) with N ∈ R , whereeach dimension of N , i.e., N i , may depend on y ( t ) and y ( t ). The diffusion coefficient takes a matrix g ∈ R × .The two Wiener processes w = ( w , w ) act as inde-pendent Brownian noises for the state variables y ( t )and y ( t ). The diagonal elements of g comprise the dif-fusion coefficients of self-contained stochastic diffusiveprocesses, and the off-diagonal elements represent in-terdependencies between the two Wiener processes, i.e.,they result from an interaction between the two pro-cesses. Each single-dimensional stochastic process ele-ment d w i of d w is an increment of a Wiener process,with (cid:104) d w i (cid:105) = 0 , (cid:104) d w i (cid:105) = d t, ∀ i . The discontinuous jumpterms are contained in ξ ∈ R × and d J ∈ N , where d J represents a two-dimensional Poisson process. Theseare Poisson-distributed jumps with an average jump rate λ ∈ R in unit time t . The average expected number ofjumps of each jump process J i in a timespan t is λ i t . Thejump amplitudes ξ are Gaussian distributed with zeromean and standard deviation ξ i,j .We here consider merely autonomous systems. Non-ergodic problems are beyond the scope of this work, anda more delicate approach to both bivariate stochastic pro-cesses would be needed. III. BIVARIATE DIFFUSION PROCESSES
Let us begin with bivariate diffusion processes, forwhich the model takes the form y (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) d y ( t )d y ( t ) (cid:19) = N (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) N N (cid:19) d t (cid:124) (cid:123)(cid:122) (cid:125) drift + g (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) g , g , g , g , (cid:19) d w (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) d w d w (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) diffusion . (2)The model consists of six functions, two for the drift co-efficients and four for the diffusion coefficients. Given abivariate diffusion process, can we reconstruct the afore-mentioned parameters strictly from data? Extensive workexists on this matter [18], especially covering purely dif-fusion processes, and we will use these now as a step-ping stone to jump-diffusion processes. Understandingthe working and contingencies of reconstructing the pa-rameters of a diffusion process (Eq. (2)) will serve as agateway to understand how a similar procedure awardsequivalent measures for jump-diffusion processes. We ad-dress the aforementioned question firstly by revisitingthe mathematical foundation that allows one to recover,strictly from data, the drift N and diffusion g coeffi-cients. Subsequently, we numerically integrate diffusionprocesses with a priori fixed values of the drift N anddiffusion g coefficients and aim at retrieving these val-ues strictly from the generated data (Euler–Mayuramascheme with a time sampling of 10 − over a total of 10 time units, i.e., 10 number of data points). If the actualand retrieve values match, the reconstruction method iseffective.A stochastic process has a probabilistic descriptiongiven by the master equation [16, 19]. It does not de-scribe a specific stochastic processes in itself, but theprobabilistic evolution of the process in time. The mas-ter equation accepts an expansion in terms, the KM ex-pansion, that allows for a purely differential descriptionof the process. More importantly, the coefficients of theexpansion, known as KM coefficients, entail directly arelation to the aforementioned parameters of a stochas-tic process given by Eq. (1). The exact relation will begiven below. There is although an important detail inwhat regards the KM coefficients: they are in themselvesnot constants, but functions on the underlying space, orin other words, a scalar field, and for our purposes here,they can be understood as two-dimensional surfaces. Wewill denote these as KM surfaces.Lastly, and more familiar, the Fokker–Planck equationis a truncation of the KM expansion at second order.It is especially relevant given its connection to physicalprocesses and the Pawula theorem [38]. The Pawula the-orem ensures that the truncation is not ill-suited for theunderlying process if the fourth-order KM coefficient ap-proaches zero in the limit d t →
0. It is now crucial tounderstand that the theorem holds for a one-dimensionalprocess, and we are not aware of a proof for higher di-mensions. This contrasts the common notion that study- -0.8-0.40.00.40.8-0.4-0.20.00.20.4-1.0-0.50.00.51.0 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.4-0.20.00.20.4-0.4-0.20.00.20.4 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.4-0.20.00.20.40.00.0040.0080.0120.016 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.4-0.20.00.20.40.120.140.16 M [2 , and f [2 , y y FIG. 1. Two-dimensional Kramers–Moyal coefficients M [ (cid:96),m ] for two independent diffusion processes given by Eq. (8). Theuncovered KM surfaces match the expectation, the cubic term in the drift term N = − x + x along the first variable isvisible on the top-left panel, M [1 , , and the negative-slope surface in the top-right panel, M [0 , . The flat surfaces reproduceas well the expected form of the constant terms involved in the diffusion terms for M [2 , on the bottom-right. Moreover,the bottom-left panel M [1 , , which accounts for the couplings terms of all diffusion terms, is also zero almost everywhere, asexpected, given g , and g , are zero. In each panel, the theoretical expected surface, given by Eqs. (4) and (5), is indicatedby a grid, with f [ (cid:96),m ] denoting the respective theoretical values introduced in the model. ing only the first two KM coefficients of two- or higher-dimensional processes is sufficient (see Refs. [32, 33, 39]and references therein).The KM coefficients M [ (cid:96),m ] ( x , x ) ∈ R of orders( (cid:96), m ) are defined as: M [ (cid:96),m ] ( x , x ) =lim ∆ t → t (cid:90) ( y ( t +∆ t ) − y ( t )) (cid:96) ( y ( t +∆ t ) − y ( t )) m · P ( y , y ; t +∆ t | y , y ; t ) | y ( t )= x ,y ( t )= x d y d y , and can be obtained from bivariate time series( y ( t ) , y ( t )). Theoretically, ∆ t should take the limit-ing case of ∆ t →
0, but the restriction of any mea-suring or storing devices—or the nature of the observ-ables themselves—permits only time-sampled or discreterecordings. The relevance and importance of adequatetime-sampling was extensively studied and discussed inRef. [19, 33].In the limiting case where ∆ t is equivalent to the sam- pling rate of the data, the KM coefficients take the form M [ (cid:96),m ] ( x , x ) = 1∆ t (cid:104) ∆ y (cid:96) ∆ y m | y ( t )= x ,y ( t )= x (cid:105) , (3)∆ y i = y i ( t + ∆ t ) − y i ( t ) . The algebraic relations between the KM coefficientsand functions in Eq. (2) are given by [19, 32] M [1 , = N , M [0 , = N , (4) M [1 , = g , g , + g , g , , M [2 , = (cid:2) g , + g , (cid:3) , M [0 , = (cid:2) g , + g , (cid:3) . (5)An explicit derivation can be found in Appendix A.Evidently, this under-determined set of five equationsis insufficient to uncover the six functions of a generalstochastic diffusion process. One must bare this in mind,for the same issue will arise again when reconstructingjump-diffusion processes from data. Nonetheless, under -0.8-0.40.00.40.8-0.2-0.10.00.10.20.050.10.150.2 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.00.20.40.60.8 M [2 , and f [2 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.00.020.040.06 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.00.0020.004 M [2 , and f [2 , y y FIG. 2. Two-dimensional Kramers–Moyal (KM) coefficients M [ (cid:96),m ] for two independent diffusion processes given by Eq. (9) andthe theoretical expected functions f [ (cid:96),m ] associated with each coefficient, according to Eq. (5). KM coefficients M [1 , , M [2 , , M [0 , , and M [2 , exhibit the quadratic multiplicative dependencies of the diffusion terms. In addition, M [1 , displays bothan offset from zero as well as a quadratic shape, entailing the desired results emerging from Eq. (5), i.e., the noise-couplingterm g , with g , . M [2 , displays an offset and has a minimum close to g , / g , / .
13. We also show the higher-ordercoefficient M [2 , and the corresponding theoretically expected value (given by Eq. (5)), both of which vanish. All obtained KMsurfaces fit considerably well their theoretically expected ones ( V [1 , = 0 . V [2 , = 0 . V [0 , = 0 . V [2 , < .
01; errorvolumes estimated over the displayed domain). certain assumptions it is possible to reduced the dimen-sion of the problem and therefore obtain a system ofequation which is not under-determined. Two methodsfor these cases are present in Ref. [19] and another crite-rion will be presented later.In order to relate the results obtained from studyingthe KM coefficients against the theoretical functions, wepropose a method to assess the difference between thevalues of the theoretically expected functions and the es-timated values of the KM coefficients. Since for bivariateprocesses the KM coefficients are two-dimensional—asare the parameters of Eq. (1)—an adequate “distance”measure between the resulting two-dimensional surfacesis required.Following Ref. [20], we propose a distance measure thatallows for the variability of the density of data in someregions of the underlying space to be taken into consid-eration.Let f [ (cid:96),m ] ( y , y ) denote the theoretical value for orders( (cid:96), m ) introduced in the model, i.e., a non-linear com-bination of the various parameters of the system. The distance between each surface can be defined as (cid:90)(cid:90) U (cid:16) M [ l,m ] ( y , y ) − f [ l,m ] ( y , y ) (cid:17) d y d y =: V , (6)where U denotes the domain of M [ l,m ] ( y , y ). The least-squared distance volume V between the surfaces is zeroif M [ l,m ] ( y , y ) = f [ l,m ] ( y , y ). It is this volume thatone aims to minimize such that the reconstructed KMcoefficients match the underlying theoretical functions inthe model. Since M [ l,m ] ( y , y ) is a real-valued functionmeasured over a distribution space U , the density of datapoints is not uniform over U . This implies that a compar-ative measure on distances between M [ l,m ] ( y , y ) and f [ l,m ] ( y , y ) would be non-normalized to the density ofpoints of the space. We therefore introduce a normaliza-tion to Eq. (6) that ensures the less dense areas of U are normalized accordingly, thus mitigating the effect ofscarcity of points at the borders of U and an overestima-tion of V due to outliers in the distribution. We derivesuch a normalization by considering the zeroth-order KM -0.8-0.40.00.40.8-0.2-0.10.00.10.2-1.0-0.50.00.51.0 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.2-0.4-0.20.00.20.4 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.050.10.150.2 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.10.20.30.4 M [2 , and f [2 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.040.080.12 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.00.0010.0020.003 M [2 , and f [2 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.2-0.0020.00.0020.0040.006 M [4 , and f [4 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.2-0.0020.00.0020.0040.006 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.2-0.10.00.10.20.00.0040.0080.012 M [4 , and f [4 , y y FIG. 3. Two-dimensional Kramers–Moyal (KM) coefficients M [ (cid:96),m ] of bivariate diffusion processes given by Eq. (15) (with φ = 0 .
0) together with the theoretically expected functions f [ (cid:96),m ] associated with each coefficient according to Eqs. (10),(11), and (12). KM coefficients M [1 , , M [0 , , M [1 , , M [2 , , M [0 , , M [2 , , M [4 , , M [0 , , and M [4 , are shown. Althoughseemingly small, the higher-order moments are all present and non-zero. We find M [4 , = 0 .
012 and M [0 , = 0 . V [1 , = 0 . V [0 , = 0 . V [1 , < . V [2 , = 0 . V [0 , = 0 . V [2 , < . V [4 , = 0 . V [0 , = 0 . V [4 , < . coefficient M [0 , ( y , y ) which captures exactly the den-sity of points in U , although it is in itself not normalizedas a distribution. The resulting normalized volume er-ror measure V err between surfaces takes the form (statedependencies not explicit) (cid:90)(cid:90) U (cid:16) M [ l,m ] − f [ l,m ] (cid:17) p ( y , y )d y d y = V , (7)where p ( · ) denotes the probability density. Coinciden-tally, the numerical evaluation implemented via botha histogram or a kernel-based estimators immediatelyyields this density. (i.e., the zeroth-power of the right-hand side of Eq. (3), before applying the estimation op-erator). This makes it easy to retrieve p ( y , y ) as onenumerically evaluates data.With this at hand, it is now possible to relate theoret-ical and numerical results and to quantify the deviationof the obtained KM coefficients from the functions em-ployed.To showcase what two-dimensional KM coefficients areas well as how to identify drift and diffusion terms of bivariate diffusion processes, we present in the follow-ing two exemplary processes with a priori known coef-ficients. In this manner, by employing Eqs. (4) and (5),one can judge the outcome of the KM coefficient estima-tion procedure from discrete data in comparison with theexpected theoretical functions.We begin with two uncoupled processes, where one hasconstant diffusion and a quartic potential as the driftterm: N = (cid:18) N N (cid:19) = (cid:18) − x + x − x (cid:19) , g = (cid:18) g , g , g , g , (cid:19) = (cid:18) . . . . (cid:19) , (8)In Fig. 1, we show the corresponding KM coefficients M [1 , , M [0 , , M [1 , , and M [2 , together with the theo-retically expected functions. The per-design cubic-linearfunction ( N = − x + x ) acting as drift coefficient alongthe first dimension as well as the negatively-sloped sur-face of N = − x are evident. Likewise, the constant dif-fusion term leads to a flat constant-valued M [2 , , and theabsence of any non-diagonal elements ( g , = g , = 0)agrees with the zero-valued M [1 , . Alongside the sur-faces are plotted the theoretically expected values, whichagree well with the data recovery.We next extent Eq. (8) by adding multiplicative noiseto the diffusion term and by including a noise couplingterm g , (cid:54) = 0 N = (cid:18) N N (cid:19) = (cid:18) − x + x − x (cid:19) , g = (cid:18) g , g , g , g , (cid:19) = (cid:18) . x . . . x (cid:19) . (9)The recovered KM coefficients (see Fig. 2) of the driftterms remain unaltered, but as posited, the second-orderKM coefficients, i.e., M [1 , , M [2 , , M [0 , , and M [2 , ,clearly exhibit the influence of the multiplicative noise.The quadratic multiplicative dependencies of M [2 , and M [0 , and their offsets from zero are evident. More per-tinently, one can notice M [1 , to display the expectedshape arising from Eq. (5), i.e., this value is non-zero andexhibits the parabolic shape of g , · g , = 0 . . x ).For x = 0, the minimum of M [1 , coincides with 0 .
1, asexpected. This indicates that the presence of the multi-plicative noise does not hinder the assertion of the KMcoefficients. Again, the recovered KM coefficients matchthe theoretically ones.
IV. BIVARIATE JUMP-DIFFUSIONPROCESSES
The KM coefficients of bivariate jump-diffusion pro-cesses take the following form (under the parameter pre-scription used in the jump-diffusion model in Eq. (1)) [19,32, 33] M [1 , = N M [0 , = N (10) M [1 , = g , g , + g , g , M [2 , = (cid:2) g , + s , λ + g , + s , λ (cid:3) M [0 , = (cid:2) g , + s , λ + g , + s , λ (cid:3) (11) M [2 , = [ s , s , λ + s , s , λ ] M [4 , = 3 (cid:2) s , λ + s , λ (cid:3) M [0 , = 3 (cid:2) s , λ + s , λ (cid:3) M [4 , = 9 (cid:2) s , s , λ + s , s , λ (cid:3) (12) M [6 , = 15 (cid:2) s , λ + s , λ (cid:3) M [0 , = 15 (cid:2) s , λ + s , λ (cid:3) M [6 , = 225 (cid:2) s , s , λ + s , s , λ (cid:3) (13) M [8 , = 105 (cid:2) s , λ + s , λ (cid:3) M [0 , = 105 (cid:2) s , λ + s , λ (cid:3) M [8 , = 11025 (cid:2) s , λ + s , λ (cid:3) (14)where (cid:104) ξ (cid:96)ij (cid:105) = s (cid:96)ij are the variances of the Gaussian-distributed jump amplitudes. An extended derivation canbe found in Appendix A. The last equations here aretaken from the general form M [2 (cid:96), m ] = (2 (cid:96) )!2 (cid:96) (cid:96) ! (2 m )!2 m m ! (cid:2) s (cid:96) , s m , λ + s (cid:96) , s m , λ (cid:3) . A. Understanding the impact of jumps
As an illustrative case study, we investigate a generaljump-diffusion process that is based on Eq. (9) but ex-cludes the multiplicative diffusion terms. Taking into ac-count the effect of the jump terms, but maintaining thesystem independent in at least one of the dimensions,we extend Eq. (9) to include jumps only in the diagonalterms of ξ N = (cid:18) N N (cid:19) = (cid:18) − x + x − x (cid:19) g = (cid:18) g , g , g , g , (cid:19) = (cid:18) . . . . (cid:19) ξ = (cid:18) ξ , ξ , ξ , ξ , (cid:19) = (cid:18) . . φ . (cid:19) λ = (cid:18) λ λ (cid:19) = (cid:18) . . (cid:19) , (15)where for the present case, φ = 0 .
0. In this manner,jumps are added to the first dimension of the process,having an amplitude of ξ , = 0 . . t , given λ = 0 .
1. Similarly, jumps are addedto the second dimension, ξ , = 0 .
1, but the jumpsoccur three times more often than the aforementioned,given λ = 0 .
3. The influence of jumps can be observedacross all KM coefficients (see Fig. 3). The previouslysmooth KM surfaces become rugged from the fast vari-ations emerging due to the jumps, and the higher-orderKM coefficients—although small compared to the lower-order ones—clearly do not vanish. This indicates that thecontinuous stochastic modeling of time series of complexsystems (white-noise-driven Langevin equation) togetherwith the commonly used assumptions stemming from thePawula theorem are invalid for jump-diffusion processes.Modeling these processes with only the first two ordersof the KM expansion of the master equation is thereforeinsufficient.In order to understand further if it is possible to un-cover the jump amplitude terms of coupled processes, weuse the previous model Eq. (15) with φ = 0 .
3, therebyeffectively introducing a coupling via the off-diagonal el-ements of the jump matrix ξ . We show, in Fig 4, thecorresponding fourth-order KM coefficients. The impact -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.0030.0060.009 M [4 , and f [4 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.0030.0060.009 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.20.40.6 M [4 , and f [4 , y y × − FIG. 4. Two-dimensional Kramers–Moyal (KM) coefficients M [4 , , M [0 , , and M [4 , of bivariate jump-diffusion processes givenby Eq. (15) (with φ = 0 .
3) together with the respective theoretically expected functions f , associated with each coefficientaccording to Eqs. (10), (11), and (12). Notice that the estimated KM coefficients agree well with the theoretical expectedfunctions in all orders (error volumes estimated over the displayed domain). For further details, see Appendix B. of the coupling is visible, although small, in M [4 , , whichis no longer zero. Likewise, M [4 , and M [0 , also do notvanish. In Appendix B, we present the corresponding KMcoefficients up to order eight. B. Criteria for recovering coefficients in diffusionand jump-diffusion models
For the case of vanishing off-diagonal elements g , and ξ , , we can identify ways to recover the remaining coef-ficients of these processes.First, given that the noise d ω is Gaussian distributed, g is sign-reversal symmetric and one can thus assumethat it takes only positive values. One obtains that, if M [1 , = 0 then at least two elements of g must be zero,and if M [2 , = 0 then at least two elements of ξ mustbe zero (by assuming that λ and λ are non-vanishingrates). These findings reduce the dimensionality of theestimation procedure and ensure that the underlying pro-cesses are less complex that the full-fledged descriptionof Eq. (1), although they do not grant which coefficientsare zero-valued.Second, if one either employs a heuristic argument ofindependence of the jump processes or neglects the off-diagonal jump amplitudes ξ , and ξ , (e.g., by assumingthey are small compared to the diagonal terms of ξ ), onefinds the following approximations:15 M [6 , M [4 , = 15 153 s , λ s , λ = s , , M [0 , M [0 , = 15 153 s , λ s , λ = s , . (16)Likewise, the jump rates λ and λ can be obtained equiv-alently as1059 M [4 , M [8 , = 1059 3
105 ( s , λ ) s , λ = λ , M [0 , M [0 , = 1059 3
105 ( s , λ ) s , λ = λ . (17) Taking again model Eq. (15) with φ = 0 . s est1 , = 0 . ≈ . s , ,s est2 , = 0 . ≈ . s , . These estimated values (indicated by the superscript · est )are very close to the actual ones. The criteria and ap-proximations are especially relevant when constructingor analyzing systems which are known to have a specificunidirectional coupling form, e.g. a master-slave system,where for example the noise or the slave system are dic-tated by the driving master system. C. Factors influencing quality of recovery ofcoefficients
In order to validate the quality of the non-parametricrecovery of the KM coefficients, we now turn to two criti-cal aspects: firstly, bivariate processes may require a highnumber of data points in a time series for the estimationto be reliable; secondly, the interplay between the drift,diffusion, and jump parts of a stochastic processes mayrender the estimation incorrect.Addressing these aspects, we include a more contrivedmodel involving couplings and interactions in both thediffusion and jump terms, thus theoretically resulting inhaving all higher-order KM coefficients non-zero, and es-pecially the KM coefficients with (cid:96) (cid:54) = 0 and m (cid:54) = 0. Theparameters for the model read: N = (cid:18) N N (cid:19) = (cid:18) − x + x − x (cid:19) , g = (cid:18) g , g , g , g , (cid:19) = (cid:18) . . α . (cid:19) , ξ = (cid:18) ξ , ξ , ξ , ξ , (cid:19) = (cid:18) . . β . (cid:19) , λ = (cid:18) λ λ (cid:19) = (cid:18) . . (cid:19) . (18)The parameters α (diffusion term) and β (jump term)can be freely varied. . . . V e rr V [1 , V [0 , . . . V e rr V [2 , V [0 , . . . V e rr V [4 , V [0 , . . . V e rr V [6 , V [0 , n . . V e rr V [1 , V [2 , n . . V e rr V [4 , V [6 , FIG. 5. [Abscissa in logarithmic scale] Error volume V err for abivariate jump-diffusion process (Eq. (15)) depending on thenumber of data points n in the time series. Each process is nu-merically integrated with random initial conditions, for vary-ing number of data points n ∈ [10 , ] and over 50 times.The time sampling used was of 10 − . The average value of V err and one standard deviation (shaded area) are displayed. No-tice the clear decrease on all KM coefficients with either (cid:96) = 0or m = 0, e.g. M [2 , or M [0 , , as the number of data points n increases. This can be seen since the volume between the the-oretically expected values and the KM coefficients decreasesconsistently, i.e., V [ (cid:96),m ]err decreases for an increasing number ofdata points. The KM coefficients with (cid:96) (cid:54) = 0 and m (cid:54) = 0, suchas M [4 , or M [6 , present themselves as non-decreasing, butthe error volume is overall considerable small in value (cf.Fig. 6). It is important to notice that V [1 , does not convergeto zero since the KM coefficient is associated with the quarticpotential (i.e., the term N = − x + x ). Due to its shape, theprocess has two preferred states, either at x = − x = 1,thus spends little time at any intermediary point, like x = 0,damaging the statistics of the recovery. Let us focus firstly on the number of data points in atime series. We utilize the models Eq. (15) with φ = 0 . α = β = 0 .
3, and show in Figs. 5and 6, respectively, the error volumes V err for the KMcoefficients for an increasing number of data points. Thereliability of the recovery of KM coefficient is valid forhigher amount of data ( n ≥ ), as expected, althoughthe scarcity of data posits no extensive problem for thecalculation. It is especially important to notice that atime series with a lower amount of data entails naturallyfewer jumps in the process, hindering the possibility ofaccurately recovering the jump terms from such shorttime series. For n ≥ , the estimation seems reliable,the standard deviations become minute, and most errorvalues approach zero, i.e., the theoretical and estimatedKM surfaces are close.One remark is necessary on the recovery of the driftterms. The presence of noise and jumps in the processtakes its toll on the recovery of the exact form of the KMcoefficients as well as the explicit dependence of the statevariables, i.e., the quartic potential in both Eqs. (15) V e rr V [1 , V [0 , . . . V e rr V [2 , V [0 , V e rr V [4 , V [0 , V e rr V [6 , V [0 , n . . . V e rr V [1 , V [2 , n . . V e rr V [4 , V [6 , FIG. 6. [Abscissa in logarithmic scale] Same as Fig. 5 but forthe bivariate jump-diffusion process (Eq. (18)) with α = β =0 .
3. Integration parameters as in Fig. 5. and (18). A finer time sampling can help to improve theresults.To further test the limitation of retrieving the KMcoefficients from data, we utilize model Eq. (18) oncemore and investigate the influence of parameters α (dif-fusion term) and β (jump term). For increasing α ( α ∈ [10 − , ]) and β = 0 .
3, we observe a considerable im-pact on the error volume V err after the order of magnitudeon the diffusion parameter α is ten-fold bigger in com-parison to the diffusion parameter g , (Fig. 7). Similarly,for increasing β ( β ∈ [10 − , ]) and α = 0 .
3, the errorvolume V err is considerably impacted already when β isof similar size as the other parameters, namely ξ , = 0 . M [1 , , M [0 , , . . . , M [0 , as presented in Fig. 9takes about 5.5 s on a desktop computer (quad-core withclock speed 2.20 GHz) for a bi-dimensional time seriesof size 2 × data points (and 22 s for 2 × datapoints). Our approach might thus be advantageous forfield applications that aim at an investigation of interac-tions between complex systems with poorly understooddynamics. V. CONCLUSION
We introduced the bivariate jump-diffusion processwhich comprises two-dimensional diffusion and two-dimensional jumps that can be coupled to one an- V e rr V [1 , V [0 , V e rr V [2 , V [0 , V e rr V [4 , V [0 , V e rr V [6 , V [0 , − − α − V e rr V [1 , V [2 , − − α V e rr V [4 , V [6 , FIG. 7. [Double logarithmic scale] Error volume V err for thebivariate jump-diffusion process Eq. (18) for a varying diffu-sion term α ∈ [10 − , ]. The vertical dotted line at α = 0 . α = g , = 0 . g , = 0 .
5. A small value of the diffusion term α ,in comparison to the diffusion parameters g , and g , , en-sures a good reconstruction, i.e., a small error volume V err .The average and one-standard deviations (shaded area) aredisplayed. For each point 50 iterations are taken, each with atotal number of data points of 5 × and a time samplingof 10 − . − V e rr V [1 , V [0 , − V e rr V [2 , V [0 , V e rr V [4 , V [0 , V e rr V [6 , V [0 , − − β − − V e rr V [1 , V [2 , − − β V e rr V [4 , V [6 , FIG. 8. [Double logarithmic scale] Error volume V err for thebivariate jump-diffusion process Eq. (18) for a varying jumpterm β ∈ [10 − , ]. The vertical dotted line indicates thebiggest jump term value ξ , = 0 . β . Asin direct analogy to Fig. 7, a small jump term β ensures agood reconstruction, i.e., a small error volume V err . Increas-ing values of the jump term β in comparison to the otherparameters in the system make the reconstruction unreliable.Iteration scheme identical to the one in Fig. 7. other. For such a process we presented a data-driven,non-parametric estimation procedure of higher-orderKramers–Moyal coefficients and investigated its pros andcons using synthetic bivariate time series from continuousand discontinuous processes. The procedure allows oneto reconstruct relevant aspects of the underlying jump-diffusion processes and to recover the underlying param-eters.Having now a traceable mathematical framework,the model can be extended to embody other noiseand jump properties. An extension from the underly-ing Wiener process to include e.g. fractional Brownianmotion is straightforward but nevertheless requires fur-ther investigations to derive an explicit forward Kol-mogorov equation [19]. Also, a generalization to con-tinuous jump processes—originating from alpha-stableor other heavy-tailed distributions (L´evy noise-drivenLangevin equation)—is possible, however, with the draw-back that calculating the conditional moments may notalways be mathematically possible [19]. On the otherhand, a numerical estimation of generalised momentsshould be possible but these still require a physical in-terpretation.We are confident that our novel approach provides ageneral avenue to further understanding of interactingcomplex systems (e.g. brain or power grids [33, 41–43])whose dynamics exhibit nontrivial noise contributions. ACKNOWLEDGEMENTS
The authors would like to thank Thorsten Rings andMehrnaz Anvari for interesting discussions, and Fran-cisco Meirinhos for the tremendous help in devising themethodology behind the results. L.R.G. thanks Dirk Wit-thaut and Giulia di Nunno for everlasting support, andgratefully acknowledges support by the grant No. VH-NG-1025 and a scholarship from the E.ON Stipendien-fonds. [1] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, andD.-U. Hwang, “Complex networks: Structure and dynam-ics,” Phys. Rep. , 175–308 (2006).[2] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Phys.Rep. , 93–153 (2008).[3] E. T. Bullmore and D. S. Bassett, “Brain graphs: Graph-ical models of the human brain connectome,” Annu. Rev. Clin. Psychol. , 113–140 (2011).[4] M. Barth´elemy, “Spatial networks,” Phys. Rep. , 1–101 (2011).[5] M. E. J. Newman, “Communities, modules and large-scale structure in networks,” Nat. Phys. , 25–31 (2012).[6] P. Holme and J. Saram¨aki, “Temporal networks,” Phys.Rep. , 97–125 (2012).[7] M. Kivel¨a, A. Arenas, M. Barthelemy, J. P. Gleeson,Y. Moreno, and M. A. Porter, “Multilayer networks,”J. Complex Netw. , 203–271 (2014).[8] A. S. Pikovsky, M. G. Rosenblum, and J. Kurths, Syn-chronization: A universal concept in nonlinear sciences (Cambridge University Press, Cambridge, UK, 2001).[9] H. Kantz and T. Schreiber,
Nonlinear Time Series Anal-ysis , 2nd ed. (Cambridge University Press, Cambridge,UK, 2003).[10] E. Pereda, R. Quian Quiroga, and J. Bhattacharya,“Nonlinear multivariate analysis of neurophysiologicalsignals,” Prog. Neurobiol. , 1–37 (2005).[11] K. Hlav´aˇckov´a-Schindler, M. Paluˇs, M. Vejmelka,and J. Bhattacharya, “Causality detection based oninformation-theoretic approaches in time series analysis,”Phys. Rep. , 1–46 (2007).[12] N. Marwan, M. C. Romano, M. Thiel, and J. Kurths,“Recurrence plots for the analysis of complex systems,”Phys. Rep. , 237–329 (2007).[13] K. Lehnertz, S. Bialonski, M.-T. Horstmann, D. Krug,A. Rothkegel, M. Staniek, and T. Wagner, “Synchro-nization phenomena in human epileptic brain networks,”J. Neurosci. Methods , 42–48 (2009).[14] K. Lehnertz, “Assessing directed interactions from neu-rophysiological signals – an overview,” Physiol. Meas. ,1715–1724 (2011).[15] T. Stankovski, T. Pereira, P. V. E. McClintock, andA. Stefanovska, “Coupling functions: Universal insightsinto dynamical interaction mechanisms,” Rev. Mod.Phys. , 045001 (2017).[16] H. Risken, The Fokker-Planck Equation , 2nd ed.,Springer Series in Synergetics (Springer, Berlin, 1996).[17] N. G. Van Kampen,
Stochastic processes in physics andchemistry (North Holland, Amsterdam, 1981).[18] R. Friedrich, J. Peinke, M. Sahimi, and M. R. R. Tabar,“Approaching complexity by stochastic methods: Frombiological systems to turbulence,” Phys. Rep. , 87–162 (2011).[19] M. R. R. Tabar,
Analysis and Data-Based Reconstruc-tion of Complex Nonlinear Dynamical Systems: Using theMethods of Stochastic Processes (Springer-Nature, 2019).[20] J. Prusseit and K. Lehnertz, “Measuring interdepen-dences in dissipative dynamical systems with estimatedFokker-Planck coefficients,” Phys. Rev. E , 041914(2008).[21] B. Lehle, “Stochastic time series with strong, correlatedmeasurement noise: Markov analysis in n dimensions,” J.Stat. Phys. , 1145–1169 (2013).[22] T. Scholz, F. Raischel, V. V. Lopes, B. Lehle,M. W¨achter, J. Peinke, and P. G. Lind, “Parameter-free resolution of the superposition of stochastic signals,”Phys. Lett. A , 194–206 (2017).[23] B. Wahl, U. Feudel, J. Hlinka, M. W¨achter, J. Peinke,and J. A. Freund, “Granger-causality maps of diffusionprocesses,” Phys. Rev. E , 022213 (2016).[24] P. G. Lind, M. Haase, F. B¨ottcher, J. Peinke, D. Klein-hans, and R. Friedrich, “Extracting strong measurement noise from stochastic time series: Applications to empir-ical data,” Phys. Rev. E , 041125 (2010).[25] M. B. Weissman, “ f noise and other slow, nonexponen-tial kinetics in condensed matter,” Rev. Mod. Phys. ,537–571 (1988).[26] G. Bakshi, C. Cao, and Z. Chen, “Empirical perfor-mance of alternative option pricing models,” J. Finance , 2003–2049 (1997).[27] D. Duffie, J. Pan, and K. Singleton, “Transform analysisand asset pricing for affine jump-diffusions,” Economet-rica , 1343–1376 (2000).[28] T. G. Andersen, L. Benzoni, and J. Lund, “An empiricalinvestigation of continuous-time equity return models,”J. Finance , 1239–1284 (2002).[29] S. R. Das, “The surprise element: jumps in interestrates,” J. Econometrics , 27–65 (2002).[30] M. Johannes, “The statistical and economic role of jumpsin continuous-time interest rate models,” J. Finance ,227–260 (2004).[31] Z. Cai and Y. Hong, “Some recent developments innonparametric finance,” in Nonparametric EconometricMethods , Advances in Econometrics, Volume 25, editedby Q. Li and J. S. Racine (Emerald Group PublishingLimited, 2009) pp. 379–432.[32] M. Anvari, M. R. R. Tabar, J. Peinke, and K. Lehnertz,“Disentangling the stochastic behavior of complex timeseries,” Sci. Rep. , 35435 (2016).[33] K. Lehnertz, L. Zabawa, and M. R. R. Tabar, “Charac-terizing abrupt transitions in stochastic dynamics,” NewJ. Physics , 113043 (2018).[34] C. T. Chudley and R. J. Elliott, “Neutron scattering froma liquid on a jump diffusion model,” Proc. Phys. Soc. ,353 (1961).[35] R. C. Merton, “Option pricing when underlying stockreturns are discontinuous,” J. Financ. Econ. , 125–144(1976).[36] P. L. Hall and D. K. Ross, “Incoherent neutron scatter-ing functions for random jump diffusion in bounded andinfinite media,” Mol. Phys. , 673–682 (1981).[37] M. T. Giraudo and L. Sacerdote, “Jump-diffusion pro-cesses as models for neuronal activity,” Biosystems ,75–82 (1997).[38] R. F. Pawula, “Approximation of the linear Boltzmannequation by the Fokker-Planck equation,” Phys. Rev. , 186–188 (1967).[39] J. Prusseit and K. Lehnertz, “Stochastic qualifiers ofepileptic brain dynamics,” Phys. Rev. Lett. , 138103(2007).[40] L. R. Gorj˜ao and F. Meirinhos, “Python KM: Kramers–Moyal coefficients for stochastic processes,” J. OpenSource Softw. (under review). Preview available at https://github.com/openjournals/joss-papers/blob/joss.01679/joss.01679/10.21105.joss.01679.pdf .[41] B. Sch¨afer, C. Beck, K. Aihara, D. Witthaut, andM. Timme, “Non-gaussian power grid frequency fluctua-tions characterized by L´evy-stable laws and superstatis-tics,” Nat. Energy , 119 (2018).[42] L. R. Gorj˜ao, M. Anvari, H. Kantz, D. Witthaut,M. Timme, and B. Sch¨afer, “Data-driven model ofthe power-grid frequency dynamics,” arXiv preprintarXiv:1909.08346 (2019).[43] M. Anvari, L. R. Gorj˜ao, M. Timme, B. Sch¨afer, D. Wit-thaut, and H. Kantz, “Stochastic properties of the fre- quency dynamics in real and synthetic power grids,”arXiv preprint arXiv:1909.09110 (2019). Appendices
A. EXTENDED DERIVATION OF THE TWO-DIMENSIONAL KRAMERS–MOYAL COEFFICIENTSFOR A JUMP-DIFFUSION PROCESS
The following derivations stem from Eq. (3) and apply to the two-dimensional jump-diffusion process ( y , y ), as inEq. (1). All orders of the Kramers–Moyal coefficients M [ (cid:96),m ] are ( (cid:96), m ) ∈ N + .
1. Kramers–Moyal coefficients M [1 , and M [0 , M [1 , ( x , x ) = lim d t → t (cid:104) (d y ) (d x ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) d y (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) N d t + g , d w + g , d w + ξ , d J + ξ , d J (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t [ N d t + g , (cid:104) d w (cid:105) + g , (cid:104) d w (cid:105) + (cid:104) ξ , (cid:105)(cid:104) d J (cid:105) + (cid:104) ξ , (cid:105)(cid:104) d J (cid:105) ]= N , where (cid:104) g i,j d W j (cid:105) = (cid:104) g i,j (cid:105)(cid:104) d W j (cid:105) = 0, because a Wiener process has the property (cid:104) d W j (cid:105) = 0. Further (cid:104) ξ i,j d J j (cid:105) = (cid:104) ξ i,j (cid:105)(cid:104) d J j (cid:105) = 0, since ξ i,j is a Gaussian with zero mean, i.e., (cid:104) ξ i,j (cid:105) = 0.Mutatis mutandis for M [0 , .
2. Kramers–Moyal coefficient M [1 , M [1 , = lim d t → t (cid:104) (d y ) (d y ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) · ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → (cid:20) N N d t + g , g , (cid:104) (d w ) (cid:105) t + g , g , (cid:104) (d w ) (cid:105) dt + O (d t ) (cid:21) = g , g , + g , g , , where higher-order terms O (d t ) (cid:15) , with (cid:15) >
0, vanish in the limit d t →
0. Recall as well (cid:104) (d w i ) (cid:105) = d t .
3. Kramers–Moyal coefficients M [2 , and M [0 , M [2 , = lim d t → t (cid:104) (d y ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → (cid:20) N d t + g , (cid:104) (d w ) (cid:105) t + g , (cid:104) (d w ) (cid:105) t + (cid:104) ξ , (cid:105)(cid:104) (d J ) (cid:105) t + (cid:104) ξ , (cid:105)(cid:104) (d J ) (cid:105) t + O (d t ) (cid:21) = (cid:2) g , + s , λ + g , + s , λ (cid:3) , (cid:104) ξ ij (cid:105) = σ ξ ij = s ij as well as (cid:104) (d J i ) (cid:105) = λ i d t .Mutatis mutandis for M [0 , M [0 , = (cid:2) g , + s , λ + g , + s , λ (cid:3) ,
4. Kramers–Moyal coefficient M [2 , M [2 , = lim d t → t (cid:104) (d y ) (d y ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) · ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:20) terms( N , N , O (d t )) + terms( g ij , O (d t )) + terms(mixing ξ ij )+ (cid:104) ξ , (cid:105)(cid:104) ξ , (cid:105)(cid:104) (d J ) (cid:105) + (cid:104) ξ , (cid:105)(cid:104) ξ , (cid:105)(cid:104) (d J ) (cid:105) + (cid:104) ξ , (cid:105)(cid:104) ξ , (cid:105)(cid:104) (d J ) (cid:105)(cid:104) (d J ) (cid:105) + (cid:104) ξ , (cid:105)(cid:104) ξ , (cid:105)(cid:104) (d J ) (cid:105)(cid:104) (d J ) (cid:105) (cid:21) = [ s , s , λ + s , s , λ ] . Terms including d t on the right-hand side of the above equation vanish for d t →
0, where as well (cid:104) ξ , ξ , (cid:105) = (cid:104) ξ , (cid:105)(cid:104) ξ , (cid:105) = 0, and t [ (cid:104) (d J ) (cid:105)(cid:104) (d J ) (cid:105) ] = t [ λ d tλ d t ] ∝ d t vanishes in the limit d t →
5. Kramers–Moyal coefficients M [ (cid:96),m ] , for × ( (cid:96), m ) ≥ For (2 (cid:96), m ), with ( (cid:96), m ) ≥
4, the Kramers–Moyal coefficients M [2 (cid:96), m ] are as follows M [2 (cid:96), m ] = lim d t → t (cid:104) (d y ) (cid:96) (d y ) m (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) (cid:96) · ( N d t + g , d w + g , d w + ξ , d J + ξ , d J ) m (cid:105)| y ( t )= x ,y ( t )= x = lim d t → t (cid:104) (cid:104) ξ (cid:96) , (cid:105)(cid:104) ξ m , (cid:105)(cid:104) (d J ) (cid:96) + m ) (cid:105) + (cid:104) ξ (cid:96) , (cid:105)(cid:104) ξ m , (cid:105)(cid:104) (d J ) (cid:96) + m ) (cid:105) (cid:105) = (cid:2) (cid:104) ξ (cid:96) , (cid:105)(cid:104) ξ m , (cid:105) λ + (cid:104) ξ (cid:96) , (cid:105)(cid:104) ξ m , (cid:105) λ (cid:3) = (2 (cid:96) )!2 (cid:96) (cid:96) ! (2 m )!2 m m ! (cid:2) s (cid:96) , s m , λ + s (cid:96) , s m , λ (cid:3) . In the last step, take the fact that the jump amplitudes ξ i,j are Gaussian distributed, thus (cid:104) ξ (cid:96)i,j (cid:105) ∝ σ (cid:96)ξ i,j = s (cid:96)i,j . In thismanner, all Kramers–Moyal coefficients M [2 (cid:96), m ] , with ( (cid:96), m ) ≥ B. EXTENDED RESULTS FOR MODELED DATA BY EQ. (15) -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.3-1.0-0.50.00.51.0 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.3-0.4-0.20.00.20.4 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.050.10.150.2 M [1 , and f [1 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.10.20.30.4 M [2 , and f [2 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.3-0.10.00.10.2 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.3-0.0050.00.0050.010.015 M [2 , and f [2 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.0030.0060.009 M [4 , and f [4 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.0030.0060.009 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.20.40.6 M [4 , and f [4 , y y × − -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.3-0.00040.00.00040.00080.0012 M [6 , and f [6 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.3-0.00040.00.00040.00080.0012 M [0 , and f [0 , y y -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.00.20.40.60.8 M [6 , and f [6 , y y × − -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.30.51.01.52.02.53.0 M [8 , and f [8 , y y × − -0.8-0.40.00.40.8-0.3-0.2-0.10.00.10.20.32468 M [0 , and f [0 , y y × − FIG. 9. Two-dimensional Kramers–Moyal (KM) coefficients M [ (cid:96),m ] of bivariate jump-diffusion processes given by Eq. (15) andall theoretical expected functions f [ (cid:96),m ] associated with each KM coefficient according to Eqs. (10), (11), and (12). Shown arethe KM coefficients M [1 , , M [0 , , M [1 , , M [2 , , M [0 , , M [2 , , M [4 , , M [0 , , M [4 , , M [6 , , M [0 , , M [6 , , M [8 , , and M [0 , . The respective error volumes read: V [1 , = 0 . V [0 , = 0 . V [1 , = 0 . V [2 , = 0 . V [0 , = 0 . V [2 , = 0 . V [4 , = 0 . V [0 , = 0 . V [4 , < . V [6 , = 0 . V [0 , = 0 . V [6 , = 0 . V [8 , = 0 . V [0 , = 0 ..