Transfer Entropy for Coupled Autoregressive Processes

Abstract

A method is shown for computing transfer entropy over multiple time lags for coupled autoregressive processes using formulas for the differential entropy of multivariate Gaussian processes. Two examples are provided: (1) a first-order filtered noise process whose state is measured with additive noise, and (2) two first-order coupled processes each of which is driven by white process noise. We found that, for the first example, increasing the first-order AR coefficient while keeping the correlation coefficient between filtered and measured process fixed, transfer entropy increased since the entropy of the measured process was itself increased. For the second example, the minimum correlation coefficient occurs when the process noise variances match. It was seen that matching of these variances results in minimum information flow, expressed as the sum of transfer entropies in both directions. Without a match, the transfer entropy is larger in the direction away from the process having the larger process noise. Fixing the process noise variances, transfer entropies in both directions increase with the coupling strength. Finally, we note that the method can be generally employed to compute other information theoretic quantities as well.

1. Introduction

Transfer entropy [1] quantifies the information flow between two processes. Information is defined to be flowing from system X to system Y whenever knowing the past states of X reduces the uncertainty of one or more of the current states of Y above and beyond what uncertainty reduction is achieved by only knowing the past Y states. Transfer entropy is the mutual information between the current state of system Y and one or more past states of system X, conditioned on one or more past states of system Y. We will employ the following notation. Assume that data from two systems X and Y are simultaneously available at k timestamps: $t_{n-k+2:n+1}\equiv \left\{t_{n-k+2},t_{n-k+2},\mathrm{...},t_n,t_{n+1}\right\}$. Then we express transfer entropies as:

$$T{E}_{x\to y}^{\left(k\right)}=I\left(y_{n+1};x_{n-k+2:n}|y_{n-k+2:n}\right)=H\left(y_{n+1}|y_{n-k+2:n}\right)-H\left(y_{n+1}|y_{n-k+2:n},x_{n-k+2:n}\right)$$

(1)

$$T{E}_{y\to x}^{\left(k\right)}=I\left(x_{n+1};y_{n-k+2:n}|x_{n-k+2:n}\right)=H\left(x_{n+1}|x_{n-k+2:n}\right)-H\left(x_{n+1}|x_{n-k+2:n},y_{n-k+2:n}\right).$$

(2)

Each of the two transfer entropy values TE_x→y and TE_y→x is nonnegative and both will be positive (and not necessarily equal) when information flow is bi-directional. Because of these properties, transfer entropy is useful for detecting causal relationships between systems generating measurement time series. Indeed, transfer entropy has been shown to be equivalent, for Gaussian variables, to Granger causality [2]. Reasons for caution about making causal inferences in some situations using transfer entropy, however, are discussed in [3,4,5,6]. A formula for normalized transfer entropy is provided in [7].

The contribution of this paper is to explicitly show how to compute transfer entropy over a variable number of time lags for autoregressive (AR) processes driven by Gaussian noise and to gain insight into the meaning of transfer entropy in such processes by way of two example systems: (1) a first-order AR process X = {x_n} with its noisy measurement process Y = {y_n}, and (2) a set of two mutually-coupled AR processes. Computation of transfer entropies for these systems is a worthwhile demonstration since they are simple models that admit intuitive understanding. In what follows we first show how to compute the covariance matrix for successive iterates of the example AR processes and then use these matrices to compute transfer entropy quantities based on the differential entropy expression for multivariate Gaussian random variables. Plots of transfer entropies versus various system parameters are provided to illustrate various relationships of interest.

Note that Kaiser and Schreiber [8] have previously shown how to compute information transfer metrics for continuous-time processes. In their paper they provide an explicit example, computing transfer entropy for two linear stochastic processes where one of the processes is autonomous and the other is coupled to it. To perform the calculation for the Gaussian processes the authors utilize expressions for the differential entropy of multivariate Gaussian noise. In our work, we add to this understanding by showing how to compute these quantities analytically for higher time lags. We now provide a discussion of differential entropy, the formulation of entropy appropriate to continuous-valued processes as we are considering.

2. Differential Entropy

The entropy of a continuous-valued process is given by its differential entropy. Recall that the entropy of a discrete-valued random variable is given by the Shannon entropy $H=-{\displaystyle \sum _i{p}_i\log }{p}_i$ (we shall always choose log base 2 so that entropy will be expressed in units of bits) where p_i is the probability of the i^th outcome and the sum is over all possible outcomes.

Following [9] we derive the appropriate expression for differential entropies for conditioned and unconditioned continuous-valued random variables. When a process X is continuous-valued we may approximate it as a discrete-value process by identifying p_i = f_iΔx where f_i is the value of the pdf at the i^th partition point and Δx is the refinement of the partition. We then obtain:

$$\begin{array}{l}H\left(X\right)=-{\displaystyle \sum _i{p}_i\log }{p}_i\\ =-{\displaystyle \sum _i{f}_i\Delta x\log }{f}_i\Delta x\\ =-{\displaystyle \sum _i{f}_i\Delta x\left(\log f_i+\log \Delta x\right)}\\ =-{\displaystyle \sum _i{f}_i\log f_i\Delta x}-{\displaystyle \sum _i\log \Delta x{f}_i\Delta x}\\ =-{\displaystyle \int f\log fdx}-\log \Delta x{\displaystyle \int fdx}\\ =h\left(X\right)-\log \Delta x\end{array}$$

(3)

Note that since the X process is continuous-valued, then, as Δx → 0, we have H(X) → + infinity. Thus, for continuous-valued processes, the quantity h(X), when itself defined and finite, is used to represent the entropy of the process. This quantity is known as the differential entropy of random process X.

Closed-form expressions for the differential entropy of many distributions are known. For our purposes, the key expression is the one for the (unconditional) multivariate normal distribution [10]. Let the probability density function of the n-dimensional random vector x be denoted f(x), then the relevant expressions are:

$$\begin{array}{l}f\left(\overline{x}\right)=\frac{\exp \left[-\frac{1}{2}{\left(\overline{x}-\overline{\mu }\right)}^T{C}^{-1}\left(\overline{x}-\overline{\mu }\right)\right]}{{\left(2\pi \right)}^{\frac{n}{2}}{\left[\det C\right]}^{\frac{1}{2}}}\\ h\left(\overline{x}\right)=-{\displaystyle \int f\left(\overline{x}\right)\log \left[f\left(\overline{x}\right)\right]dx}\\ =\frac{1}{2}\log \left[{\left(2\pi e\right)}^n\det C\right]\end{array}$$

(4)

where detC is the determinant of matrix C, the covariance of x. In what follows, this expression will be used to compute differential entropy of unconditional and conditional normal probability density functions. The case for conditional density functions warrants a little more discussion.

Recall that the relationships between the joint and conditional covariance matrices, C_XY and C_Y|X, respectively, of two random variables X and Y (having dimensions n_x and n_y, respectively) are given by:

$$\begin{array}{l}{C}_{XY}=\mathrm{cov}\left(\left[\begin{array}{c}X\\ Y\end{array}\right]\right)=\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]\\ \\ \mathrm{cov}\left[Y|X=x\right]=C_{Y|X}={\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{12}.\end{array}$$

(5)

Here blocks Σ₁₁ and Σ₂₂ have dimensions n_x by n_x and n_y by n_y, respectively. Now, using Leibniz’s formula, we have that:

$$\det C_{XY}=\det \left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]=\det {\Sigma }_{11}\det \left({\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{12}\right)=\det C_X\det C_{Y|X}.$$

(6)

Hence the conditional differential entropy of Y, given X, may be conveniently computed using:

$$\begin{array}{l}h\left(Y|X\right)=\frac{1}{2}\log \left[{\left(2\pi e\right)}^{n_y}\det C_{Y|X}\right]\\ =\frac{1}{2}\log \left[{\left(2\pi e\right)}^{n_y}\frac{\det C_{XY}}{\det C_X}\right]\\ =\frac{1}{2}\log \left[{\left(2\pi e\right)}^{n_x+n_y}\det C_{XY}\right]-\frac{1}{2}\log \left[{\left(2\pi e\right)}^{n_x}\det C_X\right]\\ =h\left(X,Y\right)-h\left(X\right).\end{array}$$

(7)

This formulation is very handy as it allows us to compute many information-theoretic quantities with ease. The strategy is as follows. We define C^(k) to be the covariance of two random processes sampled at k consecutive timestamps {t_n−k+2, t_{n− k+1}, …, t_n, t_n+1}. We then compute transfer entropies for values of k up to k sufficiently large to ensure that their valuations do not change significantly if k is further increased. For our examples, we have found k = 10 to be more than sufficient. A discussion of the importance of considering this sufficiency is provided in [11].

3. Transfer Entropy Computation Using Variable Number of Timestamps

We wish to consider two example processes each of which conforms to one of the two model systems having the general expressions:

$$(1)\{\begin{array}{c}{x}_{n+1}=a_0{x}_n+a_1{x}_{n-1}+\cdots +a_m{x}_{n-m}+w_n\\ y_{n+1}=c_{-1}{x}_{n+1}+v_n\\ v_n~N\left(0,R\right),w_n~N\left(0,Q\right)\end{array}$$

(8)

and:

$$(2)\{\begin{array}{c}{x}_{n+1}=a_0{x}_n+a_1{x}_{n-1}+\cdots +a_m{x}_{n-m}+b_0{y}_n+b_1{y}_{n-1}+\cdots +b_j{y}_{n-j}+w_n\\ y_{n+1}=c_0{x}_n+c_1{x}_{n-1}+\cdots +c_m{x}_{n-m}+d_0{y}_n+d_1{y}_{n-1}+\cdots +d_j{y}_{n-j}+v_n\\ v_n~N\left(0,R\right),w_n~N\left(0,Q\right).\end{array}$$

(9)

Here, v_n and w_n are zero mean uncorrelated Gaussian noise processes having variances R and Q, respectively. For system stability, we require the model poles to lie within the unit circle.The first model is of a filtered process noise X one-way coupled to an instantaneous, but noisy measurement process Y. The second model is a two-way coupled pair of processes, X and Y.

Transfer entropy (as defined by Schreiber [1]) considers the flow of information from past states (i.e., state values having, timetags $t_{n-k+2:n}\equiv \left\{t_{n-k+2},t_{n-k+2},\mathrm{...},t_n\right\}$) of one process to the present $\left(t_{n+1}\right)$ state of another process. However, note that in the first general model (measurement process) there is an explicit flow of information from the present state of the X process; x_n+1 determines the present state of the Y process y_n+1 (assuming c_{− 1} is not zero). To fully capture the information transfer from the X process to the current state of the Y process we must identify the correct causal states [4]. For the measurement system, the causal states include the current (present) state. This state is not included in the definition of transfer entropy, being a mutual information quantity conditioned on only past states. Hence, for the purpose of this paper, we will temporarily define a quantity, “information transfer,” similar to transfer entropy, except that the present of the driving process, x_n+1, will be lumped in with the past values of the X process: x_n−k+2:x_n. For the first general model there is no information transferred from the Y to the X process. We define the (non-zero) information transfer from the X to the Y process (based on data from k timetags) as:

$$I{T}_{x\to y}^{\left(k\right)}=I\left(y_{n+1};x_{n-k+2:n+1}|y_{n-k+2:n}\right)=H\left(y_{n+1}|y_{n-k+2:n}\right)-H\left(y_{n+1}|y_{n-k+2:n},x_{n-k+2:n+1}\right).$$

(10)

The major contribution of this paper is to show how to analytically compute transfer entropy for AR Gaussian processes using an iterative method for computing the required covariance matrices. Computation of information transfer is additionally presented to elucidate the power of the method when similar information quantities are of interest and to make the measurement example more interesting. We now present a general method for computing the covariance matrices required to compute information-theoretic quantities for the AR models above. Two numerical examples follow.

To compute transfer entropy over a variable number of multiple time lags for AR processes of the general types shown above, we compute its block entropy components over multiple time lags. By virtue of the fact that the processes are Gaussian we can avail ourselves of analytical entropy expressions that depend only on the covariance of the processes. In this section we show how to analytically obtain the required covariance expressions starting with the covariance for a single time instance. Taking expectations, using the AR equations, we obtain the necessary statistics to characterize the process. Representing these expectation results in general, the process covariance matrix C⁽¹⁾(t_n) corresponding to a single timestamp, t_n, is:

$$C^{\left(1\right)}\left(t_n\right)\equiv \mathrm{cov}\left(\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]\right)=\left[\begin{array}{cc}E\left[x_n^2\right]& E\left[x_n{y}_n\right]\\ E\left[y_n{x}_n\right]& E\left[y_n^2\right]\end{array}\right].$$

(11)

To obtain an expanded covariance matrix, accounting for two time instances (t_n and t_n+1), we compute the additional expectations required to fill in the matrix C ⁽²⁾(t_n):

$$C^{\left(2\right)}\left(t_n\right)\equiv \mathrm{cov}\left(\left[\begin{array}{c}{x}_n\\ y_n\\ x_{n+1}\\ y_{n+1}\end{array}\right]\right)=\left[\begin{array}{cccc}E\left[x_n^2\right]& E\left[x_n{y}_n\right]& E\left[x_n{x}_{n+1}\right]& E\left[x_n{y}_{n+1}\right]\\ E\left[x_n{y}_n\right]& E\left[y_n^2\right]& E\left[x_{n+1}{y}_n\right]& E\left[y_n{y}_{n+1}\right]\\ E\left[x_n{x}_{n+1}\right]& E\left[x_{n+1}{y}_n\right]& E\left[x_{n+1}^2\right]& E\left[x_{n+1}{y}_{n+1}\right]\\ E\left[x_n{y}_{n+1}\right]& E\left[y_n{y}_{n+1}\right]& E\left[x_{n+1}{y}_{n+1}\right]& E\left[y_{n+1}^2\right]\end{array}\right].$$

(12)

Because the process is stationary, we may write:

$$C^{\left(2\right)}\left(t_n\right)=C^{\left(2\right)}=\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]$$

(13)

where:

$$\begin{array}{l}{\Sigma }_{11}\equiv \left[\begin{array}{cc}E\left[x_n^2\right]& E\left[x_n{y}_n\right]\\ E\left[x_n{y}_n\right]& E\left[y_n^2\right]\end{array}\right]\\ {\Sigma }_{12}\equiv \left[\begin{array}{cc}E\left[x_n{x}_{n+1}\right]& E\left[x_n{y}_{n+1}\right]\\ E\left[x_{n+1}{y}_n\right]& E\left[y_n{y}_{n+1}\right]\end{array}\right]\\ {\Sigma }_{21}={\Sigma }_{12}^T\\ {\Sigma }_{22}={\Sigma }_{11}.\end{array}$$

(14)

Thus we have found the covariance matrix C⁽²⁾ required to compute block entropies based on two timetags or, equivalently, one time lag. Using this matrix the single-lag transfer entropies may be computed.

We now show how to compute the covariance matrices corresponding to any finite number of time stamps. Define vector ${\overline{z}}_n=\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]$. Using the definitions above, write the matrix C⁽²⁾ as a block matrix and, using standard formulas, compute the conditional mean and covariance C_c of ${\overline{z}}_{n+1}$ given ${\overline{z}}_n$:

$$\begin{array}{l}{C}^{\left(2\right)}=\mathrm{cov}\left(\left[\begin{array}{c}{\overline{z}}_n\\ {\overline{z}}_{n+1}\end{array}\right]\right)=E\left[\left(\left[\begin{array}{c}{\overline{z}}_n\\ {\overline{z}}_{n+1}\end{array}\right]\left[\begin{array}{cc}{\overline{z}}_n& {\overline{z}}_{n+1}\end{array}\right]\right)\right]=\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]\\ E\left[{\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right]=E\left[{\overline{z}}_n\right]+{\Sigma }_{21}{\Sigma }_{11}^{-1}\left[\overline{z}-E\left[{\overline{z}}_n\right]\right]\\ ={\overline{\mu }}_{\overline{z}}+{\Sigma }_{21}{\Sigma }_{11}^{-1}\left[\overline{z}-{\overline{\mu }}_{\overline{z}}\right]\\ C_c\equiv \mathrm{cov}\left[{\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right]={\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{12}.\end{array}$$

(15)

Note that the expected value of the conditional mean is zero since the mean of the ${\overline{z}}_n$ process, ${\overline{\mu }}_{\overline{z}}$, is itself zero.

With these expressions in hand, we note that we may view propagation of the state ${\overline{z}}_n$ to its value ${\overline{z}}_{n+1}$ at the next timestamp as accomplished by the recursion:

$$\begin{array}{l}{\overline{z}}_{n+1}={\overline{\mu }}_z+D\left({\overline{z}}_n-{\overline{\mu }}_z\right)+S{\overline{u}}_n:{\overline{u}}_n~N\left(0_2,I_2\right)\\ D\equiv {\Sigma }_{21}{\Sigma }_{11}^{-1}\\ C_c\equiv S{S}^T\equiv {\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{12}.\end{array}$$

(16)

Here S is the principal square root of the matrix C_c. It is conveniently computed using the inbuilt Matlab function sqrtm. To see analytically that the recursion works, note that using it we recover at each timestamp a process having the correct mean and covariance:

$$E\left\{{\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right\}=E\left\{{\overline{\mu }}_z+D\left({\overline{z}}_n-{\overline{\mu }}_z\right)+S{\overline{u}}_n|{\overline{z}}_n=\overline{z}\right\}={\overline{\mu }}_z+D\left(\overline{z}-{\overline{\mu }}_z\right)$$

(17)

and:

$$\begin{array}{l}{\overline{z}}_{n+1}-E\left\{{\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right\}={\overline{\mu }}_z+D\left({\overline{z}}_n-{\overline{\mu }}_z\right)+S{\overline{u}}_n-\left({\overline{\mu }}_z+D\left(\overline{z}-{\overline{\mu }}_z\right)\right)=S{\overline{u}}_n+D\left({\overline{z}}_n-\overline{z}\right)\\ \\ \mathrm{cov}\left({\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right)=E\left\{\left[{\overline{z}}_{n+1}-E\left\{{\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right\}\right]{\left[{\overline{z}}_{n+1}-E\left\{{\overline{z}}_{n+1}|{\overline{z}}_n=\overline{z}\right\}\right]}^T|{\overline{z}}_n=\overline{z}\right\}\\ =E\left\{\left[S{\overline{u}}_n+D\left({\overline{z}}_n-\overline{z}\right)\right]{\left[S{\overline{u}}_n+D\left({\overline{z}}_n-\overline{z}\right)\right]}^T|{\overline{z}}_n=\overline{z}\right\}\\ =E\left\{\left[S{\overline{u}}_n\right]{\left[S{\overline{u}}_n\right]}^T\right\}=SE\left\{{\overline{u}}_n{\overline{u}}_n^T\right\}{S}^T=S{S}^T.\end{array}$$

(18)

Thus, because the process is Gaussian and fully specified by its mean and covariance, we have verified that the recursive representation yields consistent statistics for the stationary AR system. Using the above insights, we may now recursively compute the covariance matrix C^(k) for a variable number (k) of timestamps. Note that C^(k) has dimensions of 2k × 2k. We denote 2 × 2 blocks of C^(k) as C^(k)_ij for i, j = 1,2, ..., k , where C^(k)_ij is the 2-by-2 block of C^(k) consisting of the four elements of C^(k) that are individually located in row 2i − 1 or 2i and column 2j − 1 or 2j.

The above recursion is now used to compute the block elements of C⁽³⁾. Then each of these block elements is, in turn, expressed in terms of block elements of C⁽²⁾. These calculations are shown in detail below where we have also used the fact that the mean of the z_n vector is zero:

$$\begin{array}{l}{C}_{ij}^{\left(3\right)}=C_{ij}^{\left(2\right)}:i=1,2;j=1,2\\ {\overline{z}}_{n+2}=D{\overline{z}}_{n+1}+S{\overline{u}}_{n+1}\\ =D\left[D{\overline{z}}_n+S{\overline{u}}_n\right]+S{\overline{u}}_{n+1}=D^2{\overline{z}}_n+DS{\overline{u}}_n+S{\overline{u}}_{n+1}\end{array}$$

(19)

$$\begin{array}{l}{C}_{13}^{\left(3\right)}=E\left[{\overline{z}}_n{\overline{z}}_{n+2}^T\right]=E\left[{\overline{z}}_n{\left(D^2{\overline{z}}_n+DS{\overline{u}}_n+S{\overline{u}}_{n+1}\right)}^T\right]={\Sigma }_{11}{\left[D^2\right]}^T\\ C_{31}^{\left(3\right)}={\left[C_{13}^{\left(3\right)}\right]}^T\end{array}$$

(20)

$$\begin{array}{l}{C}_{23}^{\left(3\right)}=E\left[{\overline{z}}_{n+1}{\overline{z}}_{n+2}^T\right]=E\left[\left(D{\overline{z}}_n+S{\overline{u}}_n\right){\left(D^2{\overline{z}}_n+DS{\overline{u}}_n+S{\overline{u}}_{n+1}\right)}^T\right]\\ =D{\Sigma }_{11}{\left[D^2\right]}^T+C_c{D}^T=D{C}_{13}^{\left(3\right)}+C_c{D}^T\\ C_{32}^{\left(3\right)}={\left[C_{23}^{\left(3\right)}\right]}^T\end{array}$$

(21)

$$\begin{array}{l}{C}_{33}^{\left(3\right)}=E\left[{\overline{z}}_{n+2}{\overline{z}}_{n+2}^T\right]=E\left[\left(D^2{\overline{z}}_n+DS{\overline{u}}_n+S{\overline{u}}_{n+1}\right){\left(D^2{\overline{z}}_n+DS{\overline{u}}_n+S{\overline{u}}_{n+1}\right)}^T\right]\\ =D^2{\Sigma }_{11}{\left[D^2\right]}^T+D{C}_c{D}^T+C_c=D{C}_{23}^{\left(3\right)}+C_c.\end{array}$$

(22)

By continuation of this calculation to larger timestamp blocks (k > 3), we find the following pattern that can be used to extend (augment) C^(k−1) to yield C^(k). The pattern consists of setting most of the augmented matrix equal to that of the previous one, and then computing two additional rows and columns for C^(k), k > 2, to fill out the remaining elements. The general expressions are:

$$\begin{array}{l}{C}_{m,n}^{\left(k\right)}=C_{m,n}^{\left(k-1\right)}:m,n=1,2,\mathrm{...},k-1\\ C_{1k}^{\left(k\right)}={\Sigma }_{11}{\left[D^{k-1}\right]}^T\\ C_{ik}^{\left(k\right)}=D{C}_{i-1,k}^{\left(k\right)}+C_c{\left[D^{k-i}\right]}^T:i=2,3,\dots k\\ C_{ki}^{\left(k\right)}={\left[C_{ik}^{\left(k\right)}\right]}^T:i=1,2,\mathrm{...},k.\end{array}$$

(23)

At this point in the development we have shown how to compute the covariance matrix:

$$C^{\left(k\right)}=\mathrm{cov}\left({\overline{z}}^{\left(k\right)}\right)=\mathrm{cov}\left({\left[\begin{array}{ccccccc}{x}_n& y_n& x_{n+1}& y_{n+1}& \cdots & x_{n+k-1}& y_{n+k-1}\end{array}\right]}^T\right)$$

(24)

Since the system is linear and the process noise w_n and measurement noise v_n are white zero-mean Gaussian noise processes, we may express the joint probability density function for the 2k variates as:

$$f\left({\overline{z}}^{\left(k\right)}\right)=pdf\left({\overline{z}}^{\left(k\right)}\right)=pdf\left(\left[\begin{array}{ccccccc}{x}_n& y_n& x_{n+1}& y_{n+1}& \cdots & x_{n+k-1}& y_{n+k-1}\end{array}\right]\right)=\frac{\exp \left\{-\frac{1}{2}{\left[{\overline{z}}^{\left(k\right)}\right]}^T{\left[C^{\left(k\right)}\right]}^{-1}\left[{\overline{z}}^{\left(k\right)}\right]\right\}}{{\left(2\pi \right)}^{\frac{n}{2}}{\left(\det \left[C^{\left(k\right)}\right]\right)}^{\frac{1}{2}}}$$

(25)

Note that the mean of all 2k variates is zero.

Finally, to obtain empirical confirmation of the equivalence of the covariance terms obtained using the original AR system and its recursive representation, numerical simulations were conducted. Using the example 1 system (below) 500 sequences were generated each of length one million. For each sequence the C⁽³⁾ covariance was computed. The error for all C⁽³⁾ matrices was then averaged, assuming that the C⁽³⁾ matrix calculated using the method based on the recursive representation was the true value. The result was that for each of the matrix elements, the error was less than 0.0071% of its true value. We are now in position to compute transfer entropies for a couple of illustrative examples.

4. Example 1: A One-Way Coupled System

For this example we consider the following system:

$$\begin{array}{l}{x}_{n+1}=a{x}_n+w_n:w_n~N\left(0,Q\right)\\ y_{n+1}=h_c{x}_{n+1}+v_n:v_n~N\left(0,R\right)\end{array}$$

(26)

Parameter hc specifies the coupling strength of the Y process to the first-order AR process X, and R and Q are their respective (wn and vn) zero-mean Gaussian process noise variances. For stability, we require |a| <1. Comparing to the first general representation given above, we have m = 0, $a_0=a,$ and $c_{-1}=h_{c}a$. The system models filtered noise x_n and a noisy measurement, y_n, of x_n. Thus the x_n sequence represents a hidden process (or model) which is observable by way of another sequence, y_n. We wish to examine the behavior of transfer entropy as a function of the correlation ρ between x_n and y_n. One might expect that the correlation ρ between x_n and y_n to be proportional of the degree of information flow; however, we will see that the relationship between transfer entropy and correlation is not quite that simple.

Both the X and Y processes have zero mean. Computing the joint covariance matrix C⁽¹⁾ for x_n and y_n and their correlation we obtain:

$$\begin{array}{l}Var\left(x_n\right)=\frac{Q}{1-a^2}\\ Var\left(y_n\right)=h_c^{2}Var\left(x_n\right)+R\\ E\left(x_n{y}_n\right)=h_{c}Var\left(x_n\right)\\ \rho \equiv \frac{E\left(x_n{y}_n\right)}{\sqrt{Var\left(x_n\right)Var\left(y_n\right)}}\end{array}$$

(27)

Hence the process covariance matrix C⁽¹⁾ corresponding to a single timestamp, t_n is:

$$C^{\left(1\right)}\equiv \mathrm{cov}\left(\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]\right)=\left[\begin{array}{cc}Var\left(x_n\right)& hVar\left(x_n\right)\\ h_{c}Var\left(x_n\right)& h_c^{2}Var\left(y_n\right)+R\end{array}\right].$$

(28)

In order to obtain an expanded covariance matrix, accounting for two time instances (t_n and t_n+1) we compute the additional expectations required to fill in the matrix C⁽²⁾:

$$C^{\left(2\right)}\equiv \mathrm{cov}\left(\left[\begin{array}{c}{x}_n\\ y_n\\ x_{n+1}\\ y_{n+1}\end{array}\right]\right)=\left[\begin{array}{cccc}Var\left(x_n\right)& h_{c}Var\left(x_n\right)& aVar\left(x_n\right)& h_{c}aVar\left(x_n\right)\\ h_{c}Var\left(x_n\right)& h_c^{2}Var\left(x_n\right)+R& h_{c}aVar\left(x_n\right)& h_c^{2}aVar\left(x_n\right)\\ aVar\left(x_n\right)& h_{c}aVar\left(x_n\right)& Var\left(x_n\right)& h_{c}Var\left(x_n\right)\\ h_{c}aVar\left(x_n\right)& h_c^{2}aVar\left(x_n\right)& h_{c}Var\left(x_n\right)& h_c^{2}Var\left(x_n\right)+R\end{array}\right].$$

(29)

Thus we have found the covariance matrix C⁽²⁾ required to compute block entropies based on a single time lag. Using this matrix the single-lag transfer entropies may be computed. Using the recursive process described in the previous section we can compute C⁽¹°⁾. We have found that using higher lags does not change the entropy values significantly.

To aid the reader in understanding the calculations required to compute transfer entropies using higher time lags, it is worthwhile to compute transfer entropy for a single lag. We first define transfer entropy using general notation indicating the partitioning of the X and Y sequences in to past and future $\left(\overleftarrow{x},\overrightarrow{x}\right)$ and $\left(\overleftarrow{y},\overrightarrow{y}\right)$, respectively. We then compute transfer entropy as a sum of block entropies:

$$\begin{array}{l}T{E}_{x->y}=I\left(\overleftarrow{x};\overrightarrow{y}|\overleftarrow{y}\right)=h\left(\overleftarrow{x}|\overleftarrow{y}\right)+h\left(\overrightarrow{y}|\overleftarrow{y}\right)-h\left(\overleftarrow{x};\overrightarrow{y}|\overleftarrow{y}\right)\\ =\left[h\left(\overleftarrow{x},\overleftarrow{y}\right)-h\left(\overleftarrow{y}\right)\right]+\left[h\left(\overleftarrow{y},\overrightarrow{y}\right)-h\left(\overleftarrow{y}\right)\right]-\left[h\left(\overleftarrow{x},\overleftarrow{y},\overrightarrow{y}\right)-h\left(\overleftarrow{y}\right)\right]\\ =h\left(\overleftarrow{x},\overleftarrow{y}\right)+h\left(\overleftarrow{y},\overrightarrow{y}\right)-h\left(\overleftarrow{y}\right)-h\left(\overleftarrow{x},\overleftarrow{y},\overrightarrow{y}\right).\end{array}$$

(30)

Similarly:

$$T{E}_{y->x}=I\left(\overleftarrow{y};\overrightarrow{x}|\overleftarrow{x}\right)=h\left(\overleftarrow{x},\overleftarrow{y}\right)+h\left(\overleftarrow{x},\overrightarrow{x}\right)-h\left(\overleftarrow{x}\right)-h\left(\overleftarrow{x},\overleftarrow{y},\overrightarrow{x}\right)$$

(31)

The Y states have no influence on the X sequence in this example. Hence TE_y→x = 0. Since we are here computing transfer entropy for a single lag (i.e., two time tags t_n and t_n+1) we have:

$$T{E}_{x->y}^{\left(2\right)}=I\left(x_n;y_{n+1}|y_n\right)=h\left(x_n,y_n\right)+h\left(y_n,y_{n+1}\right)-h\left(y_n\right)-h\left(x_n,y_n,y_{n+1}\right)$$

(32)

By substitution of the expression for the differential entropy of each block we obtain:

$$\begin{array}{l}T{E}_{x->y}^{\left(2\right)}=\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^2\det C_{[1,2],[1,2]}^{\left(2\right)}\right]+\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^2\det C_{[2,4],[2,4]}^{\left(2\right)}\right]-\\ \frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^1\det C_{[2],[2]}^{\left(2\right)}\right]-\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^3\det C_{[1,2,4],[1:,2,4]}^{\left(2\right)}\right]\\ =\frac{\text{1}}{\text{2}}\log \left[\frac{\det C_{[1,2],[1,2]}^{\left(2\right)}\det C_{[2,4],[2,4]}^{\left(2\right)}}{\det C_{[2],[2]}^{\left(2\right)}\det C_{[1,2,4],[1,2,4]}^{\left(2\right)}}\right].\end{array}$$

(33)

For this example, note from the equation for y_n+1 that state x_n+1 is a causal state of X influencing the value of y_n+1. In fact, it is the most important such state. To capture the full information that is transferred from the X process to the Y process over the course of two time tags we need to include state x_n+1. Hence we compute the information transfer from x → y as:

$$I{T}_{x->y}^{\left(2\right)}=I\left(x_n,x_{n+1};y_{n+1}|y_n\right)=h\left(x_n,x_{n+1},y_n\right)+h\left(y_n,y_{n+1}\right)-h\left(y_n\right)-h\left(x_n,x_{n+1},y_n,y_{n+1}\right)$$

(34)

$$\begin{array}{l}I{T}_{x->y}^{\left(2\right)}=\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^3\det C_{[1,2,3],[1,2,3]}^{\left(2\right)}\right]+\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^2\det C_{[2,4],[2,4]}^{\left(2\right)}\right]-\\ \frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^1\det C_{[2],[2]}^{\left(2\right)}\right]-\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^4\det C_{[1:4],[1:4]}^{\left(2\right)}\right]\\ =\frac{\text{1}}{\text{2}}\log \left[\frac{\det C_{[1:,2,3],[1,2,3]}^{\left(2\right)}\det C_{[2,4],[2,4]}^{\left(2\right)}}{\det C_{[2],[2]}^{\left(2\right)}\det C_{[1:4],[1:4]}^{\left(2\right)}}\right].\end{array}$$

(35)

Here the notation $\det C_{[i],[i]}^{\left(2\right)}$ indicates the determinant of the matrix composed of the rows and columns of C⁽²⁾ indicated by the list of indices i shown in the subscripted brackets. For example, $\det C_{[1:4],[1:4]}^{\left(2\right)}$ is the determinant of the matrix formed by extracting columns {1, 2, 3, 4} and rows {1, 2, 3, 4} from matrix C⁽²⁾. In later calculations we will use slightly more complicated-looking notation. For example, $\det C_{[2:2:20],[2:2:20]}^{\left(10\right)}$ is the determinant of the matrix formed by extracting columns {2, 4 ,…, 18, 20} and the same-numbered rows from matrix C⁽¹°⁾. (Note C^(k)_[i],[i] is not the same as C^(k)_ii as used in Section 3).

It is interesting to note that a simplification in the expression for information transfer can be obtained by writing the expression for it in terms of conditional entropies:

$$I{T}_{x->y}^{\left(2\right)}=I\left(x_n,x_{n+1};y_{n+1}|y_n\right)=h\left(y_{n+1}|y_n\right)-h\left(y_{n+1}|x_n,y_n,x_{n+1}\right)$$

(36)

From the fact that y_n+1 = x_n+1 + v_n+1 we see immediately that:

$$h\left(y_{n+1}|x_n,y_n,x_{n+1}\right)=h\left(v_{n+1}\right)=\frac{1}{2}\log \left(2\pi eR\right).$$

(37)

Hence we may write:

$$\begin{array}{l}I{T}_{x->y}^{\left(2\right)}=h\left(y_{n+1}|y_n\right)-h\left(y_{n+1}|x_n,y_n,x_{n+1}\right)\\ =\frac{1}{2}\log \left[\frac{2\pi e\det C_{[2,4],[2,4]}^{\left(2\right)}}{\det C_{[2],[2]}^{\left(2\right)}}\right]-\frac{1}{2}\log \left[2\pi eR\right]\\ =\frac{1}{2}\log \left[\frac{\det C_{[2,4],[2,4]}^{\left(2\right)}}{R\det C_{[2],[2]}^{\left(2\right)}}\right].\end{array}$$

(38)

To compute transfer entropy using nine lags (ten timestamps) assume that we have already computed C⁽¹⁰⁾ as defined above. We partition the sequence ${\left\{{\overline{z}}_{n+i}^T\right\}}_{i=0}^9=\left\{x_n,y_n,x_{n+1},y_{n+1},x_{n+2},y_{n+2},x_{n+3},y_{n+3},x_{n+4},y_{n+4},x_{n+5},y_{n+5},x_{n+6},y_{n+6},x_{n+7},y_{n+7},x_{n+8},y_{n+8},x_{n+9},y_{n+9}\right\}$ into three subsets:

$$\begin{array}{l}\overleftarrow{x}\equiv \left\{x_n,x_{n+1},\dots ,x_{n+8}\right\}\\ \overleftarrow{y}\equiv \left\{y_n,y_{n+1},,\dots ,y_{n+8}\right\}\\ \overrightarrow{y}\equiv \left\{y_{n+9}\right\}.\end{array}$$

(39)

Now, using these definitions, and substituting in expressions for differential block entropies we obtain:

$$\begin{array}{l}T{E}_{x->y}^{\left(10\right)}=I\left(\overleftarrow{x};\overrightarrow{y}|\overleftarrow{y}\right)=h\left(\overleftarrow{x},\overleftarrow{y}\right)+h\left(\overleftarrow{y},\overrightarrow{y}\right)-h\left(\overleftarrow{y}\right)-h\left(\overleftarrow{x},\overleftarrow{y},\overrightarrow{y}\right)\\ =\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^{18}\det C_{[1:18],[1:18]}^{\left(10\right)}\right]+\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^{10}\det C_{[2:2:20],[2:2:20]}^{\left(10\right)}\right]-\\ \frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^9\det C_{[2:2:18],[2:2:18]}^{\left(10\right)}\right]-\frac{\text{1}}{\text{2}}\log \left[{\left(2\pi e\right)}^{19}\det C_{[1:18,20],[1:18,20]}^{\left(10\right)}\right]\\ =\frac{\text{1}}{\text{2}}\log \left[\frac{\det C_{[1:18],[1:18]}^{\left(10\right)}\det C_{[2:2:20],[2:2:20]}^{\left(10\right)}}{\det C_{[2:2:18],[2:2:18]}^{\left(10\right)}\det C_{[1:18,20],[1:18,20]}^{\left(10\right)}}\right].\end{array}$$

(40)

Similarly:

$$I{T}_{x->y}^{\left(10\right)}=h\left(\overleftarrow{y}|\overleftarrow{y}\right)-h\left(\overrightarrow{y}|\overleftarrow{y},\overleftarrow{x},x_{n+1}\right)=\frac{\text{1}}{\text{2}}\log \left[\frac{\det C_{[2:2:20],[2:2:20]}^{\left(10\right)}}{R\det C_{[2:2:18],[2:2:18]}^{\left(10\right)}}\right].$$

(41)

As a numerical example we set h_c = 1, Q = 1, and for three different values of a (0.5, 0.7 and 0.9) we vary R so as to scan the correlation ρ between the x and y processes between the values of 0 and 1.

In Figure 1 it is seen that for each value of parameter a there is a peak in the transfer entropy TE^(k)_x→y. As the correlation ρ between x_n and y_n increases from a low value the transfer entropy increases since the amount of information shared between y_n+1 and x_n is increasing. At a critical value of ρ transfer entropy peaks and then starts to decrease. This decrease is due to the fact that at high values of ρ the measurement noise variance R is small. Hence y_n becomes very close to equaling x_n so that the amount of information gained (about y_n+1) by learning x_n, given y_n, becomes small. Hence h(y_n+1 | y_n) ‑ h(y_n+1 | y_n, x_n) is small. This difference is TE⁽²⁾_x→y.

Figure 1. Example 1: Transfer entropy TE^(k)_x→y versus correlation coefficient ρ for three values of parameter a (see legend). Solid trace: k = 10, dotted trace: k = 2.

Figure 1. Example 1: Transfer entropy TE^(k)_x→y versus correlation coefficient ρ for three values of parameter a (see legend). Solid trace: k = 10, dotted trace: k = 2.

The relationship between ρ and R is shown in Figure 2. Note that when parameter a is increased, a larger value of R is required to maintain ρ at a fixed value. Also, in Figure 1 we see the effect of including more timetags in the analysis. When k is increased from 2 to 10 transfer entropy values fall, particularly for the largest value of parameter a. It is known that entropies decline when conditioned on additional variables. Here, transfer entropy is acting similarly. In general, however, transfer entropy, being a mutual information quantity, has the property that conditioning could make it increase as well [12].

Figure 2. Example 1: Logarithm of R versus ρ for three values of parameter a (see legend).

Figure 2. Example 1: Logarithm of R versus ρ for three values of parameter a (see legend).

The observation that the transfer entropy decrease is greatest for the largest value of parameter a is perhaps due to the fact that the entropy of the X process is itself greatest for the largest a value and therefore has more sensitivity to an increase in X data availability (Figure 3).

From Figure 1 it is seen that as the value of parameter a is increased, transfer entropy is increased for a fixed value of ρ. The reason for this increase may be gleaned from Figure 3 where it is clear that the amount of information contained in the x process, H_X, is greater for larger values of a. Hence more information is available to be transferred at the fixed value of ρ when a is larger. In the lower half of Figure 3 we see that as ρ increases the entropy of the Y process, H_Y, approaches the value of H_X. This result is due to the fact that the mechanism being used to increase ρ is to decrease R. Hence as R drops close to zero y_n looks increasingly identical to x_n (since h_c = 1).

Figure 3. Example 1: Process entropies H_X and H_Y versus correlation coefficient ρ for three values of parameter a (see legend).

Figure 3. Example 1: Process entropies H_X and H_Y versus correlation coefficient ρ for three values of parameter a (see legend).

Figure 4 shows information transfer IT^(k)_x→y plotted versus correlation coefficient ρ. Now note that the trend is for information transfer to increase as ρ is increased over its full range of values. °

Figure 4. Example 1: Information transfer IT^(k)_x→y versus correlation coefficient ρ for three different values of parameter a (see legend) for k = 10 (solid trace) and k = 2 (dotted trace).

Figure 4. Example 1: Information transfer IT^(k)_x→y versus correlation coefficient ρ for three different values of parameter a (see legend) for k = 10 (solid trace) and k = 2 (dotted trace).

This result is obtained since as ρ is increased y_n+1 becomes increasingly correlated with x_n+1. Also, for a fixed ρ, the lowest information transfer occurs for the largest value of parameter a. We obtain this result since at the higher a values x_n and x_n+1 are more correlated. Thus the benefit of learning the value of y_n+1 through knowledge of x_n+1 is relatively reduced, given that y_n (itself correlated with x_n) is presumed known. Finally, we have IT⁽¹⁰⁾_x→y < IT⁽²⁾_x→y since conditioning the entropy quantities comprising the expression for information transfer with more state data acts to reduce their difference. Also, by comparison of Figure 2 and Figure 4, it is seen that information transfer is much greater than transfer entropy. This relationship is expected since information transfer as defined herein (for k = 2) is the amount of information that is gained about y_n+1 from learning x_n+1 and x_n, given that y_n is already known. Whereas transfer entropy (for k = 2) is the information gained about y_n+1 from learning only x_n, given that y_n is known. Since the state y_n+1 in fact equals x_n+1, plus noise, learning x_n+1 is highly informative, especially when the noise variance is small (corresponding to high values of ρ). The difference between transfer entropy and information transfer therefore quantifies the benefit of learning x_n+1, given that x_n and y_n are known (when the goal is to determine y_n+1).

Figure 5 shows how information transfer varies with measurement noise variance R. As R increases the information transfer decreases since measurement noise makes determination of the value of y_n+1 from knowledge of x_n and x_n+1 less accurate. Now, for a fixed R, the greatest value for information transfer occurs for the greatest value of parameter a. This is the opposite of what we obtained for a fixed value of ρ as shown in Figure 4. The way to see the rationale for this is to note that, for a fixed value of information transfer, R is highest for the largest value of parameter a. This result is obtained since larger values of a yield the most correlation between states x_n and x_n+1. Hence, even though the measurement y_n+1 of x_n+1 is more corrupted by noise (due to higher R), the same information transfer is achieved nevertheless, because x_n provides a good estimate of x_n+1 and, thus, of y_n+1.

Figure 5. Example 1: Information transfer IT⁽¹⁰⁾_x→y versus measurement error variance R for three different values of parameter a (see legend).

Figure 5. Example 1: Information transfer IT⁽¹⁰⁾_x→y versus measurement error variance R for three different values of parameter a (see legend).

5. Example 2: Information-theoretic Analysis of Two Coupled AR Processes.

In example 1 the information flow was unidirectional. We now consider a bidirectional example achieved by coupling two AR processes. One question we may ask in such a system is how transfer entropies change with variations in correlation and coupling coefficient parameters. It might be anticipated that increasing either of these quantities will have the effect of increasing information flow and thus transfer entropies will increase.

The system is defined by the equations:

$$\begin{array}{l}{x}_{n+1}=a{x}_n+b{y}_n+w_n:w_n~N\left(0,Q\right)\\ y_{n+1}=c{x}_n+d{y}_n+v_n:v_n~N\left(0,R\right).\end{array}$$

(42)

For stability, we require that the eigenvalues of the constant matrix $\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$ lie in the unit circle. The means of processes X and Y are zero. The terms w_n and v_n are the X and Y processes noise terms respectively. Using the following definitions:

$$\begin{array}{l}{\lambda }_0\equiv 1+ad-bc\\ {\lambda }_1\equiv 1-ad-bc\\ {\psi }_a\equiv \left(1-ad\right)\left(1-a^2\right)-bc\left(1+a^2\right)\\ {\psi }_d\equiv \left(1-ad\right)\left(1-d^2\right)-bc\left(1+d^2\right)\\ \tau \equiv {\psi }_a{\psi }_d-b^2{c}^2{\lambda }_0^2\\ {\eta }_{x1}\equiv {\lambda }_1{\psi }_d/\tau \\ {\eta }_{x2}\equiv b^2{\lambda }_0{\lambda }_1/\tau \\ {\eta }_{y1}\equiv c^2{\lambda }_0{\lambda }_1/\tau \\ {\eta }_{y2}\equiv {\lambda }_1{\psi }_a/\tau \end{array}$$

(43)

we may solve for the correlation coefficient ρ between x_n and y_n to obtain:

$$\left[\begin{array}{c}Var\left(x_n\right)\\ Var\left(y_n\right)\end{array}\right]=\left[\begin{array}{cc}{\eta }_{x1}& {\eta }_{x2}\\ {\eta }_{y1}& {\eta }_{y2}\end{array}\right]\left[\begin{array}{c}Q\\ R\end{array}\right].$$

(44)

$$\begin{array}{l}{C}_{\left[xy\right]}\equiv \mathrm{cov}\left(\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]\right)=\left[\begin{array}{cc}Var\left(x_n\right)& \xi \\ \xi & Var\left(y_n\right)\end{array}\right]\\ \xi \equiv E\left[x_n{y}_n\right]=\frac{b\left(d{\psi }_a+abc{\lambda }_0\right)R+c\left(a{\psi }_d+bcd{\lambda }_0\right)Q}{{\psi }_a{\psi }_d-b^2{c}^2{\lambda }_0^2}\\ \rho =\frac{\xi }{\sqrt{Var\left(x_n\right)Var\left(y_n\right)}}.\end{array}$$

(45)

Now, as we did previously in example 1 above, compute the covariance C⁽²⁾ of the variates obtained at two consecutive timestamps to yield:

$$C^{\left(2\right)}\equiv \mathrm{cov}\left(\left[\begin{array}{c}{x}_n\\ y_n\\ x_{n+1}\\ y_{n+1}\end{array}\right]\right)=\left[\begin{array}{cccc}Var\left(x_n\right)& \xi & aVar\left(x_n\right)+b\xi & cVar\left(x_n\right)+d\xi \\ \xi & Var\left(y_n\right)& bVar\left(y_n\right)+a\xi & dVar\left(y_n\right)+c\xi \\ aVar\left(x_n\right)+b\xi & bVar\left(y_n\right)+a\xi & Var\left(x_n\right)& \xi \\ cVar\left(x_n\right)+d\xi & dVar\left(y_n\right)+c\xi & \xi & Var\left(y_n\right)\end{array}\right].$$

(46)

At this point the difficult part is done and the same calculations can be made as in example 1 to obtain C^(k); k = 3,4, …, 10 and transfer entropies. For illustration purposes, we define the parameters of the system as shown below, yielding a symmetrically coupled pair of processes. To generate a family of curves for each transfer entropy we choose a fixed coupling term ε from a set of four values. We set Q = 1000 and vary R so that ρ varies from about 0 to 1. For each ρ value we compute the transfer entropies. The relevant system equations and parameters are:

$$\begin{array}{l}{x}_{n+1}=\left(\frac{1}{2}-\epsilon \right)x_n+\epsilon y_n+w_n:w_n~N\left(0,Q\right)\\ y_{n+1}=\epsilon x_n+\left(\frac{1}{2}-\epsilon \right)y_n+v_n:w_n~N\left(0,R\right)\\ \epsilon \in \left\{0.1,0.2,0.3,0.4\right\}\\ Q=1000.\end{array}$$

(47)

Hence, we make the following substitutions to compute C⁽²⁾:

$$\begin{array}{l}a=\left(\frac{1}{2}-\epsilon \right)\\ b=\epsilon \\ c=\epsilon \\ d=\left(\frac{1}{2}-\epsilon \right).\end{array}$$

(48)

For each parameter set {ε, Q, R} there is a maximum possible ρ, ${\rho }_{\infty }$ obtained by taking the limit as R→ ∞ of the expression for ρ given above. Doing so, we obtain:

$${\rho }_{\infty }=\frac{{\varphi }_1{\varphi }_2+{\varphi }_3}{\sqrt{{\varphi }_1\left({\varphi }_1{\mu }_1+1\right)}}$$

(49)

where:

$$\begin{array}{l}{\varphi }_1\equiv \frac{2a{b}^{2}d+b^2{\lambda }_1}{\left(1-a^2-b^2{\mu }_1\right){\lambda }_1-2ab\left(ac+bd{\mu }_1\right)}\\ {\varphi }_2\equiv \frac{ac+bd{\mu }_1}{{\lambda }_1}\\ {\varphi }_2\equiv \frac{bd}{{\lambda }_1}\end{array}$$

(50)

$$\begin{array}{l}{\lambda }_1\equiv 1-ad-bc\\ {\mu }_1\equiv \frac{c^2{\lambda }_1+2a{c}^{2}d}{\left(1-d^2\right){\lambda }_1-2bc{d}^2}.\end{array}$$

(51)

There is a minimum value of ρ also. The corresponding value for R, R_min, was found by means of the inbuilt Matlab program fminbnd. This program is designed to find the minimum of a function in this case ρ(a, b, c, d, R, Q)) with respect to one parameter (in this case R) starting from an initial guess (here, R = 500). The program returns the minimum functional value (ρ_min) and the value of the parameter at which the minimum is achieved (R_min). After identifying R_min a set of R values were computed so that the corresponding set of ρ values spanned from ρ_min to the maximum ${\rho }_{\infty }$ in fixed increments of Δρ (here equal to 0.002). This set of R values was generated using the iteration:

$$R_{new}=R_{old}+\Delta R=R_{old}+{{\left(\frac{\partial \rho }{\partial R}\right)}^{-1}|}_{R=R_{old}}\Delta \rho$$

(52)

For the four selections of parameter ε we obtain the functional relationships shown in Figure 6.

From Figure 6 we see that for a fixed ε, increasing R increases (or decreases) ρ depending on whether R is less than (or greater than) Q (Q = 1000). Note that large increases in R > Q are required to marginally increase ρ when ρ nears its maximum value. The reason that the minimum ρ value occurs when Q equals R is because whenever they are unequal one of the processes dominates the other, leading to increased correlation. Also, note that if R << Q, then increasing ε will cause ρ to decrease since increasing the coupling will cause the variance of the y process Var(y_n), a term appearing in the denominator of the expression for ρ, to increase. If Q << R, a similar result is obtained when ε is increased.

Figure 6. Example 2: Process noise variance R versus correlation coefficient ρ for a set of ε parameter values (see figure legend).

Figure 6. Example 2: Process noise variance R versus correlation coefficient ρ for a set of ε parameter values (see figure legend).

Transfer entropies in both directions are shown in Figure 7. Fixing ε, we note that as R is increased from a low value both ρ and TE_{x− >y} initially decrease while TE_{y− >x} increases. Then for further increases of R, ρ reaches a minimum value then begins to increase, while TE_x→y continues to decrease and TE_y→x continues to increase.

Figure 7. Example 2: Transfer entropy values versus correlation ρ for a set of ε parameter values (see figure legend). Arrows indicate direction of increasing R values.

Figure 7. Example 2: Transfer entropy values versus correlation ρ for a set of ε parameter values (see figure legend). Arrows indicate direction of increasing R values.

Figure 8. Example 2: Transfer entropies difference (TE_{x− >y} – TE_{y − > x}) and sum (TE_{x− > y} + TE_{y− > x}) versus correlation ρ for a set of ε parameter values (see figure legend). Arrow indicates direction of increasing R values.

Figure 8. Example 2: Transfer entropies difference (TE_{x− >y} – TE_{y − > x}) and sum (TE_{x− > y} + TE_{y− > x}) versus correlation ρ for a set of ε parameter values (see figure legend). Arrow indicates direction of increasing R values.

By plotting the difference TE_x→y – TE_y→x in Figure 8 we see the symmetry that arises as R increases from a low value to a high value. What is happening is that when R is low, the X process dominates the Y process so that TE_x→y > TE_y→x. As R increases, the two entropies equilibrate. Then, as R rises above Q, the Y process dominates giving TE_x→y < TE_y→x. The sum of the transfer entropies shown in Figure 8 reveal that the total information transfer is minimal at the minimum value of ρ and increases monotonically with ρ. The minimum value for ρ in this example occurs when the process noise variances Q and R are equal (matched). Figure 9 shows the changes in the transfer entropy values explicitly as a function of R. Clearly, when R is small (as compared to Q = 1000), TE_x→y > TE_y→x. Also it is clear that at every fixed value of R, both transfer entropies are higher at the larger values for the coupling term ε.

Figure 9. Example 2: Transfer entropies TE_x→y and TE_{y→x versus} process noise variance R for a set of ε parameter values (see figure legend).

Figure 9. Example 2: Transfer entropies TE_x→y and TE_{y→x versus} process noise variance R for a set of ε parameter values (see figure legend).

Another informative view is obtained by plotting one transfer entropy value versus the other as shown in Figure 10.

Figure 10. Example 2: Transfer entropy TE_{x− > y} plotted versus TE_{y− > x} for a set of ε parameter values (see figure legend). The black diagonal line indicates locations where equality obtains. Arrow indicates direction of increasing R values.

Figure 10. Example 2: Transfer entropy TE_{x− > y} plotted versus TE_{y− > x} for a set of ε parameter values (see figure legend). The black diagonal line indicates locations where equality obtains. Arrow indicates direction of increasing R values.

Here it is evident how TE_y→x increases from a value less than TE_x→y to a value greater than TE_x→y as R increases. Note that for higher coupling values ε this relative increase is more abrupt.

Finally, we consider the sensitivity of the transfer entropies to the coupling term ε. We reprise example system 2 where now ε is varied in the interval (0, ½) and three values of R (somewhat arbitrarily selected to provide visually appealing figures to follow) are considered:

$$\begin{array}{l}{x}_{n+1}=\left(\frac{1}{2}-{\epsilon }_x\right)x_n+{\epsilon }_x{y}_n+w_n:w_n~N\left(0,Q\right)\\ y_{n+1}={\epsilon }_y{x}_n+\left(\frac{1}{2}-{\epsilon }_y\right)y_n+v_n:w_n~N\left(0,R\right)\\ R\in \left\{{10}^0,{10}^3,{10}^4\right\}\\ Q={10}^3.\end{array}$$

(53)

Figure 11 shows the relationship between ρ and ε, where ε_x = ε_y = ε for the three R values. Note that for the case R = Q the relationship is symmetric around ε = ¼. As R departs from equality more correlation between x_n and y_n is obtained.

Figure 11. Example 2: Correlation coefficient ρ vs coupling coefficient ε for a set of R values (see figure legend).

Figure 11. Example 2: Correlation coefficient ρ vs coupling coefficient ε for a set of R values (see figure legend).

The reason for this increase is that when the noise driving one process is greater in amplitude than the amplitude of the noise driving the other process, the first process becomes dominant over the other. This domination increases as the disparity between the process noise variances increases (R versus Q). Note also that as the disparity increases, the maximum correlation occurs at increasingly lower values of the coupling term ε. As the disparity increases at fixed ε = ¼ the correlation coefficient ρ increases. However, the variance in the denominator of ρ can be made smaller and thus ρ larger, if the variance of either of the two processes can be reduced. This can be accomplished by reducing ε.

The sensitivities of the transfer entropies to changes in coupling term ε are shown in Figure 12. Consistent with intuition, all entropies increase with increasing ε. Also, when R < Q (blue trace) we have TE_x‑>y > TE_y‑>x and the reverse for R > Q. (red). For R = Q, TE_x‑>y = TE_y‑>x (green).

Figure 12. Example 2: Transfer entropies TE_x→y (solid lines) vs TE_y→x (dashed lines) vs coupling coefficient ε for a set of R values (see figure legend).

Figure 12. Example 2: Transfer entropies TE_x→y (solid lines) vs TE_y→x (dashed lines) vs coupling coefficient ε for a set of R values (see figure legend).

Finally, it is interesting to note that whenever we define three cases by fixing Q and varying the setting for R ( one of R₁, R₂ and R₃ for each case) such that R₁ < Q, R₂ = Q and R₃ = Q²/R₁ (so that R_i+1 = QR_i/R₁ for i = 1 and i = 2) we then obtain the symmetric relationships TE_{x ‑ >y}(R₁) = TE_{y ‑ >x}(R₃) and TE_{x ‑ >y}(R₃) = TE_{y ‑ >x}(R₁) for all ε in the interval (1, ½). For these cases we also obtain ρ(R₁) = ρ(R₃) on the same ε interval.

6. Conclusions

It has been shown how to compute transfer entropy values for Gaussian autoregressive processes for multiple timetags. The approach is based on the iterative computation of covariance matrices. Two examples were investigated: (1) a first-order filtered noise process whose state is measured with additive noise, and (2) two first-order symmetrically coupled processes each of which is driven by independent process noise. We found that, for the first example, increasing the first-order AR coefficient at a fixed correlation coefficient, transfer entropy increased since the entropy of the measured process was itself increased.

For the second example, it was discovered that the relationships between the coupling and correlation coefficients and the transfer entropies is more complicated. The minimum correlation coefficient occurs when the process noise variances match. It was seen that matching of these variances results in minimum information flow, expressed as the sum of both transfer entropies. Without a match, the transfer entropy is larger in the direction away from the process having the larger process noise. Fixing the process noise variances, transfer entropies in both directions increase with coupling strength ε.

Finally, it is worth noting that the method for computing covariance matrices for a variable number of timetags as presented here facilitates the calculation of many other information-theoretic quantities of interest. To this purpose, the authors have computed such quantities as crypticity [13] and normalized transfer entropy using the reported approach.