# EvolveMT: an Ensemble MT Engine Improving Itself with Usage Only

Kamer Ali Yuksel, Ahmet Gunduz, Mohamed Al-Badrashiny,
Shreyas Sharma, and Hassan Sawaf
aiXplain Inc., 16535 Grant Bishop Lane, Los Gatos, CA 95032, US
{kamer, ahmet, mohamed, shreyas<EMAIL_ADDRESS>
###### Abstract

This paper presents EvolveMT for efficiently combining multiple machine
translation (MT) engines. The proposed system selects the output from a single
engine for each segment by utilizing online learning techniques to predict the
most suitable system for every translation request. A neural quality
estimation metric supervises the method without requiring reference
translations. The online learning capability of this system allows for dynamic
adaptation to alterations in the domain or machine translation engines,
thereby obviating the necessity for additional training. EvolveMT selects a
subset of translation engines to be called based on the source sentence
features. The degree of exploration is configurable according to the desired
quality-cost trade-off. Results from custom datasets demonstrate that EvolveMT
achieves similar translation accuracy at a lower cost than selecting the best
translation of each segment from all translations using an MT quality
estimator. To our knowledge, EvolveMT is the first meta MT system that adapts
itself after deployment to incoming translation requests from the production
environment without needing costly retraining on human feedback.

##  1 Introduction

Machine Translation (MT) has experienced substantial progress in recent years,
resulting in improving accuracy and more human-like translation output.
Despite these advancements, challenges remain, particularly in ensemble
modeling. Ensemble models integrate predictions from multiple individual
models to achieve a more accurate final output. However, the effective
combination of these models is often a complex task that requires thoughtful
consideration of factors such as the model architecture, training data, and
prediction combination methods. One of the significant challenges in MT
ensembling is that the training data used to train the ensemble model, may not
be fully representative of the data to be translated later, leading to a
mismatch between the model and the data. This paper presents EvolveMT, a
method that addresses data drift in ensemble models by continual self-
adaptation for optimal performance during usage.

In the subsequent section, we review existing machine translation (MT) quality
estimation metrics in the literature, which have been trained on human
evaluation or post-editing datasets. In the Approach section, we present a
comprehensive explanation of the proposed method. In the Experiments section,
we describe our experimental design and provide quantifiable results
demonstrating the enhancement resulting from the application of the proposed
method, as compared to state-of-the-art quality estimation metrics. Finally,
we discuss the obtained results and present our conclusions.

##  2 Related Work

In the WMT20 Metrics Shared Task (Mathur et al., 2020), four reference-free
metrics were submitted to evaluate machine translation outputs in the news
translation task. These metrics use bilingual mapping of contextual embeddings
from language models such as XLM-RoBERTa (Conneau et al., 2019) to assess
cross-lingual semantic similarity. However, they often struggle to accurately
differentiate between human and machine translations, except for COMET-QE (Rei
et al., 2020), the only reference-free metric capable of doing so.

The study by Freitag et al. (2021a) evaluated top MT systems from WMT 2020
using Multidimensional Quality Metrics (MQM) and professional translator
annotations. Their results showed a low correlation between crowd worker
evaluations and MQM, leading to different rankings and questioning previous
conclusions. The study also found that automatic metrics based on pre-trained
embeddings can outperform human crowd workers, suggesting that models trained
with crowd-sourced human evaluations may have higher accuracy.

The WMT21 Metrics Shared Task (Freitag et al., 2021b), used MQM expert-based
human evaluation to acquire reliable ratings, and evaluate metrics on news and
TED talk translations produced by MT systems. Results showed reference-free
metrics COMET-QE and OpenKiwi performed well in scoring human translations but
not as well with MT outputs, and were strong at segment-level human
translation evaluation while competitive with reference-based metrics in
system-level evaluation.

REGEMT (Štefánik et al., 2021) is a reference-free metric in WMT21 that uses
an ensemble of surface, syntactic, and semantic similarity metrics as input to
a regression model. As demonstrated by CushLEPOR (Han et al., 2021), it allows
customization, outperforming lexical semantic similarity-based metrics with a
higher computational cost.

Onception (Mendonça et al., 2022) used active learning to converge an MT
ensemble in a production environment to the best MT with evaluations acquired
online.

(Naradowsky et al., 2020) used bandit-learning to adapt MT policies based on
simulated user feedback, outperforming the best single MT in mixed-domain
settings. A contextual bandit strategy was proposed to make instance-specific
decisions, but the system still required a human-in-the-loop (HITL) process.

##  3 Approach

EvolveMT is a quasi MT ensemble technique. In contrast to the traditional
multi-system MT approach, which combines outputs from multiple MT systems to
enhance translation accuracy and fluency, EvolveMT prioritizes the selection
of the most optimal translation from a finite set of MT systems, as we
demonstrate in this section. Figure 1 below shows the system architecture of
EvolveMT. The system is centered around a multi-class classification model
that drives multiple processes to select the best MT model for translation
requests.

For each incoming machine translation request, we use SpaCy (Honnibal and
Montani, 2017) and Stanza (Qi et al., 2020) frameworks to extract
morphological and lexical features. These features include the count of
tokens, characters, and the average word length, as well as the frequency of
Part-of-Speech labels (such as nouns, verbs, adjectives, etc.), the frequency
of Named Entity Recognition labels (including entities such as persons,
locations, organizations, etc.), and the frequency of morphological features
(e.g. gender and aspect). These features are combined with the
1024-dimensional embedding vector generated by the XLM-RoBERTa encoder of the
COMET-QE model and stored alongside the input sentence in the _Ranked Batch
Requests Queue_. This queue serves the purpose of prioritizing translation
requests that necessitate precedence in processing.

At the outset, requests in the Ranked Batch Requests Queue are ranked based on
the order in which they are added. The highest-ranked Machine Translation (MT)
request is selected for translation. The _Multi-class MT Classifier_ employs
the extracted features of the selected MT request to determine the MT systems
to be utilized. The classifier prioritizes MT systems with a higher
probability of having a higher COMET-QE value. Exploration of additional MT
systems becomes more likely only if the probabilities from the classifier’s
prediction exhibit high entropy. This enables EvolveMT to minimize the cost of
exploration when the best MT is predicted with high certainty.

Finally, the selected MT systems are utilized to translate the MT request, and
the COMET-QE score is calculated for each translation. The translation with
the highest score is chosen and returned in response to the MT request. The
Multi-class MT Classifier is then updated online with the best MT system, as
determined by the COMET-QE score, serving as a label for the extracted
features of the MT request. The online machine learning (ML) functionality of
FLAML AutoML framework (Wang et al., 2021) is utilized for online learning.
This capability enables the optimization of model hyper-parameters during the
iterative course of ML, facilitating continual ML without repetitively hyper-
tuning the classifier from scratch.

Subsequently, the classifier is employed to re-rank the Ranked Batch Requests
Queue, based on the uncertainty of the classifier, with requests having higher
entropy being placed at the top of the queue for prioritized translation.
Getting the MT request with the maximum entropy from the queue after each
iteration, helps prioritize the most informative sample for the iterative
training of the classifier. As the classifier improves its ability to predict
the best MT model for MT requests via learning, it reduces the likelihood of
exploring other MT(s).

The driving algorithm in EvolveMT, which outlines the primary process of the
proposed method, is presented in pseudo-code form in Algorithm 1.

Algorithm 1 EvolveMT with online active-learning

0: $MTQueue$: list of tuples (source text, features),
and $MaxMTs$: maximum number of MTs to sample

while $len(MTQueue)$ do

$Classifier.rankByUncertainty(MTQueue)$

$source,feats\leftarrow MTQueue.popMaxEntropyItem()$

$predMT,classProbs\leftarrow Classifier.predict(feats)$

$predTrans\leftarrow Translate(source,predMT)$

$randMTs\leftarrow sampleMTs(classProbs,MaxMTs)$

$maxEnt\leftarrow normalizedEntropy(classProbs)$

if $randMTs_{0}=predMT$ and $maxEnt<\alpha$ then

$Classifier.learn(feats,predMT)$

else

$sampled\leftarrow Translate(source,randMTs)$

$randScores\leftarrow CometQE(source,sampled)$

$predMTScore\leftarrow CometQE(source,predTrans)$

if $max(randScores)>predMTScore$ then

$Classifier.learn(feats,randMTs)$

$IndexOfBestMT\leftarrow randScores.argMax()$

$predTrans\leftarrow sampled_{IndexOfBestMT}$

else

$Classifier.learn(feats,predMT)$

end if

end if

$respondMTRequest((source,predTrans))$

end while

Figure 1: EvolveMT System Architecture

##  4 Experiments

###  4.1 Data

A multi-lingual corpus of 37,500 human-translated sentences in Czech, German,
and Russian, along with their corresponding English source-texts, was
collected for the OPUS repository (Lison and Tiedemann, 2016; Aulamo and
Tiedemann, 2019) using stratified random sampling for each language and
dataset. To evaluate EvolveMT, translations for each sentence in the corpus
were obtained from one open-source machine translation system (Tiedemann et
al., 2022) and five major machine translation service providers in the
industry (Google, Azure, AWS, ModernMT, and AppTek).

###  4.2 Setup

Experiments were conducted on a 64-bit Ubuntu 22.04 LTS computer system with
an AMD Ryzten 5950X CPU (16 processors, 32 threads) and 64GB of memory. An
Nvidia GeForce RTX GPU was used for XLM-RoBERTa embedding extraction from the
fine-tuned COMET-QE encoder. The results showed an average 2.88 ($\pm$0.06)
millisecond response time for the EvolveMT system to return an MT output and
update its classifier when the MT request queue contained a single item.
Depending on GPU usage, the feature extraction time was (183.43 -244.80)
milliseconds. It’s worth noting that in a production setting, feature
extraction can be performed in parallel for multiple MT requests.

###  4.3 Evaluation criteria and baselines

In this paper, grid search is used to evaluate the impact of two
hyperparameters, $maxMTs$ and $\alpha$, on the classifier’s performance.
$maxMTs$ refers to the maximum number of machine translation systems the
classifier can select, and $\alpha$ is the maximum entropy threshold (as
described in Algorithm 1). The grid search involves varying $maxMTs$ from 1 to
6 (the maximum number of the MT systems we are using), and $\alpha$ from 0.1
to 1.0 with increments of 0.1. The experimental results are obtained by
averaging 100 repetitions to account for the method’s inherent stochasticity.
For clarity in the results section, we present the results over the $maxMTs$
range while setting $\alpha$ to its optimal value of 0.2, determined through
the grid search.

For the evaluation, we adopt the reference-based quality score COMET-DA, as
detailed in (Rei et al., 2020), as the evaluation metric for our ensemble
output. This choice is motivated by the results of prior research which have
demonstrated that COMET-DA exhibits a higher correlation with human evaluation
scores compared to other widely used machine translation metrics, such as BLEU
and METEOR. The evaluation of EvolveMT is being conducted against the
following baselines:

  * •

COMET-QE Ensemble: translation is performed utilizing the six MT systems. The
translation with the highest COMET-QE score is selected for each input
sentence as the ensemble translation. Then, the COMET-DA score is calculated
using selected translations.

  * •

Best MT: involves translating the entire data using all six MT systems. The MT
system that produces the highest overall COMET-DA is then selected as the Best
MT to employ.

###  4.4 Results

The comparison of COMET-DA scores of the COMET-QE ensemble and Best MT
concerning various variants of EvolveMT with varying $MaxMTs$ values are
presented in Table 1. The results are depicted for the three language pairs of
English-to-Czech, English-to-German, and English-to-Russian. In addition, the
average translation cost of each system across the three languages is also
documented in the table. The findings indicate that EvolveMT approximates the
COMET-QE ensemble’s quality while incurring significantly lower costs.

Furthermore, the results in the table reveal that the optimal cost-quality
trade-off for EvolveMT varies depending on the target language. Specifically,
for all three language pairs, it can be observed that EvolveMT with $MaxMTs=3$
and $MaxMTs=4$ provide the best balance between cost and quality when compared
to other individual and ensemble MTs. As $MaxMTs$ increases, EvolveMT can
achieve higher MT quality by exploring a larger pool of MT systems from which
the best translation can be selected. Hence, the $MaxMTs$ parameter can be
adjusted to achieve the desired cost-quality trade-off.

Notably, after only a few hundred Machine Translation (MT) requests from the
total dataset, the EvolveMT algorithm demonstrates convergence towards an
upper limit of its weighted F1-score, which depends on the parameter $maxMTs$.
Figure 2 shows the confusion matrix between the outputs of EvolveMT (with
$MaxMTs$ = 4) and the COMET-QE ensemble after 100 translations requests. The
swift convergence of EvolveMT with a limited number of requests is mainly due
to the utilization of XLM-RoBERTa embeddings that have been fine-tuned
specifically for the COMET-QE task. This exemplifies the model’s
effectiveness, as it begins with no prior knowledge, and within a few hundred
requests, it can converge and approach the performance of the COMET-QE
ensemble. It is crucial to mention that the results presented in Table 1
encompass the "warm up" phase where EvolveMT starts from zero knowledge until
full convergence is achieved. If this phase were excluded, the COMET-DA scores
of EvolveMT would likely be even higher.

Figure 2: The normalized confusion matrix between EvolveMT ($MaxMTs$ = 4) and COMET-QE after 100 translations requests. |  | COMET-DA
---|---|---
Model | Cost ($) | English-to-Czech | English-to-German | English-to-Russian
Best MT (1) | 20.000 | 0.867 | 0.586 | 0.617
COMET-QE (6) | 77.000 | 0.900 | 0.605 | 0.658
EvolveMT (1) | 12.312 | 0.851 | 0.567 | 0.605
EvolveMT (2) | 23.442 | 0.870 | 0.586 | 0.627
EvolveMT (3) | 32.358 | 0.878 | 0.591 | 0.637
EvolveMT (4) | 39.905 | 0.882 | 0.596 | 0.643
EvolveMT (5) | 46.067 | 0.887 | 0.598 | 0.647
EvolveMT (6) | 51.095 | 0.887 | 0.599 | 0.651
Table 1: The Cost and COMET-DA comparison of the Best MT system, COMET-QE and
EvolveMT ensembles for various $MaxMTs$ parameters (indicated in parentheses).
The MT quality increases as COMET-DA scores increase

##  5 Discussion

The cost-benefit analysis of EvolveMT highlights the trade-off between run-
time efficiency and training expenses. While the run-time cost of EvolveMT may
be higher than that of Best MT, it does not require the extensive and time-
consuming training process required for traditional MT ensemble methods. This
training process involves obtaining translations from all MTs and scoring them
using references generated by annotators.

However, the increased run-time cost of EvolveMT is offset by its ability to
achieve superior production quality and adapt to changes in the data domain
with a minimum amount of overhead. As the data domain changes, traditional MT
ensemble techniques require costly retraining to accommodate the new domain,
whereas EvolveMT can adapt to changes with a few hundred MT requests.

This versatility and adaptability of EvolveMT make it a robust solution for
machine translation tasks that may be subject to data variation, as it can
adjust to these changes with minimal effort. The cost-benefit analysis results
clearly demonstrate that the increased run-time cost of EvolveMT is outweighed
by its high performance and adaptability in the face of changing data domains.

##  6 Limitations

The performance of the EvolveMT system is contingent upon the reliability of
the COMET-QE model in providing accurate labels for the MT requests. Utilizing
the encoder’s embeddings as features necessitates that the COMET-QE model
performs effectively on blind MT requests. The batch re-ranking of MT requests
after each learning step may result in a computational bottleneck if the queue
size is substantial. To mitigate this issue, an asynchronous re-ranking
process could be implemented, whereby the queue is only reorganized once the
re-ranking is completed. Additionally, before the re-ranking process, a
diverse subset of the queue can be selected based on the XLM-RoBERTa
embeddings, which reflect the novelty of the requests relative to previously
processed MT requests. The source embeddings from the XLM-RoBERTa model can be
cached in parallel during the batch feature extraction process utilizing GPU
capabilities, thus facilitating efficient COMET-QE inference. EvolveMT could
also be optimized for cost-effectiveness by incorporating the cost of each MT
in the ensemble into the algorithm.

##  7 Conclusion and Future Work

This study presents a novel approach called EvolveMT for ensembling machine
translation (MT) engines, focusing on minimizing the number of engines
required to be queried to achieve optimal quality. To evaluate the efficacy of
the proposed method, a series of experiments were conducted, wherein EvolveMT
was implemented with varying levels of granularity in terms of the maximum
number of engines permitted for each individual MT request. The quantitative
results of the experiments indicate that, compared to the traditional method
of querying all available MT engines, EvolveMT offers a more cost-effective
solution for the ensembling process without compromising the quality of the
resulting translations.

EvolveMT presents a unique advantage in terms of cost efficiency compared to
COMET-QE Ensemble. This is achieved by utilizing a stochastic exploration
approach that selectively queries additional MT engines based on predicted
probabilities, which are also employed in an active-learning framework by re-
ranking MT requests after each learning step. Furthermore, unlike traditional
MT ensemble techniques, EvolveMT can adapt in real-time to changes in
customers’ translation requests, without incurring the cost of acquiring human
references or undergoing costly re-training or fine-tuning.

In conclusion, this paper presents four significant contributions to the field
of machine translation: (1) the introduction of the first self-improving MT
system that operates without the need for human feedback; (2) the capability
of adaptively optimizing the MT ensemble in response to production environment
translation requests through online machine-learning; (3) the development of a
novel approach for selectively querying MT engines rather than relying on
translations from all available engines; and (4) the implementation of an
active-learning framework that leverages uncertainties from the ensemble for
batch translation.

##  8 Ethics and Impact Statement

EvolveMT is a high-quality machine translation (MT) system for individuals or
organizations. It can improve translation accuracy if validated on a specific
MT corpus. EvolveMT is trained from scratch for each customer or project,
eliminating biases in the algorithm but may still present biases in the
quality estimation metric or training dataset. The system is self-adaptable,
secure, and protects user privacy by deleting data immediately after
translation. EvolveMT eliminates the need for re-training and re-hypertuning,
reducing computational costs and being environmentally friendly. The only
potential harm is to linguists who perform post-editing as it reduces their
dependence on references or evaluations.

## References

  * Aulamo and Tiedemann (2019) Mikko Aulamo and Jörg Tiedemann. 2019.  The OPUS resource repository: An open package for creating parallel corpora and machine translation services.  In _Proceedings of the 22nd Nordic Conference on Computational Linguistics_ , pages 389–394, Turku, Finland. Linköping University Electronic Press.
  * Conneau et al. (2019) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Unsupervised cross-lingual representation learning at scale.  _arXiv preprint arXiv:1911.02116_.
  * Freitag et al. (2021a) Markus Freitag, George Foster, David Grangier, Viresh Ratnakar, Qijun Tan, and Wolfgang Macherey. 2021a.  Experts, errors, and context: A large-scale study of human evaluation for machine translation.  _Transactions of the Association for Computational Linguistics_ , 9:1460–1474.
  * Freitag et al. (2021b) Markus Freitag, Ricardo Rei, Nitika Mathur, Chi-kiu Lo, Craig Stewart, George Foster, Alon Lavie, and Ondřej Bojar. 2021b.  Results of the wmt21 metrics shared task: Evaluating metrics with expert-based human evaluations on ted and news domain.  In _Proceedings of the Sixth Conference on Machine Translation_ , pages 733–774.
  * Han et al. (2021) Lifeng Han, Irina Sorokina, Gleb Erofeev, and Serge Gladkoff. 2021.  cushlepor: customising hlepor metric using optuna for higher agreement with human judgments or pre-trained language model labse.  In _Proceedings of the Sixth Conference on Machine Translation_ , pages 1014–1023.
  * Honnibal and Montani (2017) Matthew Honnibal and Ines Montani. 2017.  spaCy 2: Natural language understanding with Bloom embeddings, convolutional neural networks and incremental parsing.  To appear.
  * Lison and Tiedemann (2016) Pierre Lison and Jörg Tiedemann. 2016.  OpenSubtitles2016: Extracting large parallel corpora from movie and TV subtitles.  In _Proceedings of the Tenth International Conference on Language Resources and Evaluation (LREC’16)_ , pages 923–929, Portorož, Slovenia. European Language Resources Association (ELRA).
  * Mathur et al. (2020) Nitika Mathur, Johnny Wei, Markus Freitag, Qingsong Ma, and Ondřej Bojar. 2020\.  Results of the wmt20 metrics shared task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 688–725.
  * Mendonça et al. (2022) Vânia Mendonça, Ricardo Rei, Luisa Coheur, and Alberto Sardinha. 2022\.  Onception: Active learning with expert advice for real world machine translation.  _arXiv preprint arXiv:2203.04507_.
  * Naradowsky et al. (2020) Jason Naradowsky, Xuan Zhang, and Kevin Duh. 2020.  Machine translation system selection from bandit feedback.  _arXiv preprint arXiv:2002.09646_.
  * Qi et al. (2020) Peng Qi, Yuhao Zhang, Yuhui Zhang, Jason Bolton, and Christopher D. Manning. 2020\.  Stanza: A python natural language processing toolkit for many human languages.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ , pages 101–108, Online. Association for Computational Linguistics.
  * Rei et al. (2020) Ricardo Rei, Craig Stewart, Ana C Farinha, and Alon Lavie. 2020.  COMET: A neural framework for MT evaluation.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 2685–2702, Online. Association for Computational Linguistics.
  * Štefánik et al. (2021) Michal Štefánik, Vít Novotnỳ, and Petr Sojka. 2021.  Regressive ensemble for machine translation quality evaluation.  _arXiv preprint arXiv:2109.07242_.
  * Tiedemann et al. (2022) Jörg Tiedemann, Mikko Aulamo, Daria Bakshandaeva, Michele Boggia, Stig-Arne Grönroos, Tommi Nieminen, Alessandro Raganato, Yves Scherrer, Raul Vazquez, and Sami Virpioja. 2022.  Democratizing machine translation with opus-mt.
  * Wang et al. (2021) Chi Wang, Qingyun Wu, Markus Weimer, and Erkang Zhu. 2021.  Flaml: a fast and lightweight automl library.  _Proceedings of Machine Learning and Systems_ , 3:434–447.


# Enhancing the accuracy of a data-driven reconstruction of bivariate jump-
diffusion models with corrections for higher orders of the sampling interval

Esra Aslim1,2, Thorsten Rings1,2, Lina Zabawa1,2 and Klaus Lehnertz1,2,3 1
Department of Epileptology, University of Bonn Medical Centre, Venusberg
Campus 1, 53127 Bonn, Germany  2 Helmholtz Institute for Radiation and Nuclear
Physics, University of Bonn, Nussallee 14–16, 53115 Bonn, Germany  3
Interdisciplinary Center for Complex Systems, University of Bonn, Brühler
Straße 7, 53175 Bonn, Germany<EMAIL_ADDRESS>
###### Abstract

We evaluate the significance of a recently proposed bivariate jump-diffusion
model for a data-driven characterization of interactions between complex
dynamical systems. For various coupled and non-coupled jump-diffusion
processes, we find that the inevitably finite sampling interval of time-series
data negatively affects the reconstruction accuracy of higher-order
conditional moments that are required to reconstruct the underlying jump-
diffusion equations. We derive correction terms for conditional moments in
higher orders of the sampling interval and demonstrate their suitability to
strongly enhance the data-driven reconstruction accuracy.

  *

Keywords: nonlinear dynamics, stochastic processes, interactions, higher-order
corrections

##  1 Introduction

The problem of reliably characterizing interactions between complex dynamical
systems pervades many scientific fields. Since real-world systems quite often
impose restrictions to standard approaches, linear and nonlinear time-series-
analysis techniques have been developed that allow one to estimate the
strength, the direction, and the functional form of an interaction from pairs
of time series of appropriate system observables. Given that interactions can
manifest themselves in various aspects of the dynamics, analysis techniques
have been developed in diverse fields such as statistics, synchronization
theory, nonlinear dynamics, information theory, and statistical physics (for
an overview, see [1, 2, 3, 4, 5, 6]). Most of these techniques specifically
concentrate on the (low-dimensional) deterministic part of the dynamics and
have shown interactions to impact on amplitudes, phases, frequencies – or even
combinations thereof – as well as on trajectories in the respective phase
spaces.

Many natural systems, however, exhibit both deterministic and stochastic
features, even if the underlying dynamic is deterministic. Stochastic features
in a system’s dynamics may arise from a high number of degrees of freedom,
from random forcing, and/or from nonlinear couplings. If the systems’ dynamics
can be described sufficiently by the Langevin equation, the first- and second-
order Kramers-Moyal (KM) coefficients estimated from time-series data [7, 8]
can serve as an indicator for interactions between stochastic processes [9,
10, 11]. A more general ansatz that accounts for the presence of additive as
well as multiplicative diffusive fluctuations and discontinuous jump
contributions in the time-series data [12, 13, 14] has been proposed only
recently [15], namely a bivariate jump-diffusion model. It consists of two-
dimensional diffusion and two-dimensional jumps that can be coupled to one
another. Such a model could improve theoretical modeling of time-series data
e.g., in the neurosciences [16, 17], condensed matter physics [18, 19, 20],
ecology [21, 22], or finance [23].

When analyzing empirical time-series data, one is faced with the issue of an
inevitably finite sampling interval $\mathop{}\\!\mathrm{d}\hskip 1.0ptt$,
which not only influences the first- and second-order KM coefficients [24] but
also causes non-vanishing higher-order ($>2$) ones. For (one-dimensional)
jump-diffusion processes, additional influences need to be taken into account
[13]: jump events induce terms of order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ in the conditional moments
of even orders and the jump rate and amplitude induce terms of order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ in all conditional
moments. We here extend these studies and investigate the data-driven
reconstruction of stochastic dynamical equations underlying interacting jump-
diffusion processes with finite sampling interval. We will show that in these
cases, corrections for higher orders of the sampling interval strongly enhance
reconstruction accuracy.

The outline of this paper is as follows. In section 2, we recall the
definition of a bivariate jump-diffusion model, and we define our scale-
independent measure to assess the accuracy of the reconstruction of
conditional moments from time-series data. In section 3, we first illustrate
the reconstruction of conditional moments of various jump-diffusion models
with and without couplings from time-series data with finite sampling
interval, thereby emphasizing the necessity for corrections for higher orders
of the sampling interval. We then derive these corrections and demonstrate
their suitability to enhance reconstruction accuracy. Finally, in section 4 we
draw our conclusions.

##  2 Methods

###  2.1 Bivariate jump-diffusion model

A bivariate jump-diffusion process consists of two-dimensional diffusion and
two-dimensional jumps, that can be coupled to one another. It can be modeled
via [12, 8, 15]

$\displaystyle\qquad\quad\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptx_{1}(t)\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptx_{2}(t)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}h_{1}\\\
h_{2}\end{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptt+\begin{pmatrix}g_{1,1}&g_{1,2}\\\
g_{2,1}&g_{2,2}\end{pmatrix}\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptW_{1}\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptW_{2}\end{pmatrix}$
$\displaystyle+\begin{pmatrix}\xi_{1,1}&\xi_{1,2}\\\
\xi_{2,1}&\xi_{2,2}\end{pmatrix}\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptJ_{1}\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptJ_{2}\end{pmatrix}.$ (1)

The drift is a two-dimensional vector $\boldsymbol{h}=(h_{1},h_{2})$ with
$\boldsymbol{h}\in\mathbb{R}^{2}$, where each dimension of $\boldsymbol{h}$,
i.e., $h_{i}$, may depend on state variables $x_{1}(t)$ and $x_{2}(t)$. The
diffusion takes a matrix ${\boldsymbol{g}\in\mathbb{R}^{2\times 2}}$ with the
diagonal elements of $\boldsymbol{g}$ comprise the diffusion coefficients of
self-contained stochastic diffusive processes. The off-diagonal elements
represent interdependencies between the two Wiener processes
$\boldsymbol{W}=(W_{1},W_{2})$, i.e., they result from an interaction between
the two processes. The Wiener processes act as independent Brownian noises for
the state variables with
$\mathopen{}\mathclose{{}\left\langle\mathop{}\\!\mathrm{d}\hskip
1.0ptW_{i}}\right\rangle=0,\mathopen{}\mathclose{{}\left\langle\mathop{}\\!\mathrm{d}\hskip
1.0ptW_{i}^{2}}\right\rangle=\mathop{}\\!\mathrm{d}\hskip 1.0ptt,\forall i$.
The discontinuous jump terms are contained in
$\boldsymbol{\xi}\in\mathbb{R}^{2\times 2}$ and
$\boldsymbol{\mathop{}\\!\mathrm{d}\hskip 1.0ptJ}\in\mathbb{N}^{2}$, where
$\boldsymbol{\mathop{}\\!\mathrm{d}\hskip 1.0ptJ}$ 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 (or size)
$s_{i,j}$. We note that elements of vectors $\boldsymbol{h}$ and
$\boldsymbol{\mathop{}\\!\mathrm{d}\hskip 1.0ptJ}$ as well as of matrices
$\boldsymbol{g}$ and $\boldsymbol{\xi}$ may, in general, be state- and time-
dependent; for convenience of notation, we omit these dependencies.

The two-dimensional KM coefficients of orders $(\ell,m)$ of a bivariate
process read

$\displaystyle D^{(\ell,m)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)$ $\displaystyle=\frac{1}{(\ell+m)!}\lim_{\mathop{}\\!\mathrm{d}\hskip
1.0ptt\to 0}\\!\frac{1}{\mathop{}\\!\mathrm{d}\hskip
1.0ptt}K^{(\ell,m)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ (2)

with the conditional moments

$\displaystyle K^{(\ell,m)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=\mathopen{}\mathclose{{}\left\langle[x_{1}(t+\mathop{}\\!\mathrm{d}\hskip
1.0ptt)-x_{1}(t)]^{\ell}[x_{2}(t+\mathop{}\\!\mathrm{d}\hskip
1.0ptt)-x_{2}(t)]^{m}}\right\rangle\Big{|}_{\begin{subarray}{c}x_{1}(t)=x_{1}\\\
x_{2}(t)=x_{2}\end{subarray}}$ (3)

which can be directly estimated from time-series data [15]. They are related
to the elements of the drift vector, the elements of the diffusion matrix, and
the jump components via:

$\displaystyle K^{(1,0)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)$ $\displaystyle=h_{1}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$ (4)
$\displaystyle K^{(0,1)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)$ $\displaystyle=h_{2}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle K^{(1,1)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)$
$\displaystyle=\Big{[}g_{11}g_{21}+g_{12}g_{22}\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt$ $\displaystyle K^{(2,0)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)$
$\displaystyle=\Big{[}g_{11}^{2}+s_{11}\lambda_{1}+g_{12}^{2}+s_{12}\lambda_{2}\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt$ $\displaystyle K^{(0,2)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)$
$\displaystyle=\Big{[}g_{21}^{2}+s_{21}\lambda_{1}+g_{22}^{2}+s_{22}\lambda_{2}\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt$ $\displaystyle
K^{(2\ell,2m)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$
$\displaystyle=\Big{[}s_{11}^{\ell}s_{21}^{m}\lambda_{1}+s_{12}^{\ell}s_{22}^{m}\lambda_{2}\Big{]}\frac{(2\ell)!}{2^{\ell}\ell!}\frac{(2m)!}{2^{m}m!}\mathop{}\\!\mathrm{d}\hskip
1.0ptt,$

with $(\ell,m)\in\mathbb{N}^{+}$. We omitted the state- and time-dependencies
in the drift and diffusion functions and jump components to enhance
readability. If such processes are coupled, the elements of the drift vector,
the elements of the diffusion matrix, and the jump components may depend on
both state variables and thus information about the coupling may be contained
in the respective conditional moments. However, apart from the drift
functions, there is – to our knowledge – not yet an analytical way to estimate
all diffusion functions and jump components in order to characterize
interactions in the diffusion and jump part of the dynamics. As an
alternative, one can use some optimization techniques to estimate these
unknown functions [8].

###  2.2 Assessing the accuracy of a data-driven reconstruction of conditional
moments

Following Refs. [9, 15], we employ a distance measure to relate theoretical
and numerical results and to quantify the deviation of the obtained
conditional moments from the functions employed. In order to allow a
comparison of the accuracy of the reconstruction of conditional moments of
order $(\ell,m)$, we here use a scale-independent distance measure,
$\mathcal{U}^{(\ell,m)}$, which is based on the bounded relative error of the
difference between estimated and theoretical conditional moments (see A for
details). $\mathcal{U}^{(\ell,m)}<1$ indicates a sufficient accuracy of the
reconstruction of the conditional moment of order $(\ell,m)$.

We estimate conditional moments up to orders $\ell=m=6$ from normalized time
series (zero mean and unit variance) that we obtain from numerically
integrating bivariate jump-diffusion equations (Euler-Maruyama scheme [25]
with a sampling interval $\mathop{}\\!\mathrm{d}\hskip 1.0ptt=10^{-3}$). Time
series consist of $n=10^{7}$ data points (after eliminating $5\times 10^{3}$
transients), and we ensure that individual jump numbers $n_{j}\simeq\lambda n$
varied by at most 10 % for a constant jump rate.

##  3 Results

For our investigations, we consider various interacting jump-diffusion
processes. We begin with reconstructing conditional moments of a bivariate
jump-diffusion model with uni-directional couplings in the drift and in the
diffusion from time-series data. Since interactions between jump-diffusion
processes may lead to over- and underrepresented parts of the dynamics [15],
we then reconstruct conditional moments of a bivariate jump-diffusion model
with a disproportionally weighted drift, diffusion, and jump part from time-
series data.

Figure 1:  a) Reconstruction accuracy ($\mathcal{U}^{(\ell,m)}$) for
conditional moments up to order $\ell=m=6$ of a bivariate jump-diffusion model
with a uni-directional coupling in the diffusion part (see equation (3.1)) for
various values of coupling strength $c_{2}$. $\mathcal{U}^{(\ell,m)}<1$
(horizontal dotted line) indicates a sufficient accuracy. Medians and
interquartile ranges (shaded area) derived from 50 time series generated with
the respective models using random initial conditions. Lines are for eye
guidance only. b) Two-dimensional theoretical conditional moments (equation
(4)) up to order $\ell=m=6$ (black grid) and reconstructed ones (gray surface)
from exemplary time series generated using equation (3.1) with $c_{1}=0$ and
$c_{2}=100$. Reconstruction accuracy is insufficient for conditional moments
of orders $(2,2)$, $(0,4)$, $(4,4)$, $(0,6)$, and $(6,6)$
($\mathcal{U}^{(2,2)}\approx 23$, $\mathcal{U}^{(0,4)}\approx 2\,639$,
$\mathcal{U}^{(4,4)}\approx 9$, $\mathcal{U}^{(0,6)}\approx 682$, and
$\mathcal{U}^{(6,6)}\approx 7$).

###  3.1 Jump-diffusion model with uni-directional couplings

We examine uni-directional couplings in the drift and in the diffusion of an
exemplary bivariate jump-diffusion process modeled via:

$\displaystyle\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip 1.0ptx_{1}(t)\\\
\mathop{}\\!\mathrm{d}\hskip 1.0ptx_{2}(t)\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-x_{1}^{3}+x_{1}\\\
-x_{2}+f_{1}(x_{1})\end{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptt+\begin{pmatrix}0.1&0.5\\\
0.3&0.2+f_{2}(x_{1})\end{pmatrix}\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptW_{1}\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptW_{2}\end{pmatrix}$
$\displaystyle+$ $\displaystyle\begin{pmatrix}\xi_{1,1}&\xi_{1,2}\\\
\xi_{2,1}&\xi_{2,2}\end{pmatrix}\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptJ_{1}\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptJ_{2}\end{pmatrix},$

with

$\displaystyle\boldsymbol{s}=\begin{pmatrix}0.2&0.5\\\
0.3&0.1\end{pmatrix},\,\,\boldsymbol{\lambda}=\begin{pmatrix}0.1\\\
0.3\end{pmatrix}$

and with the coupling terms $f_{1}(x_{1})=c_{1}x_{1}$ and
$f_{2}(x_{1})=c_{2}x_{1}$. We vary the coupling strengths $c_{1}$ and $c_{2}$
each over four orders of magnitude and reconstruct conditional moments from 50
time series with random initial conditions and with finite sampling interval.
For couplings in the drift part ($c_{1}\in[0.01,100]$; $c_{2}=0$), almost all
conditional moments up to order $\ell=m=6$ can be reconstructed with
sufficient accuracy (data not shown). An exception builds the conditional
moment $K^{(0,2)}$ with diffusion contributions of process $x_{2}$. For this
moment, we observe an insufficient reconstruction accuracy
($\mathcal{U}^{(0,2)}>1$) in case of strong couplings ($c_{1}>50$). For
couplings in the diffusion part ($c_{2}\in[0.01,100]$; $c_{1}=0$; see figure
1), we observe at large values of the coupling strength ($c_{2}>1$) rather
strong inaccuracies for some conditional moments. These include conditional
moments with jump contributions of process $x_{2}$, namely $K^{(0,4)}$,
$K^{(0,6)}$, and $K^{(i,i)}$ with $i\in\\{2,4,6\\}$ (figure 1a). A visual
inspection of these conditional moments for a coupling in the diffusion part
with $c_{2}=100$ (figure 1b) indicates that the theoretical moments (equation
(4)) fail to account for a dependency on process $x_{1}$. This dependency on
process $x_{1}$ appears similar to the dependency of $K^{(0,2)}$ on $x_{1}$,
where we note that $K^{(0,2)}$ contains information about the coupling in the
diffusion part (see equation (4)).

###  3.2 Jump-diffusion model with disproportionally weighted parts

Figure 2:  Reconstruction accuracy ($\mathcal{U}^{(\ell,m)}$) of conditional
moments up to order $\ell=m=6$ of a bivariate jump-diffusion model with
disproportionally weighted parts (see equation (6)) for various values of
drift-, diffusion-, and jump-scaling parameters ($\alpha$, $\beta$, $\gamma$).
The vertical dotted line at $\beta=g_{2,1}=0.3$ and $\gamma=0.3$ indicates the
point where the diffusion-scaling parameter is equal to the jump-scaling
parameter. The horizontal dotted line indicates a sufficient accuracy
($\mathcal{U}^{(\ell,m)}<1$). Medians and interquartile ranges (shaded area)
derived from 50 time series generated with the respective models using random
initial conditions. Lines are for eye guidance only.

Next, we examine processes derived from an exemplary bivariate jump-diffusion
model with disproportionally weighted parts, which we obtain by rescaling
drift, diffusion, and jump dynamics. For the latter two, we allow for a mixing
of the Wiener processes and a mixing of the Poisson processes (non-vanishing
off-diagonal elements in diffusion and jump size matrix):

$\displaystyle\qquad\quad\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptx_{1}(t)\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptx_{2}(t)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}-x_{1}^{3}+x_{1}\\\ -\alpha
x_{2}\end{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptt+\begin{pmatrix}0.1&0.5\\\
\beta&0.2\end{pmatrix}\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptW_{1}\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptW_{2}\end{pmatrix}$
$\displaystyle+\begin{pmatrix}\xi_{1,1}&\xi_{1,2}\\\
\xi_{2,1}&\xi_{2,2}\end{pmatrix}\begin{pmatrix}\mathop{}\\!\mathrm{d}\hskip
1.0ptJ_{1}\\\ \mathop{}\\!\mathrm{d}\hskip 1.0ptJ_{2}\end{pmatrix},$ (6)

with

$\displaystyle\boldsymbol{s}=\begin{pmatrix}0.2&0.5\\\
\gamma&0.1\end{pmatrix},\,\,\boldsymbol{\lambda}=\begin{pmatrix}0.1\\\
0.3\end{pmatrix}.$

Here $\alpha$ denotes the drift-, $\beta$ the diffusion- and $\gamma$ the
jump-scaling parameter. We rescale the jump size since it equals the variance
of the Gaussian-distributed jump amplitudes and thus weighs the jump part of
the dynamics. We vary each scaling parameter over four orders of magnitude
(while keeping the others fixed) and reconstruct conditional moments of the
model from 50 time series with random initial conditions and with finite
sampling interval.

When rescaling the drift part (with $\beta=0.3$ and $\gamma=0.3$), almost all
conditional moments up to order $\ell=m=6$ can be reconstructed with
sufficient accuracy (figure 2). The inaccuracy seen for $K^{(0,1)}$ for
$\alpha<1$ is to be expected given our chosen bivariate jump-diffusion model.
If we rescale the diffusion part (with $\alpha=1$ and $\gamma=0.3$), we
observe at large values of the scaling parameter ($\beta>1$) rather strong
inaccuracies for conditional moments $K^{(0,4)}$, $K^{(0,6)}$, and $K^{(i,i)}$
with $i\in\\{2,4,6\\}$ (figure 2). As with jump-diffusion models with uni-
directional couplings (section 3.1), these conditional moments contain jump
contributions of process $x_{2}$. Rescaling the jump part (with $\alpha=1$ and
$\beta=0.3$) has no effect for the considered range of values here
($\gamma\in[0.01,100]$; see figure 2).

###  3.3 Corrections for higher orders of the sampling interval

Our findings presented above demonstrate that uni-directional couplings in the
diffusion may have a similar impact on the accuracy of the reconstruction of
conditional moments as rescaling the diffusion part. Estimated conditional
moments differ from the respective theoretical moments with jump contributions
of process $x_{2}$ at large values of either the coupling strength or the
scaling parameter. Given these observations and since we know that a finite
sampling interval may have a non-negligible impact on the reconstruction of
higher-order conditional moments of one-dimensional jump-diffusion models [13]
from time-series data, we conjecture that a similar impact can be expected for
bivariate jump-diffusion models.

To test this conjecture, we derive the theoretical conditional moments for
different orders of the sampling interval $\mathop{}\\!\mathrm{d}\hskip
1.0ptt$ of a bivariate jump-diffusion model (the derivation and further
expressions of conditional moments can be found in B). With the abbreviations

$\displaystyle A^{(1,0)}$ $\displaystyle=h_{1}$ $\displaystyle A^{(0,1)}$
$\displaystyle=h_{2}$ $\displaystyle B^{(1,1)}$
$\displaystyle=\frac{1}{2}\Big{[}g_{11}g_{21}+g_{12}g_{22}\Big{]}$
$\displaystyle B^{(2,0)}$
$\displaystyle=\frac{1}{2}\Big{[}g_{11}^{2}+s_{11}\lambda_{1}+g_{12}^{2}+s_{12}\lambda_{2}\Big{]}$
$\displaystyle B^{(0,2)}$
$\displaystyle=\frac{1}{2}\Big{[}g_{21}^{2}+s_{21}\lambda_{1}+g_{22}^{2}+s_{22}\lambda_{2}\Big{]}$
$\displaystyle C^{(2\ell,2m)}$
$\displaystyle=\frac{1}{(2\ell+2m)!}\Big{[}s_{11}^{\ell}s_{21}^{m}\lambda_{1}+s_{12}^{\ell}s_{22}^{m}\lambda_{2}\Big{]}\frac{(2\ell)!}{2^{\ell}\ell!}\frac{(2m)!}{2^{m}m!},$

where $(\ell,m)\in\mathbb{N}^{+}$, the theoretical conditional moments of
orders $(0,4)$ and $(2,2)$ with correction terms up to order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ read

$\displaystyle K^{(0,4)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 4!C^{(0,4)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}4!\big{(}B^{(0,2)}\big{)}^{2}+4!\big{(}A^{(0,1)}\partial_{x_{2}}C^{(0,4)}+A^{(1,0)}\partial_{x_{1}}C^{(0,4)}\big{)}$
$\displaystyle+4\cdot
4!C^{(0,4)}\partial_{x_{2}}A^{(0,1)}+4!\big{(}B^{(0,2)}\partial_{x_{2}}^{2}C^{(0,4)}+B^{(2,0)}\partial_{x_{1}}^{2}C^{(0,4)}$
$\displaystyle+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(0,4)}\big{)}+6\cdot
4!\big{(}C^{(0,4)}\partial_{x_{2}}^{2}B^{(0,2)}+C^{(2,2)}\partial_{x_{1}}^{2}B^{(0,2)}\big{)}$
$\displaystyle+4\cdot
5!\big{(}C^{(0,6)}\partial_{x_{2}}^{3}A^{(0,1)}+3C^{(2,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}A^{(0,1)}\big{)}$
$\displaystyle+6\cdot
4!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(0,4)}+3\cdot
6!C^{(2,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(0,2)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$

and

$\displaystyle K^{(2,2)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 4!C^{(2,2)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}\frac{4!}{3}\big{(}B^{(2,0)}B^{(0,2)}+2\big{(}B^{(1,1)}\big{)}^{2}\big{)}$
$\displaystyle+4!\big{(}A^{(1,0)}\partial_{x_{1}}C^{(2,2)}+A^{(0,1)}\partial_{x_{2}}C^{(2,2)}\big{)}$
$\displaystyle+2\cdot
4!\big{(}C^{(2,2)}\partial_{x_{1}}A^{(1,0)}+C^{(2,2)}\partial_{x_{2}}A^{(0,1)}\big{)}$
$\displaystyle+4!\big{(}B^{(2,0)}\partial_{x_{1}}^{2}C^{(2,2)}+B^{(0,2)}\partial_{x_{2}}^{2}C^{(2,2)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(2,2)}\big{)}$
$\displaystyle+4!\big{(}C^{(4,0)}\partial_{x_{1}}^{2}B^{(0,2)}+C^{(2,2)}\partial_{x_{2}}^{2}B^{(0,2)}+C^{(2,2)}\partial_{x_{1}}^{2}B^{(2,0)}+C^{(0,4)}\partial_{x_{2}}^{2}B^{(2,0)}$
$\displaystyle+8C^{(2,2)}\partial_{x_{1}}\partial_{x_{2}}B^{(1,1)}\big{)}+2\cdot\frac{6!}{3!}\big{(}C^{(4,2)}\partial_{x_{1}}^{3}A^{(1,0)}+3C^{(2,4)}\partial_{x_{1}}\partial_{x_{2}}^{2}A^{(1,0)}$
$\displaystyle+C^{(2,4)}\partial_{x_{2}}^{3}A^{(0,1)}+3C^{(4,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}A^{(0,1)}\big{)}+6\cdot
4!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(2,2)}$
$\displaystyle+\frac{5!}{2!}\big{(}6\big{(}C^{(4,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(0,2)}+C^{(2,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(2,0)}\big{)}$
$\displaystyle+16\big{(}C^{(4,2)}\partial_{x_{1}}^{3}\partial_{x_{2}}B^{(1,1)}+C^{(2,4)}\partial_{x_{1}}\partial_{x_{2}}^{3}B^{(1,1)}\big{)}\big{)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3}).$

For the differential operator, we use the short notation
$\partial_{x_{i}}=\frac{\partial}{\partial x_{i}}$. With
$\mathcal{O}(\delta)$, we indicate all terms that contain $C^{(\ell,m)}$ of
higher-order or derivatives $\partial_{x_{i}}^{j}$, $j>3$, and with
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$ all terms that contain
higher orders $(\geq 3)$ of the sampling interval
$\mathop{}\\!\mathrm{d}\hskip 1.0ptt$. We note that the terms of order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ can introduce drift,
diffusion and jump contributions to each conditional moment.

Figure 3: a) Same as figure 1b) but in addition we present theoretical
conditional moments, for which we considered terms of
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ (red grid). In this
case, we obtain a sufficient quality of the reconstruction of conditional
moments of all shown orders except of $(0,6)$, $(4,4)$ and $(6,6)$
($\mathcal{U}_{\text{corr}}^{(0,6)}\approx 37$,
$\mathcal{U}_{\text{corr}}^{(4,4)}\approx 9$ and
$\mathcal{U}_{\text{corr}}^{(6,6)}\approx 7$). b) Same as a) but for a
bivariate jump-diffusion model with a more weighted diffusion part (see
equation (6) with $\alpha=1$, $\beta=100$ and $\gamma=0.3$). By considering
terms of order $\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ in the
theoretical conditional moments, we obtain a sufficient quality of the
reconstruction of conditional moments of all orders except of $(0,6)$, $(4,4)$
and $(6,6)$ ($\mathcal{U}_{\text{corr}}^{(0,6)}\approx 2\,736\,249$,
$\mathcal{U}_{\text{corr}}^{(4,4)}\approx 41$,
$\mathcal{U}_{\text{corr}}^{(6,6)}\approx 8$).

With these correction terms, we now focus on conditional moments, which are
affected by a uni-directional coupling in the diffusion or a more weighted
diffusion part. Already a visual inspection reveals that considering terms of
order $\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ can clearly
improve the reconstruction of some conditional moments (see figure 3;
$\mathcal{U}_{\text{corr}}^{(\ell,m)}$ indicates the distance measure for
which correction terms were considered).

Figure 4:  Reconstruction accuracies with
($\mathcal{U}_{\text{corr}}^{(\ell,m)}$; red) and without
($\mathcal{U}^{(\ell,m)}$; black) integrating correction terms of order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ for various conditional
moments (cf. figure 1a and figure 2). a) Bivariate jump-diffusion model with
uni-directional couplings (see equation (3.1)) for various values of coupling
strength $c_{2}$. b) Bivariate jump-diffusion model (see equation (6)) for
various values of the diffusion-scaling parameter $\beta$. The horizontal
dotted line indicates a sufficient accuracy. Medians and interquartile ranges
(shaded area) derived from 50 time series generated with the respective models
using random initial conditions. Lines are for eye guidance only.

For the data shown in figure 3a, the correction terms of order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ can be of the same or
even greater magnitude than terms of order
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ and thus have a non-
negligible effect on the accuracy of the reconstruction of conditional
moments. Particularly the accuracy of the reconstruction of moments
$K^{(0,4)}$ and $K^{(2,2)}$ is considerably improved
($\mathcal{U}_{\text{corr}}^{(0,4)}<1$ and
$\mathcal{U}_{\text{corr}}^{(2,2)}<1$; see figure 4) even at large values of
the coupling strength or at large values of the diffusion-scaling parameter.
However, we still observe inaccuracies in the reconstruction of conditional
moments $K^{(0,6)}$, $K^{(4,4)}$, and $K^{(6,6)}$, and we expect that
considering terms of higher order of $\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
($\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{i})$, $i\geq 3$) will
further improve the accuracy of the reconstruction of these conditional
moments.

##  4 Concluding remarks

We evaluate the significance of a bivariate jump-diffusion model for a data-
driven characterization of interactions between complex dynamical systems.
Investigating various coupled and non-coupled jump-diffusion processes, we
observed strong deviations between conditional moments of the underlying jump-
diffusion model and those estimated from time-series data and conjectured that
these deviations result from the finiteness of the sampling interval. We
derived correction terms for conditional moments in higher orders of the
sampling interval and could demonstrate that these corrections strongly
enhance the accuracy of a data-driven reconstruction of stochastic evolution
equations from time-series data in terms of bivariate jump-diffusion models.

Our findings demonstrate that the drift, diffusion and jumps induce terms of
order $\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{2})$ in all
conditional moments and are most pronounced in conditional moments with jump
contributions (orders $\geq 4$). A blending of all parts of the dynamics
should thus be taken into account when investigating interacting jump-
diffusion processes. To further enhance the significance of the bivariate
jump-diffusion model for the analysis of empirical data, future studies should
investigate other possible influencing factors such as measurement noise,
limited observation time (finite number of data points), or the impact of
indirect interactions mediated by observed/unobserved additional processes.

We are grateful to M. Reza Rahimi Tabar for constructive discussions and
valuable comments.

## Appendix

##  Appendix A Scale-independent measure to assess the accuracy of a data-
driven reconstruction of conditional moments

For each two-dimensional conditional moments of order $(\ell,m)$, we consider
the unscaled mean bounded relative absolute error [26]

$\displaystyle\mathcal{U}^{(\ell,m)}=\frac{\mathcal{R}_{b}^{(\ell,m)}}{1-\mathcal{R}_{b}^{(\ell,m)}},$

with the weighted average of bounded relative errors

$\displaystyle\mathcal{R}_{b}^{(\ell,m)}=\sum\limits_{i,j=1}^{B}p(x_{1,i},x_{2,j})\frac{\Big{|}\Delta_{ij}^{(\ell,m)}\Big{|}}{\Big{|}\Delta_{ij}^{(\ell,m)}\Big{|}+\Big{|}K^{(\ell,m)}(x_{1,i},x_{2,j},\mathop{}\\!\mathrm{d}\hskip
1.0ptt)\Big{|}},$

where
$\Delta_{ij}^{(\ell,m)}=\hat{K}^{(\ell,m)}(x_{1,i},x_{2,j},\mathop{}\\!\mathrm{d}\hskip
1.0ptt)-K^{(\ell,m)}(x_{1,i},x_{2,j},\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ is
the difference between the estimated and theoretical conditional moment,
$p(\cdot)$ the estimated probability density that is used as normalized
weight, and $B$ the number of bins for each dimension.

$\mathcal{U}^{(\ell,m)}<1$ indicates a sufficient accuracy of the
reconstruction of conditional moments of order $(\ell,m)$ in the sense that
$\Delta_{ij}^{(\ell,m)}$ is smaller than
$K^{(\ell,m)}(x_{1,i},x_{2,j},\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ on
average. We note that the value of $\mathcal{U}^{(\ell,m)}$ becomes large or
undefined if the value of the theoretical conditional moments tends to $0$,
which can lead to a misinterpretation.

In our histogram-based investigations, we used $B=20$ bins for each dimension
and considered a range of $\pm\sigma_{i}$ for each $x_{i}$. We refer to Ref.
[27] for a discussion on the optimal choice of the number of bins and to Refs.
[28, 29] for other, e.g., kernel-based estimation techniques.

##  Appendix B Derivation of conditional moments of bivariate jump-diffusion
models for different orders of $\mathop{}\\!\mathrm{d}\hskip 1.0ptt$

We follow Refs. [13, 8] to derive conditional moments of bivariate jump-
diffusion models for different orders of $\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
using the Kramers-Moyal adjoint operator. With the abbreviations

$\displaystyle A^{(1,0)}$ $\displaystyle=h_{1}$ $\displaystyle A^{(0,1)}$
$\displaystyle=h_{2}$ $\displaystyle B^{(1,1)}$
$\displaystyle=\frac{1}{2}\Big{[}g_{11}g_{21}+g_{12}g_{22}\Big{]}$
$\displaystyle B^{(2,0)}$
$\displaystyle=\frac{1}{2}\Big{[}g_{11}^{2}+s_{11}\lambda_{1}+g_{12}^{2}+s_{12}\lambda_{2}\Big{]}$
$\displaystyle B^{(0,2)}$
$\displaystyle=\frac{1}{2}\Big{[}g_{21}^{2}+s_{21}\lambda_{1}+g_{22}^{2}+s_{22}\lambda_{2}\Big{]}$
$\displaystyle C^{(2\ell,2m)}$
$\displaystyle=\frac{1}{(2\ell+2m)!}\Big{[}s_{11}^{\ell}s_{21}^{m}\lambda_{1}+s_{12}^{\ell}s_{22}^{m}\lambda_{2}\Big{]}\frac{(2\ell)!}{2^{\ell}\ell!}\frac{(2m)!}{2^{m}m!},$

where $(\ell,m)\in\mathbb{N}^{+}$, one can find the following corrections for
conditional moments of orders $\ell=m=6$ of a bivariate jump-diffusion model:

$\displaystyle K^{(1,0)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle A^{(1,0)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}A^{(1,0)}\partial_{x_{1}}A^{(1,0)}+A^{(0,1)}\partial_{x_{2}}A^{(1,0)}$
$\displaystyle+B^{(2,0)}\partial_{x_{1}}^{2}A^{(1,0)}+B^{(0,2)}\partial_{x_{2}}^{2}A^{(1,0)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}A^{(1,0)}$
$\displaystyle+6C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}A^{(1,0)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(2,0)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 2B^{(2,0)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}2\mathopen{}\mathclose{{}\left(A^{(1,0)}}\right)^{2}+2\mathopen{}\mathclose{{}\left(A^{(1,0)}\partial_{x_{1}}B^{(2,0)}+A^{(0,1)}\partial_{x_{2}}B^{(2,0)}}\right)$
$\displaystyle+4\mathopen{}\mathclose{{}\left(B^{(2,0)}\partial_{x_{1}}A^{(1,0)}+B^{(1,1)}\partial_{x_{2}}A^{(1,0)}}\right)$
$\displaystyle+2\mathopen{}\mathclose{{}\left(B^{(2,0)}\partial_{x_{1}}^{2}B^{(2,0)}+B^{(0,2)}\partial_{x_{2}}^{2}B^{(2,0)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}B^{(2,0)}}\right)$
$\displaystyle+8\big{(}C^{(4,0)}\partial_{x_{1}}^{3}A^{(1,0)}+3C^{(2,2)}\partial_{x_{1}}\partial_{x_{2}}^{2}A^{(1,0)}\big{)}$
$\displaystyle+12C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(2,0)}+5!C^{(4,2)}\partial_{x_{1}}^{3}\partial_{x_{2}}^{2}A^{(1,0)}$
$\displaystyle+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(4,0)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 4!C^{(4,0)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}4!\big{(}B^{(2,0)}\big{)}^{2}+4!\big{(}A^{(1,0)}\partial_{x_{1}}C^{(4,0)}+A^{(0,1)}\partial_{x_{2}}C^{(4,0)}\big{)}$
$\displaystyle+4\cdot
4!C^{(4,0)}\partial_{x_{1}}A^{(1,0)}+4!\big{(}B^{(2,0)}\partial_{x_{1}}^{2}C^{(4,0)}+B^{(0,2)}\partial_{x_{2}}^{2}C^{(4,0)}$
$\displaystyle+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(4,0)}\big{)}+6\cdot
4!\big{(}C^{(4,0)}\partial_{x_{1}}^{2}B^{(2,0)}+C^{(2,2)}\partial_{x_{2}}^{2}B^{(2,0)}\big{)}$
$\displaystyle+4\cdot
5!\big{(}C^{(6,0)}\partial_{x_{1}}^{3}A^{(1,0)}+3C^{(4,2)}\partial_{x_{1}}\partial_{x_{2}}^{2}A^{(1,0)}\big{)}$
$\displaystyle+6\cdot
4!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,0)}+3\cdot
6!C^{(4,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(2,0)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(6,0)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 6!C^{(6,0)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}2\cdot
6!B^{(2,0)}C^{(4,0)}+6!\big{(}A^{(1,0)}\partial_{x_{1}}C^{(6,0)}+A^{(0,1)}\partial_{x_{2}}C^{(6,0)}\big{)}$
$\displaystyle+6\cdot
6!C^{(6,0)}\partial_{x_{1}}A^{(1,0)}+6!\big{(}B^{(2,0)}\partial_{x_{1}}^{2}C^{(6,0)}+B^{(0,2)}\partial_{x_{2}}^{2}C^{(6,0)}$
$\displaystyle+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(6,0)}\big{)}+6\cdot
6!\big{(}C^{(4,0)}\partial_{x_{1}}^{2}C^{(4,0)}+C^{(2,2)}\partial_{x_{2}}^{2}C^{(4,0)}\big{)}$
$\displaystyle+\frac{(6!)^{2}}{2!4!}\big{(}C^{(6,0)}\partial_{x_{1}}^{2}B^{(2,0)}+C^{(4,2)}\partial_{x_{2}}^{2}B^{(2,0)}\big{)}$
$\displaystyle+6\cdot\frac{8!}{3!}\big{(}C^{(8,0)}\partial_{x_{1}}^{3}A^{(1,0)}+3C^{(6,2)}\partial_{x_{1}}\partial_{x_{2}}^{2}A^{(1,0)}\big{)}$
$\displaystyle+6\cdot
6!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(6,0)}+6\cdot\frac{(6!)^{2}}{4!2!}C^{(4,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,0)}$
$\displaystyle+30\cdot\frac{8!}{2!2!}C^{(6,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(2,0)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(1,1)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 2B^{(1,1)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}2A^{(1,0)}A^{(0,1)}+2\big{(}A^{(1,0)}\partial_{x_{1}}B^{(1,1)}+A^{(0,1)}\partial_{x_{2}}B^{(1,1)}\big{)}$
$\displaystyle+2\big{(}B^{(2,0)}\partial_{x_{1}}A^{(0,1)}+B^{(1,1)}\partial_{x_{2}}A^{(0,1)}+B^{(1,1)}\partial_{x_{1}}A^{(1,0)}+B^{(0,2)}\partial_{x_{2}}A^{(1,0)}\big{)}$
$\displaystyle+2\big{(}B^{(2,0)}\partial_{x_{1}}^{2}B^{(1,1)}+B^{(0,2)}\partial_{x_{2}}^{2}B^{(1,1)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}B^{(1,1)}\big{)}$
$\displaystyle+4\big{(}C^{(4,0)}\partial_{x_{1}}^{3}A^{(0,1)}+3C^{(2,2)}\partial_{x_{1}}\partial_{x_{2}}^{2}A^{(0,1)}+C^{(0,4)}\partial_{x_{2}}^{3}A^{(1,0)}$
$\displaystyle+3C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}A^{(1,0)}\big{)}+12C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(1,1)}$
$\displaystyle+\frac{5!}{2!}\big{(}C^{(4,2)}\partial_{x_{1}}^{3}\partial_{x_{2}}^{2}A^{(0,1)}+C^{(2,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{3}A^{(1,0)}\big{)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(2,2)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 4!C^{(2,2)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}\frac{4!}{3}\big{(}B^{(2,0)}B^{(0,2)}+2\big{(}B^{(1,1)}\big{)}^{2}\big{)}$
$\displaystyle+4!\big{(}A^{(1,0)}\partial_{x_{1}}C^{(2,2)}+A^{(0,1)}\partial_{x_{2}}C^{(2,2)}\big{)}$
$\displaystyle+2\cdot
4!\big{(}C^{(2,2)}\partial_{x_{1}}A^{(1,0)}+C^{(2,2)}\partial_{x_{2}}A^{(0,1)}\big{)}$
$\displaystyle+4!\big{(}B^{(2,0)}\partial_{x_{1}}^{2}C^{(2,2)}+B^{(0,2)}\partial_{x_{2}}^{2}C^{(2,2)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(2,2)}\big{)}$
$\displaystyle+4!\big{(}C^{(4,0)}\partial_{x_{1}}^{2}B^{(0,2)}+C^{(2,2)}\partial_{x_{2}}^{2}B^{(0,2)}+C^{(2,2)}\partial_{x_{1}}^{2}B^{(2,0)}+C^{(0,4)}\partial_{x_{2}}^{2}B^{(2,0)}$
$\displaystyle+8C^{(2,2)}\partial_{x_{1}}\partial_{x_{2}}B^{(1,1)}\big{)}+2\cdot\frac{6!}{3!}\big{(}C^{(4,2)}\partial_{x_{1}}^{3}A^{(1,0)}+3C^{(2,4)}\partial_{x_{1}}\partial_{x_{2}}^{2}A^{(1,0)}$
$\displaystyle+C^{(2,4)}\partial_{x_{2}}^{3}A^{(0,1)}+3C^{(4,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}A^{(0,1)}\big{)}+6\cdot
4!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(2,2)}$
$\displaystyle+\frac{5!}{2!}\big{(}6\big{(}C^{(4,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(0,2)}+C^{(2,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}B^{(2,0)}\big{)}$
$\displaystyle+16\big{(}C^{(4,2)}\partial_{x_{1}}^{3}\partial_{x_{2}}B^{(1,1)}+C^{(2,4)}\partial_{x_{1}}\partial_{x_{2}}^{3}B^{(1,1)}\big{)}\big{)}+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(4,4)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 8!C^{(4,4)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}12\cdot
6!\big{(}B^{(2,0)}C^{(2,4)}+B^{(0,2)}C^{(4,2)}\big{)}$ $\displaystyle+2\cdot
4!4!\big{(}C^{(4,0)}C^{(0,4)}+18\big{(}C^{(2,2)}\big{)}^{2}\big{)}$
$\displaystyle+12\cdot
6!\big{(}C^{(2,4)}B^{(2,0)}+C^{(4,2)}B^{(0,2)}\big{)}+8!\big{(}A^{(1,0)}\partial_{x_{1}}C^{(4,4)}+A^{(0,1)}\partial_{x_{2}}C^{(4,4)}\big{)}$
$\displaystyle+4\cdot
8!\big{(}C^{(4,4)}\partial_{x_{1}}A^{(1,0)}+C^{(4,4)}\partial_{x_{2}}A^{(0,1)}\big{)}$
$\displaystyle+8!\big{(}B^{(2,0)}\partial_{x_{1}}^{2}C^{(4,4)}+B^{(0,2)}\partial_{x_{2}}^{2}C^{(4,4)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(4,4)}\big{)}$
$\displaystyle+3\cdot
4!6!\big{(}C^{(4,0)}\partial_{x_{1}}^{2}C^{(2,4)}+C^{(2,2)}\partial_{x_{2}}^{2}C^{(2,4)}+C^{(0,4)}\partial_{x_{2}}^{2}C^{(4,2)}$
$\displaystyle+C^{(2,2)}\partial_{x_{1}}^{2}C^{(4,2)}\big{)}$
$\displaystyle+12\cdot
6!\big{(}C^{(6,0)}\partial_{x_{1}}^{2}C^{(0,4)}+C^{(4,2)}\partial_{x_{2}}^{2}C^{(0,4)}+C^{(2,4)}\partial_{x_{1}}^{2}C^{(4,0)}+C^{(0,6)}\partial_{x_{2}}^{2}C^{(4,0)}$
$\displaystyle+36\big{(}C^{(4,2)}\partial_{x_{1}}^{2}C^{(2,2)}+C^{(2,4)}\partial_{x_{2}}^{2}C^{(2,2)}\big{)}\big{)}$
$\displaystyle+6\cdot
8!\big{(}C^{(4,4)}\partial_{x_{1}}^{2}B^{(2,0)}+C^{(2,6)}\partial_{x_{2}}^{2}B^{(2,0)}+C^{(6,2)}\partial_{x_{1}}^{2}B^{(0,2)}+C^{(4,4)}\partial_{x_{2}}^{2}B^{(0,2)}$
$\displaystyle+16\cdot\frac{2!}{3!}C^{(4,4)}\partial_{x_{1}}\partial_{x_{2}}B^{(1,1)}\big{)}+6\cdot
8!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,4)}$
$\displaystyle+6\cdot
8!\big{(}C^{(6,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(0,4)}+C^{(2,6)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,0)}+36C^{(4,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(2,2)}\big{)}$
$\displaystyle+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$
$\displaystyle K^{(6,6)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip
1.0ptt)=$ $\displaystyle 12!C^{(6,6)}\mathop{}\\!\mathrm{d}\hskip 1.0ptt$
$\displaystyle+\frac{1}{2}\Big{[}30\cdot
10!\big{(}B^{(2,0)}C^{(4,6)}+B^{(0,2)}C^{(6,4)}\big{)}$
$\displaystyle+\frac{6!8!}{2!}\big{(}C^{(4,0)}C^{(2,6)}+C^{(0,4)}C^{(6,2)}+15C^{(2,2)}C^{(4,4)}\big{)}$
$\displaystyle+2\cdot
6!6!\big{(}C^{(6,0)}C^{(0,6)}+\frac{(6!)^{2}}{(2!4!)^{2}}C^{(2,4)}C^{(4,2)}\big{)}$
$\displaystyle+\frac{6!8!}{2!}\big{(}C^{(2,6)}C^{(4,0)}+C^{(6,2)}C^{(0,4)}+15C^{(4,4)}C^{(2,2)}\big{)}$
$\displaystyle+30\cdot 10!\big{(}C^{(4,6)}B^{(2,0)}+C^{(6,4)}B^{(0,2)}\big{)}$
$\displaystyle+12!\big{(}A^{(1,0)}\partial_{x_{1}}C^{(6,6)}+A^{(0,1)}\partial_{x_{2}}C^{(6,6)}\big{)}+6\cdot
12!\big{(}C^{(6,6)}\partial_{x_{1}}A^{(1,0)}+C^{(6,6)}\partial_{x_{2}}A^{(0,1)}\big{)}$
$\displaystyle+12!\big{(}B^{(2,0)}\partial_{x_{1}}^{2}C^{(6,6)}+B^{(0,2)}\partial_{x_{2}}^{2}C^{(6,6)}+2B^{(1,1)}\partial_{x_{1}}\partial_{x_{2}}C^{(6,6)}\big{)}$
$\displaystyle+\frac{6!10!}{2!2!}\big{(}C^{(4,0)}\partial_{x_{1}}^{2}C^{(4,6)}+C^{(2,2)}\partial_{x_{2}}^{2}C^{(4,6)}+C^{(2,2)}\partial_{x_{1}}^{2}C^{(6,4)}+C^{(0,4)}\partial_{x_{2}}^{2}C^{(6,4)}\big{)}$
$\displaystyle+30\cdot\frac{6!8!}{2!}\big{(}C^{(6,0)}\partial_{x_{1}}^{2}C^{(2,6)}+C^{(4,2)}\partial_{x_{2}}^{2}C^{(2,6)}+C^{(2,4)}\partial_{x_{1}}^{2}C^{(6,2)}+C^{(0,6)}\partial_{x_{2}}^{2}C^{(6,2)}$
$\displaystyle+15\big{(}C^{(4,2)}\partial_{x_{1}}^{2}C^{(4,4)}+C^{(2,4)}\partial_{x_{2}}^{2}C^{(4,4)}\big{)}\big{)}$
$\displaystyle+\frac{6!8!}{2!}\Big{(}C^{(8,0)}\partial_{x_{1}}^{2}C^{(0,6)}+C^{(6,2)}\partial_{x_{2}}^{2}C^{(0,6)}+C^{(2,6)}\partial_{x_{1}}^{2}C^{(6,0)}+C^{(0,8)}\partial_{x_{2}}^{2}C^{(6,0)}$
$\displaystyle+\frac{(6!)^{2}}{(4!2!)^{2}}\big{(}C^{(4,4)}\partial_{x_{1}}^{2}C^{(4,2)}+C^{(2,6)}\partial_{x_{2}}^{2}C^{(4,2)}+C^{(6,2)}\partial_{x_{1}}^{2}C^{(2,4)}+C^{(4,4)}\partial_{x_{2}}^{2}C^{(2,4)}\big{)}\Big{)}$
$\displaystyle+\frac{6!10!}{2!2!}\big{(}C^{(4,6)}\partial_{x_{1}}^{2}C^{(4,0)}+C^{(2,8)}\partial_{x_{2}}^{2}C^{(4,0)}+C^{(8,2)}\partial_{x_{1}}^{2}C^{(0,4)}+C^{(6,4)}\partial_{x_{2}}^{2}C^{(0,4)}$
$\displaystyle+15\big{(}C^{(6,4)}\partial_{x_{1}}^{2}C^{(2,2)}+C^{(4,6)}\partial_{x_{2}}^{2}C^{(2,2)}\big{)}\big{)}$
$\displaystyle+15\cdot
12!\big{(}C^{(8,4)}\partial_{x_{1}}^{2}B^{(0,2)}+C^{(6,6)}\partial_{x_{2}}^{2}B^{(0,2)}+C^{(6,6)}\partial_{x_{1}}^{2}B^{(2,0)}+C^{(4,8)}\partial_{x_{2}}^{2}B^{(2,0)}$
$\displaystyle+\frac{4!4!}{5!}C^{(6,6)}\partial_{x_{1}}\partial_{x_{2}}B^{(1,1)}\big{)}$
$\displaystyle+6\cdot
12!C^{(2,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(6,6)}$
$\displaystyle+30\cdot
7!8!\big{(}C^{(6,2)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(2,6)}+C^{(2,6)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(6,2)}+15C^{(4,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,4)}\big{)}$
$\displaystyle+\frac{(6!)^{2}}{(4!2)^{2}}\frac{6!10!}{2!2!}\big{(}C^{(4,6)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,2)}+C^{(6,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(2,4)}\big{)}$
$\displaystyle+90\cdot
12!\big{(}C^{(4,8)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(4,0)}+C^{(8,4)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(0,4)}+15C^{(6,6)}\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}C^{(2,2)}\big{)}$
$\displaystyle+\mathcal{O}(\delta)\Big{]}\mathop{}\\!\mathrm{d}\hskip
1.0ptt^{2}+\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$

For the differential operator, we use the short notation
$\partial_{x_{i}}=\frac{\partial}{\partial x_{i}}$. With
$\mathcal{O}(\delta)$, we indicate all terms that contain $C^{(\ell,m)}$ of
higher-order or derivatives $\partial_{x_{i}}^{j}$, $j>3$, and with
$\mathcal{O}(\mathop{}\\!\mathrm{d}\hskip 1.0ptt^{3})$ all terms that contain
higher orders $(\geq 3)$ of the sampling interval
$\mathop{}\\!\mathrm{d}\hskip 1.0ptt$. One obtains
$K^{(m,\ell)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ from
$K^{(\ell,m)}(\boldsymbol{x},t,\mathop{}\\!\mathrm{d}\hskip 1.0ptt)$ by
interchanging the indices $i\in\\{1,2\\}$ of the differential operator
$\partial_{x_{i}}\equiv\frac{\partial}{\partial_{x_{i}}}$ and by interchanging
$\ell$ with $m$ in the orders of $A^{(\ell,m)}$, $B^{(\ell,m)}$, and
$C^{(\ell,m)}$.

## References

## References

  * [1] Pikovsky A S, Rosenblum M G and Kurths J 2001 Synchronization: A universal concept in nonlinear sciences (Cambridge, UK: Cambridge University Press)
  * [2] Kantz H and Schreiber T 2003 Nonlinear Time Series Analysis 2nd ed (Cambridge, UK: Cambridge University Press)
  * [3] Reinsel G C 2003 Elements of multivariate time series analysis 2nd ed (New York: Springer)
  * [4] Hlaváčková-Schindler K, Paluš M, Vejmelka M and Bhattacharya J 2007 Phys. Rep. 441 1–46
  * [5] Marwan N, Romano M C, Thiel M and Kurths J 2007 Phys. Rep. 438 237–329
  * [6] Stankovski T, Pereira T, McClintock P V E and Stefanovska A 2017 Rev. Mod. Phys. 89 045001
  * [7] Friedrich R, Peinke J, Sahimi M and Tabar M R R 2011 Phys. Rep. 506 87–162
  * [8] Tabar M R R 2019 Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems: Using the Methods of Stochastic Processes (Cham-Switzerland: Springer)
  * [9] Prusseit J and Lehnertz K 2008 Phys. Rev. E 77 041914
  * [10] Lehle B 2013 J. Stat. Phys. 152 1145–1169
  * [11] Scholz T, Raischel F, Lopes V V, Lehle B, Wächter M, Peinke J and Lind P G 2017 Phys. Lett. A 381 194–206
  * [12] Anvari M, Tabar M R R, Peinke J and Lehnertz K 2016 Sci. Rep. 6 35435
  * [13] Lehnertz K, Zabawa L and Tabar M R R 2018 New J. Physics 20 113043
  * [14] Hashtroud A M, Mirzahossein E, Zarei F and Tabar M R R 2019 J. Stat. Mech.: Theory Exp. 2019 083213
  * [15] Rydin Gorjão L, Heysel J, Lehnertz K and Tabar M R R 2019 Phys. Rev. E 100 062127
  * [16] Ditlevsen S and Löcherbach E 2017 Stoch. Process. Their Appl. 127 1840–1869
  * [17] Lombardi F, Gómez-Extremera M, Bernaola-Galván P, Vetrivelan R, Saper C B, Scammell T E and Ivanov P C 2020 J. Neurosci 40 171–190
  * [18] Scalliet C, Gnoli A, Puglisi A and Vulpiani A 2015 Phys. Rev. Lett. 114 198001
  * [19] Plati A, Baldassarri A, Gnoli A, Gradenigo G and Puglisi A 2019 Phys. Rev. Lett. 123 038002
  * [20] Plati A and Puglisi A 2020 Phys. Rev. E 102 012908
  * [21] Carpenter S R and Brock W A 2011 Ecology 92 2196–2201
  * [22] Li D, Cui J and Song G 2015 J. Math. Anal. Appl. 430 438–464
  * [23] Aït-Sahalia Y, Cacho-Diaz J and Laeven R J 2015 J. Financ. Econ. 117 585–606
  * [24] see, e.g., M. Ragwitz and H. Kantz, Phys. Rev. Lett. 87, 254501 (2001); R. Friedrich, C. Renner, M. Siefert, and J. Peinke, Phys. Rev. Lett. 89, 217 (2002); P. Sura and J. Barsugli, Phys. Lett. A 305, 304 (2002); D. Kleinhans, R. Friedrich, A. Nawroth, and J. Peinke, Phys. Lett. A 346, 42 (2005); J. Gottschall and J. Peinke, New J. Physics 10, 083034 (2008); S. J. Lade, Phys. Lett. A 373, 3705 (2009); C. Honisch and R. Friedrich, Phys. Rev. E 83, 066701 (2011); C. Honisch, R. Friedrich, F. Hörner, and C. Denz, Phys. Rev. E 86, 026702 (2012); K. Tang, P. Ao, and B. Yuan, EPL (Europhysics Letters) 102, 40003 (2013); A. Vulpiani, and M. Baldovin, J. Stat. Mech.: Theory Exp. 2020, 014003 (2020)
  * [25] Kloeden P E and Platen E 1999 Numerical Solution of Stochastic Differential Equations (Berlin, Heidelberg: Springer)
  * [26] Chen C, Twycross J and Garibaldi J M 2017 PloS one 12 e0174202
  * [27] Knuth K H 2019 Digit. Signal Process. 95 102581
  * [28] Lamouroux D and Lehnertz K 2009 Phys. Lett. A 373 35073512
  * [29] Rydin Gorjão L and Meirinhos F 2019 J. Open Source Softw. 4 1693