Transfer Entropy Rate Through Lempel-Ziv Complexity

In this article we present a methodology to estimate the Transfer Entropy Rate between two systems through the Lempel-Ziv complexity. This methodology carries a set of practical advantages: it can be estimated from two single discrete series of measures, it is not computationally expensive and it does not assume any model for the data. The results of simulations over three different unidirectional coupled systems, suggest that this methodology can be used to assess the direction and strength of the information flow between systems.


Introduction
The Transfer Entropy (TE) and the Transfer Entropy Rate (TER) are closely related concepts that measure information transport. The former was proposed by Schreiber in [1] and independently by Paluš in [2]. The later was described by Amblard et al. in [3,4]. They are able to quantify the strength and direction of the coupling between simultaneously observed systems [5]. Moreover, they have become of general interest since they can be used to study complex interaction phenomena found in many disciplines [6].
Since our intent is to investigate a possible causality connection between two dynamical systems, we need to analyze the signals that they produce. We will assume the existence of ergodic probability measures that describe the density of trajectories in phase space, such us it can be treated as probability densities. This allows us to analyze the dynamics of the systems through the construction of random processes from their signals.
Consider a system X that produces a time series x t = x 1 . . . x T . We can compose samples of an m-dimensional time-embedded process X (m) = {X 1 , . . . , X m } by sampling x t with a frequency of 1/τ [10,11]: We can define the m-order entropy rate as [11]: where H X (m) is the entropy of the joint distribution p x The m-order entropy rate measures the variation of the total information in the time-embedded process when the embedding dimension m is increased by 1. From this definition we can calculate the entropy rate of the system X as [12,13]: Equations (1) and (2) relate two different interpretations of the entropy rate. The first one tells that h (X) is a measure of our uncertainty about the present state of the system under the assumption that its entire past is observed. The second one states that the entropy rate is the average information gained by observing the system. In this respect, systems with a higher entropy rates generate information at a higher rate, and this make their dynamics more difficult to predict.

Lempel-Ziv Complexity
The concepts of entropy rate and Lempel-Ziv complexity are closely related since systems with higher entropy rate tend to generate more complex sequences (time series). In that context, the entropy rate of an ergodic system can be estimated by measuring its capacity to generate new patterns [9]. Estimating the entropy rate of a system using the Lempel-Ziv algorithm carries a set of practical advantages: it can be estimated from a single discrete series of measures, the algorithm is fast and it does not assume any model for the data.
Suppose a stationary stochastic process {X t } that produces a sequence x t of length T , where for a fixed t, the random variable X t can take values from an alphabet Ω x of size α. To estimate the complexity of this process we will use the Lempel and Ziv's scheme proposed in 1976 [14]. In this approach, a sequence x t is parsed into a number C xt of words, by considering as a new word any subsequence that has not yet been encountered. For example the sequence 100110111001010001011 is parsed in 7 words: 1 · 0 · 01 · 101 · 1100 · 1010 · 001011. Then, the entropy rate can be computed as [8]: This approach can be easily generalized to multivariate processes by extending the alphabet size [7]. Consider an m-dimensional binary stationary process X (m) , that produces the sequences x t,i = x 1,i , . . . , x T,i with i = 1, . . . , m. Let z t = z 1 , . . . , z T be a new sequence defined over an extended alphabet of size α = 2 m [7]: then the joint Lempel-Ziv complexity C xt,i = C zt and the m-order entropy rate can be calculated as [7,8]:

Transfer Entropy and Transfer Entropy Rate
The Transfer Entropy is able to assess the amount of information transferred from process Y (driver/source) to process X (driven/target). It is defined as [1,5,6,15]: The parameter m is commonly called the history length or embedding dimension and τ is the lag or embedding lag [6]. T t−τ ) about the current state of the process X (X t ), that is not already explained by the m-past states of process X (X In [3] Amblard et al. suggest that under stationarity conditions the TE can be considered as an information flow rate. This idea leads to the definition of TER [3,4,16]: where h (X) is the entropy rate of X and h (X|Y ) is the conditional entropy rate [16]:

Transfer Entropy Rate Through Lempel-Ziv Complexity
A preprint  Figure 1: Diagram to obtain the sequence z n to calculate the entropy rate h Y The TER lies between zero and the entropy rate of the target X, being equal to zero if X and Y are independent [16].
If the processes X and Y had no relationship then t (m) Y →X should be equal to zero. However, in practical applications, the estimation of t (m) Y →X could present a bias due to the finite length of the data. Some authors [6] proposed to correct this bias by empirically finding the distribution of the surrogate measurementt (m) Y →X . The surrogate data must be generated in such a way that the temporal correlation between the source and the target is destroyed but statistical properties and the temporal structure of both processes are preserved [5,6]. Note that only the second term in equation (5) depends on the source, so the surrogate transfer entropy from whereŶ t−τ , and < · > K is the mean value over the k = 1, 2, . . . , K surrogate realizations.
In order to assess the directionality of the information transport we need to analyze the global TER estimator: A positive value of T suggests that the information flow goes from system Y to system X, meanwhile a negative value suggests the contrary. Finally, if T = 0, there is no information flow between systems.

Transfer Entropy Rate Based on Lempel-Ziv Complexity
In this section we formalize our approach to estimate the Transfer Entropy Rate using the Lempel-Ziv complexity. The idea is to estimate the two joint entropy rates on the right-hand side of equation (5) by means of their associated joint Lempel-Ziv complexities. To this end, we propose a methodology based on the construction of binary delayed embedding vectors from time series. Consider two binarized time series Transfer Entropy Rate Through Lempel-Ziv Complexity A preprint from a coupled system: let x t = x 1 . . . x T be the target and y t = y 1 . . . y T be the source. In this document each time series is binarized with its own median value. Set the parameter m (embedding dimension) and τ (embedding lag) and create the collection of embedding vectors (see Fig. 1): By construction V n is a collection of (2m + 1)-uples and we can define the sequence z n = 2m+1 i=1 2 i−1 v n,i , over an extended alphabet of size α = 2 2m+1 (see Fig. 1). Then, the joint entropy rate can be calculated as [7]: where C zn is the LZC of the sequence z n . This procedure can be followed to estimate the first term of equation (5), considering the collection of embedding vectors V n = X t−τ , X t tuples.

Implementation
As we have mentioned before, the here proposed methodology is based on the construction of embedding spaces from time series. In this direction, our algorithm has two parameters: the embedding dimension (m) and the embedding lag (τ ). As well as other embedding based algorithms, we have found that the best results are achieved when a good reconstruction of the state space is guaranteed [17]. In other words, when m is bigger than the minimum embedding dimension of the system and τ is large enough so that the various coordinates of the embedding vectors contains as much new information as possible, without being entirely independent. In this sense a good choice of embedding dimension is m = m x + m y + 1, here m x and m y are estimations of the minimum embedding dimension of X and Y , respectively. On the other hand, we propose to use an embedding lag value τ = max (τ x , τ y ), where τ x and τ y are the lags that minimize the mutual information function between x t and x t−τ , and between y t and y t−τ , respectively.
The Algorithm 1 describes the steps to calculate the global Transfer Entropy Rate (T ), between processes X and Y . In the first step both time series, x t = x 1 . . . x T and y t = y 1 . . . y T , must be binarized. This can be done using measures as the mean value, or the median value, among other options. Consider x t as the target/driven series and y t as the source/driver series and follow steps 3-6 to obtain the TER from Y to X (t (m) Y →X ) and step 7 to obtain its surrogate estimation (t (m) Y →X ). As it was previously described, to estimate each term in equation (5)  t−τ , X t , can be computed using equation (8). Regarding the second term in equation (5), given that the space X t−τ , X t , the sequence z n needed to estimate h X (m) t−τ , X t can be computed by taking the product of the last m + 1 columns of V with 2 0 , 2 1 , . . . , 2 m T . Finally, use equation (5) to obtain t (m) Y →X . In the step 7 we must generate the surrogate data. Choose a number K of surrogate data sets to be generated and build the surrogate data matrices V k = Ŷ (m) t−τ , X X→Y just set y t as the target series, x t as the source and repeat the procedure described in the steps 3-8. Finally, the global TER must be estimated using equation (7). Algorithm 1 LeZTER Algorithm. MATLAB code: https://bitbucket.org/jrinckoar/tentropyrate-lzc/src/master/ 1: Binarize the temporal series x t and y t using each median value. 2: Set x t as target/driven series and y t as source/driver series.
3: Given a value of m and τ , set the matrix of embedding vectors: t−τ , X t , and obtain the sequence z n (see Fig. 1).

Results
We have conducted three simulations using different unidirectional coupled systems: the Henon-Henon, the Lorenz-Lorenz and the Lorenz driven by Rössler. The results presented in Figs. 2, 3 and 4 were computed using values of m † and τ that met the conditions mentioned in Subsection 3.1. † The minimum embedding dimension for the Henon-Henon system is mx + my = 4 and, for the Lorenz-Lorenz and Rössler-Lorenz is mx + my = 6.

Transfer Entropy Rate Through Lempel-Ziv Complexity
A preprint The first system was the coupled Henon-Henon [2,18]: The results are shown in the Fig. 2. Each plot presents the global TER (T ), calculated with m = 5 and τ = 1, as a function of the coupling parameter . It can be observed in Fig. 2a (N = 2000) that the median value of estimator T is zero for = 0. This is an expected result since there is no information flow between the two systems. Moreover, the median value of T increases along with the coupling parameter until = 0.6. The positivity of T points out the correct direction of coupling and its increasing magnitude indicates a rising strength of the coupling. On the other hand, for ≥ 0.7 the median value of T is zero. For these values of the coupling parameter the Henon-Henon system is synchronized in such a way that both systems are statistically indistinguishable. In this kind of situations, T is unable to point out any information flow. This behaviour has been already observed on other transfer entropy estimators [2,15,18]. It can be seen in Figs. 2b and 2c (N = 5000 and N = 10000, respectively) that the variance of T decreases as long as the data length is increased.
Transfer Entropy Rate Through Lempel-Ziv Complexity For the second simulation we have chosen the Lorenz-Lorenz system:  (Fig. 3a) it can be observed that for uncoupled systems = 0 the median value of T ≈ 0. As the coupling parameter increases, T is positive and grows until the synchronization threshold is reached ≈ 12 [18]. From this point, the value of T goes toward zero despite the systems are coupled. The same behavior was displayed by the symbolic transfer entropy [19] and other information flow estimators calculated over the same system [18]. As well as in the case of the Henon-Henon coupled system, the variance of T decreases as the number of data points is increased (figures 3b and 3c).
The third system is the Lorenz driven by Rössler (Rössler-Lorenz) system: where α = 6, β = 2 and ∈ {0, 0.2, . . . , 5}. For each value of the coupling parameter, 200 realization using random initial conditions were computed. The numerical integration was performed using the same methodology described for the Lorenz-Lorenz system, but ∆t = 0.03 [19]. The TER was calculated using the same parameter's values of the afore simulations.
In Fig. 4 the behavior of T as a function of the coupling parameter for m = 7 and τ = 10 is shown. For this coupled system the synchronization threshold is ≈ 2.8 [18,20]. It can be observed in Fig. 4a (N = 2000) that the median value of T is always positive, even for = 0. This means that the T estimator is detecting false coupling for = 0. This phenomenon has been also observed in the Symbolic Transfer Entropy [19]. However, the T estimator is detecting the correct coupling direction. Figs. 4b and 4c display a similar behaviour but notice that the variance of T decreases.

Discussion
Based on the results shown in the last Section, we can conclude that the estimator here proposed T (equation (7)) is able to detect the direction and strength of the information flow between two coupled ergodic systems. This methodology, based in the LZC, is computationally fast and it does not assume any model for the data. Our results are comparable with those obtained by the Symbolic Transfer Entropy [19] and with the ones reported by Krakovská et al. in [18].
The TER depends on two parameters: m (the history length) and τ (the lag). For simplicity, we considered that these parameters should be the same for the embedding of both time series, although they can be different [6]. We have observed that the best results are achieved when the values of m and τ ensure good reconstruction of the state space.
In future studies, we will address the implementation of the of our methodology using different embedding parameters for x t and y t .

Conclusions
In this article we have presented a new methodology to calculate the Transfer Entropy Rate between two systems based on the Lempel-Ziv's complexity. Because of the properties of the Lempel-Ziv's algorithm, we were able to propose a computationally fast methodology to estimate the information flow between two systems. This methodology have been assessed using three unidirectional coupled systems: Henon-Henon system, Lorenz-Lorenz system and Rössler-Lorenz system. The results show that our estimator is able to detect the direction and strength of the information flow.