## Abstract

We present an approach for the quantification of directional relations in multiple time series exhibiting significant zero-lag interactions. To overcome the limitations of the traditional multivariate autoregressive (MVAR) modelling of multiple series, we introduce an extended MVAR (eMVAR) framework allowing either exclusive consideration of time-lagged effects according to the classic notion of Granger causality, or consideration of combined instantaneous and lagged effects according to an extended causality definition. The spectral representation of the eMVAR model is exploited to derive novel frequency domain causality measures that generalize to the case of instantaneous effects the known directed coherence (DC) and partial DC measures. The new measures are illustrated in theoretical examples showing that they reduce to the known measures in the absence of instantaneous causality, and describe peculiar aspects of directional interaction among multiple series when instantaneous causality is non-negligible. Then, the issue of estimating eMVAR models from time-series data is faced, proposing two approaches for model identification and discussing problems related to the underlying model assumptions. Finally, applications of the framework on cardiovascular variability series and multichannel EEG recordings are presented, showing how it allows one to highlight patterns of frequency domain causality consistent with well-interpretable physiological interaction mechanisms.

## 1. Introduction

The assessment of causality, i.e. the presence of directional interactions among a set of measured variables, is becoming of paramount importance in the study of physiological systems. Application of methods aimed at detecting causality from the analysis of experimental multivariate (MV) time series, indeed, ranges from neurophysiology [1] to cardiovascular physiology [2]. While a universally accepted definition of causality is still lacking, a very popular and practically useful notion is that first proposed by Wiener [3] and then formalized by Granger [4] in the context of linear regression modelling of MV stochastic processes. Wiener–Granger causality is based on the principles that the cause occurs in time before the effect, and that the causing process contains information about the caused process that is unique, i.e. not present in any other relevant process. While its original definition was set in the time domain [4], the interest in applying this concept to physiological processes that exhibit important and meaningful oscillatory content has fostered the development of methods operating in the frequency domain. A pioneering approach was that proposed by Akaike [5], who introduced a causality measure based on spectral decomposition, the so-called noise contribution ratio, which was named directed coherence (DC) by Saito & Harashima [6] for bivariate time series, and was used in the context of MV time series later on [7]. As a further step, partial directed coherence (PDC) was introduced by Baccala & Sameshima [8,9] as a straightforward frequency domain descriptor of Granger causality. DC and PDC are derived as factors in the decomposition of ordinary coherence and partial coherence and, as such, are able to elicit information on directionality from the spectral representation of multiple time series [10]. While, in bivariate systems, DC and PDC provide essentially the same information, their behaviour differs when MV series are considered: PDC exclusively detects direct interactions between two series in the MV representation; DC captures both direct and indirect effects (i.e. effects mediated by one or more other series) [10].

The computation of DC and PDC is based on fitting the observed set of time series with a parametric model, interpreting the model coefficients as causal effects, and finally elaborating the coefficients in the spectral domain to get frequency-dependent information on causality. Such an interpretation presupposes that the developed model is able to capture the whole interaction structure of the observed time series; if this is not the case, the causal interpretation provided by the model coefficients may not be justified and, ultimately, then DC/PDC may yield misleading frequency domain patterns of causality. The model underlying the definition of DC and PDC is a strictly causal MV autoregressive (MVAR) model that describes, for each time series, the linear causal effects coming from its own past and from the past of all other time series. Strict causality of the model means that only lagged effects are modelled, whereas instantaneous (i.e. not lagged) effects are not described by any model coefficient. Therefore, a strictly causal MVAR model can fully describe a set of multiple time series, and thus yield unambiguous causal information based on its coefficients in either time or frequency domains, only when instantaneous effects are negligible. Despite this basic requirement, strictly causal MVAR models are ubiquitously used in the experimental frequency domain causality analysis without verifying the absence of instantaneous effects among the modelled time series [11–16]. It might be questioned that instantaneous causality is not a practical issue as zero-lag causal effects are unattainable in real-world physical systems where interactions take time to occur. Actually, instantaneous causality shows up in practical time-series modelling whenever the time resolution of the measurements is lower than the time scale of the lagged causal influences occurring among the processes under analysis. In the study of experimental time series, instantaneous causality may reflect either fast (within-sample) physiologically meaningful interactions or non-physiological effects (e.g. due to unobserved confounders).

The adverse impact of instantaneous effects on MVAR-based causality analysis has been made explicit in recent studies [17,18], which have provided theoretical arguments showing that the omission of zero-lag correlations from the model can change quite drastically the values of the time-lagged coefficients, and thus of the frequency domain causality measures. As a possible solution, these studies proposed the utilization of an extended MVAR (eMVAR) model accounting for both instantaneous and time-lagged effects in such a way that the whole data covariance could be explained in terms of the model coefficients. In this study, we build on this idea, presenting a complete framework for the analysis of causality in the frequency domain based on eMVAR models. After providing a thorough categorization of the various forms of causality that may be defined for an MV process, we particularize them to eMVAR processes and derive frequency domain measures that generalize the known DC and PDC. The correspondence of these novel measures with the time domain causality definitions is illustrated in theoretical examples. Moreover, we elaborate on the practical estimation of eMVAR processes on time-series data. In fact, while strictly causal MVAR identification is easily performed by classic regression methods, the inclusion of instantaneous effects in the model makes eMVAR identification much less straightforward. We show that the problem can only be solved using more information than that provided by the data covariance alone, and propose two alternative identification procedures based on imposing the direction of the instantaneous transfer paths, or on assuming non-Gaussianity of the model innovations. Finally, we present an application of the framework to two datasets of physiological time series that entail utilization of the two identification procedures and investigation of different aspects of causal information transfer based on DC and PDC, i.e. cardiovascular variability series and multichannel electroencephalography (EEG) recordings.

We distribute a Matlab toolbox for the identification of eMVAR processes and estimation of related frequency domain causality measures; the toolbox is available online at http://www.science.unitn.it/biophysicslab/research/sigpro/eMVAR.html.

## 2. Causality in linear multivariate processes

### (a) Definitions of causality in the time domain

According to the framework first introduced by Granger [4], the concept of causality is defined in terms of prediction improvement, i.e. reduction of the variance of the prediction error resulting from modifying the information set on which the prediction model is conditioned. Let us consider the MV vector stochastic process **Y** composed of *M* scalar processes of zero mean, **Y**=[*y*_{1},…,*y*_{M}]^{T}. To introduce the notation with regard to the *m*th scalar process, *m*=1,…,*M*, we let *y*_{m}(*n*) denote the current value of *y*_{m}, *Y* _{m}={*y*_{m}(*n*−1),…,*y*_{m}(*n*−*p*)} the set of the *p* past values of *y*_{m}, and the set of the *p* past values enlarged with the current value. Moreover, *Z* is an information set consisting of an ensemble of values properly chosen from the process, and *σ*^{2}(*y*_{m}|*Z*) is the prediction error variance of the optimal linear predictor of *y*_{m}(*n*) based on *Z*. Given this notation, different definitions of causality from the process *y*_{j} to the process *y*_{i} (*i*,*j*=1,…,*M*;*i*≠*j*) are formally provided in table 1 and are interpreted as follows. *Instantaneous causality* from *y*_{j} to *y*_{i} exists if the knowledge of *y*_{j}(*n*) improves the prediction of *y*_{i}(*n*); *lagged direct causality* from *y*_{j} to *y*_{i}, , exists if the knowledge of *Y* _{j} improves the prediction of *y*_{i}(*n*); *extended direct causality* from *y*_{j} to *y*_{i}, , exists if the knowledge of improves the prediction of *y*_{i}(*n*); *lagged causality* from *y*_{j} to *y*_{i}, *y*_{j}⇒*y*_{i}, exists if a cascade of *L* direct causality relations occurs such that ; *extended* *causality* from *y*_{j} to *y*_{i}, , exists if a cascade of *L* extended direct causality relations occurs such that .

A first differentiation among the definitions provided above may be made on the basis of the role played by instantaneous effects, i.e. effects from one series to another occurring within the same time lag. These effects are the basis of instantaneous causality, and are excluded from the definitions of direct causality and causality, which consider only time-lagged effects. On the contrary, instantaneous effects are explicitly accounted for in the extended causality definitions. In the absence of instantaneous causality, extended direct causality reduces to lagged direct causality, and extended causality reduces to lagged causality. Another distinction among the proposed causality definitions involves the comparison between direct and indirect effects. In fact, lagged causality and extended causality may be viewed as generalizations of lagged direct causality and extended direct causality, respectively, which consider the indirect effects between two processes (i.e. effects mediated by one or more other processes), in addition to the direct effects. Note that for a bivariate process (*M*=2) lagged causality reduces to lagged direct causality and extended causality reduces to extended direct causality. Note also that the original definition of causality provided by Granger [4] corresponds to our definition of lagged direct causality stated for bivariate processes. For MV processes (*M*≥3), our lagged direct causality definition agrees with the notion of *prima facie* cause introduced later on [19], whereas our lagged causality definition may be viewed as a generalization incorporating indirect effects. As to the extended definitions, they may be viewed as further generalizations that combine instantaneous effects with the lagged effects traditionally studied for the estimation of causality, in accordance with recently proposed ideas [17,18].

### (b) Extended multivariate autoregressive processes

The causality definitions provided above may be expressed in a parametric form describing the considered process by means of an eMVAR model:
2.1where **Y**(*n*)=[*y*_{1}(*n*),…,*y*_{M}(*n*)]^{T} is the vector variable obtained by sampling the stochastic process **Y** at time instant *n*, **B**(*k*) are *M*×*M* coefficient matrices in which the element *b*_{ij}(*k*) describes the dependence of *y*_{i}(*n*) on *y*_{j}(*n*−*k*) (*i*,*j*=1,…,*M*;*k*=0,1,…,*p*), and **W**=[*w*_{1},…,*w*_{M}]^{T} is an innovation process formed by white and independent scalar processes with diagonal covariance matrix . Note that the formulation in (2.1) represents an extension of classic strictly causal modelling where present observations are linearly predicted from past ones [20]. The difference with respect to the strictly causal representation (see also §3*a*) consists of the fact that the extended model describes explicitly instantaneous effects from one scalar process to another in the form of the matrix **B**(0). As a consequence, the extended representation in (2.1) allows us to express all the definitions of causality provided above in terms of the off-diagonal elements of the coefficient matrices **B**(*k*). Specifically, instantaneous causality, lagged direct causality and extended direct causality occur from *y*_{j} to *y*_{i} when *b*_{ij}(0)≠0, *b*_{ij}(*k*)≠0 for at least one *k*≥1, and *b*_{ij}(*k*)≠0 for at least one *k*≥0, respectively. Correspondingly, causality and extended causality are detected when non-zero coefficients are present, for at least one relevant lag *k*, in the positions of **B**(*k*) that identify each direct connection of the cascade linking *y*_{j} to *y*_{i}.

### (c) Measures of causality in the frequency domain

The spectral representation of an eMVAR process is obtained taking the Fourier transform (FT) of (2.1) to yield the equation **Y**(*f*)=**B**(*f*)**Y**(*f*)+**W**(*f*), where **Y**(*f*) and **W**(*f*) are the FTs of **Y**(*n*) and **W**(*n*), and the *M*×*M* frequency domain coefficient matrix **B**(*f*) results as . If we want to consider the transfer function from **W**(*n*) to **Y**(*n*), then the spectral representation may be rewritten as **Y**(*f*)=**G**(*f*)**W**(*f*), where is the *M*×*M* transfer matrix in the frequency domain. The elements of the coefficient and transfer matrices may be suitably combined to define the extended directed coherence (eDC)
2.2and the extended partial directed coherence (ePDC)
2.3The functions defined in (2.2) and (2.3) constitute directional frequency domain measures of connectivity between the processes *y*_{i} and *y*_{j}, because *ξ*_{ij}(*f*) and *χ*_{ij}(*f*) quantify the influence of *y*_{j} on *y*_{i} at frequency *f*, as opposed to *ξ*_{ji}(*f*) and *χ*_{ji}(*f*), which quantify the influence in the opposite direction from *y*_{i} to *y*_{j}. eDC and ePDC constitute an extension of the well-known DC and PDC functions [6,7,8,9], the extension being in the fact that they are derived from the extended model (2.1), which describes instantaneous effects in addition to the commonly studied lagged effects.

When we use the extended measures (2.2) and (2.3), the information that flows from one process to another is both lagged (*k*>0) and instantaneous (*k*=0), because it is measured in the frequency domain by the function , which incorporates both **B**(0) and **B**(*k*) with *k*>0, and by its inverse **G**(*f*). If we want to explore lagged causality in the presence of zero-lag interactions, we have to exclude the coefficients related to the instantaneous effects from the desired spectral causality measure. Hence, we set
2.4and then we define the normalized lagged directed coherence (nDC)
2.5and the normalized lagged partial directed coherence (nPDC)
2.6nDC and nPDC constitute a variant of the known DC and PDC functions [6,7,8,9], the variation being in the fact that they are derived from the MVAR model (2.1), including instantaneous effects, rather than from a classic strictly causal MVAR model. It can be shown that, because, in the absence of instantaneous causality, the eMVAR model reduces to a strictly causal model (see §3*a*), in this case, eDC and nDC amount to the DC, whereas ePDC and nPDC amount to the PDC.

### (d) Interpretation of frequency domain causality measures

A straightforward interpretation of the four causality measures defined in §2*c* may be obtained considering that they reflect in the frequency domain the different time domain definitions of causality given in §2*a*. First, we note that the numerator terms of (2.3) and (2.6) are non-zero, with *i*≠*j*, when *b*_{ij}(*k*)≠0 for at least one *k*≥0 and at least one *k*≥1, respectively. Therefore, ePDC and nPDC measure in the frequency domain the concepts of extended direct causality and lagged direct causality, respectively. As to the eDC and nDC, one can show that, expanding the numerator terms of (2.2) or (2.5) as a geometric series, a sum of terms is obtained in which each term contains the FT of the model coefficients describing one of the (direct or indirect) pathways that connect *y*_{j} to *y*_{i} [21]. Therefore, eDC and nDC measure in the frequency domain the concepts of extended causality and lagged causality, respectively. The correspondence between time domain definitions and frequency domain measures is summarized in table 1.

Because the causality measures defined in §2*c* are complex-valued, the squared modulus of nPDC, ePDC, nDC and eDC is commonly used to quantify directional interactions in the frequency domain. Thus, , |*χ*_{ij}(*f*)|^{2}, and |*ξ*_{ij}(*f*)|^{2} are computed to quantify, respectively, lagged direct causality, extended direct causality, lagged causality and extended causality as a function of frequency. All these squared measures are normalized, so that they take values between 0, representing absence of causal coupling, and 1, representing full causal coupling from the process *y*_{j} to the process *y*_{i} at the frequency *f*. The concept of ‘causal coupling’ takes different meanings depending on the specific causality definition to which it is associated. In fact, it is well known that DC and PDC constitute terms into which the coherence and partial coherence functions may be decomposed [8,18], and that they may be related to the decomposition of diagonal elements of the spectral density matrix of the MV process and of its inverse, respectively [10]. As a consequence, the squared modulus of eDC (or nDC) from *y*_{j} to *y*_{i} quantifies the causal coupling intended as the fraction of the power spectrum of *y*_{i} (or of its part excluding instantaneous effects) that is due to *y*_{j} at each given frequency. The interpretation of the squared ePDC and nPDC functions is less meaningful, because they are related to functions of the inverse spectral domain, which does not provide clear physical interpretation of the measured quantities. On the other hand, the measures of (lagged or extended) direct causality determine the frequency domain connectivity structure of the considered MV process, in the sense that showing the presence of significant non-zero values for the squared nPDC or ePDC at a given frequency may be taken as an indication that a direct connection exists between two of the processes. From this point of view, nDC and eDC are less straightforward because they measure the ‘total’ (direct and indirect) effects from one process to another, so that the structural information cannot be elicited simply by checking whether nDC or eDC are non-zero.

## 3. Practical analysis of multivariate autoregressive processes

### (a) Model identification

The practical application of the theoretical concepts described in §2 on a set of *M* time series of *N* samples, *y*_{m}(*n*), *m*=1,…,*M*, *n*=1,…,*N*, measured from a physical system, is based on considering the series as a finite length realization of the MV process **Y**=[*y*_{1},…,*y*_{M}]^{T}, which describes the evolution of the system over time. Thus, the descriptive equation (2.1) is seen as a model of how the observed data have been generated, and a model identification procedure is applied to provide estimates of the coefficients and innovation variances to be used in equations (2.2)–(2.6) for computing the various frequency domain causality functions. The procedure for identification of the eMVAR model (2.1) is based on the prior identification of the corresponding strictly causal model,
3.1where **U**(*n*)=[*u*_{1}(*n*),…,*u*_{M}(*n*)]^{T} is a vector of white innovations with covariance matrix , and the coefficients of **A**(*k*) describe lagged interactions among the observed time series. The model (3.1) is traditionally used for the parametric frequency domain analysis of MV time series [20]. The basic difference from the extended model (2.1) is in the fact that the lag variable *k* starts from 1, so that possible instantaneous interactions among the *y*_{i} are not accounted for by the model coefficients. As a consequence, instantaneous interactions among the *y*_{i}, when present, result in a correlation among the innovations *u*_{i} such that the covariance **Σ** is not diagonal. In order to evaluate the relation between the two model representations, we note that, defining the matrix **L**=**I**−**B**(0)^{−1}, (2.1) can be rewritten as , which, compared with (3.1), leads to
3.2Hence, provided that an estimate of **L** (or equivalent of **B**(0)) is known, the identification of the extended model follows intuitively from the identification of the strictly causal model. Therefore, we propose the following identification algorithm:

— estimate

**A**(*k*),**U**(*n*) and**Σ**for the strictly causal MVAR model (3.1);— solve the instantaneous model

**U**(*n*)=**LW**(*n*) to estimate**B**(0)=**I**−**L**^{−1}and**W**(*n*);— exploit (3.2) to estimate coefficients and innovation covariance of the extended model as

**B**(*k*)=**L**^{−1}**A**(*k*) (*k*≥1) and**Λ**=**L**^{−1}**Σ**(**L**^{−1})^{T}.

The first step is easy to perform by means of classic regression methods (in the appendix, we describe the standard least-squares identification). As to model order selection, in this study, the order *p* was set at the value yielding the minimum of the Akaike figure of merit, AIC(*p*)=*N* [22]. The second identification step is less straightforward, because the instantaneous model may suffer from lack of identifiability, being related to the zero-lag covariance structure of the observed data, which is, *per se*, non-directional. In other words, using only the covariance information, one may find several combinations of **L** and **W**(*n*) that result in the same **U**(*n*), and thus describe the observed data **Y**(*n*) equally well. We propose two approaches to solve the instantaneous model, both based on using more information than the covariance information alone. The additional information consists of imposing *a priori* the structure of instantaneous causality for the first approach, and in exploiting the non-Gaussian distribution of the model innovations for the second approach. Accordingly, the overall identification procedures will be denoted, respectively, as eMVARis and eMVARng, to indicate that identification is performed assuming knowledge of instantaneous structure or non-Gaussianity of the innovation process **W**. The two approaches are described briefly in the following, and with more details in the appendix.

The eMVAR is approach for estimation of the instantaneous model **U**(*n*)=**LW**(*n*) is based on setting *a priori* the direction (though not the strength) of the instantaneous effects among the observed scalar processes [18]. This is achieved by imposing a causal ordering for the observed time series *y*_{m}, *m*=1,…,*M*, denoted as **κ**=[*κ*(1),…,*κ*(*M*)], such that no later series causes instantaneously any earlier series. The identification proceeds by: (i) permuting the estimated strictly causal residuals **U**(*n*) and their covariance **Σ** according to the imposed causal order; (ii) applying the Cholesky decomposition to estimate the permuted version of the extended innovation process **W**(*n*) and of the diagonal covariance **Λ**; and (iii) applying inverse permutation to get estimates of **W**(*n*) and of the instantaneous effects **B**(0). The eMVARng identification approach exploits the non-Gaussian structure of the observed data to overcome the identifiability problem of the instantaneous model [17]. It makes use of an algorithm recently proposed for the identification of instantaneous models [23], which is based on departing from the hypothesis of Gaussianity commonly assumed for the distribution of the independent processes **W**(*n*). The identification proceeds in two steps: first, independent component analysis (ICA) is performed on the estimated strictly causal innovations **U**(*n*), finding a mixing matrix **M** that represents an unordered and non-normalized version of **L**; second, appropriate row permutation followed by normalization is applied to **M**^{−1} to get an estimate of **L**^{−1}, and thus of **B**(0)=**I**−**L**^{−1}.

### (b) Model validation

Model validation refers to the use of a range of diagnostic tools that are available to check the adequacy of the estimated structure. The basic assumptions that need to be satisfied for the extended model (2.1) regard *whiteness* of the innovation process **W** and *independence* of its scalar components. These assumptions entail, respectively, that the stochastic variables *w*_{i}(*n*−*l*) and *w*_{j}(*n*−*m*) are independent for each *i*,*j*=1,…,*M* and for each *m*≠*l*, and that the stochastic variables *w*_{i}(*n*) and *w*_{j}(*n*) are independent for each *i*≠*j*. Note that whiteness of the extended innovations in **W** always corresponds to whiteness of the strictly causal innovations in **U**, because **W** and **U** differ only in the instantaneous structure. On the contrary, independence of **W** corresponds to independence of **U** only in the absence of instantaneous causality; in fact, violation of the assumption of independence of the strictly causal innovations is the reason leading to the introduction of the extended model in place of the strictly causal one. When the causal innovations are not strictly independent, the conditions underlying identifiability of the instantaneous model **U**=**LW** have to be verified in addition to the basic assumptions. These conditions consist of checking the adequacy of the imposed instantaneous structure, or testing non-Gaussianity of the extended innovations in **W**, respectively, when the eMVARis approach or the eMVARng approach is followed to identify the instantaneous model.

To perform model validation in practical applications, we propose the use of standard statistical tests, which are reviewed in Lutkepohl [20]. Whiteness of the model residuals representing estimates of the innovations **U** and **W** was tested using the Ljung–Box portmanteau test, which checks for the overall significance of groups of residual correlations. Independence of strictly causal or extended model residuals was tested by means of Spearman's rho test for independence; in particular, rejection of the test for the **U** residuals was taken as an indication of the significance of instantaneous effects among the observed series. If this was the case and if the eMVARng approach was used to estimate the instantaneous model, then non-Gaussianity of the extended residuals was tested by means of the Jarque–Bera test for non-normality of a vector stochastic process.

## 4. Numerical example

### (a) Theoretical formulation

In order to describe the four frequency domain causality measures defined in this study, and to illustrate their relationship with the different causality definitions, let us consider the MVAR process of order *p*=2, composed by *M*=4 scalar processes, generated by the equations:
4.1where the innovations *w*_{i} are independent white noises of unit variance (**Λ**=**I**). The diagonal elements of the coefficient matrix **B**(*k*) (i.e. *b*_{11}(1)=0.8√2,*b*_{11}(2)=−0.64, *b*_{22}(2)=−0.64) generate complex conjugate poles (with modulus 0.8 and phases ±*π*/4 and ±*π*/2) for the processes *y*_{1} and *y*_{2}, determining autonomous oscillations at 0.125 and 0.25 Hz, respectively. The off-diagonal elements of **B**(*k*) determine direct causal effects from one process to another, which are set according to the value of the parameter *δ*.

We consider the conditions *δ*=1 and *δ*=0, entailing the absence and presence of instantaneous causality, respectively. The corresponding causality patterns occurring among the processes in the two conditions are shown in figures 1*a*,*b* and 2*a*,*b*. When *δ*=1, all the imposed effects from one process to another are lagged. In this case, lagged direct causality is set over the directions (*b*_{21}(1)=1, *b*_{21}(2)=−0.5), (*b*_{32}(1)=0.5), (*b*_{42}(1)=0.5) and (*b*_{13}(2)=0.7) (figure 1*a*). The corresponding lagged causality relations are *y*_{1}⇒*y*_{2}, *y*_{2}⇒*y*_{3}, *y*_{2}⇒*y*_{4}, *y*_{3}⇒*y*_{1} (direct effects), and *y*_{1}⇒*y*_{3}, *y*_{1}⇒*y*_{4}, *y*_{2}⇒*y*_{1}, *y*_{3}⇒*y*_{2}, *y*_{3}⇒*y*_{4} (indirect effects) (figure 1*b*). Because instantaneous causality is absent, extended direct causality is equivalent to lagged direct causality (figure 1*a*), and extended causality is equivalent to lagged causality (figure 1*b*). When *δ*=0, the imposed causal effects are exclusively instantaneous from *y*_{2} to *y*_{3} and to *y*_{4} (*b*_{32}(0)=*b*_{42}(0)=0.5), exclusively lagged from *y*_{3} to *y*_{1} (*b*_{13}(2)=0.7), and mixed instantaneous and lagged from *y*_{1} to *y*_{2} (*b*_{21}(0)=1,*b*_{21}(2)=−0.5). The resulting causality relations are illustrated in figure 2*a*,*b*; note that, owing to the presence of instantaneous causality, extended direct causality differs from lagged direct causality (figure 2*a*), and extended causality differs from lagged causality (figure 2*b*).

The theoretical profiles of the four squared causality measures computed from the coefficient values according to equations (2.2), (2.3), (2.5) and (2.6) are depicted in figures 1*c*,*d* and 2*c*,*d*. Results show the close correspondence between each time domain causality definition and the associated frequency domain causality measure. Specifically, comparing (*a*) with (*c*) and (*b*) with (*d*) in figures 1 and 2, we note that, with if and only if , |*χ*_{ij}(*f*)|^{2}≠0 if and only if , if and only if *y*_{j}⇒*y*_{i}, and |*ξ*_{ij}(*f*)|^{2}≠0 if and only if . In the absence of instantaneous causality (*δ*=1), the extended measures provide the same quantitative information as the lagged measures, as seen in figure 1*c*, where ePDC overlaps with nPDC, and in figure 1*d*, where eDC overlaps with nDC. When instantaneous causality is significant (*δ*=0), lagged measures differ from extended ones, the former reflecting exclusively lagged causal effects and the latter reflecting either instantaneous and/or lagged effects (figure 2*c*,*d*). We see that the only non-zero nPDC profiles outside the diagonal of the matrix layout plot of figure 2*c* are and , reflecting the two lagged direct causality relations and shown in figure 2*a*; the larger number of extended direct relations is reflected by the non-zero squared ePDCs |*χ*_{21}(*f*)|^{2}, |*χ*_{13}(*f*)|^{2}, |*χ*_{32}(*f*)|^{2} and |*χ*_{42}(*f*)|^{2}. Besides the two direct causal effects, indirect lagged causality occurs also from *y*_{3} to *y*_{2} (figure 2*b*) and is reflected by the non-zero profile of the nDC (figure 2*d*); finally, the squared eDC is non-zero in a large number of connections (figure 2*d*), reflecting the broad nature (direct and/or indirect, instantaneous and/or lagged) of directional interactions described by this measure of extended causality.

### (b) Implementation on simulations

The feasibility of the two model identification approaches presented in §3 (i.e. eMVARis and eMVARng) was tested on practical realizations of the eMVAR process illustrated in §4*a*. Realizations of (4.1) lasting *N*=500 points were obtained by generating realizations of the extended innovations **W**(*n*), estimating the corresponding strictly causal innovations as **U**(*n*)=[**I**−**B**(0)]^{−1}**W**(*n*), and finally feeding a strictly causal model in the form of (3.1) with the estimated **U**(*n*) to get the series **Y**(*n*), *n*=1,…,*N*. The innovations, *w*_{m}(*n*), *m*=1,…,*M*, were generated as independent Gaussian white noises for testing the eMVARis identification approach, and as independent non-Gaussian white noises for testing the eMVARng approach. In the second case, non-Gaussianity was achieved by first generating *z*_{m}(*n*) as independent Gaussian white noises and then applying the nonlinear transformation *w*_{m}(*n*)=sign(*z*_{m}(*n*))|*z*_{m}(*n*)|^{q}, with exponent *q* chosen in the range [0.5, 0.8] or [1.2, 2.0] to yield, respectively, a sub-Gaussian or super-Gaussian distribution for *w*_{i}. For each of the two identification approaches, the analysis was performed under three different scenarios: (i) absence of instantaneous effects (*δ*=1 in (4.1)); (ii) presence of instantaneous effects (*δ*=0) with satisfied model assumptions; and (iii) presence of instantaneous effects (*δ*=0) with non-satisfied model assumptions. Scenarios (i) and (ii) were designed to test the ability of the two identification approaches to estimate frequency domain causality when instantaneous effects are missing and significant, respectively, whereas scenario (iii) was conceived to investigate the effects of wrong model assumptions on the various causality measures. Wrong assumptions were simulated in scenario (iii) by imposing, in eMVARis identification, an incorrect causal ordering for the observed series, i.e. any ordering different from the two correct orderings set in (4.1) for instantaneous causal effects (**κ**=[1,2,3,4] and **κ**=[1,2,4,3]), and imposing, in eMVARng identification, a Gaussian distribution for the extended innovations **W**.

Figures 3 and 4 depict simulation results obtained for 100 realizations of the MVAR process (4.1) using the two identification approaches under the three scenarios presented above. In the two figures, the top two rows show representative trends of expected and estimated frequency domain causality functions (one selected eDC and ePDC in figure 3, and one selected nDC and nPDC in figure 4); the third row shows the percentage of rejection of the four performed validation tests (whiteness and independence of **U** and **W** residuals, and Gaussianity of **W** residuals); and the fourth row shows the distribution of the error, estimated as the mean absolute value of the difference between estimated and true causal coupling, computed for each measure (eDC, nDC, ePDC and nPDC) and averaged over the full frequency range. When instantaneous effects were not present in the simulations (*δ*=1 in (4.1); figures 3*a* and 4*a*), both eMVARis and eMVARng identification approaches estimated correctly the imposed causal relations, as documented by the close adherence of theoretical and estimated measures and by the low error values. In this case, the assumptions of whiteness and independence of the residuals, and of non-Gaussianity of the extended residuals for the eMVARng approach, were satisfied in almost all simulations. In the presence of significant instantaneous effects (*δ*=0 in (4.1); figures 3*b*,*c* and 4*b*,*c*), the strictly causal residuals **U** were never independent, indicating the necessity of moving to the extended representation. When the additional assumptions required for identification of the instantaneous model (i.e. correct imposition of the direction of instantaneous effects for the eMVARis approach, and presence of non-Gaussian innovations for the eMVARng approach) were set, the estimation was successful, as demonstrated respectively in figures 3*b* and 4*b* showing low estimation errors and fulfilment of the validation tests. On the contrary, the identification was unsuccessful when the additional assumptions were not met, resulting in estimated frequency domain causality measures that markedly deviated from the expected profiles and in substantial estimation errors (figures 3*c* and 4*c*). In the case of eMVARng identification, the presence of Gaussian innovations led to failure of validation tests performed on the extended residuals (*iW* is high and *g* is very low in figure 4*c*). In the case of eMVARis identification, the imposition of the wrong direction for the instantaneous effects was not associated with unfulfilment of the validation tests (*iW* and *g* are near zero in figure 3*c*). This result is due to the fact that the models estimated with different instantaneous causal orderings are equally admissible, in terms of description of the data structure, as the considered one, and thus pass the validation tests. In such a case, selection of the model to be used should rely on prior knowledge of which is the most appropriate instantaneous structure, for example, on the basis of physical considerations.

## 5. Real data applications

This section describes two applications of the proposed framework for estimation of frequency domain causality in the presence of instantaneous causality, considering cardiovascular and neurophysiological MV time series, respectively. In the first application, the analysis is devoted to investigate the source mechanisms underlying cardiovascular oscillations, whereas the second application deals with determining the propagation directions of the alpha EEG rhythm. These two applications are chosen because their different goals entail utilization of different identification procedures and causality measures.

### (a) Cardiovascular variability signals

Short-term cardiovascular and cardiorespiratory variability analysis was performed on 15 young healthy subjects (25±3 years old) standing in the 60^{°} upright position after passive head-up tilt [24]. All subjects were in sinus rhythm and breathing spontaneously. The surface ECG, the finger photopletismographic arterial pressure signal and the respiratory nasal flow signal were acquired simultaneously and digitized with 1 kHz sampling rate and 12-bit resolution. Starting from these signals, the beat-to-beat variability series of heart period (HP, *y*_{1}), systolic arterial pressure (SAP, *y*_{2}) and respiratory flow (RF, *y*_{3}) were measured as follows (*M*=3): the *n*th HP value, *y*_{1}(*n*), was taken as the time interval between the *n*th and the (*n*+1)th R peaks of the ECG; the *n*th SAP value, *y*_{2}(*n*), was taken as the local maximum of the arterial pressure signal measured inside *y*_{1}(*n*); the *n*th RF value, *y*_{3}(*n*), was taken as the value of the nasal flow signal sampled at the onset of *y*_{1}(*n*) (i.e. at the *n*th R peak). A representative example is reported in figure 5*a*. According to these measurement settings, instantaneous effects may occur from RF to HP and to SAP, as well as from SAP to HP, but not over the reverse directions (see also Porta *et al.* [25], for the setting of instantaneous causality in cardiovascular variability analysis). This prior information about zero-lag effects was exploited to set the order ** κ**=[3,2,1] for instantaneous causality, thus allowing utilization of the eMVARis approach for model identification.

For each subject, synchronous stationary windows of the three time series, of length *N*=300 beats, were selected for the analysis. An example is reported in figure 5*b*; in this case, the optimal model order, selected with the Akaike criterion, was *p*=7. In model identification, whiteness of the model residuals was verified by means of the Ljung–Box test; independence (Spearman's rho test) was not satisfied between the strictly causal residuals of SAP and RF (*u*_{2} and *u*_{3}), indicating the significance of instantaneous effects. Utilization of the extended model led to independence between all pairs of extended residuals *w*_{1}, *w*_{2} and *w*_{3}. After validation, frequency domain analysis was performed computing the eDC and nDC from the model coefficients as in (2.2) and (2.5). The results in figure 5*c* also report the power spectral density function of each single time series (diagonal plots), which can be simply computed from the coefficients as the denominator of the squared eDC. The RF power spectrum (panel with *i*,*j*=3,3) shows a clear peak at the frequency of respiration (0.27 Hz for this subject); this activity is reflected by synchronous small peaks in the so-called high-frequency band (HF, ±0.04 Hz around respiration frequency) of the spectrum of HP (*i*,*j*=1,1) and SAP (*i*,*j*=2,2) series, which present also a prominent peak in the low-frequency band (LF, from 0.04 to 0.15 Hz) related to the so-called Mayer waves [26]. According to this scenario, we observe that lagged and extended causality are set from respiration to the cardiovascular variables in a unidirectional way, as the eDC and nDC evaluated in the HF band are substantial from RF to HP (*i*,*j*=1,3) and to SAP (*i*,*j*=2,3), and negligible over the opposite directions from HP to RF (*i*,*j*=3,1) and from SAP to RF (*i*,*j*=2,3). In the LF band, the analysis of causality shows a bidirectional interaction between SAP and HP, with eDC comparable to nDC from RR to SAP (*i*,*j*=2,1), and eDC higher than nDC from SAP to RR (*i*,*j*=1,2).

Figure 6 summarizes the results of the analysis extended to the 15 considered subjects, which display well-interpretable patterns of causality related to the generation of cardiovascular and cardiorespiratory oscillations. Because the RF time series does not present any meaningful oscillation within the LF band, only the analysis relevant to the interactions between HP and SAP are reported for this frequency band (figure 6*a*). The results suggest the existence of a closed-loop interaction between RR and SAP, as both the squared eDC and the squared nDC are non-negligible along the two pathways of the loop between *y*_{1} and *y*_{2}. This finding agrees with the known coexistence of regulatory mechanisms such as the so-called baroreflex operating from SAP to HP and the mechanical interaction exerted from HP to SAP [27,28]. The analysis performed in the HF band (figure 6*b*) evidences the key role of respiration in generating the oscillations observed in the two cardiovascular variables within this frequency range. Indeed, a strong unidirectional causal coupling is detected from RF to SAP and from RF to HP, whereas very low values are observed—both for the squared eDC and for the squared nDC—along all other possible directions of interaction. This finding is probably reflecting known phenomena, such as respiratory sinus arrhythmia and respiratory arterial pressure variability, which occur through the (direct or indirect) driving of respiration onto cardiovascular variability [29]. Our analysis further shows that instantaneous causality might play a role in determining the regulatory mechanisms observed in both frequency bands, as the squared eDC results are significantly higher, according to a paired Student's *t*-test, than the squared nDC when evaluated from SAP to HP at LF, and from RF to SAP and to HP at HF. As a consequence, we suggest the utilization of extended causality measures such as the eDC for the frequency domain analysis of cardiovascular interactions, to account for non-negligible zero-lag effects, which are representative of fast, likely of vagal origin, cardiovascular and cardiorespiratory mechanisms.

### (b) Electroencephalography signals

Analysis of neurophysiological signals was performed on 25 young healthy subjects (25±3 years old) lying with eyes closed in the relaxed awake state [30]. Multichannel EEG recordings were acquired through the standard 10–20 system (Fpz common reference) at 256 Hz sampling rate, then truncated to artefact-free synchronous stationary windows of 8 s and digitally filtered (0.3–40 Hz, zero phase-shift FFT band-pass filter). To reduce redundancy and favour MVAR model identifiability, the pre-processed signals were downsampled to 128 Hz and divided into five subsets describing front-to-back directions in the scalp (left-lateral: electrodes O1, T5, T3, F7; left-central: P3, C3, F3, Fp1; central: Pz, Cz, Fz; right-central: P4, C4, F4, Fp2; right-lateral: O2, T6, T4, F8) (figure 7*a*). For reducing the effects of reference electrode location, a modified common reference averaging procedure was applied for each subset, consisting of subtracting from each signal of the considered subset the average of all signals from the other subsets. An example of pre-processed signals is shown in figure 7*b* for the subset {Pz,Cz,Fz} of a representative subject. We used the convention of assigning suffix number to signals according to back-to-front scalp locations (i.e. *y*_{i} is measured at a more anterior location than *y*_{j} when *i*>*j*).

In this application, eMVAR model identification was performed according to the eMVARng approach, because no prior information about the direction of instantaneous effects can be inferred from synchronously sampled EEG signals. Hence, validation checked whiteness, independence and non-Gaussianity of the model residuals for each identified model. For each subject, we proceeded with frequency domain analysis only for the subsets for which all validation tests were passed. As the aim of this application was to infer the preferential directions of propagation for the alpha rhythm, we focused the analysis on computation of the direct causality spectral functions, which are useful for structure determination. Specifically, ePDC and nPDC were computed as in (2.3) and (2.6) from the eMVAR coefficients estimated from a given subset. Then, statistical significance inside the alpha band of the frequency spectrum (8–13 Hz) was assessed for each estimated function (squared modulus) through comparison with its corresponding significance threshold derived by means of the so-called causal FT surrogate data approach (the threshold was taken as the 95th percentile of the distribution of each squared ePDC or nPDC computed over 100 surrogate sets with absence of direct causality over the investigated causal direction; see Faes *et al.* [31] for details). An example of the analysis is reported in figure 7*c*, showing that the power spectral densities (diagonal plots) highlight the presence of dominant alpha activity at each channel, and how ePDC and nPDC analyses (off-diagonal plots) help in eliciting propagation patterns for this activity. In particular, the squared nPDC (figure 7*c*, dashed line) exceeds its significance threshold, inside the alpha band, over all three back-to-front directions of propagation (i.e. for *i*,*j*=2,1, *i*,*j*=3,1, and *i*,*j*=3,1), and only over one front-to-back direction (*i*,*j*=2,3). On the contrary, the extended direct causality relations resulting as statistically significant in the alpha band according to the squared ePDC (figure 7*c*, solid line) are in this case and .

In the overall analysis involving all subsets passing model validation and all subjects, results were collected counting the number of back-to-front connections and front-to-back connections for which the squared ePDC and nPDC resulted as statistically significant within the alpha band. As shown in figure 8*a*, the percentage of significant direct connections was balanced over the back-to-front and front-to-back directions when assessed by ePDC, and was prevalent over the back-to-front direction when assessed by nPDC. Furthermore, a significantly lower percentage of active front-to-back connections was detected using the nPDC (statistical analysis was assessed by means of the paired Student's *t*-test). Moreover, a prevalence of back-to front over front-to back active connections was observed in 19 of 25 subjects using the nPDC, and in only 12 of 25 subjects using the ePDC; the difference between the two paired proportions was found to be statistically significant according to the McNemar test. These results suggest that instantaneous causality, which in this application is likely to be related in major part to non-physiological phenomena such as volume conduction [32], may mask the detection of the direction of propagation of the alpha EEG waves. Indeed, a preferential direction of propagation in the alpha band, and specifically that going from the posterior towards the anterior scalp regions, could be evidenced only by looking at the patterns of lagged direct causality measured by the nPDC. The detection of a preferential back-to-front direction for the alpha EEG activity is in agreement with previous findings, and is explained by the occipital nature of the alpha oscillations, which are supposed to originate in the posterior visual cortex and then spread towards the central and frontal regions in the brain [13,31].

## 6. Discussion and conclusions

The framework for the assessment of frequency domain causality proposed in this study has been devised to deal with the adverse effects that instantaneous correlations among multiple time series may have on the estimation of time-lagged causal effects performed by means of the very popular MVAR approach. The impact of zero-lag correlations on the classic, strictly causal MVAR representation of multiple series has been recognized for a long time [4,20]. Recent studies have shown that neglecting instantaneous causality when it is non-trivial can lead to a misleading interpretation of causal effects, performed either in terms of the model coefficients [17], or in terms of their frequency domain representation [18]. The proposed solution consisted of moving to the eMVAR representation in which instantaneous effects are included in the model as zero-lag coefficients. This combination of a standard MVAR approach representing time-lagged influences with a structural equation modelling approach representing instantaneous influences allows one to describe simultaneously both the temporal precedence among variables that is actually detected as lagged in time, and the one that is hidden in the zero-lag correlations. Here, we have shown that explicit consideration of instantaneous causality leads to a distinction between the concepts of lagged causality and extended causality. Within the proposed framework, these two concepts are reflected in the frequency domain respectively by the nDC/nPDC functions measuring exclusive lagged effects, and by the eDC/ePDC functions measuring combined instantaneous and lagged effects. The reported analytical formulations and theoretical examples demonstrate that there exists a one-to-one correspondence between non-zero coefficients of the time-lagged causal effects and non-zero spectral profiles of the nPDC. As a generalization of this correspondence, the presence of extended causality (i.e. non-zero coefficients of instantaneous and/or lagged effects) entails non-zero spectral profiles for the ePDC. In a similar way, the nDC and eDC measure exclusive lagged effects and extended (instantaneous and/or lagged) effects, respectively, but accounting for both the direct pathway and all possible indirect pathways from one series to another in the MV representation. We note also that these functions reduce to the well-known DC and PDC functions [7,8], so that the proposed eMVAR framework reduces to the classic strictly causal MVAR framework in the absence of instantaneous correlations among the considered time series.

A major problem with the eMVAR approach, which has so far limited its utilization in practical time-series analysis, is that concerning model identification. In fact, unambiguous estimation of the coefficients describing instantaneous effects is a daunting task that cannot be solved by relying on the data covariance alone [20,23]. This issue is highly relevant for causality analysis, as taking instantaneous effects into account changes the estimation procedure for all autoregressive matrices. Therefore, any ambiguity in the estimation of the instantaneous model determines a lack of identifiability even for the time-lagged portion of the eMVAR model. The only way to guarantee identifiability is to incorporate in the identification procedure additional information with respect to the data covariance, so that uniqueness of the solution is guaranteed by exploiting such information. Note that closed form estimation methods such as the ordinary least squares, though recently proposed for eMVAR identification [33], are not feasible because they force identification of the instantaneous model to an arbitrary solution. In this study, we considered two possible approaches to deal with the identifiability problem: the eMVARis procedure exploits prior knowledge about the direction of instantaneous effects and makes use of the Cholesky decomposition [34]; and the eMVARng procedure exploits non-Gaussianity of the model innovations and makes use of ICA [23]. Whereas the first procedure is straightforward but does not solve the arbitrariness issue when the direction of instantaneous transfer paths cannot be deduced from the observed measurements, the second procedure is more objective but adds the requirement of non-Gaussian innovations. The feasibility of the two procedures was tested on realizations of the proposed theoretical examples, showing that, when the relevant additional identification assumptions are met, they reproduce to good approximation the expected frequency domain connectivity patterns.

As practical applications of the proposed framework, we considered biological MV time series where the autonomic nervous system is involved in signal interactions. The chosen applications reflect spheres of applicability for which the two presented eMVAR identification procedures are eligible. In cardiovascular time-series analysis, prior determination of the direction of instantaneous effects may be in agreement with physiological considerations: SAP might determine HP changes within the same cardiac beat as an effect of the fast vagal arm of the baroreflex (i.e. ); RF might produce modifications of the intrathoracic pressure that can be transferred within the same beat to SAP (i.e. ); respiratory centres are rapidly linked to vagal outflow, which, in turn, affects HP within the same beat (i.e. ). On the contrary, in EEG analysis, where no prior information about instantaneous causality may be exploited, the eMVARng procedure is recommended. Moreover, these two different applications suggested the use of different causality estimators. In cardiovascular analysis, where the interest was focused on determining the source of a given oscillation (LF or HF rhythm), we exploited the nDC/eDC measures that provide information on the full (i.e. over both direct and indirect pathways) information transferred as a function of frequency. In EEG analysis, where the focus was on determining the preferential propagation direction of the alpha waves, we used the nPDC/ePDC measures that are most indicated for determining the interaction structure within networks. Finally, the reported results suggest a different usefulness of lagged versus extended causality measures in the two applications. In cardiovascular analysis, where instantaneous causality is expected to be physiologically meaningful [29], the results were more interpretable, according to the known physiology, using the eDC than using the nDC. In EEG analysis, instantaneous causality is probably including confounding, non-physiological effects such as the signal cross-talk induced by volume conduction [32]; in this application, the detection of a propagation direction for the alpha waves was possible using the nPDC but not using the ePDC.

In conclusion, the proposed tool for the assessment of frequency domain causality has been shown to be useful, both theoretically and in practical analysis, for counteracting the problems arising from the presence of significant instantaneous effects in multiple interacting time series. The versatility of the tool has been demonstrated, showing that different identification procedures and causality measures may be chosen depending on the specific application of interest. We showed that unambiguous model identification in practical analysis is enabled by additional assumptions such as knowledge of the instantaneous causal structure or non-Gaussianity of the model innovations. However, we also showed that failure to satisfy these assumptions may lead to arbitrary estimates of model coefficients and connectivity patterns. Ongoing research in machine learning is witnessing several efforts aimed at proposing new methods for the identification of the instantaneous model that relax, at least partly, the model assumptions [35,36], as well as new variants that improve estimation accuracy and robustness against near-Gaussianity, scale invariance and outliers [37]. Because these modified estimators can be easily incorporated in our framework for eMVAR identification, future developments of this study will be directed to devising more accurate and general procedures for the estimation of frequency domain causality in the presence of instantaneous effects.

## Appendix A

Here, we formalize the aspects of model identification introduced in §3*a*. A simple, consistent and asymptotically efficient estimator for the strictly causal MVAR model (3.1) is the MV least-squares method. It is based first on representing (3.1) through the compact representation **Y**=**AZ**+**U**, where **A**=[**A**(1)⋯**A**(*p*)] is the *M*×*pM* matrix of the unknown coefficients, **Y**=[**Y**(*p*+1)⋯**Y**(*N*)] and **U**=[**U**(*p*+1)⋯**U**(*N*)] are *M*×(*N*−*p*) matrices and **Z**=[**Z**_{1}⋯**Z**_{p}] is a *pM*×(*N*−*p*) matrix having **Z**_{i}=[**Y**(*p*−*i*+1)⋯**Y**(*N*−*i*)] as *i*th row block (*i*=1,…,*p*). Given this representation, the method estimates the coefficient matrices through the well-known least-squares formula , the innovation process as the residual time series , and the innovation covariance as .

Once the strictly causal model is estimated, the second crucial step in the estimation of the eMVAR model (2.1) consists of the identification of the instantaneous model **U**=**LW**. The eMVARis approach is based on the initial knowledge of a causal ordering of the observed processes *y*_{1},…,*y*_{M}, i.e. of a vector of indices ** κ**=[

*κ*(1),…,

*κ*(

*M*)], such that instantaneous causality is allowed, for each

*j*<

*i*, from

*y*

_{κ(j)}to

*y*

_{κ(i)}but not from

*y*

_{κ(i)}to

*y*

_{κ(j)}. To identify the model, first an

*M*×

*M*permutation matrix

**P**is defined such that its elements are equal to one in the position (

*m*,

*κ*(

*m*)) and are equal to zero elsewhere (

*m*=1,…,

*M*). Second, the strictly causal residual covariance estimated as indicated above is permuted to obtain . The permuted covariance is then decomposed in Cholesky factors as , where the resulting matrix

**L**

_{p}is lower triangular by construction and represents the permuted version of the mixing matrix estimate. Returning back to the original ordering of the time series yields , from which instantaneous effects and extended residuals are easily estimated as and . Note that, because application of the Cholesky decomposition always leads to a lower triangular matrix, this approach is feasible only as a consequence of the initial imposition of a causal ordering for the instantaneous effects (so that the permuted instantaneous effects matrix is lower triangular).

The second considered approach, i.e. the eMVARng approach, identifies the instantaneous model **U**=**LW** according to the algorithm proposed in Shimizu *et al.* [23]. As a first step, the algorithm applies ICA on the residuals . ICA finds a mixing matrix **M** and a set of source processes **S** such that **U**=**MS**. The mixing matrix **M** can be identified as long as the distribution of **U** is a linear, invertible mixture of independent and non-Gaussian components. In fact, the problem is solved looking for the transformation **Q**=**M**^{−1} that makes the sources **S**=**QU** maximally non-Gaussian, and thus maximally independent. In principle, identification of the instantaneous model could result directly from ICA, i.e. we could take and . Unfortunately, the order and scaling of the independent components resulting from ICA are completely arbitrary; this means that the matrices **M** and **Q** are found up to permutation and scaling. Thus, the second step of the algorithm exploits the desired relationship between **L** and **B**(0), i.e. **L**^{−1}=**I**−**B**(0), to solve the permutation and scaling indeterminacies. Specifically, as **B**(0) is expected to have zero diagonal, the estimate of **L**^{−1} must have all ones on the main diagonal, and thus the row permutation of **Q**=**M**^{−1} corresponding to the correct model cannot have a zero on the diagonal. Hence, the matrix **Q** resulting from ICA is row-permuted to get a matrix with the largest possible elements on the diagonal; each row of is then divided by its diagonal element to get rescaling towards a new matrix with all ones on the diagonal. Finally, the model is solved by taking (or, equivalently, ), and .

The conditions for identifiability of the instantaneous model are that the variables in **W** are mutually independent and, if the eMVARng approach is used, then non-normally distributed. The assumption of non-normality follows directly from the utilization of ICA, which is known to find a mixture of independent components, and is successful only when these components have non-Gaussian distribution. Another condition for the algorithm proposed in Shimizu *et al.* [23] is that the matrix of the instantaneous effects **B**(0) corresponds to a directed acyclic graph, that is, a causal ordering exists for the variables in **W** such that no later variable causes any earlier variable. Under this condition, the row permutation performed on **Q** is unique, i.e. there is no other admissible permutation yielding a correct representation of the instantaneous model [23]. If the graph depicting **B**(0) is not acyclic, the obtained solution is not improper, but is rather the representation of one of the (many) admissible models fitting the observed data distribution [17]. In order to check the ‘degree of acyclicity’ of the estimated instantaneous model, one might proceed as follows [23]: find the identical row and column permutation of in terms of minimizing the sum of squares of elements on and above the diagonal, that is , where are elements of the permuted matrix; then, take as a score and issue a warning, indicating the probability of a non-unique representation of the instantaneous effects, if this score is larger than a predefined threshold (say, 5%). Note that the assumptions of non-Gaussianity and acyclicity are not required when the instantaneous model is identified through the eMVARis approach, because imposition of a causal ordering entails acyclicity of **B**(0) and avoids the utilization of ICA from which the requirement of non-Gaussianity stems.

## Footnotes

One contribution of 13 to a Theme Issue ‘Assessing causality in brain dynamics and cardiovascular control’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.