{'fraction_non_alphanumeric': 0.04153031633277463, 'fraction_numerical': 0.008513581823502807, 'mean_word_length': 4.140590809628009, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variable-lag Granger Causality for Time Series Analysis

Granger causality is a fundamental technique for causal inference in time series data, commonly used in the social and biological sciences. Typical operationalizations of Granger causality make a strong assumption that every time point of the effect time series is influenced by a combination of other time series with a fixed time delay. However, the assumption of the fixed time delay does not hold in many applications, such as collective behavior, financial markets, and many natural phenomena. To address this issue, we develop variable-lag Granger causality, a generalization of Granger causality that relaxes the assumption of the fixed time delay and allows causes to influence effects with arbitrary time delays. In addition, we propose a method for inferring variable-lag Granger causality relations. We demonstrate our approach on an application for studying coordinated collective behavior and show that it performs better than several existing methods in both simulated and real-world datasets. Our approach can be applied in any domain of time series analysis.


I. INTRODUCTION
Inferring causal relationships from data is a fundamental problem in statistics, economics, and science in general. The gold standard for assessing causal effects is running randomized controlled trials which randomly assign a treatment (e.g., a drug or a specific user interface) to a subset of a population of interest, and randomly select another subset as a control group which is not given the treatment, thus attributing the outcome difference between the two groups to the treatment. However, in many cases, running such trials may be unethical, expensive, or simply impossible [1]. To address this issue, several methods have been developed to estimate causal effects from observational data [2], [3].
In the context of time series data, a well-known method that defines a causal relation in terms of predictability is Granger causality [4]. X Granger-causes Y if past information on X predicts the behavior of Y better than Y 's past information alone [5]. In this work, when we refer to causality, we mean specifically the predictive causality defined by Granger causality. The key assumptions of Granger causality are that 1) the process of effect generation can be explained by a set of structural equations, and 2) the current realization of the effect at any time point is influenced by a set of causes in the past. Similar to other causal inference methods, Granger causality assumes unconfoundedness and that all relevant variables are included in the analysis. There are several studies that have been developed based on Granger causality [6]- [8]. The typical operational definitions [7] and inference methods for inferring Granger causality, including the common software implementation packages [9], [10], assume that the effect is influenced by the cause with a fixed and constant time delay.
This assumption of a fixed and constant time delay between the cause and effect is, in fact, too strong for many applications of understanding natural world and social phenomena. In such domains, data is often in the form of a set of time series and a common question of interest is which time series are the (causal) initiators of patterns of behaviors captured by another set of time series. For example, who are the individuals who influence a group's direction in collective movement? What are the sectors that influence the stock market dynamics right now? Which part of the brain is critical in activating a response to a given action? In all of these cases, effects follow the causal time series with delays that can vary over time.
To address the remaining gap, we introduce the concept Variable-lag Granger causality and a method to infer it in time series data. We prove that our definition and the proposed inference method can address the arbitrary-timelag influence between cause and effect, while the traditional operationalizations of Granger causality and the corresponding inference methods cannot. We show that the traditional definition is indeed a special case of the new relation we define. We demonstrate the applicability of the newly defined causal inference framework by inferring initiators of collective coordinated movement, a problem proposed in [11].
We use Dynamic Time Warping (DTW) [12] to align the cause X to the effect time series Y while leveraging the power of Granger causality. In the literature, there are many clustering-based Granger causality methods that use DTW to cluster time series and perform Granger causality only for time series within the same clusters [13], [14]. Previous work on inferring causal relations using both Granger causality and DTW has the assumption that the smaller warping distance between two time series, the stronger the causal relation is [15]. If the minimum distance of elements within the DTW optimal warping path is below a given distance threshold, then the method considers that there is a causal relation between the two time series. However, their work assumes that Granger causality and DTW should run independently.In contrast, our method formalizes the integration of Granger causality and DTW by generalizing the definition of Granger causality itself and using DTW as an instantiation of the optimal alignment requirement of the time series.
In addition to the standard uses of Granger causality, our method is capable of: • Inferring arbitrary-lag causal relations: our method can infer Granger causal relation where a cause influences an effect with arbitrary delays that can change dynamically; • Quantifying variable-lag emulation: our method can report the similarity of time series patterns between the cause and the delayed effect, for arbitrary delays; We also prove that when multiple time series cause the behavioral convergence of a set of time series then we can treat the set of these initiating causes in the aggregate and there is a causal relation between this aggregate cause (of the set of initiating time series) and the aggregate of the rest of the time series. We provide many experiments and examples using both simulated and real-world datasets to measure the performance of our approach in various causality settings and discuss the resulting domain insights. Our framework is highly general and can be used to analyze time series from any domain.

II. RELATED WORK
Many causal inference methods assume that the data is i.i.d. and rely on knowing a mechanism that generates that data, e.g., expressed through causal graphs or structural equations [2]. In time series data, the values of the consecutive time steps violate the i.i.d. assumption. Another set of causal inference methods relax the strong i.i.d assumption, and instead assume independence between the cause and the mechanism generating the effect [16]- [18]. Specifically, knowing the cause X = x never reveals information about the structural function f (X) and vice versa. This idea has been used in the context of times series data [18] by relying on the concept of Spectral Independence Criterion (SIC). If a cause X is a stationary process that generates the effect Y via linear time invariance filter h (mechanism), then X and h should not contain any information about each other but dependency between them and Y exists in spectral sense.
Granger causality has inspired a lot of research since its introduction in 1969 [4]. Recent work on Granger causality has focused on various generalizations for it, including ones based on information theory, such as transfer entropy [19], [20] and directed information graphs [21]. Recent inference methods are able to deal with missing data [22] and enable feature selection [23]. Granger causality has even been explored as a method to offer explainability of machine learning models [24]. However, none of them study tests for variablelag Granger causality, as we propose in this work. Besides, no method studies a causal structure that is unstable overtime [25]. In our work, we also relax the stationary assumption of time series.

III. GRANGER CAUSALITY AND FIXED LAG LIMITATION
Let X = (X(1), . . . , X(t), . . . ) be a time series. We will use X(t) ∈ R to denote an element of X at time t. Given two time series X and Y , it is said that X Granger-causes [4] Y if the information of X in the past helps improve the prediction of the behavior of Y , over Y 's past information alone [5]. The typical way to operationalize this general definition of Granger causality [7] is to define it as follows: Definition 1 (Granger causal relation): Let X and Y be time series, and δ max ∈ N be a maximum time lag. We define two residuals of regressions of X and Y , r Y , r Y X , below: where a i and b i are constants that optimally minimize the residual from the regression. Then X Granger-causes Y if the variance of r Y X is less than the variance of r Y . This definition assumes that, for all t > 0, Y (t) can be predicted by the fixed linear combination of with some fixed ∆ > 0 and every a i , b i is a fixed constant over time [5], [7]. However, in reality, two time series might influence each other with a sequence of arbitrary, non-fixed time lags. For example, Fig. 1(a2.) has X as a cause time series and Y as the effect time series that imitates the values of X with arbitrary lags. Because Y is not affected by X with a fixed lags and the linear combination above can change over time, the standard Granger causality tests cannot appropriately infer Granger-causal relation between X and Y even if Y is just a slightly distorted version of X with some lags. For a concrete example, consider a movement context where time series represent trajectories. Two people follow each other if they move in the same trajectory. Assuming the followers follow leaders with a fixed lag means the followers walk lockstep with the leader, which is not the natural way we walk. Imagine two people embarking on a walk. The first starts walking, the second catches up a little later. They may walk together for a bit, then the second stops to tie the shoe and catches up again. The delay between the first and the second person keeps changing, yet there is no question the first sets the course and is the cause of the second's choices where to go. Fig. 1 illustrates this example.

IV. VARIABLE-LAG GRANGER CAUSALITY
Here, we propose the concept of variable-lag Granger causality, VL-Granger causality for short, which generalizes the Granger causal relation of Definition 1 in a way that addresses the fixed-lag limitation. We demonstrate the application of the new causality relation for a specific application of inferring initiators and followers of collective behavior.
Definition 2 (Alignment of time series ): An alignment between two time series X and Y is a sequence of pairs of indices (t i , t j ), aligning X(t i ) to Y (t j ), such that for any two pairs in the alignment (t i , t j ) and (t i , t j ), if t i < t i then t j < t j (non-crossing condition). The alignment defines a sequence of delays P = (∆ 1 , . . . , ∆ t , . . . ), where ∆ t ∈ Z and X(t − ∆ t ) aligns to Y (t).
Definition 3 (VL-Granger causal relation): Let X and Y be time series, and δ max ∈ N be a maximum time lag (this is an upper bound on the time lag between any two pairs of time series values to be considered as causal). We define residual r * Y X of the regression: is a time delay constant in the optimal alignment sequence P * of X and Y that minimizes the residual of the regression. The constants a i , b i , and c i optimally minimize the residuals r Y , r Y X , and r * Y X , respectively. The terms b i and c i can be combined but we keep them separate to clearly denote the difference between the original and proposed VL-Granger causality. We say that X VL-Granger-causes Y if the variance of r * Y X is less than the variances of both r Y and r Y X . In order to make Definition 3 fully operational for this more general case (and to find the optimal constants values), we need a similarity function between two sequences which will define the optimal alignment. We propose such a similaritybased approach in Definition 4. Before defining this approach, we show that VL-Granger causality is the proper generalization of the traditional operational definition of Granger causality stated in Definition 1. Clearly, if all the delays are constant then r * Y X (t) = r Y X (t). Proposition 4.1: Let X and Y be time series and P be their alignment sequence. If ∀t, ∆ t = ∆, then r * Y X (t) = r Y X (t).
We must also show that the variance of r * Y X (t) is no greater than the variance of r Y X (t).
Proposition 4.2: Let X and Y be time series, P = (∆ 1 , . . . , ∆ t , . . . ) be their alignment sequence such that According to Propositions 4.1 and 4.2, VL-Granger causality is the generalization of the Def. 1 and always has lower or equal variance.
Of a particular interest is the case when there is an explicit similarity relation defined over the domain of the input time series. The underlying alignment of VL-Granger causality then should incorporate that similarity measure and the methods for inferring the optimal alignment for the given similarity measure.
Note, here the sim similarity function does not have to be a distance function that obeys, among others, a triangle inequality. It can be any function that quantitatively compares the two time series. For example, it may be that when one time series increases the other decreases. We provide a more concrete and realistic example in the application setting below.
Adding this similarity measure to Definition 3 allows us to instantiate the notion of the optimal alignment P * as the one that maximizes the similarity between X and Y : where X (t) = X(t − ∆ t ) for any given P and ∆ t ∈ P . With that addition, if X Y , then X VL-Granger-causes Y . This allows us to operationalize VL-Granger causality by checking for variable-lag emulation, as we describe in the next section.

A. Example application: Initiators and followers
In this section, we demonstrate an application of the VL-Granger causal relation to finding initiators of collective behavior. The Variable-lag emulation concept corresponds to a relation of following in the leadership literature [11]. That is, X ≺ Y if Y is a follower of X. We are interested in the phenomenon of group convergence to a consensus behavior and answering the question of which subset of individuals, if any, initiated that collective consensus behavior. With that in mind, we now define the concept of an initiator and provide a set of subsidiary definitions that allow us to formally show (in Proposition 4.3) that initiators of collective behavior are indeed the time series that VL-Granger-cause the collective pattern in the set of the time series. In order to do this, we generalize our two-time series definitions to the case of multiple time series by defining the notion of an aggregate time series, which is consistent with previous Granger causality generalizations to multiple time series [26]- [28].
Definition 5 (Initiators): Let U = {U 1 , . . . , U n } be a set of time series. We say that X ⊆ U is a set of initiators if ∀U ∈ U \ X , ∃X ∈ X , s.t.X ≺ U , and, conversely, ∀X ∈ X ∃U ∈ U \ X , s.t.X ≺ U . That is, every time series follows some initiator and every initiator has at least one follower.
Given a set of time series U = {U 1 , . . . , U n } , and a set of time series X ⊆ U, we can define an aggregate time series as a time series of means at each step: In order to identify the state of reaching a collective consensus of a time series, while allowing for some noise, we adopt the concept of -convergence from [29].
Definition 6 ( -convergence): Let Q and U be time series, dist : R 2 × [0, 1] be a distance function, and 0 < ≤ 1/2. If for all time t ∈ [t 0 , t 1 ], dist(Q(t), U (t)) ≤ , then Q and U -converge toward each other in the interval [t 0 , t 1 ]. If t 1 = ∞ then we say that Q and U -converge at time t 0 .
Definition 7 ( -convergence coordination set): Given a set of time series U = {U 1 , . . . , U n }, if all time series in Uconverge toward agg(U), then we say that the set U is an -convergence coordination set. We are finally ready to state the main connection between initiation of collective behavior and VL-Granger causality. Proposition 4.3: Let dist : R 2 ×[0, 1] be a distance function, U be a set of time series, and X ⊆ U be a set of initiators, which is an -convergence coordination set converging towards T .
If for any U, U ∈ U their similarity sim(U, U ) ≥ 1 − in the interval [t 0 , t 1 ], then agg(X ) VL-Granger-causes agg(U \ X ) in that interval.
Proof Suppose ∀X ∈ X , X and agg(X ) -converge toward each other in the interval [t 0 , t 1 ], then, by definition, for all the times t ∈ [t 0 , t 1 ], dist(agg(X )(t), X(t)) ≤ . By the definition of initiators, ∀U ∈ U \ X , ∃X ∈ X , such that X ≺ U , from some time t 2 > t 0 . Thus, we have ∀t, We have now shown that a subset of time series are initiators of a pattern of collective behavior of an entire set if that subset VL-Granger-causes the behavior of the set. Thus, VL-Granger causality can solve the COORDINATION INITIATOR INFERENCE PROBLEM [11], which is a problem of determining whether a pattern of collective behavior was spurious or instigated by some subset of initiators and, if so, finding those initiators who initiate collective patterns that everyone follows.

V. VL-GRANGER CAUSALITY INFERENCE
Given a target time series Y , a candidate causing time series X, a similarity threshold σ, a significance level α, and the max lag δ max , our framework evaluates whether X VL-Granger causes Y (with a variable lag), X Granger causes Y (with a fixed lag) or no conclusion of causation between X and Y .
In Algorithm 1 line 1-2, we have a fix-lag parameter F ixLag that controls whether we choose to compute the normal Granger causality (F ixLag = true) or VL-Granger causality (F ixLag = f alse). We present the high level logic of the algorithm. However, the actual implementation is more efficient by removing the redundancies of the presented logic.
First, we compute Granger causality (line 1 in Algorithm 1). The flag GrangerResult = true if X Granger-causes Y , otherwise GrangerResult = f alse. Second, we compute VL-Granger causality (line 2 in Algorithm 1). The flag Based on the work in [7], we use the Bayesian Information Criteria (BIC) to compare the residual of regressing Y on Y past information, r Y , with the residual of regressing Y on Y and X past information r Y X . We use v 1 v 2 to represent that v 1 is less than v 2 with statistical significance by using some statistical test(s). If BIC(r Y ) BIC(r Y X ), then we conclude that the prediction of Y using Y, X past information is better than the prediction of Y using Y past information alone. After we got the results of both GrangerResult and V LGrangerResult, then we proceed to report the conclusion of causal relation between X and Y w.r.t. the following four conditions. • If both GrangerResult and V LGrangerResult are true, then we compare the residual of variable-lag regression r DT W with bothr Y and r Y X . If BIC(r DT W ) < min(BIC(r Y X ), BIC(r Y )), then we conclude that X causes Y with variable lags, otherwise, X causes Y with a fix lag (line 3 in Algorithm 1). • If GrangerResult is true but V LGrangerResult is false, then we conclude that X causes Y with a fix lag (line 4 in Algorithm 1). • If GrangerResult is false but V LGrangerResult is true, then we conclude that X causes Y with variable lags (line 5 in Algorithm 1). • If both GrangerResult and V LGrangerResult are false, then we cannot conclude whether X causes Y (line 4 in Algorithm 1). Algorithm 1: Time-lag test function input : X, Y , σ, α, δ max output: XGrangerCausesY Note that we assume the maximum lag value δ max is given as input, as it is for all definitions of Granger causality. For practical purposes, a value of a large fraction (e.g., half) of the length of the time series can be used. However, there is, of course, a computational trade-off between the magnitude of δ max and the time it takes to compute Granger causality by almost all methods.
In the next section, we describe the details of VL-Granger function that we use in Algorithm 1: line 1-2.
A. VL-Granger causality operationalization . As shown in the previous section, we can operationalize VL-Granger causality by checking for variable-lag emulation. Given time series X, Y , a similarity threshold σ, a significance level α, the maximum possible lag δ max , and whether we want to check for variable or fixed lag F ixLag, Algorithm 2 reports whether X causes Y by setting GrangerResult to true or false, and by reporting on two residuals r Y and r Y X .
First, we compute the residual r Y of regressing of Y on Y 's information in the past (line 1). Then, we regress Y (t) on Y and X past information to compute the residual r Y X (line 2). If BIC(r Y X ) BIC(r Y ), then X Granger-causes Y and we set GrangerResult = true (line 3). If F ixLag is true, then we report the result of typical Granger causality. Otherwise, we consider VL-Granger causality (lines 3-7) by computing the emulation relation between X DT W and Y where X DT W is a version of X that is reconstructed through DTW and is most similar to Y , captured by DT W reconsF unc(X, Y ) which we explain in Section V-B.
Afterwards, we do the regression of Y on X DT W 's past information to compute residual r DT W (line 4). Finally, we check whether BIC(r DT W ) BIC(r Y ) and the similar value simV alue ≥ σ (line [8][9][10][11]. If so, X VL-Granger-cause Y . In the next section, we describe the details of how to construct X DT W and how to estimate the emulation similarity value simV alue. In this work, we propose to use Dynamic Time Warping (DTW) [12], which is a standard distance measure between two time series. DTW calculates the distance between two time series by aligning sufficiently similar patterns between two time series, while allowing for local stretching (see Figure 1). Thus, it is particularly well suited for calculating the variable lag alignment.
Given time series X and Y , Algorithm 3 reports reconstructed time series X DT W based on X that is most similar to Y , as well as the emulation similarity simV alue between the two series. First, we use DT W (X, Y ) to find the optimal alignment sequenceP = (∆ 1 , . . . , ∆ t , . . . ) between X and Y , as defined in Definition 2. Efficient algorithms for computing DT W (X, Y ) exist and they can incorporate various kernels between points [12], [30]. Then, we useP to construct X DT W where X DT W (t) = X(t − ∆ t ). Afterwards, we use X DT W to predict Y instead of using only X information in the past in order to infer a VL-Granger causal relation in Definition 3. The benefit of using DTW is that it can match time points of Y and X with non-fixed lags (see Figure 1). LetP = (∆ 1 , . . . , ∆ t , . . . ) be the DTW optimal warping path of X, Y such that for any ∆ t ∈P , Y (t) is most similar to X(t − ∆ t ).
In addition to finding X DT W , DT W ReconstructionF unction estimates the emulation similarity simV alue between X, Y in line 3. For that, we adopt the measure from [11] below: Since the sign(∆ t ) represents whether Y is similar to X in the past (sign(∆ t ) > 0) or X is similar to Y in the past (sign(∆ t ) < 0), by comparing the sign of sign(∆ t ), we can infer whether Y emulates X. The function s(P ) computes the average sign of sign(∆ t ) for the entire time series. If s(P ) is positive, then, on average, the number of times that Y is similar to X in the past is greater than the number of times that X is similar to some values of Y in the past. Hence, s(P ) can be used as a proxy to determine whether Y emulates X or vice versa.

VI. EXPERIMENTS
We measured our framework performance on the task of inferring causal relations using both simulated and real-world datasets. The notations and symbols we use in this section are in Table I.

A. Experimental setup
We tested the performance of our method on synthetic datasets, where we explicitly embedded a variable-lag causal relation, as well as on biological datasets in the context of the application of identifying initiators of collective behavior.
We compared our method, VL-Granger causality (VG) with several existing methods: Granger causality with F-test (G) [7], Copula-Granger method (CG) [6], and Spectral Independence Criterion method (SIC) [18]. We have two different versions of VG, VG All and VG No F-test. VG All consider a relation to be a causal relation only if it has a strong predictability property (verified by F-test to check BIC(r Y X ) BIC(r Y )), strong emulation relation (simV alue ≥ σ), and a strong dependency between the effect and the past information of the cause (verified by the independence test that will be discussed in Section VI-A). VG No F-test is the same as VG All except that a causal relation can have a weak predictability property (verified by the condition BIC(r Y X ) < BIC(r Y ) without using F-test).
In this paper, we explore the choice of δ max in {0.1T, 0.2T, 0.3T, 0.4T } for all methods to analyze the sensitivity of each method, where T is the length of time series, and set σ = 0.5 as default unless explicitly stated otherwise.
Independence test. To infer Granger causality, we estimate the regression between Y and the past part of X, to check whether X predicts Y . However, if X causes Y , then, in addition to predicting, Y must have dependency with some pattern of X in the past. By checking for dependency between Y and an aligned version of X, we can have higher confidence that X might indeed be the cause of Y .
We deploy Hilbert-Schmidt Independence Criterion (HSIC) test [31], which outperforms contingency table and correlationbased tests [31], for testing the dependency between Y and the past of X. The null hypothesis H 0 is that the given variables are independent. We use α as a significance level needed to reject H 0 . We check the dependency between Y starting at time t + 1 with all the shifted versions of X and X DT W , starting from time 1 to t. If Y is simply a shifted version of X or a non-fixed-lag distortion of X, then Y must have a dependency with some version of X and HSIC must reject H 0 . We use HSIC with bootstrap resampling to estimate the test threshold.

B. Datasets
Synthetic data: pairwise level. The main purpose of the synthetic data is to generate settings that explicitly illustrate the difference between the original Granger causality method and the proposed variable-lag approach. We generate pairs of time series for which the fixed-lag Granger causality methods would fail to find a relationship but the variable-lag approach would find the intended relationships.
We generated a set of synthetic time series of 200 time steps. We generated two sets of pairs of time series X and Y . First, we generated X either by drawing the value of each time step from a normal distribution with zero mean and a constant variance (X(t) ∼ N ) or by Auto-Regressive Moving Average model (ARMA or A.) with X(t) = 0.9µ + 0.1X(t − 1).
The first set we generated was of explicitly related pairs of time series X and Y , where Y emulates X with some time lag ∆ ∈ [1, 20] (X ≺ Y ). One way to ensure lag variability is to "turn off" the emulation for some time. In our data, Y remains constant between 110th and 150th time steps. This makes Y a variable-lag follower of X. Figure 3 shows examples of the generated time series.
The second set of time series pairs X and Y were generated independently and as a result have no causal relation. We used these pairs to ensure that our method does not infer spurious relations.
We set the significance level for both F-test and independence test at α = 0.01. We considered there to be a causal relation only if simV alue ≥ σ for our method.  2. The causal graph where the edges represent causal directions from the cause time series (e.g. X 1 ) to the effect time series (e.g. Y 1 ). Y ij represents a time series generated by agg({X i , X j }), where X i ≺ X i with some fixed lag ∆.
Synthetic data: group level. This experiment explores the ability of causal inference methods to retrieve multiple causes of a time series Y ij , which is generated from multiple time series X i , X j . Fig. 2 shows the ground truth causal graph we used to generate simulated datasets. The edges represent causal directions from the cause time series (e.g. X 1 ) to the effect time series (e.g. Y 1 ). Y ij represents the time series generated by agg({X i , X j }), where X i ≺ X i and X j ≺ X j with some fixed lag ∆ ∈ [1,20]. The task is to infer edges of this causal graph from the time series. We generated time series for each generator model 100 times.
Schools of fish. We used the dataset of golden shiners (Notemigonus crysoleucas) that is publiclly available. The dataset has been collected for the study of information propagation over the visual fields of fish [32]. There are 24 Fig. 3. The comparison between the original time series X, variable-lag follower Y , fixed-lag time series modified from X to match Y , and variablelag time series modified from X to match Y . The traditional Granger causality uses only fixed-lag version of X to infer whether X causes Y , while our approach uses both versions of X to determine the causality between X, Y . Both X, Y are generated from N . Y remains constant from time 110 to 150, which makes it a variable-lag follower of X. coordination events in this dataset. Each coordination event consists of two-dimensional time series of fish movement that are recorded by video. The time series of fish movement are around 600 time steps. The number of fish in each dataset is around 70 individuals, of which 10 individuals are "informed" fish who have been trained to go to a feeding site. Trained fish lead the group to feeding sites while the rest of the fish just follow the group. We represent the dataset as a pair of aggregated time series: X being the aggregated time series of the directions of trained fish and Y being the aggregated   [33], [34]. For the analysis, we selected the 16 members of the troop that have GPS information available for 10 consecutive days, with no missing data. We selected a set of trajectories of lat-long coordinates from a highly coordinated event that has the length of 600 time steps (seconds) for each baboon. This known coordination event is on August 02, 2012 in the morning, with the baboon ID3 initiating the movement, followed by the rest of the troop [11]. Again, the goal is to infer ID3 (time series X) as the cause of the movement of the rest of the group (aggregate time series Y ). In this experiment we set δ max = 120 seconds, informed by the biological meaning and expertise.

C. Time complexity and running time
The main cost of computation in our approach is DTW. We used the "Windowing technique" for the search area of warping [35]. The main parameter for windowing technique is the maximum time delay δ max . Hence, the time complexity of our approach is O(T δ max ). Table II

VII. RESULTS
We report the results of our proposed approach and the standard Granger causality method on both synthetic and biological data. We also explore how the performance of the methods depends on the basic parameter, δ max .
A. Synthetic data: pairwise level Table III (1st-5th rows) shows the results of the accuracy of inferring causal relations and directions. For each row, we repeated the experiment 100 times on different simulated datasets and computed the accuracy and reported the mean. The result shows that our methods, VG (all) and (No F-test), performed better than the rest of other methods. Moreover, we also investigate the sensitivity of varying the value of the x-axis represents the value of δmax as a fraction of the time series length T and y-axis is the average accuracy.
δ max parameter for all methods. We aggregated the accuracy of inferring causal direction from various cases that have the same δ max value and report the result. The result in Fig. 4 shows that our approaches: VG (ALL) and VG (No F-test), can maintain the high accuracy throughout the range of the values of δ max .
B. Synthetic data: group level Table IV shows the result of causal graph inference. The VG (All), VG (No F-test), and G (1st-3rd rows) performed the best overall with the highest F1 score in the normal generator model (N ). However, VG (All) and VG (No Ftest) have higher F1 score than G in ARMA generator model. This result reflects the fact that our approaches can handle complicated time series in the Granger causal inference task better than the rest of other methods. In addition, we aggregated X = agg({X 1 , X 2 , X 3 }) and the rest of time series Y = agg({Y 1 , Y 2 , . . . , Y 123 }), then we measured the ability of methods to infer that X is a cause of Y . The results, which are in the "Group N " and "Group A." rows in Table III, show that VG (All), VG (No F-test), and G performed well in this task, while CG and SIC failed to infer causal relations. C. Animal data: schools of fish Table V shows the result of inferring causal direction that X (trained fish) cause Y (untrained fish) group behavior, as observed from 24 traces of coordinated movements. Granger (G) performed poorly while VG, CG, and SIC performed well on this task. This result reflects the fact that the causal relation between X and Y have highly variable lags of influence between causes and effects. Granger causality typically dismisses the existence of this type of causation and, therefore, misses the behavior.

D. Animal data: a troop of baboons
Recall, that in this dataset ID3 is the true initiator of the group's behavior during the observed interval of coordinated behavior. Fig. 5 shows the result of initiator inference in this coordinated interval. Only VG (all) in the top-left of the Fig. 5 gave the correct initiator (ID3) the highest weight. G did not differentiate among a large group of individuals as the potential cause. CG gave a week evidence of 0.1 to ID4 (wrong) as the initiator. Finally, SIC gave equal weak support of less than 0.3 to both ID3 and ID21 (false positive). This result shows the robustness of our approach over other methods in the real-world noisy datasets. According to [18], this implies that ID3 has a spectral independence from the collective time series Y = agg(U \ {ID3}), while Y has spectral dependence with ID3. This result shows that not only can we gain insight into whether X VG-causes Y from our method, but we can use SIC (and perhaps other methods) to confirm the causal properties between X and Y . ID3 is the only true initiator in the dataset. Each plot represents the inferred level of each individual (x-axis) being the cause of the group movement (y-axis: higher value is better). As shown, VG (all) is able to uncover the true initiator as the cause of the group's movement, unlike G, CG and SIC (not the y-axis scale).

SUGGESTIONS
One of the main input of Granger causality is the max lag variable δ max . There is a report in [36] that Granger causality might face the overfitting issue when δ max is inappropriately chosen. The VL-Granger causality also has the same issue. The cross-validation technique might be used to combat this issue when there is no ground-truth of δ max is available. Moreover, since VL-Granger causality is the generalization of traditional Granger method, it can easily face the issue of overfittingfalse positive discovery of causality relation. This situation is similar to the case of fitting linear vs. non-linear regression. The more generalized model seems to suffer from the issue of overfitting easier than a simpler model. We suggest that if the variance of VL-Granger is not different from traditional method, the conclusion of traditional Granger approach should be used.
Regarding the scalability aspect, even though VL-Granger causality is a generalization of the traditional method, it requires more resources in both time and space due to the need of inferring variable lags. The efficient approach of inferring variable lags should be developed in future.

IX. CONCLUSIONS
In this work, we proposed a method to infer Granger causal relations in time series where the causes influence effects with arbitrary time delays, which can change dynamically. We formalized a new Granger causal relation, proving that it is a true generalization of the traditional Granger causality. We demonstrated on both carefully designed synthetic datasets and noisy biological data that the new causal relation can address the arbitrary-time-lag influence between cause and effect, while the traditional Granger causality cannot. Moreover, in addition to improving and extending Granger causality, our approach can be applied to infer leader-follower relations, as well as the dependency property between cause and effect. We have shown that, in many situations, the causal relations between time series do not have a lock-step connection of a fixed lag that the traditional Granger causality assumes. Hence, traditional Granger causality misses true existing causal relations in such cases, while our method correctly infers them. Our approach can be applied in any domain of study where the causal relations between time series is of interest. In future work, we plan to explore non-linear relations between time series using non-linear kernels. We also provide the source code at [37].