{'fraction_non_alphanumeric': 0.039983952928590534, 'fraction_numerical': 0.011589551573504503, 'mean_word_length': 4.281812831077104, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discovery and Visualization of Nonstationary Causal Models

It is commonplace to encounter nonstationary data, of which the underlying generating process may change over time or across domains. The nonstationarity presents both challenges and opportunities for causal discovery. In this paper we propose a principled framework to handle nonstationarity, and develop some methods to address three important questions. First, we propose an enhanced constraint-based method to detect variables whose local mechanisms are nonstationary and recover the skeleton of the causal structure over observed variables. Second, we present a way to determine some causal directions by taking advantage of information carried by changing distributions. Third, we develop a method for visualizing the nonstationarity of causal modules. Experimental results on various synthetic and real-world data sets are presented to demonstrate the efficacy of our methods.


Introduction
In many fields of empirical sciences and engineering, we would like to obtain causal knowledge for many purposes. As it is often difficult if not impossible to carry out randomized experiments, inferring causal relations from purely observational data, known as the task of causal discovery, has drawn much attention in several fields including computer science, statistics, philosophy, economics, and neuroscience. With the rapid accumulation of huge volumes of data of various types, causal discovery is facing exciting opportunities but also great challenges. One phenomenon such data often feature is that of distribution shift. Distribution shift may occur across domains or over time. For an example of the former kind, consider the problem of remote sensing image classification, which aims to derive land use and land cover information through the process of interpreting and classifying remote sensing imagery. The data collected in different areas and at different times usually have different distributions due to different physical factors related to ground, vegetation, illumination conditions, etc. As an example of the latter kind, the fMRI recordings are usually nonstationary: the causal connections in the brain may change with stimuli, tasks, states, the attention of the subject, etc. More specifically, it is believed that one of the basic properties of the neural connections in the brain is its time-dependence [1]. To these situations many existing approaches to causal discovery fail to apply, as they assume a fixed causal model and hence a fixed joint distribution underlying the observed data.
In this paper we assume that the underlying causal structure is a directed acyclic graph (DAG), but the mechanisms or parameters associated with the causal structure, or in other words the causal model, may change across domains or over time (we allow mechanisms to change in such a way that some causal links in the structure become vacuous or vanish over some time periods or domains). We aim to develop a principled framework to model such situations as well as practical methods to address these questions: • How to efficiently identify the variables whose local causal mechanisms are nonstationary and recover the skeleton of the causal structure over the observed variables?
• How to take advantage of the information carried by distribution shifts for the purpose of identifying causal directions?
• How to visualize the nonstationarity of those causal mechanisms that change over time or across domain?
This paper is organized as follows. In Section 2 we define and motivate the problem in more detail and review related work. Section 3 proposes an enhanced constraint-based method for recovering the skeleton of the causal structure over the observed variables and identify those variables whose generating processes are nonstationary. Section 4 develops a method for determining some causal directions by exploiting nonstationarity. Section 5 proposes a way to visualize nonstationarity. Section 6 reports simulations results to test the performance of the proposed causal discovery approach when the ground truth is known. Finally, we apply the method to some real-world datasets, including financial date and fMRI data, in Section 7.

Causal Discovery of Fixed Causal Models
Most causal discovery methods assume that there is a fixed causal model underlying the observed data and aim to estimate it from the data. Classic ap-proaches to causal discovery divide roughly into two types. In late 1980's and early 1990's, it was noted that under appropriate assumptions, one could recover a Markov equivalence class of the underlying causal structure based on conditional independence relationships among the variables [2,3]. This gives rise to the constraint-based approach to causal discovery, and the resulting equivalence class may contain multiple DAGs (or other related graphical objects to represent causal structures), which entail the same conditional independence relationships. The required assumptions include the causal Markov condition and the faithfulness assumption, which entail a correspondence between separation properties in the underlying causal structure and statistical independence properties in the data. The so-called score-based approach (see, e.g., [4,5]) searches for the equivalence class which gives the highest score under some scoring criterion, such as the Bayesian Information Criterion (BIC) or the posterior of the graph given the data.
Another set of approaches is based on restricted functional causal models, which represent the effect as a function of the direct causes together with an independent noise term [3]. Under appropriate assumptions, these approaches are able to identify the whole causal model. More specifically, the causal direction implied by the restricted functional causal model is generically identifiable, in that the model assumptions, such as the independence between the noise and cause, hold only for the true causal direction and are violated for the wrong direction. Examples of such restricted functional causal models include the Linear, Non-Gaussian, Acyclic Model (LiNGAM [6]), the additive noise model [7,8], and the post-nonlinear causal model [9]. The method presented in [10] makes use of a certain type of smoothness of the function in the correct causal direction to distinguish cause from effect, though it does not give explicit identifiability conditions.

With Nonstationary Causal Models
Suppose we are working with a set of observed variables V = {V i } n i=1 and the underlying causal structure over V is represented by a DAG G. For each V i , let P A i denote the set of parents of V i in G. Suppose at each point in time or in each domain, the joint probability distribution of V factorizes according to G: We call each P (V i | P A i ) a causal module. If there are distribution shifts (i.e., P (V) changes over time or across domains), at least some causal modules P (V k | P A k ), k ∈ N must change. We call those causal modules nonstationary causal modules. Their changes may be due to changes of the involved functional models, causal strengths, noise levels, etc. We assume that those quantities that change over time or cross domains can be written as functions of a time or domain index, and denote by C such an index. Figure 1: An illustration on how ignoring changes in the causal model may lead to spurious connections by the constraint-based method. (a) The true causal graph (including confounder g(C)). (b) The estimated conditional independence graph on the observed data in the asymptotic case.
If the changes in some modules are related, one can treat the situation as if there exists some unobserved quantity (confounder) which influences those modules and, as a consequence, the conditional independence relationships in the distribution-shifted data will be different from those implied by the true causal structure. Therefore, standard constraint-based algorithms such as PC [2,3] may not be able to reveal the true causal structure. As an illustration, suppose that the observed data were generated according to Fig. 1(a), where g(C), a function of C, is involved in the generating processes for both V 2 and V 4 ; the conditional independence graph for the observed data then contains spurious connections V 1 − V 4 and V 2 − V 4 , as shown in Fig. 1  Moreover, when one fits a fixed functional causal model (e.g., a linear, non-Gaussian model [6]) to distribution-shifted data, the estimated noise may not be independent from the cause any more. Consequently, the approach based on restricted functional causal models in general cannot infer the correct causal structure either. Fig. 2 gives an illustration of this point. Suppose we have two datasets for variables V 1 and V 2 : V 2 is generated from V 1 according to V 2 = 0.3V 1 + E in the first and according to V 2 = 0.7V 1 + E in the second, and in both datasets V 1 and E are mutually independent and follow a uniform distribution. Fig. 2(a -c) show the scatter plots of V 1 and V 2 on dataset 1, on dataset 2, and on merged data, respectively. (d) then shows the scatter plot of V 1 , the cause, and the estimated regression residual on both datasets; they are not independent any more, although on either dataset the regression residual is independent from V 1 .
To tackle the issue of changing causal models, one may try to find causal models on sliding windows [11] (for nonstationary data) or for different domains (for data from multiple domains) separately, and then compare them. Improved versions include the online changepoint detection method [12], the online undirected graph learning [13], the locally stationary structure tracker algorithm [14]. Such methods may suffer from high estimation variance due to sample scarcity, large type II errors, and a large number of statistical tests. Some methods aim to estimate the time-varying causal model by making use of certain types of smoothness of the change [15], but they do not explicitly locate the nonstationary causal modules. Several methods aim to model time-varying time-delayed causal relations [16,17], which can be reduced to online parameter learning because the direction of the causal relations is given (i.e., the past influences the future). Compared to them, learning changing instantaneous causal relations, with which we are concerned in this paper, is generally more difficult. Moreover, most of these methods assume linear causal models, limiting their applicability to complex problems with nonlinear causal relations.
In contrast, we will develop a nonparametric and computationally efficient method that can identify nonstationary causal modules and recover the causal skeleton. We will also show that distribution shifts actually contain useful information for the purpose of determining causal directions and develop practical algorithms accordingly.
3 Enhanced Constraint-Based Procedure

Assumptions
As already mentioned, we allow changes in some causal modules to be related, which may be explained by positing unobserved confounders. Intuitively, such confounders may refer to some high-level background variables. For instance, for fMRI data, they may be the subject's attention or unmeasured background stimuli impinging on a subject-scanner noise, random thoughts, physical sensations, etc.; for the stock market, they may be related to economic policies and changes in the ownership among the companies, etc. Thus we do not assume causal sufficiency for the set of observed variables. However, we assume that the confounders, if any, can be written as smooth functions of time or domain index. It follows that at each time or in each domain, the values of these confounders are fixed. We call this a weak causal sufficiency assumption.
Denote by {g l (C)} L l=1 the set of such confounders (which may be empty). We further assume that for each V i the local causal process for V i can be represented by the following structural equation model (SEM): where g i (C) ⊆ {g l (C)} L l=1 denotes the set of confounders that influence V i , θ i (C) denotes the effective parameters in the model that are also assumed to be functions of C, and i is a disturbance term that is independent of C and has a non-zero variance (i.e., the model is not deterministic). We also assume that the 's are mutually independent.
Note that {g l (C)} L l=1 are introduced to account for changes in different causal modules that are not independent. As a result, although θ i (C) may also contribute to a change in the causal module for V i , changes to the module for V i due to θ i (C) are independent of changes to the module for V j due to θ j (C), i = j. In other words, θ i (C) is specific to V i and is independent of θ j (C) for i = j. Note that g i (C) and θ i (C) can be constant, corresponding to stationary causal modules.
In this paper we treat C as a random variable, and so there is a joint distribution over V ∪ {g l (C)} L l=1 ∪ {θ m (C)} n m=1 . We assume that this distribution is Markov and faithful to the graph resulting from the following additions to G (which, recall, is the causal structure over V): to G, and for each i, add an arrow from each variable in g i (C) to V i and add an arrow from θ i (C) to V i . We will refer to this augmented graph as G aug .
Obviously G is simply the induced subgraph of G aug over V.

Detecting Changing Modules and Recovering Causal Skeleton
In this section we propose a method to detect variables whose modules change and infer the skeleton of G. The basic idea is simple: we use the (observed) variable C as a surrogate for the unobserved V ∪ {g l (C)} L l=1 , or in other words, we take C to capture C-specific information. 1 We now show that given the assumptions in 3.1, we can apply a constraint-based algorithm to V ∪ {C} to detect variables with changing modules and recover the skeleton of G.

(Detection of changing modules) For every i, test for the marginal and
conditional independence between V i and C. If they are independent given a subset of {V k | k = i}, remove the edge between V i and C in U C .
3. (Reovery of causal skeleton) For every i = j, test for the marginal and conditional independence between V i and V j . If they are independent given a subset of The procedure is briefly described in Algorithm 1. It outputs an undirected graph, U C , that contains C as well as V. In Step 2, whether a variable V i has a changing module is decided by whether V i and C are independent conditional on some subset of other variables. The justification for one side of this decision is trivial. If V i 's module does not change, that means P (V i | P A i ) remains the same for every value of C, and so V i ⊥ ⊥ C | P A i . Thus, if V i and C are not independent conditional on any subset of other variables, V i 's module changes with C, which is represented by an edge between V i and C. Conversely, we assume that if V i 's module changes, which entails that V i and C are not independent given P A i , then V i and C are not independent given any other subset of V\{V i }. If this assumption does not hold, then we only claim to detect some (but not necessarily all) variables with changing modules.
Step 3 aims to discover the skeleton of the causal structure over V. Its (asymptotic) correctness is justified by the following theorem: Theorem 1. Given the assumptions made in Section 3.1, for every V i , V j ∈ V, V i and V j are not adjacent in G if and only if they are independent conditional on some subset of Proof. Before getting to the main argument, let us establish some implications of the SEMs Eq. 2 and the assumptions in Section 3.1. Since the structure is assumed to be acyclic or recursive, according to Eq. 2, all variables V i can be written as a function of As a consequence, the probability distribution of V at each value of C is determined by the distribution of 1 , ..., n , and the values of {g l (C)} L l=1 ∪ {θ m (C)} n m=1 . In other words, p(V|C) is determined by n i=1 p( i ) (for 1 , ..., n are mutually independent), and {g l (C)} L l=1 ∪ {θ m (C)} n m=1 , where p(·) denotes the probability density or mass function. For any That is, By the weak union property of conditional independence, it follows that We are now ready to prove the theorem. Let V i , V j be any two variables in V. First, suppose that V i and V j are not adjacent in G. Then they are not adjacent in G aug , which recall is the graph that incorporates Because all g l (c) and θ m (C) are deterministic functions of C, . According to [18] , Eqs. 6 and 4 imply By the weak union property of conditional independence, it fol- In other words, V i and V j are conditionally independent given a subset of Conversely, suppose V i and V j are conditionally independent given a subset We show that V i and V j are not adjacent in G, or equivalently, that they are not adjacent in G aug . There are two possible cases to consider: • Suppose S does not contain C. Then since the joint distribution over is assumed to be Faithful to G aug , V i and V j are not adjacent in G aug , and hence not adjacent in G.
According to Eq. 3, and also noting that which also implies Substituting Eqs. 8 -10 into Eq. 7 gives That is, Again, by the Faithfulness assumption on G aug , this implies that V i and V j are not adjacent in G aug and hence are not adjacent in G.
Therefore, V i are V j are not adjacent in G if and only if they are conditionally independent given some subset of In the above procedure, it is crucial to use a general, nonparametric conditional independence test, for how variables depend on C is unkown and usually very nonlinear. In this work, we use the kernel-based conditional independence test (KCI-test [19]) to capture the dependence on C in a nonparametric way. By contrast, if we use, for example, tests of vanishing partial correlations, as is widely used in the neuroscience community, the proposed method will not work well.

An Advantage of Nonstationarity in Determination of Causal Direction
We now show that using the additional variable C as a surrogate not only allows us to infer the skeleton of the causal structure, but also facilitates the determination of some causal directions. Let us call those variables that are adjacent to C in the output of Algorithm 1 "C-specific variables", which are actually the effects of nonstationary causal modules. For each C-specific variable V k , it is possible to determine the direction of every edge incident to V k , or in other words, it is possible to infer P A k . Let V l be any variable adjacent to V k in the output of Algorithm 1. There are two possible cases to consider: 1. V l is not adjacent to C. Then C − V k − V l forms an unshielded triple in the skeleton. For practical purposes, we can take the direction between C and V k as C → V k (though we do not claim C to be a cause in any substantial sense). Then we can use the standard orientation rules for unshielded triples to orient the edge between V k and V l [2,3]: if V l and C are independent given a set of variables excluding V k , then the triple is a V-structure, and we have V k ← V l . Otherwise, if V l and C are independent given a set of variables including V k , then the triple is not a V-structure, and we have V k → V l .
2. V l is also adjacent to C. This case is more complex than Case 1, but it is still possible to identify the causal direction between V k and V l , based on the principle that P (cause) and P (effect | cause) change independently; a heuristic method is given in Section 4.1.
The procedure in Case 1 contains the methods proposed in [20,21] for causal discovery from changes as special cases, which may also be interpreted as special cases of the principle underlying the method for Case 2: if one of P (cause) and P (effect | cause) changes while the other remains invariant, they are clearly independent.

Inference of the Causal Direction between Variables with Changing Modules
We now develop a heuristic method to deal with Case 2 above. For simplicity, let us start with the two variable case: suppose V 1 and V 2 are adjacent and are both adjacent to C (and not adjacent to any other variable). We aim to identify the causal direction between them, which, without loss of generality, we suppose to be V 1 → V 2 . The guiding idea is that nonstationarity may carry information that confirms "independence" of causal modules, which, in the simple case we are considering, is the "independence" between P (V 1 ) and P (V 2 |V 1 ). If P (V 1 ) and P (V 2 |V 1 ) are "independent" but P (V 2 ) and P (V 1 |V 2 ) are not, then the causal direction is inferred to be from V 1 to V 2 . The idea that causal modules are "independent" is not new, but in a stationary situation where each module is fixed, such independence is very difficult, if not impossible, to test. By contrast, in the situation we are considering presently, both P (V 1 ) and P (V 2 |V 1 ) are nonstationary, and we can try to measure the extent to which variation in P (V 1 ) and variation in P (V 2 ) are dependent (and similarly for P (V 2 ) and P (V 1 |V 2 )). This is the sense in which nonstationarity actually helps in the inference of causal directions, and as far as we know, this is the first time that such an advantage is exploited in the case where both P (cause) and P (effect | cause) change.
We now derive a method along this line. Note that although both of V 1 and V 2 are adjacent to C, there does not necessarily exist a confounder. Fig. 3(a) shows the case where the involved changing parameters, θ 1 (C) and θ 2 (C) are independent, i.e., P (V 1; θ 1 ) and P (V 2 | V 1 ; θ 2 ) change independently. (We dropped the argument C in θ 1 and θ 2 to simplify notations.) For the reverse direction, one can decompose the joint distribution of (V 1 , V 2 ) according to where θ 1 and θ 2 are assumed to be sufficient for the corresponding distribution terms. Generally speaking, θ 1 and θ 2 are not independent, because they are determined jointly by both θ 1 and θ 2 . We assume that this is the case, and identify the direction between V 1 and V 2 based on this assumption. Now we face two problems. First, how can we compare the dependence between θ 1 and θ 2 and that between between θ 1 and θ 2 ? Second, in our nonparametric setting, we do not really have such parameters. How can we compare the dependence based on the given data?
For the first problem, we make use of the following measures of contributions from the parameters. The total contribution (in a way analogous to causal effect; see [22]) from θ 1 and θ 2 to (V 1 , V 2 ) can be measured with mutual information: where the second equality holds because of the chain rule, and the last one because the sufficiency of θ 1 for P involves the regular mutual information and conditional mutual information.
Since θ 1 and θ 2 are dependent, their individual contributions to (V 1 , V 2 ) are redundant. Below we calculate the individual contributions. The contribution from θ 2 to V 2 is S θ 2 →V2 = I(θ 2 ; V 2 ). The contribution from θ 1 to V 1 has been derived in [22]: , whereθ 1 is an independent copy of θ 1 . As a consequence, the redundancy in the contributions from θ 1 and θ 2 is ∆ V2→V1 is always non-negative because it is a Kullback-Leibler divergence. One can verify that if θ 1 ⊥ ⊥ θ 2 , which implies , leading to ∆ V2→V1 = 0. ∆ V2→V1 provides a way to measure the dependence between θ 1 and θ 2 . Regarding the second problem mentioned above, since we do not have parametric models, we propose to estimate ∆ V2→V1 from the data by: where · denotes the sample average,P (V 1 | V 2 ) is the empirical estimate of P (V 1 | V 2 ) on all data points, and P (V 1 | V 2 ) denotes the sample average of P (V 1 | V 2 ), the estimate of P (V 1 | V 2 ) at each time (or in each domain). In our implementation, we used kernel density estimation (KDE) on all data points to estimateP (V 1 | V 2 ), and used KDE on sliding windows (or in each domain) to estimateP (V 1 | V 2 ). We take the direction for which∆ is smaller to be the causal direction. If there is a confounder g 1 (C) underlying V 1 and V 2 , as shown in Fig. 3(b), we conjecture that the above approach still works if the influences from g 1 (C) is not very strong, for the following reason: for the correct direction,∆ measures the influence from the confounder; for the wrong direction, it measures the influence from the confounder and the dependence in the "parameters" caused by the wrong causal direction. A future line of research is to seek a more rigorous theoretical justification of this method. When there are more than two variables which are connected to C and inter-connected, we try all possible causal structures and choose the one that minimizes the total∆ value, i.e., i:P A i =∅∆ P A i →Vi .

Kernel Nonstationarity Visualization of Causal Modules
It is informative to determine for which variable the causal model (data-generating process), or P (V i | P A i ), changes. But usually it is not enough -one often wants to interpret the pattern of the changes, find what causes the changes, and understand the causal process in more detail. To achieve so, it is necessary to discover how the causal model changes, i.e., where the changes occur and how fast it changes, and visualize the changes. Although the changes occur in the conditional distribution P (V i | P A i ), usually it is not straightforward to see the properties of the changes by directly looking at the distribution itself. A lowdimensional representation of the changes is needed.
In the parametric case, if we know which parameters of the causal model P A i → V i are changing, which could be the mean of a root cause, the coefficients in a linear SEM, etc., then we can estimate such parameters for different values of C and see how they change. However, such knowledge is usually not available, and more importantly, for the sake of flexibility we often model the causal processes nonparametrically. Therefore, it is desirable to develop a general nonparametric procedure for nonstationarity visualization of changing causal modules.
Note that changes in P (V i | P A i ) are irrelevant to changes in P (P A i ), and accordingly, they are not necessarily the same as changes in the joint distribution P (V i , P A i ). (If V i is a root cause, P A i is an empty set, and P (V i | P A i ) reduces to the marginal distribution P (V i ).) We aim to find a mapping of P (V i | P A i ) which captures its nonstationarity: We call λ i (C) the nonstationarity encapsulator for P (V i | P A i , C). This formulation is rather general: any identifiable parameters in P (V i | P A i , C) can be expressed this way, and in the nonparametric case, λ i (C) can be seen as a statistic to summarize changes in P (V i | P A i , C) along with C. If P (V i | P A i , C) does not change along with C, then λ i (C) remains constant. Otherwise, λ i (C) is intended to capture the variability of P (V i | P A i , C) across different values of C. Now there are two problems to solve. One is given only observed data, not the conditional distribution, how to represent λ i (C) in Eq. 15 conveniently.
The other is what criterion and method to use to enable λ i (C) to capture the variability in the conditional distribution along with C. We tackle the above two problems by making use of kernels [23], and accordingly propose a method called kernel nonstationarity visualization (KNV) of causal modules.

Using Kernel Embedding of Conditional Probabilities
We use the kernel embedding of conditional distributions [24] instead of the original conditional distributions. Suppose we have kernels k XX , where C Y X and C XX denote the (uncentered) crosscovariance and covariance operators, respectively [25]. The empirical estimate of U Y |X isÛ Y |X = Ψ Y (K X + βI) −1 Ψ X , where β is a regularization parameter (set to 0.05 in our experiments), and Ψ Y , Ψ X , and K X are the feature matrix on Y , feature matrix on X, and the kernel matrix on X, respectively [24]. We use the Gaussian kernel for k (1) X and k (1) Y with kernel width σ 1 , andÛ Y |X encodes the information of P (Y | X) on the given data.
In our problem, we need consider the kernel conditional distribution embedding of P (V i | P A i ) for each value of C. If C is a domain index, for each value of C we have a dataset of (V i | P A i ). If C is a time index, we use a sliding window to find the data corresponding to C = c, by using the data of (V i , P A i ) in the window of length L centered at c. As we shall see later, It is possible to avoid directly calculating the empirical estimate of the embedding, but we need the following ("cross") kernel (or Gram) matrices: K Vi (c, c ) is the "cross" kernel matrix between the values of V i corresponding to C = c and those corresponding to C = c , and similarly for K P A i (c, c ).

Nonstationary Encapsulator Extraction by Eigenvalue Decomposition
Next, in principle, we use the estimated kernel embedding of conditional distributions,Û Vi | P A i ,C=c , as input, and aim to findλ i (c) as a (nonlinear) mapping ofÛ Vi | P A i ,C=c , to capture its variability across different c. This can be readily achieved by exploiting some nonlinear principle component analysis (PCA) techniques, and here we adopted kernel principal component analysis problem (KPCA) [26], for its nice formulation and computational efficiency. KPCA computes principal components in high-dimensional feature spaces of the input. In our case, for each c the input,Û Vi | P A i ,C=c , is a matrix. We can stack it into a long vector, and then representλ i (c) by making use of a second kernel, k (2) (which is usually different from k (1) ), as required by KPCA. Denote by the corresponding Gram matrix by M , whose (c, c )th entry is, the empirical kernel maps of V i and P A i ; below we show that we can directly find M without explicitly making use of empirical kernel maps. If we use a linear kernel for k (2) , the (c, c )th entry of M is 2 If k (2) is a Gaussian kernel with kernel width σ 2 , we have where || · || F denotes the Frobenius norm. Finally,λ i (C) can be found by performing eigenvalue decomposition on the above Gram matrix, M l or M g ; for details please see [26]. Algorithm 2 summarizes the proposed KNR method. There are several hyperparameters to set. In our experiments, we set the kernel width σ 2 1 (for k (1) ) and σ 2 2 (for k (2) ) to the median distance between points in the sample, as in [27]. We kept the window length L = 100. k (1) . If C is a time index, they can be obtained by extracting corresponding entries of the kernel matrices K Vi and K P A i on the whole data.
2. Calculate Gram matrix M with kernel k (2) (see Eq. 16 for linear kernels and Eq. 17 for Gaussian kernels).
3. Findλ i (C) by directly feeding Gram matrix M to KPCA. That is, perform eigenvalue decomposition on M to find the nonlinear principal componentŝ λ i (C), as in Section 4.1 of [26]. 2 When P A i is an empty set, P (V i | P A i ) reduces to P (V i ). In this case we use the em- , whose empirical estimate is the sample mean of ψ(V i ) on the sample. Here ψ(·) denotes the feature map. Accordingly, M l (c, c ) reduces to 1 ncn c 1 n c K V i (c , c)1n c , where nc and n c are the sizes of the data corresponding to C = c and C = c , respectively, and 1n c is the vector of 1's of length nc.

A Toy Example
We generated synthetic data according to the SEMs specified in Fig. 4. More specifically, the exogenous input to V 1 , the causal strength from V 3 to V 5 (the coefficient f 3 in the structural equation for V 5 ), and the noise variance in the equation for V 4 are time varying; the changing parameters were represented by sinusoid or cosine functions of T . We used different periodic levels (w = 5, 10, 20, 30) of the varying components, as well as different sample sizes (N = 600, 1000). In each setting, we run 10 replications, with both our enhanced constraint-based method (Algorithm 1, with the time index for C) and the original constraint-based method; we used the SGS search procedure [28] and kernel-based conditional independence test [19].
with i = 1, · · · , 6.    of the discovered adjacencies between the V 's at significance level 0.05. It is obvious that compared to the original method, our method effectively reduces the number of spurious connections, i.e., edges (V 1 , V 4 ), (V 1 , V 5 ) and (V 4 , V 5 ), in all the settings. The FN rate only very slightly increases. As w increases, the FP rate stays stable, and the FN rate slightly increases for both methods; as N increases, the FN rate is greatly reduced. In addition, from the augmented causal graph, we can identify causal directions by the procedure in Section 4. In this simulation, the whole causal DAG is correctly identified. However, with the original SGS method, we can only identify two causal directions: 5 → 6 and 2 → 6, and there are spurious edges (V 1 , V 4 ), (V 1 , V 5 ) and (V 4 , V 5 ). Furthermore, we visualized the nonstationarity of causal modules, P (V 1 ), (V 2, V 3) → V 4 , and V 3 → V 5 , by KNV (Algorithm 2). We tried both the linear kernel and Gaussian kernel for k (2) . Figure 6 shows the first component of the extracted nonstationarity encapsulatorsλ i , i = 1, 4, 5, corresponding to the three nonstationary causal models; see the blue solid lines. Panels (a) and (b) correspond to the setting w = 5, N = 600 and w = 30, N = 600, respectively. The red dashed lines show the changing parameters f 1 , f 2 , and f 3 in the respective causal models. Note that they have been rescaled to match with the nonstationarity encapsulatorsλ i . We can see that KNV successfully recovers the variability in the causal models (as represented by changing parameters f 1 , f 2 , f 3 , corresponding to changes in the causal strength or noise variance). In addition, the Gaussian kernel gives better results especially in the case where w = 30. To summarize, we found that when there is only one changing parameter in the causal model P (V i | P A i ), which may be the linear coefficient, the mean of the noise, or its variance, with the Gaussian kernel for k (2) , one component of λ i (C) is usually enough to capture the changes -this component is close to a nonlinear transformation of the changing parameter, and its corresponding eigenvalue is at least five times bigger than the remaining ones. However, if the functional form of the causal model changes, say, if the SEM changes from a linear one to a quadratic one, more than one component of λ i (C) has relatively large eigenvalues, and they jointly capture the change in P (V i | P A i ) (results are not included here).

Experimental Results on Simulated fMRI
In recent years the brain effective connectivity study from fMRI has received much attention. The fMRI experiments may last for a relatively long time period, during which the causal influences are likely to change along with certain unmeasured states (e.g., the attention) of the subject and ignoring the timedependence may lead to spurious connections. Likewise, the causal influences may also vary as a function of the experimental condition (e.g., health, disease, and behavior) [11].
Currently little is known for the causal connectivity in our brain, so firstly we applied our approach on simulated fMRI data which enables us to evaluate the robustness of our approach with known ground truth. We generated the simulated fMRI signal according to the DCM forward model [29]. Fig. 7 shows a basic setting of the network topologies, where we modeled the external input u 1 to the nodes as random square wave [30], and the external input to the connections with different kinds of functions, e.g., exponential decay, square wave, and log functions. Since the study on how causal connections between brain regions are changing is very limited, we tried to represent them with different functions to model different possible scenarios. In addition, in practice we may analyze the fMRI signal concatenated from different scans (different subjects or different instruments), so in order to model this situation, we concatenated two generated BOLD signals to derive the final signal. Figure 7: The basic setting of the network topologies.
We tested our enhanced constraint-based method on 50 realizations, where the time information T is included into the system to capture smooth varying causal relations and the influences from smooth varying latent confounders. Fig.  8(a) gives the False Positive (FP) rate and False Negative (FN) rate at significance level 0.03. We compared our enhanced constraint-based method with the original one (both with SGS search and KCI test), and we also compared with partial correlation test since it is widely used in fMRI analysis [30]. It is obvious that our approach greatly reduces the FP rate, that is, it effectively reduces spurious connections which are induced by the time-varying connections, while at the same time increases the FN rate in a reasonable range. The partial correlation test gives the worst results, with the FP rate 1 and the FN rate 0.1016 in a small-sample-size case. Since there is a certain amount of variation across realizations, we give a causal connection if it exists in more than 80% of all the realizations. Fig. 8(b-c) show the causal structures estimated by our approach and the original constraint-based method with KCI-test. The partial correlation test produces fully connected graph. (a) Estimation error.
(c) Estimated graph by the original method Figure 8: (a) The estimation error, FP rate and FN rate, derived from our enhanced constraint-based method, the original constraint-based method with KCI and the partial correlation test. (b,c) The estimated causal graph by the our approach and the original constraint-based one. Seng Utilities Sub-index (HSU). It is believed that during the financial crisis around 2008, the causal relations in Hong Kong stock market have changed. Fig. 9 shows the estimated causal structure by our method, where 2, 3, 4, 5, 6 and 7 are found to be time-dependent as indicated by red cycles. In contrast, the original constraint-based method produces four more edges, which are (2,3), (3,6), (5,7) and (6,8). We found that all time-dependent returns are in HSF, HSP, and HSU, which are directly affected by some unconsidered factors, e.g. policy changes. Furthermore, we inferred the causal directions by the procedure given in Section 4, and we found that all the inferred directions are reasonable. In particular, the within sub-index causal directions tend to satisfy the ownermember relationship. For example, 4 → 1 because 1 partially owns 4, and similarly for 5 → 3 and 9 → 8. Those stocks in HSF are the major causes to those in HSC and in HSP, and the stocks in HSP and HSU impact those in HSC. These causal relations match with the fact that financial institutions are in the leading position to impact other fields, and industries are usually affected by financial institutions, companies in properties, and companies in utilities. One exception is that, 10, Bank of China Hong Kong in HSF, is affected by 2 in HSP; it is perhaps because of Bank of China Hong Kong's close relation with Bank of China in mainland China. Figure 10 (bottom panels) visualizes the nonstationarity of (changing) causal modules, for stocks 2, 3, 4, 5, 6, and 7. We can see that the nonstationary en- The nonstationarity encapsulators, especially those of 2, 4, 5, and 3, seem to be consistent with the change of the TED spread, 4 which is an indicator of perceived credit risk in the general economy and shown in Figure 10 (top panel) for comparison; 7 and 6 seem to be directly influenced by the change in the underlying unmeasured factor, which may be related to the credit risk, mainly from 2008.

On fMRI Hippocampus
This fMRI Hippocampus dataset [31] was recorded from six separate brain regions: perirhinal cortex (PRC), parahippocampal cortex (PHC), entorhinal cortex (ERC), subiculum (Sub), CA1, and CA3/Dentate Gyrus (CA3) in the resting states on the same person in 64 successive days. We are interested in investigating causal connections between these six regions in the resting states. The anatomical connections between them reported in the literature are shown in Fig. 11. We used the anatomical connections as a reference, because in theory a direct causal connection between two areas should not exist if there is no anatomical connection between them.
PHC ERC CA3 CA1 PRC Sub Figure 11: The anatomical connections between the six separate brain regions.
We applied our enhanced constraint-based method on 10 successive days separately, with time information T as an additional variable in the system. We assumed that the underlying causal graph is acyclic, although the anatomical structure gives cycles. We found that our method effectively reduces the FP rate, from 62.86% to 17.14%, compared to the original constraint-based method with SGS search and KCI-test. Here we regard those connections that do not exist in the anatomical structure as spurious; however, with the lack of ground truth, we are not able to compare the FN rate. We found that the causal structure varies across days, but the connections between CA1 and CA3, and between CA1 and SUB are robust, which coincides with the current findings in neuroscience [32]. In addition, on most datasets the causal graphs we derived are acyclic, which validates the use of constraint-based method. Furthermore, we applied the procedure in Section 4 to infer the causal directions. We successfully recovered the following causal directions: CA3 → CA1, CA1 → Sub, Sub → ERC, ERC → CA1 and PRC → ERC, and the accuracy of direction determination is 85.71% (we consider the anatomical connections, shown in Fig. 11, as ground truth for the directions).

On WiFi Dataset
The WiFi dataset has been seen as a benchmark dataset to test the performance of domain adaptation algorithms. The indoor WiFi localization data can be easily outdated since the WiFi signal strength (features) may vary with time periods, devices, space and usage of the WiFi [33]. Therefore, it is important to detect the domain-varying features in domain adaptation. In this dataset, the data were collected from three different time periods in the same locations.
We added the domain information (D = 1, 2, 3) as an additional variable in the causal system to capture the domain-varying features. Here we set the significance level as 0.05 and we found that only a small subsets of features (8/67) vary across domains, with feature index = {1, 2, 3, 4, 5, 6, 12, 44}, which provides benefits for further analysis in domain adaptation. We also found that compared to the original constraint-based method, our method gives much sparser connections between the features (the number of connections between the features is reduced from 52 to 26).

Conclusion
This paper is concerned with discovery and visualization of nonstationary models, where causal modules may change over time or across datasets. We assume a weak causal sufficiency condition, which states that all confounders can be written as smooth functions of time or the domain index. We proposed (1) an enhanced constraint-based method for locating variables whose causal modules are nonstationary and estimating the skeleton of the causal structure over the observed variables, (2) a method for causal direction determination that takes advantage of the nonstationarity, and (3) a technique for visualizing nonstationary causal modules.
In this paper we only considered instantaneous or contemporaneous causal relations, as indicated by the assumption that the observed data are independently but not identically distributed; the strength (or model, or even existence) of the causal relations is allowed to change over time. We did not explicitly consider time-delayed causal relations and in particular did not engage autoregressive models. However, we note that it is natural to generalize our framework to incorporate time-delayed causal relations, just in the way that constraint-based causal discovery was adapted to handle time-series data (see, e.g., [34]).
There are several open questions we aim to answer in future work. First, in this paper we assumed that causal directions do not flip despite of nonstationarity. But what if some causal directions also change over time or across domains? Can we develop a general approach to detect causal direction changes? Second, to fully determine the causal structure, one might need to combine the proposed framework with other approaches, such as those based on restricted functional causal models. How can this be efficiently accomplished? Third, the issue of distribution shift may decrease the power of statistical (conditional) independence tests. Is it possible to mitigate this problem?