Causal Discovery in Semi-Stationary Time Series
Authors: Shanyun Gao
Raghavendra Addanki
Tong Yu
Ryan A. Rossi
Murat Kocaoglu
Venue: NeurIPS 2023 poster (2023)
Keywords: time-series causal discovery
constraint-based causal discovery
TL;DR: We propose a constraint-based, non-parametric approach to discover the causal relations for a specific type of time series without assuming stationarity.

Abstract: Discovering causal relations from observational time series without making the stationary assumption is a significant challenge. In practice, this challenge is common in many areas, such as retail sales, transportation systems, and medical science. Here, we consider this problem for a class of non-stationary time series. The structural causal model (SCM) of this type of time series, called the semi-stationary time series, exhibits that a finite number of different causal mechanisms occur sequentially and periodically across time. This model holds considerable practical utility because it can represent periodicity, including common occurrences such as seasonality and diurnal variation. We propose a constraint-based, non-parametric algorithm for discovering causal relations in this setting. The resulting algorithm, PCMCI$_{\Omega}$, can capture the alternating and recurring changes in the causal mechanisms and then identify the underlying causal graph with conditional independence (CI) tests. We show that this algorithm is sound in identifying causal relations on discrete time series. We validate the algorithm with extensive experiments on continuous and discrete simulated data. We also apply our algorithm to a real-world climate dataset.

Causal Discovery in Semi-Stationary Time Series
Shanyun Gao Purdue University gao565@purdue.. edu
Raghavendra Addanki Adobe Research raddenk1@adobe. com com
Tong Yu Adobe Research tyu@adobe.o. com
Ryan A. Rossi Adobe Research
Syrossi@adobe.. com
Murat Kocaoglu
Purdue University
mkocaoglu@purdue.. edu
Abstract
Discovering causal relations from observational time series without making the stationary assumption is significant challenge. In practice, this challenge is common many areas, such as retail sales, transportation systems, and medical science. Here, consider this problem for class of non-stationary time series. The structural causal model (SCM) this type of time series, called the semi- stationary time series, exhibits that a finite number of different causal mechanisms occur sequentially and periodically across time. This model holds considerable practical utility because such as seasonality and diurnal variation. We propose a constraint- parametric algorithm for discovering causal relations in this setting. The resulting algorithm, PCMCI2, can capture the alternating and recurring changes in the causal mechanisms and then identify the underlying causal graph with conditional independence (CI) tests. We show that this algorithm is sound in identifying causal relations on discrete-valued time series. We validate the algorithm with extensive experiments on continuous and discrete simulated data. We also apply our algorithm real-world climate dataset.
represent periodicity, including common occurrences
non-
Introduction
In modern sciences, causal discovery aims to identify the collection of causal relations from obser-
vational data, in Pearl (1980, Peters al. [2017] and (Spirtes et al. (2000) One of the most popular causal discovery approaches is the so-called constraint-based method. Constraint-based approaches assume that the probability distribution of variables causal Markov and faithful
directed acyclic graph called the causal graph. Given large enough data, they can then recover the corresponding Markov equivalence class by exploiting conditional independence relationships variables. See Peters et al. [2017]. There are many constraint-bas
the algorithms such as PC and FCI algorithms Spirtes al. (2000) The standard assumption of these approaches is that data samples are independent and identically distributed, which makes it possible to perform CI tests. See Bergsma
2004 Zhang et al. 2012 and Shah and Peters 2020
Recently, there have been numerous efforts to extend such constraint-ban -based algorithms to accommo-
date time series data. For instance, PCMCI in Ruugee al. [2019] and LPCMCI in Gerhardus and Runge [2020] are the PC-based algorithms for time series. Inspired by FCI algorithms, approaches designed for time series include ANLSTM in Chu and Glymour 2008 tsFCI Entner and Hoyer [2010] and SVAR-FCI in Malinsky and Spirtes] (2018). This setup is relevant in several industrial applications since many data points have an associated time-point, such root-cause analysis in Ikram et al. 2022 Most of the existing causal discovery algorithms make the stationary assumption.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

See Chu and Glymour [2008]. Hyvärinen et al. 2010 Entner and Hoyer 2010 Peters et al. 2013
Mahinsky and Spiries 20188 Runge et al. [0019, Pamhl al. 2020 and Assaad et al. 2022
Non-stationary temporal data makes causal discovery more challenging since the statistics are time-variant, and it is unreasonable expect that the underlying causal structure is time-invariant. Identifying causal relations from non-stationary time series without imposing any restriction on the data is difficult. Here, we focus on specific class of non-stationary time series, called the semi- stationary time series, whose structural causal model (SCM) exhibits that a finite number of different causal mechanisms occur sequentially and periodically across time. One example is illustrated in 1 where the time series X has three different causal mechanisms across time, shown as red edges, green edges, and blue edges, respectively. Similarly, time series X has two alternative causal mechanisms. This setting holds considerable practical utility. Periodic nature is commonly observed in many real-world time series data. See Han al. [2002], Nakamura et al. 2003). Carskadon et al. [2055] and Komarzynski et al. [2018]. Here are a few additional intuitive examples: poor traffic conditions often coincide with commute time and weekends; household electricity consumption
typically follows a pattern of being higher at night and lower during the daytime. Consequently, it is reasonable to expect periodic changes in the causal relations underlying this type of time series without assuming stationarity. Here, the constraint-based methods in Chu and Glymour (2008). Entner and Hoyer 2010]] Malinsky and Spirtes (018 and Runge et al. (2019), designed for stationary time series, may fail. Given observational data with periodically changing causal structures, it is hard to apply CI tests directly. Most of the other algorithms designed for non-stationary time series rely
heavily on model assumptions, as in Gong et al. [2015]. Pamfil et al. [2020] and Huang et al. 22199] These algorithms are discussed further in the related work.
In this paper, we propose an algorithm to address this problem, namely non- non-parametric causal discovery in time series data with semi-stationary SCMs. The key contributions of our work are:
We develop algorithm discover the causal structure from semi-stationary time series data where the underlying causal structures change periodically. Our algorithm systemati- cally uses the PCMCI algorithm proposed for the stationary setting in Runge et al. (2019] The resulting algorithm is hence named PCMCIa where Q denotes periodicity.
We validate our method with synthetic simulations on both continuous- valued and discrete- valued time series, showing that our method can correctly learn the periodicity and causal mechanism of the synthetic time series.
We utilize our method in real-world climate application. The result reveals the potential existence of periodicity in those time series, and the stationary assumption made by previous works could be relaxed in some practical situations.
1.1 Related Work
PCMCI has been applied in diverse domains to investigate atmospheric interactions in the biosphere,
as demonstrated in Krich et al. [2020] global wildfires as explored in Qu et al.] [2021], J. water usage
as studied in Zou et al. [022]. ultra-processed food manufacturing as examined in Menegozzo
le al. [2020] and causal feature selection as discussed in Peterson [2022]. among other applications.
See Arvind et al. 2021]. Gerhardus and Runge (2020). Castri et al. 2023a bj. While PCMCI has
achieved considerable success, it not without its limitations. One notable assumption that can be challenged is the concept of causal slationary, that is, causal relations are time-invariant. PCMCI
exhibits robustness when applied to linear models with an added non-stationary trend. See also Runge
et al. [019]. However, there is an ongoing exploration to enhance its performance in wider range of non-stationary settings.
Although not as extensively as the stationary case, causal discovery in non-stationary time series has been studied by some authors. However, many of those algorithms rely on parametric assumptions such as the vector autoregressive model in Gong et al. (2015) and Malinsky and Spirtes (2019) linear
and nonlinear state-space model in [Huang et al. [2019] One non- ponammraetttic algorithm the literature is CD-NOD proposed by |Huang et al. 2020). which has been extended to recover time-varying instantaneous and lagged causal relationships. In very recent work, Fujiwara et al. 20033 proposed an algorithm JTT-LiNGAM to obtain a local approximated linear causal model combining algorithm
LiNGAM and JIT framework for non-linear and non-stationary data. To the best of our knowledge,
no other non non-parametric approaches can discover causal relations underlying time series without

---

21
x7
X
X
C:
11 12
13
14
15
with Semi-Stationary SCM where Tmax = 3, W1 = 3, w22 2, W3 1, Ω = and 8 = 6. The first 3(=Tmax) time slices color edges represent the same causal mechanism. E.g. for X there are 3 (= wi) time partition subsets [ Jisggs. The time points & of nodes X; sharing the same filling color are in the same time partition subsets. The time points of nodes X; sharing both the same filling color and the same outline shape are in the same homogenous time partition subsets (the definitions are the supplementary material). There are (= 5) different Markov chains in this multivariate time series V. and the first element of these 6 Markov chains is shown as {44) sssse and are tinted with a gradient of blue hues. The superscript of 2 is the index of different Markov chains, whereas the subscript denotes the running index of that specific Markov chain. For instance, i and Z3 denote the first two elements of the first Markov chain, while 23 and Z denote the first elements of the second Markov chain.
Figure 1: Partial causal graph for 3-variate time series V =
X2 X (X
[X:) Sisa are the starting points. The
assuming stationarity and can also allow for sudden changes in causal mechanisms, Our proposed approach does not directly enforce the stationary assumption the time series. The SCM also integrates finite set of causal mechanisms that exhibit periodic variations over time.
2 PCMCI:: Capturing Periodicity of the Causal Structure
In this section, formulate the problem of learning the causal graph multi-variat time series data when the SCM exhibits periodicity in causal mechanisms. In section 2. we present the necessary
definitions and provide an overview of the problem setting. In section 2.2 we introduce the required assumption.
2.1 Preliminaries
Let G(V, F.) denote the underlying causal graph, and for each variable X E V, we denote the set of all incoming neighbors as parents, denoted by Pa(X).
For any two variables X, V and C V, we denote the CI relation X is independent Y conditioned on S, by X IL Y S.
For simplicity's sake, define sets: [b] = {1, 21.,..... b} and [a, := {a, 1,..., b), where a, N.
In the time series setting, let X R² denote the variable of jth time series time t, X (X/hez IRT denote a univariate time series and X, = [Xx])en E R" denote a slice of variables at time point t. V = (X))ee = (Xthery E == denotes a n-variate time series. By default, we assume n > 1 and hence X Ç V, and p(V) o 0, where p(.) denotes the probability or probability density.
Definition 2.1 (Non-Stationany SCM), Stationai Structural Causal Model (SCM) is tuple M (V, F, EP where there exists Tmax N+ defined
Tmax arg max Xi Pa(X)), E [m]].
such that with this Tmax each variable X E V is deterministic function of its parent set Pa(Xt>rmal E V and an unobserved (exogenous) variable 2/> E:
x² = (Pa(()))) d), E [n], t [Tmax 1, 1,7].
3

---

and there exist at least two different time points to, t1 E [Tmax + 1] satisfying
S1 + Jan 3je E [n), 41 C [Tmax 1, 1.
where fj Jo Sys, F and {E) heery are jointly independent with probability measure P. Tmax is the finite maximal lag in terms of the causal graph G.
Definition 2.2 Remi-Stationary SCM). A
SCM is Non-Stationary SCM that additionally satisfies the following conditions. For each j E [n], there exists an w E such that:
a) fi fia+NWw
(4)
b) Pa(X++NN) (X++N=: X E E [n]},
(5)
c) c e NN are i.i.d.
(6)
are satisfied for all [Tmax + 1,77,, number of causal mechanisms are repeated periodically for every univariate time series X in V. One example of this model illustrated in Fig I For X in Fig three causal mechanisms are reiterated periodically with w₁ = 3, represented by red, green, and blue edges, respectively.
{N N N, S T). This means that a finite
The minimum value that satísfies the above conditions for X is defined as the periodicity of X', denoted by wj. Furthermore, for an n-variate time series V, Ω2 denotes the minimum periodicity across all time series X j E [n.] The number of causal mechanisms occurring sequentially and periodically of univariate time series X is wj. and that number of causal mechanisms of multivariate time series Q. For X the causal mechanisms are associated with each variable X However, for V, the causal mechanisms are related to each time slice vector X, as a whole. The relationship between Q and W can be captured by:
f2 = LC=LCM(14 : [n]})
where LCM(.) is an operation to find the least common multiple between any two or more numbers. Here, Z is the smallest common multiple among {wi,... In Fig Ω LCM 2, 1) 6.
Definition 2.3 (Time Partition). A time partition IT) of univariate time series X in Semi- Stationary SCM with periodicity W is way of dividing all time points E [2'] into collection of non-overlapping non- non-empty subsets {IR& (T))k((ss) such that:
(T) == {t: Tmax + 1 S < T, (t mod wj) + =k.
kj.
(8)
where mod denotes the modulo operation. For instance, 5 mod 3 2
We can observe that the variables in hent share the same causal mechanism. Since the number of potentially different causal mechanisms of variables in X is Wj. the number of such time partition subsets is wy. For simplicity, notations II' and II are used instead of IT'(T) and TH, In Fig
I 7, 7, 10, 13, 4 + -.N,...,, IT} (5, 8, 11, 14, 3.....) and TTg {6,9. 12, 15, 6 1++33 .......> where N N The nodes X are classified into their associated time partition subsets by the matching colors.
Definition 2.4 (Illusory Parent Sets). For univariate time series X V with Semi-Stationary
SCM having periodicity wj 1, parent set index plndkele,) is defined as:
pInd = (1i, he(w") given Pa(X)) (X Xi X J. AI E II
(9)
where n' = (aa()))) Ta is the time lag and ys is the variable index. Given pInd Illusory Parent Sets are defined as:
Pau (x)( (x (vi, Ti) E plndj), Ak E (k: a IT))
(10)
Put simply, the illusory parent sets of X/ are the time-shift version of the parent set of other variables in X that have different causal mechanism from X Note that the illusory parent sets constructed specifically in the Semi-Stationary SCM. For stationary SCM, there is no illusory parent set needed. To maintain consistency in notation, for time points E IT/ the notation Pay (X)) can also be extended to encompass the true parent set of X?:
Pat(X)) Pa(X() Vt E III
(11)

---

By doing so, Pa(X)) c Use[u,) (x)( In Fig I by observing Pa(X;), we have pind =
11, (2, 2)); by observing Pa(X)), we have pind {(1,1), (3, 1)} and finally by observ-
ing Pa(Xg), plnd {(L, 1), (1, 2)}. Based on those indexes, Pa (X²) {
(X
{Xg, X3), (a2(X))
{ {Xg. ). The first one is the true parent set of X and the latter two are the illusory parent
{ X and Pas (aa(X))
sets. The order of those parent sets is not important.
At last, we need to further define a series of Markov chains that are associated tightly with Semi- Stationary SCM. The presence and characteristics of these Markov chains are thoroughly examined in the supplementary materials. The motivation behind creating such Markov chains is to introduce assumptions on them rather than directly on the original data V.
Definition 2.5 (Time Series as a Markov Chain). For time series V with Semi-Stationary SCM, there (potentially) 8 different Markov chains { MENN E [8]:
2% = (X Tmax Tmax +-+
1)6 X Tmax +-xxx++++++++ +n
Tmox tmmx)))
where N == {n E N+; Tmax +9 19-1- 1 + nd < 7), 5 = [magt1] 12. Note that in { {Z is used to indicate a specific Markov chain, while serves as the running index for that particular Markov chain. Such Markov chain { / [o] exists as long as Pa(Zğ) C 7g U 7 for all n. This is finite state Markov Chain if all time series in V are discrete-valued time series. The state space of {2A} the set containing all possible realization of +g+ (( hel.))n The transition probabilities between the states are the product of associated causal mechanisms based on Assumption A2, which elaborated by an example section C.1 (1q.(10)-(11)) of the supplementary material.
2.2 Assumptions for PCMCI2
A1. Sufficiency: There are no unobserved confounders.
A2. Causal Markov Condition: Each variable X is independent of all non-descendants, given its
parents Pa(X) in G.
A3. Faithfulness Condition (Pearl (9800 ): Let P be a probability distribution generated by G. (G, P) satisfies the Faithfulness Condition if and only if every conditional independence relation true in P is entailed by the Causal Markov Condition applied to G.
A4. No Contemporaneous Causal Effects: Edges between variables at the same time are not allowed.
A5. Temporal Priority: Causal relations that point from the future the past are not permitted.
A6. Hard Mechanism Change: If at time points t and ig, the causal mechanisms X and x are different, then their corresponding parent sets not be transformed to each other by time shifts:
fan, * f3.02 = Pa(X),) # {X] (x2-(ts-1)) X E Pa(X))), [n]].
A7. Irreducible and Aperiodic Markov Chain: The Markov chains {Zn} of V are assumed to be irreducible (Serfozo 2009 ): for all states i and j of [Zn J. In so that
Pij = p(Zn+1 = j = i) >
(12)
and aperiodic(Karlin 2014 ): for every state (Zn), d(i) = where the period d(i) of the state i is the greatest common divisor of all integers for which pi (x) 0.
Assumptions A1-A5 are conventional and commonly employed in causal discovery methods for time series data. On the other hand, our approach requires additional Assumptions A6-A7 to be in place. To clarify, A6 is essential because our method may encounter challenges in distinguishing distinct causal mechanisms for variables (X% (xtr if they share identical parent sets after time shifts. As for A7, it serves crucial role in establishing the soundness of our algorithm.
3 PCMCI, Algorithm
In this section, we propose an algorithm called PCMCIa, and in section 3.1 we present a theorem demonstrating the soundness of PCMCIa and its ability to recover the causal graph. Our algorithm
5

---

Algorithm 1 PCMCI
1: Input: A n-variate time series = ,..,,"" X"), periodicity upper bound wibs time
lag upper bound Tub By default, we assume Tuh and wab are larger than their true value.
2: A superset of parent set is obtained using PCMCI with Tub and denote by SPa(X?) Vj,t.
3: for X where j E [n] do
4: for guess w E Wub]] of wy do
5: Let .t == - lk: E [w]] where nT.. = (2Tub + k; 2T + w + k. 2Tub + 2w + ...)
6: for l E [w] do
Initialize the parent set for X2, ( E {t 2 L N (with guess w) denoted by (X,)) SPa(X;).
7:
8:
Consider X (X/). Remove X from Pa (X?) if X -T X
UPP(X)) \ Xx_r using a CI Test with samples & E ( 2 2Tuh, / E m...
9:
Store Pa (X²) for Xi, t (t:t 2Tbb,& t 11
end for 11: end for
10:
12:
w f arg c(ww) maxkelw]
(X ( )l. Set Pa(X() + Pas, (X²) for x2, , E {t : 2 2nbb}
13:
14:
end for
15: return W, and ( Vj E [n]. 2T
PCMCIa builds the Algorithm PCMCI in Runge et al. 2019 Additional details about PCMCI are provided in the supplementary material.
Overview of Algorithm PCMCIn- We assume that the periodicity and time lag upper bounded by Wab and Tuh respectively. Using PCMCI Runge et al. (2019). we obtain a superset of parents for every variable X denoted by SPa(X() (line 2). Our goal is to identify the correct set of parents along with its periodicity for every variable in V. For variable x² guess its periodicity w by iterating over all possible values in [wabb]. Next, construct time partition subsets II; k E [w] based on the guess of periodicity w. In each time partition subset, we maintain parent set, denoted by Paw (X²), initializing it with the superset SPa(X?). Then test the causal relations between X;_, E (X?) and X! using a CI test on the sample t fi (lines 6-10).
For each guess w, every variable in X should have its estimated parent set (line 9), and there are total potentially different parent index pludketa) X3, We return an estimate W, that maximizes the sparsity of the causal graph (Lemma 3.4. Therefore, we select the value w E [wab]] that minimizes the maximum value of (P/) [T] as the estimator of wy (line 12).
3.1 Theoretical Guarantees
Our main theorem shows that PCMCIa recovers the true causal graph on discrete data. There are three important lemmas. We provide all the detailed proof in the supplementary material.
Theorem 3.1. Let & be the estimated graph using the Algorithm PCMCIg- Under assumptions AI-A7 and with an oracle (infinite sample size limit), we have that:
δ Ō=0 G
(13)
almost surely.
Lemma and Lemma jointly state that if CI tests are conducted on samples generated by different causal mechanisms, the obtained parent sets SPa(X:) should be the superset of the union of the true and illusory parent sets Ug Pay (X That is, the estimated graph should be denser than the correct graph. The true parent set can then obtained by directly testing the independent relations between the target variable X and the variables in SPa( (P))), assuming a consistent CI test. Note that the CI tests in algorithm are assumed to be consistent given i..dd. samples. We do not assume

---

the consistency of CI tests with respect to semi-stationary data. Therefore, any CI tests that maintain consistency with i..d. samples can be seamlessly integrated into our algorithm.
w,)) kei,,) contain the true and illusory parent sets, where wy is the true periodicity of X'. For any random variable X? with large enough t, under assumptions Al-A7 and with an oracle (infinite sample size limit), we have:
Lemma 3.2. Denote that (Xax(()) kEE
p p(x22 Pa (XF)) + p(x)) UiLn Pax \u)) = 1, by E Pan())) (14)
Here, p(X - (x?)) - p -
Pak (x)() where p(X;1 U Pat(X))) is
an estimated conditional distribution using all samples / E [maax + 1,7::
E, 1 U Pa (Xx( E, (
x - Pax (xF))
(15)
Proof sketch. We argue that the estimated conditional distribution in Eq. 15 can be written as a linear combination of p(X? [Pa(X?)) where & E I1 k E [wj]. i.e., as a mixture of conditional distributions. The coefficients in the linear function, say Ouk, k [way], can be further decomposed based finer time partition called the homogenous time purtition, which consists of subsets constructed according the Markov chains ge corresponding to the time series. Based on Assumption A7, the Markov chains are stationary and ergodic. Therefore, after sufficiently large time steps, the it achieves unique equilibrium. With this type of stationary sample, we can express Og by joint distributions instead of the indicators. Then, we can complete the proof of our inequality claim in Eq. 14 using Assumption A2 and Bayes theorem.
distribution of {2i] ge()) will be invariant across
Lemma 3.3. Let SPa(X() denote the estimated superset of parent set for X E V obtained from the Algorithm (line 2). {Pau (Xi)keEw,) contain the true and illusory parent sets, where wi is the true periodicity of X'. Under assumptions A A1-A7 and with an oracle (infinite sample size limit), we have that:
J R=1
Pax(X() C SPa(X)), W [Tmax 1,T]
almost surely.
Proof sketch. Assume the contrary, i.e., there exists s E Up. Pag (X²) \ SPa(X?). From Lemma 3.2 we have X Pax \ s. By the Definition 2.4 we have Pa(X() Pax If 8 & Pa (X)), by the causal Markov property (Assumption A2), the dependence relation can not be true because s is a non-descendant of X/. If s E Pa(X/), our Algorithm would have concluded that X (Pa())) (line 2), evident from the causal Markov property, contradicting our assumption. Hence, the lemma.
Based on Lemma 3.2 and Lemma 3.3 we can identify the true wj for X through Lemma 3.4
Lemma 3.4. Let Wj denote the true periodicity for X E V and Pa, (X. (an)) denote the estimated
purent set for X3 obtained from Algorithm 1 line where E ITC. Define:
B arg min max (X tang wE(wa)) kE[w]
(16)
Under assumptions A1-A7 and with an oracle (infinite sample limit), we have that رo = wjs Vj E [m] almost surely.
Proof sketch. Assume the contrary that Wj # Wj. then in the Algorithm ! we have an incorrect time partition ÍP' Hence, CI tests that are performed use samples with different causal mechanisms.
7

---

1.0
5. Times
5, Tiue
0.8
Time (sounnits)
0.6
T=1000
Arrumay
100
T-2000
0.4
T-1000
50
0.2
T-2000
T=4000
0.0
1
2
3
4
6
7
8
9
10
Wimax
(a) Accuracy of
(b) Runtime (in sec.) when Wmax varies
Figure 2: PCMCIQ is tested on 5-variate time series with Tmax 5. Set Tub 15, with 15 for all variables. Every line corresponds to different time series length. Every marker corresponds to the average accuracy or average running time 100 trials. a) The accuracy rate of à for different time series lengths and different Wmax- b) Illustration of Runtime (in sec.) when wmax varies.
5 (a)(X in Eq 13 is estimated from mixture of two more time partition subsets, say 111 and 1zz. We can apply Lemmal 3.2 with Pa With Lemma 3.3 L Pag. (X/) S
Pazy (X²). Hence, for K.j. [Payy (x2) 2 J (Xy)) using mixture samples / E U 1 IT
For true wj, we have PPa, (X/)) = [Pa(X With Assumption A6 the Hard Mechanism Change,
L (aa(X()) > (P((X)) so that w, always leads to smaller size of estimated parent sets than 6j. contrary to the definition of aj. Hence, ش, = W/.
Experiments
4.1 Experiments on Continuous-valued Time Series
To validate the correctness and effectiveness our algorithm, we perform a series of experiments. The Python code is provided at we test four algorithms PCMCIg. PCMCI Runge et al. 2019]. VARLiNGAM Hyvärinen et al. 2010 and DYNOTEARS Pamfil et al. 2020 on continuous valued time series with Gaussian noise. The experiments for continuous- time series with exponential noise and binary-valued time series in the supplementary material.
In this section,
Following Runge et al. 2019), we generate the caatinuous-valluee time series in three steps:
Construct an n-variate time series V with length T using independent and identical (Standard Gaussian Exponential) noise temporarily. Determine Tmax and Wmax where Wmax max(w;) (e After making sure that one univariate time series, say X2, has periodicity Wmaxx the periodicity the remaining time series X', * j is randomly selected from {1,... (muxx respectively.
2. Randomly generate wi binary edge matrices with shape (n, Tmax) for each time series jj E [n]. denotes an edge and 0 denotes no edge. Each binary matrix represent one parent set index pInd k E [wy]. Randomly generate Wj coefficient matrices with shape (n, Tmax) for each time series x E [n]. One binary edge matrix and one coefficient matrix jointly determine one causal mechanism. Hence, total W/ causal mechanisms are constructed. Here, make sure that V satisfies Assumption A6.
3. Starting from time point > Tmaxx generate vector X; over time according to all the causal mechanisms of V. until t achieves T.
We did not conduct experiments on JIT JIT-LiNGAM because this is from recent paper Fujiwara et al 2023 and is considered concurrent per NeurIPS policy.

---

500
T 2000
8000
T 2000, =5
1.0
1.0
0.8
0.44
A. Score
0.6
Predsion
Recall
2 0.2
0.0
0.4
3 9
5
1
0.2
1.0
Trediion
0.8
10
12
5
10
12
10
15
Traa
0.2
1.0
2000,
Wlear
5, Trax
0.0
3
1351 9 5351 9
1.0
S.oore
0.6
Precision
Recall
0.8
0.6
0.4
PCMOla
0.2
PCMC
POMCI
PCMC
0.0
5
5 7 6hnxx
VARL INGAM
0.0
2
10
15
10
15
20
VARLINEAM:
DYNOTEARS
DYNOTEARS
Gaussian noise
(b) Performance when Tmax or varies
Figure 3: 4 algorithms tested on 5-variate time series. Set Tub 15, wis 15 for all variables. Every line corresponds to a different algorithm. Every marker corresponds to the average performance over 100 trials.
Following the previous work Huang al. (2020]. Fi score, Adjacency Precision, and Adjacency Recall are used to measure the performance of the algorithms. The details of calculating these metrics are described the Appendix. All the performance statistics are averaged over 100 trials. The standard error of the averaged statistics is displayed either by color filling or by error bars.
correct estimator accuracy rate of @ for different time lengths T'. Here, elements in {Nw;} } where N Wab are all treated correct estimations. By Definition 2. and 2
the prerequisite for obtaining the correct causal graph. Fig 2(a) shows the
the multiple of w, is still associated with correct causal graph. However, it leads to finer time point partition IT' decreasing the sample size used in each CI test from approximately T/ww to approximately T/(Nw;;) The accuracy rate is sensitive Wmax for small T. This result verifies that algorithm PCMCIa has the capacity to detect the true periodicity of each X E with large enough time length.
We evaluate the performances of PCMCIa on continuous-valueee time series with Gaussian noise shown as Fig As T increases, The sub-figures show that all three evaluation metrics decrease when Wmax gets larger. The precision of PCMCIa is always far better than other algorithms when wmax not equal
natural to see a continuous improvement in performance.
Given the fact that the parent sets Pa(X)) Vj, obtained from PCMCIg are subsets of the parent set SPa(X/) Vj. estimated from PCMCI, the recall rate of PCMCI should be the upper bound of the recall rate of PCMCIp. This assertion has been verified as the red recall line of PCMCI is always below the blue recall line of PCMCI 7 increases.
In Fig ((a) the recall of PCMCIq worse than PCMCI for T 500. In this regime, the accuracy rate of low, shown as the dark blue line in Fig 2 a sparser causal graph. Hence the number of true positive edges may decrease. This a common problem for many constraint-bast based algorithms, but it hurts PCMCIp the most because in PCMCIp, the sample sizes in each CI test are approximate T / instead of T. As T increases, the red recall line of PCMCIg push forward to the blue recall line of PCMCI. The high value of both adjacent precision and recall rate with large T verify that PCMCI, can identify the correct causal graph.
Small sample sizes in CI tests may result in
We also observe the performance of our algorithm as Tmax and N varies in Fig 3(b) As the performance of PCMCIa consistent over n-variate time series with different n, large Tmax may lead to smaller precision and recall rate.

---

The parent sets for each variable obtained from PCMCIa algorithms are shown in Table | Sets of true and illusory parents of variable time are separated by curly braces. For instance, variable X with Wwp 1 means that the causal mechanism of the surface pressure in the West Pacific remains invariant over time with the estimated parent set {X X 1-11 x1=1]- Only time 1-1 series X has three different parent sets, including one true parent set and two illusory parent sets, which appear periodically over time. The three parent sets of X imply that the causal effect from the tropical Atlantic air temperature X the Central Pacific temperature X;" would disappear every quarter of year. Note that we do not have a ground truth in this case, and we do not possess Case Study
4.2
Here, we construct an experiment with real-world climate time series dataset. In Runge et al. 2019]. the authors tested dependencies among monthly surface pressure anomalies in the West Pacific and surface air temperature anomalies in the Central Pacific, East Pacific, and tropical Atlantic from 1948 to 2012. Our application explores the causal relations among the monthly mean of the same set of variables from 1948-2022 with 900 months. Let X, denote the monthly mean of surface pressure in the West Pacific, X X and X denote the monthly mean of air temperature in the Central Pacific, East Pacific, and tropical Atlantic, respectively.
the necessary knowledge in this area, so the significance of these results under-explored. More discussion about this application can be found in the supplementary materials.
Table I: Climate application results estimated from PCMCIa.
PCMCL
X
w
(Pax kelp true and illusory parent sets
X
X X, X { 1 +-2:
1-1
X.P
3
}:
{X, X; x_ X {x-
X
X
X
{X.P X X/a 2, X,
X
{X] 12-11-0 X
Conclusions
In this paper, propose non-parametric, constraint-based causal discovery algorithm PCMCI designed for semi-stationary time-series data, in which a finite number of causal mechanisms are repeated periodically. We establish the soundness of our algorithm and assess its effectiveness continuous-valle valued and discrete-valued time series data. The algorithm PCMCIg has the capacity to reveal the existence of periodicity of causal mechanisms in real-world datasets.
Acknowledgements
This research has been supported in part by NSF CAREER 2239375 and Adobe Research. We wish to convey our heartfelt gratitude to the anonymous reviewers for their invaluable and constructive feedback, which significantly contributed to enhancing the quality of the manuscript.
10

---

References
Judea Pearl. Causality: models, reasoning, and inference, 1980.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithmss The MIT Press, 2017.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Wicher Pieter Bergsma. Testing conditional independence for continuous random variables. Citeseer, 2004.
Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kemnel-based conditional independence test and application in causal discovery. arXiv preprint arXiv:1202.3775, 2012.
Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. 2020.
Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science advances, (1 :1)aaau9999, 2019.
Andreas Gerhardus and Jakob Runge. High-recall causal discovery for autocorrelated time series with latent confounders. Advances in Neural Information Processing Systems, 12655-1552,,, 2020.
Tianjiao Chu and Clark Glymour. Search for additive nonlinear time series causal models. Journal of Machine Learning Research, 9(5), 2008.
Doris Entner and Patrik O Hoyer. On causal discovery from time series data using fci. Probubilistic graphical models, pages 121-128, 2010.
Daniel Malinsky and Peter Spirtes. Causal structure learning from multivariate time series in settings with unmeasured confounding. In Proceedings of 2018 ACM SIGKDD workshop on causal discovery, pages 23-47. PMLR, 2018.
Azam Ikram, Sarthak Chakraborty, Subrata Mitra, Shiv Saini, Saurabh Bagchi, and Murat Kocaoglu. Root cause analysis of failures in microservices through causal discovery. Advances in Neural Information Processing Systems, 35:31158 151111180,1,,,, 2022.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer. Estimation of structural vector autoregression model using non non-gaussiamity. Journal of Machine Learning Research, 11(5), 2010.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Causal inference time series using restricted structural equation models. Advances in neural information processing systems, 26, 2013.
Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Geor- gatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from -eme-s data. In International Conference on Artificial Intelligence and Statistics, pages 1595-16000 1605. PMLR, 2020.
Charles K Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. Journal of Artificial Intelligence Research, 73:767-819, 2022.
Fang Han, Shyam Subramanian, Edwin R Price, Joseph Nadeau, and Kingman P Strohl. Periodic breathing in the mouse. Journal of Applied Physiology, 9333)1133-1110, 2002.
Akira Nakamura, Yasuichiro Fukuda, and Tomoyuki Kuwaki. Sleep apnea and effect of chemostimu- lation on breathing instability in mice. Journal of Applied Physiology, 94(2):25-533, 2003.
Mary Carskadon, William C Dement, et al. Normal human sleep: an overview. Principles and practice of sleep medicine, 4(1):13-23, 2005.

---

Sandra Komarzynski, Qi Huang, Pasquale F Innominato, Monique Maurice, Alexandre Arbaud, Jacques Beau, Mohamed Bouchahda, Ayhan Ulusakarya, Nicolas Beaumatin. Gabrièle Breda, et al. Relevance of mobile internet platform for capturing inter-and intrasubject variabilities in circadian coordination during daily routine: pilot study, Journal of Medical Internet Research, 20 (6):e204, 2018.
Mingming Gong, Kun Zhang, Bernhard Schoelkopf, Dacheng Tao, and Philipp Geiger. Discovering temporal causal relations from subsampled data. In International Conference on Machine Learning, pages 1898-1906. PMLR, 2015.
Biwei Huang, Kun Zhang, Mingming Gong, and Clark Glymour. Causal discovery and forecasting in nonstationary environments with state-space models. In International conference on machine learning, pages 2901-2910. PMLR, 2019.
Christopher Krich, Jakob Runge, Diego G Miralles, Mirco Migliavacca, Oscar Perez- PerezoPriego, Tarek El-Madany, Arnaud Carrara, and Miguel D Mahecha. Estimating causal networks in biosphere- atmosphere interaction with the pemei approach. Biogeoxciences, 1734): 2020.
Yuquan Qu, Carsten Montzka, and Harry Vereecken. Causation discovery of weather and vegetation condition global wildfire using the pemei approach. In 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS, pages
IEEE, 2021.
Liangfeng Zou, Yuanyuan Zha, Yuqing Diao, Chi Tang, Wenquan Gu, and Dongguo Shao. Coupling the causal inference and informer networks for short-term forecasting in irrigation water usage. Water Resources Management, pages 1-23, 2022.
Giovanni Menegozzo. Diego Dall' Alba, and Paolo Fiorini. Causal interaction modeling on ultra- processed food manufacturing. In 2020 IEEE 16th International Conference on Automation Science and Engineering (CASE), pages 200-205 IEEE, 2020.
Sayeh Peterson. Comparison of Lasso Granger and PCMCI for Causal Feature Selection in Multi- variate Time Series. PhD thesis, The University of Arizona, 2022.
DK Arvind, Sharan Maiya, and P Andreu Sedeño. Identifying causal relationships in time-series data from pair of wearable sensors. In 2021 IEEE 17th International Conference Wearable and Implantable Body Sensor Networks (BSN), pages 1-4. IEEE, 2021.
Luca Castri, Sariah Mghames, Marc Hanheide, and Nicola Bellotto. Causal discovery of dynamic models for predicting human spatial interactions. In Social Roboties: 14th International Conference, ICSR 2022, Florence, Italy, December 13-16, 2022, Proceedings, Part 1, pages 154-164. Springer, 2023a.
Luca Castri, Sariah Mghames, Mare Hanheide, Nicola Bellotto, et al. Enhancing causal discovery from robot sensor data in dynamic scenarios. 2023b.
Daniel Malinsky and Peter Spirtes. Learning the structure of nonstationary vector autoregression. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2986-2994. PMLR, 2019.
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph D Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bernhard Schölkopf. Causal discovery from hecerogencouu/ioontationnayy data. J. Mach. Learn. Res., 21(89): 1-53, 2020.
Daigo Fujiwara, Kazuki Koyama, Keisuke Kiritoshi, Tomomi Okawachi, Tomonori Izumitani, and Shohei Shimizu. Causal discovery for non-stationary non-linear time series data using just-in-time modeling. In 2nd Conference on Causal Learning and Reasoning. 2023.
Richard Serfozo. Basics of applied stochastic processes. Springer Science & Business Media, 2009.
Samuel Karlin. A first course in stochastic processes. Academic press, 2014.
12

---

Appendix
A PCMCI Algorithm
The PCMCI algorithm is proposed by Runge al. 2019 aiming to detect time-lagged causal relations in a window causal graph. There are two stages of PCMCI: the condition-selectie stage and the causal discovery stage. In the first stage, unnecessary edges removed based the conditional independencies from an initialized partially connected graph where Assumption A4-A5 should be satisfied. In the second stage, Momentary Conditional Independence tests (MCI) used to further remove the false positive edges caused by autocorrelations time series data. More specifically, these steps can be briefly formalized as follows:
PC₁ in Algorithm Al part of the PC algorithm in a more robust version named stable-PC al. (2016] The goal
Condition selection stage. PC1 is variant of the skeleton-discovery
in this stage is to obtain superset of the parents Pa( (a()) for all variables (E[Tamx E V.
X
Initialize Pa(X)) (Xi r/ Pa x will remove X if
x 11 X Pa(x))
MCI in Algorithm A2 Causal discovery stage. In this stage, do MCI tests for all variable pairs xx)+,X(] with j E [n]] and time delays [Tmax]:
where Pa(X() and Pa(X are estimated from the PC1 stage.
Note that Tinax in this section the same as Tub in the main paper, serving as the upper bound for the time lag that exhibits causal effects. On the other hand, Tmax in the main paper denotes the maximum time lag observed within the multivariate time series. Essentially, in the main paper, Tub is parameter that must be fed into the algorithm, and Tmax is observed from the true causal graph. As default, assume Tuh is configured with a value greater than Tmax, ensuring that the algorithm uncovers the correct causal relations. See Fig for more detail.
X
x²
x²
10
11
12
13
14
15
Tmax max(t; 2,,, 3, =
Tab
Figure 1: Set Tub to be 5, then all parent candidates of variables at f = 15 are included in the large orange box, ranging from 10 14. Consequently, the algorithm will only examine causal effects with time lag not exceeding 5. In the causal graph, Tmax is 3, representing the maximum time lag observed among the 3-variate time series. Specifically, the maximum time lag for each component time series is 72 3, T3 1, respectively, and Tmax represents the largest value among these three maximum lags.

---

Algorithm A1 PC
1: Input: A n-variate time series V = xx',,, x"),, target time series X maximum
time lag 7max, significance threshold apc, maximum condition dimension Pmaxx (default Puax (Tmax). maximum number of combinations qmax (default qmax 1), conditional independence
test function C'I
2: function CI (X, (X,Y.Z) Y.Z)
3: Test X IL Y using test statistic measure
4: return p-value, test statistic value /
5: Initialize preliminary set of parents Pa(X)) {Xi_, :iE fl,...,,,,77 fi,... Tmax }
6: Initialize dictionary of test statistic values I (Xi-r + X?) oc VXi-r E Pa(X!)
1..2..... .... Pmax do 8: if ([ < then
7: for p 00 1,2,
9: Break for-loop
end
11:
for all Xi_, in Pa(X:) do --
12:
13:
for all lexicographically chosen subsets S C Pa(X)) (x-_,) with [S| = do
14:
q =
15:
if' q 2 4max then
Break from inner for-loop
17:
end if
18:
Run CI test obtain (p-value, 1) C((xi_x:X
19:
if [ < m (x: r- + x then
D~Store min. / of parent among all tests
20:
7min (X x |7
21:
end
22:
if p-value arC then
Removed only after all X have been tested
23:
Mark X;_ for removal from Pa(X/)
Break from inner for-loop
25:
end if end for end for
27:
Remove sinniinnif parents from Pa(X?)
29:
Sort parents in Pa(X/) by pXi - X;) from largest to smallest
end for
31:
return Pa(X()
Algorithm A2 MCI
1: Input: A n-variate time series V
X2 X X"), sorted parents Pa(X)) for all
variables X estimated with Algorithm
A1 maximum time lag Tmax, maximum number px of parents variable X and conditional independence test function C'I
2: for all (Xi_r, X!) with i.j E (1,...,n),,, E Tima) excluding (X/, x!) do
3: Remove X;_, from Pa(X)) if necessary
4: Define Papx (X:--) as the first parents from Pa(Xi). shifted by
5: Run MCI test to obtain (p-value, I) + CI(Xi X² = (Pa(X)), Papx (X/))))
6: end for
7: Optionally adjust p-value of all links by False Discovery Rate-approach (FDR)
8: return p-vauue and MCI test statistic values

---

B pCMCI
For simplicity's sake, define sets: [t] == b) and [a, = {a, a + 1, I,... bj, where a, b E N.
Algorithm B1 PCMCIa
1: Input: A n-variate time series V = (X x² x3 ....., )"), periodicity upper bound wab, time
lag upper bound Tub- By default, we assume Tub and Wab are larger than their true value.
2: A superset of parent set is obtained using PCMCI with Tub and denote by SPa(X!) Vj,
3: for X where j E [n] do
4: for a guess w E [wbb]] of do
5: Let II := [ E [w]] where I.ç (2 + k; 2Tu + RTbbub + 2w + k.
6: for E [w)] do
Initialize the parent set for X t E {t: t 2 2Tub, t E ( (with guess w) denoted by
7:
Paw < SPa(X?).
8:
Consider Xi E (X)( Remove X from Pa)(X if Xi 1T X! | ( \ [X,,, using CI Test with samples t E {t: t 2 2Tb,, t m.).
9: Store (X!) for X2, E (t 1 2 Tub. t f.
10: end for
11: end for
12:
13:
if there exists turning points Sj, S E [wib then
14:
@3 min
min Sy
15:
else
16:
17:
3g + arg we(wa) maxkelw) Paw (X') reñt! )l-
18: end if
19: Set Pa(X)) Pas, (XF) for X/ {t: 2
20:
end for
21:
return رj and Pa( [m], 2 2T
C Soundness of PCMCI
C.1 Stationary Markoy Chain
Any discrete-valued time series V with Semi-Stationary Structural Causal Model (SCM) satisfying assumption A1, A2, A4, A5 can be written chain satisfies Pa( Z,,) C Zn U Zn-1 for all n, where Z,, is set of variables in V. This Markov chain has a finite number of states if all time series in V are saccete-valued time series.
Claim:
Markov chain {Z,) as long as this Markov
Note that when the notation n is related to Markov chain 2n, it means the running index. In the context of xx?E(m), represents the index of component time series within the n-variate time series.
To simplify, assume that one associated Markov chain of V {Xe, Y,, Xx-1, Y--1) with t E {t E N+ :,t < T} satisfying Pa Pa(Zn) C Zn U Zn-1- Here, the notation for the time points of variables is simplified as and function of n. the running index of the Markov chain. Note that Zn-1 (X+-2, Y,- X,-3, Y,-3) rather than (X,-1, Y,-1, X, Y--2). 2), as the simplified notation could erroneously suggest the latter sequence. A simple proof is shown below through Markov assumption (A2).
{X. Y) has Zn
even though it should be

---

21
x
X
X
C:
11
12
13
14
15
Figure 2: Partial causal graph for 3-variate time series V -= -X X2 X with Semi-Stationary
SCM where Tmax = 3, W1 = 3, w2 = 2, W3 1, 6 and 8 = 6. The first 3(=Tmax) time slices IXsss <isa are the starting points. The same color edges denote the same causal mechanism. E.g. for X there are 3 (= wy) time partition subsets [ The time points / of nodes X; sharing the same filling color are in the same time partition subsets. The time points of nodes X; sharing both the same filling color and the same outline shape are in the same homogenous time partition subsets. There are 6 (= 8) different Markov chains in this multivariate time series V, and the first element of these 6 Markov chains is shown as { {Z7) S5555 and tinted with a gradient of blue hues.
[X:)
Zi and 2 denote the first two elements of the first Markov chain while Z and Z denote the first two elements of the second Markov chain.
Proof.
pl Zn 1ii--1,Z22220 Zn- p(x,, Ye, X1-1, Y 1 2n-
(3)
Zn, Zn- |x,, -n- )p YelZa U Un-1 (X,
(5)
= -)(x//a(Xx)) Z, U. - (X
(6)
xP Y.)Pa(Y)), 2, J2 1000,1XXX,,,10 (X,UX,
xp XP-iPPa(X,-)) Zn - /X..uXXXUUUUU.U.UU U UPA(X
p,,,/,, - Xx,p
Y2 U2-7 1
=((,, V1, X Y (2-11 = (m
(8)
Assume that the space of both X, and Y, with < T are (1,2). There are total 2 16 states of Markov Chain { ((X,, Y,, Y-1)) The transition probability P for this Markov Chain
illustrated as a 16 x 16 matrix:
(1,1,1,1)
12,1,1,1)
(11111
P1,1
P1,2
P (2,1,1,1)
P2.1
P2.2
16X16
where (1, 1,, 1, means X, 1,Y, 1, 1, Each in this transition probability
matrix is conditional distribution of Z,, given one realization of Zn Each entry is probability
of having one specific realization Z,, given one realization of 7n-1. This probability can be
decomposed by conditional distributions based on Markov assumption (A2). Take P1,1 an example:
Pl,1 = p(X, 1, 1,, = 1|X,-2 Y, 1, X,-3 1,
(10)
P(X, -
1\Pa((,-))))) (11)

---

where Pa(.) here are realizations, not random variables.
For time series V with Semi-Stationary SCM, there are (potentially) å different Markov chains {2%), E [o]:
2 = (X Tmax Tmax tg+ 1)6, X moxx + (n 1)8,
X to
where n E {n: n E N+, -max + q 1 + nō Σ T), [mmst 12. As proved in the claim, such a Markov chain exists as long as Pa(Za) C 2% - for all n. The value of can guarantee the existence of such Markov chain because 8 is larger than Tmax + 1 and is multiple of Ω2 that is, multiple of all {wy je(n)- By doing so, Pa(Zy C Zg U - is satisfied; for any variable X7. there exists E [6] and n E N+ such that variable X and its parent set Pa(X;) can be included in 7g: and the causal mechanism generating Z is invariant for different The state space of (Zx) is the set containing all possible realizations of (X, Tmax t0
1)8 ) (ej)).eNN-- The transition
probabilities between the states are the product of associated causal mechanisms based on Markov assumption (A2).
Determined by the starting slice X, where Tmax <t I Tmax + 6, there should be 8 potentially different Markov chains {2X} where S S 8. To be more specific, those Markov chains are:
Markov Chain 1: Z, {X T +1+ 1)6 X Tmax +2+ (n-1)6- 136, X mnd
(12)
where {n N Tmax nō S T)
Markov Chain 2: 2 {X T +2 136 X +3 186, X +1m66)
(13)
where n (n: No, Tmax + +1+ nd S T).
Markov Chain 8: Z
[X TTmax +n8 X Tmaxe +max X Tams - 1+++1)66
(14)
where {n N+ -max - 11 (n + 1)8 - T).
Given Irreducible and Aperiodic Markov Chain assumption (A7), discrete-time Markov chain {7% Jocn; 9 E [8] with finite states should be a stationary and ergodic Markov chain, and there is a unique stationary distribution Fy (Bertsekas and Tsitsiklis (2088), Karlin (0044). Additionally, the large power of the associate transition matrix P, will eventually converge to a matrix in which each row is the stationary distribution TT. Equivalently,
lim pZZ alz? b) = p(Z = Va.bs.. 200
(15)
where S is the state space of Z
In other words, after sufficiently long time, equivalently, is large enough, the distribution of (2A) does not change with increasing n. That is, for large enough n::
p(Z = p( (2 Vn1, n2 > n.
(16)
Returning from the stationary and ergodic Markov chains {2g}, E [8] back to the original data
V through Eq. 12) to Eq. 14, the distribution of the original data must adhere to the following
condition:
p( X mon X Troex +g+ 1 +oo
pl X maax tog X
Twaa + 4nggs vnry
Tmax +
(17)
for any q E [6] and A1, n2 >> n.
Given these clarifications, yet finer than, the time partition defined Definition 2.3 the main paper.
naturally introduce
refined time partition that is based on,
Definition C.1 (Homogenous Time Partition). For univariate time series X in a Semi-Stationary SCM with periodicity wj. the time partition 112 of X can be further divided into series of non- overlapping and non-empty subsets (T(K..) 4. For each t E [Tmax + 1,7], there exists
E [wس] so that & IT. and further there exists E [,,) so that E ((k.s)" Tik.o) can be written as:
- ))+ 11 s). (18)
3 (k,s) = {t: Tmax + 1 S T, (t mod w/) + 1 = k, (t mod
Wj
5

---

With this definition, we have - (k,,)) IT. While time partition IT, guarantees that all variables in {X/) share the same causal mechanism, homogenous time partition T(k,s) guarantees that all variables in {x/' /eseseeee share the same distribution where t' represent the steps needed by the associated Markov chain achieve equilibrium.
2 shows a partial causal graph for a 3-variate time series with Semi-Stationary SCM. Tmax = 3 means that the causal mechanisms start from /. 4, and the random variables with / 2,,2,3) are 2,3)
random noises. For the first time series X the periodicity wi is 3. And the periodicity of the time series X² and X is 2 and 1, respectively, The periodicity of the whole time series V obtained by LCM(3, 2, 1) = 6. 8 = -m S2
S+11 x 6 = 1 x 6 = 6. 6
In Fig 2 periodicity W1 = 3 means that the causal mechanisms repeat every three time points and hence there are three time partition subsets 11 E [3]. More specifically, II = (4, 7, 10, 13, ..... 4 3.V,...., {5,8, 8, 11, 14,..., 5 + ..., ITg = {6, 9, 12, 15,..., 3N, where E N+ Random variables {X²) with L in the same time partition subset share the same causal mechanism. However, they may not share the same marginal distribution.
Still in Fig Q based on the definition of homogenous time partition, time partition subset 11; for X can be further decomposed as T/K "(k=1,s=1)> = {4, 10, .... + 6N, 6..,..., m (k==1.==)) = {7, 13,
J....... where E ]. After long run n, Z, and Z,+1 will eventually share the same distri- bution, that is, all the variables inside 2 will share the same joint or marginal distribution as the corresponding variables inside Zh+1- To illustrate this, we assume that this Markov chain has already achieved its equilibrium at time point t 4. Based Eq 12 and Eq we have:
p(X4: X X9) p( 10, XI1, .... X16) P(X X
X,,, ...
(19)
From the identical joint distribution, we can further have:
x( p(Xio) - p(Xis)
(20)
as X1 E X1o E Xio and X1 E X16-
Therefore, for sufficiently large values of ensuring that 2,, has reached its stationary distribution, all variables within (X-) (eEms)))) will share the same distribution.
In Z there are 6(= 5) potentially different Markov chains {2%} E [o] in V. For any time window with length 8, (X,,.. .... X(+6-1), there exists q E =10), E N so that this time window can be completely included in Zg. For instance, set {X,, X Xe...., X is in Z7, which is the first
element of Markov chain {
Constructing Markov chains and applying the Irreducible and Aperiodic Markov Chain assumption (A7) enable us obtain a consistent estimator for the conditional and joint distributions of interest.
C.2 Consistent Estimator
The conditional distributions for variables in (x:) are the same, that is, p(1/,, (Pa(x),))
p(xy, (P(((())) Vh,t2 E TT. For simplicity, denote
Prem (a/(Pa(x=)) = p(a:/|Pa(()) Vi ITŞ
(21)
Consider an indicator function such that 1 (w/ Pa(x))) = if configuration (x/ Pa(a:)) has realized, otherwise 1 (x/, Pa(x))) 0.
Since every t E (k,,s) is apart from each other with N8 steps where N E N+, and there must exist q E [رس]] and E N+ so that {xf, Pa(a()} E Z%, then for the same 9. there must exist another n2 so that {at Pa(r/4N5)) E Z% Hence, we have {1(x], Pa(z))))
(eer with some function f R2 + IR satisfying (Z ( < co. Since the value of deter- mines q and n2, use qt) and (t) to emphasize their relations. For large enough t > t,.
= (f(2 ) (k..)

---

Pa(x)))())ssse are identical samples where t' is the time point needed by the associate
Markov chain to achieve its equilibrium after ra ((') steps.
Without loss of generality, we assume T is multiple of á all the time.
We can construct an estimator of p(z/, Pa(x))) with large enough t as:
5
p(((,)((x))) T
Σ 1( Pa(x)) tE#
(22)
where k, s is determined by and there must exist one and only one k. satisfying & E (kke)) Now, we are going to show this estimator is consistent.
We first decompose the estimator into two parts: time point S t' and t t', where t' represents the time point when the equilibrium of the associated Markov chain is achieved.
Σ
1, Pa(x())
(23)
(EST
8
Σ
(al,
Σ
(24)
T
(St (Sti,tEF().))
8
8
Σ
((tl,
Σ
(a(a))
(25)
T
T
1S1 (E
(E
8
8
T 8
Σ
(a(,)((z))) +
Σ
(a(x))
(26)
T
T t'
5 T
1>t1,EE (k..)
6
8
T
1'
Σ
(al,Pa(z)) +
Σ
Pa(x))
(27)
T
T t'
6 T
(SU, (K.))
(>V (h.a)
8
T - t'
Σ
(a(x)) +
Σ
(28)
T
T
T - t
(SH (B.))
(29)
Take a limit of Eq. 23
we have:
lim T
Σ
1(2), Pa(x())
(30)
T-t t' 8
= lim
Σ
(a(,)) + lim
Σ 1(2 Pa(x)))
(31)
T0
T=20
T
T-
(SU,EEF(K
0 lim
T --t t'
1
Σ
r(2a(()) where i 1,, E T(P.0)
(32)
-00
T
(T) - (1())
(1())21(())
Birkhoff's Ergodic Theorem
0+ E ((ZZ(()))
(33)
=E( [x/] Pa(x))))))) where t > t', E T(k
(34)
= p(o), Pa(a))), where / 1', / E <<(k.))
(35)
Denote
PiEm (v,, Pa(x()) = p(at, Pa(a!)), where ( 1>1' (k.))
(36)
7

---

Based on the definition of homogenous time partition and time partition, PremiE.)
(./(Pa(x)))
Preng (a/[Pa(z)) Vs
Similar to Eq. one estimator of Preug (x) (a(aa)))) vk = [w)]] is
Pre (x/(Pa(x:)) = Eveng 1 (z Pa(a)))
(37)
Eieng I(Pa(z)))
Σ EiER (k,)) ((t),
(38)
EieRy 1(Pa(1!))
Σ EIEFiO.,
1(2/ Pa(x!))
(39)
E $ EIEnt 1(Pa(1]))
Take a limit of Eq. we have:
lim P Prem: (1) (a(()))) T-000
(40)
4 Σ PIE ()
(of Pa(x()
Eq
(41)
E PrET(RE)) (Pa(z)))
(xf (a(( (Pa(x!))
s= PIE
(42)
E (Pa(z!))
Pier) (e/(Pa(x))) same for all (t (h:0)
Preng (oy Pa(x()) E b PIET(A.)) (Pa(x)))
(43)
E PrEm/en) (Pa(x?))
Prent (r((Pa())))
(44)
Hence, Preng (x. Pa(æ()) is consistent estimator of Preng (x/(Pa(zz))
Similarly, construct estimator of p(z/1 Un Pan (a))) where t E [T]: p(x)1 Un Pay (xt)) = Σ 1x0/ Un Par (xx))
(45)
E. Il (zf Un
(46)
Σ, (an,
We will prove that this estimator converged as goes to infinity in Lemma D.2 Hence, it is consistent estimator.
In this section, we have proved that p(æ Pa(x))) in Eq. is consistent estimator of p(a?,, Pa(w:)) using samples with in the same homogenous time partition subset and p(x² (Pax=))) Eq. is consistent estimator of (P(x:)) using samples with t the same time partition subset.
D Theorem
Theorem D.1. Let be the estimated graph using the Algorithm PCMCIg. Under assumptions Al-A7 and with an oracle (infinite sample size limit), we have that:
G G-0
(47)
almost surely.

---

Lemma D.2. Denote that {Pak (X)) contain the true and illusory parent sets, where wy is the true periodicity of X3 For any random variable X with large enough i, under assumptions Al-A7 and with an oracle (infinite sample size limit), we have:
p [ U (x)( * p(X?1 L#L1 (xx(X)) v)) = 1, Vy E
(48)
Here, p(X (x)(X lim p(X (x)()
Proof. We first prove that there exist a sequence of coefficients {axtkelus) satistying Σx= ag = 1 so that:
configuration Up Pas (x/),
W3 Σ anPs (a(xx))
p(z!l Un Pan (x:))
(49)
If this is correct, then p(ag) Un Pan (x))) would be consistent estimator of p(0 Un Pas (x()).
Based on Eq. (46. we have:
p( Un (x)(
(50)
Uar
(51)
E. (((nao)n(()))
E.Ex 142 Pan (t)))
(52)
E, A(UnPan
E, (r)) mI))
(53)
E, E.L(UPPan(tt)
Ereny
(54)
Energ 1(2/ Un (nPa)(T))
Exeng I(UnPan(z())
(55)
Exeny ((UnPandTt)
E, ((UnPan(z))
h
Exemy 1(1% Un (n)))
Exeng 1(U,Pa (xf))
(56)
ExETR ((?))an(x))
((((
k
Exeug ((UnPan(xt))
Σ (Drm (of) Un (v!(
(57)
E, 1(Uh (xt))
Eveny I(U ((UnPan(s))
Σ (Drent UnPan (xt))
(58)
E. PU,,
Exerg I(UnPan(=())
Σ (Premz (x/ (x<(Pa(x())
(59)
E, ((UnPan(27))
Σ (n (ir) Pa(xt)),
(60)
k
Event (an))
where as (T)
(61)
E, I(UnPan

---

we can decompose the numerator and denominator of ag with homogenous time partition until each component converges a stationary distribution.
Using the same logic in Eq. 30
BS
Eien; (UnPan(z())
(T)
(62)
Σ EnETy, 1(U, Pan (112))
T Σ ==1 8 EVER(EN ((U)Pan(z)))
(63)
T Σx ==1 8 EIER/A (an))a)()))
(k))
Σ 8= P (NPPa((z))
lim (T) T-000
(64)
Σ Exs ma PrEF(M.))
(UhPas (x())
Without loss of generality, assume y Pa(xz), where E M and y & Pa(az), where & ITJ Then we have
E, Uha,,,) p(x/1 Un Pan (af( (v)
(65)
Σ, I(un (1)
Exeng ((U))an(T)))))) (w)
Σ (Prem (io Un Pan \.)
(66)
E, [(((unaantt)) (1)
Exemy (x/ UnP 1) Erem ((UnPan \8)
Exeny A(UnPan(c)) (x{) v) ((L((Paan(T)) (1)
Be (T(T (o) (a(x))) + (T) (xf Pa(x)) v))
(67)
2=
Σ Einng P(UnPan (z)) vv)
where (c( =
(68)
Σ, L(UnPan (1/) \2)
Similarly, we have:
Enery I(UnPan(at) (1)
Bx(T) Σn Eneng 1(U)Pan (x7) (1)
(69)
Σmh I 8 EiET/
(UnPan( \4)
(70)
Σn Ein x= EiEn(an P(Uh \ v)
Σ POET(U.,, (UnPan(Zt) v)
lim (T)
(71)
T-000
Σn Eesr PrEsie., (UnPan (x (v)
Proving p(a?) Un Pan (x/1) # p(a?f Un Pan (x!)) \ v) is equal to proving:
U Pan()))) - p( Un Pan Pxt( \ v) + 0
(72)
10

---

Substitutes Eq. 60 and Eq. 67 in Eq 72 we have the following derivation:
p(x.] Un (a)())) p(12l Un Pan(x?) \ v))
(73)
= lim T0000
Σ a (a (o) (r/ (a(Paa)))))
(74)
lim
Σ (p (a(()(()))) + B1 Pren (of (a((()))))) (11)
T-
k=
Σ s=1 PIEN (UnPan(z!))
Preng (x/(P((z)))
(75)
(UnPan
Σ == PrEn(e..)
(UnPan (xt) \y)
Prent (x/ (a((Pa(()))
PrEatu.. (UnPan(at) (n)
PLETT (Un (UnPan(xt) \ v)
Σx EE P (z)) Prenf (o/ (((laa()))) ww)
After equating the denominators, the numerator is:
45
(UnPan (1/ pa(x)))))
[Σ> k=1s=1
(UnPan(x7) \w)
Premiu.))
Preat
1m s=1
wy
PrEF(K.)) Pap (xt) MP (xg) Pa(x)))
My
(UnPan(Tt) \ W (r/()a((t!) (v)
PiEI
My
wj
Pie
(UnPay, (xP))())))
(76)
(k..)
For the sake of simplicity, denote
0k := Σ Prem (k.))
(UnPan(x!)
(77)
bk = Σ PrEm
(78)
#=1
C Premç (a((x))
(79)
d Prent (or (P((z)) \v)
(80)

---

After substituting the simple notations in Eq.76
C&CR)
bk)
back + bidl)
as)
(81)
k=1
2=2
==1
Ck
cr. caab +
cc)aabb
(82)
k=1
k= i> 1,itt
W/
Wy
Wy
Wy
Wy
b
CLUl - cy
ak
CkOK
bi
Ok
G
(83)
k=1
k=1
k=1
k=1
1>1
Define
V (X, < <t}
(84)
That is, V2 contains all the nodes before time point t.
Denote {br, hi]n)] = Par (z() and assume {be, hissh <n (x/),, where E II
(k.s)
We express Prem (Un, Pau (Ir)) by marginalizing all other random variables occurring before the latest variables Pa (xt) and utilizing the Causal Markov assumption (A2):
PrEm (anPan(z))
(8:5)
Pau (x/))t m())
(86)
E
plbe
Vn \ (butiemm)l
E IT (M))))
(87)
Σ
(((e,)a((be.)
Σ 3 Σ
P( FMon
(88)
(Pa(bi)))e(m))
Vimox Tmaxe KCC' Sh af Pa(x))())))
Note that 2 E Vn Jein))
This joint distribution is now represented by conditional distributions of one related variable given its parents.
Similarly, assume y = br,, we have
Prem (U, (an()))) \ (1)
(89)
= POET (K.))
(Un Pan be.)
(90)
p(UupPa (ot w)(t E 1
(91)
= E
bibes....b benn Ven \ (br, m n))
(92)
Ve, Vol, Fige
Σ
Σ Σ Σ
( Vinne)
(93)
(Pat
Vemax Tmax C Cl'Shefl at Pa(x))
Note that oi E Ven \ (be, J
12

---

Prem (a: Pa((x)) \ v) can also be represented by those conditional distributions based on Bayes rule.
Prenf (r](Pa(x))) b
(94)
Prent (et, beg beg Deng)
(95)
Preng Deag)
Edi, PREmT (Ite bege bear)
(96)
Dmn, PrEmt ( be,))
Emi, PrEmT (o7. Par (x?))
(97)
Dm,, PrEm (Pai (x())
Eny p(az [Pai (xt)) (be beng)
(98)
Eny PLEUf (be, .., blm))
Ene,
Eva Enalereetnn Prent (be, Vn (Dess 3 = max(tie(n]]))
Emn Eva Whlleeemm PLEnt
Vn (breetor) H/ = max((se(m.13)
(99)
Ein AB
(100)
Emn CD
where
A = Par(())))
E
p(b,, (Pa(b:)))
(101)
(Pu(te,)hellnl
B = Σ E
Σ
Σ
p( (a(x)) rmax)
(102)
Vimax Tmax <t'
Sh (E() Pa(t))
C= Σ
p(br, (Pa(br.))
(103)
(Pa(bi)))eenn
D Σ E
(104)
Vanax Tiax <V af) Pa(,)
Note that & E II, for distributions in above section from Eq. 9q. to Eq. 104 and that o/. Vh. {b,, (e(m))
Hence, every term in Eq.
be expressed function of those conditional distributions. Substituting Eq. 88) Eq. and Eq. (100) in Eq.833). we have a polynomial equation only composed
of conditional distributions (x Pa(x)))) he except the joint distribution of the starting
points p(V...... Note that the conditional distributions of variables in {X/) {X/heen [n], [wy]]
are the same. Since sets do not allow duplicate values, set (p(=/]Pa(<>) [Pa(=)))) (e contains only different conditional distributions. There should be potentially total Σ wi different causal mechanisms. The total number of conditional probabilities should be jointly determined by the number of causal mechanisms and also the number of realizations that variables can take. After adjusting those conditional distributions by the linear restriction E, p(xly) 1, all components in the set of all the causal mechanisms because the first starting points are random noises. That is, upon adjustments, all the terms in Eq should be rendered independent of each other, without any imposed constraints across them.
are mutually independent, and p(V is also independent
After expanding all the summations in Eq the coefficients of this polynomial equation either or - 1. Each coefficient is accompanied by one unique monomial as index (k, s) in the joint distribution Prexî_ determined unique product of conditional distributions, i.e., with different pair of (k, s). the product should be different. Considering all random and independent conditional
13

---

distributions in [p(x//|Pa(x)) hie fjefn),'s. the polynomial is not identically zero, and the probability of choosing a root of this polynomial is zero.
Denote the polynomial equation in (83) as A, we have:
p = =0
(105)
Back to the original E3.773. we finally have p p(x//1 Un Par (x()) & p(æf1 Un Pan(2)) \w)
Vy UaPan (x:).
Lemma D.3. Let SPa(X.) denote the estimated superset of parent set for X E V obtained from the Algorithm B/ (line 2). (Pax(X))) contain the true and illusory parent sets, where wj is
the true periodicity of X3. Under assumptions AI-A7 and with an oracle (infinite sample size limit),
we have:
1 Pax S SPa(X)), W [mmax 1,2]
almost surely.
Proof. Assume the contrary, i.e., there exists 8 E U&Pai(Xy)) \ (SP((X)) From Lemma D.2 have Xi 1 Pax \ 8. By the Definition 2.4, have Pa(X/) U (X/).
& Pa(X))), by the causal Markov property (A2), the dependence relation can not be true, because is a non-descendant of X/. If E Pa(X)), our Algorithm would have concluded that < SPa(X;) (line 2) with a consistent CI test, evident from the causal Markov property, contradicting our assumption. Hence, the lemma.
Lemma D.4. Let Pa(X)) denote the estimated parent set for X E V obtained the Algorithm Bl (line 19) assuming that true wy has obtained (line 17). (Pak (Pax(H)) (X/)kke(ws) ke(wsy] contain the true and illusory parent sets. Under assumptions Al-A7 and with an oracle (infinite sample size limit), we have:
Fa(x:) = Pa(X)), Vt [max 1,
(106)
almost surely.
Proof. From Lemma
D.3
Pa(X)) Pax())) S SPa(X))) Vt [maa + 1,T7 E [m]]
(107)
In Runge al. 2019]. the author proved Pa(X Pa(X;) if we run PCMCI stationary time series. Using the same logic, we have the following proof.
Suppose X & Pa(X)) but X E Pa(X)). With consistent conditional independence test and correct time partition, the MCT test (line 8 in Algorithm B1 will remove X; from Paw, (X?) and only
x 1L x SPa(X)) \ (
(108)
Based on Eq. (107) the rule is equivalent to removing X from Pan, (X/) if and only if:
Pa(x)) X (axx))
t_ 11 X/
SPa(x
() \Pa(ii
(109)
*Xi_ 11 X! Pa(X)) \ XiL.).P(Xii.
(110)
14

---

Based on Causal Markov Condition assumption (A2) and Faithfulness Condition (A3), from Eq. have Xi & Pa(X))- In other words, if X & Pa(()) then X & Pa(Xi). That is, Pa( (2) _ Pa(X?)
Suppose Xi_, E Pa(X)) but _ $ Pa(X/). By the contraposition of Faithfulness (A1), we know that Xx X Pa(X)) (Xi_r), Pa( (,,) Denote W = (Pu(X)) \ (Pa(X)), } U (Pa(Xi-)) \ Pa(Xi_))) Since Xi & Pa(X!), based on Causal Markov Condition assumption (A2),
WUX 11 X Pa(x?)
WUx 11 x Pa(x)), Pa(xi_)
Weak Union Xi 11 xi (Pa(X))) Pa(xi=))UW
X 11 X Pa(Xt) Pu(x))
This contrary to the assumption so that there is no such Xi_, satisfying Xi E Pa(X?) but Xi_, & Pa(X)). In other words, if Xi-r Pa(X)), then Xi_, Pa((X). That is, Pa(X)) C Pa(X?). Combined with the previous conclusion that Pa(X?) S Pa(X?), we have Pa(X)) Pa(X;).
Based on Lemma D.2 Lemma D.3 and Lemma Lemma
D.4 we can identify the true wy for X through
Lemma D.5. Let wij denote the true periodicity for X E E and Pa(x reng
Yeug) denote the estimated
parent set for X2 obtained from Algorithm
B1 where t E IP Define:
W5 arg min max [Pa(X ( wE(wal kelw)
(111)
Under assumptions A1-A7 and with an oracle (infinite sample limit), we have that ر = Wjs Vj E [M] almost surely.
Proof. Assume the contrary that Wj * wy. then in the Algorithm we have an incorrect time partition ft. Hence, CI tests that are performed use samples with different causal mechanisms. P(X - (ax(X)) (X?)) ) in 50. is estimated from mixture of two more time partition subsets, say II and IPz. We can apply Lemma D.2 where Pax is replaced by _2__ Pa (X/) and then in Lemma D.3 SPa(X;) is replaced by Pas, (X3) and hence _ Pax S Pas, (x))
where Pas, (X/) is obtained from samples with é from the mixture of two different partition subsets (line 8). Hence, with مj. |Pas, (XP) 2 1U Pak (Xx3) using mixture samples U²_ IT! However, with true Wy, we have Paw, (X/) ) = [Pa(X/)] based on Lemma D.4 With Assumption
A6 the Hard Mechanism Change, U = Pag (X()) > Pa(X?)] so that wj always leads to smaller size of estimated parent sets than Z contrary to the definition of wy. Hence, W w/.
With those lemmas, we can prove Theorem 1.
Proof. Assuming that a correct wy has already been obtained, from Lemma D.4 we have
Pa(X;) Pa(X)), Vi 2 27.h. E [n]]
From Lemma D.5
know that correct w; must be obtained with consistent CI tests, that is.
4y = Wj, Vj E n Therefore from Algorithm B we have
Pa(x?) Pa(X?), Vt 2 2T j E [n]
15

---

[Pa, (x22(x
[Pag((51
@ 1
30
0 2
26
21
>
3
11
10
0
21
12
15
9
&
13
Figure 3: In the above illustration of the "turning point," the sizes of parent sets for different estimates Wy depicted as [Xa)(x))], parent set or the illusory parent set of X In this context, we are interested the sizes of these parent sets. The first occurrence of the "turning point" happens Wj 3 since the sizes of parent sets obtained when W3 and w = are larger than the corresponding size when W = 3, respectively, The term "turning point" denotes that as Wy increases, the size of the parent set initially decreases and then starts increasing once the local minimum reached. The corresponding relations exist because as long Wj is not multiple of the true wij, the estimated time partition subsets with wi must be a mixture of some correct time partition subsets with wj. Therefore, is reasonable to use this trick rather than looking at the maximum size of the parent sets Pa,(X/) k E [(y]] (line 17 in Algorithm BTp.
[.a]] It is worth noting that Pak (X²( represents either the true
If the causal mechanism is fixed across time, i.e., Wy 1,j E [n], the proof of PCMCI Runge et al. 2019 showed that for all X E V.
X + X p 6= x X く X + xX و x + X و
Therefore 9 G.
3w, > we can simply separate the whole graph into sub graphs (9% heewe,) consisting of only target variable X? with corresponding ( and parent variables X e Pa(X/). Focusing only on one time partition subset IT. E [Wr]]. we have
gh 9%
(112)
any [wy]] and E [n]] based on the proof of Proposition 1 in the supplementary materials of
Runge et al. 2019
Each sub-graph 9% includes only variable X the edges entering X for time points / E IT and the
corresponding parent variables X Pa(X)). Given IP'
IT and V X we have:
kE[wy]
jE(n)
g
(113)
jEin], KE|wy)
U g jEln|, KElwyl
guus
(114)
On the basis of Eq. 112 we finally have:
G
E
Turning Points
Given infinite samples, our estimate 4; (line 17 Algorithm Bi
the exact value w, (see Lemma D.5 However, for finite samples, estimating wy by the equation in line 17 in Algorithm Bl does
16

---

not yield good performance when T is small. While searching, larger guesses w lead to finer time partitions in I1,. resulting in smaller sizes for II (see Line in Algorithm B1 Due the power limit of CI tests on smaller sample given by 11 the number of false negative edges increases. In order solve this issue, introduce turning points. A turning point is guess i satisfying:
max Pag (X()1 < min (Pas-1 max (Pa)+1(X(11
where PPa (X²) is the estimated parent set for X/ with periodicity guess W. See line 19 in Algorithm
We illustrate it with a special example in Fig If there are several turning points, then ش, is the first turning point. If there is no turning point, then we obtain Wj using Line of Algorithm BT
The concept of the turning point is not based on any formal theorem but rather on experimental observations. In experiments, the turning point often corresponds when 7 is not large. This occurs due to the limitations of CI tests on finite samples. In such cases, the causal graph can still be correct because the estimated time partition remains accurate. these experiments, the accuracy rate is calculated by considering (Nay NE, m
multiple of the periodicity
as correct estimations.
F
Computational Complexity
Executing the PCMCI algorithm on the entire time series constitutes the initial phase of the proposed
approach (Algorithm BT line 2). The algorithm's worst-case overall computational complexity is O(n ab + On Tub). discussed in Runge et al. [2019) Here, the symbol n denotes n n-variate time series and Tub represents the upper boundary for time lags.
The subsequent computational load stemming from the remaining components of our algorithm follows complexity of Ouw "(a)) This encompasses the O(n²Tab) complexity associated with conducting Momentary Conditional Independence (MCI) tests on all n univariate time series. The parameter wib here arises due to the search procedure involving w, iterating through values from 1 to Web for all n univariate time series.
The runtime of the computation is further influenced by the scaling behavior of the CI test concerning the dimensionality of the conditioning set and the temporal series length T. For further details, see section 5.1 in Runge et al. 2019
G Experiments
All experiments, including those detailed in the main paper, are conducted on single node with one core, utilizing 512 GB of memory in the Gilbreth cluster at Purdue University.
Here, we describe how to calculate the metrics (Fi score, Adjacency Precision, and Adjacency Recall) in our setting. In stationary time series, the output of the causal discovery algorithm is typically an adjacency matrix with dimensions [n, n, Tmax + 1]. Within the three-dimensional binary array, edge pointing from one variable another with specific time lag, while 0
the value 1 signifies indicates the absence of an edge. For instance, if element [ k] in the matrix is 1, then there is an edge pointing from X; to X². In semi cemi-stationary time series, due to the presence of multiple causal mechanisms, the binary edge matrix is four- oour-dimensional array with dimensions [n. Ω. n, Tmax + 1]. where Q is defined as (7) in the main paper. This expanded binary matrix is constructed based on the edge matrix of each variable X E [n], through repetition. For instance, if S2 2wj, setting the third dimension of the large binary matrix to j should yield Wj potentially different parent sets (including illusory and true parent sets), each appearing twice.
We should have two such binary arrays, one representing the ground truth with dimensions
[n, 22 n. Tmax + 1] and one obtained from the algorithm with dimensions [nz ②2 Tmax + 1). If the estimator Ω is incorrect, those two binary arrays will have different sizes, compare them. To solve this problem, we do the same operation and calculate the least common multi- ple of Ω and Ω. Denoting this least common multiple as LCM(9, we create two four- niu-ddimensionll binary arrays with dimensions of [r, LCM(S, Q), Tmax 1] based on the true edge array and the
not directly
17

---

estimated edge array, respectively, through repetition. The metrics are then computed by comparing the values in these two arrays.
G.1 More Discussion regarding the Case Study
As stated in the main paper, we express our inability to comment on the significance of the case study results. We open door for the related experts; if assumptions A1-A7 are satisfied, the stationary assumption may not hold in this real-world dataset, and such periodicity exists. However, the finding not correct from expert's 's viewpoint, the following assumptions may be violated:
Assumption A4 No Contemporaneous Causal Effects: There is a possibility of potential causal effects from X to X that the algorithm is unable to capture.
Assumption A6 Hard Mechanism Change combined with limited power CI tests: If there is soft mechanism change in the variables, the reliability of the CI test of variables given their parents will be influenced by the skewed distribution of the parent variables. This effect will be exacerbated by the fact that the sample size will be shrunk by w.
We provide sound and robust algorithm for experts in various fields who interested in validating the presence of periodicity within the causal mechanisms specific to their domain.
G.2 Experiments on Continuous-valued Time Series with Exponential Noise
Considering that VARLiNGAM a temporal extension of LiNGAM and LiNGAM is an algorithm designed for non-Gaussian data, following the work in Pamfil et al. 2020 also construct experiments on contnuuus--vlluee time series data with Exponential noise. Shown as Fig [a] the performance of PCMCIa, PCMCI and VARLiNGAM, are quite similar with their performance on Gaussian noise. The recall rate of DYNOTEARS, however, gets worse with Exponential noise.
G.3 Experiments Binary Time Series
Similar to the process of generating continuous-valuee time series, the generation of binary time series also involves three steps. However, the main difference lies the last two steps. the third step, we simulate the conditional distributions of each child variable based on all possible combinations of parent variable values. Subsequently, we randomly generate the value of the child variable by considering the corresponding conditional distribution given its parent sets.
For discrete-valie valued time series, longer time length is required. To evaluate performance, conduct a series of experiments following the same methodology as described in section 4.1. Fig 4bb) illustrates the variation in comprehensive performance with respect Wmax. PCMCI demonstrates similar performance to PCMCI in terms of the F'l score, indicating well-balanced trade-off between precision and recall. This outcome expected since discrete-valiee
time series demand larger sample sizes, and the increases in Wmax negatively impact the power of MCI tests. This observation is further supported by Fig ((a] where an increase time length T from 4000 to 12000 does not lead significant improvement in the accuracy rate of W, while the accuracy decreases rapidly with higher values of winax-
Comparing these results to the experiments conducted on continuous-valuud time series, becomes evident that the demand for efficient samples is even more substantial for binary time series, and the influence of increasing Wmax on performance becomes more pronounced.
Fig ( shows how the performance of the algorithm varies across Tmax and the same trade-off between recall and precision has been shown.
G.4 More experiments on Continuous-valuu time series
In this section, we conduct more experiments with continuous- time series with Gaussian noises.
In Fig 6(a) we test our algorithm with and without utilizing the furning point rule. See lines 13-14 in Algorithm B] and section E for more information about the turning point rule. Let PCMCIg TP
18

---

Tmur
n m,Tmac
1.0
1.0
T 4000
T 8000
12000
0.8
0.8
0.6
Score
0.6
0.4
0.4
0.2
0.2
PCMC PCMC
0.0
0.0
3
7
9
3
5
7
9
1
3
9
1.0
1.0
0.8
0.8
0.6
0.4
0.4
0.2
0.0
0.0
3 7
3
3
2
3
5
1
2
3
4
5
1
2
3
4 5
1.0
0.8
0.8
0.6
0.6
0.4
PCMCI
0.4
0.2
POMCI
3
5
DYNOTEARS
VARLINGAM
4
(a) Exponential noise
(b) Binary time series
Figure 4: F) Score, Adjacency Precision, and Adjacency Recall when Wmax varies for experiments on continuous-val valuue time series with Exponential noise, length T {500, 2000,8000), Tmax and n = 5. b) F1 Score, Adjacency Precision, and Adjacency Recall when Wmax varies for experiments binary time series with length T {4000, 8000, 12000], Tmax and 3.
I1 3, Umax
T 4000, Omex
1.0
1.0
I
T=4000
F, Score
+
T=8000
0.5
0.8
+
T=12000
0.0
2
Accuracy
0.6
1.0
Precision
0.5
0.4
0.0
2
4
5
0.2
1.0
PCMC
Recall
0.5
0.0
2
4
0.0
( Wmax
2
(a) Accuracy of W
(b) Performance when Tiax varies
Figure 5: PCMCIp is tested 3-variate binary time series. Every marker corresponds to the average accuracy rate or average running time over 100 trials. a) The accuracy rate of á for different time series lengths and different Wmax. b) F1 Score, Adjacency Precision, and Adjacency Recall when Tmax varies for experiments with time series length T 4000, Wmax = and n 3.
denote the version of PCMCIg that the turning point rule is utilized in choosing w. PCMCIg non-TP TP means that the turning point rule is not applied and w is chosen directly according to Lemma D.5
Fig 6(a) shows that the algorithm PCMCI, non-TP and PCMCIa have similar performance with various T and wmax- With T = 500, PCMCIa non-TP yields slightly larger standard errors for those metrics, compared to PCMCIa TP. As time length 7 increases, the performance of the
algorithm PCMCIp non-TP TP has consistently increased and is even slightly better than PCMCIa TP.
19

---

The consistent performance of PCMCIg under different chosen rules of w supports our theoretical result; that is, the correct periodicity leads the most sparse causal graph.
In Fig 6(b) non-stationary time series are produced instead of semi-stationary ones. Consequently, the causal mechanisms for each univariate time series longer appear sequentially and periodically, The proposed method performs slightly better in terms of Fl score and precision. However, the recall rate is the worst compared other baselines.
Fig 6(c) we conduct experiments in the nonlinear setting. The proposed algorithms PCMCIg TP and PCMCIa non-TP perform the best.
In Fig F(d) with wub Wmax- the performance of the proposed algorithm is significantly worse compared to the scenario where wab > Wmax. However, with < Wmaxx the proposed algorithm can still detect a less dense graph in comparison other baselines. Based on these outcomes, is essential to maintain slightly higher wab without significantly impacting the number of efficient samples utilized in each CI test.
20

---

n 5,Tra -5
S,Tmax m
500
1.0
1.0
0.8
0.6
0.4
0.4
0.2
0.0
0.0
9
3
5
9
1.0
Srecisio
0.8
0.6
0.6
0.2
0.2
13579
7
1.0
0.8
0.8
0.4
0.4
PCMC
0.2
PCMCI TP PCMCI non-JP PCMC/
0.2
PCMCIa non-TP
5
PCMCY
Ginxe
VARLWSAN
VARLWGAM
(a) Performance with and without turning point (b) Performance in non-stationary setting with- out periodicity
Trox - 2000
T.ae
T 500
4000
500
2000
8000
1.0
0.8
0.8
0.6
0.4
0.2
0.0
3 5 7 91 3
3
1
1.0
Precision
0.8
Prectoie
W
0.6
0.4
0.4
0.2
0.0
3
5
7 3 1357
1.0
0.8
Recait
0.6
0.6
0.4
0.4
PCMCI, 4hp Wleer
0.2
PCMCI
RCMCI,
PCMC nop-TP
0.0
PCMCI
PCMCI
tsFC/
DYNOTEARS
(c) Performance in non-linear SCM
(d) Performance when With Wmaax: and With
Wmar
Figure 6: Multiple algorithms are tested 5-variate time series with different time lengths T. Every line corresponds to a different algorithm. Every marker corresponds to the average performance over 50 trials. In (a), the consistent performance of PCMCI under different chosen rules of supports theoretical result; that is, the correct periodicity w leads to the most sparse causal graph. (b), data sets are in non-stationary setting without periodicity. In (c), the structural causal model (SCM) is non-linear. (d), algorithm PCMCIg are tested under conditions that Wab > Wmax and Wab < Wmax respectively.
21

---

References
Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science advances, (1))aaa44999, 2019.
Thue Duy Le, Tao Hoang, Jiuyong Li, Lin Liu, Huawen Liu, and Shu Hu. fast pe algorithm for high dimensional causal discovery with multi-core pes, IEEE/ACM transactions on computational biology and bioinformatics, 16(5):1483-1495, 2016.
Dimitri Bertsekas and John N Tsitsiklis. Introduction to probability, volume 1. Athena Scientific, 2008.
Samuel Karlin. A first course in stochastic processes. Academic press, 2014.
Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Geor- gatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, pages 1595-1605 PMLR, 2020.
22


Reviews:
[{"review_id": "aAEegbt218", "reviewer": "Reviewer_nsZ6", "summary": "The authors propose a causal discovery algorithm for a class of periodic random processes called semi-stationary process. The algorithm is built upon an existing algorithm assuming stationarity. Consistency results of the algorithm is provided. ", "questions": "Since PCMCI assumes stationarity, in Algorithm 1, why PCMCI is used to estimate the superset $\\widehat{SPa}(X_{t}^{j})$ for semi-stationary processes?", "limitations": "No potential negative societal impact.\n\n", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. Consistency results of the algorithm was shown even though the algorithm is relatively complicated.\n\n2. The algorithm shows superior performance against the baselines.   \n", "weaknesses": "1. The motivation for generalizing PCMCI to handle semi-stationary processes is not clear, especially from a theoretical perspective. The description of PCMCI is completely moved to the supplementary. I think including a brief description and the key procedures (e.g., how the superset is estimated) is helpful for understanding why PCMCI is good fit for the considered setting. \n\n2. As mentioned in Section F. from Appendix, the proposed algorithm does not yield good performance in finite samples until the heuristic method $turning\\text{ } points$ is introduced. To support the theoretical results, it should be shown that the algorithm starts to work when the sample size is sufficiently large. \n\ntypo: \nDefinition 2.3 is not correct: $\\Pi_{k}^{j}(T)$ depends on $k$, but its definition in (8) does not depend on $k$.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We really appreciate the feedback. Regarding your questions,\n\n## Q1. Description of PCMCI\n\nThank you for the suggestion. In the final version of the paper, we will provide a concise overview of how PCMCI functions, utilizing the additional page available in the camera-ready version.\n\nHere is a brief overview of PCMCI. There are two stages of PCMCI: the condition-selection stage and the causal discovery stage. In the first stage, unnecessary edges are removed based on the conditional independencies from an initialized partially connected graph. In the second stage, Momentary Conditional Independence tests (MCI) are used to further remove the false positive edges caused by autocorrelations in data.\n\n## Q2. Relevance of PCMCI\n\nAfter running PCMCI on the whole multi-variate time series with Semi-Stationary SCM, the parent sets  $\\widehat{SPA}(X^{j}_{t}),j\\in [n]$ we obtained (line 2 in **Algorithm 1**) for each variable should be a superset of its true parent set for any time point $t\\in T$, proved by **Lemma 3.2-3.3**.\n\nFor a guess of $\\hat{\\omega}_{j}$, we can construct a corresponding series of time partition subsets. A series of MCI tests are conducted on samples whose time points are from the same time partition subset.\n\nBased on the result of the MCI tests, variables in $\\widehat{SPA}(X^{j}_{t})$ will be removed based on the conditional independencies.\n\nIf $\\hat\\omega_{j}=\\omega_{j}$, the samples we used in each MCI tests are from the same causal mechanism and $\\widehat{SPA}(X^{j}_{t})$ will shrink to the true parent set under assumption **A1-A7** with an oracle (infinite sample size limit).\n\nIf $\\hat\\omega_{j}\\neq\\omega_{j}$, however, the samples in MCI tests are from different causal mechanisms, bringing in more causal relations under assumption  **A6**. Hence according to the MCI test results, $\\widehat{SPA}(X^{j}_{t})$ will shrink (or not shrink) to a wrong parent set whose size is larger than the true parent set, leading to a denser causal graph, proved by **Lemma 3.4**. \n\nHence, guess $\\hat \\omega_{j}$ resulting in the most sparse causal graph equal to the correct estimation $ \\omega_{j}$. Based on the consistency of PCMCI algorithm and our proof of **Lemma 3.2-3.4**, the proposed algorithm PCMCI$_\\Omega$ is sound under assumption **A1-A7** with an oracle (infinite sample size limit).\n\n## Q3. PCMCI$_\\Omega$ without *turning point*\nThank you for this insightful comment, and we will add the related discussion in the camera-ready paper.\n\nWe have conducted new experiments without using *turning point* and therefore $\\hat\\omega_{j}$ has been chosen according to the original rule stated in line 12 in **Algorithm 1**. The results show that algorithm shows similar performance with and without the *turning point*. \n\nShown as **Fig.1(a)** in the Rebuttal pdf file, PCMCI$_\\Omega$ without *turning point* results in a slightly larger standard error with a smaller sample size. \n\nAs time length $T$ increases, the performance of the algorithm without *turning point* has consistently increased and is even slightly better than PCMCI$_\\Omega$ with *turning point*. \n\nThe consistent performance of PCMCI$_\\Omega$ under different chosen rules of $\\omega$ supports our theoretical result, that is, the correct periodicity $\\omega$ leads to the most sparse causal graph.\n\n## Q4. Updated Definition 2.3 *Time Partition*\nThank you for pointing this out. This definition in the main paper is not precise. The updated definition should be as follows.\n\nA time partition $\\Pi^{j}(T)$ of a univariate time series $X^{j}$ in a Semi-Stationary SCM with periodicity $\\omega_{j}$ is a way of dividing all time points $t\\in[T]$ into a collection of non-overlapping non-empty subsets $\\Pi^{j}_{k} (T),k\\in [\\omega_j]$ such that\n\n$\\Pi^j_k(T):=\\\\{t:\\tau_{\\max}+1 \\leq t \\leq T, (t \\bmod \\omega_j)+1=k \\\\} $\n\nwhere $\\bmod$ denotes the modulo operation. For instance, $5\\bmod3=2$."}, {"review_id": "UuWEh7yA2v", "reviewer": "Reviewer_YnE6", "summary": "## Summary\nThe paper addresses the problem of discovering causal relations in semi-stationary time series data, where a finite number of causal mechanisms occur sequentially and periodically across time. The authors propose a constraint-based, non-parametric algorithm called $PCMCI_\\Omega$ to handle this type of time series data. \n\nThe main contributions are \n1. The proposed algorithm, PCMCIΩ, is an extension of the PCMCI algorithm designed for stationary time series data. It leverages the PCMCI algorithm to systematically discover the superset \nand use CI test to find out the correct periodicity and remove unnecessary parents.\n\n2. The soundness of the PCMCIΩ algorithm is demonstrated through a theoretical analysis, showing that it can recover the true causal graph under specific assumptions.\n\n3. The authors validate the effectiveness of the PCMCIΩ algorithm through experiments on both continuous and discrete synthetic time series data. Although I think the baselines can be stronger.\n\n", "questions": "1. I am a bit confused on line 123, why minimum periodicity has the same symbol as the number of causal mechanisms in $V$?\n\n2. In definition 2.3, eq.8, the subscript $k$ is not showing in the latter equations. What does this k means?\n\n3. Can you make the definition of Illusory parent sets easier to understand and explain the intuition behind it?\n\n4. For figure 1, why $Z_2^1$ is identical to $Z_1^7$? I thought $Z_2$ defines the Markov chain for second time series? Also, you use superscript in $X$ to indicate time series number, I think $Z$ should follow the same pattern to avoid confusion. \n", "limitations": "The author does not explicitly provide a discussion on the limitations. Hence, I suggest the author to include a limitation section. ", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Strengths of the paper:\n1. Originality and significance: The paper proposes a novel algorithm, PCMCIΩ, to address the challenging problem of causal discovery in semi-stationary time series data. Fewer work have considered such settings and most of the previous method assume stationarity. This is a step towards generalising causal discovery to non-stationary time series. \n\n2. Soundness: This work provides a theoretical guarantees on the soundness of the proposed method, and validate it with experiments.\n", "weaknesses": "The main weakness of this paper are\n 1. Limited scope of experiments: Although the authors have performed experiments on synthetic continuous and discrete time series data, as well as a real-world climate dataset. However, the baseline selection is too simple. For example, VARLinGaM and DYNOTEARS are both linear models with stationarity assumptions. It is expected to outperform those two methods. I suggests including more stronger baselines, such as Granger causality baselines [1], and SCM-based method [2]. \n\n2. The paper is hard to understand when I first read it. There are many definitions defined in section 2, without knowing why they are important. The graphical illustration helps a little bit, but still confuses me on why we have to define them. For example, I am not sure I understand why Illusory Parent sets requires a subscript. Like why we need $Pa_2(X_7^1)$ and $Pa_3(X_7^1)$. Isn't it the same as $Pa_1(X_8^1)$? Since the following theoretical analysis extensively uses these definitions, I recommend the author to be very clear on these definitions and should make them easier to understand for the readers. \n\n3. How does your method performs when there are mismatches in the assumptions? For example, the ground-truth is stationary/non-stationary. When the lag upper bound and period upper bound are shorter than the truth. etc. More ablation study would be helpful.\n\n4. What are the potential limitations of the proposed method? Some discussion should be added. For example, instantaneous effect can be important when aggregation happens. Computational complexity with dimensionality, etc. \n\n[1] Khanna, S., & Tan, V. Y. (2019). Economy statistical recurrent units for inferring nonlinear granger causality. arXiv preprint arXiv:1911.09879.\n\n[2] Gong, W., Jennings, J., Zhang, C., & Pawlowski, N. (2022). Rhino: Deep causal temporal relationship learning with history-dependent noise. arXiv preprint arXiv:2210.14706.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Your insights are greatly appreciated as we work towards enhancing the robustness of our method. Regarding your questions,\n\n## Q1. Scope of experiments\n\nTo our best knowledge, there is no such non-parametric algorithm that can handle the periodicity in causal mechanisms and discover the time-lag effect with a window causal graph. \n\nFor the Granger causality baselines [1], the estimated causal graph is a summary causal graph, while the output of our proposed algorithm is an estimated window causal graph. They highlight distinct facets. Summary graphs are less informative than the window graphs about the underlying causal structure in our case.\n\nFor the SCM-based method [2], we have contacted the authors and will conduct the experiments as soon as we get the code. \n\nBased on your valuable suggestions, we conduct more experiments in the nonlinear setting with the baseline method \"tsFCI.\" In the new experiments, shown as **Fig.1(b)** in the Rebuttal pdf file, the SCM are non-linear and the proposed algorithm performs well.\n\nAs for the \"non-stationary\" assumption, to showcase our method's capacity to handle periodicity in causal mechanisms with time-lag effect, a capability lacking in other approaches, we will introduce a related baseline called \"Regime-PCMCI\" in our new experiment. It's worth noting that this algorithm necessitates a linear model and additionally mandates precise foreknowledge of the number of causal mechanisms inherent within the time series. Remarkably, based on several experiment results we have, our algorithm's superior performance is sustained even when the accurate count of causal mechanisms is supplied to the \"Regime-PCMCI\" method. Given the considerable time requirement associated with the \"Regime-PCMCI\" method, we regret to say that the result plot is currently unavailable. However, we intend to incorporate the said plot into the paper at a later stage.\n\n## Q2. Visualization of definitions\n\nWe appreciate your suggestion, and we will incorporate additional figures that illustrate each definition separately into the supplementary material. \n\n## Q3. Ablation study\nWe have conducted new experiments for the ablation study. \n\nOur proposed algorithm can handle stationary SCM as stationary SCM is a special case of semi-stationary SCM where $\\Omega=1$.\n\nIn a non-stationary setting without periodicity, shown as **Fig.1(c)** in the Rebuttal pdf file, the proposed method performs slightly better in terms of F1 score and precision. However, the recall rate is the worst. \n\nIf $\\omega_{ub}<\\omega_{max}$, the performance of the proposed algorithm decreases, but still can detect a sparser graph, shown as **Fig.1(d)** in the pdf file.\n\n## Q4. Limitations\n\n**Assumption A4** No Contemporaneous Causal Effects is needed in the underlying Markov chains of the time series. See line 39 in the Appendix. \n\nThe worst case complexity of our algorithm is $O(n^{3}\\tau_{ub}^{2})$+$O(\\omega_{ub}^2n^{2}\\tau_{ub})$. The runtime of the computation is further influenced by the dimensionality of the conditioning set in CI tests. This set's maximal size is $2+|Pa(X^j_t)|+|Pa(X^i_{t-\\tau})|$, coupled with the temporal series length $T$. Refer to section 5.1 in the work by Runge et al.[2019] for more details.\n\n## Q5. Minimum Periodicity is the number of causal mechanisms.\nAs defined in Eq.4 in **Definition 2.2**, the causal mechanisms of $\\mathbf{X}^{j}$ follow a sequential and periodic pattern, occurring every $\\omega$ time points. The smallest value of $\\omega$ that satisfies this condition is termed the periodicity of $X^{j}$, which necessitates $\\omega$ causal mechanisms.\n\nTo illustrate, considering each number as representing a distinct causal mechanism, the causal mechanism progression of $\\mathbf{X}^{j}$ over time could be [1,2,3,1,2,3,1,2,3,...]. In this instance, $\\mathbf{X}^{j}$ comprises 3 distinct causal mechanisms, which cyclically repeat every 3,or 6, or 9,..., or $3N,N\\in \\mathbf{N}^{+}$ time points. The smallest value of $3N$ here is 3, which is the number of causal mechanisms.\n\nHowever, it's important to note that these $\\omega$ causal mechanisms need not to be entirely unique. For instance, the causal mechanism progression of $\\mathbf{X}^{j}$ over time could be [1,1,2,1,1,2,1,1,2,...]. In this scenario, the algorithm treats the sequence [1,1,2] as a single unit, effectively resulting in a periodicity of 3 while containing only 2 distinct causal mechanisms. To maintain simplicity, we still refer to this as having 3 causal mechanisms.\n## Q6. Updated Definition 2.3 *Time Partition*\nThe updated definition should be as follows.\n\nA time partition $\\Pi^{j}(T)$ of a univariate time series $X^{j}$ in a Semi-Stationary SCM with periodicity $\\omega_{j}$ is a way of dividing all time points $t\\in[T]$ into a collection of non-overlapping non-empty subsets $\\Pi^{j}_{k} (T),k\\in [\\omega_j]$ such that\n\n$\\Pi^j_k(T):=\\\\{t:\\tau_{\\max}+1 \\leq t \\leq T, (t \\bmod \\omega_j)+1=k \\\\} $\n\nwhere $\\bmod$ denotes the modulo operation. For instance, $5\\bmod3=2$.\n\nIn this context, by gathering all time points $t$ where the corresponding variable $X^j_t$ share the same causal mechanism, we can form $\\omega_{j}$ distinct subsets of time points, denoted as $\\Pi^j_k(T)$, where $k$ ranges from 1 to $\\omega_{j}$. \n\n## Q7. Markov chains\nThe Discrete-time Markov Chains $\\\\{Z^q_n\\\\}$, where $q\\in [\\delta]$ denotes the index of Markov Chains and $n$ denotes the state index of each Markov Chain, are defined for the entire multivariate time series $V$ and are not defined for individual univariate components $\\textbf{X}^{j}\\in V$. When a specific value of $q\\in [\\delta]$ is chosen, the sequence $\\\\{Z^q_n\\\\}_\\{n\\in \\mathcal{N}\\}$ in **Definition 2.5** constitutes a Markov Chain, with $n$ serving as the progressing index. \n\nIn Fig.1, there is no $Z^7_1$ because there are only six Markov Chains associated with $V$, and $Z^1_2$ denotes the second state within the first Markov Chain."}, {"review_id": "SX6egxrzy9", "reviewer": "Reviewer_ZpAB", "summary": "The paper describes a constraint-based causal discovery method for semi-stationary SCMs, where  the causal structure is periodically changing over time.\nThis type of structure is relevant in real-life scenarios where seasonal or diurnal variation is present for example.\nThe problem is rigorously formalized and and algorithm is give to reconstruct the SCM. As a first step the algorithm identifies a superset of parents using PCMCI naively, \nthen searches for the correct periodicity for all variables by minimizing parent set size. The author provides theoretical guarantees that the reconstructed causal graph \nis indeed the generating graph in the infinite sample size limit, if all assumptions hold.", "questions": "l257: **Following the previous work in Huang et al. [2020], F1 score, Adjacency Precision, and Adjacency Recall are used to measure the performance**\nPlease specify in details how you compute these metrics. Classical methods reconstruct a single causal mechanism per variable, while the ground truth has multiple.\n\nl259: **The standard error of the averaged statistics, displayed either by color filling or by error bars is usually too small to be observed.**\nYou mean here the standard error of the mean estimate? You can use the sample variance instead.\n\nFig 2a, In horizontal axis $\\omega_{max}$ is shown. Has the code been executed with the same $\\omega_{ub}$ for all cases, or the $\\omega_{ub}$ was smaller for smaller $\\omega_{max}$ trials?\nIn general what hyper parameters (upper bounds, significance thresholds) were selected? What CI test was used?\n\nIn the Case Study, if you have monthly data, and Central Pacific has 3 parent sets, it seems to me the causal effect is different in every month say the time partitions are {January,April, July, .. } {February,May, August, ..} {March, June, September, ...}, how you come to the conclusion that: **\"the causal effect from the tropical Atlantic air temperature to the Central Pacific air temperature would disappear every quarter of a year\"*** ? I am not sure I understand this point.\n\nMinor comment:\nI would suggest to use \"discrete valued\" and \"continuous valued\" terminology instead of \"discrete time series\" / \"continuous time series\", as it can be mistaken with discrete and continuous time.\n\nHow the  sampling frequency influences the reconstruction of the SCM? (See weaknesses part.)\n\n**Rebuttal**: I have read the author's rebuttal and we exchanged further comments. I still see the paper having potential interest for the audience of this conference. I had one main open question, the result of the case study.  Its result is unfortunately quite implausible, therefore I cannot further increase my score.", "limitations": "The author points out some limitation of the method. For example the higher sample complexity of the method relative to vanilla PCMCI.\nIt is not clear how realistic the present semi-stationarity assumption.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper describes a novel extension of PCMCI for semi-stationary SCMs, and gives theoretical guarantees. The probrem formulation is sufficiently rigorous, and the algorithm is described in sufficient details. This makes the works reproducible.\n\nThe method relaxes a quite frequent assumption of causal discovery methods in the temporal domain: stationarity. Even if this new relaxed assumption also likely violated by many real life dataset, it is clearly an improvement  in cases when stationarity is clearly a wrong assumption.", "weaknesses": "The main limitation I see is the following: For example, if the causal effect changes by the day of the week but the samples are collected in every hour, this periodicity is hard to represent in the presented framework. Given the process inherent timescale is hourly, we cannot just subsample the data as this would induce transitive dependencies. However in the current form the method will try to fit a new causal mechanism in every hour (until the period is reached of course).\nIt would be good to know what happens if this type of data was fed to the algorithm, as this is hard to avoid. To avoid the identification of the transitive dependencies as genuine we are motivated to work as fine timescale as available.\n\nSimilar question applies for the \"sudden\" mechanism change. In real life most likely the sudden change will transform to a smooth change as the sampling frequency increases (for example seasonal effect in a daily scale). Most likely the effect of the Atlantic air temperature on the Central Pacific air temperature does not disappears suddenly.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We really appreciate the feedback. Regarding your questions, please see below for our responses.\n\n## Q1. Sampling issue\n\nThank you for this great question. \n\nOur proposed method can handle the \"unchanged\" causal mechanism as long as there exists a periodicity in how the causal mechanisms change over time. \n\nFor instance, if the causal effect changes by the end of the day but the samples are collected in every hour, considering each number as representing a distinct causal mechanism, the causal mechanism progression of $\\mathbf{X}^{j}$ over time could be [1,1,...,1,2,2,...,2,3,3,...,3,4,4,...,4,...,7,7,...,7,1,1,...,1,2,2,...,2,3,3,...,3,...4,4,...,4,...,7,7,...,7,...]. Even though the causal mechanism doesn't change every time point (every hour for example), the causal mechanism will eventually repeat itself every week and the algorithm will treat the sequence [1,1,...,1,2,2,...,2,3,3,...,3,4,4,...,4,...,7,7,...,7] for one week as a single unit, eventually leading to a periodicity of 24*7=168 hours=7 days. \n\nIn this situation, our proposed algorithm will have 168 correspondingly distinct time point subsets where the 1th, 169th, and 337th causal mechanism \"1\" will be included in the same time point subsets instead of having all causal mechanism \"1\" in the same time point subsets like [1,2,3,4,5,6,7,1,2,3,4,5,6,7,...,1,2,3,4,5,6,7]. Hence, the drawback is that the effective sample used in each MCI test for the first case will be 1/24 of the samples used for the second case. That is, the sample size used in MCI tests is shrunk. \n\nIn our ongoing research, our goal is to overcome this limitation and use the samples more efficiently in the first case.\n\n## Q2. \"sudden\" change\n\nAs for the \"sudden\" change in the causal mechanism, we can handle the \"soft\" change as long as the periodicity of the \"soft\" change mechanism is not a multiple of the periodicity of the \"hard\" change mechanism. \n\nFor instance, suppose the causal mechanism progression of $\\textbf{X}^{j}$ over time is [1,2,3,1,2,3,1,2,3,...] and the change between causal mechanism \"1\" and \"2\" is \"soft\". More specifically, suppose the only incoming edge of $X^j_t$ in causal mechanism \"1\" and \"2\" is $X^i_{t-1}$, but the causal effect between $X^i_{t-1}$ and $X^j_t$ in causal mechanism \"1\" is stronger than that in causal mechanism \"2\". This is a violation of the **Assumption A5** Hard Mechanism Change. \n\nSuppose in causal mechanism \"3\", **Assumption A5**  still holds. Our algorithm will treat this case as the same as [1,1,3,1,1,3,1,1,3,...]. Because of the causal mechanism \"3\", we will still have a periodicity of 3, and the time partition will still be correct. Even though there are no sudden changes between causal mechanisms \"1\" and \"2\", the corresponding samples will still be partitioned according to their causal mechanism, and then the causal effect will be estimated correctly. \n\nHowever, without causal mechanism \"3\", the algorithm will fail. \n\nOn the other hand, if the periodicity in the \"soft\" change mechanisms is 2 with [1,2,1,2,...], then the periodicity of the \"hard\" mechanism change is 1 with [1,1,1,1,...], therefore the algorithm will tell us that the periodicity is 1, instead of 2. Hence, the algorithm can handle the \"soft\" change as long as the periodicity of the \"soft\" change mechanism is not a multiple of the periodicity of the \"hard\" change mechanism. \n\nThank you for this great question that inspires us to relax the **Assumption A5** to handle \"soft\" change in some cases. We will update this in the camera-ready version.\n\n\n## Q3: The metrics.\n\nThank you for the comment, and we will include the explanation in the paper later. \n\nDuring our metric evaluation process, we will compute the least common multiple of $\\Omega$ and $\\hat{\\Omega}$, where $\\Omega$ itself represents the least common multiple of all $\\omega_{j}, j\\in[n]$, as defined in Eq. 7. Designating this least common multiple as \"lcm,\" we construct two binary arrays of edges with dimensions $[N, lcm, N, \\tau_{max}+1]$ for both the true causal graph and the estimated graph. In these arrays, the value 1 signifies an edge connecting one variable to another with a specific time lag, while 0 indicates the absence of an edge. The metrics are subsequently computed based on the difference between these two arrays of edges.\n\n## Q4: Standard error\n\nWe appreciate your suggestion, and we will make the necessary updates accordingly. To maintain consistency, we will continue to utilize the standard error for the new experiments. These updates will be synchronized with the previous experiments in the camera-ready version.\n\n## Q5: hyperparameters and CI tests.\n\nWe thank you for your comment, and we will incorporate these specific details into the main paper. \n\nFor the continuous-valued time series, $\\omega_{ub}=15$ for all cases no matter what $\\omega_{max}$ is, and $\\tau_{ub}=20$ for all cases. CI test used in continuous valued time series is a partial correlation test. For the discrete valued time series, $\\omega_{ub}=7$ and $\\tau_{ub}=7$ for all cases. CI test here is a conditional mutual information test.\n\n## Q6: Conclusion in the case study\n\nThe casual edge from $X^{ta}_{t-1} \\text{to} X^{cp}_t$ disappears in Jan, April, June, and October, while in the remaining months, the edge exists. Therefore we have an initial conclusion that the causal effect of tropical Atlantic air temperature on the Central Pacific air temperature disappears every quarter of a year.\n\n\n## Q7: \"discrete valued\" and \"continuous valued\" time series\n\nWe greatly appreciate your supportive comment, and we will revise the terminology accordingly to prevent any potential misunderstandings."}, {"review_id": "Oavf4GZBAl", "reviewer": "Reviewer_ztgZ", "summary": "The problem of causal discovery is tackled for time series data in those cases where stationarity of multivariate time series data can not be assumed. A problem formulation based on structural causal model is given by considering the case of semi-stationary time series. Under this assumptin, the paper considers the case where a finite number of different causal mechanisms happen in a sequential manner and periodically across time. A constraint based, and non parametric, structural learning algorithm is designed and developed by leveraging on the PCMCI, described in the specialized literature. Such an algorithm is extended and adapted to cope with the considered semi-stationary setting. The causal discovery problem is tackled by standard conditional independence tests. The paper proves that the proposd algorithm is sound when asked to discover and identify causal relationships between variables of discrete time series. Numerical experiments, both from synthetic and real world data, are reported and described to comment on the performance of the proposed approach. The main contributions of the paper are in my humble opinion the following:\n- designed and developing a new causal discovery algorithm for semi-stationary time series where the causal graph changes periodically.\n- validation of the proposed approach through syntethic data", "questions": "- which is the coputational complexity of the algorithm you proposed?\n- how does the proposed algrithm scales with respect to the number of variables of the multivariate time series?\n- could you please explain why the following works had not been taken into account? maybe they are slighlty out of the scope but still linked to the problem you tackled (Learning Continuous Time Bayesian Networks in Non-stationary Domains), (On-homogeneous  dynamic bayesian networks  with  bayesian  regularization  for  inferring  gene  regulatory  networks  with gradually time-varying structure) (Learning non-stationary dynamic bayesian networks)\n- line 104: could you please explain the meaning of setting P(V) different from 0?\n- line 105: the notation used for the pair of variables is not properly coherent, indeed Xj is a univariate timeseries contained in V, which has not been defined. Furthermore, what do you mean by variables here? two time serires? \n- line 115: what do you mean by finite maximum lag?\n- would you be that kind to comment explicitly and not only formally on the assumptions made about the time when nonstationarity occurs, i.e., the casual graph changes its structure? synchronous? asynchronous for each time series? could you discusse more on this aspect? what about parameters non stationarity? I guess you should consider this is excluded to happen, or not?\n- line 156: the definition of n is a little confusing. Indeed, n belongs to a set that is defined as a function of n again, would you be that kind to help me understand?\n- formula (14), the expernal p is the same probability ad the p inside brackets?\n- line 243, I read \"continuous time series in four steps\", but I found only three of them\n- would you also consider the structrual hamming distance to measure the performance of the proposed algorithm\n\n", "limitations": "- computational complexity of the proposed is not mentioned at all\n- comparison has not been performed with algorithm natively designed for semi-stationary time series data. Maybe no such algorithms exist?\n- results on real world data are of no help because no ground thruth is available and no domain expert comments are given to corroborate it", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- the tackled problem which is of both theoretical and practical interest\n- the theoretical framework of the considered problem\n", "weaknesses": "- the proposed algorithm wisely exploits the PCMCI algorithm while not being original, it looks a little as an incremental result\n- numerical experiments seem to compare the proposed algorithm to algortihms which are not conceived to tackle the considered problem\n- numerical experiments on real world data are disappointing. Indeed, the paper itself recognizes that not having knowledge about the considered domain the results obtained by the proposed algorithm can not be validated nor commented on.  Therefore, I ask myself which is the value of presenting such experiments?\n- computational complexity is not mentioned at all\n- the paper structure and organization makes it not trivial to follow it\n- I found many typos", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We really appreciate the feedback. Regarding your valuable questions, please see below for our responses.\n\n## Q1. Computational Complexity\nExecuting the PCMCI algorithm on the entire time series constitutes the initial phase of the proposed approach (designated as \"line 1\" in **Algorithm 1**). The algorithm's worst-case overall computational complexity is delineated by $O(n^{3}\\tau_{ub}^{2})$. Here, the symbol $n$ denotes the multivariate nature of the temporal data, and $\\tau_{ub}$ represents the upper boundary for time lags. \n\nThe subsequent computational load stemming from the remaining components of our algorithm follows a complexity of $O(\\omega_{ub}^2n^{2}\\tau_{ub})$. This encompasses the $O(n^{2}\\tau_{ub})$ complexity associated with conducting Momentary Conditional Independence (MCI) tests on all $n$ univariate time series. The parameter $\\omega_{ub}^2$ here arises due to the search procedure involving $\\omega_{j}$, iterating through values from $1$ to $\\omega_{ub}$ for all $n$ univariate time series.\n\nThe runtime of the computation is further influenced by the scaling behavior of the conditional independence test concerning the dimensionality of the conditioning set and the temporal series length $T$. For further details, see section 5.1 in the work by Runge et al.[2019]. \n## Q2. Novelty\nThank you for your comment.  We agree that our solution leverages PCMCI and builds on it. However, we believe that the novelty of our approach lies in simplicity by carefully leveraging it using algorithmic ideas. Please note that prior to our work, it was not clear if we needed to rethink the causal discovery problem from scratch for non-stationary settings or if some of the existing work could be leveraged in a clever way. Our algorithmic contribution demonstrates the latter. \n## Q3.  Lack of appropriate Baselines\nWe agree with the reviewer. To the best of our knowledge, there are no non-parametric baselines that can handle non-stationarity or semi-stationarity in the data with a window causal graph. \n## Q4. Experiment in the case study\n\nWhile we understand that our results are based on simulated data, it is fairly common in the area to benchmark against the simulated settings that we considered. \nAs one of the main contributions of our work is the introduction of the problem, we believe that our empirical contributions along with the problem formulation are significant.  Moreover, many results in causal discovery and inference rely on synthetic simulations, which is necessary in the absence of ground truth.\n\n## Q5. Related Literature\nWe appreciate the valuable literature you've shared; it bears significant relevance to our ongoing project. Regarding our proposed algorithm here, we emphasize a non-parametric approach capable of uncovering causal relationships with time-lag effects and adapting to both abrupt changes in causal mechanisms at each time point while also addressing consistent causal mechanisms as long as periodicity is present. In contrast, the works mentioned are not distribution-free and generate summary causal graphs as outputs. \n\n## Q6. $P(V)\\neq 0$\nThe probability of any given realization involving all variables within $V$ should not be zero. For instance, consider $V=\\\\{X,Y\\\\}$; when computing $P(X|Y)=\\frac{P(X,Y)}{P(Y)}$, we aim to avoid scenarios where the denominator or numerator becomes zero. This is needed in the proof, and this prerequisite is already met by **Assumption A7**. $P(V) \\neq 0$ is also called the Positivity assumption.\n\n## Q7. Clarification on some notations in Line 105\nThe notation $X, Y\\in V$ signifies that any two random variables, denoted by $X$ and $Y$, belong to the set $V$.  $X, Y\\in \\mathbb{R}^{1}$. Meanwhile, $S$ represents a subset of $V$. As indicated in line 101, $X_{t}^{j} \\in \\mathbb{R}^{1}$ designates an individual variable, whereas $\\mathbf{X}^{j}_{t\\in[T]} \\in \\mathbb{R}^{T}$ represents a time series.\n## Q8. Finite maximum lag\nThis implies that as $t$ approaches infinity, $X^i_1$ will not serve as a cause for $X^j_t$. The maximum time lag $\\tau_{\\text{max}}$, as defined in Eq.1, must remain finite.\n## Q9. Non-stationarity in the Semi-Stationary setting\nEach univariate time series $\\mathbf{X}^j$ in the n-variate time series $V$ can have its own periodicity $\\omega_j, j\\in[n]$ in the causal mechanisms. Within each cycle of $\\omega_j$ time points, the same causal mechanism is reiterated, resulting in $\\omega_j$ distinct causal mechanisms for the time series $\\mathbf{X}^j$.\n\nIn accordance with **Assumption A6**, when the causal mechanisms governing two variables, $X^j_{t1}$ and $X^j_{t2}$, have undergone alteration, it follows that their respective parent sets cannot remain identical. In essence, when the parameters within the structural causal model shift, there is a corresponding change in the parent set as well.\n## Q10. Line 156, definition of $n$\n$n$ belongs to a subset of $\\mathbb{N}$ that should satisfying $\\tau_{max}+q-1+n\\delta \\leq T$, i.e, $n\\leq (T-(\\tau_{max}+q-1))/\\delta$\n## Q11. Formula 14, two $p$\nYes, they both denote the probability. That means the chance of the two probabilities enclosed within the brackets being identical is zero.\n## Q12. Line 243 typo\nThere should be only three steps.\n## Q13. Structrual hamming distance\nWe did not use structural hamming distance because of the semi-stationary context, wherein the estimated binary edge array's dimensionality grows with rising estimated $\\hat\\omega_{j}$. As a result, larger $\\hat\\omega_{j}$ values inherently amplify the dissimilarity in SHD between the true and estimated edge arrays. This complexity makes direct comparisons challenging. Precision, recall, and F1 score, on the other hand, are intrinsically normalized within the range of [0,1], rendering them suitable metrics for the semi-stationary context."}]

Author rebuttal:
Dear reviewers,

We extend our sincere gratitude for dedicating your time to meticulously assess our submission and offer invaluable insights.

With the initial rebuttal phase comes to an end, we believe that our response effectively addresses the raised points of concern.

We have attached a one-page pdf with new experiment results.

We remain fully open to continued discussions throughout the forthcoming discussion phase.

In light of our responses and the implemented revisions, we kindly request you to reconsider your evaluation.

Warm regards, Authors

Meta-review: This paper is concerned with causal discovery method for semi-stationary time series, where causal structure may periodically change over time. A constraint-based, non-parametric, structural learning algorithm is proposed based on the PCMCI approach. It is extended and adapted to cope with the considered semi-stationary setting. The causal discovery problem is tackled by standard conditional independence tests. The paper shows that the proposed algorithm is sound in the discovered causal relationships. Experiments on both synthetic and real-world data are reported to illustrate the good performance of the proposed approach. Overall, the considered setting is practical and theoretical and empirical results in the paper demonstrate that the proposed method is a nice contribution to the community.
Decision: Accept (poster)