{'fraction_non_alphanumeric': 0.04103868013529994, 'fraction_numerical': 0.009698796566879128, 'mean_word_length': 4.32233176168024, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Network inference via process motifs for lagged correlation in linear stochastic processes

A major challenge for causal inference from time-series data is the trade-off between computational feasibility and accuracy. Motivated by process motifs for lagged covariance in an autoregressive model with slow mean-reversion, we propose to infer networks of causal relations via pairwise edge measure (PEMs) that one can easily compute from lagged correlation matrices. Motivated by contributions of process motifs to covariance and lagged variance, we formulate two PEMs that correct for confounding factors and for reverse causation. To demonstrate the performance of our PEMs, we consider network interference from simulations of linear stochastic processes, and we show that our proposed PEMs can infer networks accurately and efficiently. Specifically, for slightly autocorrelated time-series data, our approach achieves accuracies higher than or similar to Granger causality, transfer entropy, and convergent crossmapping -- but with much shorter computation time than possible with any of these methods. Our fast and accurate PEMs are easy-to-implement methods for network inference with a clear theoretical underpinning. They provide promising alternatives to current paradigms for the inference of linear models from time-series data, including Granger causality, vector-autoregression, and sparse inverse covariance estimation.


I. INTRODUCTION
The task of deducing a network's structure from observations of dynamics on the network is a challenging problem with many important applications. For example, many researchers infer gene-interaction networks and/or metabolic networks from data from biological experiments to visualize and understand complex biological systems [1]. The evaluation of clinical trials and marketing tests often involves the inference of causal relationships among a set of observed variables, which one can interpret as a network of causal relationships [2,3]. Neuroscientists infer structural and functional connections among neurons or brain regions from neural recordings to inform their understanding of the brain [4]. Financial networks inferred from recordings of stock and/or currency prices or market transactions are important tools for assessing national and global risks of financial crises [5,6].
In these and other domains, researchers increasingly use the structural analysis of inferred networks to draw conclusions with wide-reaching impact, including implications for design of medical advances, corporate strategies, and public policy. These researchers come from different scientific backgrounds, and some may lack access to advanced computing resources. We are therefore motivated to develop network inference methods (1) that reconstruct networks from observational data with high * Corresponding author: alice.c.schwarze@dartmouth.edu accuracy, (2) whose assumptions and limitations are easy to understand, and (3) that are easy to implement and use, even without access to high-performance computing systems.
The relationship between the network structure and the correlations in a coupled stochastic dynamical system has been the subject of numerous studies across many disciplines [25,27,29,30]. Inverse covariance estimation is a popular network-inference strategy that assumes that a time-series data set is a set of observations of a linear stochastic process [29]. The method aims to identify a arXiv:2208.08871v2 [stat.ML] 22 Aug 2022 plausible coupling structure via inverting the data's sample covariance matrix. Sparse inverse covariance estimation includes additional steps to ensure that the resulting inferred networks are sparse. Graphical LASSO is an example of a widely-used algorithm for performing sparse inverse covariance estimation [30]. For possibly nonlinear systems with linearizable coupling functions, Lai (2017) proposed to infer a network's structure via products of lagged covariance matrices and their inverses [27]. In our paper, we derive process motifs and process-motif contributions for covariance and lagged covariance in a linear stochastic process, and we use these theoretical results to develop new PEMs. Effects of motifs and other aspects of local network structure on inference results has also been the subject of a recent study on network inference from oscillatory dynamics [25].
We propose to use linear combinations of lagged sample correlation matrices as PEMs to infer the structure of unweighted, directed networks with high accuracy at low computational cost from longitudinal observations of a linear stochastic process when the sampling rate of observations is equal to or larger than the characteristic timescale τ of the stochastic process. Similar PEMs (e.g., Granger causality and transfer entropy) are associated with much higher computational cost than our proposed approach, and tend to yield inaccurate results when large sampling rates and/or small characteristic timescales lead to time-series data sets with strong correlations between consecutive observations. In particular, we propose two new PEMs: lagged correlation with a correction for confounding factors (LCCF) and lagged correlation with a correction for reverse causation (LCRC). One can compute these PEMs very fast and easily; an implementation is available as part of our code repository [31].
Our paper proceeds as follows. In Section II, we give a brief introduction to PEM-based network inference. In Section III, we compare two models for linear stochastic processes -the Ornstein-Uhlenbeck process (OUP) and a vector-autoregressive (VAR) model -and explain their connection to process motifs. We introduce our stochastic delay-difference (SDD) model and derive its process motifs for covariance and lagged covariance in Section IV. In Section V, we use the contributions of these process motifs to develop our proposed PEMs, LCCF and LCRC. We demonstrate in Section VI that, for time-series data that we generate using the SDD model, one or both of our PEMs lead to similarly accurate or more accurate inferred networks than several widely-used PEMs including Granger causality, transfer entropy, and convergent crossmapping. We discuss limitations and possible future directions in Section VII.  Figure 1. Inferring networks via pairwise edge measures (PEMs). Our starting point is a directed, unweighted network with unknown structure and observable nodes. Observing node states over time leads to a time-series data set (xt)t=1,2,3,... from which one computes a matrix of PEM values (i.e., "PEM matrix"). For our proposed inference method, the PEM function is a linear combination of lag-k covariance matrices S (k) with coefficient α ∈ R. Thresholding the PEM matrix leads to the adjacency matrix of the inferred network.

II. NETWORK INFERENCE VIA PAIRWISE EDGE MEASURES
In this section, we define the network-inference problem that we consider in this paper. We show a schematic pipeline for PEM-based network inference in Fig. 1. Here our starting point is a network with known nodes and unknown structure (i.e., unknown edges), and we assume that this network is static (as opposed to changing with time), directed, and unweighted. Further, we consider a scenario in which we are able to observe all nodes in the network.
The first step of network inference via PEMs is to collect observations of all node states over a period of time, yielding a time-series data set. Here we assume that one can model the dynamics that determine the node states by a stable linear stochastic system (e.g., a stable Ornstein-Uhlenbeck process or a stable vectorautoregressive model; see Section III) and that the sampling rate 1/∆t at which one observes the node states is constant and smaller than or equal to the inverse characteristic time τ [32] of the linear stochastic system. The next step is to compute PEM values for each node pair from the time-series data. This step yields a "PEM matrix", which is a real-valued square matrix with the same number of rows and columns as the adjacency matrix of the network to be inferred. We reason that node pairs with large PEM values are more likely to be connected by an edge than node pairs with low PEM values, so that the PEM values yield a ranking of node pairs by their likelihood of being connected by an edge.
Either from a priori knowledge of the number of edges in a network or via a statistical test, one determines a threshold for the PEM values. The final step is to apply that threshold to the PEM matrix. One interprets the thresholded PEM matrix as the adjacency matrix of the inferred network.
In our paper, we assume that one knows the number of edges in the network a priori, and we use that knowledge to determine the threshold for the PEM values. This approach obviates the need to choose one or several statistical tests to determine PEM thresholds, but such tests can be used if prior knowledge is unavailable. It thus allows us to compare the ability of different PEMs to rank node pairs correctly by their edge likelihood regardless of the choice of any statistical test.

III. LINEAR STOCHASTIC SYSTEMS AND PROCESS MOTIFS
In this section, we give brief introductions to a continuous-time and a discrete-time model for linear stochastic dynamics on networks. These models are two limiting cases of the dynamical system that is the theoretical foundation for our proposed PEMs. We also give a brief introduction to process motifs, which we use to derive PEMs (see Section V).
In Subsection III A, we briefly review the multivariate OUP, which is a model for continuous-time linear stochastic dynamics on networks. In Subsection III B, we briefly review VAR models, which model discrete-time linear stochastic dynamics on networks. We compare the OUP with VAR models in Subsection III C. Finally, in Subsection III D, we give a brief introduction to process motifs for linear stochastic processes.
A. The multivariate Ornstein-Uhlenbeck process (OUP) The OUP is a linear stochastic differential equation that is commonly used to model the transmission of signals on networks. Applications of the OUP include the coupled dynamics of neurons [33], stock prices [34], and genes [35]. In these studies, a neuron, stock, or gene corresponds to a node i in a network and its state (e.g., a neuron's action potential, a stock's price, or a gene's expression level) is described by a variable x . . of all nodes in a state vector x t , the OUP takes the form where A is the adjacency matrix of a network of coupled variables and w t is a vector-valued Wiener process that has an L 2 -norm with an expectation dw 2 t = dt. The parameters of the OUP include its characteristic time τ , its coupling strength , and its intrinsic noise strength σ. In Fig. 2 (a), we show an example of an OUP with τ = 1, = 0.9, and σ = 0.1 on a unidirectional ring with three nodes.
The existence of steady-state distributions of the stochastic differential equation in Eq. (1) depends on the eigenvalues of A and the coupling parameter . If all eigenvalues of A − I are negative, the OUP has a single stable steady state. At that steady state, x has a multivariate Gaussian distribution with steady-state covariance matrix S, which solves the continuous-time Lyapunov equation [33,36] One can solve Eq. (2) for S and represent its solution as a sum of powers of the adjacency matrix and its transpose [33,36], where x (c) SDD 0.2 (1)) on a unidirectional ring with three nodes. In (b), we show a realization of an order-1 vector-autoregressive model (VAR(1), see Eq. (5)) on the same network. In (c) and (d), we show realizations of the stochastic delay-difference model (SDD, see Eq. (12)) for two different values of the sampling period ∆t. For small ∆t (panel (c)), the SDD model is similar to the OUP. For large ∆t, the SDD model is similar to the VAR(1) model. We show the respective autocorrelation functions for a single node in panels (e) and (f). For all models, we used coupling strength = 0.9 and noise strength σ = 0.1.
An example of a VAR model is the discrete-time stochastic process of the form where the matrices A (k) describe the influence of the states x t−k on x t . To facilitate a direct comparison with the OUP, we include a coupling-strength parameter in Eq. (5). It is more common to use VAR models to fit time-series data, rather than to simulate synthetic data. When using VAR models to fit time-series data, the error term, e t , is the residual of the fit. When one can assume that e t has a Gaussian distribution, one can use the sample covariance matrix e t e T t of the stochastic term to measure a VAR model's goodness of fit [47].
A common way to measure the unconditional linear multivariate Granger causality f (GC) y→x from a variable y to a variable x is to fit (1) the VAR model that regresses over past values of x and (2) the VAR model that regresses over past values of x and y to observed data. The unconditional linear multivariate Granger causality is given by the difference of log-determinants of e t e T t for the two models [23] . When using the model in Eq. (5), one assumes that only the p most recent measurements have a nonnegligible influence on x t . The integer p is called the order of the VAR model, and it is common to denote a VAR model of order p by VAR(p). For example, the VAR(1) model from Eq. (6) can have the form (For the remainder of our paper, we refer to the model in Eq. (6) as "the VAR(1) model".) In Fig. 2 (b), we show an example realization of the VAR(1) model for which A (1) is the adjacency matrix of a unidirectional ring with 3 nodes. If the error e t follows a multivariate Gaussian distribution with covariance matrix e t e T t = τ σ 2 n I and all eigenvalues λ i of A (1) have negative real parts and absolute values |λ i | < 1, the VAR(1) model in Eq. (6) has a stable steady-state distribution that is Gaussian. This stable distribution has an associated steady-state covariance matrix S (0) that satisfies the discrete-time Lyapunov equation [33] One can write the solution of Eq. (8) as a sum of powers of the adjacency matrix and its transpose, In Eq. (8) and Eq. (9), we denote the steady-state covariance matrix of a discrete-time process by S (0) instead of S (which is our notation for covariance matrices of continuous-time processes). Denoting covariance matrices of a discrete-time processes by S (0) makes it possible to denote lag-k covariance matrices of discrete-time processes in Subsection IV C and thereafter by S (k) .

C. Autocorrelation in the OUP and VAR models
One can view the VAR(1) model as a discrete-time equivalent of an OUP [48]. However, when fitting Eqs. (1) and (5) to data, one can expect that the estimated matrices A for the OUP and A (1) models are structurally different. The OUP with A = 0, has an exponentially decaying autocorrelation function, In contrast, the VAR(1) model with A (1) = 0 has the autocorrelation function, The autocorrelation function of the OUP with A = 0 indicates that the OUP is a continuous-time limit of a first-order Markov process. In contrast, the autocorrelation function of the VAR(1) model A (1) = 0 indicates that the VAR(1) model describes a process in which the value of a state variable is independent of the variable's past values. (Such processes are colloquially called "zeroth"-order Markov processes [49,50].) The VAR(1) model describes a first-order Markov process only if A (1) = 0.
These differences between the OUP and the VAR(1) model influence the interpretation of diagonal elements of A and A (1) . A non-zero diagonal element A (1) ii in A (1) indicates that the state variable x (i) t is positively or negatively correlated with its most recent past value x (i) t−1 . By contrast, a non-zero diagonal element A ii in A indicates the state variable x (i) t depends more strongly or more weakly on its immediate past x (i) t−dt than one would expect given the exponentially decaying autocorrelation function in Eq. (10). A comparison of panels (a) and (b) in Fig. 2 illustrates the qualitative differences between the OUP and the VAR(1) model with A = A (1) . In panels (e) and (f) of Fig. 2, we use a dark blue curve and a gray curve to indicate the normalized autocorrelation func- for the OUP with A = 0 and the VAR(1) model with A (1) = 0. The normalized autocorrelation function for the OUP gradually decays as a function of t − t. The autocorrelation function for the VAR(1) model is only defined where t − t has integer values. It is equal to 1 for t − t = 0 and zero for all other integer values of t − t.

Matrix
Edge graph Walk graph i j i j i j Figure 3. Representations of process motifs. On the left, we show matrix functions whose elements indicate counts of process motifs for covariance in a network. In the center, we show the examples of small network structures on which the corresponding process motifs can occur. White nodes indicate source nodes, which are the starting points of the walks that form a process motif. On the right, we show the corresponding process motifs as walk graphs, in which each curved blue edge with a label "k" indicates a walk of length k on the network.

D. Process motifs for linear stochastic processes
Many researchers have studied motifs in network structure. A structure motif G is a (small) network that (1) occurs as a subnetwork of a network G ⊇ G and (2) is considered to be important for a function of the network G [51]. By contrast, a process motif is a structured set of walks that (1) occurs on a network G and (2) is considered to be important for a function of G [36].
For linear dynamical systems, one can express many system properties as a sum of contributions of process motifs [36]. An example are the OUP's steady-state covariances S ij , for which we show several representations of their process motifs in Fig. 3. Equation (3) indicates that one can express the steady-state covariance S ij between the nodes i and j as If A is an unweighted adjacency matrix, the matrix element A F (A T ) L− F ij corresponds to the number of occurrences of a process motif (i.e., a structured set of walks). The value of c L− F , F indicates the contribution of each of these process-motif occurrences to S ij .
In Fig. 3, we show the examples of process motifs for steady-state covariance in the OUP. For each example, we show its matrix representation, its representation in a graph of edges, and its representation as a structured set of walks (i.e., walk graphs). If one allows walks to have a 0 length (i.e., one allows the source node v to coincide with i or j), one can consider all process motifs for steady-state covariance in the OUP to be a structured set of two walks that start at the same node v and end at nodes i and j, respectively. A process motif for steadystate covariance is then fully characterized by the lengths of the two walks. We call the walk from v to j the walk in the "forward" direction, and we denote its length by F . We call the walk from v to i the walk in the "backward" direction, and we denote its length by B . The total length of a process motif is L := B + F . We use the pair ( B , F ) to denote the process motif whose number of occurrences between i and j is given by

IV. A STOCHASTIC DELAY-DIFFERENCE MODEL FOR A SAMPLED PROCESS WITH SLOW MEAN-REVERSION
In this section, we introduce a stochastic delaydifference (SDD) model for sampled time-series data, and we derive its process motifs and process-motif contributions for steady-state covariance and lagged steady-state covariance. In Subsection IV A, we introduce the stochastic difference equation that we use to interpolate between the OUP in Eq. (1) and the VAR model in Eq. (5). In Subsections IV B and IV C, we derive and discuss the process motifs for steady-state covariance and lagged steadystate covariance for the stochastic delay-difference equation with Markovian dynamics. In Subsection IV D, we show how one can extend these results to non-Markovian dynamics.

A. Stochastic delay-difference equation
Time-series data consists of measurements or observations at a sequence of discrete time points. Whether consecutive observations in the data are correlated depends on the ratio between the characteristic time τ of the observed dynamics and the timestep ∆t between consecutive observations. As a dynamical model for such data, we consider a stochastic delay-difference (SDD) model, where ∆t τ := ∆t/τ is the sampling period in units of τ and every element of the vector ∆w t is a n-variate 0mean Gaussian random variable with ∆w t 2 = ∆t. In analogy to the VAR model, we call p the order of the SDD model. When ∆t = τ = 1, the SDD model is equivalent to the VAR model in Eq. (5).
The SDD model of order 1 has the form (12) which is formally identical to the OUP in Eq. (1) when one replaces ∆t with a differential dt. We consider the sampling period ∆t to be a parameter of the SDD model, and it is not infinitesimally small. The SDD model with order p > 1 allows edges in a network to have associated transmission lags that are multiples of ∆t. In contrast to the OUP, one can thus use the SDD model for non-Markovian dynamics.
In Fig. 2 (c), we show a realization of the SDD model of order 1 with a small sampling period (i.e., ∆t = 0.2) on a unidirectional ring. In panel (d), we show a realization of the same SDD model with a large sampling period (i.e., ∆t = 0.8). One can observe that the SDD model with small ∆t is visually similar to the realization of the OUP shown in Fig. 2 (a). In Fig. 2 (e), one can observe that the SDD model with small ∆t and the OUP have similar autocorrelation functions. The SDD model with large ∆t is visually similar to a realization of the VAR(1) (see Fig. 2 (b)) and has a similar autocorrelation function to the VAR(1) (see Fig. 2 (f)).

B. Process motifs for covariance
In this section, we derive the process motifs and process-motif contributions for steady-state covariance for the SDD model of order 1. The steady-state covariance matrix of the SDD model satisfies the discrete-time Lyapunov equation where K := (1 − ∆t τ )I + ∆t τ A (1) . One can express the solution of Eq. (13) as The contributions of process motifs to steady-state covariance in the SDD model are where we define and use the shorthand ∧ p,q := max{p, q} + 1 .
We derive Eq. (14) in Appendix A. The function ψ p,q (z) captures the effect of nodes "remembering" their previous states (i.e., the effect of non-zero correlations between x (i) t and x (i) t−∆t ) on the contributions c (0) B , F . This effect depends on the autocorrelation function of the process, which depends on the ratio of the sampling period ∆t and the characteristic time τ (see Section III C). The function ψ p,q (z) thus depends on ∆t τ , and it is independent of the coupling parameter and the noise strength σ.
As ∆t τ approaches 0, the steady-state covariance matrix S (0) of the SDD model converges to the steady-state covariance matrices of the OUP. We derive this limit in Appendix B. When ∆t τ = 1 in Eq. (14), the steadystate covariance matrix of the SDD model is equal to the steady-state covariance matrix of a VAR(1) model with a time-independent error covariance matrix given by Eq. (7).
One can interpret the elements of the matrix functions (A (1) ) ((A (1) ) T ) L− of the adjacency matrix as counts of process-motif occurrences (see Section III D) and the (i, j)-th element of the SDD model's steady-state covariance matrix as a sum of these counts -each weighted by c (0) B , F . Comparing Eqs. (14) and (3), one can see that both steady-state covariance matrices are sums of the same matrix functions of the adjacency matrix. The OUP in Eq. (1) and the SDD model therefore have the same process motifs for covariance, albeit with different contributions -c B , F for the OUP, and c (0) B , F for the SDD model. Comparing Eqs. (14) and (9), one can see that the steady-state covariance matrix of the VAR(1) model is a sum of a subset of the matrix functions that contribute to the steady-state covariance matrix of the SDD model. The process motifs for covariance in the VAR(1) model in Eq. (6) are therefore a subset of the process motifs for covariance in the SDD model.
In Fig. 4, we show the contributions of process-motif occurrences to steady-state covariance for the OUP (which corresponds to the SDD model with ∆t τ → 0), the SDD model with ∆t τ = 0.2, the SDD model with ∆t τ = 0.8, and the VAR(1) model (which corresponds to the SDD model with ∆t τ = 1. One observes that many different process motifs contribute to the steadystate covariances of the OUP. The contributions tend to decrease with increasing difference between the length of the walk to node i and the length of the walk to node j (i.e., | B − F |). The contributions c (0) B , F also tend to decrease with increasing overall length of the process motif (i.e., B + F ). All shown process motifs have larger contributions to the steady-state covariance of the OUP than for the SDD model with a finite sampling rate. As ∆t τ approaches 1, the SDD model corresponds with the VAR(1) model and the only process motifs with positive contributions to steady-state covariance are process motifs for which the walk to node i and the walk to node j have the same length (i.e., B = F ). For the three discrete-time processes that we show in Fig. 4, the process motif (1, 1) has the largest contribution to steady- Figure 4. Contributions of process motifs to covariance in linear stochastic processes. We show the contributions of process motifs to covariance for the Ornstein-Uhlenbeck process (OUP), the stochastic delay-difference model with ∆tτ = 0.2 (SDD0.2), the stochastic delay-difference model with ∆tτ = 0.8 (SDD0.8), and an order-1 vector autoregressive model (VAR (1)). Bars with black outlines indicate the largest contributions to covariance for each stochastic process. We observe that the process motifs (0,1) and (1,0) -which correspond to true causation and reverse causation [52,53], respectively -contribute most to covariance in the OUP. In the VAR(1) model, the process motif (1,1) -which corresponds to existence of a confounding factor -contributes most to covariance. The contributions to covariance in the SDD model interpolate between the contributions to covariance in the OUP and the VAR(1) model. For all models, we used = 0.9, σ = 1.0. state covariance. For the OUP, the process motifs (0, 1) and (1, 0) have the largest contribution to steady-state covariance.
To explain why the process motif (1, 1) is not the process motif with the largest contribution to steady-state covariance in the OUP, we remark that one can interpret each process motif as a way by which a signal on any node in the network can reach the nodes i and j at the same time. On its way to i and j, a signal in the OUP can traverse at most one edge per timestep ∆t. If it does, its amplitude is scaled down by a factor that depends on B and F and is proportional to ∆t τ . For a process in which node states x (i) t are correlated with their past x (i) t for t < t, signals do not need to traverse an edge in each timestep. One can think of signals "lingering" at a node over some number of timesteps. If a signal lingers at a node, its amplitude is scaled down by a factor that depends on B and F and is proportional to (1−∆t τ ). In the process motif (1, 1), a signal does not have to linger at any node to contribute to steady-state covariance. Its contribution c 1,0 . One thus has larger contributions of the process motifs (0, 1) and (1, 0) than the contribution of the motif (1, 1) in the OUP.

C. Process motifs for lagged covariance
In this section, we derive the process motifs and process-motif contributions for lagged steady-state covariance in the SDD model.
At steady state, the lagged covariance functions s ij (t, t ) := x j,t x T i,t of a stochastic process for given nodes i and j depend only on the difference ζ := t −t between time points. We can thus write the matrix-valued lagged-covariance function For a discrete-time process, the function S(ζ) is defined on ζ ∈ {0, ±∆t, ±2∆t, . . . } . For any integer k, we call the matrix the lag-k steady-state covariance matrix of a discrete stochastic process. For positive k, the lag-k steady-state covariance matrix of the SDD model of order 1 is connected to the k -lag steady-state covariance matrices with k ∈ {0, . . . , k − 1} via a recursive relationship, From the recursive relationship in Eq. (15) and the order-1 SDD model's covariance matrix in Eq. (14), it follows that the lag-k steady-state covariance matrix has the form where c (k) From comparing Eqs. (14) and (16), one can see that the same process motifs contribute to all lag-k steadystate covariance matrices for, k = 0, 1, 2, . . . , of the SDD model. Their contributions c (k) B , F to the lag-k steadystate covariance matrix are different for different k. The differences between them depend on the parameters and τ of the SDD model, as well as its sampling period ∆t.
In Fig. 5, we show the contributions of process motifs to the lag-k steady-state covariance matrices with k ∈ {0, 1, 2, 3}. We show results for an order-1 SDD model with ∆t τ = 0.2 and for an order-1 SDD model with ∆t τ = 0.8. One observes that for ∆t τ = 0.8 and k ≥ 1, the largest contributions to elements of the lagk steady-state covariance matrix are associated with klength walks from node j to node i. For ∆t τ = 0.2 and k ≥ 1, the largest contributions are associated with a length-1 walk from j to i.
A comparison of the contributions of process motifs to the lag-1 steady-state covariance matrix (see second column of Fig. 5) for different ∆t τ suggests that using the elements of the lag-1 sample covariance matrix, S (1) , as a PEM tends to yield a higher accuracy for network inference when node states at subsequent timesteps, x (i) t and x (i) t+∆t , are weakly correlated (i.e., ∆t τ is large) than when node states at subsequent timesteps are strongly correlated (i.e., ∆t τ is small). When ∆t τ is large, the contribution of the process motif (0, 1) is much larger than the contributions of other process motifs. It follows that a large lagged covariance from node j to node i either suggests a direct connection from j to i or many indirect connections corresponding to other process motifs to covariance. By contrast, for the SDD model with small ∆t τ , the process motifs (0, 1), (1, 0), and (1, 1) have similar contributions to the lagged covariance from node j to node i. A large lagged covariance from node j to node i can thus result from many network structures, including an edge from j to i, or edges from a third node to both i and j. For the SDD model with order p > 1, the SDD equation (11) includes matrices A (k) with k > 1. The nonzero elements of these matrices indicate edges with lagged signal transmission (i.e., it takes a signal more than one Figure 5. Contributions of process motifs to lag-k covariance. We show results for the order-1 stochastic delaydifference (SDD) model with ∆tτ = 0.2 and ∆tτ = 0.8. Bars with black outlines indicate the largest contributions to lag-k steady-state covariance for each stochastic process. For small ∆tτ , the process motif (0,1) tends to contribute the most to lagged steady-state covariance. For large ∆tτ , the process motifs (0, k) tend to contribute the most to the lag-k steadystate covariance. Results are for the SDD model with = 0.9, σ = 1.0.
timestep to traverse the edge). We define the transmission delay associated with an edge e to be δ e = k e − 1, where k e is the number of timesteps that it takes for a signal to traverse e. There are many different reasons why one can expect to observe lagged signal transmission in time-series data. For example, networks that are embedded in a physical space (e.g., road networks, transport networks, vascular networks, natural neural networks, and so on) can include edges with different lengths. When objects or signals travel at a constant speed, they traverse some edges faster than others in these networks. Unobserved nodes are another possible cause for observing delayed signal transmission. For example, when a network includes a path (j, u, i) through an unobserved node u, signals traversing this path can contribute positively to the lag-2 correlation from node j to node i.
One can model delayed signal transmission in the SDD model as unobserved nodes in a weighted network. Consider a network G in which each edge (j, i) has a transmission delay δ i,j . One can construct another network G with unobservable nodes such that each edge (j, i) in G corresponds to a path (j, u 1 , u 2 , . . . , u δij , i) of length δ ij + 1 from j via δ ij unobservable nodes to i. To model the decay of a signal's amplitude along edges, one can associate each edge in G with a weight ∆t τ . In G , each path that corresponds to an edge in G includes one edge with weight ∆t τ . All other edges on the path have a weight of 1. This choice of edge weights for G ensures that the amplitude of a signal that traverses an edge in G always decreases by a factor ∆t τ , regardless of the transmission delay associated with that edge. In Fig.6, we give an example of a network G and its associated network G with unobserved nodes.
When one constructs G as described above, the lag-k steady-state covariance matrix of the observable nodes in G is equal to the lag-k steady-state covariance matrix of G for all integers k. A process motif ( B , F ) that has transmission delays on edges has a contribution to the lag-k steady-state covariance in the SDD model of order p > 1. In Eq. (17), we useδ B to denote the sum transmission delays of the edges traversed in the "backward" direction (i.e., from the source node v to node i) andδ F to denote the sum transmission delays of the edges traversed in the "forward" direction (i.e., from v to j). The data shown in Fig. 5 illustrates that, for large ∆t τ , the process motif (0, k) has the largest contribution to the elements of the lag-k steady-state covariance matrix. This observation motivates us to use lag-k correlation matrices to infer edges with transmission delays k − 1 in the same way that we use lag-1 correlation matrices to infer edges without transmission delay.

V. FROM PROCESS MOTIFS TO PAIRWISE EDGE MEASURES
Previously, we discussed that the process-motif contributions to lagged steady-state covariances in the SDD model can suggest whether lag-k steady-state covariances are likely to be good PEMs (see Section IV D). In the current section, we use process-motif contributions to derive new PEMs that are linear combinations of lag-k steadystate covariances. In Subsection V A, we discuss the combination of process-motif contributions that would be an Representation of a stochastic delaydifference (SDD) model of order k as the SDD model of order 1. In (a), we show a network G of nodes and edges. Edge labels indicate the transmission delays δe and weights associated with the SDD model on the network. In (b), we show a corresponding network G , in which we replaced each edge e of G by a path of length δe + 1. Gray nodes are nodes that are in G but not in G. We consider them to be "unobservable nodes". We choose edge weights in G such that the weight of each edge in G is equal to the product of edge weights along the corresponding path in G .
ideal PEM for network inference from the SDD model. In Subsection V B, we propose two PEMs for network inference from the SDD model; one PEM that corrects lag-1 steady-state covariance for confounding factors and one PEM that corrects lag-1 steady-state covariance for reverse causation. In Subsection V C, we explain how we compute these PEMs from time-series data.
A. An ideal edge-likelihood score function An ideal PEM would indicate edges with certainty. A hypothetical example of such an ideal score would be as a PEM, one would be able to infer networks from observations of the SDD model of order 1 with an accuracy of 1.
In the left column of Fig. 7, we show the contributions of process motifs to the ideal PEM f (ideal) ij . The only process motif that has a nonzero contribution to f (ideal) ij is the process motif (0, 1). A comparison of the processmotif contributions for f (ideal) ij to the process-motif contribution for S (1) (see Fig. 5) illustrates why lag-k covariance is not an ideal PEM. Many process motifs contribute positively to lag-k covariance. The process motif (0, 1) has the largest contribution to lag-1 covariance, but the motifs (1, 0) and (1, 1) can have almost equally large contributions for small ∆t τ .
One can think of process motifs for covariance as defining a basis of an infinite-dimensional "process-motif contribution space" and use the contributions c (k) B , F of process motifs to a PEM to construct a "contribution vector" for the PEM in that space. We use such contribution vectors to compare PEMs. For example, the contribution vector for the ideal PEM f (ideal) ij is only supported on the process motif (0, 1). For ∆t τ < 1, the contribution vectors for lag-k steady-state covariance of the SDD model of order 1 are supported on all process motifs. As ∆t τ approaches one, the contribution vector for lag-k steadystate covariance localizes on process motifs ( B , F ) with F = B + k. We are not aware of any function that one can compute in practice from a time-series data set and that fulfills our criterion for being an ideal PEM (i.e., it has a contribution vector with support only on the process motif (0, 1)). However, one can combine the contribution vectors of lag-k steady-state covariances to construct PEMs whose contribution vectors localize more strongly on (0, 1) than the contribution vectors of lag-k steady-state covariance for any k = 0, 1, 2, . . . . We propose two such constructions in Subsection V B.

B. Correcting lag-1 covariance for confounding factors and reverse causation
To construct a PEM that is localized on the process motif (0, 1), we use linear combinations of contribution vectors for lag-k steady-state covariances. In principle, one can use any number of lag-k steady-state covariances with any values of k to construct a PEM. However, when using lagged sample covariance matrices to approximate lagged steady-state covariance matrices, the sampling errors of each lagged sample covariance matrix can lead to errors of computed PEM values, especially when k is large.
To minimize these sources of error, we construct PEMs that are linear combinations of few lag-k correlation matrices that have small k. Specifically, we consider candidate PEMs of the form for some non-negative α.

Correction for confounding factors
Most process motifs for covariance and lag-1 covariance in the SDD model include two walks of positive length and thus three distinct nodes; a node i, a node j, and a source node v from which signals can reach both i and j. The state of v is a confounding factor for the node pair (j, i) because v can lead to an increase in covariance and lagged covariance between the states of the nodes i and j even when there is no path from i to j and no path from j to i. We thus refer to the process motifs for covariance and lagged covariance with three distinct nodes as confounder motifs. To construct a PEM with a contribution vector that is localized on the process motif (0, 1), we use the steady-state covariance matrix S (0) to correct for the confounder motif (1, 1), which is the confounder motif that contributes most to lag-1 steady-state covariances. To do so, define the PEM with a correction term [54] α (LCCF) := c The PEM in Eq. (18) combines lag-1 steady-state covariance with a correction term that cancels the contributions of the confounder motif (1, 1) to lag-1 covariance. We call this PEM lagged correlation corrected for confounding factors (LCCF), and we show the contributions of several process motifs in the third column of Fig. 7.

Correction for reverse causation
Some process motifs for lag-1 steady-state covariance from j to i include only two distinct nodes and a single walk from i to j. This means that signal transmission from i to j can lead to an increase of the lag-1 steadystate covariance from j to i and thus to a false-positive detection of an edge from j to i. Some researchers refer to this phenomenon as reverse causation [52,53]. When ∆t τ is small, reverse causation can impact the accuracy of network inference because the difference between the contributions of the process motif (0, 1) and (1, 0) to lag-1 steady-state covariance is small. This small difference between c 1,0 makes it difficult to infer the direction of an edge in a network accurately from the lag-1 steady-state covariance matrix of the SDD model with small ∆t τ . A large lag-1 steady-state covariance from a node j to a node i can be caused by an edge (j, i), but for small ∆t τ , it is almost equally likely that it is caused by an edge in the reverse direction (i.e., an edge (i, j)).
To mitigate false-positive edge detections due to reverse causation, we define the PEM with a correction term The correction term α (LCRC) cancels the contributions of reverse causation to lag-1 steady-state covariance. We call the PEM in Eq. (20) lagged correlation corrected for reverse causation (LCRC), and we show the contributions of several process motifs in the fourth column of Fig. 7.

C. Computing pairwise edge measures from time-series data
In the OUP, the SDD process, and the VAR model, the process motifs that contribute to covariance are a subset of the process motifs that contribute to correlations for k = 0, 1, 2, . . . . In the OUP, the process motifs for covariance tend to be the process motifs that contribute the most to correlation [36]. We thus expect that the correction factors that are needed to correct lag-1 steadystate correlation matrices for confounding factors or reverse causation are very similar to the correction factors that we proposed in Eqs. (19) and (21). Using sample covariance matrices as estimators of steady-state covariance matrices and sample correlation matrices as estimators of steady-state correlation matrices for synthetic time-series data (see Section VI), we found in initial tests that we can infer networks with higher accuracy using correlations than when using covariances. We thus use the following estimators for f (LCCF) and f (LCRC) to infer networks from time-series data; whereR (k) is the k-lag sample correlation matrix of a time-series data set.
To compute the PEM estimators in Eqs. (22) and (23), one needs to know the value of ∆t τ := ∆t/τ . We assume that the sampling period ∆t of the time-series data is known. In Appendix C, We propose a simple procedure for estimating τ −1 from a time-series data set. In Appendix, D, we explain how to generalize Eqs. (22) and (23) when one uses the SDD model of order p > 1 as a model for time-series data.

VI. NETWORK INFERENCE VIA CORRELATION-CORRECTED CROSS-CORRELATION
In this section, we evaluate the accuracy of our proposed PEMs and several commonly used PEMs for network inference from simulated time-series data. We compare the accuracy of inferred networks from different PEMs as well as the computation time associated with different PEMs. For all PEMs that we consider here, the accuracy of inferred networks depends on parameters of the SDD model that we use to generate time-series data and on structural properties of the ground-truth networks. In Subsection VI A, we investigate how the mean accuracy of inferred networks changes as one varies several parameters of the SDD model. These parameters are the coupling strength , the characteristic time τ , and the sampling period ∆t. In Subsection VI B, we present results on the accuracy of inferred networks as one varies structural properties of ground-truth networks. We consider variations in the network size n, the edge density d e , and the edge reciprocity r e . We also compare the mean accuracy of inferred networks when ground-truth networks are realizations of one of several random-graph models. In Subsection VI C, we investigate the effect of system noise σ and measurement noise η on the mean accuracy of inferred networks. Since most of the PEMs formally rely on time-series data being stationary, we also study how nonstationarity of time-series data affects the mean accuracy of inferred networks. We show our results of this study in Subsection VI D. In Subsection VI E, we compare the computation time for different PEMs as one varies the size n of the ground-truth network, the number N of observations, and the maximum edge transmission lagδ.
For generating time-series data, we simulate the SDD model with coupling strength , characteristic time τ , sampling period ∆t, and system noise strength σ on a ground-truth network with n nodes, an edge density d e , and an edge reciprocity r e . Unless we state otherwise, we use the same default parameters for all numerical experiments presented in this section. For the network generation, our default network model is a G(n, m) randomgraph model and with n = 10, d e = 0.5, and r e = 0.5. To ensure convergence of simulations and facilitate a comparison of processes with different coupling strengths, we normalize adjacency matrices so that their largest eigenvalue has absolute value |λ max | = 1. Whenever a sample network from a random-graph model has a nilpotent adjacency matrix, we reject that network and sample a new network from the random-graph model. Our default parameters for simulating the SDD model are = 0.9, τ = 1, ∆t = 0.5, and σ = 0.2 (see [58] for complete list of default simulation parameters).
A. Dependence on parameters of the stochastic delay-difference model In Figure 8, we show the mean inference accuracy Φ for inferring realizations of a G(n, m) random-graph model with n = 10 nodes and m = 45 (i.e., an edge density d e = 0.5) as we vary parameters of the SDD model. In Figure 9, we show differences in mean inference accuracy between several PEMs and our proposed PEM f (LCRC) as one varies and τ .

Sampling period
In Figure 8 (a), we show our results for the mean inference accuracy Φ as one varies the sampling period ∆t of a SDD model with no transmission lags on edges (i.e., the SDD model of order 1). We observe that the mean inference accuracy for most of the considered PEMs improves as ∆t increases. An exception is the PEM for OUI. This PEM relies on the assumption that the timeseries data arises from a continuous-time process. This assumption is approximately correct when ∆t is small (e.g., ∆t 0.3 in Fig. 8 (a)). When ∆t is large, the SDD model is very different from to a continuous-time process. Accordingly, assuming that the time-series data is the result of a continuous-time process leads to lower mean inference accuracy for OUI.
When ∆t = τ , the correction factors for LCCF (see Eq. (19)) and LCRC (see Eq. (21)) have the same value. We thus observe equal mean inference accuracy for LCCF and LCRC at ∆t = 1. The PEM LCRC leads to a higher mean inference accuracy than LCCF at small ∆t (i.e., ∆t 0.2) and to a lower mean inference accuracy than LCCF at larger ∆t (i.e., ∆t 0.3).
Among the PEMs that we compare in Fig. 8 (a), at least one of our proposed PEMs leads to the highest or second-highest mean inference accuracy for all considered values of ∆t. For ∆t ≈ 0.3 and for ∆t 0.8, we obtain the highest mean inference accuracy with GC among our considered network-inference methods. When using linear Granger causality to infer networks from time-series data, one uses a VAR(p) model for the data. For ∆t close to 1, the SDD model is similar to a VAR(1) model. The VAR(p) model that one uses for the calculation of Granger causality is thus a good model for the time-series data, and using GC as a PEM leads to high inference accuracy.
To compare the results for network inference from Markovian dynamics to network inference from non-Markovian dynamics, we also show our results for mean  23) and (22)) for inferring networks. The baseline "base" corresponds to inferring a network by selecting the correct number of edges uniformly at random. In (a) and (b), we vary the sampling period ∆t. In (c) and (d), we vary the characteristic time τ . In (e) and (f), we vary the coupling strength . In (a), (c), and (e), we show the results for dynamics with no transmission lag on edges (i.e., δ = 0). In (b), (d), and (f), we show results for dynamics with a maximum transmission lag δ = 5 on edges. We sampled networks from a directed G(n, m) random-graph model. Unless specified otherwise, we use the default parameter values given in [58].
inference accuracy as one varies ∆t of a SDD model of order 6 (see Fig. 8 Figure 9. Comparison of inference methods with varying parameters of the dynamical system. We show the difference between the mean accuracy ΦLCRC of inferring networks with LCRC and the mean accuracy of inferring networks with other PEMs as we vary the coupling strength , the characteristic time τ , and the maximum edge transmission lag δ. (For values of other parameters, see [58].) Red pixels indicate that a PEM leads to higher mean accuracy than LCRC. Blue pixels indicate that a PEM leads to lower mean accuracy than LCRC. For PEM acronyms, see caption of Fig. 8. from the SDD model without transmission lags on edges (i.e., δ = 0) to the SDD model with δ = 5, one observes that the mean inference accuracy decreases for LC and OUI, which are PEMs that do not account for edge transmission lags. The mean inference accuracy for CM is larger for δ = 5 than for δ = 0. The other considered PEMs lead to similar or slightly lower mean inference accuracy for δ = 5 than for δ = 0. This negative effect of edge transmission lags on mean inference accuracy is larger for GC than for LCRC or LCCF. In Fig. 8 (b), either LCRC or LCCF yields the highest mean inference accuracy among all considered PEMs for all considered values of ∆t ∈ (0, 1].

Characteristic time
In Figure 8 (c), we show our results for mean inference accuracy as one varies the characteristic time τ of the SDD model with no transmission lags on edges (i.e., a SDD model of order 1).
We noted in Section IV B that the process motifs for covariance and lagged covariance of the SDD model depend on ∆t and τ only through the ratio ∆t τ := ∆t/τ . Accordingly, one observes similar effects when one increases τ as when one decreases ∆t. The mean inference accuracy for all considered PEMs decreases with increasing τ . In Fig. 8 (c), the PEMs LCRC, and GC tend to lead to higher mean inference accuracy than all other PEMs that we considered. For small τ (i.e., τ ≈ 1), one obtains a higher mean inference accuracy for LCCF than for LCRC or GC. As one increases τ , the ranking of the three PEMs changes first such that GC ranks best (for 2 τ 3), then such that LCRC ranks best (for τ 4).
In Figure 8 (d), we show our results for mean inference accuracy as one varies the τ of a SDD model with a maximum transmission lag δ = 5 on edges. At τ = 1, the PEM LCCF leads to the highest mean inference accuracy among our considered PEMs. However, as one increases the characteristic time to τ = 2, the mean inference accuracy for LCCF drops by almost 15%. For τ 2, the PEM LCRC leads to the highest mean inference accuracy in our comparison of PEMs.

Coupling strength
In Figure 8 (e), we show our results for mean inference accuracy as one varies the coupling strength of the SDD model with no transmission lags on edges (i.e., δ = 0).
For small to intermediate values of (i.e., 0.7), the mean inference accuracy for all considered PEMs increases monotonically with . This monotonic increase of Φ with is intuitive. At = 0, the SDD model produces only Gaussian white noise, which carries no information on a network's structure. When is very small, a signal's amplitude decreases strongly when it traverses an edge. The strong decrease of a signal's amplitude along a single edge makes it difficult to infer edges correctly. As increases, signals can traverse edges without a strong decrease of amplitude, which makes it possible to infer networks with an increasingly high accuracy.
For 0.7, the mean inference accuracy for all considered PEMs except LCCF and OUI decreases with increasing coupling strength. A possible explanation for this change in monotonicity with increasing is that, as approaches 1, there is almost no change in signal amplitude when a signal traverses an edge. Signals thus tend to persist on the network for a very long time, and the assumption that the observed time-series data is stationary becomes inaccurate. Network inference from nonstationary time-series data via PEMs tends to be more challenging than network inference from stationary timeseries data and lead to a decreased mean accuracy of inferred networks (see Subsection VI E).
In Fig. 8 (e), the PEM f (LCRC) leads to the highest mean accuracy of network inference among all considered PEMs when 0.6. For ∈ [0.7, 0.8], Granger causality is the best-performing PEM in our comparison. For = 0.9, LCCF leads to higher mean inference accuracy than any of the other considered methods.
In Figure 8 (f), we show our results for mean inference accuracy as one varies for the SDD model with transmission lags on edges. The results are qualitatively similar to the results for the SDD model without transmission lags on edges. Notably, we observe that both of our proposed PEMs have larger mean inference accuracy than GC at low coupling strengths (i.e., 0.4) for networks with and without transmission lags on edges.

Interplay of characteristic time and coupling strength
From the results that we presented in Subsections VI A 1-VI A 3, one can see that the ranking of PEMs based on their associated mean inference accuracy changes with variations in the parameters of the dynamical system that we use to generate time-series data. To give an overview of how changes in the parameters of the SDD model affect the performance of our PEMs, we show heat maps of the difference in mean inference accuracy between LCRC and several PEMs in Figure 9. Positive values (i.e., red-shaded areas) indicate parameter pairs (τ, ) for which another PEM leads to higher mean inference accuracy than LCRC. Negative values (i.e., blueshaded areas) indicate parameter pairs (τ, ) for which LCRC leads to higher mean inference accuracy than another PEM.
We find that for most combinations of ∈ (0, 1) and τ ∈ [1,10], the PEM LCRC leads to better or comparable mean inference accuracy when compared to LCCF and all other PEMs that we considered. For dynamics with small to intermediate τ (i.e., τ 5), large (i.e., 0.8) and no transmission lags on edges (i.e., δ = 0), the mean inference accuracy with LCCF or GC tends to be higher than with LCRC. For δ = 5, the mean inference accuracy associated with LCRC is either higher or less than 1% lower than the mean inference accuracy for GC.

Effect of lag-parameter mismatch
Our comparisons of mean accuracy for network inference from dynamics with and without transmission lags on edges in Fig. 8 suggest that for PEMs that do not account for transmission lags on edges (e.g., LC and OUI), the mean inference accuracy decreases as the maximum transmission lag δ increases. Most of our considered PEMs accounting for transmission lags on edges have a parameter that determines the maximum edge transmission lag that one assumes when inferring a network. Our proposed PEMs, LCRC and LCCF, include the assumed maximum transmission lagδ as a parameter (see Appendix D). For GC, the parameter setting the assumed maximum transmission time is the maximum orderp of the VAR(p) model used for regressing the time-series data. For TE, the maximum edge transmission timep is given by the parameter that sets the number of data points that one uses to compute transfer entropy. The assumed maximum transmission lag on edges for GC and TE isδ =p − 1.
For the results shown in Figs. 8 and 9, we chose the assumed maximum transmission lagδ for the network inference to be equal to the actual maximum transmission lag δ that we used in the simulation of the SDD model.
To demonstrate how a mismatch of δ andδ can affect the accuracy of inferred networks, we show the mean inference accuracy for LCRC, LCCF, GC, and TE with independently varying δ andδ in Fig. 10.
In Fig. 10 (a), we show results for increasing δ when the assumed maximum transmission lag on edges is equal to the actual maximum transmission lag on edges. We find that an increased transmission lag on edges does not have a negative effect on mean inference accuracy with LCRC, LCCF, or GC as a PEM when settingδ = δ. The mean inference accuracy for TE decreases with increasing δ even if one setsδ = δ.
We show results for mean inference accuracy for a mismatch of the assumed and actual lag parameter (i.e., δ = δ) in Fig. 10 (b). For the four considered PEMs, one obtains the best inference results when settingδ equal to or larger than δ. When one underestimates the actual maximum edge transmission lag by settingδ < δ, the mean inference accuracy for all four PEMs tends to be much lower than when the assumed maximum edge transmission lag and the actual maximum edge transmission lag match (i.e.,δ = δ). For all four PEMs with a parameter that setsδ, overestimating the maximum edge transmission lag (i.e.,δ > δ) tends to have a small negative effect or no effect on the mean inference accuracy. To demonstrate that these observations regarding a possible mismatch ofδ and δ also hold for δ = 5, we show heat maps of the mean inference accuracy for LCRC, LCCF, GC, and TE for varying δ andδ in Fig. 10 (c)-(f).

B. Inferring different network structures
For the results that we presented in Subsection VI A, we simulated dynamical systems on networks that are realizations of a G(n, m) random-graph model with n = 10 nodes, an edge density d e = 0.5, and an edge reciprocity of r e = 0.5. In the current subsection, we investigate how changes in structural properties of ground-truth networks affect the mean inference accuracy for different PEMs.
To that end, we show the mean inference accuracy for several PEMs as we vary the edge density d e , the edge reciprocity r e , and the network size n of ground-truth networks in Fig. 11.

Edge density
In Fig. 11 (a), we show our results for mean inference accuracy as one varies edge density. If the number of edges is known, the mean baseline accuracy Φ base (i.e., the accuracy achieved by guessing the positions of edges in a network uniformly at random) is close to 1 for very sparse and very dense networks and Φ base ≈ 0.5 for networks with edge density d e = 0.5. For all PEMs that we consider in this paper, the mean inference accuracy changes with increasing edge density in a similar way to  Figure 11. Accuracy with varying structural properties of networks. We compare inference results from using lag-1 correlation (LC), transfer entropy (TE), linear Granger causality (GC), convergent cross mapping (CM), Ornstein-Uhlenbeck inference (OUI), and our proposed methods LCCF and LCRC (see Eqs. (22) and (23)) to infer networks that we sampled from a directed G(n, m) random-graph model. The baseline "base" corresponds to inferring a network by selecting the correct number of edges uniformly at random. In (a) and (b), we vary the edge density de. In (c) and (d), we vary the edge reciprocity re. In (e) and (f), we vary the network size n. In (a), (c), and (e), we show the results for dynamics with no transmission lag on edges (i.e., δ = 0). In (b), (d), and (f), we show results for dynamics with a maximum transmission lag δ = 5 on edges. (For default parameters, see [58].) Each data point corresponds to the mean over the accuracy of inferred networks for 100 ground-truth networks.
Φ base . The mean inference accuracy for all considered PEMs is high when networks are very sparse or very dense and low when networks have an edge density of d e ≈ 0.5. In contrast to Φ base however, the mean infer-ence accuracy for the considered PEMs is not symmetric about d e = 0.5. Instead, all considered PEMs tend to infer networks with a small edge density d e = d with higher accuracy than networks with a correspondingly large edge density 1 − d. This asymmetry is more pronounced for LCCF and GC than for the other considered PEMs and becomes stronger for LCCF, LCRC, and GC when the SDD model includes transmission lags on edges (see Fig. 11 (b)).
For LCCF and LCRC, we hypothesize that this asymmetry is connected to the existence of many long process motifs in very dense networks. When setting a value for the correction factor α for LCCF and LCRC in Section V B, we aimed to nullify the contributions of short process motifs (i.e., the motifs (1,1) and (1,0)) to lag-1 covariance. We did not consider how the resulting correction terms affect the contributions of long process motifs because long process motifs tend to have very small contributions to lag-1 covariance per motif occurrence. However, networks that are very dense include many long process motifs whose cumulative contribution to lag-1 covariance can be large and thus can negatively affect the mean inference accuracy when inferring very dense networks with LCCF or LCRC.

Edge reciprocity
In Fig. 11 (c) and (d), we show our results for mean inference accuracy as one varies edge reciprocity of groundtruth networks with no edge transmission lags and with edge transmission lags of up to 5 timesteps, respectively. We find for both cases that, with most of the considered PEMs, one can infer networks with high edge reciprocity more accurately than networks in which most edges are not reciprocated. The increase of mean inference accuracy with increasing edge reciprocity is most pronounced when using OUI as a PEM. This PEM can only infer undirected networks accurately; edges in networks inferred via OUI are always reciprocated. The only PEM for which we find that mean inference accuracy is larger for low edge reciprocity than for high edge reciprocity is LCCF. An possible explanation for this observation is that a process motif (0, 1) (i.e., a length-1 walk from j to i) in an undirected network is always accompanied by a process motif (1, 0) (i.e., a length-1 walk from i to j), which has a negative contribution to LCCF and thus makes the detection of an undirected edge via LCCF less likely than that of a directed edge. The high mean inference accuracy for LCCF at low edge reciprocities suggests an interesting use case for LCCF as a PEM for inferring networks that one knows to have a low edge reciprocity.

Network size
In Fig. 11 (e) and (f), we show our results for mean inference accuracy as one varies the network size. Whether the dynamics include transmission lags on edges or not, the mean inference accuracy decreases for all considered PEMs with increasing network size.
In Appendix E, we compare several inferred networks to their corresponding ground-truth networks. For a set of small networks, which we call "shooting-star networks", we find that the coexistence of (1) nodes with low degree and (2) nodes with high degree and small local clustering coefficient makes network inference with LCRC difficult. This observation suggests a possible explanation for the decreasing mean inference accuracy with increasing network size: with increasing network size, the probability of nodes with high degree and small local clustering coefficient in the Erdős-Rényi model increases, and the increase of this probability makes it harder to infer these networks accurately.
Regardless of the cause of the low mean inference accuracy for large networks, our results point to a serious problem for network inference via PEMs. The PEMs that we have considered here include several PEMs that researchers have used to infer large real-world networks with hundreds of nodes from time-series data [17,59]. We find that the mean inference accuracy for these PEMs approaches Φ base with increasing network size for our simulated time-series data. It is already close to (i.e., within 10% of) baseline accuracy for networks with 50 nodes or more. Inferring large networks accurately via PEMs thus seems infeasible in practice without improved PEMs or a carefully designed preprocessing routine for time-series data.
For our default parameter values [58], we find that LCCF, LCRC, and GC tend to lead to the highest mean inference accuracies. As we have seen in Section VI A, a ranking of these three PEMs by mean inference accuracy depends on one's choice of the parameters τ , , and ∆t. We revisit the interplay between parameters of randomgraph models and the parameters of the dynamics on networks in Subsection VI B 5.

Other variations of network structure
To investigate how other aspects of the structure of ground-truth networks can affect mean inference accuracy, we compare the mean inference accuracy Φ LCRC for LCRC for four graph models: a regular ring, the Watts-Strogatz small-world model [60], the Barabási-Albert model [61], and the Erdős-Rényi model [62]. We use model variants for which one can specify the network size, edge density, and edge reciprocity as model parameters (see implmentations in our accompanying code library [31]).
In Fig. 12, we show the mean inference accuracy for LCRC, Φ LCRC , for realizations of the four considered graph models. We find that the relationships between Φ LCRC and edge density d e , edge reciprocity r e , and network size n are qualitatively similar for realizations of all four considered graph models; the mean infer-  Figure 12. Inference of realizations of graph models. We show the mean accuracy Φ of inferring networks via lag-1 correlation with a correction for reverse causation (LCRC) for variants of the Erdős-Reényi (ER) model, the Barabási-Albert (BA) model, a regular-ring (RR) model, and the Watts-Strogatz small-world (SW) model. In (a), we vary the edge density de. In (b), we vary the network size n. In (c), we vary the edge reciprocity re. In (d), we show the mean anti-clustering coefficient c with increasing network size. We observe that the inference accuracy tends to be large when c is small and vice versa. ence accuracy decreases with increasing edge density d e ∈ (0, 0.5] (see Fig. 12 (a)) or increasing network size (see Fig. 12 (b)), and it increases with increasing edge reciprocity (see Fig. 12 (c)).
For almost all parameter combinations that we considered, we obtain much higher Φ LCRC for the regularring model and the Watts-Strogatz model than for the Barabási-Albert model and the Erdős-Rényi model. A possible explanation for the differences in mean inference accuracy between the four graph models is that nodes with a high degree and small local clustering coefficient -whose neighborhoods tend to be difficult to infer (see Appendix E) -are more frequent in realizations of the Erdős-Rényi model and the Barabási-Albert model than in realizations of the regular-ring model or the Watts-Strogatz model.
To characterize nodes with a high degree and a small local clustering coefficient, we define the anti-clustering coefficient of a node i with degree k i and local clustering coefficient c i as In Fig. 12 (d), we show the mean anti-clustering coefficient c := 1 n i c i over all nodes i in a network as we increase the network size. Comparing Fig. 12 (b) with Fig. 12 (d), we observe that the mean inference accuracy tends to be high when the mean anti-clustering coefficient is small (Spearman correlation r s = −0.93 with p < 10 −100 ).

Interplay of the characteristic time of the stochastic delay-difference model and structural properties of networks
In this subsection, we demonstrate that results on the ranking of PEMs based on their associated mean inference accuracy can depend on an interplay of parameters of the dynamics on a network and the structural properties of the network.
In Fig. 13, we show the difference Φ GC −Φ LCRC of mean inference accuracy for GC and LCRC with varying network size and edge reciprocity for four different values of the characteristic time τ . We observe that Φ GC − Φ LCRC decreases and even changes from positive to negative values as one increases the network size or edge reciprocity of ground-truth networks, or the characteristic time of the SDD model. When using mean inference accuracy as the only criterion for deciding which PEM to use for network inference, one thus has to take into account structural properties of the networks that one wants to infer as well as parameters of the dynamics that give rise to the observed time-series data. We show the difference of the mean accuracy ΦGC for inferring networks with Granger causality (GC) and the mean accuracy ΦLCRC for inferring networks with lag-1 correlation with a correction for reverse causation (LCRC). We vary the network size n, the edge reciprocity re, and the characteristic time τ .
(For values of other parameters, see [58].) Red pixels indicate that GC leads to higher mean accuracy than LCRC. Blue pixels indicate that GC leads to lower mean accuracy than LCRC.

C. Effect of noisy dynamics and noisy measurements
In this subsection, we study the effect of noise on mean inference accuracy. We consider noise in the dynamics and noise in the measurement process. In the SDD model, we model noise in the dynamical system as additive Gaussian white noise in the delay-difference equation (see Eq. (12)). The parameter that tunes the strength of this system noise is σ. We model measurement noise as additive Gaussian white noise added to each observation in a time-series data set. The parameter η tunes the strength of the measurement noise.
In Fig. 14, we show our results for mean inference accuracy as one varies the strength σ of the system noise and the strength η of the measurement noise. In panels (a) and (b), we show the mean inference accuracy for several PEMs as one increases the strength of system noise for the SDD model with maximum edge transmission lag on edges of δ = 0 and δ = 5, respectively. For all considered PEMs, our results indicate that the strength of system noise does not affect mean inference accuracy. In panels (c) and (d), we show our results for increasing the strength of measurement noise. We observe that the mean inference accuracy for all considered PEMs is unaffected by measurement noise when η σ. When the strength of the measurement noise exceeds the strength of the system noise, it leads to a decrease in mean inference accuracy for all considered PEMs.

D. Nonstationarity and data efficiency
When calculating LC, GC, TE, OUI, LCRC, and LCCF, one estimates the steady-state covariance matrices or transfer entropy of a stationary stochastic process from time-series data. To ensure that the sample covariance matrix of a set of time-series data is a good approximation of the steady-state covariance matrix of the underlying stochastic process, one needs a large number of independent observations. What happens when this is not the case?
In Fig. 15, we investigate the effect of the number N of samples on the mean inference accuracy. As one increases the number of observations, one either also increases the observation time T while keeping the sampling rate 1/∆t fixed (see Fig. 15 (a)) or one increases the sampling rate while keeping T fixed (see Fig. 15 (b-c)) .

Fixed sampling rate vs. fixed observation time
When one increases N and T while keeping 1/∆t τ fixed, one successively includes more and more independent observations in the data set. This increase of independent observations tends to lead to an increase of the mean inference accuracy for all considered PEMs (see Fig. 15 (a)).
When one increases N while keeping T fixed, one adds more samples to the data set but increases the statistical dependence between subsequent data points. The the ratio ∆t τ := ∆t/τ changes, and this change can have a strong effect on mean inference accuracy Φ (see Fig. 15 (b-c)). At ∆t τ = 1, the SDD model is equivalent to a VAR model, and we thus expect PEMs that perform well on VAR models (i.e., GC, TE, LC) to achieve their highest mean inference accuracies when ∆t τ = 1. As ∆t τ decreases, the SDD model becomes more and more similar to the OUP, and we thus expect the mean inference accuracy of PEMs for continuous-time processes (i.e., OUI and CM) to increase. Our results in Fig. 15 (bc) are consistent with these hypotheses.

Varying ∆tτ in directed and undirected networks
For networks with edge reciprocity r e = 0.5, network inference via GC tends to be more accurate than network inference via OUI, even for large N and small ∆t τ . This is because it is impossible to accurately infer directed networks via OUI, because OUI is symmetric (i.e., f ). For undirected networks (i.e., r e = 1), the mean inference accuracy Φ (OUI) is close to Φ (base) when ∆t τ = 1, but it rapidly increases in with increasing N and decreasing ∆t τ . At N ≈ 250 (i.e., ∆t τ ≈ 0.4), we observe similar mean inference accuracy for OUI and GC. At N ≈ 650 (i.e., ∆t τ ≈ 0.15), the mean inference accuracy of CM is comparable to GC.
In accordance with results in previous subsections, we find that (i) LCCF and LCRC infer networks with similar or higher accuracy than all other considered PEMs when ∆t τ ∈ [0.5, 1] (see Section VI A 1), (ii) LCRC tends to infer directed networks with higher accuracy than the other considered PEMs when ∆t τ > 0.5 (see Section VI A 2), and (iii) LCCF performs much better on directed networks than on undirected networks (see Section VI B 2).
Our results demonstrate that collecting more data on a stochastic process does not necessarily improve the accuracy of inferred networks. For many PEMs, including our proposed PEMs LCRC and LCCF, one obtains the most accurately inferred networks when ∆t ≈ τ . We conclude that it is important to consider the characteristic time of a stochastic process on a network when setting the sampling rate to achieve a high accuracy of inferred networks.

E. Computation time
In Subsections VI A-VI D, we investigated how various parameters of graph models and stochastic processes on networks affect the mean inference accuracy for several PEMs. In addition to accuracy, the computational costs associated with performing network inference may be a relevant consideration when selecting a PEM. In Fig. 16, we show our results for the mean computation time for inferring a network via one of several PEMs as one varies the size of the network n, the number of samples in a time-series data set N , and the anticipated maximum edge transmission lagδ.
Our results for increasing the network size (see Fig. 16 (a)) demonstrate that the mean computation time increases fast with increasing network size for all considered PEMs. However, one can group the seven PEMs that we considered into two groups -one of computationally expensive PEMs and one of computationally inexpensive PEMs. The group of computationally expensive PEMs includes GC, TE, and CM. The group of computationally inexpensive PEMs includes LC, OUI, LCCF, and LCRC. For all considered network sizes, the mean computation time associated with any of the computationally expensive PEMs is at least 100 times longer than the mean computation time associated with any of the computationally inexpensive PEMs. The only exception to this statement is the mean computation time associated with CM at n = 1000. At this network size, a time-series data set with N = 1000 observations is too small to compute CM. In that case, the algorithm stops without returning an inferred network.
In Fig. 16 (b), we show results for increasing the number N of observations. We find that the mean computation time for all considered PEMs increases with N . The mean computation time for GC, TE, and CM increases more strongly with N than the mean computation time for the other considered PEMs. Low mean computation times for CM for N 50 occur because the algorithm aborts due to an insufficient number of observations.
For PEMs that include a parameter to set the anticipated maximum edge transmission lagδ, we show how the mean computation time changes withδ in Fig. 16 (c). We again find that we can group the considered PEMs into a group of computationally expensive PEMs, which includes GC and TE, and a group of computationally inexpensive PEMs, which includes LCCF and LCRC.
When comparing the mean inference accuracy for several PEMs in Subsections VI A-VI D, we consistently obatined the most accurate network inferences for many parameter pairs with GC, LCCF, or LCRC. When computational costs are a relevant consideration for selecting a network inference method, we believe that LCCF and LCRC offer a very good trade-off between highly accurate network inference and low computational costs. We do not find a comparably good combination of high mean inference accuracy and low computation costs for any other PEM that we considered in our study.

VII. DISCUSSION
Inferring networks from time-series data is a challenging problem with many applications. Many researchers have used vector-autoregressive (VAR) models to fit time-series data and infer networks from it via pairwise edge measures (PEMs), such as lag-1 correlation, Granger causality, and transfer entropy. We proposed a stochastic delay-difference (SDD) model (see Eq. (12)) that is a generalization of a VAR model (see Eq. (5)) and that allows for a mismatch of the characteristic time τ of the dynamics and the sampling rate 1/∆t. We then derived the process motifs for covariance and lagged covariance for the SDD model, and we used the contributions of process motifs to covariance and lagged covariance to formulate a family of PEMs that are linear combinations of the sample correlation matrix and the lag-1 correlation matrix of a time-series data set. From this family of PEMs, we selected two members that have a particularly strong theoretical motivation.
Our first proposed PEM is lag-1 correlation with a correction for confounding factors (LCCF, see Eq. (22)), which one can use to infer networks with high accuracy from mildly autocorrelated time-series data (i.e., 0.3 ∆t/τ ≤ 1). Our second proposed PEM is lag-1 correlation with a correction for reverse causation (LCRC, see Eq. (23)), which leads to accurate network-inference results for strongly autocorrelated time-series data (i.e., ∆t/τ 0.3). Importantly, one can compute both of our proposed PEMs more easily and with less computational resources than many widely-used PEMs, while obtaining similar or better inference results than the widely-used PEMs for many different network structures and for large regions of the parameter space of the SDD model. When considering inference accuracy and computational costs in the selection of a PEM for network inference, we suggest that our proposed PEMs present a better trade-off between high inference accuracy and low computational cost for the inference of networks from linear stochastic dynamics than any other PEM that we are currently aware of.
We anticipate the strong performance of our proposed PEMs is limited to the inference of linear stochastic dynamics with additive Gaussian white noise. Non-linear dynamics and dynamics with system noise that is not additive or not Gaussian are not well modelled by our SDD model, which is the theoretical base of our proposed PEMs. We do not anticipate, however, that linear Granger causality leads to much better inference results than our proposed PEMs in either of these settings. The theoretical base for linear Granger causality is a VAR model (see Eq. (5)) and any process that is not well modelled by our SDD model is also not well modelled by this VAR model.
In future work, we would like to conduct a systematic comparison of network inference via our proposed PEMs to alternative approaches to causal inference from timeseries data sets with strong autocorrelation. For example, many researchers have used a combination of preprocessing strategies (e.g. de-trending [63,64] and subsampling [65,66]) with PEMs to circumvent problems arising from strong autocorrelation in stochastic processes. A possi-ble future direction for extending the applicability of our PEMs would be to develop significance tests for LCCF and LCRC. Such tests would allow researchers to use LCCF and LCRC for network inference even when the number of edges in a network is not known a priori. Finally, our results demonstrated that the mean inference accuracy for many PEMs decreases with increasing network size. While our current paper includes results that suggest that a cause of the decline of mean inference accuracy with network size is connected to an increase of a network's mean anticlustering coefficient (see Section VI B), we anticipate that further research on the challenges posed by large networks to inference methods is a particularly interesting avenue of future research. One can compute our proposed PEMs at low computational cost for very large networks with hundreds or thousands of nodes, but we do not expect inference results for such large networks to be more accurate than guessing edge positions uniformly at random. An improved understanding of the problems that large network sizes pose for network inference via PEMs has the potential to lead to improved PEMs that can infer large networks accurately with low computational cost.
Using the binomial theorem, one can express the k-th power of K as One can thus write for which we define the 3-index sequence, g k, B , F := ζ k, B , F φ B , F , which has a scalar-valued part, which depends on k and ∆t τ , and a matrix-valued part which is independent of k and ∆t τ . The discrete-time Lyapunov equation Eq. (13) has the solution [67, p. 108] When all eigenvalues of A (1) have absolute values less than 1, the 3-index sum in Eq. (A1) is absolutely convergent, and one can thus apply Fubini's theorem for infinite series and exchange the summation indices, such that, where we use the shorthand notationˆ := max{ B , F }.
which proves Eq. (B3). One can use Eq. (B3) to derive the limit Substituting F by and B by L − F in Eq. (B6), one can see that the limit of the contributions c (0) B , F of process motifs to the steady-state covariances in the SDD model is identical to the contributions of process motifs to the steady-state covariances in the OUP (see Eq. (4)) in the limit ∆t τ → 0.
for some non-negative α. When one aims to infer edges of any transmission lag (including edges with no transmission lag) up to a maximum transmission lag δ, a suitable PEM estimator is the maximum of allf (δ e ) for δ e ∈ {0, 1, . . . , δ}, i.e., Following the same rationale that we used to motivate the LCCF and LCRC for the order-1 SDD model (see Subsection V B), we obtain the PEM estimatorŝ for the SDD model of order δ + 1.

Appendix E: Challenges in inferring large networks
To understand some of the challenges that arise in inferring large networks, we characterize structural features that tend to lead to false-postive and/or false-negative edge inferences.
We define a shooting-star network to be a graph that one can construct from connecting a path network to a star network by an edge. When each pair of adjacent nodes is connected by a bidirectional edge, one can fully characterize a shooting star network by its number n of nodes and its maximum node degree K ∈ {2, . . . , n − 1}. In Fig. 17, we show several examples of small (i.e., n = 10) shooting-star networks, and corresponding examples of networks inferred via LCRC from SDD simulations on the shooting-star networks. We also show the mean inference accuracy for shooting-star networks with increasing K/n. A comparison of the inferred networks in Fig. 17 with their corresponding ground-truth networks illustrates that (1) false-positive edge inferences tend to occur for pairs of neighbors of a high-degree node, and (2) false-negative edge inferences tend to occur for pairs of nodes with low degree. Accordingly, the mean inference accuracy for shooting-star networks is high when K/n is close to 0 or close to 1, and is low when K/n ≈ 1/2.
We explain these observations as follows. When the number of edges in a network is known, node pairs are in competition for being assigned an edge in the inferred network because a false-positive edge inference necessitates a false-negative edge inference somewhere else. One can interpret a pair of a false-positive and a false-negative edge inference as an edge in the ground-truth network having lost to a non-edge in the ground-truth network in this competition. Edges between low-degree nodes tend  At the top, we show four exemplary pairs of a ground-truth shooting-star network and a corresponding network inferred from the simulation of the SDD model on the ground-truth network. At the bottom, we show the mean accuracy of network inference via LCRC from SDD simulations on shootingstar networks. Each data point corresponds to the mean over the accuracy of inferred networks for 100 simulations with inverse sampling rate ∆t = 0.5, characteristic time τ = 1, coupling strength = 0.9, system noise σ = 0.2, measurement noise η = 0, and N = 1000 observations. to have the lowest associated PEM values among edges are thus weak competitors in this competition. Conversely, pairs of neighbors of a high-degree node tend to be associated with the highest PEM values among nonedges and thus are strong competitors. A ground-truth network that includes edges between low-degree nodes as well as non-edges between neighbors of a high-degree node (e.g., a shooting-star network with K/n ≈ 1/2) is thus hard to infer with high accuracy, when compared with a ground-truth a network that does not include edges between low-degree nodes (e.g., shooting-star networks with K/n ≈ 1) or a ground-truth network that does not include non-edges between neighbors of high-degree nodes (e.g., shooting-star networks with K/n ≈ 0). In conclusion, we observe that the coexistince of (1) edges between low-degree nodes and (2) non-edges between pairs of neighbors of a high-degree node pose a challenge for PEM-based inference methods. This observation has implications for the accuracy of PEM-based network inference with increasing network size. A network that includes high-degree nodes with a small local clustering coefficient necessarily also includes non-edges between neighbors of a high-degree node. One can use a node's anticlustering coefficient c i (see Eq. (24)) as a measure of its tendency to have a high degree and low clustering coefficient. For several graph models, the mean anticlustering coefficient increases with network size (see Fig. 12 (d)). It follows that with increasing network size the number of non-edges that can out-compete edges between low-degree nodes also increases. We expect that this is one of possibly several aspects that make the inference of large networks challenging.