PHYSICAL REVIEW E 77, 036217 共2008兲

Comparison of nonlinear Granger causality extensions for low-dimensional systems
1

Katsuhiko Ishiguro,1,* Nobuyuki Otsu,1,2,† Max Lungarella,1,‡ and Yasuo Kuniyoshi1,§

Graduate School of Information Science and Technology, University of Tokyo, 113-8656 Tokyo, Japan
2
National Institute of Advanced Industrial Science and Technology, 305-8568 Tsukuba, Japan
共Received 5 April 2007; revised manuscript received 4 November 2007; published 25 March 2008兲

The identification of hidden interdependences among the parts of a complex system is a fundamental issue.
Typically, the objective is, given only a sequence of scalar measurements, to infer as much as possible about
the internal dynamics of the system and about the interactions between its subsystems. In general, suchinteractions are not only nonlinear but also asymmetric. Constraints on the estimation of hidden relationships
are further posed by noise and by the length of signals sampled from real world systems. The focus of thispaper is causal dependences between bivariate time series. We especially focus on the nonlinear extension ofGranger causality with polynomial terms of the conventional embedding vector. In this paper, we study theperformance of this measure in comparison with three alternative methods proposed recently in lowdimensional and low-order-nonlinearity systems. Those methods are tested with three different artificial chaotic
maps with several noise contamination setups. As a result, we find that the polynomial embedding techniquesuccessfully detects asymmetric 共causal兲 dependences between bivariate time series in many low-dimensional
cases.
DOI: 10.1103/PhysRevE.77.036217

PACS number共s兲: 05.45.Tp, 02.50.Sk

I. INTRODUCTION

introduced 关1–5兴. Although all these measures are asymmetric by construction and can be shown to work in specificinstances, their general application is not a straightforward
subject 共for a review, see 关6兴兲. Real world signals 共e.g., neurophysiological signals兲 are typically noisy, nonstationary,and of finite sample size; moreover, their dependences areoften characterized by nonlinearities.
Here, we focus on another type of extension of Granger
causality. We augment the time-delayed embedding vector
by incorporating multivariate polynomial terms, and use theresulting nonlinear embedding vector to extend linear
Granger causality 关7兴. Polynomial extension of Granger causality was reported in 关8兴, but intensive comparison withother nonlinear extensions of Granger causality was not
made. In this paper, we provide a quantitative study of several methods’ performance on some numerical data sets. As afirst step, we test low-dimensional and low-order-nonlinear
systems in this paper.
In what follows, we first describe three previously introduced methods. We then introduce the notion of the polynomial embedding vector and show how it can be used to extend Granger causality. We proceed by exposing the resultsof a set of numerical experiments with chaotic maps. Finally,
we discuss possible extensions and alternatives to the proposed method.

The detection and characterization of dependences amonginteracting components of a complex system can give novel
insights into its functioning and lead to a better understanding of its dynamics. In most natural as opposed to laboratory
settings, however, we are not granted direct access to thesystem’s components, but possess instead only a set of strategically selected and typically simultaneously recorded variables. The problem of inferring and quantifying relationships
between the parts of the system 共e.g., how they contribute to
the generation of information and at what rate they exchange
information兲 is thus mapped onto the one of studying therelationships between multivariate time series corresponding
to the recorded variables.
A popular way of evaluating the statistical dependency
between time series is cross correlation. This method is not
only computationally efficient but also has an appealing andnatural interpretation. Cross correlation, however, is asecond-order statistics and handles merely linear dependences. Mutual information in the information theory provides an attractive way of circumventing the restrictions ofcross correlation because it is sensitive to higher-order relationships. Yet mutual information is a symmetric measure,and fails to detect directional 共asymmetric兲 flow of information unless one of the time series is delayed. Measures ofinterdependence or “causality structure” aimed at overcoming the limitations of mutual information have recently been

*Corresponding author. Present address: NTT Communication
ScienceLaboratories,
619-0237
Kyoto,
Japan.
ishiguro@cslab.kecl.ntt.co.jp†
otsu.n@aist.go.jp‡
Present address: Artificial Intelligence Laboratory, University of
Zurich, Switzerland. lunga@ifi.uzh.ch§kuniyosh@isi.imi.i.u-tokyo.ac.jp1539-3755/2008/77共3兲/036217共9兲

II. EXISTING METHODS

Given a complex system, it is most likely the case that its
internal dynamics is a priori unknown. There is also noknowledge of whether the system is deterministic or stochastic. Here, we assume that all we can do is to measure two or
more variables observed separately from different parts ofthe complex system.
A. Embedding vector

For simplicity, let us consider the bivariate case, and let us
assume that X = 兵x1 , x2 , . . . , xN其 and Y = 兵y 1 , y 2 , . . . , y N其 are

036217-1©2008 The American Physical Society

PHYSICAL REVIEW E 77, 036217 共2008兲

ISHIGURO et al.

two simultaneously measured time series each consisting of
N scalar quantities. The time-delayed embedding vector reconstructing the state 共or phase兲 space X of time series X isexpressed as␶Txm,t = 共xt,xt−␶,xt−2␶, . . . ,xt−共m−1兲␶兲 ,共1兲

where m is the embedding dimension and ␶ is the time delay
共or lag兲 between successive elements of the state vector
关9,10兴. The choice of both parameters depends on the dynamics underlying the data. Here, first guesses for the embedding dimension and the time lag were obtained by thefalse nearest neighbor technique 关11兴 and as the first minimum of the mutual information function 关12兴, respectively.
The time-delayed embedding vector reconstructing the state
space Y of time series Y can be defined analogously. If nototherwise stated in this paper, ␶ = 1, and we use the following
expression:T
xt = xmt = 共xt,xt−1,xt−2, . . . ,xt−m+1兲 .

autoregression to estimate the causal influence of the twotime series on each other. The autoregression models of Eqs.
共3兲 and 共4兲 are replaced by the following expressions:

xt+1 = ␣Txt + ␲ + ⑀共x兲
t ,

共3兲

,
xt+1 = aTxt + bTyt + p + ⑀共兩x兩y兲
t

共4兲

共兩x兩y兲
represent the prediction errors, ␣, a, andwhere ⑀共x兲t and ⑀tb denote regression coefficient vectors, and ␲ and p are constants. The coefficient vectors and constants are determined
共兩x兩y兲
. Once theso as to minimize the variance of ⑀共x兲t and ⑀tcoefficient vectors have been calculated, the causal influenceof Y on X can be expressed as兲
var共⑀共兩x兩y兲
t
GY⇒X = 1 −共x兲 .
var共⑀t 兲

共5兲Intuitively, the larger GY⇒X the stronger is the estimatedcausal influence of Y on X. In other words, if the predictionof xt improves by incorporating the past values of y t, then y t
is said to “Granger cause” xt. In a similar way, it is possible
共兩y兩x兲
, and GX⇒Y .
to define ⑀共y兲t , ⑀tGranger causality is formulated for linear models, and its
application to nonlinear systems may not be appropriate in
the general case. Here, we describe two recently introduced
nonlinear extensions of Granger causality. The first one isnonlinear Granger causality 共NLGC兲 关1兴. NLGC employs anonlinear kernel autoregression scheme instead of a linear,
xt+1 = aT⌽共xt兲 + bT⌽共yt兲 + p + ␩共兩x兩y兲
t

共7兲

冢 冣
兩xt − x␳1兩␴2
兩xt − x␳2兩
−
␴2
]
−

⌽共xt兲 −共8兲

兩xt − x␳ P兩␴2with ␴ a parameter defined a priori, and the center vectors
chosen using fuzzy c-means clustering 共other clustering algorithms can also be applied兲. The nonlinear Granger causal
influence of Y on X or X on Y can be expressed as thedifference between the variances of the two error terms

B. Granger causality and its extensions

The idea of Granger causality applied to bivariate timeseries 关7兴 relies on the notion that the prediction of one timeseries through linear regression can be improved by incorporating information about the past values of the other timeseries. To assess the dependence of X on Y 共subsequently
also denoted as Y ⇒ X兲, consider the following set of autoregressive predictions:共6兲

where ⌽ is a P-dimensional vector whose elements are nonlinear radial based functions 共RBFs兲 centered at x␳1 , . . . , x␳ P,

共2兲

Similarly, we can define yt.
Note that many of the existing causality measures 共including the ones described in this paper兲 are based on the notionof the embedding vector.

xt+1 = ␣T⌽共xt兲 + ␲ + ␩共x兲t ,共兩x兩y兲
兲,cY⇒X = var共␩共x兲
t 兲 − var共␩t

共9兲

共兩y兩x兲
兲.
cX⇒Y = var共␩共y兲t 兲 − var共␩t

共10兲

If cY⇒X ⬎ 0 then the prediction of xt improves by including
y t, and it can be said that Y causally affects X. Analogously,
if cX⇒Y ⬎ 0 then X has a causal influence on Y.
Another extension of Granger causality, called the extended Granger causality index 共EGCI兲, was proposed in 关4兴.
The rationale of this technique is to divide the phase space
into a set of small neighborhoods and approximate the globally nonlinear dynamics by local linear regression models.
Let us express the joint dynamics of state space X and state
space Y as
zt =

冉冊

xt
.
yt

共11兲

The resulting state vector zt 僆 R2m is a point in the
2m-dimensional reconstructed phase space Z. Given L ⬍ N
points in Z, we define L clusters of local neighbors ⌰k= 兵zt : 兩zt − zt共k兲兩 ⬍ d其, where d is the size of the neighborhoodand t共k兲 共k = 1 , . . . , L兲 is the index of an arbitrarily selected
centroid vector zt共k兲. The parameter d is the size of the neighborhood and is chosen of the order of the resolution of themeasurements. If the d neighborhood contains no points, d isincreased.
Once the entire reconstructed phase space Z has been segmented, the linear Granger “causalities” GY⇒X共k兲 and
GX⇒Y 共k兲 are calculated for each local neighborhood ⌰k. By
averaging over all neighborhoods sampling Z 共i.e., L neighboring clusters兲, we obtain two extended Granger causality
indices

036217-2

PHYSICAL REVIEW E 77, 036217 共2008兲

COMPARISON OF NONLINEAR GRANGER CAUSALITY…
L

1
␦Y⇒X = 兺 GY⇒X共k兲,
L k=1

共12兲

joint probability p共xt+1 , xkt , ylt兲 which, in turn, is estimatedusing the histograms of the embedding vectors. The only
parameter is the number of bins, r.

L

1
␦X⇒Y = 兺 GX⇒Y 共k兲.
L k=1

共13兲

One relevant problem of both nonlinear extensions ofGranger causality is their dependency on at least two parameters which need to be estimated and tuned to the time series
at hand 共e.g., size and number of the local neighborhoods,
number of radial basis functions, and their parameter ␴兲.
Another drawback is their respective computational costs. In
the case of NLGC, the computational cost of the least-meansquare regression is polynomial in the number of RBFs. For
the EGCI, it is necessary to exhaustively search the state
space to determine the local neighborhood of every centroid
vector.

III. POLYNOMIAL EMBEDDING VECTOR ANDNONLINEAR GRANGER CAUSALITY

Here, we present an alternative way of extending linear
Granger causality, which not only performs better than
NLGC and the EGCI but solves the problem of parametertuning. We first extend the time-delayed embedding vector
by incorporating multivariate polynomial terms. We thengeneralize Granger’s idea by using the resulting “polynomial” embedding vector.
A. Polynomial embedding vector

Let us consider all possible combinations of order R ofdifferent monomial terms generated from the elements of xt,
k=j

C. Transfer entropy

Another approach to obtaining knowledge about asymmetric 共causal兲 dependences between coupled systems is bymeasuring to what extent its individual components contribute to information production and at what rate they exchange
information among each other. Transfer entropy is an information theoretic measure which not only shares some of theproperties of mutual information but also takes the dynamics
of information transport into account 关3兴. It is a specific version of the mutual information for conditional probabilities,
but unlike mutual information it is designed to detect thedirected exchange of information between two systems,
separate for both directions, and conditional to common history and input signals.
Transfer entropy measures the deviation from the following conditional independency assumption about the Markov
process:p共兩xt+1兩xkt 兲 = p共兩xt+1兩xkt ,ylt兲

TY⇒X = 兺共15兲

A similar expression exists for TX⇒Y . Transfer entropy represents the information about a future observation of variable
xt obtained from the simultaneous observation of past valuesof both xt and y t, after discarding the information about the
future of xt obtained from the past of xt alone. If not statedotherwise, in the subsequent part of this paper k = l = 1. The
conditional probabilities of Eq. 共15兲 are calculated from the

共16兲

k=i

where 0 ⱕ i ⱕ j ⬍ m, Rk ⱖ 0, and 1 ⱕ Ri + ¯ + R j ⱕ R. By arranging all possible x̃t in an array, we can define a nonlinear
extension of the time-delayed embedding vector, subsequently called the polynomial embedding vector 共PEV兲:R
x̃t  共xt, . . . ,xt−m+1,x2t ,xtxt−1, . . . ,xt−m+1
兲T .

共17兲

B. PEV-extended Granger causality

The PEV can be used to extend Granger causality. One
naive formulation is to replace the conventional delayed embedding vectors xt and yt with the PEVs x̃t and ỹt, that is,

共14兲

where p denotes the transition 共conditional兲 probability. If
the deviation is small, then we can safely assume that thestate of space Y has no 共or little兲 relevance on the transition
probabilities of the state vectors of space X. If the deviationis large, however, then the conditional independency assumption is not valid. The incorrectness of the assumption
can be quantified by the transfer entropy, which is formulated
as the Kullback-Leibler entropy between p共xt+1 兩 xkt 兲 andp共xt+1 兩 xkt , ylt兲,p共xt+1兩xkt ,ylt兲
.
p共xt+1,xkt ,ylt兲log
p共xt+1兩xkt 兲t

Rk,x̃t共i, j,R兲 = 兿 xt−k

xt+1 = ␣Tx̃t + ␲ + ␰共x兲
t ,

共18兲

y t+1 = ␤Tỹt + ␪ + ␰共y兲
t ,

共19兲

with ␣, ␲, ␤, and ␪ determined by autoregression so as to共y兲minimize the variance of ␰共x兲t and ␰t . The “mixed” autoregression equations are also modified according to,xt+1 = aTx̃t + bTỹt + p + ␰共x兩y兲
t

共20兲

.
y t+1 = cTx̃t + dTỹt + q + ␰共y兩x兲t

共21兲In this paper, however, we consider the polynomials of themixed-state embedding vector zt as in 关8兴. The polynomiallyextended mixed-state embedding vector z̃t incorporates thehigher-order cross-correlation terms of xt and y t as well as
the elements of x̃t and ỹt 关8兴.
In this case, the mixed information regression would be,
xt+1 = aTz̃t + p + ␰共x兩y兲
t

共22兲

,
y t+1 = cTz̃t + q + ␰共y兩x兲t

共23兲

and the resulting causality indices are

036217-3

PHYSICAL REVIEW E 77, 036217 共2008兲

ISHIGURO et al.

共24兲

共y兩x兲
兲.
CX⇒Y = var共␰共y兲t 兲 − var共␰t

共25兲

(a)

(b)

0.6

T
0.2

Hereafter, we call this method PEV-extended Granger causality 共PEV-GC兲.
C. Computational properties

IV. NUMERICAL EXPERIMENTS

To evaluate quantitatively the behavior of the PEV-GC
and to illustrate its effectiveness, we conducted numericalexperiments with two noiseless and two noisy model systems
with built-in causality structure. In all experiments we compared the proposed measure with transfer entropy 共TE兲, nonlinear Granger causality, and extended Granger causality. All
simulations were done in MATLAB .

0.4
0

0
0 0.2 0.4 0.6 0.8
e

(c)

1

(d)

0.16

δ

0.12

c

In this section, we describe the computational propertiesof this polynomial approach. Compared to existing nonlinear
extensions of Granger causality this method has two advantages. First, the only parameter is the order R of the largest
polynomial. It follows that the search for parameters is simplified in relation to the EGCI and NLGC, which both depend on the choice of at least two critical parameters. Second, due to its close relationship to the “classical” delayed
embedding vector, the PEV is straightforward to implementand does not require any additional computational effort suchas fuzzy c means clustering, evaluation of RBFs 共as in
NLGC兲, or search of local neighborhoods 共as in the EGCI兲.
To compute the PEV of embedding dimension m and order R, the size of the nonlinear embedding vector is the sumof the number of different monomials of degree l, that is,
l=R 1
兺l=1
l! m共m + 1兲 ¯ 共m + l − 1兲, where m is the size of the original linear embedding vector. For example, in our experiments we used 共m , R兲 = 共2 , 2兲 or 共m , R兲 = 共3 , 3兲 which correspond to a five- and a 19-dimensional vector, respectively
共the linear vector is two and three dimensional, respectively兲.
It is clear that by using a polynomial embedding vector thecomputational cost due to the evaluation of the autoregressive equations increases. However, it is still reasonably small
if compared to the additional operations required by the
NLGC and EGCI.

1.2
0.8

0.4

C

共x兩y兲
CY⇒X = var共␰共x兲兲,t 兲 − var共␰t

0.08
0.04
0

0 0.2 0.4 0.6 0.8
e1

0 0.2 0.4 0.6 0.8
e

1

2
1
0
−1
−2
−3
0

0.2 0.4 0.6 0.8
e1

FIG. 1. Unidirectionally coupled Hénon maps 共N = 103; nonoise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger Causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table I
for parameters.

By construction the two time series are unidirectionallycoupled with time series xt influencing time series y t. It follows that the causal influence of xt on y t should increase if
the strength of the coupling e increases. On the other hand,causal influence of y t on xt should be negligible because
there is no interaction term in Eq. 共26兲. From other studies
共e.g., 关3兴兲, it is known that xt and y t will enter a 共generally兲
synchronized state for e ⬎ 0.7. In that area, the two maps areindistinguishable from each other and the flow of information is zero 共despite a nonzero coupling between the maps兲.
The results of our experiments are shown in Fig. 1 共N
= 103兲 and in Fig. 2 共N = 104兲. The parameter settings used forthe analysis are summarized in Table I. Those parameters areexperimentally selected so as to give the best results. As
evident from the graphs, PEV-GC gives noise-free and consistent results even for small sample sizes 共N = 103兲. We also
observe that only PEV-GC successfully indicates absence ofa causal influence of y t on xt 共CX⇒Y = 0兲 at any e.

A. Noiseless polynomial chaotic map

We first studied a model system consisting of unidirectionally coupled noiseless chaotic 共Hénon兲 maps 关13兴:2
xt+2 = 1.4 − xt+1
+ 0.3xt ,共26兲

y t+2 = 1.4 − 关ext+1 + 共1 − e兲y t+1兴y t+1 + 0.3y t ,共27兲

where the parameter e regulates the coupling strength between the two maps. We varied e from 0 to 1 with increments of 0.025. For each coupling and for sample sizes of
N = 103 and 104, the time series xt and y t were generated 共theinitial 105 samples were discarded as transients兲. Each experiment was repeated five times.

B. Noisy polynomial chaotic map

In order to assess the robustness of PEV-GC against noise,in the second set of experiments we considered noiseaffected unidirectionally coupled Hénon maps. After generating 105 + N samples 共as in the previous experiment, the first
105 samples were discarded as transients兲, independent identically distributed white noise was added to the time series.
We tested the causality measures for N = 104 and noise levelsof 10%, 30%, and 50%. Each experiment was conducted fivetimes.
The results are presented in Fig. 3 共noise level 10%兲, Fig.
4 共noise level 30%兲, and Fig. 5 共noise level 50%兲. The pa-

036217-4

PHYSICAL REVIEW E 77, 036217 共2008兲

COMPARISON OF NONLINEAR GRANGER CAUSALITY…

0 0.2 0.4 0.6 0.8
e

0

1

0

δ
0 0.2 0.4 0.6 0.8
e

1

T
0

0.2 0.4 0.6 0.8
e

0 0 0.2 0.4 0.6 0.8 1
e

1

(d)

10

0.12

0

0.08

0.4
0

0.2
0.1

(c)

0.2

0.04

0.2

1

0.6

0.08

0

0 0.2 0.4 0.6 0.8
e

0.8

0.12

0.3

0.1

(d)

0.16

c

0.4
0.2

0

0.4

0.3

C

T

C
0.2

(b)

0.4

0.6

0.4

(c)

(a)

0.8

δ

(b)

0.6

c

(a)

0.04
0 0.2 0.4 0.6 0.8
e−20

0 0 0.2 0.4 0.6 0.8
e

1

−10

1

−300

0.2 0.4 0.6 0.8
e1

FIG. 2. Unidirectionally coupled Hénon maps 共N = 104, nonoise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table I
for parameters.

FIG. 4. Unidirectionally coupled Hénon map 共N = 104; 30%
noise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table I
for parameters.

rameter settings for each method are shown in Table I. We
observe that the PEV-GC and TE methods are robust against
the influence of fairly high levels of noise for all values of
the coupling parameter e. PEV-GC performs better because it
successfully shows that there is no information flow between
the two maps. NLGC, on the other hand, manages to disclose
the unidirectional causal dependency between the two maps

only for values of e 僆 关0.5, 0.7兴. Finally, the EGCI, being
sensitive to both sample size and noise level, fails to give asatisfying result. From these results, we can safely conclude
that only PEV-GC and TE give consistent results for thisnoisy chaotic map data.

(b)

0.4

0.4

0

(c)

0.2 0.4 0.6 0.8
e

(d)

C

0 0.2 0.4 0.6 0.8
e

01−0.4 0

0.1

0.06
0

0.2 0.4 0.6 0.8
e1c

0 0 0.2 0.4 0.6 0.8
e1

FIG. 3. Unidirectionally coupled Hénon maps 共N = 104; 10%noise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table I
for parameters.

1

0
−20

0.02

0.2 0.4 0.6 0.8
e0.02 0 0.2 0.4 0.6 0.8
e

(d)

20

0.04

−0.2
0 0.2 0.4 0.6 0.8
e

0.14

0.06

δ
0

0.18

0.2

(c)

0

0.1

0.22

0.1

1

0.2

0.2

c

0

1

(b)
0.3

0.2

0.2
0

(a)

T

0.6

δ

0.6

T

C

(a)

1

−40 0 0.2 0.4 0.6 0.8
e1

FIG. 5. Unidirectionally coupled Hénon maps 共N = 104; 50%
noise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table I
for parameters.

036217-5

PHYSICAL REVIEW E 77, 036217 共2008兲

ISHIGURO et al.
C. Noiseless sine-based chaotic map

(b)

(a)

In the third numerical experiment, we applied all fourmeasures to coupled chaotic maps featuring nonpolynomial
components:2
y t+2 = 1.7 − 0.4y t+1+ 1.1 sin共y t − ext兲,

共29兲

where the parameter e 僆 关0 , 1兴 denotes the coupling strengthof the causal interaction from X to Y. The embedding of thecoupling parameter in the phase term of the sine function
makes the detection of the causal relationship more challenging. The parameter e was incremented by 0.025 steps, andfor each e, we collected N = 103 and 104 samples after discarding 105 samples as transients. For each experiment, we
performed five runs.
By construction the two time series were unidirectionallycoupled with time series xt influencing time series y t. It follows that xt causally affects y t 共the strength of the causal
interaction should increase if the coupling parameter e increases兲. On the other hand, the causal influence of y t on xt
should be negligible due to the absence of an interaction
term in Eq. 共28兲.
The results are presented in Fig. 6 共N = 103兲 and in Fig. 7
共N = 104兲. The parameters used for the analysis are summarized in Table II. Those parameters are experimentally selected so as to gave the best results. In this experiment, allmethods except the EGCI fail to measure the unidirectionalcausal structure in the weak coupling condition 共e ⬍ 0.5兲. Another rather surprising result is that TE is unstable in the
weak coupling area. This result differs from the previous two
experiments, in which TE performed quite well.

T

C

0.4
0.2

0.1

0

0
0

0.2 0.4 0.6 0.8
e−0.2 0 0.2 0.4 0.6 0.8 1
1
e

(c)

(d)

1.2

0.08

0.8

0.04

δ

共28兲

0.6

0.2

0.4

c

2
xt+2 = 1.7 − 0.4xt+1+ 1.1 sin共xt兲,

0.8

0.3

0

0
0

0.2 0.4 0.6 0.8
e1

0 0.2 0.4 0.6 0.8
e1

FIG. 6. Unidirectionally coupled sine-based maps 共N = 103; nonoise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table IIfor parameters.

not even affected by heavy noise 共30%兲. On the other hand,a deterioration of performance is notable in the TE and
EGCI. Especially for the TE, the results are far from thetheoretical predictions; for weaker couplings, the estimatedcoupling direction is even inverted. NLGC displays a strange
behavior for values of e close to 0.6 where it suddenly falls
and rises again.

D. Noisy sine-based chaotic map

In the fourth experiment, we studied a noisy sine-based
map. As in the case of the second experiment, independent
identically distributed white noise was added to the extracted
bivariate time series. We tested the causality measures for
N = 104 and noise levels of 10%, 20%, and 30%; each experiment was conducted five times.
The results of our numerical experiments are presented inFig. 8 共noise level 10%兲, Fig. 9 共noise level 20%兲, and Fig.
10 共noise level 30%兲. The parameters for each method arepresented in Table II. As evident from the graphs, PEV-GC is
robust against noise, and the estimated coupling strengths are
TABLE I. Parameters for experiments on unidirectionallycoupled Hénon maps. If there are two values, they correspond to thecases of N = 103 , 104, respectively. The embedding dimension ism = 2 except for TE.
MethodParameters

Transfer entropy 共TE兲Extended GC index 共EGCI兲
Nonlinear GC 共NLGC兲
PEV-GC

r=8, m=1
L = 50, d = 0.5, 0.3
P = 50, ␴ = 0.05
R=2

(b)
0.3

0.6

0.2

0.4

T

C

(a)

0.2

0.1

0

0
0 0.2 0.4 0.6 0.8
e

(c)

1

−0.2

(d)

0 0.2 0.4 0.6 0.8
e

1

0.8

0.08

0.6

c

δ

0.04

0.4
0.2

0
0 0.2 0.4 0.6 0.8
e1

0
0 0.2 0.4 0.6 0.8 1eFIG. 7. Unidirectionally coupled sine-based maps 共N = 104, nonoise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table IIfor parameters.

036217-6

PHYSICAL REVIEW E 77, 036217 共2008兲

COMPARISON OF NONLINEAR GRANGER CAUSALITY…
TABLE II. Parameters for experiments on unidirectionallycoupled sine-based maps. If there are two values, they correspondto the cases of N = 103 , 104, respectively. The embedding dimension is m = 3 except for TE.

(b)

(a)
0.4

0.3

0.3

T

C

0.2

0.1

0.1

Parameters

Transfer entropy 共TE兲Extended GC index 共EGCI兲
Nonlinear GC 共NLGC兲
PEV-GC

r=8, m=1
L = 50, d = 0.8, 0.6
P = 50, ␴ = 1.0
R=3

0

(c)

0

0.16

0.2 0.4 0.6 0.8
e

1

0

(d)

0.12

c

0.04
0

E. Exponential function coupling data

25xt
4共x2t + 1兲

+ 2 cos共1.2t兲 + N共0,1兲, 共30兲

y t+1 = 0.1y 2t − e关f共xt兲2 − 0.3兴 + N共0,1兲,

共31兲

where
1.0
f共xt兲 =
1.0 + exp共− xt兲

共32兲

and N共0 , 1兲 is a random value from a normal distribution ofzero mean and unit variance. One big difference from previous chaotic maps is that X drives Y through an exponential
共sigmoid兲 function. Apparently, exponential bases of NLGC
will suit for estimating the X ⇒ Y causal couplings. Note also

C

T

C

0.1
0 0.2 0.4 0.6 0.8 1
e

0.12

0

(d)

δ

c
0.04

0 0.2 0.4 0.6 0.8
e

0.2 0.4 0.6 0.8
e

1

0.3

0.3

0.2

0

1

(c)

0.1
0 0.2 0.4 0.6 0.8
e

0

1

0

0.12

−10

0.08
0.04

0

0.2 0.4 0.6 0.8
e01

FIG. 8. Unidirectionally coupled sine-based maps 共N = 104; 10%noise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table IIfor parameters.

1

10

0.16

0.1

0 0.2 0.4 0.6 0.8
e

(d)

0.2

0.3

0
0

0.4

0.1

0.2

0.08

(b)

0.5

0.2

0.2

0

0

(a)

0.4

0.2

(c)

that the noise term is incorporated in each dynamic model.
The results are presented in Fig. 11 for N = 104. The parameters for each method are presented in Table III. Thoseparameters are experimentally selected so as to gave the best
results. The NLGC presents excellent estimation results aspredicted, reflecting the coupling structure of the data.
Though the bases are not polynomials, the PEV-GC still performs quite well.

0.6

0.3

1

FIG. 9. Unidirectionally coupled sine-based maps 共N = 104; 20%
noise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table IIfor parameters.

(b)

0.4

1

−2
−4

0 0.2 0.4 0.6 0.8
e

c

(a)

0 0.2 0.4 0.6 0.8e2

δxt+1 = 0.125xt +1

TSo far, we have examined the performances of causality
measures on several chaotic maps. As you can see, PEV-GC
shows preferable results against the others.
In this experiment, we use the following dynamic models
for data generation:

0 0.2 0.4 0.6 0.8
e

0

0.08

δ

Method0.2−20
0

0.2 0.4 0.6 0.8
e1

0 0.2 0.4 0.6 0.8 1eFIG. 10. Unidirectionally coupled sine-based maps 共N = 104;
30% noise兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy. 共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality index. The horizontal axis denotes the strength of coupling
and the vertical axis denotes the causality measure. Thick lines
represent X ⇒ Y causality and thin lines represent Y ⇒ X causality.
The plots are averages over five experimental runs. The error bars
denote standard deviations of causality measures for five runs. See
Table II for parameters.

036217-7

PHYSICAL REVIEW E 77, 036217 共2008兲

ISHIGURO et al.

(a)

(b)
0.12

0.12

0.08

C

T

0.08
0.04
0

(c)

0.04
0 0.2 0.4 0.6 0.8 1
e

0.16
0.08
0.04
0

(d)

0 0.2 0.4 0.6 0.8 1
e
-3

x10

20
15
10
5
0
−5
0 0.2 0.4 0.6 0.8 1
e

δ

C

0.12

0

0 0.2 0.4 0.6 0.8 1eFIG. 11. Unidirectional coupling data with exponential function
共N = 104兲. 共a兲 PEV-extended Granger causality. 共b兲 Transfer entropy.
共c兲 Nonlinear Granger causality. 共d兲 Extended Granger causality
index. The horizontal axis denotes the strength of coupling and thevertical axis denotes the causality measure. Thick lines represent
X ⇒ Y causality and thin lines represent Y ⇒ X causality. The plots
are averages over five experimental runs. The error bars denotestandard deviations of causality measures for five runs. See Table
III for parameters.
V. DISCUSSIONAs the experimental results show, the PEV-GC successfully detects the underlying causal relationship in manycases, regardless of the number of samples and the heavinessof noise contamination. One possible reason for this is that
the PEV may be interpreted as a truncated Taylor expansionof order R. This view is interesting because any elementary
and well-behaved function can be expanded into polynomials
by Taylor expansion. This indicates that the PEV is potentially applicable to many different kinds of bivariate timeseries.
It is also appealing for users that the PEV is easy to
implement and has only one control parameter 共R兲, compared to other alternatives presented in this paper. Since tuning of multiple parameters is a difficult task, the smaller
number of parameters is attractive in practice.
We also point out a substantial difference between thepolynomial embedding vector and the “kernel method” usedin nonlinear Granger causality 关1兴. Kernel methods are rather
popular nonlinear techniques in statistical pattern recognition
共e.g., 关14兴兲. They rely on the use of “kernel functions,” deTABLE III. Parameters for experiments on data unidirectionallycoupled through exponential function. The embedding dimension ism = 3 except for TE.
MethodParameters

Transfer entropy 共TE兲Extended GC index 共EGCI兲
Nonlinear GC 共NLGC兲
PEV-GC

r=8, m=1
L = 100, d = 2.0
P = 100, ␴ = 1.0
R=3signed as nonlinear distance functions 共e.g., Gaussian, sigmoid兲 between samples. One characteristic of kernel techniques is that the size of the nonlinear 共embedding兲 vector isvery dependent on the number of samples, not on the dimensionality of each sample. On the other hand, the size of the
PEV is basically dependent on the dimension of a sample
共embedding dimension兲, not on the sample size.
This point makes a difference between the two techniquesin computational costs. We can see this by comparing two
research tasks: image processing and time series analysis. In
image processing, an image has very high dimensionality
共e.g., even a very small 64⫻ 64 pixel image has 4096 pixels兲
and the number of samples 共images兲 is rather small 共in manycases, some hundreds兲. Clearly, this is not an ideal setting for
the PEV—small number of samples for high-dimensional
data. Next, let us consider time series analysis. Time series
analysis techniques need to handle large sample pools 共N⬎ 103兲 for which, generally speaking, the embedding dimension is not too large 共in many cases, m ⬍ 10兲. In this case,kernel techniques may face a computational problem—manysamples for low-dimensional data. The number of samples
can be reduced by introducing a heuristic trick such as clustering preprocessing in nonlinear Granger causality 关1兴, but it
may lead to a deterioration of performance unless the heuristics properly reflect the nature of the data. Therefore, intime series analysis we have two reasons to use the PEV: one
is computational advantage and the other is that the PEV
does not need heuristics which may deteriorate the performance. As the results show, the PEV-GC performs efficiently
in the low-dimensional data set. However, we should notethat this may not be true for highly nonlinear or highdimensional data.
The main rationale behind using a nonlinear embedding
vector in combination with Granger causality was to extend
the traditional linear analysis and take into account also nonlinear structure. In principle it should be possible to extend
other measures of causality also, such as transfer entropy,
using the polynomial embedding vector. One caveat is thatits formulation could become rather complex. Specifically, inthe case of transfer entropy, we would need to estimate theprobability density function of the time series formed by theembedding vectors. If the PEV is used, then such estimation
is rather difficult and requires a substantial computational
effort and a very large data pool 共e.g., for m = 3 , R = 3 thePEV is a 19-dimensional vector for which we would have toestimate probability distribution functions with up to 19 dimensions兲.
VI. CONCLUSION

In this paper, we survey several nonlinear extensions ofGranger causality and assess their performances through numerical experiments using a low-dimensional and low-ordernonlinearity data set. In particular, we focus on the extensionof the time-delayed embedding vector 共the polynomial embedding vector兲 which in combination with Granger causality leads to a computationally efficient nonlinear extension ofGranger causality. We studied the effectiveness of this typeof extension in numerical experiments with some artificial

036217-8

PHYSICAL REVIEW E 77, 036217 共2008兲

COMPARISON OF NONLINEAR GRANGER CAUSALITY…

chaotic maps. Intensive experiments concerning the type ofthe maps, the number of samples, and the heaviness of noisecontamination were conducted, and the results showed thatthe polynomial embedding vector based technique is superior
or comparable to other alternative measures in quantifying

the asymmetric causal interdependence between bivariate
time series data, at least in low-complexity systems. In future
work, we will apply the proposed PEV-GC to data gatheredfrom real world systems such as climate data, neurophysiological signals, and sensory and motor variables of robots.

关1兴 N. Ancona, D. Marinazzo, and S. Stramaglia, Phys. Rev. E 70,
056221 共2004兲.
关2兴 J. Arnhold, P. Grassberger, K. Lehnertz, and C. E. Elger,
Physica D 134, 419 共1999兲.
关3兴 T. Schreiber, Phys. Rev. Lett. 85, 461 共2000兲.
关4兴 Y. Chen, G. Rangarajan, J. Feng, and M. Ding, Phys. Lett. A
324, 26 共2004兲.
关5兴 U. Feldmann and J. Bhattacharya, Int. J. Bifurcation Chaos
Appl. Sci. Eng. 14, 505 共2004兲.
关6兴 M. Lungarella, K. Ishiguro, Y. Kuniyoshi, and N. Otsu, Int. J.
Bifurcation Chaos Appl. Sci. Eng. 17, 903 共2007兲.
关7兴 C. W. J. Granger, Econometrica 37, 424 共1969兲.
关8兴 I. I. Mokhov and D. A. Smirnov, Geophys. Res. Lett. 33,

L03708 共2006兲.
关9兴 N. H. Packard, J. P. Crutchfield, J. D. Farmer, and R. S. Shaw,
Phys. Rev. Lett. 45, 712 共1980兲.
关10兴 M. Takens, in Dynamical Systems and Turbulence, edited by
D. Rand and L. Young, Lecture Notes in Mathematics Vol. 898
共Springer Verlag, Berlin, 1981兲, pp. 366–381.
关11兴 M. B. Kennel, R. Brown, and H. D. I. Abarbanel, Phys. Rev. A
45, 3403 共1992兲.
关12兴 A. M. Fraser and H. L. Swinney, Phys. Rev. A 33, 1134
共1986兲.
关13兴 M. Henon, Commun. Math. Phys. 50, 69 共1976兲.
关14兴 V. N. Vapnik, The Nature of Statistical Learning Theory
共Springer, Berlin, 1995兲.

036217-9