Analysis and data-driven reconstruction of bivariate jump-diffusion processes

TL;DR

This paper introduces a data-driven, non-parametric method to reconstruct and analyze bivariate jump-diffusion processes, enabling parameter estimation and process understanding from synthetic data.

Contribution

It presents a novel approach for reconstructing bivariate jump-diffusion processes using higher-order Kramers-Moyal coefficients, validated with synthetic data.

Findings

Successful reconstruction of jump-diffusion processes from synthetic data

Effective parameter estimation using higher-order Kramers-Moyal coefficients

Identification of limitations due to data scarcity and parameter imbalance

Abstract

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-parametric estimation procedure of higher-order (up to 8) Kramers-Moyal coefficients that allows one to reconstruct relevant aspects of the underlying jump-diffusion processes and to recover the underlying parameters. The procedure is validated with numerically integrated data using synthetic bivariate time series from continuous and discontinuous processes. We further evaluate the possibility of estimating the parameters of the jump-diffusion model via data-driven analyses of the higher-order Kramers-Moyal coefficients, and the limitations arising from the scarcity of points in the data or disproportionate parameters in the system.

Equations

$$\stackrel{\mathbf{y}}{\overbrace{\left(\begin{array}{c}\mathrm{d}{y}_1(t)\\ \mathrm{d}{y}_2(t)\end{array}\right)}}$$

$$\underset{\text{Poissonian jumps}}{\underbrace{\stackrel{\mathbf{\xi }}{\overbrace{\left(\begin{array}{cc}{\xi }_{1,1}& {\xi }_{1,2}\\ {\xi }_{2,1}& {\xi }_{2,2}\end{array}\right)}}\stackrel{\mathrm{d}\mathbf{J}}{\overbrace{\left(\begin{array}{c}\mathrm{d}{J}_1\\ \mathrm{d}{J}_2\end{array}\right)}}}},$$

$$\begin{array}{r}\stackrel{\mathbf{y}}{\overbrace{\left(\begin{array}{c}\mathrm{d}{y}_1(t)\\ \mathrm{d}{y}_2(t)\end{array}\right)}}=\underset{\text{drift}}{\underbrace{\stackrel{\mathbf{N}}{\overbrace{\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)}}\mathrm{d}t}}+\underset{\text{diffusion}}{\underbrace{\stackrel{\mathbf{g}}{\overbrace{\left(\begin{array}{cc}{g}_{1,1}& g_{1,2}\\ g_{2,1}& g_{2,2}\end{array}\right)}}\stackrel{\mathrm{d}\mathbf{w}}{\overbrace{\left(\begin{array}{c}\mathrm{d}{w}_1\\ \mathrm{d}{w}_2\end{array}\right)}}}}\end{array}.$$

$${\mathcal{M}}^{[\ell ,m]}(x_1,x_2)=$$

$$\underset{\Delta t\to 0}{\lim }\frac{1}{\Delta t}\int (y_1(t+\Delta t)-y_1(t){)}^{\ell }(y_2(t+\Delta t)-y_2(t){)}^m\cdot$$

$$P(y_1,y_2;t+\Delta t|y_1,y_2;t){|}_{y_1(t)=x_1,y_2(t)=x_2}\mathrm{d}{y}_1\mathrm{d}{y}_2,$$

$${\mathcal{M}}^{[\ell ,m]}(x_1,x_2)$$

$$\Delta y_i$$

$$\begin{array}{rl}{\mathcal{M}}^{[1,0]}& =N_1,\\ {\mathcal{M}}^{[0,1]}& =N_2,\end{array}$$

$$\begin{array}{rl}{\mathcal{M}}^{[1,1]}& =g_{1,1}{g}_{2,1}+g_{1,2}{g}_{2,2},\\ {\mathcal{M}}^{[2,0]}& =\left[g_{1,1}^2+g_{1,2}^2\right],\\ {\mathcal{M}}^{[0,2]}& =\left[g_{2,1}^2+g_{2,2}^2\right].\end{array}$$

$$\int {\int }_U{\left({\mathcal{M}}^{[l,m]}(y_1,y_2)-f^{[l,m]}(y_1,y_2)\right)}^2\mathrm{d}{y}_1\mathrm{d}{y}_2=:V^2,$$

$$\int {\int }_U{\left({\mathcal{M}}^{[l,m]}-f^{[l,m]}\right)}^{2}p(y_1,y_2)\mathrm{d}{y}_1\mathrm{d}{y}_2=V_{\text{err}}^2,$$

$$\mathbf{N}=\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)=\left(\begin{array}{c}-x_1^3+x_1\\ -x_2\end{array}\right),$$

$$\mathbf{g}=\left(\begin{array}{cc}{g}_{1,1}& g_{1,2}\\ g_{2,1}& g_{2,2}\end{array}\right)=\left(\begin{array}{cc}0.5& 0.0\\ 0.0& 0.5\end{array}\right),$$

$$\mathbf{N}=\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)=\left(\begin{array}{c}-x_1^3+x_1\\ -x_2\end{array}\right),$$

$$\mathbf{g}=\left(\begin{array}{cc}{g}_{1,1}& g_{1,2}\\ g_{2,1}& g_{2,2}\end{array}\right)=\left(\begin{array}{cc}0.1+x_1^2& 0.5\\ 0.0& 0.2+2x_2^2\end{array}\right).$$

$$\begin{array}{rl}{\mathcal{M}}^{[1,0]}& =N_1\\ {\mathcal{M}}^{[0,1]}& =N_2\end{array}$$

$$\begin{array}{rl}{\mathcal{M}}^{[1,1]}& =g_{1,1}{g}_{2,1}+g_{1,2}{g}_{2,2}\\ {\mathcal{M}}^{[2,0]}& =\left[g_{1,1}^2+s_{1,1}{\lambda }_1+g_{1,2}^2+s_{1,2}{\lambda }_2\right]\\ {\mathcal{M}}^{[0,2]}& =\left[g_{2,1}^2+s_{2,1}{\lambda }_1+g_{2,2}^2+s_{2,2}{\lambda }_2\right]\end{array}$$

$$\begin{array}{rl}{\mathcal{M}}^{[2,2]}& =\left[s_{1,1}{s}_{2,1}{\lambda }_1+s_{1,2}{s}_{2,2}{\lambda }_2\right]\\ {\mathcal{M}}^{[4,0]}& =3\left[s_{1,1}^2{\lambda }_1+s_{1,2}^2{\lambda }_2\right]\\ {\mathcal{M}}^{[0,4]}& =3\left[s_{2,1}^2{\lambda }_1+s_{2,2}^2{\lambda }_2\right]\\ {\mathcal{M}}^{[4,4]}& =9\left[s_{1,1}^2{s}_{2,1}^2{\lambda }_1+s_{1,2}^2{s}_{2,2}^2{\lambda }_2\right]\end{array}$$

$$\begin{array}{rl}{\mathcal{M}}^{[6,0]}& =15\left[s_{1,1}^3{\lambda }_1+s_{1,2}^3{\lambda }_2\right]\\ {\mathcal{M}}^{[0,6]}& =15\left[s_{2,1}^3{\lambda }_1+s_{2,2}^3{\lambda }_2\right]\\ {\mathcal{M}}^{[6,6]}& =225\left[s_{1,1}^3{s}_{2,1}^3{\lambda }_1+s_{1,2}^3{s}_{2,2}^3{\lambda }_2\right]\end{array}$$

$$\begin{array}{rl}{\mathcal{M}}^{[8,0]}& =105\left[s_{1,1}^4{\lambda }_1+s_{1,2}^4{\lambda }_2\right]\\ {\mathcal{M}}^{[0,8]}& =105\left[s_{2,1}^4{\lambda }_1+s_{2,2}^4{\lambda }_2\right]\\ {\mathcal{M}}^{[8,8]}& =11025\left[s_{2,1}^4{\lambda }_1+s_{2,2}^4{\lambda }_2\right]\end{array}$$

$$M^{[2\ell ,2m]}=\frac{(2\ell )!}{2^{\ell }\ell !}\frac{(2m)!}{2^{m}m!}\left[s_{1,1}^{\ell }{s}_{2,1}^m{\lambda }_1+s_{1,2}^{\ell }{s}_{2,2}^m{\lambda }_2\right].$$

$$\mathbf{N}=\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)=\left(\begin{array}{c}-x_1^3+x_1\\ -x_2\end{array}\right)$$

$$\mathbf{g}=\left(\begin{array}{cc}{g}_{1,1}& g_{1,2}\\ g_{2,1}& g_{2,2}\end{array}\right)=\left(\begin{array}{cc}0.1& 0.5\\ 0.0& 0.2\end{array}\right)$$

$$\mathbf{\xi }=\left(\begin{array}{cc}{\xi }_{1,1}& {\xi }_{1,2}\\ {\xi }_{2,1}& {\xi }_{2,2}\end{array}\right)=\left(\begin{array}{cc}0.2& 0.0\\ \varphi & 0.1\end{array}\right)$$

$$\mathbf{\lambda }=\left(\begin{array}{c}{\lambda }_1\\ {\lambda }_2\end{array}\right)=\left(\begin{array}{c}0.1\\ 0.3\end{array}\right),$$

$$\frac{1}{5}\frac{{\mathcal{M}}^{[6,0]}}{{\mathcal{M}}^{[4,0]}}=\frac{1}{5}\frac{15}{3}\frac{s_{1,1}^3{\lambda }_1}{s_{1,1}^2{\lambda }_1}$$

$$\frac{1}{5}\frac{{\mathcal{M}}^{[0,6]}}{{\mathcal{M}}^{[0,4]}}=\frac{1}{5}\frac{15}{3}\frac{s_{2,2}^3{\lambda }_1}{s_{2,2}^2{\lambda }_1}$$

$$\frac{105}{9}\frac{{{\mathcal{M}}^{[4,0]}}^2}{{\mathcal{M}}^{[8,0]}}=\frac{105}{9}\frac{3^2}{105}\frac{(s_{1,1}^2{\lambda }_1{)}^2}{s_{1,1}^4{\lambda }_1}$$

$$\frac{105}{9}\frac{{{\mathcal{M}}^{[0,4]}}^2}{{\mathcal{M}}^{[0,8]}}=\frac{105}{9}\frac{3^2}{105}\frac{(s_{2,2}^2{\lambda }_2{)}^2}{s_{2,2}^4{\lambda }_2}$$

$$s_{1,1}^{\mathrm{est}}=0.16\approx 0.2=s_{1,1},$$

$$s_{2,2}^{\mathrm{est}}=0.09\approx 0.1=s_{2,2}.$$

$$\mathbf{N}=\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)=\left(\begin{array}{c}-x_1^3+x_1\\ -x_2\end{array}\right),$$

$$\mathbf{g}=\left(\begin{array}{cc}{g}_{1,1}& g_{1,2}\\ g_{2,1}& g_{2,2}\end{array}\right)=\left(\begin{array}{cc}0.1& 0.5\\ \alpha & 0.2\end{array}\right),$$

$$\mathbf{\xi }=\left(\begin{array}{cc}{\xi }_{1,1}& {\xi }_{1,2}\\ {\xi }_{2,1}& {\xi }_{2,2}\end{array}\right)=\left(\begin{array}{cc}0.2& 0.5\\ \beta & 0.1\end{array}\right),$$

$$\mathbf{\lambda }=\left(\begin{array}{c}{\lambda }_1\\ {\lambda }_2\end{array}\right)=\left(\begin{array}{c}0.1\\ 0.3\end{array}\right).$$

$${\mathcal{M}}^{[1,0]}(x_1,x_2)$$

$$=\underset{\mathrm{d}t\to 0}{\lim }\frac{1}{\mathrm{d}t}⟨\mathrm{d}{y}_1⟩{|}_{y_1(t)=x_1,y_2(t)=x_2}$$

$$=\underset{\mathrm{d}t\to 0}{\lim }\frac{1}{\mathrm{d}t}⟨N_1\mathrm{d}t+g_{1,1}\mathrm{d}{w}_1+g_{1,2}\mathrm{d}{w}_2+{\xi }_{1,1}\mathrm{d}{J}_1+{\xi }_{1,2}\mathrm{d}{J}_2⟩{|}_{y_1(t)=x_1,y_2(t)=x_2}$$

$$=\underset{\mathrm{d}t\to 0}{\lim }\frac{1}{\mathrm{d}t}\left[N_1\mathrm{d}t+g_{1,1}⟨\mathrm{d}{w}_1⟩+g_{1,2}⟨\mathrm{d}{w}_2⟩+⟨{\xi }_{1,1}⟩⟨\mathrm{d}{J}_1⟩+⟨{\xi }_{1,2}⟩⟨\mathrm{d}{J}_2⟩\right]$$

$$=N_1,$$

$${\mathcal{M}}^{[1,1]}$$

Full text

I Introduction

Research over the last two decades has demonstrated the high suitability of the network paradigm in advancing our understanding of natural and man-made complex dynamical systems Boccaletti et al. (2006); Arenas et al. (2008); Bullmore and Bassett (2011); Barthélemy (2011); Newman (2012); Holme and Saramäki (2012); Kivelä et al. (2014). With this paradigm, a system component is represented by a vertex and interactions between components are conveyed by edges connecting vertices, and graph theory provides a large repertoire of methods to characterize networks on various scales.

Characterizing properties of interactions using the knowledge of the dynamics of each of the components is key to understanding real-world systems. To achieve this goal, a large number of time-series-analysis methods has been developed that originate from synchronization theory, nonlinear dynamics, information theory, and from statistical physics (for an overview, see Refs. Pikovsky et al. (2001); Kantz and Schreiber (2003); Pereda et al. (2005); Hlaváčková-Schindler et al. (2007); Marwan et al. (2007); Lehnertz et al. (2009); Lehnertz (2011); Stankovski et al. (2017)). Some of these methods make rather strict assumptions about the dynamics of network components generating the time series and many approaches preferentially focus on the low-dimensional deterministic part of the dynamics.

Real-world systems, however, are typically influenced by random forcing and interactions between constituents are highly non-linear, which results in very complex, stochastic, and non-stationary system behavior that exhibits both deterministic and stochastic features. Aiming at determining characteristics and strength of fluctuating forces as well as at assessing properties of non-linear interactions, the analysis of such systems is associated with the problem of retrieving a stochastic dynamical system from measured time series. There is a substantial existing literature Risken (1996); Van Kampen (1981); Friedrich et al. (2011); Tabar (2019) for the modeling of complex dynamical systems which employs the conventional Langevin equation that is based on the first- and second-order Kramers–Moyal (KM) coefficients, known as drift and diffusion terms. All functions and parameters of this modeling can be found directly from the measured time series employing a widely used non-parametric approach. There are by now only few studies that make use of this ansatz to characterize interactions between stochastic processes Prusseit and Lehnertz (2008); Lehle (2013); Scholz et al. (2017); Wahl et al. (2016); Lind et al. (2010).

Despite its successful application in diverse scientific fields, growing evidence indicates that the continuous stochastic modeling of time series of complex systems (white-noise-driven Langevin equation) should account for the presence of discontinuous jump components Weissman (1988); Bakshi et al. (1997); Duffie et al. (2000); Andersen et al. (2002); Das (2002); Johannes (2004); Cai and Hong (2009); Anvari et al. (2016); Lehnertz et al. (2018); Tabar (2019). In this context, the jump-diffusion model Chudley and Elliott (1961); Merton (1976); Hall and Ross (1981); Giraudo and Sacerdote (1997) was shown to provide a theoretical tool to study processes of known and unknown nature that exhibit jumps. It allows one to separate the deterministic drift term as well as different stochastic behaviors, namely diffusive and jumpy behavior Anvari et al. (2016); Lehnertz et al. (2018); Tabar (2019). Moreover, all of the unknown functions and coefficients of a dynamical stochastic equation that describe a jump-diffusion process can be derived directly from measured time series. This approach involves estimating higher-order ($\geq 3$) KM coefficients and it provides an intuitive physical meaning of these coefficients.

The focus of this paper is to introduce a method to investigate bivariate time-series with discontinuous jump components. We begin with an overview of bivariate diffusion processes that exhibit the known relation between the parameters and functions in stochastic modeling and the KM coefficients. Exemplary processes are portrayed, and we propose a measure to judge the quality of our reconstruction procedure. We then present bivariate jump-diffusion processes alongside with the associated KM expansion Anvari et al. (2016), and investigate the suitability of our reconstruction procedure using various examples. We conclude this paper by summarizing our findings.

II State of the art: Bivariate jump-diffusion model

A bivariate jump-diffusion process can be modeled via Anvari et al. (2016); Tabar (2019)

[TABLE]

where all the elements of vectors $\boldsymbol{N}$, $\boldsymbol{\mathrm{d}J}$, and $\boldsymbol{\mathrm{d}w}$ as well as of matrices $\boldsymbol{g}$ and $\boldsymbol{\xi}$ may, in general, be state- and time-dependent (dependencies not shown for convenience of notation). The drift coefficient is a two-dimensional vector $\boldsymbol{N}=(N_{1},N_{2})$ with $\boldsymbol{N}\in\mathbb{R}^{2}$, where each dimension of $\boldsymbol{N}$, i.e., $N_{i}$, may depend on $y_{1}(t)$ and $y_{2}(t)$. The diffusion coefficient takes a matrix $\boldsymbol{g}\in\mathbb{R}^{2\times 2}$. The two Wiener processes $\boldsymbol{w}=(w_{1},w_{2})$ act as independent Brownian noises for the state variables $y_{1}(t)$ and $y_{2}(t)$. The diagonal elements of $\boldsymbol{g}$ comprise the diffusion coefficients of self-contained stochastic diffusive processes, and the off-diagonal elements represent interdependencies between the two Wiener processes, i.e., they result from an interaction between the two processes. Each single-dimensional stochastic process element $\mathrm{d}w_{i}$ of $\boldsymbol{\mathrm{d}w}$ is an increment of a Wiener process, with $\langle\mathrm{d}w_{i}\rangle=0,\langle\mathrm{d}w_{i}^{2}\rangle=\mathrm{d}t,\forall i$. The discontinuous jump terms are contained in $\boldsymbol{\xi}\in\mathbb{R}^{2\times 2}$ and $\boldsymbol{\mathrm{d}J}\in\mathbb{N}^{2}$, where $\boldsymbol{\mathrm{d}J}$ represents a two-dimensional Poisson process. These are Poisson-distributed jumps with an average jump rate $\boldsymbol{\lambda}\in\mathbb{R}^{2}$ in unit time $t$. The average expected number of jumps of each jump process $J_{i}$ in a timespan $t$ is $\lambda_{i}t$. The jump amplitudes $\boldsymbol{\xi}$ are Gaussian distributed with zero mean and standard deviation $\xi_{i,j}$.

We here consider merely autonomous systems. Non-ergodic problems are beyond the scope of this work, and a more delicate approach to both bivariate stochastic processes would be needed.

III Bivariate diffusion processes

Let us begin with bivariate diffusion processes, for which the model takes the form

[TABLE]

The model consists of six functions, two for the drift coefficients and four for the diffusion coefficients. Given a bivariate diffusion process, can we reconstruct the aforementioned parameters strictly from data? Extensive work exists on this matter Friedrich et al. (2011), especially covering purely diffusion processes, and we will use these now as a stepping stone to jump-diffusion processes. Understanding the working and contingencies of reconstructing the parameters of a diffusion process (Eq. (2)) will serve as a gateway to understand how a similar procedure awards equivalent measures for jump-diffusion processes. We address the aforementioned question firstly by revisiting the mathematical foundation that allows one to recover, strictly from data, the drift $\boldsymbol{N}$ and diffusion $\boldsymbol{g}$ coefficients. Subsequently, we numerically integrate diffusion processes with a priori fixed values of the drift $\boldsymbol{N}$ and diffusion $\boldsymbol{g}$ coefficients and aim at retrieving these values strictly from the generated data (Euler–Mayurama scheme with a time sampling of $10^{-3}$ over a total of $10^{5}$ time units, i.e., $10^{8}$ number of data points). If the actual and retrieve values match, the reconstruction method is effective.

A stochastic process has a probabilistic description given by the master equation Risken (1996); Tabar (2019). It does not describe a specific stochastic processes in itself, but the probabilistic evolution of the process in time. The master equation accepts an expansion in terms, the KM expansion, that allows for a purely differential description of the process. More importantly, the coefficients of the expansion, known as KM coefficients, entail directly a relation to the aforementioned parameters of a stochastic process given by Eq. (1). The exact relation will be given below.

There is although an important detail in what regards the KM coefficients: they are in themselves not constants, but functions on the underlying space, or in other words, a scalar field, and for our purposes here, they can be understood as two-dimensional surfaces. We will denote these as KM surfaces.

Lastly, and more familiar, the Fokker–Planck equation is a truncation of the KM expansion at second order. It is especially relevant given its connection to physical processes and the Pawula theorem Pawula (1967). The Pawula theorem ensures that the truncation is not ill-suited for the underlying process if the fourth-order KM coefficient approaches zero in the limit $\mathrm{d}t\to 0$. It is now crucial to understand that the theorem holds for a one-dimensional process, and we are not aware of a proof for higher dimensions. This contrasts the common notion that studying only the first two KM coefficients of two- or higher-dimensional processes is sufficient (see Refs. Prusseit and Lehnertz (2007); Anvari et al. (2016); Lehnertz et al. (2018) and references therein).

The KM coefficients $\mathcal{M}^{[\ell,m]}(x_{1},x_{2})\in\mathbb{R}^{2}$ of orders $(\ell,m)$ are defined as:

[TABLE]

and can be obtained from bivariate time series $(y_{1}(t),y_{2}(t))$. Theoretically, $\Delta t$ should take the limiting case of $\Delta t\to 0$, but the restriction of any measuring or storing devices—or the nature of the observables themselves—permits only time-sampled or discrete recordings. The relevance and importance of adequate time-sampling was extensively studied and discussed in Ref. Lehnertz et al. (2018); Tabar (2019).

In the limiting case where $\Delta t$ is equivalent to the sampling rate of the data, the KM coefficients take the form

[TABLE]

The algebraic relations between the KM coefficients and functions in Eq. (2) are given by Anvari et al. (2016); Tabar (2019)

[TABLE]

An explicit derivation can be found in Appendix A. Evidently, this under-determined set of five equations is insufficient to uncover the six functions of a general stochastic diffusion process. One must bare this in mind, for the same issue will arise again when reconstructing jump-diffusion processes from data. Nonetheless, under certain assumptions it is possible to reduced the dimension of the problem and therefore obtain a system of equation which is not under-determined. Two methods for these cases are present in Ref. Tabar (2019) and another criterion will be presented later.

In order to relate the results obtained from studying the KM coefficients against the theoretical functions, we propose a method to assess the difference between the values of the theoretically expected functions and the estimated values of the KM coefficients. Since for bivariate processes the KM coefficients are two-dimensional—as are the parameters of Eq. (1)—an adequate “distance” measure between the resulting two-dimensional surfaces is required.

Following Ref. Prusseit and Lehnertz (2008), we propose a distance measure that allows for the variability of the density of data in some regions of the underlying space to be taken into consideration.

Let $f^{[\ell,m]}(y_{1},y_{2})$ denote the theoretical value for orders $(\ell,m)$ introduced in the model, i.e., a non-linear combination of the various parameters of the system. The distance between each surface can be defined as

[TABLE]

where $U$ denotes the domain of $\mathcal{M}^{[l,m]}(y_{1},y_{2})$. The least-squared distance volume $V$ between the surfaces is zero if $\mathcal{M}^{[l,m]}(y_{1},y_{2})=f^{[l,m]}(y_{1},y_{2})$. It is this volume that one aims to minimize such that the reconstructed KM coefficients match the underlying theoretical functions in the model. Since $\mathcal{M}^{[l,m]}(y_{1},y_{2})$ is a real-valued function measured over a distribution space $U$, the density of data points is not uniform over $U$. This implies that a comparative measure on distances between $\mathcal{M}^{[l,m]}(y_{1},y_{2})$ and $f^{[l,m]}(y_{1},y_{2})$ would be non-normalized to the density of points of the space. We therefore introduce a normalization to Eq. (6) that ensures the less dense areas of $U$ are normalized accordingly, thus mitigating the effect of scarcity of points at the borders of $U$ and an overestimation of $V$ due to outliers in the distribution. We derive such a normalization by considering the zeroth-order KM coefficient $\mathcal{M}^{[0,0]}(y_{1},y_{2})$ which captures exactly the density of points in $U$, although it is in itself not normalized as a distribution. The resulting normalized volume error measure $V_{\text{err}}$ between surfaces takes the form (state dependencies not explicit)

[TABLE]

where $p(\cdot)$ denotes the probability density. Coincidentally, the numerical evaluation implemented via both a histogram or a kernel-based estimators immediately yields this density. (i.e., the zeroth-power of the right-hand side of Eq. (3), before applying the estimation operator). This makes it easy to retrieve $p(y_{1},y_{2})$ as one numerically evaluates data.

With this at hand, it is now possible to relate theoretical and numerical results and to quantify the deviation of the obtained KM coefficients from the functions employed.

To showcase what two-dimensional KM coefficients are as well as how to identify drift and diffusion terms of bivariate diffusion processes, we present in the following two exemplary processes with a priori known coefficients. In this manner, by employing Eqs. (4) and (5), one can judge the outcome of the KM coefficient estimation procedure from discrete data in comparison with the expected theoretical functions.

We begin with two uncoupled processes, where one has constant diffusion and a quartic potential as the drift term:

[TABLE]

In Fig. 1, we show the corresponding KM coefficients $\mathcal{M}^{[1,0]}$, $\mathcal{M}^{[0,1]}$, $\mathcal{M}^{[1,1]}$, and $\mathcal{M}^{[2,0]}$ together with the theoretically expected functions. The per-design cubic-linear function ($N_{1}=-x_{1}^{3}+x_{1}$) acting as drift coefficient along the first dimension as well as the negatively-sloped surface of $N_{2}=-x_{2}$ are evident. Likewise, the constant diffusion term leads to a flat constant-valued $\mathcal{M}^{[2,0]}$, and the absence of any non-diagonal elements ($g_{1,2}=g_{2,1}=0$) agrees with the zero-valued $\mathcal{M}^{[1,1]}$. Alongside the surfaces are plotted the theoretically expected values, which agree well with the data recovery.

We next extent Eq. (8) by adding multiplicative noise to the diffusion term and by including a noise coupling term $g_{1,2}\neq 0$

[TABLE]

The recovered KM coefficients (see Fig. 2) of the drift terms remain unaltered, but as posited, the second-order KM coefficients, i.e., $\mathcal{M}^{[1,1]}$, $\mathcal{M}^{[2,0]}$, $\mathcal{M}^{[0,2]}$, and $\mathcal{M}^{[2,2]}$, clearly exhibit the influence of the multiplicative noise. The quadratic multiplicative dependencies of $\mathcal{M}^{[2,0]}$ and $\mathcal{M}^{[0,2]}$ and their offsets from zero are evident. More pertinently, one can notice $\mathcal{M}^{[1,1]}$ to display the expected shape arising from Eq. (5), i.e., this value is non-zero and exhibits the parabolic shape of $g_{1,2}\cdot g_{2,2}=0.5(0.2+2x_{2}^{2})$. For $x_{2}=0$, the minimum of $\mathcal{M}^{[1,1]}$ coincides with $0.1$, as expected. This indicates that the presence of the multiplicative noise does not hinder the assertion of the KM coefficients. Again, the recovered KM coefficients match the theoretically ones.

IV Bivariate jump-diffusion processes

The KM coefficients of bivariate jump-diffusion processes take the following form (under the parameter prescription used in the jump-diffusion model in Eq. (1)) Anvari et al. (2016); Lehnertz et al. (2018); Tabar (2019)

[TABLE]

[TABLE]

where $\langle\xi_{ij}^{2\ell}\rangle=s_{ij}^{\ell}$ are the variances of the Gaussian-distributed jump amplitudes. An extended derivation can be found in Appendix A. The last equations here are taken from the general form

[TABLE]

IV.1 Understanding the impact of jumps

As an illustrative case study, we investigate a general jump-diffusion process that is based on Eq. (9) but excludes the multiplicative diffusion terms. Taking into account the effect of the jump terms, but maintaining the system independent in at least one of the dimensions, we extend Eq. (9) to include jumps only in the diagonal terms of $\boldsymbol{\xi}$

[TABLE]

where for the present case, $\phi=0.0$. In this manner, jumps are added to the first dimension of the process, having an amplitude of $\xi_{1,1}=0.2$ and occurring every $0.1t$, given $\lambda_{1}=0.1$. Similarly, jumps are added to the second dimension, $\xi_{2,2}=0.1$, but the jumps occur three times more often than the aforementioned, given $\lambda_{2}=0.3$. The influence of jumps can be observed across all KM coefficients (see Fig. 3). The previously smooth KM surfaces become rugged from the fast variations emerging due to the jumps, and the higher-order KM coefficients—although small compared to the lower-order ones—clearly do not vanish. This indicates that the continuous stochastic modeling of time series of complex systems (white-noise-driven Langevin equation) together with the commonly used assumptions stemming from the Pawula theorem are invalid for jump-diffusion processes. Modeling these processes with only the first two orders of the KM expansion of the master equation is therefore insufficient.

In order to understand further if it is possible to uncover the jump amplitude terms of coupled processes, we use the previous model Eq. (15) with $\phi=0.3$, thereby effectively introducing a coupling via the off-diagonal elements of the jump matrix $\boldsymbol{\xi}$. We show, in Fig 4, the corresponding fourth-order KM coefficients. The impact of the coupling is visible, although small, in $\mathcal{M}^{[4,4]}$, which is no longer zero. Likewise, $\mathcal{M}^{[4,0]}$ and $\mathcal{M}^{[0,4]}$ also do not vanish. In Appendix B, we present the corresponding KM coefficients up to order eight.

IV.2 Criteria for recovering coefficients in diffusion and jump-diffusion models

For the case of vanishing off-diagonal elements $g_{2,1}$ and $\xi_{1,2}$, we can identify ways to recover the remaining coefficients of these processes.

First, given that the noise $\boldsymbol{\mathrm{d}}\boldsymbol{\omega}$ is Gaussian distributed, $\boldsymbol{g}$ is sign-reversal symmetric and one can thus assume that it takes only positive values. One obtains that, if $\mathcal{M}^{[1,1]}=0$ then at least two elements of $\boldsymbol{g}$ must be zero, and if $\mathcal{M}^{[2,2]}=0$ then at least two elements of $\boldsymbol{\xi}$ must be zero (by assuming that $\lambda_{1}$ and $\lambda_{2}$ are non-vanishing rates). These findings reduce the dimensionality of the estimation procedure and ensure that the underlying processes are less complex that the full-fledged description of Eq. (1), although they do not grant which coefficients are zero-valued.

Second, if one either employs a heuristic argument of independence of the jump processes or neglects the off-diagonal jump amplitudes $\xi_{1,2}$ and $\xi_{2,1}$ (e.g., by assuming they are small compared to the diagonal terms of $\boldsymbol{\xi}$), one finds the following approximations:

[TABLE]

Likewise, the jump rates $\lambda_{1}$ and $\lambda_{2}$ can be obtained equivalently as

[TABLE]

Taking again model Eq. (15) with $\phi=0.0$ and following Eq. (16), we obtain

[TABLE]

These estimated values (indicated by the superscript $\cdot^{\mathrm{est}}$) are very close to the actual ones. The criteria and approximations are especially relevant when constructing or analyzing systems which are known to have a specific unidirectional coupling form, e.g. a master-slave system, where for example the noise or the slave system are dictated by the driving master system.

IV.3 Factors influencing quality of recovery of coefficients

In order to validate the quality of the non-parametric recovery of the KM coefficients, we now turn to two critical aspects: firstly, bivariate processes may require a high number of data points in a time series for the estimation to be reliable; secondly, the interplay between the drift, diffusion, and jump parts of a stochastic processes may render the estimation incorrect.

Addressing these aspects, we include a more contrived model involving couplings and interactions in both the diffusion and jump terms, thus theoretically resulting in having all higher-order KM coefficients non-zero, and especially the KM coefficients with $\ell\neq 0$ and $m\neq 0$. The parameters for the model read:

[TABLE]

The parameters $\alpha$ (diffusion term) and $\beta$ (jump term) can be freely varied.

Let us focus firstly on the number of data points in a time series. We utilize the models Eq. (15) with $\phi=0.3$ and Eq. (18) with $\alpha=\beta=0.3$, and show in Figs. 5 and 6, respectively, the error volumes $V_{\text{err}}$ for the KM coefficients for an increasing number of data points.

The reliability of the recovery of KM coefficient is valid for higher amount of data ($n\geq 10^{5}$), as expected, although the scarcity of data posits no extensive problem for the calculation. It is especially important to notice that a time series with a lower amount of data entails naturally fewer jumps in the process, hindering the possibility of accurately recovering the jump terms from such short time series. For $n\geq 10^{6}$, the estimation seems reliable, the standard deviations become minute, and most error values approach zero, i.e., the theoretical and estimated KM surfaces are close.

One remark is necessary on the recovery of the drift terms. The presence of noise and jumps in the process takes its toll on the recovery of the exact form of the KM coefficients as well as the explicit dependence of the state variables, i.e., the quartic potential in both Eqs. (15) and (18). A finer time sampling can help to improve the results.

To further test the limitation of retrieving the KM coefficients from data, we utilize model Eq. (18) once more and investigate the influence of parameters $\alpha$ (diffusion term) and $\beta$ (jump term). For increasing $\alpha$ ($\alpha\in[10^{-2},10^{2}]$) and $\beta=0.3$, we observe a considerable impact on the error volume $V_{\text{err}}$ after the order of magnitude on the diffusion parameter $\alpha$ is ten-fold bigger in comparison to the diffusion parameter $g_{1,2}$ (Fig. 7). Similarly, for increasing $\beta$ ($\beta\in[10^{-2},10^{2}]$) and $\alpha=0.3$, the error volume $V_{\text{err}}$ is considerably impacted already when $\beta$ is of similar size as the other parameters, namely $\xi_{1,2}=0.5$ (Fig. 8).

These findings point to the difficulty of recovering the KM coefficients in the presence of jumps. Nonetheless, our findings indicate that the current understanding, modeling, and numerical recovery of KM surfaces, for the case of jumps of comparable size to the diffusion terms, is possible and reliable.

Before closing this section, we note that the calculations of the KM coefficients are computationally inexpensive Gorjão and Meirinhos . E.g., estimating all the 14 KM coefficients, $\mathcal{M}^{[1,0]},\mathcal{M}^{[0,1]},\dots,\mathcal{M}^{[0,8]}$ as presented in Fig. 9 takes about 5.5 s on a desktop computer (quad-core with clock speed 2.20 GHz) for a bi-dimensional time series of size $2\times 10^{6}$ data points (and 22 s for $2\times 10^{7}$ data points). Our approach might thus be advantageous for field applications that aim at an investigation of interactions between complex systems with poorly understood dynamics.

V Conclusion

We introduced the bivariate jump-diffusion process which comprises two-dimensional diffusion and two-dimensional jumps that can be coupled to one another. For such a process we presented a data-driven, non-parametric estimation procedure of higher-order Kramers–Moyal coefficients and investigated its pros and cons using synthetic bivariate time series from continuous and discontinuous processes. The procedure allows one to reconstruct relevant aspects of the underlying jump-diffusion processes and to recover the underlying parameters.

Having now a traceable mathematical framework, the model can be extended to embody other noise and jump properties. An extension from the underlying Wiener process to include e.g. fractional Brownian motion is straightforward but nevertheless requires further investigations to derive an explicit forward Kolmogorov equation Tabar (2019). Also, a generalization to continuous jump processes—originating from alpha-stable or other heavy-tailed distributions (Lévy noise-driven Langevin equation)—is possible, however, with the drawback that calculating the conditional moments may not always be mathematically possible Tabar (2019). On the other hand, a numerical estimation of generalised moments should be possible but these still require a physical interpretation.

We are confident that our novel approach provides a general avenue to further understanding of interacting complex systems (e.g. brain or power grids Lehnertz et al. (2018); Schäfer et al. (2018); Gorjão et al. (2019); Anvari et al. (2019)) whose dynamics exhibit nontrivial noise contributions.

Acknowledgements

The authors would like to thank Thorsten Rings and Mehrnaz Anvari for interesting discussions, and Francisco Meirinhos for the tremendous help in devising the methodology behind the results. L.R.G. thanks Dirk Witthaut and Giulia di Nunno for everlasting support, and gratefully acknowledges support by the grant No. VH-NG-1025 and a scholarship from the E.ON Stipendienfonds.

Appendix A Extended derivation of the two-dimensional Kramers–Moyal coefficients for a jump-diffusion process

The following derivations stem from Eq. (3) and apply to the two-dimensional jump-diffusion process $(y_{1},y_{2})$, as in Eq. (1). All orders of the Kramers–Moyal coefficients $\mathcal{M}^{[\ell,m]}$ are $(\ell,m)\in\mathbb{N}_{+}$.

A.0.1 Kramers–Moyal coefficients $\mathcal{M}^{[1,0]}$ and $\mathcal{M}^{[0,1]}$

[TABLE]

where $\langle g_{i,j}\mathrm{d}W_{j}\rangle=\langle g_{i,j}\rangle\langle\mathrm{d}W_{j}\rangle=0$, because a Wiener process has the property $\langle\mathrm{d}W_{j}\rangle=0$. Further $\langle\xi_{i,j}\mathrm{d}J_{j}\rangle=\langle\xi_{i,j}\rangle\langle\mathrm{d}J_{j}\rangle=0$, since $\xi_{i,j}$ is a Gaussian with zero mean, i.e., $\langle\xi_{i,j}\rangle=0$.

Mutatis mutandis for $\mathcal{M}^{[0,1]}$.

A.0.2 Kramers–Moyal coefficient $\mathcal{M}^{[1,1]}$

[TABLE]

where higher-order terms $\mathcal{O}(\mathrm{d}t)^{\epsilon}$, with $\epsilon>0$, vanish in the limit $\mathrm{d}t\to 0$. Recall as well $\langle(\mathrm{d}w_{i})^{2}\rangle=\mathrm{d}t$.

A.0.3 Kramers–Moyal coefficients $\mathcal{M}^{[2,0]}$ and $\mathcal{M}^{[0,2]}$

[TABLE]

using the previously employed nomenclature $\langle\xi_{ij}^{2}\rangle=\sigma_{\xi_{ij}}^{2}=s_{ij}$ as well as $\langle(\mathrm{d}J_{i})^{2}\rangle=\lambda_{i}\mathrm{d}t$.

Mutatis mutandis for $\mathcal{M}^{[0,2]}$

[TABLE]

A.0.4 Kramers–Moyal coefficient $\mathcal{M}^{[2,2]}$

[TABLE]

Terms including $\mathrm{d}t$ on the right-hand side of the above equation vanish for $\mathrm{d}t\to 0$, where as well $\langle\xi_{1,1}\xi_{1,2}\rangle=\langle\xi_{1,1}\rangle\langle\xi_{1,2}\rangle=0$, and $\frac{1}{\mathrm{d}t}[\langle(\mathrm{d}J_{1})^{2}\rangle\langle(\mathrm{d}J_{2})^{2}\rangle]=\frac{1}{\mathrm{d}t}[\lambda_{1}\mathrm{d}t\lambda_{2}\mathrm{d}t]\propto\mathrm{d}t$ vanishes in the limit $\mathrm{d}t\to 0$.

A.0.5 Kramers–Moyal coefficients $\mathcal{M}^{[\ell,m]}$, for $2\times(\ell,m)\geq 2$

For $(2\ell,2m)$, with $(\ell,m)\geq 4$, the Kramers–Moyal coefficients $\mathcal{M}^{[2\ell,2m]}$ are as follows

[TABLE]

In the last step, take the fact that the jump amplitudes $\xi_{i,j}$ are Gaussian distributed, thus $\langle\xi_{i,j}^{2\ell}\rangle\propto\sigma_{\xi_{i,j}}^{2\ell}=s_{i,j}^{\ell}$. In this manner, all Kramers–Moyal coefficients $\mathcal{M}^{[2\ell,2m]}$, with $(\ell,m)\geq 1$ are obtained.

Appendix B Extended results for modeled data by Eq. (15)