PHYSICAL REVIEW E 77, 026216 共2008兲

Detecting direction of causal interactions between dynamically coupled signals
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 February 2008兲

The problem of temporal localization and directional mapping of the dynamic interdependencies between
parts of a complex system is addressed. We present a technique that weights the sampled values so as tominimize the mutual prediction error between pairs of measured signals. The reliability of the detected intermittent causal interactions is maximized by 共a兲 smoothing the weight landscape through regularization, and 共b兲
using a nonlinear 共polynomial兲 variant of the conventional embedding vector. The effectiveness of the proposed technique is demonstrated by studying three numerical examples of dynamically coupled chaotic maps
and by comparing it with two other measures of causal dependency.
DOI: 10.1103/PhysRevE.77.026216

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

The identification and quantification of hidden interdependencies between the parts of a complex system promotes a
deeper understanding of the system’s dynamics. One way oftackling this important but challenging issue is by estimating
to what extent a number of carefully selected observables
共e.g., pairs of measured signals兲 are mutually predictive ofeach other. Causality measures based on this idea have been
proposed which quantify the strength and the direction of thecausal interaction between bivariate time series 共e.g., 关1–9兴兲.
The basic assumption underlying all these methods is that thecausal relationship does not change too rapidly over time.
Signals sampled from the real world, however, are rarely
stationary and well behaved, and couplings typically appear,
disappear, and reappear over time, becoming weaker or
growing stronger 共e.g., causal interactions between the stock
market exchanges of Tokyo and New York兲. Here, we address the question of how to detect and localize the directionof such temporally intermittent couplings, that is, of how to
detect the waxing and waning of causal interactions.
A naive approach is to repeatedly estimate causal relationships within a moving data window of fixed size. In thistechnique, the moving window slides across the data set inincrements of the data interval, and the causal analysis isconducted using the samples contained in the window. The
major problem of this approach is how to find an appropriate
window size. Too wide a window decreases the temporal
resolution of the analysis; too narrow a window reduces thestatistical reliability. A less naive technique recently proposed 关10兴 relies on embeddings across multiple realizationsof time-dependent dynamics of a system. This technique isquite powerful and can resolve causal relations based on only
200 samples, its major drawback being that it requires atleast a few tens of realizations of sample sequences, independently generated from the analyzed system.

The purpose of this paper is to motivate and derive analternative method for measuring changes in the couplings
between interacting systems. To this end, we introduce thenotion of a causality marker: at each time step the causal
interaction is weighted based on its strength. For nonstationary data sets, such weighting is unequal and changes overtime as a result of changes in the temporal causal relation
between the data. The concept of the causality marker isgeneral and potentially applicable to various existing methods used to detect causal structure 关1,5,7,11兴. In this paper,
we focus on its formulation in the context of Granger causality 关4兴.
For simplicity, we consider the bivariate case. Let X= 兵x1 , x2 , . . . , xN其 and Y = 兵y 1 , y 2 , . . . , y N其 denote two simultaneously measured time series each consisting of N scalarquantities. As in most previous approaches, we start by reconstructing the state spaces corresponding to both time series by the method of time delay embedding, in which thestate 共or embedding兲 vectors are formed by delayed pastscalar observations of the time series 关12兴. The embedding vector of dimension m of X is xm,dt
= 共xt , xt−d , xt−2d , . . . , xt−共m−1兲d兲T, where d is the time delay 共orlag兲 between successive elements of the vector. If not otherwise stated, in this paper d = 1, and we use the expression
xt = 共xt , xt−1 , xt−2 , . . . , xt−m+1兲T. Similarly, we define yt. For the
choice of adequate values for m and d, we refer to 关13,14兴.
The notion of Granger causality between two time series
is based on the idea that the prediction of one time series
through linear autoregression can be improved by incorporating information about the past values of the other timeseries. To estimate the causal effect of Y ⇒ X 共dependency of
X on Y兲, the two following autoregression equations areused:*Corresponding author. Present address: NTT Communication

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

共1兲

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

共2兲

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

where the terms ␣, ␲, a, b, and p are determined so as to共兩x兩y兲
,minimize the variance of the prediction errors ⑀共x兲t and ⑀t

026216-1©2008 The American Physical Society

PHYSICAL REVIEW E 77, 026216 共2008兲

ISHIGURO et al.

min

1
共x兲 2
兺 共⑀共x兲t −␩ 兲 ,␣,␲ N t共3兲1
共⑀共兩x兩y兲− ␩共兩x兩y兲兲2 ,
兺
t
N
a,b,p
t

共4兲

min

with ␩ denoting the time-averaged error. After calculating
the coefficients, we can express the Granger causality as兲
var共⑀共兩x兩y兲
t
CY⇒X = 1 −共x兲 = 1 −var共⑀t 兲− ␩共兩x兩y兲兲2
兺t 共⑀共兩x兩y兲
t共x兲 2
兺t 共⑀共x兲t −␩ 兲

.

Intuitively, the larger CY⇒X the more X is causally dependent
on Y. The dependence in the opposite direction, CX⇒Y , canbe defined in a similar way.
The idea of the causality marker is based on the assumption that, for improving the detection of causal relationships,
the analysis should place more emphasis 共weight兲 on observables having a larger causal influence on other observables.
This weight, called the causality marker, should be large if a
causal dependency exists, but should take a small value otherwise. In other words, causality markers indicate 共“mark”兲
the locations of potential causal relationships.
Given two bivariate time series X and Y consisting of
N samples, we define as causality markers the two
其 and WY⇒X
N-dimensional weight vectors WX⇒Y = 兵wx→y
t
y→x
x→y
= 兵wty→x其 for which 兺twx→y
=
兺
w
=
N,
w
ⱖ 0, and
t t
t
t
y→x
wt ⱖ 0. Using these weight vectors, we can reformulate
Eqs. 共3兲 and 共4兲 as
min

1
共x兲 2
兺 wty→x共⑀共x兲
t −␩ 兲 ,␣,␲ N t共5兲1
− ␩共兩x兩y兲兲2 .
兺 wty→x共⑀共兩x兩y兲ta,b,p N t

共6兲

min

其 and 兵wty→x其What one seeks are two sets of weights 兵wx→ytthat maximize CX⇒Y and CY⇒X. These weights can be calculated by substituting the results of Eqs. 共5兲 and 共6兲 in thefollowing “weighted” Granger causality measure:冢max CY⇒X = max 1 −

WY⇒X

WY⇒X

− ␩共兩x兩y兲
兲2
兺t wty→x共⑀共兩x兩y兲
t
t
共x兲 2
兺t wty→x共⑀共x兲
t − ␩t 兲冣. 共7兲

Analogous expressions exist for CX⇒Y and WX⇒Y .
Although this maximization problem cannot be solved
analytically, it is possible to analytically derive the partial
derivatives of C with respect to W, and numerically determine the local maximum through the method of steepest descent.
Our method has several interesting properties. First, it
does not rely on parameters such as window size or theamount of overlap between adjacent windows, as is the casefor the moving window technique. Second, the number ofprocessed samples is not limited to the ones contained in aparticular window centered at a specific time step. The absence of windowing enables us to avoid the dilemma between statistical reliability and granularity of analysis causedby an inadequate window size.
One major drawback of the analysis based on causality
markers is the irregularities 共roughness兲 originating from theapplication of the steepest descent algorithm 共which computes a gradient兲. The causality marker is calculated locally
based on the improvement of the regression at each time
step, making the estimation dependent on the quality of theobserved values. This property is not desirable if the method
is applied to data sampled from chaotic systems, which in thegeneral case can display rather complex dynamics, or to datacollected from real world systems, which are typically contaminated by noise or inconsistencies. An interesting denoising technique for noisy bivariate data sets relying on Kalmanfiltering was recently proposed in 关15兴. This technique, however, assumes that the time series are generated from a system with stationary dynamics, and not characterized by temporally intermittent couplings.
We cope with both issues by penalizing Eq. 共7兲 with aregularization term冢1Nmax CY⇒X + ␭

WY⇒X

where
l

冉

兺t 共wty→x − w̄兲2
1
N兺t ␴2共t兲
l

冣

1
1
y→x
␴ 共t兲 =wt+i−兺
兺 wy→x2l + 1 i=−l2l + 1 j=−l t+j
2

共8兲

,

冊

2

,

w̄ is the average of the causality markers over all time steps,␭ ⱖ 0 controls the effectiveness of the regularizing second
term, and l is the size of the window used to calculate the“local” variance ␴2共t兲 共in all our experiments, l = 2兲. To maximize Eq. 共8兲, it is necessary to maximize CY⇒X as well as its
second term. To achieve this, the denominator ␴2共t兲 should
be small and the numerator 共“global” variance兲 should belarge. By satisfying these two requirements, the causality
marker is “smoothed” and the difference of coupling strength
between the weakly and strongly causally related periods iswidened, keeping the global variance large. In this paper, we
optimize ␭ to yield the best results in our experiments. Thereis no decisive approach, however, for the choice of ␭ if novalidation information is provided.
To more effectively handle data with nonlinear dynamics,
we use a polynomial embedding vector 共PEV兲 关11,16兴. The
PEV consists of the polynomial combinations of the monoRkjxt−k, where 0
mial components of xt, namely, x̃t共i , j , R兲 = 兿k=ij
ⱕ i ⱕ j ⬍ m, Rk ⱖ 0, and 1 ⱕ 兺k=iRk ⱕ R. By arranging allx̃t共i , j , R兲 in an array, we can define the PEV asRx̃t  共xt , . . . , xt−m+1 , x2t , xtxt−1 , . . . , xt−m+1
兲T. We define ỹt analogously. PEV-based Granger causality 共PEVGC兲 extends
Granger causality by replacing xt with x̃t and yt with ỹt in
Eqs. 共1兲 and 共2兲. Instead of xt and yt, in Eq. 共2兲 we can also
employ a joint-state PEV z̃t, derived from zt = 共xTt yTt 兲T. Such a
joint-state PEV has been recently applied to the analysis of026216-2

DETECTING DIRECTION OF CAUSAL INTERACTIONS …

PHYSICAL REVIEW E 77, 026216 共2008兲

TABLE I. Coupling strength parameters for artificial data
sets.
Sine-based

eyx

exy

eyx

exy

0.0
0.0
0.5
0.0
0.0

0.0
0.5
0.5
0.0
0.25

0.0
0.0
0.7
0.0
0.35

0.0
0.7
0.0
0.0
0.7

C

0–300
300–600
600–900
900–1200
1200–1500

Hénon

C

Time step5

the dynamics of climatic modes 关16兴. Selecting R is also notan easy problem, but in general a small value of R 共0–3兲 issufficient and preferable. There are two reasons for choosing
a smaller value of R: one is that higher-order correlations are
not effective for regression in many cases because of thevolatility, and the other is that larger R reduces computational efficiency exponentially.
To illustrate the performance of the proposed technique,
we analyzed the bivariate time series extracted from three
dynamical systems and compared the results with two otherapproaches: 共i兲 a combination of a moving window analysis
and naive 共linear兲 Granger causality; 共ii兲 a combination of amoving window analysis and PEVGC. In what follows, wewill show that the causality-marker-based method outperforms the other two.
The first system of interest was composed of two bidirectionally coupled identical Hénon maps

5

FIG. 1. Expected 共true兲 and estimated coupling directions 共causality markers兲 for coupled Hénon maps. The horizontal axis denotes time; the vertical axis represents the value of the markers
WX⇒Y 共left兲 and WY⇒X 共right兲. The dotted lines represent the threshold values. 共a兲,共b兲 Expected 共true兲 coupling: 共a兲 X ⇒ Y and 共b兲
Y ⇒ X. 共c兲,共d兲 Estimated coupling without regularization 共␭ = 0.0兲:
共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Estimated coupling with regularization 共␭ = 10.0兲: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

FIG. 2. Moving window analysis of temporal causality measuresfor coupled Hénon maps. The horizontal axis denotes time; thevertical axis represents the value of the measures CX⇒Y 共left兲 andCY⇒X 共right兲. The dotted lines represent the threshold values. 共a兲,共b兲
Expected 共true兲 coupling borrowed from Fig. 1: 共a兲 X ⇒ Y and 共b兲
Y ⇒ X. 共c兲,共d兲 Estimated causality measures with naive Granger
causality: 共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Estimated causality measures with PEVGC: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

xt+1 = 1.4 − 关eyxy t + 共1 − eyx兲xt兴xt + 0.3xt−1 ,y t+1 = 1.4 − 关exyxt + 共1 − exy兲y t兴y t + 0.3y t−1 ,where exy and eyx denote the coupling strengths of xt on y t
and y t on xt, respectively. After discarding 104 samples astransients, we collected N = 1500 samples while the couplingstrengths were altered according to Table I. Due to the sudden changes of the coupling constants, the extracted bivariate
data were rendered manifestly nonstationary. The embedding
dimension m and order R of the PEV were optimized insample to yield the best result, and were set to m = 3 and R= 2, respectively. All causality markers were initialized to 1,
and the steepest descent algorithm maximizing CX⇒Y andCY⇒X was iterated 20 times to yield the estimated couplings
between the time series. To eliminate edge effects, we cut off
the last few causality marker values. We show the results in
Fig. 1. As shown in Figs. 1共e兲 and 1共f兲, the regularized causality markers are sufficiently similar to the expected ones
关Figs. 1共a兲 and 1共b兲兴. The smoothing effect of the regularization term is also noteworthy.
For a comparison, we present the results of a moving
window analysis with naive and PEV Granger causality. To
estimate the amount of causal interaction, a moving windowof 300 samples was slid across the time axis one time step at
a time, and the causality measures were evaluated for every
window. The estimated couplings are shown in Fig. 2. As
evident from the figures, each method 关Figs. 2共c兲 and 2共d兲for naive Granger causality, and Figs. 2共e兲 and 2共f兲 for PEVGC兴 fails to detect the X ⇒ Y couplings in the interval
t = 600– 900. The limited sample size makes true directional
coupling Y ⇒ X difficult to estimate. We can also observe

026216-3

PHYSICAL REVIEW E 77, 026216 共2008兲

ISHIGURO et al.

X

-0.15

Activation
Repression-0.2

Controlled

A
0.1

-0.05

B

C

D

-0.05

-0.01

Y

FIG. 3. Expected 共true兲 and estimated coupling directions 共causality markers兲 for coupled sine-based chaotic maps. The horizontal
axis denotes time; the vertical axis represents the value of the markers WX⇒Y 共left兲 and WY⇒X 共right兲. The dotted lines represent thethreshold values. 共a兲,共b兲 Expected 共true兲 coupling: 共a兲 X ⇒ Y and 共b兲
Y ⇒ X. 共c兲,共d兲 Estimated coupling without regularization 共␭ = 0.0兲:
共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Estimated coupling with regularization 共␭ = 10.0兲: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

that, by contrast with the causality-marker-based method, theother two measures are characterized by more gentle trends.
In our second numerical experiment, we studied a sinebased chaotic map

C

C

(c)

FIG. 5. Studied gene regulatory network structure. The network
consists of six nodes; the time series are extracted from nodes X and
Y. The numbers labeling the lines are the strength of excitory 共⬎0兲or inhibitory 共⬍0兲 connections.

xt+1 = 1.7 − 0.4x2t + 1.1 sin共xt−1 − eyxy t−1兲,
y t+1 = 1.7 − 0.4y 2t + 1.1 sin共y t−1 − exyxt−1兲.
All parameters were the same as in the previous experiment,
except for the coupling strengths 共Table I兲 and the PEV order
共R = 3兲. The coupling constants 共exy and eyx兲 were varied according to the pattern shown in Figs. 3共a兲 and 3共b兲. As can beobserved in Fig. 3, despite the peakiness of the curves, thecausal dependency between the two maps is well reflected
for the X ⇒ Y relationship 关Figs. 3共c兲 and 3共e兲兴. For the caseof Y ⇒ X 关Figs. 3共d兲 and 3共f兲兴, the calculated causality markers for t = 1200– 1500 are relatively small. It is known thatfor this map for e ⬍ 0.5 the detection of the coupling is nontrivial 关11兴.
As in the previous experiment, here also we present theresults of a moving window analysis with naive and PEVGranger causality. The width of the moving window was 300
samples. The estimated coupling constants are shown in Fig.
4. As evident from the figure, the obtained causality measures are less similar to the true coupling status.
In our third experiment, we tested the proposed causality
marker with time series extracted from a genetic regulatory
network model 关17兴. The studied network is depicted in Fig.
5. The network has several nodes 共“genes”兲 linked byweighted connections. The connections are directional and
TABLE II. Coupling strength parameters for gene regulatory
network data.

FIG. 4. Moving window analysis of temporal causality measures
for sine-based chaotic maps. The horizontal axis denotes time; thevertical axis represents the value of the measures CX⇒Y 共left兲 andCY⇒X 共right兲. The dotted lines represent the threshold values. 共a兲,共b兲
Expected 共true兲 coupling borrowed from Fig. 3: 共a兲 X ⇒ Y and 共b兲
Y ⇒ X. 共c兲,共d兲 Estimated causality measures with naive Granger
causality: 共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Estimated causality measures with PEVGC: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

Time step

Y → X共Hxy兲

X → Y共Hyx兲

0–200
200–400
400–600
600–800
800–1000

−0.15
0.0
−0.15
0.0
−0.15

0.1
0.1
0.0
0.0
0.1

026216-4

DETECTING DIRECTION OF CAUSAL INTERACTIONS …

PHYSICAL REVIEW E 77, 026216 共2008兲

FIG. 6. Expected 共true兲 and estimated coupling directions 共causality markers兲 for gene regulatory network data. The horizontal
axis denotes time; the vertical axis represents the value of the markers: WX⇒Y 共left兲 and WY⇒X 共right兲. The dotted lines represent thethreshold values. 共a兲,共b兲 Expected 共true兲 coupling: 共a兲 X ⇒ Y and 共b兲
Y ⇒ X. 共c兲,共d兲 Estimated coupling without regularization 共␭ = 0.0兲:
共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Estimated coupling with regularization 共␭ = 10.0兲: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

CC

have weights as shown in the figure. All the connections are
inhibitory except the one between the nodes A and C, whichis excitatory. Specifically, in our experiments we controlled
two connection weights to realize temporal changes in thecausal relationships between the nodes X and Y. The connection weight between the nodes Y and X can be understood as

FIG. 7. Moving window analysis of temporal causality measures
for gene regulatory network data. The horizontal axis denotes time;the vertical axis represents the value of the measures CX⇒Y 共left兲
and CY⇒X 共right兲. The dotted lines represent the threshold values.
共a兲,共b兲 Expected 共true兲 coupling borrowed from Fig. 6: 共a兲 X ⇒ Y
and 共b兲 Y ⇒ X. 共c兲,共d兲 Estimated causality measures with naive
Granger causality: 共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Estimated causality measures with PEVGC: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

FIG. 8. Binary classified causality markers. Solid lines denotethe classification results of expected causality markers ŴX⇒Y 共left兲
and ŴY⇒X 共right兲. Dashed lines represent the classification resultsof estimated causality markers W̃X⇒Y 共left兲 and W̃Y⇒X 共right兲.
共a兲,共b兲 Results of analysis on Hénon map data 共␭ = 10.0兲: 共a兲 X ⇒ Y
and 共b兲 Y ⇒ X. 共c兲,共d兲 Results of analysis on sine-based chaotic map
data 共␭ = 10.0兲: 共c兲 X ⇒ Y and 共d兲 Y ⇒ X. 共e兲,共f兲 Results of analysis
on gene regulatory network data 共␭ = 10.0兲: 共e兲 X ⇒ Y and 共f兲 Y ⇒ X.

the strength of the Y ⇒ X causal relation; whereas the one
between the nodes A and C is used as an indicator of the
X ⇒ Y causal relationship. The X ⇒ Y link between A and C
does not directly connect X and Y. We thus assume that thecausal relationship in the X ⇒ Y direction is more difficult to
detect.
The expression levels of all M nodes composing the network are described as 关17兴 gt+1 = gt + H共gt − ⌫兲 + ⑀, where git isthe expression level of the ith gene, H is a connectivity
strength matrix whose elements are zero except for the designated connections in Fig. 5, ⌫ is a constant vector whose
elements are all set to 50, and ⑀ is a uniform noise drawn
from 关−10, 10兴. A ceiling function is used to truncate thevalues of git at 0 for the lowest and 100 for the largest, yielding signals with highly nonlinear properties. We collectedgene expression values for 1000 time steps, and estimatedthe causal relationship between the nodes X and Y. The coupling strengths Hxy and Hyx are presented in Table II. We
selected the embedding dimension m = 2 and R = 3. Experimental results are presented in Fig. 6 for the causality marker
TABLE III. Misclassification rates for Hénon map data set.
MethodX⇒YY ⇒X

Moving window+ naive GC
Moving window+ PEVGC
Causality marker+ PEVGC

0.21270.23000.02270.24870.17670.0647026216-5

PHYSICAL REVIEW E 77, 026216 共2008兲

ISHIGURO et al.
TABLE IV. Misclassification rates for sine-based map data
set.

TABLE
network.

V.

Misclassification

rates

for

gene

regulatory

MethodX⇒YY ⇒X

MethodX⇒YY ⇒X

Moving window+ naive GC
Moving window+ PEVGC
Causality marker+ PEVGC

0.38730.12070.0727

0.46270.3453
0.118

Moving window+ naive GC
Moving window+ PEVGC
Causality marker+ PEVGC

0.66600.43300.417

0.29100.32000.17

By comparing the classification results W̃X⇒Y and W̃Y⇒X withthe classification results of the expected causality markers
关Figs. 1共a兲, 1共b兲, 3共a兲, 3共b兲, 6共a兲, and 6共b兲兴, we can computethe misclassification rate with P = 兺t储w̃t − ŵt储 / N where ŵt is
the “true” classification result, and N is the number ofsamples. The classified causality markers for all data sets areshown in Fig. 8 and the misclassification rates are reproduced in Tables III–V. The misclassification rates for themoving window analysis are calculated analogously. At eachtime step 共window兲, the computed causality measure isthresholded by the average of causality measures of all thewindows. The comparison with the moving window analysis
reveals that the misclassificaation rates of causality markers
are quite satisfactory.

In conclusion, we addressed an important issue in the context of complex dynamical systems analysis, that is, the assessment of dynamic interdependencies between the constituent parts of a complex system based on measured
variables. To automatically localize and map the direction ofsuch dependencies, we proposed the notion of causality
marker, which evaluates the strength of the causal relationship at every time step rather than an average throughout alltime steps. The advantages of our approach are twofold: 共a兲
it is bottom up, that is, there is no need for a priori assumptions on the analyzed data or on the employed causality measure; 共b兲 it can be easily applied to all nonlinear extensionsof Granger causality—at least to the ones known to the authors 关1,7,11,18兴. Our experimental results show that our
method in combination with nonlinear extended Granger
causality and a simple classifier is able to detect changes in
the coupling between two interacting system with a hightemporal resolution. As one of the directions for further
study, we mention a formulation of the causality marker applicable to other measures of causal dependency, such asprobabilistic and information-theoretic approaches 共e.g., 关6兴兲.
We point out that, in contrast to the iterative 共batch兲 calculations of the steepest descent, an incremental variant of thecausality marker algorithm might be particularly beneficial
for very long time series data. We finally note that largedimensional or continuous-time problems remain an open
issue, and that our experiments serve only the purpose ofdemonstrating the potential of the proposed method.

关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兴 U. Feldmann and J. Bhattacharya, Int. J. Bifurcation Chaos
Appl. Sci. Eng. 14, 505 共2004兲.
关4兴 C. W. J. Granger, Econometrica 37, 424 共1969兲.
关5兴 M. Lungarella, K. Ishiguro, Y. Kuniyoshi, and N. Otsu, Int. J.
Bifurcation Chaos Appl. Aci. Eng. 17, 903 共2007兲.
关6兴 T. Schreiber, Phys. Rev. Lett. 85, 461 共2000兲.
关7兴 D. Marinazzo, M. Pellicoro, and S. Stramaglia, Phys. Rev. E
73, 066216 共2006兲.
关8兴 M. G. Rosenblum and A. S. Pikovsky, Phys. Rev. E 65,
041909 共2002兲.
关9兴 D. A. Smirnov and B. P. Bezruchko, Phys. Rev. E 68, 046209
共2003兲.
关10兴 R. Andrzejak, A. Ledberg, and G. Deco, New J. Phys. 8, 6
共2006兲.

关11兴 K. Ishiguro, N. Otsu, M. Lungarella, and Y. Kuniyoshi, Phys.
Rev. E 共to be published兲.
关12兴 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.
关13兴 M. B. Kennel, R. Brown, and H. D. I. Abarbanel, Phys. Rev. A
45, 3403 共1992兲.
关14兴 A. M. Fraser and H. L. Swinney, Phys. Rev. A 33, 1134
共1986兲.
关15兴 H. Nalatore, M. Ding, and G. Rangarajan, Phys. Rev. E 75,
031123 共2007兲.
关16兴 I. I. Mokhov and D. A. Smirnov, Geophys. Res. Lett. 33,
L03708 共2006兲.
关17兴 J. Yu, V. A. Smith, P. P. Wang, A. J. Hartemink, and E. D.
Jarvis, Bioinformatics 20, 3594 共2004兲.
关18兴 Y. Chen, G. Rangarajan, J. Feng, and M. Ding, Phys. Lett. A
324, 26 共2004兲.

based approach and in Fig. 7 for the moving window analysis with a window size of 200 samples. The detection of the
X ⇒ Y coupling is difficult and fails for all methods. However, the causality marker performs better in estimating theY ⇒ X causal relationship than the moving window method.
To obtain a better estimate of the coupling direction, theresulting causality markers 关Figs. 1共e兲, 1共f兲, 3共e兲, 3共f兲, 6共e兲,and 6共f兲兴 were partitioned into two classes by applying thefollowing classification rule 共thresholded with the averagevalue of wt兲:wt ⇒再w̃t = 1 if wt ⱖ 1 共causally related兲,w̃t = 0 if wt ⬍ 1 共causally unrelated兲.

冎

共9兲

026216-6