Characterizing time series: when Granger causality triggers complex networks

In this paper, we propose a new approach to characterize time series with noise perturbations in both the time and frequency domains by combining Granger causality and complex networks. We construct directed and weighted complex networks from time series and use representative network measures to describe their physical and topological properties. Through analyzing the typical dynamical behaviors of some physical models and the MIT-BIH human electrocardiogram data sets, we show that the proposed approach is able to capture and characterize various dynamics and has much potential for analyzing real-world time series of rather short length.

3 because those interactions even in different frequency bands could be clearly figured out instead of a single number [12]. Therefore, all these important advantages of the GC allow us to have a deeper insight into the characterization of the dynamics of the time series through the established Cau-CNs. More specifically, we will use three frequently used measures, the average path length (APL), the average vertex strength (AVS) and the average clustering coefficient (ACC) [13,14], computed for the so constructed Cau-CNs. Interestingly, these measures in both the time and frequency domains are verified to be closely connected to typical dynamical systems. In particular, the APL and ACC show more potential in distinguishing dynamical behaviors of given time series in all frequencies. In addition, to exemplify the broad perspective of the proposed approach in real applications, this paper investigates human electrocardiogram (ECG) recordings from the MIT-BIH data sets.

Methods
In this section, we describe the method to construct Cau-CNs by computing the GC in the context of the linear autoregressive (AR) model in both the time and frequency domains and utilize a static as well as a dynamic approach to characterizing dynamics of time series through the network to be established.

Causal complex network construction
Given a time series consisting of n observations, denoted by {x[q]} n q=1 , the time series-induced Cau-CN is designed as follows. Firstly, we divide the given time series equally into m pieces, which are called sub-time series, denoted by x j (t) = x[t + ( j − 1) n/m ] with j = 1, . . . , m and t = 1, . . . , n/m , where · is the floor function. Later, we will discuss the impact of the piece number as well as the length of each sub-time series. Secondly, we suppose that each sub-time series corresponds to one vertex in the finally obtained Cau-CN, denoted by (V, W) J . Here, V is a set consisting of all m vertices in the network, and W = {w i j } m×m is a directed weight matrix, in which each of its elements describes the GC from the vertices i to j. The subscript J represents an index indicating the network constructed in the time or frequency domains. We include in appendix the exact definitions of GC in both the time and the frequency domains. In particular, when J = T , a Cau-CN (V, W) T in the time domain is constructed, namely each weight connectivity w i j of the network is computed through the formula of the causal influence F i→ j (see equation (A.5)); meanwhile, J = f theoretically invites an infinite number of networks (V, W) f in the frequency domain, namely the weight connectivities with frequency f are obtained by the causal influence I i→ j ( f ) (see equation (A.6)), where the frequency f ∈ [0, +∞).

Network measures-a static characterization
Once the Cau-CN (V, W) J is constructed in either the time or frequency domains, not only the distributions of the weights but also the three aforementioned measures are calculated.
1. AVS describes the average of the strengths for all the network vertices [13,14]. The strength of vertex i is defined as s tot i = s in i + s out i , where the in-strength s in i = j =i w ji and 4 the out-strength s out i = j =i w i j . Also the degree of vertex i is defined as d tot = j =i a i j and a i j is the element of the adjacency matrix A corresponding to the network (V, W) J . 2. APL describes the average of the shortest geodesic path lengths between all pairs of network vertices, where the shortest geodesic path refers to the path between two vertices with the sum of the weights of its constituent edges minimized [13,14]. 3. ACC represents the average fraction of the connections between the topological neighbors of a vertex with respect to the maximum possible for all network vertices [13,14]. In particular, the clustering coefficient (CC) of vertex i is defined as It should be noted that, in the following discussion, we use only the above three measures to characterize the constructed networks in both domains. However, there are some other typical measures that are frequently used in the characterization of complex networks [13,14]; they have been included in our ongoing works.

Artificial neural networks-a dynamic characterization
In addition to the above static measures, we design another distinguishing method by observing the dynamics of a discrete-time artificial neural network (DANN) endowed with the weighted connection matrix obtained above. The DANN is a mathematical model inspired by the structure and functional aspects of biological neural networks which defines a discrete map to describe the evolution of the states of neurons over time. Concretely, these neurons are connected and we endow the weighted connection matrix {w i j } of a DANN:  [15]. The dynamics of the time series can then be characterized by the evolution of the DANN over time.
All these characterization methods are illustrated and validated in the following examples.

Examples: from toy models to real-world data
Example 1 (Non-chaos versus chaos). Consider a noised-perturbed time series: where ε[q] is a zero mean independent Gaussian noise with variance 0.01, and each x[q] is generated by the standard Logistic model: Here, with an increase of the parameter λ in the interval [2,4], the dynamics of the Logistic model show qualitative changes from stable fixed points to cascades of periodic orbits and chaos, as shown in figure 1. For each λ and, correspondingly, the noise-perturbed series y[q] that contains 20 sub-time series, we establish Cau-CNs by the approach we proposed above.
According to the definition of the GC, the inclusion of the knowledge about one random sub-time series cannot significantly reduce the variance of the prediction error of another random one. Nevertheless, the sub-time series in the periodic series, with stronger mutual causal relations, do bring about the reduction. It is thus expected that the variation of the APL, AVS and ACC of the causal network in the time domain is closely connected to the dynamics of the Logistic model with different λ. As shown in figure 1, when λ ∈ [2, 3] that corresponds to the appearance of stable fixed points, y[q] that consists of stable fixed point dynamics and noise perturbation is seen as a random series, so that the weak causal internal relations between its sub-time series naturally result in small values of the network measures. The measures become larger when λ is taken from the region where stable periodic dynamics appear and again become smaller as λ is in the parameter region of chaos. Note also that, although the measures for the chaotic dynamics are relatively small, most of them are still larger than those for the random series as λ ∈ [2,3]. This somewhat implies the difference between chaos and randomness from a casual network point of view. In addition, the peaks of the measurement curves in figure 1 correspond to the islands of stable periodic dynamics embedded in the chaotic region, which further illustrates the reliability of the proposed approach in differentiating between chaos and non-chaotic dynamics. To be candid, this kind of differentiation can be analogously done by some existing techniques [5][6][7]; however, the proposed approach is further capable of showing the frequency bands where the differences appear between the causal networks corresponding to different dynamics. As shown in figure 1, the three measures become larger in both lower and some higher frequency bands when λ takes the value by which the Logistic model has periodic dynamics. This observation indeed implies some concrete frequency bands at which causal relations among the sub-time series of the periodic time series are stronger.
The time series generated from the Logistic model with different parameters are also characterized by the dynamics of their corresponding DANNs. As interestingly shown in figure 1, all the states of the neurons converge to the fixed points as λ corresponds to random or chaotic behavior in the noise-perturbed series y[q]; however, the states converge to periodic or quasi-periodic orbits, forming thin stripe patterns, as λ takes a value by which the Logistic model generates periodic dynamics. Hence, DANNs, when endowed with Cau-CNs, can preserve and even be used to differentiate some dynamical properties of time series for rather short data sets.
To validate the robustness of the proposed approach with regard to the length of the subtime series, the measures for the networks in the time domain are, respectively, plotted versus length in figure 1, where the curves become flatter as the length increases. This observation is consistent with the result, reported in [16], on the length requirement (at least 100 points) of the time series for computing the GC. Also this illustrates that our method is applicable for analyzing a rather short data set.
In addition, random phase surrogates and standard significance tests [17] are utilized to further validate the proposed approach. Specifically, for each λ, the generated time series are randomly permuted for 100 times. Surrogates of the measure APL are then computed based on the Cau-CNs constructed from these permuted time series. We denote the mean and standard deviation (SD) of these surrogates by m APL and σ APL , respectively. The statistic S APL is defined as S APL = (APL − m APL )/σ APL , which measures the deviation of APL from its surrogates. S AVS and S ACC are defined similarly. It can be seen from figure 1 that when λ corresponds to periodic dynamics, the three network measures are significantly different from their surrogates, while when λ corresponds to random or most of the chaotic time series, the difference, although not very outstanding, can still be identified. If we need to identify the difference significantly between random and chaotic dynamics, then in addition to the AR model some nonlinear models should also be utilized in the computation of the GC. Example 2 (Low-frequency versus high-frequency). Next we consider neural dynamics, described by the Rulkov model [18,19], with sinusoidal stimuli in two frequencies. We generate two time series {x i [q]} (i = 1, 2), corresponding to the membrane potential of the neuron, as follows: , Here, frequencies (see figure 2). These frequencies clearly correspond to either the intrinsic frequency of the chaotic bursting or the high frequency of the stimuli. More importantly, there is no overlap between the curves of APL and ACC across different frequencies in figure 2. In fact, these two quantities, rather than the AVS, reflect the global property of the underlying Cau-CNs. They therefore show more potential in differentiating the different time series in all frequencies.
Besides, after the first-order difference of the time series, that is, , the maxima positions of the measures do not change much (figure 2). This implies that the proposed approach could be suitable for characterizing some non-stationary time series, since the difference procedure is commonly used to render time series stationary and thus meet the requirements of classical Granger causal analysis.

Example 3 (Health versus arrhythmia)
. Now we analyze real-world data. Two MIT-BIH ECG data sets 8 are used in the analysis: one is from the Normal Sinus Rhythm Database (NSRD) and the other is from the Arrhythmia Database (AD). We select the first 18 subjects from each of the two data sets. ECG recorded the electrical activity of the heart and produced a time series for each subject covering about 20 min. 500 sub-time series are used to construct the corresponding Cau-CNs based on the above-proposed approach in the time domain. Here, after selecting a sub-time series that contains 500 points, we set the next R-peak as the start point of the following sub-time series and calculate the three network quantities mentioned above. For the NSRD, the means and SDs of the APL, AVS and ACC are 0.000 67 ± 0.000 49, 26.7 ± 9.5 and 0.015 ± 0.006, respectively, whereas those for the AD are 0.028 ± 0.020, 110.3 ± 32.5 and 0.081 ± 0.030, respectively. More striking differences are shown in figure 3, which characterizes the in-and out-strengths and CCs of the vertices in the constructed Cau-CNs. All these 9 differences suggest that stronger causal relations exist between the sub-time series for the arrhythmia subjects, but weaker causal relations for the normal sinus rhythm subjects. This also indicates that the heart rhythm of healthy people adapts more sensitively to internal or environmental fluctuations. These somewhat random fluctuations reasonably suppress the relations quantified by the GC.
To more profoundly study the difference between both data sets, we also establish and analyze the causal networks in the frequency domain. Here, we use 20 sub-time series, each of which contains 200 points. As shown in figure 3, the averages of the three quantities for the two data sets are different in all frequency bands. More importantly, the proposed approach enables us to figure out the frequency band (1-2 Hz) where the differences between the healthy subjects and arrhythmia patients are most prominent. This band interestingly corresponds to a human's normal heart rate, approximately 60-120 beats per minute. Besides, in figure 3, more overlaps of the SDs of the two datasets are observed for the AVS rather than the other two measures. The reason for this has been stressed in example 2.

Conclusion
Altogether, in this paper we have proposed a new approach to characterizing time series by combining the GC and the concept of complex networks. Directed and weighted complex networks constructed by the GC contain more details of the topological perspective of the interactions among sub-time series, and provide more possibilities to capture the various dynamics of the time series when compared with undirected and unweighted networks. More importantly, decomposition of the GC in the frequency domain enables the construction of frequency-specific complex networks, making it possible to characterize the time series at a particular frequency or frequency band. This is particularly useful in analyzing realworld time series when different frequency bands convey different signals and have distinct functions (e.g. local field potential or electroencephalography data) or when some frequencies are contaminated by noise and contain little information (e.g. functional magnetic-resonanceimaging data). Furthermore, due to the robustness of GC, the proposed approach can be also applied to non-stationary time series. By computing three representative network characteristics of the established causal networks, the proposed approach has been applied to characterize the time series produced by toy dynamical models and the recordings from two MIT-BIH data sets, in both the time and frequency domains. The results show the potential of our approach to capture various dynamics and to discriminate healthy subjects from arrhythmia patients.
Future works involve applying the proposed approach to analyze evolving time series by using a sliding window and constructing Cau-CNs, which may thus be applicable to rather short data sets. Moreover, since only the AR model is used for computing the GC, integrating the other nonlinear models [20,21] as well as more measures of complex networks into the present approach might provide a more accurate characterization, which also deserves further investigation.
for i = j. Hence, setting the transfer function matrix H ( f ) = {a i j ( f )} −1 2×2 leads the spectral model above to Then, the spectral density matrix S( f ) is given by S( f ) = H ( f ) H * ( f ). To eliminate the cross terms in S( f ), Geweke's normalization technique [12] is further adopted. Consequently, the GC from x j to x i at the frequency f is