Causal Modelling Combining Instantaneous and Lagged Effects:
an Identifiable Model Based on Non-Gaussianity

Aapo Hyvärinen
Dept of Computer Science and HIIT, University of Helsinki, Finlandaapo.hyvarinen@helsinki.fi

Shohei Shimizu
The Institute of Scientific and Industrial Research, Osaka University

sshimizu@ar.sanken.osaka-u.ac.jp

Patrik O. Hoyer
Dept of Computer Science and HIIT, University of Helsinki, Finland

patrik.hoyer@helsinki.fi

Abstract
Causal analysis of continuous-valued variables typically uses either autoregressivemodels or linear Gaussian Bayesian networks
with instantaneous effects. Estimation ofGaussian Bayesian networks poses serious
identifiability problems, which is why it wasrecently proposed to use non-Gaussian models. Here, we show how to combine the nonGaussian instantaneous model with autoregressive models. We show that such a nonGaussian model is identifiable without priorknowledge of network structure, and we propose an estimation method shown to be consistent. This approach also points out how
neglecting instantaneous effects can lead tocompletely wrong estimates of the autoregressive coefficients.

1. Introduction
Analysis of causal influences or effects has become
an important topic in machine learning (Pearl, 2000;
Spirtes et al., 1993), and has numerous applicationsin, for example, neuroinformatics (Roebroeck et al.,2005; Kim et al., 2007) and bioinformatics (OpgenRhein & Strimmer, 2007). For continuous-valued variables, such an analysis can basically be performed in
two different ways. First, if the time-resolution of themeasurements is higher than the time-scale of causal
influences, one can estimate a classic autoregressivemodel with time-lagged variables and interpret the authAppearing in Proceedings of the 25 International Conference on Machine Learning, Helsinki, Finland, 2008. Copyright 2008 by the author(s)/owner(s).

toregressive coefficients as causal effects. Second, if themeasurements have a lower time resolution than thecausal influences, or if the data has no temporal structure at all, one can use a model in which the causal
influences are instantaneous, leading to Bayesian networks or structural equation models (Bollen, 1989).
While estimation of autoregressive methods can besolved by classic regression methods, the case of instantaneous effects is much more difficult. Most methods suffer from lack of identifiability,1 because covariance information alone is not sufficient to uniquely
characterize the model parameters. Prior knowledge ofthe structure (fixing some of the connections to zero) of
the Bayesian network is then necessary for most practical applications. However, a method was recently
proposed which uses the non-Gaussian structure of thedata to overcome the identifiability problem (Shimizu
et al., 2006): If the disturbance variables (external influences) are non-Gaussian, no prior knowledge on the
network structure (other than the ubiquitous assumption of a directed acyclic graph (DAG)) is needed toestimate the model.
Here, we consider the general case where causal influences can occur either instantaneously or with considerable time lags. Such a model is called the structuralvector autoregressive (SVAR) model in econometrictheory, in which numerous attempts have been madefor its estimation, see e.g. (Swanson & Granger, 1997;
Demiralp & Hoover, 2003; Moneta & Spirtes, 2006).
We propose to use non-Gaussianity to estimate themodel. We show that this variant of the model is iden1

Identifiability is here used in the classic statisticalsense: a model is identifiable if no two different values ofthe parameter vector give the same distribution for theobserved data.

Causal Modelling Combining Instantaneous and Lagged Effects

tifiable without any other restrictions than acyclicity.
To our knowledge, no model proposed for this problem
has been shown to be fully identifiable without priorknowledge of network structure. We further propose acomputational method for estimating the model basedon the theory of independent component analysis or
ICA (Hyvärinen et al., 2001).

First, they are mutually independent, and temporally
uncorrelated, which are typical assumptions in autoregressive models. Second, they are assumed to be nonGaussian, which is an important assumption which
distinguishes our model from classic models, whether
autoregressive models, structural-equation models, orBayesian networks.

The proposed non-Gaussian model not only allows estimation of both instantaneous and lagged effects; italso shows that taking instantaneous influences into
account can change the values of the time-lagged coefficients quite drastically. Thus, we see that neglecting
instantaneous influences can lead to misleading interpretations of causal effects. The framework further
leads to a generalization of the well-known Granger
causality measure.

Further, we assume that the matrix modelling instantaneous effects, B0 , corresponds to an acyclic graph,
as is typical in causal analysis, but this may not bestrictly necessary as will be discussed below. The
acyclicity is equivalent to the existence of a permutation matrix P, which corresponds to an ordering of
the variables xi , such that the matrix PB0 PT is lowertriangular (i.e. entries above the diagonal are zero).
Acyclicity also implies that the entries on the diagonal are zero, even before such a permutation.

The paper is structured as follows. We first define themodel and discuss its relation to other models in Section 2. In Section 3 we propose an estimation method,show its consistency, and discuss an intuitive interpretation of the method. Section 4 contains some theoretical examples and a theorem on how including instantaneous effects in the model changes the resulting interpretations. The resulting generalization of Granger
causality is discussed in Section 5. The validity of theestimation method is demonstrated by simulations on
artificial data in Section 6, and experiments on financial and neuroscientific data in Section 7. Section 8
concludes the paper.

2. Model Combining Lagged andInstantaneous Effects
2.1. Definition and Assumptions
Let us denote the observed time series by xi (t), i =1, . . . , n, t = 1, . . . , T where i is the index of the variables (time series) and t is the time index. All the variables are collected into a single vector x(t). Denote by
k the number of time-delays used, i.e. the order of theautoregressive model. Denote by Bτ the n × n matrixof the causal effects between the variables xi with timelag τ, τ = 0 . . . k .
The causal dynamics in our model are a combinationof autoregressive and structural-equation models. The
model is defined asx(t) =kX

Bτ x(t − τ ) + e(t)

(1)

τ =0

where the ei (t) are random processes modelling theexternal influences or “disturbances”. We make thefollowing assumptions on the external influences ei (t).

2.2. Relation to Other Models
This model is a generalization of the linear nonGaussian acyclic model (LiNGAM) proposed in(Shimizu et al., 2006). If the order of the autoregressive part is zero, i.e. k = 0, the model is nothing else than the LiNGAM model, modelling instantaneous effects only. As shown in (Shimizu et al., 2006),
the assumption of non-Gaussianity of the ei enables
estimation of the model. This is because the nonGaussian structure of the data provides information
not contained in the covariance matrix which is the
only source of information in most methods. In thissense the model is similar to independent component
analysis, which solves the unidentifiability of factor analytic models using the assumption of non-Gaussianityof the factors (Comon, 1994; Hyvärinen et al., 2001).
In fact, the estimation method in (Shimizu et al., 2006)
uses an ICA algorithm as an essential part.
On the other hand, if the matrix B0 has all zero entries, the model in Equation (1) is a classic vector
autoregressive model in which future observations arelinearly predicted from preceding ones. If we knew inadvance that B0 is zero, the model could thus be estimated by classic regression techniques since we do
not have the same variables on the left and right-handsides of Equation (1).
We emphasize that our model is different from classic
autoregressive models two important ways: First, theexternal influences ei (t) are non-Gaussian. Second, thelag variable τ takes the value 0 as well, which brings
instantaneous effects into the model in the form of thematrix B0 . A coefficient B0 (i, j) models the instantaneous effect of xj (t) on xi (t) as in a linear Bayesian
network, or a structural equation model.

Causal Modelling Combining Instantaneous and Lagged Effects

2.3. Causality vs. Prediction
An autoregressive model can serve two different goals:
prediction and analysis of causality. Our goal here is
the latter: We estimate the parameter matrices Bτ inorder to interpret them as causal effects between thevariables. This goal is distinct from simply predicting
future outcomes when passively observing the time series, as has been extensively discussed in the literature
on causality (Pearl, 2000; Spirtes et al., 1993). Thus,
we emphasize that our model is not intended to reduce
prediction errors if we want to predict xi (t) using (passively) observed values of the past x(t−1), x(t−2), . . .;for such prediction, an ordinary autoregressive model
is likely to be just as good.
Our model is intended to be superior in causal modelling. Causality has an obvious intuitive interpretation, which is typically formalized as the ability to predict the effect of possible new interventions on the system (Pearl, 2000). Thus, our model should be betterin predicting effects of interventions, which is different
from conventional time series prediction.

4. Finally, compute the estimates of the causal effectmatrices Bτ for τ > 0 as
B̂τ = (I − B̂0 )M̂τ for τ > 0(5)This estimation method is consistent,2 as will beshown in Section 3.3. First, however, we show thederivation of Equation (5) and discuss its deep meaning.
3.2. Why Autoregressive Matrices Change due
to Instantaneous Influences
Equation (5) shows a remarkable fact already mentioned in the Introduction: Consistent estimates of theBτ are not obtained by a simple AR model fit evenfor τ > 0. Taking instantaneous effects into accountchanges the estimation procedure for all the autoregressive matrices, if we want consistent estimators as
we usually do. Of course, this is only the case if there
are instantaneous effects, i.e. B0 6= 0; otherwise, theestimates are not changed.
Why do we have (5)? This is because from (1) we have

3. Estimation of the Model(I − B0 )x(t) =3.1. Combining Least-Squares Estimation and
LiNGAM
We propose the following method for estimating ourmodel defined in Section 2.1. The method combinesclassic least-squares estimation of an autoregressive(AR) model with LiNGAM estimation:
1. Estimate a classic autoregressive model for thedatakXMτ x(t − τ ) + n(t)(2)x(t) =
τ =1

using any conventional implementation of a leastsquares method. Note that here τ > 0, so it is really a classic AR model. Denote the least-squares
estimates of the autoregressive matrices by M̂τ .
2. Compute the residuals, i.e. estimates of innovations n(t)n̂(t) = x(t) −

k
X

k
X

Bτ x(t − τ ) + e(t)

(6)

τ =1

and thusx(t) =kX(I − B0 )−1 Bτ x(t − τ ) + (I − B0 )−1 e(t) (7)τ =1Comparing this with (2), we can equate the autoregressive matrices, which gives (I − B0 )−1 Bτ = Mτ forτ ≥ 1, and thus (5) is justified.
While this phenomenon is, in principle, well-known
in econometric literature (Swanson & Granger, 1997;
Demiralp & Hoover, 2003; Moneta & Spirtes, 2006),
Equation (5) is seldom applied because estimation
methods for B0 have not been well developed. To
our knowledge, no estimation method for B0 has been
proposed which is consistent without strong prior assumptions on B0 .
3.3. Consistency and IdentifiabilityM̂τ x(t − τ )

(3)

τ =1

3. Perform the LiNGAM analysis (Shimizu et al.,2006) on the residuals. This gives the estimate of
the matrix B0 as the solution of the instantaneous
causal modeln̂(t) = B0 n̂(t) + ẽ(t)(4)

The consistency of our method relies on two facts.
First, in the estimation of an AR model as in (2), it
is not necessary that the innovation vector n(t) has
independent or even uncorrelated elements (for fixed2
Consistency means classic statistical consistency, i.e.
the estimator converges in probability to the right parameter values when the data follows the model and sample
size grows infinite.

Causal Modelling Combining Instantaneous and Lagged Effectst); least-squares estimation will still be consistent, as
is well known. Thus, least-squares estimation of (2),combined with (5), gives consistent estimators of Bτ
for τ ≥ 1, provided we have a consistent estimatorof B0 . Second, comparison of (7) with (2) showsthat the residuals n̂(t) are, asymptotically, of the form
(I − B0 )−1 e(t). This meansn̂(t) = (I − B0 )−1 e(t) ⇔ (I − B0 )n̂(t) = e(t)⇔ n̂(t) = B0 n̂(t) + e(t)(8)

which is the LiNGAM model for n̂(t). This showsthat B0 is obtained as the LiNGAM analysis of theresiduals, and the consistency of our estimator of B0
follows from the consistency of LiNGAM estimation(Shimizu et al., 2006). Thus, our method is consistent
for all the Bτ . This obviously proves, by construction,
the identifiability of the model as well.
We have here assumed that B0 is acyclic, as is typicalin causal analysis. However, this assumption is only
made because we do not know very well how to estimate a linear non-Gaussian Bayesian network in thecyclic case. Future work may produce methods which
estimate cyclic models, and then we do not need theassumption of acyclicity in our combined model either.
We could just use such a new method in Step 3 of themethod instead of LiNGAM, and nothing else would
be changed. Recent work in that direction is in (Lacerda et al., 2008); see also (Richardson & Spirtes, 1999)
for older methods on Gaussian data.
3.4. Interpretation as ICA of Residuals
Another viewpoint on our model is analysis of the correlations of the innovations after estimating a classic
AR model. Suppose we just estimate an AR model asin (2), and interpret the coefficients as causal effects.
Such an interpretation more or less presupposes thatthe innovations ni are independent of each other, because otherwise there is some structure in the model
which has not been modelled by the AR model. If theinnovations are not independent, the causal interpretation may not be justified. Thus, it seems necessary
to further analyze the dependencies in the innovations
in cases where they are strongly dependent.
Analysis of the dependency structure in the residuals (which are, by definition, estimates of innovations)is precisely what leads to the present model. As afirst approach, one could consider application of something like principal component analysis or independent
component analysis on the residuals. The problem
with such an approach is that the interpretation of
the obtained results in the framework of causal analysis would be quite difficult. Our solution is to fit

a causal model like LiNGAM to the residuals, whichleads to a straightforward causal interpretation of theanalysis of residuals which is logically consistent with
the AR model.

4. Interaction of Instantaneous and
Lagged Effects
Here we present some theoretical examples of how the
instantaneous and lagged effects interact based on theformula in (5).
An instantaneous effect may seem to be lagged
Consider first the case where the instantaneous andlagged matrices are as follows:




0 1
0.9 0
B0 =
,
B1 =
(9)
0 0
0 0.9
That is, there is an instantaneous effect x2 → x1 , andno lagged effects (other than the purely autoregressive xi (t − 1) → xi (t)). Now, if an AR(1) model isestimated for data coming from this model, without
taking the instantaneous effects into account, we get
the autoregressive matrix
0.9 0.9(10)
M1 = (I − B0 )−1 B1 =0 0.9
Thus, the effect x2 → x1 seems to be lagged although
it is, actually, instantaneous.
Spurious effects appear Consider three variables
with the instantaneous effects x1 → x2 and x2 → x3 ,and no lagged effects other than xi (t − 1) → xi (t), as
given by




0 0 0
0.9 0
0
B0 = 1 0 0 , B1 =  0 0.9 0  (11)
0 1 0
0
0 0.9
If we estimate an AR(1) model for the data comingfrom this model, we obtain


0.9 0
0
M1 = (I − B0 )−1 B1 = 0.9 0.9 0 (12)
0.9 0.9 0.9This means that the estimation of the simple autoregressive model leads to the inference of a direct lagged
effect x1 → x3 , although no such direct effect exists inthe model generating the data, for any time lag.
Causal ordering is not changed A more reassuring result is the following: if the data follows the samecausal ordering for all time lags, that ordering is notcontradicted by the neglect of instantaneous effect. A
rigorous definition of this property is the following.

Causal Modelling Combining Instantaneous and Lagged Effects

Theorem 1 Assume that there is an ordering
i(j), j = 1 . . . n of the variables such that no effect goes
backward,3 i.e.
Bτ (i(j −δ), i(j)) = 0 for δ > 0, τ ≥ 0, 1 ≤ j ≤ n (13)Then, the same property applies to the Mτ , τ ≥ 1 as
well. Conversely, if there is an ordering such that (13)applies to Mτ , τ ≥ 1 and B0 , then it applies to Bτ , τ ≥1 as well.
The proof of the theorem is based on the fact thatwhen the variables are ordered in this way (assuming such an order exists), all the matrices Bτ
are lower-triangular. The same applies to I − B0 .
Now, the product of two lower-triangular matrices islower-triangular; in particular the Mτ are also lowertriangular according to (5), which proves the first partof the theorem. The converse part follows from solving for Bτ in (5) and the fact that the inverse of alower-triangular matrix is lower-triangular.
What this theorem means is that if the variables reallyfollow a single “causal ordering” for all time lags, thatordering is preserved even if instantaneous effects are
neglected and a classic AR model is estimated for thedata. Thus, there is some limit to how (5) can change
the causal interpretation of the results.

5. Towards a Generalization of Granger
Causality
The classic interpretation of causality in instantaneous
Bayesian network models would be that xi causes xj
if the (j, i)-th coefficient in B0 is non-zero. In thetime series context, this is related to Granger causality (Granger, 1969), which formalizes causality as theability to reduce prediction error. A simple operational definition of Granger causality can be based on
the autoregressive coefficients Mτ : If at least one of
the coefficients from xi (t − τ ), τ ≥ 1 to xj (t) is (significantly) non-zero, then xi Granger-causes xj . Thisis because then the variable xi reduces the prediction
error in xj in the mean-square sense if it is includedin the set of predictors, which is the very definition ofGranger causality.
In light of the results in this paper, we propose adefinition which combines the two aspects: A variable xi causes xj if at least one of the coefficientsBτ (j, i), giving the effect from xi (t − τ ) to xj (t), is(significantly) non-zero for τ ≥ 0. The condition forτ is different from Granger causality since the value
3

In the purely instantaneous case, existence of such anordering is equivalent to acyclicity of the effects as notedin Section 2.1.

τ = 0, corresponding to instantaneous effects, is included. Moreover, since estimation of the instantaneous effects changes the estimates of the lagged ones,
the lagged effects used in our definition are different
from those usually used with Granger causality.
A more general formulation of this definition, which isin line with the general formulation of Granger causality, is that the error in the “prediction” of xj (t) is
reduced when xi (t − 1), xi (t − 2), . . . and xi (t) are included in the set of predictors. Here, we use a rather
unconventional definition of the word “prediction” because we include instantaneous effects.

6. Simulations
To verify the validity of our method, we first performedexperiments with artificial data. In the experiments,
we created data in the following manner using the
LiNGAM code package4 :1. We randomly constructed a strictly lowertriangular matrix (i.e. zero entries above and on
the diagonal), B0 , for the instantaneous causal
model so that the standard deviations of the innovations ni owing to parent innovations will bein the interval [0.5, 1.5]. The number of observed
time-series was n = 10. Both fully connected(no zeros in the strictly lower triangular part)and sparse networks (many zeros) were created.
We also randomly selected the standard deviations of the external influences ei from the interval[0.5, 1.5].
2. Next, we generated data with various lengths ofthe time series (300, 500 and 1,000) by independently drawing the external influences ei from various non-Gaussian distributions with zero mean
and unit variance5 . The values of the innovations
ni were generated according to the assumed instantaneous recursive process. This is straightforward because B0 is lower-triangular, so we justgenerate the ni in the order n1 , n2 . . . as is typicalin acyclic networks, e.g. (Shimizu et al., 2006).
3. We randomly permuted the order of the innovations ni to hide the causal order with which thedata was generated. We also permuted B0 as well4

http://www.cs.helsinki.fi/group/neuroinf/lingam/
We first generated a gaussian variable z with zero mean
and unit variance and subsequently transformed it to a
non-Gaussian variable by ei = sign(z)|z|q . The nonlinear
exponent q was selected to lie in [0.5, 0.8] or [1.2, 2.0]. The
former gave a sub-gaussian variable, and the latter a supergaussian variable. Finally, the transformed variable wasstandardized to have zero mean and unit variance.
5

Causal Modelling Combining Instantaneous and Lagged Effects

For the scatterplots in the left and center columns,
the estimation worked well when the sample size grew,as evidenced by the grouping of the data points onto
the main diagonal, although for the small sample size
300 the estimation was often inaccurate. On the other
hand, the scatterplots in the right column confirmedthat the causal effects were not correctly estimated by
the ordinary autoregressive coefficients when instantaneous influences existed since the data points were notvery close to the main diagonal.

7. Experiments on Real Data
7.1. Financial Data

Figure 1. Simulations on artificial data. Left column: Scatterplots of the estimated elements of B0 versus the generating values. Center column: Scatterplots of the estimatedelements of B1 versus the generating values. Right column:
Scatterplots of the estimated elements of M1 versus thoseof B1 .The number of observed signals was 10. Five data
sets were generated for each scatterplot.

as the variances of the external influences ei tomatch the new order.
4. We randomly generated a first-order autoregressive matrix M1 so that the spectral norm of thematrix was less than 0.99 to ensure the stabilityof the autoregressive process.
5. The values of the observed signals xi (t) were generated according to the assumed first-order autoregressive process.
6. Finally, we fed the data to our estimation method.
Here we told the method that the generating autoregressive order was 1.
Figure 1 gives the scatterplots of the elements of theestimated parameters versus the generating ones. The
left column is for the scatterplots of the estimatedcausal effects in B0 versus the generating values. The
center column is for the scatterplots of the estimatedcausal effects in B1 versus the generating values. The
right column is for the scatterplots of the estimatedautoregressive coefficients in M1 versus the generating values of the causal effects in B1 (here, the estimation was invalid because instantaneous effects wereignored).

As a first illustration of the applicability of the method
on real data, we analyzed a dataset from a time series repository on the Internet.6 The data consistedof two observed signals, x1 : weekly closing price ofToyota stock and x2 : weekly closing rate of exchange
of Japanese Yen to U.S. Dollar in 2007. The numberof time points was 50. The maximum, minimum andmean of x1 were 8,230, 5,870 and 7,102 (JPY). Thoseof x2 were 123.86, 108.51 and 117.72 (JPY).
We analyzed the data using our method with autoregressive order of 1. The estimated first-order autoregressive matrix M1 and residual correlation matrix
were as follows:
0.95 −4.22
(14)
M1 =
0.00080.78

1.00 0.66
corr(n) =0.66 1.00
The relatively strong correlation between the residuals implied that there would be some dependency that
had not been modeled by the AR model. Thus, we fitted the instantaneous causal model to the residuals, as
proposed above. The estimated instantaneous causal
effect matrix B0 and resulting lagged causal effect matrix B1 were as follows:


0
56.04
B0 =
(15)
0.00270
0.91 −48.01
B1 =
(16)
−0.0018
0.79
The matrix B0 is very close to be upper-triangular,
which implied that the model was really acyclic (because switching the order of the variables would makeB0 lower-triangular). Further, the instantaneous effect x2 →x1 in B0 was one order of magnitude larger
than the lagged effect in M1 and thus the lagged co6

Yahoo! Japan Finance: http://quote.yahoo.co.jp/

Causal Modelling Combining Instantaneous and Lagged Effects

efficients in M1 are quite different from those in B1 ,due to the formula in (5).
Figure 2 shows a graphical representation of the estimated model for financial data. First, it implies
that a higher value of the yen (x2 ) had a negative
lagged effect (-48.01) on the price of Toyota stock (x1 ).
This would be reasonable since Toyota sells many carsabroad, and a higher value of the yen would increase
the cost price and decrease the earning. Interestingly,
it was also implied that a higher value of the yen had
a positive instantaneous effect (56.04) on the price ofToyota stock. In other words, for weeks where valuesof the yen one week before were the (approximately)same, if the yen got more expensive (due to some reason other than the value of the yen one week before,
perhaps a U.S. recession, for example) then the priceof Toyota stock would get more expensive. It would beinteresting to further study the economic mechanismwith more extensive data.

was sitting with eyes closed, and did not perform any
specific task nor was there any specific sensory stimulation. The channels were first linearly projected tothe signal space to reduce noise (Uusitalo & Ilmoniemi,
1997). In this illustrative experiment, we only consider
a single (gradiometer) channel in the right occipitalcortex near the midline.
We considered the interaction of about 10 Hz (alpha)and about 20 Hz (beta) oscillations commonly observed in electromagnetic recordings of spontaneous
brain activity. We first computed the amplitudes ofthe oscillations by dividing the data into windows oflength of 0.25 seconds, performing fast Fourier transform inside each of them, and computing the total
Fourier amplitudes (unweighted Euclidean norm of
the Fourier coefficients) in the frequency ranges of
8 . . . 12Hz (alpha range, denoted by x1 ) and 15 . . . 25Hz
(beta range, denoted by x2 ). Thus we obtained twotime series of 1,200 points.
We fitted our model, with autoregressive order of 1 tothe data. The obtained matrices are
0.23381 0.14551
M1 =
(17)
0.10838 0.14314

0−0.65768(18)
B0 =0.5672200.305090.23965B1 =(19)−0.024244 0.060608

Figure 2. A graphical representation of the model estimated in Section 7.1. The x1 and x2 denote weekly closing
price of Toyota stock in 2007 and weekly closing rate ofexchange of Japanese Yen to U.S. Dollar in 2007, respectively. The arrow from x1 (t − 1) to x2 (t) was omitted since
the estimated strength was very close to zero (-0.0018).

7.2. Magnetoencephalographic Data
As a second illustration of the applicability of themethod on real data, we applied it on magnetoencephalography (MEG), i.e. measurements of the electric activity in the brain. The raw data consisted of
the 306 MEG channels measured by the Vectorview
helmet-shaped neuromagnetometer (Neuromag Ltd.,
Helsinki, Finland) in a magnetically shielded room at
the Brain Research Unit, Low Temperature Laboratory, Helsinki University of Technology. The sampling
frequency was 600 Hz. The measurements consistedof 300 seconds of resting state brain activity from theexperiment of (Ramkumar et al., 2007). The subject

What we see is that the instantaneous model is far
from trivial: the effects in B0 are relatively strong.
This is also reflected in B1 which is now rather different from M1 . Thus, the interpretation of the autoregressive matrices using just the autoregressive model(i.e. M1 ) or the combined model (i.e. B1 ) are quite
different. In the classic autoregressive case (based on
M1 ), the lagged effect x1 → x2 is relatively strongly
positive whereas in the combined model it is quite
weak. In fact, that effect is now modelled as an instantaneous effect in B0 . Even more interesting is thatthe instantaneous model has a strong negative effect
x2 → x1 which is not visible at all in the purely autoregressive matrix M1 . Thus, the results illustrate
how the interpretation of causal effects (and even of
the lagged ones) can change drastically when including the instantaneous effects.
Using an autoregressive order of 2 did not change theresults. We also ran the method many times to exclude
the problem of the ICA estimation algorithm (used in
LiNGAM estimation) getting stuck in local minima(Himberg et al., 2004), and the result was found to berobust with respect to that manipulation.

Causal Modelling Combining Instantaneous and Lagged Effects

One problem with this experiment is that the causal
model estimated by LiNGAM is far from acyclic.
Here, we can justify the procedure by using the theory of cyclic model estimation proposed by (Lacerda
et al., 2008); the estimation here gives the only“stable”
model according to that theory. Performance of
LiNGAM estimation methods in the case of cyclic
models, and the possible need for new methods for estimating cyclic models are future research topics of greatpractical importance. However, as discussed above,
they are separate from the main contribution of our paper in the sense that we can use any such new method
to estimate the instantaneous model in our framework.

8. Conclusion
We showed how non-Gaussianity enables estimationof a causal discovery model in which the linear effects
can be either instantaneous or time-lagged. Like inthe purely instantaneous case (Shimizu et al., 2006),
non-Gaussianity makes the model identifiable without explicit prior assumptions on existence or nonexistence of given causal effects. The classic assumption of acyclicity is sufficient although probably notnecessary. From the practical viewpoint, an important implication is that considering instantaneous effects changes the coefficient of the time-lagged effects
as well.

References
Bollen, K. A. (1989). Structural equations with latentvariables. John Wiley & Sons.
Comon, P. (1994). Independent component analysis—a new concept? Signal Processing, 36, 287–314.
Demiralp, S., & Hoover, K. D. (2003). Searching forthe causal structure of a vector autoregression. Oxford Bulletin of Economics and Statistics, 65 (supplement), 745–767.
Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral
methods. Econometrica, 37, 424–438.
Himberg, J., Hyvärinen, A., & Esposito, F. (2004).
Validating the independent components of neuroimaging time-series via clustering and visualization. NeuroImage, 22, 1214–1222.
Hyvärinen, A., Karhunen, J., & Oja, E. (2001). Independent component analysis. Wiley Interscience.
Kim, J., Zhu, W., Chang, L., Bentler, P. M., & Ernst,
T. (2007). Unified structural equation modeling approach for the analysis of multisubject, multivariate

functional MRI data. Human Brain Mapping, 28,
85–93.
Lacerda, G., Spirtes, P., Ramsey, J., & Hoyer, P. O.
(2008). Discovering cyclic causal models by independent components analysis. Proc. 24th Conf. on
Uncertainty in Artificial Intelligence (UAI2008).
Helsinki, Finland.
Moneta, A., & Spirtes, P. (2006). Graphical models forthe identication of causal structures in multivariate
time series models. Proc. Joint Conference on Information Sciences. Kaohsiung, Taiwan.
Opgen-Rhein, R., & Strimmer, K. (2007). From correlation to causation networks: a simple approximate learning algorithm and its application to highdimensional plant gene expression data. BMC Systems Biology, 1.
Pearl, J. (2000). Causality: Models, reasoning, andinference. Cambridge University Press.
Ramkumar, P., Parkkonen, L. T., He, B. J., Raichle,
M. E., Hämäläinen, M. S., & Hari, R. (2007). Identification of stimulus-related and intrinsic networks
by spatial independent component analysis of MEG
signals. Abstract presented at the Society for Neuroscience Meeting, San Diego, California.
Richardson, T. S., & Spirtes, P. (1999). Automateddiscovery of linear feedback models. In C. Glymour
and G. Cooper (Eds.), Computation, causation anddiscovery, 253–302. The MIT Press.
Roebroeck, A., Formisano, E., & Goebel, R. (2005).
Mapping directed influence over the brain using
granger causality and fMRI. NeuroImage, 25, 230–
242.
Shimizu, S., Hoyer, P. O., Hyvärinen, A., & Kerminen,
A. (2006). A linear non-Gaussian acyclic model forcausal discovery. J. of Machine Learning Research,
7, 2003–2030.
Spirtes, P., Glymour, C., & Scheines, R. (1993). Causation, prediction, and search. Springer-Verlag.
Swanson, N. R., & Granger, C. W. J. (1997). Impulse
response functions based on a causal approach toresidual orthogonalization in vector autoregression.
J. of the Americal Statistical Association, 92, 357–
367.
Uusitalo, M. A., & Ilmoniemi, R. J. (1997). Signalspace projection method. Med. Biol. Eng., 32, 35–42.