Transfer Entropy Rate Through Lempel-Ziv Complexity

TL;DR

This paper introduces a practical, model-free method to estimate Transfer Entropy Rate using Lempel-Ziv complexity, applicable to discrete data, and effective in assessing information flow direction and strength.

Contribution

It presents a novel methodology for estimating Transfer Entropy Rate from single discrete series using Lempel-Ziv complexity, with advantages in simplicity and computational efficiency.

Findings

Effective in unidirectional coupled systems

Can be estimated from single discrete series

Assesses direction and strength of information flow

Abstract

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.

Equations

$${\mathbf{x}}_n^{\left(m\right)}=\left(x_n,x_{n-\tau },\dots ,x_{n-\left(m-1\right)\tau }\right),$$

$$p\left({\mathbf{x}}_n^{\left(m\right)}\right)=P\left\{\left(X_1,\dots ,X_m\right)={\mathbf{x}}_n^{\left(m\right)}\right\}.$$

$$h\left(X^{\left(m\right)}\right)$$

$$=H\left(X^{\left(m+1\right)}\right)-H\left(X^{\left(m\right)}\right),$$

$$H\left(X^{\left(m\right)}\right)$$

$$=-\sum _{x_1}\cdots \sum _{x_m}p\left({\mathbf{x}}_n^{\left(m\right)}\right)\ln p\left({\mathbf{x}}_n^{\left(m\right)}\right).$$

$$h\left(X\right)$$

$$=\underset{m\to \infty }{\lim }\frac{H\left(X^{\left(m\right)}\right)}{m}.$$

$$h\left(X\right)=\underset{T\to \infty }{\lim }\frac{{\mathcal{C}}_{x_t}\left[\ln \left(\alpha \right)+\ln {\mathcal{C}}_{x_t}\right]}{T}.$$

$$z_t=\sum _{i=1}^m{\alpha }^{i-1}{x}_{t,i},$$

$$h\left(X^{\left(m\right)}\right)$$

$$=\underset{T\to \infty }{\lim }\frac{{\mathcal{C}}_{z_t}\left[\ln \left({\alpha }^m\right)+\ln {\mathcal{C}}_{z_t}\right]}{T}.$$

$${\text{T}}_{Y\to X}^{\left(m\right)}=$$

$$+H\left(X_t,X_{t-\tau }^{\left(m\right)}\right)-H\left(X_{t-\tau }^{\left(m\right)}\right).$$

$${\text{t}}_{Y\to X}^{\left(m\right)}$$

$$=h\left(X_t,X_{t-\tau }^{\left(m\right)}\right)-h\left(Y_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right),$$

$$h\left(X|Y\right)$$

$$=\underset{m\to \infty }{\lim }\frac{H\left(Y_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right)}{m},$$

$$=h\left(Y_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right).$$

$${\text{tˆ}}_{\hat{Y}\to X}^{\left(m\right)}=-{⟨h_k\left({\hat{Y}}_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right)⟩}_K,$$

$$\mathcal{T}={\text{t}}_{Y\to X}^{\left(m\right)}-{\text{t}}_{X\to Y}^{\left(m\right)}-\left({\text{tˆ}}_{\hat{Y}\to X}^{\left(m\right)}-{\text{tˆ}}_{\hat{X}\to Y}^{\left(m\right)}\right).$$

$${\mathbf{V}}_n=\left\{\left(Y_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right)\right\},$$

$$n$$

$$t$$

$${\mathbf{v}}_n$$

$${\mathbf{x}}_{t-\tau }^{\left(m\right)}$$

$${\mathbf{y}}_{t-\tau }^{\left(m\right)}$$

$$h\left(Y_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right)$$

$$=h\left(Z\right),$$

$$=\underset{N\to \infty }{\lim }\frac{{\mathcal{C}}_{z_n}\left[\ln \left(2^{2m+1}\right)+\ln {\mathcal{C}}_{z_n}\right]}{N},$$

$$z_n=\mathcal{V}{\left[2^0,2^1,\dots ,2^{2m}\right]}^T,$$

$$\mathcal{V}=\left[Y_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right],$$

$${\mathcal{V}}_k=\left[{\hat{Y}}_{t-\tau }^{\left(m\right)},X_{t-\tau }^{\left(m\right)},X_t\right].$$

$$\begin{cases}y_{1}\left[n+1\right]&=1.4-y_{1}^{2}\left[n\right]-by_{2}\left[n\right],\\ y_{2}\left[n+1\right]&=y_{1}\left[n\right],\\ x_{1}\left[n+1\right]&=1.4-\Big{(}\epsilon y_{1}\left[n\right]+\left(1-\epsilon\right)x_{1}\left[n\right]\Big{)}x_{1}\left[n\right]-bx_{2}\left[n\right],\\ x_{2}\left[n+1\right]&=x_{1}\left[n\right],\\ \end{cases}$$

$$\{\begin{array}{ll}{\dot{y}}_1& =10\left(-y_1+y_2\right),\\ {\dot{y}}_2& ={\rho }_1{y}_1-y_2-y_1{y}_3,\\ {\dot{y}}_3& =y_1{y}_2-\frac{8}{3}{y}_3,\\ {\dot{x}}_1& =10\left(-x_1+x_2\right)+ϵ\left(y_1-x_1\right),\\ {\dot{x}}_2& ={\rho }_2{x}_1-x_2-x_1{x}_3,\\ {\dot{x}}_3& =x_1{x}_2-\frac{8}{3}{x}_3,\end{array}$$

$$\begin{cases}\dot{y}_{1}&=-\alpha\left(y_{2}+y_{3}\right),\\ \dot{y}_{2}&=\alpha\left(y_{1}+0.2y_{2}\right),\\ \dot{y}_{3}&=\alpha\Big{(}0.2+y_{3}\left(y_{1}-5.7\right)\Big{)},\\ \dot{x}_{1}&=10\left(-x_{1}+x_{2}\right),\\ \dot{x}_{2}&=28x_{1}-x_{2}-x_{1}x_{3}+\epsilon y^{\beta}_{2},\\ \dot{x}_{3}&=x_{1}x_{2}-\frac{8}{3}x_{3},\end{cases}$$

Full text

K****eywords Transfer Entropy, Transfer Entropy Rate, Lempel-Ziv Complexity.

1 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].

On the other hand, Lempel-Ziv’s complexity (LZC) is a classical measure that, for ergodic sources, relates the concepts of complexity (in the Kolmogorov-Chaitin sense), and entropy rate [7, 8]. For an ergodic dynamical process, the amount of new information gained per unit of time (entropy rate) can be estimated by measuring the capacity of this source to generate new patterns (LZC). Because of the simplicity of the LZC method, the entropy rate can be estimated from a single discrete sequence of measurements with a low computational cost [9].

In this article we aim to relate the concepts of Transfer Entropy Rate and Lempel-Ziv complexity. To be precise, we will exploit the advantages of the LZC methodology to calculate the TER between two ergodic dynamical systems.

The remainder of this paper is organized as follows. Section 2 begins with a brief review of the concepts of TE, TER and LZC. In Section 3 we described the proposed methodology to estimate the Transfer Entropy Rate through the Lempel-Ziv complexity. In Section 4 we present and analyze the results of the simulations carried out to evaluate the performance of our approach. Finally, Sections 5 and 6 the discussion and conclusions are presented.

2 Methodology

In this section we will briefly review some theoretical concepts related with the LZC, TE and TER. Moreover, we will introduce the notation used along the document.

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}\dots x_{T}}$. We can compose samples of an $m$-dimensional time-embedded process $\left\{X^{\left(m\right)}\right\}=\left\{X_{1},\dots,X_{m}\right\}$ by sampling $x_{t}$ with a frequency of $1/\tau$ [10, 11]:

[TABLE]

where $n=1,\dots,T-\left(m-1\right)\tau$. The process $\left\{X^{\left(m\right)}\right\}$ is characterized by the joint probability distribution:

[TABLE]

We can define the $m$-order entropy rate as [11]:

[TABLE]

where $H\left(X^{\left(m\right)}\right)$ is the entropy of the joint distribution $p\left(\bm{x}_{n}^{\left(m\right)}\right)$:

[TABLE]

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]:

[TABLE]

Equations (1) and (2) relate two different interpretations of the entropy rate. The first one tells that $h\left(X\right)$ 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.

2.1 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 $\left\{X_{t}\right\}$ 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 $\Omega_{x}$ of $\alpha$ symbols. 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 $\mathcal{C}_{x_{t}}$ 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\cdot 0\cdot 01\cdot 101\cdot 1100\cdot 1010\cdot 001011$. Then, the entropy rate can be computed as [8]:

[TABLE]

This approach can be easily generalized to multivariate processes by extending the alphabet size [7]. Consider an $m$-dimensional stationary process $\left\{X^{\left(m\right)}\right\}$, that produces the sequences $x_{t,i}=x_{1,i},\dots,x_{T,i}$ with ${i=1,\dots,m}$, each one of them from an alphabet of $\alpha$ symbols. Let ${z_{t}=z_{1},\dots,z_{T}}$ be a new sequence defined over an extended alphabet of size $\alpha^{m}$ [7]:

[TABLE]

then the joint Lempel-Ziv complexity ${\mathcal{C}_{x_{t,i}}=\mathcal{C}_{z_{t}}}$ and the $m$-order entropy rate can be calculated as [7, 8]:

[TABLE]

2.2 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 [5, 6, 15, 1]:

[TABLE]

The parameter $m$ is commonly called the history length or embedding dimension and $\tau$ is the lag or embedding lag [6]. $\text{{T}}_{Y\to X}^{\left(m\right)}$ quantifies the amount of information contained in the $m$-past states of process $Y$ (${Y}_{t-\tau}^{\left(m\right)}$) about the current state of the process $X$ ($X_{t}$), that is not already explained by the $m$-past states of process $X$ (${X}_{t-\tau}^{\left(m\right)}$). This measure is asymmetric ($\text{{T}}_{Y\to X}^{\left(m\right)}\neq\text{{T}}_{X\to Y}^{\left(m\right)}$) and increases along with the coupling level, allowing to determine the direction and strength of the information flow [5].

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 Transfer Entropy Rate [3, 4, 16]:

[TABLE]

where $h\left(X\right)$ is the entropy rate of $X$ and $h\left(X|Y\right)$ is the conditional entropy rate [16]:

[TABLE]

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 $\text{{t}}_{Y\to X}^{\left(m\right)}$ should be equal to zero. However, in practical applications, the estimation of $\text{{t}}_{Y\to X}^{\left(m\right)}$ 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 measurement $\text{{\^{t}}}_{\widehat{Y}\to X}^{\left(m\right)}$. 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 $Y$ to $X$ is defined as:

[TABLE]

where ${\hat{Y}}_{t-\tau}^{\left(m\right)}$ is obtained by redrawing with replacement samples from ${Y}_{t-\tau}^{\left(m\right)}$, and $<\cdot>_{K}$ is the mean value over the $k=1,2,\dots,K$ surrogate realizations.

In order to assess the directionality of the information transport we need to analyze the global TER estimator:

[TABLE]

A positive value of $\mathcal{T}$ suggests that the information flow goes from system $Y$ to system $X$, meanwhile a negative value suggests the contrary. Finally, if $\mathcal{T}=0$, there is no information flow between systems.

3 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 delayed embedding vectors from quantized time series. For simplicity in the description of the method we will assume binary time series, although this methodology can be extended to higher quantization levels111For example, a time series can be quantized into $\alpha$ symbols using as thresholds its own $\alpha$-quantiles.. Consider two binarized time series ($\alpha=2$) from a coupled system: let $x_{t}=x_{1}\dots x_{T}$ be the target and $y_{t}=y_{1}\dots y_{T}$ be the source. In all the simulations, each time series was binarized with its own median value. Set the parameter $m$ (embedding dimension) and $\tau$ (embedding lag) and create the collection of embedding vectors (see Fig. 1):

[TABLE]

where:

[TABLE]

By construction $\bm{V}_{n}$ is a collection of ${\left(2m+1\right)}\text{-uples}$ and we can define the sequence ${z_{n}=\sum_{i=1}^{2m+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]:

[TABLE]

where $\mathcal{C}_{z_{n}}$ 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 ${\bm{V}_{n}=\left\{\left({X}_{t-\tau}^{\left(m\right)},X_{t}\right)\right\}}$ and the sequence ${z_{n}=\sum_{i=1}^{m+1}{2^{i-1}v_{n,i}}}$. Moreover, we can use the same methodology to obtain a surrogate measurement $h_{k}\left({\hat{Y}}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right)$. In this case, ${\hat{Y}}_{t-\tau}^{\left(m\right)}$ is obtained by shuffling (or redrawing with replacement) ${Y}_{t-\tau}^{\left(m\right)}$ amongst the set of $\left\{{Y}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right\}$ tuples.

3.1 Implementation

As we have mentioned before, our 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 ($\tau$). 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 $\tau$ 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 $\tau=\text{max}\left(\tau_{x},\tau_{y}\right)$, where $\tau_{x}$ and $\tau_{y}$ are the lags that minimize the mutual information function between $x_{t}$ and $x_{t-\tau}$, and between $y_{t}$ and $y_{t-\tau}$, respectively.

The Algorithm 1 describes the steps to calculate the global Transfer Entropy Rate ($\mathcal{T}$), between processes $X$ and $Y$. In the first step both time series, $x_{t}=x_{1}\dots x_{T}$ and $y_{t}=y_{1}\dots 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$ ($\text{{t}}_{Y\to X}^{\left(m\right)}$) and step 7 to obtain its surrogate estimation ($\text{{\^{t}}}_{\widehat{Y}\to X}^{\left(m\right)}$). As it was previously described, to estimate each term in equation (5) we need to embed the binarized time series in spaces with different dimensions. Considering the one with the greatest dimension: $\left\{{Y}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right\}$, the embedding vectors in this space can be disposed in a matrix ${\mathcal{V}=\left[{Y}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right]}$ (see figure 1). The sequence $z_{n}$ can be expressed as the product:

[TABLE]

and, the entropy rate $h\left({Y}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t},\right)$ can be computed using equation (8). Regarding the second term in equation (5), given that the space $\left\{{X}_{t-\tau}^{\left(m\right)},X_{t}\right\}$ is a subspace of $\left\{{Y}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right\}$, the sequence $z_{n}$ needed to estimate $h\left({X}_{t-\tau}^{\left(m\right)},X_{t}\right)$ can be computed by taking the product of the last $m+1$ columns of $\mathcal{V}$ with $\left[2^{0},2^{1},\dots,2^{m}\right]^{T}$. Finally, use equation (5) to obtain $\text{{t}}_{Y\to X}^{\left(m\right)}$.

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 $\mathcal{V}_{k}=\left[{\hat{Y}}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right]$ for $k=1,2,\dots,K$. Here, for each $k$, the sample $\left\{{\hat{Y}}_{t-\tau}^{\left(m\right)}\right\}$ is set by redrawing with repetition from the collection $\left\{{Y}_{t-\tau}^{\left(m\right)}\right\}$. Then, the surrogate Transfer Entropy Rate $\text{{\^{t}}}_{\widehat{Y}\to X}^{\left(m\right)}$ is obtained as the mean value of $h_{k}\left({\hat{Y}}_{t-\tau}^{\left(m\right)},{X}_{t-\tau}^{\left(m\right)},X_{t}\right)$ over all surrogate realizations.

To estimate $\text{{t}}_{X\to Y}^{\left(m\right)}$ (transfer entropy from $X$ to $Y$) and its surrogate $\text{{\^{t}}}_{\widehat{X}\to Y}^{\left(m\right)}$ 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).

4 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$222The minimum embedding dimension for the Henon-Henon system is $m_{x}+m_{y}=4$ and, for the Lorenz-Lorenz and Rössler-Lorenz is $m_{x}+m_{y}=6$. and $\tau$ that met the conditions mentioned in Subsection 3.1.

The first system was the coupled Henon-Henon [18, 2]:

[TABLE]

where $b=0.3$. For the simulation, the coupling parameter $\epsilon$ was varied from zero to one in steps of $0.1$. For each $\epsilon$, $200$ realizations were computed using random initial conditions. The Transfer Entropy Rate was calculated with $m=\left\{2,3,4,5,6,7\right\}$ and $\tau=\left\{1,3,5,7,10\right\}$. This procedure was repeated for data lengths $N=\left\{3000,5000,10000\right\}$.

The results are shown in the Fig. 2. Each plot presents the global TER ($\mathcal{T}$), calculated with $m=5$ and $\tau=1$, as a function of the coupling parameter $\epsilon$. It can be observed in Fig. 2a ($N=3000$) that the median value of estimator $\mathcal{T}$ is zero for $\epsilon=0$. This is an expected result since there is no information flow between the two systems. Moreover, the median value of $\mathcal{T}$ increases along with the coupling parameter until $\epsilon=0.5$. The positivity of $\mathcal{T}$ points out the correct direction of coupling and its increasing magnitude indicates a rising strength of the coupling. On the other hand, for $\epsilon\geq 0.7$ the median value of $\mathcal{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, $\mathcal{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 $\mathcal{T}$ decreases as long as the data length is increased.

For the second simulation we have chosen the Lorenz-Lorenz system:

[TABLE]

where $\rho_{1}=28.5$, $\rho_{2}=27.5$ and $\epsilon\in\left\{0,\dots,15\right\}$. For each value of the coupling parameter, $200$ realization were computed, each one starting from a different initial condition. The numerical integration was performed using the ode45 function of Matlab (algorithm of Dormand and Prince) with step size $\Delta t=0.03$. For each realization the first $10000$ data points where discarded. Then, $\mathcal{T}$ was calculated for all the combination of the parameters: $m=\left\{2,3,4,5,6,7\right\}$ and $\tau=\left\{1,3,5,7,10\right\}$. The above procedure was applied varying the data length $N=\left\{3000,5000,10000\right\}$.

In Fig. 3 is shown estimator $\mathcal{T}$ ($m=7$ and $\tau=5$) as a function of $\epsilon$ for the Lorenz-Lorenz coupled system. For $N=3000$ (Fig. 3a) it can be observed that for uncoupled systems $\epsilon=0$ the median value of $\mathcal{T}\approx 0$. As the coupling parameter increases, $\mathcal{T}$ is positive and grows until the synchronization threshold is reached $\epsilon\approx 12$ [18]. From this point, the value of $\mathcal{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 $\mathcal{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:

[TABLE]

where $\alpha=6$, $\beta=2$ and $\epsilon\in\left\{0,0.2,\dots,5\right\}$. 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 $\Delta t=0.02617$ [15]. The TER was calculated using the same parameter’s values of the afore simulations.

In Fig. 4 the behavior of $\mathcal{T}$ as a function of the coupling parameter for $m=7$ and $\tau=10$ is shown. For this coupled system the synchronization threshold is $\epsilon\approx 2.8$ [18, 20]. It can be observed in Fig. 4a ($N=3000$) that the median value of $\mathcal{T}$ is always positive, even for $\epsilon=0$. This means that the $\mathcal{T}$ estimator is detecting false coupling for $\epsilon=0$. This phenomenon has been also observed in the Symbolic Transfer Entropy [19]. However, the $\mathcal{T}$ estimator is detecting the correct coupling direction. Figs. 4b and 4c display a similar behaviour but notice that the variance of $\mathcal{T}$ decreases.

5 Discussion

There are two methodologies to estimate the Transfer Entropy Rate that are similar to our approach. The first one is the Symbolic Transfer Entropy [19], which finds its foundations in the Permutation Entropy [21]. The second one can be found in [22] and is based on the K-nearest-neighbor (KNN) estimation method proposed by Kraskov et al. [23]. In order to compare our methodology with the ones mentioned above, we have implemented both algorithms and calculated the global TER and the computation time of each one of the $200$ realizations of the coupled Lorenz-Lorenz system. The simulation was done using $m=\left\{2,3,4,5,6\right\}$, $N=\left\{1000,3000,5000,10000\right\}$, $\epsilon=\left\{0,\dots,15\right\}$, $\tau=5$ and $K=30$ surrogate realizations (just for our approach). The simulation was performed in a cluster with 10 nodes, each one has 2 Intel Xeon E5-2670 v3 2.5GHz processors of 12 Cores.

In Fig. 5 we show the global TER as a function of the coupling parameter using the three methods. It can be observed that $\mathcal{T}$ presents a very similar behaviour for all the three methods. However, it is important to notice that for a given $\epsilon$ their values differ. This result strength the hypothesis that our methodology can be used as a information transfer measure.

Results about the comparison of computation times among methodologies can be found in Figs. 6 and 7, where we present a boxplot of the computation time of a single realization as a function of the data lenght and the embedding dimension, respectively. As it can be seen in Figs. 6a, 6b and 6c, the computation time of a single realization, increases exponentially as the data length increases, regardless the employed methodology. Moreover, it can be observed that, for $m=4$, the fastest method is the Symbolic Transfer entropy, followed by our methodology and far away is the KNN approach. On the other hand, Figs. 7a, 7b and 7c show the execution time for a single realization for different embedding dimensions. Notice that the computational cost of each method increases with $m$ in different ways. For our approach the increasing is exponential since the $m$ parameter is linked to the alphabet size in the Lempel-Ziv’s algorithm. We must point out that our methodology outperforms the Symbolic Transfer entropy approach for $m>5$ (compare the computation time in Figs 7a and 7b for $m=6$). This is because the computational cost of the Symbolic Transfer entropy increases with $m$ in a factorial way. This suggests that our methodology has an advantage over the Symbolic Transfer entropy when the analysis of high-dimensional systems is needed. Finally, note for the greatest embedding dimension here studied ($m=8$) and the longest data length ($N=10000$) tested in this simulation, our approach performs a single Transfer Entropy Rate estimation in less than three seconds.

Based on the results, we can conclude that the estimator $\mathcal{T}$ here proposed (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 $\tau$ (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 $\tau$ ensure good reconstruction of the state space.

In future studies, we will address the implementation of our methodology using different embedding parameters for $x_{t}$ and $y_{t}$ as well as a different alphabet size.

6 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.