Causal Inference in Time Series via Supervised Learning

1 Introduction

Discovering temporal causal directions is an important task in time series analysis and has key applications in various fields. For instance, finding the causal direction that indicates that the research and development (R&D) expenditure X influences the total sales Y , but not vice versa, is helpful for decision making in companies. In addition, identifying causal (regulatory) relationships between genes from time series gene expression data is one of the most important topics in bioinformatics.

As a definition of temporal causality, Granger causality [Granger, 1969] is widely used [Kar et al., 2011; Yao et al., 2015]. According to its definition, the variable X is the cause of the variable Y if the past values of X are helpful in predicting the future value of Y .

Traditional methods for identifying Granger causality use regression models [Bell et al., 1996; Cheng et al., 2014; Granger, 1969; Marinazzo et al., 2008; Sun, 2008] such as the vector autoregressive (VAR) model and the generalized additive models (GAM). With these methods, we can determine that X is the cause of Y if the prediction errors of Y based only on its past values are significantly reduced by additionally using the past values of X. When the regression model can be well fitted to the data, we can infer correct causal directions. However, in practice, selecting an appropriate regression model for each time series data is difficult and requires a deep understanding of the data analysis. Therefore, it is not easy to identify correct causal directions with these model-based methods.

The goal of this paper is to build an approach to causal inference in time series that does not require a deep understanding of the data analysis. To realize this goal, we propose a supervised learning framework that utilizes a classifier instead of regression models. Specifically, we propose solving the problem of Granger causality identification by ternary classification, in other words, by training a classifier that assigns ternary causal labels (X → Y , X ← Y , or No Causation) to time series. In fact, several methods have already been proposed that perform classification to infer causal relationships from i.i.d. data, which have worked well experimentally [Bontempi and Flauder, 2015; Guyon, 2013; Lopez-Paz et al., 2015; 2017]. To solve causal inference in time series via classification, we formulate a feature representation that provides sufficiently different feature vectors for time series with different causal relationships. The idea for obtaining such feature vectors is founded on the definition of Granger causality: X is the cause of Y if the following two conditional distributions of the future value of Y are different; one is given the past values of Y and the other is given the past values of X and Y . To build the classifier for Granger causality identification, we utilize the distance between these distributions when preparing feature vectors. To compute the distance, by using kernel mean embedding, we map each distribution to a point in the feature space called the reproducing kernel Hilbert space (RKHS) and measure the distance between the points, which is termed the maximum mean discrepancy (MMD) [Gretton et al., 2007].

In experiments, our method sufficiently outperformed the model-based Granger causality methods and the supervised learning method for i.i.d. data by using the same feature representation and the same classifier. Furthermore, we describe how our approach can be extended to multivariate time series and show experimentally that feature vectors have a sufficient difference that depends on Granger causality, which demonstrates the effectivity of our proposed framework.

2 Granger Causality

Granger causality defines X as the cause of Y if the past values of X contain helpful information for predicting the future value of Y . Formally, it is defined as follows:

Definition 1 (Granger causality[Granger, 1969]) Suppose we have a stationary sequence of random variables {(Xt, Yt)} (t ∈ N), where Xt and Yt are on X and Y, respectively. Let SX and SY be observations of {X1, · · · , Xt} and {Y1, · · · , Yt}, respectively.

Granger causality defines {Xt} as the cause of {Yt} if

and states that {Xt} is not the cause of {Yt} if

To see if the two conditional distributions P (Yt+1|SX, SY ) and P(Yt+1|SY ) are identical, traditional methods [Bell et al., 1996; Granger, 1969; Marinazzo et al., 2008; Sun, 2008] use statistical testing to determine whether or not the two conditional means E[Yt+1|SX, SY ] and E[Yt+1|SY ] are equal, which is a much simpler problem than (1). For instance, in [Granger, 1969], the conditional means are represented by using the (V)AR model to compute the test statistic based on the prediction errors.

To represent the conditional means, these methods require us to use an appropriate regression model that can be well fitted to the data; however, such a model is difficult to select in practice. For this problem, we propose a novel approach that utilizes a classifier instead of regression models.

3 ProposedMethod

3.1 Classification Setup

Suppose the training data are N pairs of bivariate time series S1 , · · · , SN , where each time series Sj with the fixed length Tj consists of the observations of random variables {(Xj 1 , Y j 1 ), · · · , (XjTj , Y j Tj )} (j ∈ {1, · · · , N}). Here, each time series Sj has a causal label lj ∈ {+1, -1, 0} that indicates Xj → Y j , Xj ← Y j , or No Causation, where Xj = (Xj 1 , · · · , XjTj ) and Y j = (Y j 1 , · · · , Y j Tj ).

Using a function ν(·) that maps Sj to a single feature vector, we first train a classifier with {(ν(Sj ), lj )}Nj=1 . Then, the task of inferring the causal relationship in another time series S′ (test data) can be rephrased as assigning the label to ν(S′ ) by using the trained classifier.

Such a classification task can be extended to multivariate time series as detailed in Section 3.3.

3.2 Classifier Design

To build a classifier that assigns causal labels to time series, we formulate the feature representation ν(·). In what follows, we describe our ideas for obtaining feature vectors that are sufficiently different depending on Granger causality.

Basic Ideas for Granger Causality Identification

Simply by using the definition of Granger causality (Definition 1) 1 , for instance, we regard the causal label as X →Y if X is the cause of Y , and if Y is not the cause of X. Formally, we regard the causal label as 2

To assign causal labels to time series based on (2), (3), and (4), it is necessary to determine whether or not the two conditional distributions are identical. To represent information about conditional distributions, instead of using regression models, we utilize kernel mean embedding. Kernel mean embedding maps a distribution to a point in the feature space called the RKHS. Interestingly, when a characteristic kernel (e.g., a Gaussian kernel) is used, the mapping is injective: different distributions are not mapped to the same point [Sriperumbudur et al., 2010].

Suppose that kernel mean embedding maps conditional distributions P(Xt+1|SX, SY ), P(Xt+1|SX) and P(Yt+1|SX, SY ), P(Yt+1|SY ) to the points µXt+1|SX,SY , µXt+1|SY ∈ HX and µYt+1|SX,SY , µYt+1|SY ∈ HY , respectively, where HX and HY are the RKHSs. Then, when using a characteristic kernel, (2), (3), and (4) can be written as

To assign labels based on (5), (6), and (7), we only have to determine whether or not two points in the RKHS are the same over time t or, equivalently, to determine whether or not the distance between the two points in the RKHS, which is termed the MMD [Gretton et al., 2007], is zero over time t.

For this reason, to develop the classifier for Granger causality identification, in our feature representation, we utilize the MMDs, which are defined and estimated as follows.

Definition: Let kX and kY be kernels on X and Y, and HX and HY be the RKHSs defined by kX and kY , respectively. The MMD for the two distributions P(Xt+1|SX, SY ) and P(Xt+1|SX) is simply defined as the distance between µXt+1|SX,SY , µXt+1|SX ∈HX as follows

1 Note that since our approach is founded on Definition 1, which cannot address the case where there are latent confounders (i.e., unobserved variables that influence both X and Y ), as with the existing methods [Bell et al., 1996; Cheng et al., 2014; Granger, 1969; Marinazzo et al., 2008; Sun, 2008], this paper does not deal with such a case.

2 Although we do not consider the case where P(Xt+1|SX, SY ) ̸= P(Xt+1|SX) and P(Yt+1|SX, SY ) ̸= P(Yt+1|SY ) (i.e., X is the cause of Y , and Y is also the cause of X), we can straightforwardly address such a case by adding an extra label.

Similarly, MMD2Yt+1 is defined as the distance between µYt+1|SX,SY , µYt+1|SY ∈HY .

Estimation: The MMD can be estimated without using regression models and without performing a density estimation. At this point, the MMD is much more attractive than the Kolmogorov-Smirnov statistic [Chen and An, 1997] and the Kullback-Leibler divergence [Kullback and Leibler, 1951] since the former requires us to select regression models and the latter requires a density estimation, which is difficult when there are insufficient samples.

To estimate the MMD (8), we estimate the kernel mean embeddings of conditional distributions µXt+1|SX,SY and µXt+1|SX . As detailed in e.g., [Muandet et al., 2017], in general, the kernel mean embedding of the distribution is estimated by taking the weighted sum of the so-called feature mapping function. Specifically, when using the existing method called the kernel Kalman filter based on a conditional embedding operator (KKF-CEO) [Zhu et al., 2014], we can estimate µXt+1|SX,SY and µXT+1|SX by the weighted sum of the feature mapping ΦX:

where ΦX(xτ ) ≡ kX(xτ , ·) is a feature mapping function 3 , and wXY = [wXY 2 , · · · , wXY t-1 ]⊤ and wX = [wX 2 , · · · , wX t-1 ]⊤ (t > 3) are the real-valued weight vectors.

To compute the weight vectors wX and wXY , we employed KKF-CEO. In fact, KKF-CEO provides the algorithm needed to estimate wX from the observations SX for time series prediction. Therefore, we can compute wX by directly employing KKF-CEO. To estimate wXY from SX and SY , we simply used KKF-CEO with the product kernel kX · kY .

Although computing weight vectors by KKF-CEO requires the setting of several hyperparameters, they can be appropriately set for each time series by minimizing the squared errors between observations and the values predicted by KKF-CEO.

Applying (9) and (10) to (8), MMD2Xt+1 is estimated as

3 For instance, when using the Gaussian kernel kX(x, x′ ) = exp (-γ‖x-x′ ‖2 ) (γ > 0 is a parameter), the feature mapping becomes ΦX(x) = exp (-γx2 ) [1, √2γ/1! x, √(2γ)2 /2! x2 , · · · ]⊤.

Feature Representation

To build a classifier for Granger causality identification, we obtain the feature vectors by using the MMD pairs, where each pair dt = [̂ MMD2Xt+1 , ̂ MMD2Yt+1 ]⊤ is estimated by (11).

By using the MMD pairs, we can expect sufficiently different feature vectors to be obtained from time series with different causal labels. This is because whether or not the MMD becomes zero depends on the causal label as indicated by (5), (6), and (7). Although each MMD in dt cannot become exactly zero since it is a finite sample estimate, we can expect sufficiently different MMD pairs to be estimated from time series with different causal labels as intuitively shown in Fig. 1, which we confirm experimentally in Section 5.2.

To prepare dt for each t, given a time series with the length T, S = {(x1, y1), · · · , (xT , yT )}, we use its subsequence with the length W (W < T), i.e., {(xt-(W-1), yt-(W-1)), · · · , (xt, yt)} (t = W, · · · , T) 4 . As a result, we obtain the MMD pairs {dW, · · · , dT}.

Although we can directly use these MMD pairs as a single feature vector, such a feature vector has the dimensionality 2(T - W + 1), which depends on the time series length T . As feature vectors whose dimensionalities are the same for time series with different lengths, we utilize the mean of the MMD pairs. However, when simply using the mean (dW + · · · dT ) / (T -W + 1) as a feature vector, the feature vectors take the same value for the two sets of the MMD pairs whose empirical means are the same, but whose empirical distributions are different.

For this reason, to avoid mapping different distributions of the MMD pairs to the same feature vector, we again utilize kernel mean embedding. By using a different kernel function kD from kX and kY , we define our feature representation as

which is the mean over the feature mappings ΦD(dt) ≡ kD(dt, ·)5 .

4 By using shorter time series, we can reduce the time complexity when computing weight vectors by KKF-CEO (i.e., O(T 3 ) [Zhu et al., 2014]).

5 When using samples that are drawn directly from a distribution,

In (12), to compute the feature mapping ΦD(·), we employed random Fourier features (RFF) [Rahimi and Recht, 2007], which approximate a feature mapping as a lowdimensional vector of random features that are sampled from the Fourier transform of the kernel function. In experiments, we set the number of features m = 100 and obtained an mdimensional feature vector for each time series, where we observed no significant improvements in the inference accuracy when using a larger m.

3.3 Extensions toMultivariate Time Series

Finally, we describe how our approach can be extended to n-variate time series (n ≥ 3). We first present the feature representation for trivariate time series (i.e., n = 3) and we then address the case where n > 3.

Trivariate Time Series

Our feature representation for trivariate time series is founded on conditional Granger causality [Geweke, 1984], which can be applied to multivariate time series unlike Definition 1.

When using Definition 1 for trivariate time series, we can wrongly identify Granger causality. For instance, when there is no causal relationship between X and Y , if the third variable Z is their common cause, we can wrongly conclude that X is the cause of Y , or Y is the cause of X. This is because P(Yt+1|SX, SY ) ̸= P(Yt+1|SY ) or P(Xt+1|SX, SY ) ̸= P(Xt+1|SX) might hold due to the influence of Z.

To address the influence of Z, conditional Granger causality compares the two conditional distributions given SZ, i.e., the observations of the random variables {Z1, · · · , Zt}, each of which is defined on Z. Formally, it defines X as the cause of Y given Z if P(Yt+1|SX, SY , SZ) ̸= P(Yt+1|SY , SZ) holds; otherwise, X is not the cause of Y given Z.

We introduce causal labels based on conditional Granger causality. For instance, similarly to (2), we regard the causal label X →Y as

which can be rewritten as

where µXt+1|SX,SY ,SZ , µXt+1|SX,SZ , µYt+1|SX,SY ,SZ , and µYt+1|SY ,SZ are the kernel mean embeddings of P(Xt+1|SX, SY , SZ), P(Xt+1|SX, SZ), P(Yt+1|SX, SY , SZ), and P(Yt+1|SY , SZ), respectively.

Motivated by these expressions, we address the case where

Z might influence X and Y by adding ̂ MMD2Xt+1|Z and ̂ MMD2Yt+1|Z to the feature representation, where they are the MMD between µXt+1|SX,SY ,SZ and µXt+1|SX,SZ and the MMD between µYt+1|SX,SY ,SZ and µYt+1|SY ,SZ , respectively. For this reason, we extend the feature representation (12) simply by modifying dt as follows

the kernel mean embedding of the distribution is estimated by using the same weight values (e.g., 1/(T -W + 1) in (12)) [Muandet et al., 2017].

n-variate Time Series (n > 3)

Although it is possible to develop the feature representation for n-variate time series (n > 3) by adding the extra MMDs to dt, it becomes very difficult to prepare enough training data to train the classifier since the number of possible combinations of the common cause variables of the variable pair X and Y grows super-exponentially in n.

For this reason, we used the feature representation for trivariate time series. From n-variate time series, we infer a causal relationship between each variable pair X and Y in three steps. First, for each v ∈ {1, · · · , n - 2}, we obtain the feature vector from the observations of the triplet of the variables (X, Y , Zv). Second, based on each feature vector, we use a trained classifier to compute the probabilities of the three labels (X → Y , X ← Y , and No Causation). Finally, we assign the label with the highest average probability.

Addressing the case where there are more than one common cause variables is left as our future work.

4 RelatedWork

It is worth noting that one of the supervised learning methods for i.i.d. data, the randomized causation coefficient (RCC) [Lopez-Paz et al., 2015], also uses kernel mean embedding to obtain the features of the distribution that differ depending on the causal relationships. However, RCC and our method are designed to obtain different features of the distribution. Specifically, via kernel mean embedding, RCC obtains information about the marginal and conditional distributions, which are known to differ depending on the causal directions according to an assumption called the independence of cause and mechanism (ICM) [Janzing and Scholkopf, 2010]. In contrast, by using kernel mean embedding, our method measures the distance between the conditional distributions given past variable values since it becomes sufficiently different depending on Granger causality.

5 Experiments

5.1 Experimental Settings

We compared the performance of our method (hereafter referred to as the supervised inference of Granger causality (SIGC)) with the supervised learning method for i.i.d. data RCC [Lopez-Paz et al., 2015] 6 , with the Granger causality methods GCV AR [Granger, 1969] 7 , GCGAM [Bell et al., 1996] 7 , and GCKER [Marinazzo et al., 2008] 8 , which identify Granger causality by using the VAR model, the GAM, and kernel regression, respectively, and with transfer entropy TE [Schreiber, 2000] 9 , which infers causal relationships not by using regression models, but by performing density estimation.

For our SIGC, we used a random forest classifier 10 to make a fair comparison with RCC, which has achieved better performance with the random forest classifier than with the SVM [Lopez-Paz et al., 2015]. To prepare feature vectors, we used the Gaussian kernel as kX, kY , and kD and set the kernel parameter using the median heuristic, which is a well-known heuristic for selecting it [Scholkopf and Smola, 2001]. We set the parameter W in our method and the parameters in the existing methods to provide the best performance for each method in our synthetic data experiments. For our method, we selected W = 12. 1.0 0.8 0.6 0.4 Test Accuracy 1.0

6 https://github.com/lopezpaz/causation learning theory

7 http://people.tuebingen.mpg.de/jpeters/onlineCodeTimino.zip

8 https://github.com/danielemarinazzo/KernelGrangerCausality

9 https://github.com/Healthcast/TransEnt

10 The number of trees is selected from {100, 200, 500, 1000, 2000} via 5-fold cross validation.

0.8

5.2 Experiments on Bivariate Time Series Data Classifier Training 20 0.2 0 0.6 Test Accuracy SIGC RCC GC VAR

40

We trained a classifier to infer causal relationships from bivariate time series data. 0.4

As with the existing supervised learning methods [Bontempi and Flauder, 2015; Lopez-Paz et al., 2015; 2017], we used synthetic training data in both synthetic and real-world data experiments since there are few real-world data where the causal relationships are known. 20 0.2 0 40 60 80 100 120 Time Series Length

We generated 15, 000 pairs of synthetic time series with the length T = 42 so that there were 5.000 instances each with causal labels X → Y , X ← Y , and No Causation. Here, we used the following linear and nonlinear time series:

• Linear time series were sampled from the VAR model:

where τ = 1, · · · , P (P ∈ {1, 2, 3}) and EXt , EYt were sampled from the Gaussian distribution N(0, 1). To obtain time series with X → Y , we used the following coefficient matrix

where aτ , dτ were drawn from the uniform distribution U(-1, 1), and cτ ∈ {-1, 1}. Similarly, we prepared time series with X ← Y , and No Causation.

• Nonlinear time series were also similarly generated by using the VAR model with a standard sigmoid function g(x) = 1/(1+exp(-x)). For instance, we prepared time series with X → Y so that Yt depended on {[g(Xt-τ ), Yt-τ ]⊤}Pτ=1 while Xt depended only on {Xt-τ }Pτ=1 .
• Finally, we scaled each time series with mean 0 and variance 1.

Synthetic Time Series

Next, we tested the performance of our method for inferring causal relationships (X → Y , X ← Y , and No Causation) from synthetic time series data. We used the following linear and nonlinear test data:

• Linear Test Data: We prepared 300 pairs of linear time series so that the numbers of time series with X → Y , X → Y , and No Causation were 100. In a similar way to the linear time series in the training data, each time series was sampled from the VAR model (13) although several parameter settings were different (e.g., the noise variance was given as p ∈ {0.5, 1.0, 1.5, 2.0}).

GC GAM

SIGC

RCC TE GC GAM GC VAR GC KER GC VAR 60 80 100 120 Time Series Length TE GC KER Figure 2: Test accuracies for 300 pairs of time series data against time series length (left: linear test data; right: nonlinear test data). Means and standard deviations (error bars) are shown for our method and RCC based on 20 runs with different training data.

• Nonlinear Test Data: We used 300 pairs of nonlinear time series, where there were 100 time series with X → Y , X → Y , and No Causation in each dataset. We generated nonlinear time series with X → Y by

where the noise variables EXt , EYt were sampled from N(0, 1). Similarly, we prepared nonlinear time series with X ← Y . To prepare time series with No Causation, we simply ignored the exponential term in (15).

Using linear and nonlinear test data, we compared the performance of our method with that of the existing methods. Fig. 2 shows the test accuracies. Note that for SIGC and RCC, we show the means and the standard deviations (error bars) in 20 experiments with different training data since these methods use randomly generated training data.

The performance of the Granger causality methods (GCVAR, GCGAM, and GCKER) depended on whether or not the regression model could be well fitted to the data. For instance, since the VAR model could be well fitted to linear test data, GCV AR performed well on linear test data although it worked badly on nonlinear test data. In addition, with nonlinear test data, GCKER was less accurate than GCGAM because the time series was too short for us to perform kernel regression. Similarly, TE worked poorly since the time series was too short for us to perform density estimation.

In contrast, our method worked sufficiently well on linear and nonlinear test data. The main reason for the good performance lies in our feature representation. This can be seen from a comparison with the supervised learning method RCC since it prepares training data in the same way as our method.

To verify our feature representation, we confirmed experimentally that feature vectors are sufficiently different depending on causal labels. To do so, we used nonlinear test data to plot a histogram of the MMD pairs {dt} that were used to compute the feature vector for each time series. Fig. 3 shows the results. Since each MMD in dt is a finite sample estimate, no MMD becomes exactly zero. However, we can see that the MMDs became sufficiently different according to the causal labels. We obtained similar results with linear test data although we have omitted them due to space limitations.

In fact, since the MMD pairs are sufficiently different, we can assign the causal label by taking another approach, i.e., an unsupervised approach that uses no training data, which outputs the causal label in two steps. First, using the MMD pairs, the two statistical tests are performed to determine if the mean of ̂ MMD2Xt+1 is zero and if the mean of̂ MMD2Yt+1 is zero. Then, by using the two p-values and some threshold value (significance level), we can assign a causal label (X → Y , X ← Y , or No Causation) to each time series.

However, we confirmed experimentally that the performance of such an unsupervised approach depended greatly on the threshold value. Furthermore, it was less accurate than our method (e.g., on nonlinear test data with the length T = 250, its test accuracy was 0.810 (not shown in Fig. 2) while our method achieved 0.966) although we selected the threshold value that provided the best performance. These results suggest the effectiveness of our supervised learning approach, which can obtain the decision boundary needed to determine the causal label by training a classifier.

Real-world Time Series

We tested our method by using real-world time series. To improve the reliability of the experiment, we used the following two test datasets:

• The first test dataset was composed of five pairs of bivariate time series downloaded from the Cause-Effect Pairs database [Jakob, ], whose true causal relationships are reported in [Jakob, ] as X → Y for three pairs and as X ←Y for the others. For instance, the River Runoff is a bivariate time series concerning average precipitation X and average river runoff Y , and the true causal relationship is regarded as X → Y .
• Using the above five real-world time series, we prepared a second test dataset that consisted of subsequences in each time series. To prepare the subsequences, we sim-

ply chopped each time series into multiple subsequences of length T = 200.

As regards training data, we used synthetic time series that we prepared in the same way as those for synthetic data experiments.

Table 1 shows the result for each test dataset. Note that we have omitted RCC from the top table in Table 1 because it showed different outputs in 20 experiments where different training data were used as in the synthetic data experiments, while our SIGC always output the same causal directions. As shown in Table 1, our SIGC outperformed the other existing methods regardless of the time series length T .

5.3 Experiments onMultivariate Time Series Data

We tested SIGCtri, which utilizes a feature representation for trivariate time series. For classifier training, we used synthetic trivariate time series that were generated in a similar way to those used in the experiments described in Section 5.2. As test data, we used the following time series gene expression data:

• We used the Saccharomyces cerevisiae (yeast) cell cycle gene expression dataset collected by [Spellman et al., 1998]. By combining four short time series that were measured in different microarray experiments, we prepared a time series with the length T = 57, where the number of genes was n = 14. To determine the true causal relationships between the genes, we used the gene network database KEGG [KEGG, 1995].

Since the number of non-causally-related gene pairs was much larger than the number of causally-related gene pairs, we evaluated the performance of each method in terms of the macro and micro-averaged F1 scores rather than test accuracy.

Table 2 shows the results. Since the data were measured in different microarray experiments, all the methods could not sufficiently work well. However, our SIGCtri worked better than the existing Granger causality methods. It also performed better than SIGCbi, which uses the feature representation for bivariate time series, thus indicating that it is important to consider the influence of the common cause variable as described in Section 3.3.

6 Conclusions

We have proposed a classification approach to Granger causality identification. Whereas the performance of the model-based methods depended hugely on whether the regression model could be well fitted to the data, our method performed sufficiently well by using the same feature representation and the same classifier (random forest classifier). Furthermore, we demonstrated experimentally the reason for such good performance by showing a sufficient difference between the feature vectors that depends on Granger causality. These results demonstrate the effectiveness of classification approaches to Granger causality identification.

Addressing complicated real-world scenarios (e.g., inferring the causal directions that change over time) constitutes our future work.