Gaussianity Measures for Detecting the Direction of Causal Time Series

1 Introduction

Identifying the true direction of a causal time series is an interesting problem for the evaluation of novel methods for causal inference [Peters et al., 2009]. Furthermore, solving this problem can provide new insights about the asymmetries between past and future in physical systems [Janzing, 2010].

We conjecture that, when a linear process is appropriate for modeling a time series with non-Gaussian innovations, the residuals of a linear fit to the time-reversed series are more Gaussian than the residuals in the chronologically (forward in time) ordered series. This property can be derived for series that follow a causal AR process, assuming that the noise is non-Gaussian and that the cumulants of the distribution of the noise are defined. Based on this property it is possible to use Gaussianity measures to detect the true direction of a time series. The proposed method works by fitting a linear model to the time series in the original ordering and in the inverted ordering and then selecting the direction in which the corresponding residuals are less Gaussian.

∗This research has been supported by the Spanish Direcci´on General de Investigaci ´on under project TIN2010-21575-C02-02. AS acknowledges partial financial support from Shizuoka University (Japan).

Different metrics can be used to measure the distance of the empirical distribution of the residuals to the Gaussian. A simple method is to use an estimate of a cumulant of order higher than two to make the decision. More sophisticated metrics, such as the Maximum Mean Discrepancy (MMD) have been recently proposed [Gretton et al., 2007]. The performance of the novel causal inference rules is evaluated in experiments with data simulated from AR and ARMA processes and with real-world data formed by time series of EEG measurements of brain activity. In these experiments, the proposed approach outperforms state-of-the-art decision methods that work by quantifying the independence of the empirical residuals with respect to the previous time series values [Peters et al., 2009].

2 Reversibility in Time Series

The time series {Xt}∞t=-∞ is said to be reversible when the vectors {Xt1 , Xt2 , . . . Xtn } and {X-t1 , X-t2 , . . . X-tn } have the same joint distribution for all n > 0 and ti, where 1 ≤ i ≤ n [Lawrance, 1991]. Identifying the true direction of a time series is only possible when the series is not reversible. Consider the autoregressive-moving average linear process

where {ϵt} is i.d. white noise [Brockwell and Davis, 1991]. The previous definition of reversibility implies stationarity [Lawrance, 1991]. If the distribution of ϵt is Gaussian then {Xt} is time-reversible. This can be derived using the fact that, when the noise is Gaussian, the joint distribution of {Xt1 ,Xt2 , . . .Xtn } for any n > 0 and ti, where 1 ≤ i ≤ n, is multivariate Gaussian. The symmetry of the covariance matrix guarantees that the process is strictly stationary and time-symmetric. Weiss showed that if an ARMA process 1 is reversible, then {ϵt} is normally distributed [Weiss, 1975].

An alternative definition of time reversibility not restricted to strictly stationary processes is proposed in [Peters et al., 2009; 2010]. Consider a causal ARMA process, in which ϵt is independent of Xt-k ∀k > 0, ∀t ∈ Z. The causal ARMA process is said to be reversible when there is an i.i.d. sequence {˜ϵt} such that Xt = ∑pi=1 φiXt+i + ∑qj=1 θj˜ ϵt+j + ˜ ϵt, with

1 The process must satisfy p > 0 or, if p = 0, the cases θn = θq-n and θn = -θq-n are excluded.

{˜ϵt} being independent of Xt+k, ∀k > 0, ∀t ∈ Z. They then show, using the Darmois-Skitovich theorem, that a causal ARMA process with i.i.d. noise and p > 0 is time-reversible if and only if the noise is normally distributed. Finally, Peters et al. propose to identify the true direction of a time series by fitting an ARMA model and testing whether the residuals of the model are independent of the previous time series values in one direction and dependent in the opposite one.

Our approach to the problem of directionality detection is different. Assuming a linear time-series model, we show that the residuals of the time-reversed series are more Gaussian than the residuals of the original series. Therefore, once we have obtained the two series of empirical residuals, we can use measures of Gaussianity to identify the correct direction.

3 Main Result

To simplify the presentation, the main result is derived for a stationary causal AR(1) process

{ϵt} is an i.d. white noise process, which is not necessarily Gaussian. The noise is independent of the delayed values of the process: ϵt ⊥ Xt-k, ∀k > 0, ∀t ∈ Z. The condition of stationarity is fulfilled if |φ| < 1. The case φ = 0 is trivially time-reversible (Xt = ϵt) and will not be considered further.

If the process is reversible, it can be expressed as

where {˜ ϵt} is a white noise sequence and ˜ ϵt ⊥ Xt+k, ∀k > 0, ∀t ∈ Z. The coefficient φ is the same as in the forward equation because it is the one-lag autocorrelation φ = cov(XtXt-1)/var(Xt) = cov(XtXt+1)/var(Xt), which is symmetric in time.

If the process is not reversible, we define the time-reversed residuals as

In this case, {˜ ϵt} is not a white noise sequence and ˜ ϵt is in fact dependent of Xt+k, for some k > 0, ∀t ∈ Z.

Assuming that the cumulants of the white noise process exist, the cumulants of ˜ ϵt can be expressed in terms of the cumulants of ϵt:

The derivation of this expression can be found in Appendix A. Figure 1 shows plots of cn(φ) as a function of φ for n = 1, . . . , 14. The following result can be stated for causal AR(1) processes with φ ̸= 0 and |φ| < 1:

Theorem 3.1. The cumulants of ˜ϵt of order higher than 2 are smaller in magnitude than the cumulants of ϵt, if these exist.

For a Gaussian distribution, the cumulants of order higher than two are zero. Therefore, Theorem 3.1 implies that the cumulants of ˜ ϵt are closer to the cumulants of a Gaussian than the cumulants of ϵt. In this precise sense, the distribution of

˜ ϵt is closer to a Gaussian than the distribution of ϵt. This Gaussianization effect can be observed even when the higher order cumulants diverge. In fact, it obtains even when none of the cumulants are defined, as in the Cauchy distribution. This effect is stronger for AR(1) processes in which |φ| is in the vicinity of the golden ratio conjugate: φ ≈ ±(√5 - 1)/2 (see appendix A and the plots in Figure 1).

Theorem 3.1 allows us to determine the correct direction of the series {Xt}Tt=1 by measuring distances to the Gaussian of the distributions of the residuals: Assuming that the AR(1) model is correct, one performs a fit to the series and calculates the empirical residuals in the original ordering, {et}Tt=2 , and in the inverted ordering of the series, {˜ et}T -1 t=1 , where

t = 2, . . . , T and ˆ φ is the empirical correlation between Xt and Xt-1. Note the difference between {et}, the empirical residuals, and {ϵt}, the i.d. white noise process that was used to generate the time-series values. The correct direction of the series is the one in which the corresponding empirical residuals are less Gaussian.

Theorem 3.1 can be readily generalized to causal AR(p) processes and is expected to hold also for causal ARMA and vector AR processes. An analysis of simulated and real-world data provides strong evidence for this.

4 GaussianityMeasures

We consider two measures of Gaussianity for detecting the direction of causal time series that follow a linear model. The first one is directly based on the inequality given by Theorem 3.1 for a particular order of the cumulants. We focus on the fourth order cumulant because κ4 is the cumulant with lowest order which is non-zero for symmetric non-Gaussian noise. In leptokurtic distributions, the estimates of κ4 have lower variance than the estimates of cumulants of higher order.

Given a sample {Yi}Ni=1 , the unbiased estimator of κ4 with lowest variance is the fourth k-statistic [Kendall et al., 1994]:

where N [4] = N (N - 1)(N - 2)(N - 3) and S1, . . . , S4 are the first 4 power sums, that is, Sr = ∑Ni=1 Y r i . The inference rule based on κ4 selects the direction for which the empirical residuals have lower values of (3).

The deviations from Gaussianity may also appear at higher order cumulants. In fact, it is possible to build non-Gaussian random variables with a finite number of zero cumulants, including κ4. Therefore, properly designed metrics on the space of probability distributions should be a better way to quantify the deviations of the residuals fromGaussianity and to provide a more effective method for detecting the direction of the time series. We focus on metrics based on a Hilbert space embedding for probability distributions [Gretton et al., 2007; Sriperumbudur et al., 2010]. Let p and q denote two density functions defined on R and let F be the unit ball in a universal reproducing kernel Hilbert space H defined on R with kernel k(·, ·). The Maximum Mean Discrepancy (MMD) between p and q within the set of functions F is

where μp = Ey∼p[k(y, ·)] and μq = Ez∼q [k(z, ·)] are the corresponding mappings of p and q onto H and ‖ · ‖H denotes the norm operator in H [Smola et al., 2007]. MMD vanishes when p = q and is positive otherwise [Gretton et al., 2007].

Given a sample {Yi}Ni=1 with zero-mean and unit standard deviation, we can compute the MMD statistic between the standard Gaussian distribution and the empirical distribution of the sample. For this, we define

where δ(y) represents a point probability mass at y = 0, and the target function q is the N(0, 1) density. Using a Gaussian kernel k(y, z) = exp[- 1 2σ2 (y - z)2 ] of width σ, we obtain the square of the MMD statistic:

where 〈·, ·〉H denotes the dot product operator in H and k′ (y) = σ √1+σ2 exp[- 1 2(1+σ2 ) y2 ]. The MMD decision rule selects the direction of the time series in which the value of (4) for the empirical residuals is lower.

5 RelatedWork

PC-LINGAM [Hoyer et al., 2008] combines conditional independence methods with ICA based approaches for causal inference. This method discriminates between models in the same d-separation-equivalence class by selecting the model that maximizes the ICA target function

where p is the empirical distribution of the standardized model residuals (zero mean and unit standard deviation). This target function is a specific case of the MMD distance in which the set F contains only the absolute value function | · |.

A related method is described in [Daniuˇsis et al., 2010]. This technique employs Gaussianity measures for inferring causal directions between two random variables. The variable whose marginal distribution is closer to a Gaussian is selected as the cause and the other one as the effect. The distance to the Gaussian is computed using an empirical estimate of the entropy. This method cannot be directly applied to time series because, if the analyzed process is stationary, both cause (Xt) and effect (Xt+1) have the same distribution.

In [Zhang and Hyv¨arinen, 2009], a class of acyclic causal models are investigated. In these models, the causal relations among the observed variables are nonlinear but the effect of the disturbances is linear. These authors show that, when the disturbances are additive, the correct causal model is the one with the lowest entropy for the disturbances. This agrees with the results of our investigation: For linear models, (2) implies that the differential entropy of the forward residuals is lower than that of the reversed residuals. However, in the problems investigated, entropy-based measures of Gaussianity, such as those used in [Daniuˇsis et al., 2010], achieve lower accuracies than the proposed methods (results not shown).

6 Experiments

The performance of the decision rules based on MMD and κ4 is investigated in experiments with simulated and real-world data. The statistical test described in [Peters et al., 2009] is used as a benchmark for comparison. This method is based on determining whether the empirical residuals are independent of the previous time series values in one direction and dependent in the opposite direction. The degree of dependence is measured using the p-value given by a statistical test based on the Hilbert-Schmidt Independence Criterion (HSIC) [Gretton et al., 2008]. The smaller this p-value is, the larger the level of dependence in the data. The direction predicted by the HSIC test is the one for which the dependencies between the model residuals and the previous time series values up to a fixed lag γ are weaker. For example, when γ = 2 the HSIC method checks for (i) dependencies between et and Xt-1,Xt-2 and (ii) dependencies between ˜ et and Xt+1,Xt+2. The input to each decision rule (MMD, κ4 and HSIC) are the empirical residuals and the time series values (the latter ones only for HSIC). These input series are standardized so that they have zero mean and unit standard deviation.

6.1 Experiments with Simulated Data

The protocol of these experiments is similar to the one used in [Peters et al., 2009]. First, we simulate time series that follow AR(1) processes, using different values of φ and with ϵt ∼ |Z|r sgn(Z ), where Z ∼ N (0, σ2 ) and σ is such that ϵt has unit standard deviation. The parameter 0.1 ≤ r ≤ 2 determines the shape of the noise distribution. For r = 1 ϵt follows a Gaussian distribution. In the region r > 1, the noise is leptokurtic. For values r < 1, the distribution of ϵt is platykurtic and bimodal. The density of ϵt is plotted in the top row of Figure 2 for r = 0.1, 0.75, 1.25 and 2. The bottom row in this figure displays the density functions of the timereversed residuals in the AR(1) model when φ = (√5 - 1)/2. The reversed residual densities are closer to a Gaussian than the original ones.

In general, both MMD and HSIC use Gaussian kernels whose width σ is σ2 = m/2, where m is the median distance between sample points. However, MMD does not perform well using this rule when the modes of the distribution for r < 1 are far apart. In this specific case, the median distance between sample points measures mainly the distance between the modes, rather than the characteristic scale of the data. To avoid this problem, we only consider positive values of the residuals when computing m in MMD for r < 1. The results of the HSIC test are similar using either the original or the modified rule. In practice, the performances of MMD and of the HSIC test do not strongly depend on the particular choice of σ as long as reasonable values are considered.

Using different values of φ and r, we sampled 1000 series of length 100 and evaluated the accuracy of each decision method. The accuracy is defined as the fraction of series for which the method predicts the correct direction. The top row in Figure 3 displays the results for r = 0.75, 1.25, 1.5 when φ is varied between -0.98 and 0.98. For r < 1, all the decision rules perform similarly. However, when r > 1 the method based on κ4 obtains worse results than MMD and HSIC. The reason for this is the larger sampling variance of (3) when the residual distribution is heavy-tailed. When r > 1 and φ is small, MMD consistently outperforms HSIC. All the methods perform similarly as random guessing for |φ| close to 0 or 1. The reason is that, for these parameter values, the process is close to being time-reversible. The highest accuracies are obtained for φ ≈ ±(√5 - 1)/2. This is the value of φ for which the Gaussianization effect is the largest.

We also simulate ARMA(2,2) models Xt = φ1Xt-1 + φ2Xt-2 + θ1ϵt-1 + θ2ϵt-2 + ϵt, where the distribution of the noise has the same functional form as before. The parameters are φ1 = 0.9, φ2 = -0.3, θ1 = -0.29 and θ2 = 0.5 and r is varied between 0.1 and 2, as in [Peters et al., 2009]. For each value of r, 1000 series of length 200 are generated. The analysis is made also for other choices of φ1: φ1 = 0.3 and

φ1 = (√5 - 1)/2. The lag used in the HSIC independence test is γ = 1, except in the experiments with φ1 = 0.3, where γ = 2 performs better. The bottom row in Figure 3 displays the results obtained by MMD, HSIC and the method based on κ4. MMD performs best, followed by the κ4 method. The performance of HSIC for r > 1 and φ1 = 0.3 is rather poor. All the methods perform as random guessing for r = 1. In this case, the noise distribution is Gaussian and the time series are time-reversible. Similar results are obtained for ARMA models of different orders and different parameter values.

6.2 Experiments with Real-world Data

The performance of the proposed rules is also evaluated in experiments with time series of EEG measurements of brain activity. Following [Peters et al., 2009], we use 1180 time series of length 500 from the publicly available dataset [EEG data, 2010]. Different ARMA models of order up to (5,5) are calibrated to each time series in each possible direction. The best model is selected using the AIC method. For each time series, a Jarque-Bera test of Gaussianity is applied on the empirical residuals of the correct direction of the series [Peters et al., 2009]. We discard the series for which this test does not reject the Gaussian hypothesis: When the noise is Gaussian the series is time-reversible and it is not possible to predict the correct ordering. Finally, we follow [Peters et al., 2009] and select the final set of series on which we make a decision. For this, we use two p-values. The first one, p, is given by an HSIC test that checks for dependencies between et and Xt-1. The second one, ˜ p, is given by an HSIC test for dependencies between ˜ et and Xt+1. We select those time series that satisfy min(p, ˜ p) < α and max(p, ˜ p) > α and additionally, max(p, ˜ p) - min(p, ˜ p) > δ, where α and δ are two parameters whose value are fixed beforehand. We make a decision only when the difference between p-values is at least larger than δ and exactly one p-value is lower than α. Typically, α is small and δ large. The first filter (the one based on α) guarantees that the causal ARMA process provides a good description for some ordering of the series. The second filter (the one based on δ) selects series for which the HSIC method is more confident about its decision. For the HSIC test, the lag used in the experiments is γ = 1, which gives the best results. Both MMD and HSIC use Gaussian kernels whose width σ is σ2 = m/2, where m is the median distance between sample points.

Figure 4 shows the accuracy of each method for α = 0.02, 0.04, 0.06 and 0.08 and values of δ in the interval [0, 0.5]. The number of series selected varies from 164 for α = 0.8 and δ = 0 to 17 series for α = 0.02 and δ = 0.5. MMD has the best overall performance. The method based on κ4 obtains rather poor results. The reason is that the noise in the EEG data usually presents heavy tails. These heavy tails mean that the variance of the estimator of κ4 given by (3) due to the finite size of the sample is rather large. The accuracy of HSIC increases with δ and in some cases HSIC outperforms MMD (for large δ). Note that, by increasing δ, we are selecting the time series in which the HSIC method is more confident about its prediction. Hence, larger values of δ favor the HSIC method. The MMD approach could be modified to generate p-values under the null hypothesis that the distribution of the residuals is Gaussian [Gretton et al., 2007]. These p-values could be used to select the time series for which MMD is more confident about its prediction. However, this possibility has not been addressed in the present investigation.

7 Concluding Remarks

Using the properties of cumulants we have shown that for a causal AR(1) process, the residuals {˜ ϵt} of a time-reversed linear stochastic process are closer to a Gaussian than the noise {ϵt} used to generate the process in the chronologically ordered direction. The proof is based on showing that all the cumulants of ˜ ϵt of order higher than 2 are smaller in absolute value than the cumulants of ϵt. Using this property, we have designed a rule to identify the correct direction of a time series that admits a linear model and has non-Gaussian noise. This rule is based on measuring the distance of the residual distribution to the Gaussian. The correct direction is identified as the one for which the model residuals are less Gaussian. Experiments on several AR and ARMA processes with different distributions for the innovations and different parameter values confirm the excellent performance of the proposed rule. Similar results are obtained using real-world time series of EEG data. The performance is good even when the cumulants of the noise are not defined. This is probably related to the fact that the cumulants of a finite sample are finite.

A Derivations

The product cumulants of the random variables ϵt and Xt-1 are defined as

with κ00(ϵt,Xt-1) = 0, and a and b arbitrary constants. Since ϵt and Xt-1 are independent, only the single-variable cumulants in the expression above are different from zero:

where δ is the Kronecker delta and κm(Xt-1) = κm(Xt) by stationarity. The characteristic function of ˜ ϵt and Xt+1 can be written in terms of the characteristic functions of ϵt and of Xt:

where a and b are purely imaginary constants. The product cumulants of ˜ ϵt and Xt+1 can also be expressed in terms of the cumulants of ϵt and Xt :

for all m > 0, n > 0, m+ n > 0. Assuming stationarity, we can use the relation

and express the product cumulants in terms of the cumulants of ϵt alone

The cumulants of ˜ ϵt and Xt are particular cases of (5) with m = 0 and n = 0, respectively. Therefore, the cumulants of ˜ ϵt can be expressed as a function of the cumulants of ϵt:

and n > 0. For causal AR(1) processes with φ ̸= 0 and |φ| < 1 we obtain c1(φ) = c2(φ) = 1, |cn(φ)| < 1, ∀n > 2, limn→∞ cn(φ) = 0 and