# Causal Discovery in Semi-Stationary Time Series

Shanyun Gao
Purdue University
<EMAIL_ADDRESS>
&Raghavendra Addanki
Adobe Research
<EMAIL_ADDRESS>
&Tong Yu
Adobe Research
<EMAIL_ADDRESS>
&Ryan A. Rossi
Adobe Research
<EMAIL_ADDRESS>
&Murat Kocaoglu
Purdue University
<EMAIL_ADDRESS>

###### 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Ω,
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 to a real-world climate dataset.

##  1 Introduction

In modern sciences, causal discovery aims to identify the collection of causal
relations from observational data, as in Pearl (1980), Peters et 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 is causal Markov and faithful
to a 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 of the variables. See Peters et al.
(2017). There are many constraint-based algorithms such as PC and FCI
algorithms Spirtes et 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-based
algorithms to accommodate time series data. For instance, PCMCI in Runge et
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 in 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 as root-cause analysis in Ikram et al. (2022).
Most of the existing causal discovery algorithms make the stationary
assumption. See Chu and Glymour (2008), Hyvärinen et al. (2010), Entner and
Hoyer (2010), Peters et al. (2013), Malinsky and Spirtes (2018), Runge et al.
(2019), Pamfil et 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 to 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 a 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
Fig.1, where the time series $\mathbf{X}^{1}$ has three different causal
mechanisms across time, shown as red edges, green edges, and blue edges,
respectively. Similarly, time series $\mathbf{X}^{2}$ 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 et
al. (2002), Nakamura et al. (2003), Carskadon et al. (2005) 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 (2018) 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. (2019). These algorithms are
discussed further in the related work.

In this paper, we propose an algorithm to address this problem, namely non-
parametric causal discovery in time series data with semi-stationary SCMs. The
key contributions of our work are:

  * •

We develop an algorithm to discover the causal structure from semi-stationary
time series data where the underlying causal structures change periodically.
Our algorithm systematically uses the PCMCI algorithm proposed for the
stationary setting in Runge et al. (2019). The resulting algorithm is hence
named PCMCIΩ where $\Omega$ 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 a 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), water usage as studied in Zou et
al. (2022), ultra-processed food manufacturing as examined in Menegozzo et 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. (2022, 2023). While PCMCI has achieved considerable success, it
is not without its limitations. One notable assumption that can be challenged
is the concept of causal stationary, 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. (2019). However, there is an
ongoing exploration to enhance its performance in a 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-
parametric algorithm in 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. (2023)
proposed an algorithm JIT-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-parametric
approaches can discover causal relations underlying time series without
assuming stationarity and can also allow for sudden changes in causal
mechanisms.

Our proposed approach does not directly enforce the stationary assumption on
the time series. The SCM also integrates a finite set of causal mechanisms
that exhibit periodic variations over time.

Figure 1: Partial causal graph for 3-variate time series
$V=\\{\mathbf{X}^{1},\mathbf{X}^{2},\mathbf{X}^{3}\\}$ with a Semi-Stationary
SCM where $\tau_{\max}=3$, $\omega_{1}=3$, $\omega_{2}=2$, $\omega_{3}=1$,
$\Omega=6$ and $\delta=6$. The first 3(=$\tau_{\max}$) time slices
$\\{\mathbf{X}_{t}\\}_{1\leq t\leq 3}$ are the starting points. The same color
edges represent the same causal mechanism. E.g. for $\mathbf{X}^{1}$: there
are 3 ($=\omega_{1}$) time partition subsets $\\{\Pi^{1}_{k}\\}_{1\leq k\leq
3}$. The time points $t$ of nodes $X^{1}_{t}$ sharing the same filling color
are in the same time partition subsets. The time points $t$ of nodes
$X^{1}_{t}$ sharing both the same filling color and the same outline shape are
in the same homogenous time partition subsets (the definitions are in the
supplementary material). There are 6 ($=\delta$) different Markov chains in
this multivariate time series $V$, and the first element of these 6 Markov
chains is shown as $\\{Z^{q}_{1}\\}_{1\leq q\leq 6}$ and are tinted with a
gradient of blue hues. The superscript $q$ of $Z^{q}_{i}$ is the index of
different Markov chains, whereas the subscript $i$ denotes the running index
of that specific Markov chain. For instance, $Z^{1}_{1}$ and $Z^{1}_{2}$
denote the first two elements of the first Markov chain, while $Z^{2}_{1}$ and
$Z^{2}_{2}$ denote the first two elements of the second Markov chain.

##  2 PCMCIΩ: Capturing Periodicity of the Causal Structure

In this section, we formulate the problem of learning the causal graph on
multi-variate time series data when the SCM exhibits periodicity in causal
mechanisms. In section 2.1, 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 $\mathcal{G}(V,E)$ denote the underlying causal graph, and for each
variable $X\in V$, we denote the set of all incoming neighbors as parents,
denoted by $\mathrm{Pa}(X)$.

For any two variables $X,Y\in V$ and $S\subset V$, we denote the CI relation
$X$ is independent of $Y$ conditioned on $S$, by $X\perp\\!\\!\\!\perp Y\mid
S$.

For simplicity’s sake, define sets: $[b]:=\\{1,2,...,b\\}$ and
$[a,b]:=\\{a,a+1,...,b\\}$, where $a,b\in\mathbb{N}$.

In the time series setting, let $X^{j}_{t}\in\mathbb{R}^{1}$ denote the
variable of $j$th time series at time $t$,
$\mathbf{X}^{j}=\\{X^{j}_{t}\\}_{t\in[T]}\in\mathbb{R}^{T}$ denote a
univariate time series and
$\mathbf{X}_{t}=\\{X^{j}_{t}\\}_{j\in[n]}\in\mathbb{R}^{n}$ denote a slice of
all variables at time point $t$.
$V=\\{\mathbf{X}^{j}\\}_{j\in[n]}=\\{\mathbf{X}_{t}\\}_{t\in[T]}\in\mathbb{R}^{n\times
T}$ denotes a $n$-variate time series. By default, we assume $n>1$ and hence
$\mathbf{X}^{j}\subsetneq V$, and $p(V)\neq 0$, where $p(.)$ denotes the
probability or probability density.

######  Definition 2.1 (Non-Stationary SCM).

A Non-Stationary Structural Causal Model (SCM) is a tuple $\mathcal{M}=\langle
V,\mathcal{F},\mathcal{E},\mathbb{P}\rangle$ where there exists a
$\tau_{\max}\in\mathbb{N}^{+}$, defined as:

$\displaystyle\tau_{\max}$
$\displaystyle\coloneqq\arg\max_{\tau}\\{\tau:X^{i}_{t-\tau}\in\mathrm{Pa}(X^{j}_{t}),i,j\in[n]\\},$
(1)

such that with this $\tau_{\max}$, each variable $X^{j}_{t>\tau_{\max}}\in V$
is a deterministic function of its parent set
$\mathrm{Pa}(X^{j}_{t>\tau_{\max}})\in V$ and an unobserved (exogenous)
variable $\epsilon^{j}_{t>\tau_{\max}}\in\mathcal{E}$:

$\displaystyle X^{j}_{t}$
$\displaystyle=f_{j,t}(\mathrm{Pa}(X^{j}_{t}),\epsilon^{j}_{t}),\;j\in[n],t\in[\tau_{\max}+1,T],$
(2)

and there exist at least two different time points
$t_{0},t_{1}\in[\tau_{\max}+1,T]$ satisfying

$\displaystyle f_{j,t_{0}}$ $\displaystyle\neq f_{j,t_{1}},\ \
\exists\;j\in[n],\\{t_{0},t_{1}\\}\subset[\tau_{\max}+1,T].$ (3)

where $f_{j,t},f_{j,t_{0}},f_{j,t_{1}}\in\mathcal{F}$ and
$\\{\epsilon^{j}_{t}\\}_{t\in[T]}$ are jointly independent with probability
measure $\mathbb{P}$. $\tau_{\max}$ is the finite maximal lag in terms of the
causal graph $\mathcal{G}$.

######  Definition 2.2 (Semi-Stationary SCM).

A Semi-Stationary SCM is a Non-Stationary SCM that additionally satisfies the
following conditions. $\text{ For each }j\in[n],\text{ there exists an
}\omega\in\mathbb{N}^{+}$such that:

$\displaystyle a)\;f_{j,t}=f_{j,t+N\omega},$ (4) $\displaystyle
b)\;\mathrm{Pa}(X^{j}_{t+N\omega})=\\{X^{i}_{s+N\omega}:X^{i}_{s}\in\mathrm{Pa}(X^{j}_{t}),i\in[n]\\},$
(5) $\displaystyle c)\;\epsilon^{j}_{t},\epsilon^{j}_{t+N\omega}\text{ are
i.i.d.}$ (6)

are satisfied for all
$t\in[\tau_{\max}+1,T],N\in\\{N:N\in\mathbb{N},t+N\omega\leq T\\}$. This means
that a finite number of causal mechanisms are repeated periodically for every
univariate time series $\mathbf{X}^{j}$ in $V$. One example of this model is
illustrated in Fig.1. For $\mathbf{X}^{1}$ in Fig.1, three causal mechanisms
are reiterated periodically with $\omega_{1}=3$, represented by red, green,
and blue edges, respectively.

The minimum value that satisfies the above conditions for $\mathbf{X}^{j}$ is
defined as the _periodicity_ of $\mathbf{X}^{j}$, denoted by $\omega_{j}$.
Furthermore, for an $n$-variate time series $V$, $\Omega$ denotes the minimum
periodicity across all time series $\mathbf{X}^{j},j\in[n]$. The number of
causal mechanisms occurring sequentially and periodically of univariate time
series $\mathbf{X}^{j}$ is $\omega_{j}$, and that number of causal mechanisms
of multivariate time series $V$ is $\Omega$. For $\mathbf{X}^{j}$, the causal
mechanisms are associated with each variable $X^{j}_{t}$. However, for $V$,
the causal mechanisms are related to each time slice vector $\mathbf{X}_{t}$
as a whole. The relationship between $\Omega$ and $\omega_{j}$ can be captured
by:

$\displaystyle\Omega=\text{LCM}(\\{\omega_{j}:j\in[n]\\})$ (7)

where $\text{LCM}(.)$ is an operation to find the least common multiple
between any two or more numbers. Here, $\Omega$ is the smallest common
multiple among $\\{\omega_{1},...\omega_{n}\\}$. In Fig.1,
$\Omega=\text{LCM}(3,2,1)=6$.

######  Definition 2.3 (Time Partition).

A time partition $\Pi^{j}(T)$ of a univariate time series $\mathbf{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_{k}^{j}(T)\\}_{k\in[\omega_{j}]}$ such that:

$\displaystyle\Pi^{j}_{k}(T):=\\{t:\tau_{\max}+1\leq t\leq
T,(t\bmod\omega_{j})+1=k\\}.$ (8)

where $\bmod$ denotes the modulo operation. For instance, $5\bmod 3=2$.

We can observe that the variables in $\\{X_{t}^{j}\\}_{t\in{\Pi^{j}_{k}(T)}}$
share the same causal mechanism. Since the number of potentially different
causal mechanisms of variables in $\mathbf{X}^{j}$ is $\omega_{j}$, the number
of such time partition subsets is $\omega_{j}$. For simplicity, notations
$\Pi^{j}$ and $\Pi^{j}_{k}$ are used instead of $\Pi^{j}(T)$ and
$\Pi^{j}_{k}(T)$. In Fig.1, $\Pi^{1}_{1}=\\{4,7,10,13,..,4+3N,...\\}$,
$\Pi^{1}_{2}=\\{5,8,11,14,..,5+3N,...\\}$ and
$\Pi^{1}_{3}=\\{6,9,12,15,..,6+3N,...\\}$ where $N\in\mathbb{N}^{+}$. The
nodes $X^{1}_{t}$ are classified into their associated time partition subsets
by the matching colors.

######  Definition 2.4 (Illusory Parent Sets).

For a univariate time series $\mathbf{X}^{j}\in V$ with Semi-Stationary SCM
having periodicity $\omega_{j}>1$, parent set index
$\mathrm{pInd}^{j}_{k\in[\omega_{j}]}$ is defined as:

$\displaystyle\mathrm{pInd}^{j}_{k}\coloneqq\\{(y_{i},\tau_{i})\\}_{i\in[n^{\prime}]}\text{,
given
}\mathrm{Pa}(X^{j}_{t})=\\{X^{y_{1}}_{t-\tau_{1}},X^{y_{2}}_{t-\tau_{2}},...,X^{y_{n^{\prime}}}_{t-\tau_{n^{\prime}}}\\},\;\forall
t\in\Pi^{j}_{k}$ (9)

where $n^{\prime}=|\mathrm{Pa}(X^{j}_{t}))|$, $\tau_{i}$ is the time lag and
$y_{i}$ is the variable index. Given $\mathrm{pInd}^{j}_{k}$, Illusory Parent
Sets are defined as:

$\displaystyle\mathrm{Pa}_{k}(X^{j}_{t})=\bigl{\\{}X^{y_{i}}_{t-\tau_{i}}:(y_{i},\tau_{i})\in\mathrm{pInd}^{j}_{k}\bigr{\\}},\;\forall
k\in\\{k:t\notin\Pi^{j}_{k}\\}$ (10)

Put simply, the illusory parent sets of $X^{j}_{t}$ are the time-shift version
of the parent set of other variables in $\mathbf{X}^{j}$ that have a different
causal mechanism from $X^{j}_{t}$. Note that the illusory parent sets are
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 $t\in\Pi_{k}^{j}$, the notation $\mathrm{Pa}_{k}(X_{t}^{j})$ can
also be extended to encompass the true parent set of $X^{j}_{t}$:

$\displaystyle\mathrm{Pa}_{k}(X_{t}^{j})$
$\displaystyle:=\mathrm{Pa}(X_{t}^{j}),\;\forall t\in\Pi_{k}^{j}$ (11)

By doing so,
$\mathrm{Pa}(X_{t}^{j})\subset\cup_{k\in[\omega_{j}]}\mathrm{Pa}_{k}(X_{t}^{j})$.
In Fig.1, by observing $\mathrm{Pa}(X_{7}^{1})$, we have
$\mathrm{pInd}^{1}_{1}=\\{(1,1),(2,2)\\}$; by observing
$\mathrm{Pa}(X_{8}^{1})$, we have $\mathrm{pInd}^{1}_{2}=\\{(1,1),(3,1)\\}$
and finally by observing $\mathrm{Pa}(X_{9}^{1})$,
$\mathrm{pInd}^{1}_{3}=\\{(1,1),(1,2)\\}$. Based on those indexes,
$Pa_{1}(X_{7}^{1})=\\{X^{1}_{7-1},X^{2}_{7-2}\\}=\\{X^{1}_{6},X^{2}_{5}\\}$,
$Pa_{2}(X_{7}^{1})=\\{X^{1}_{7-1},X^{3}_{7-1}\\}=\\{X^{1}_{6},X^{3}_{6}\\}$
and
$Pa_{3}(X_{7}^{1})=\\{X^{1}_{7-1},X^{1}_{7-2}\\}=\\{X^{1}_{6},X^{1}_{5}\\}$.
The first one is the true parent set of $X_{7}^{1}$ and the latter two are the
illusory parent 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 are (potentially) $\delta$
different Markov chains $\\{Z^{q}_{n}\\}_{n\in\mathcal{N}},q\in[\delta]$:

$\displaystyle Z^{q}_{n}$
$\displaystyle=\\{\mathbf{X}_{\tau_{\max}+q+(n-1)\delta},\mathbf{X}_{\tau_{\max}+q+1+(n-1)\delta},...,\mathbf{X}_{\tau_{\max}+q-1+n\delta}\\},$

where $\mathcal{N}\coloneqq\\{n\in\mathbb{N}^{+}:\tau_{\max}+q-1+n\delta\leq
T\\}$, $\delta=\lceil\frac{\tau_{\max}+1}{\Omega}\rceil\Omega$. Note that in
$\\{Z^{q}_{n}\\}$, $q$ is used to indicate a specific Markov chain, while $n$
serves as the running index for that particular Markov chain. Such a Markov
chain $\\{Z^{q}_{n}\\},q\in[\delta]$ exists as long as
$\mathrm{Pa}(Z^{q}_{n})\subset Z^{q}_{n}\cup Z^{q}_{n-1}$ for all $n$. This is
a finite state Markov Chain if all time series in $V$ are discrete-valued time
series. The state space of $\\{Z^{q}_{n}\\}$ is the set containing all
possible realization of
$\\{\mathbf{X}_{\tau_{\max}+q+(i-1)+(n-1)\delta}\\}_{i\in[\delta],n\in\mathbb{N}}$.
The transition probabilities between the states are the product of associated
causal mechanisms based on Assumption A2, which is elaborated by an example in
section C.1 (Eq.(10)-(11)) of the supplementary material.

###  2.2 Assumptions for PCMCIΩ

  1. A1.

Sufficiency: There are no unobserved confounders.

  2. A2.

Causal Markov Condition: Each variable $X$ is independent of all its non-
descendants, given its parents $\mathrm{Pa}(X)$ in $\mathcal{G}$.

  3. A3.

Faithfulness Condition (Pearl (1980)): Let $P$ be a probability distribution
generated by $\mathcal{G}$. $\langle\mathcal{G},P\rangle$ satisfies the
Faithfulness Condition if and only if every conditional independence relation
true in $P$ is entailed by the Causal Markov Condition applied to
$\mathcal{G}$.

  4. A4.

No Contemporaneous Causal Effects: Edges between variables at the same time
are not allowed.

  5. A5.

Temporal Priority: Causal relations that point from the future to the past are
not permitted.

  6. A6.

Hard Mechanism Change:

If at time points $t_{1}$ and $t_{2}$, the causal mechanisms of
$X^{j}_{t_{1}}$ and $X^{j}_{t_{2}}$ are different, then their corresponding
parent sets can not be transformed to each other by time shifts:

$f_{j,t_{1}}\neq
f_{j,t_{2}}\Rightarrow\mathrm{Pa}(X^{j}_{t_{2}})\neq\\{X^{i}_{s+(t_{2}-t_{1})}:X^{i}_{s}\in\mathrm{Pa}(X^{j}_{t_{1}}),i\in[n]\\}.$

  1. A7.

Irreducible and Aperiodic Markov Chain: The Markov chains $\\{Z_{n}\\}$ of $V$
are assumed to be irreducible (Serfozo (2009)): for all states $i$ and $j$ of
$\\{Z_{n}\\}$, $\exists n$ so that

$\displaystyle p_{ij}^{(n)}\coloneqq p(Z_{n+1}=j|Z_{1}=i)>0$ (12)

and aperiodic(Karlin (2014)): for every state $i$ of $\\{Z_{n}\\}$, $d(i)=1$,
where the period $d(i)$ of the state $i$ is the greatest common divisor of all
integers $n$ for which $p_{ii}^{(n)}>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 in $\\{X^{j}_{n}\\}_{n\in[T]}$ if they share
identical parent sets after time shifts. As for A7, it serves a crucial role
in establishing the soundness of our algorithm.

##  3 PCMCIΩ Algorithm

Algorithm 1 PCMCIΩ

1:Input: A $n$-variate time series
$V=(\mathbf{X}^{1},\mathbf{X}^{2},\mathbf{X}^{3},...,\mathbf{X}^{n})$,
periodicity upper bound $\omega_{\text{ub}}$, time lag upper bound
$\tau_{\text{ub}}$. By default, we assume $\tau_{\text{ub}}$ and
$\omega_{\text{ub}}$ are larger than their true value.

2:A superset of parent set is obtained using PCMCI with $\tau_{\text{ub}}$ and
denote it by $\widehat{S\mathrm{Pa}}(X^{j}_{t})\ \forall j,t$.

3:for $\mathbf{X}^{j}$ where $j\in[n]$ do

4: for a guess $\omega\in[\omega_{ub}]$ of $\omega_{j}$ do

5: Let $\widehat{\Pi}^{j}:=\\{\widehat{\Pi}^{j}_{k}|k\in[\omega]\\}$ where
$\widehat{\Pi}^{j}_{k}=\\{2\tau_{\text{ub}}+k,2\tau_{\text{ub}}+\omega+k,2\tau_{\text{ub}}+2\omega+k,\cdots\\}$.

6: for $k\in[\omega]$ do

7: Initialize the parent set for $X^{j}_{t},t\in\\{t:t\geq
2\tau_{\text{ub}},t\in\widehat{\Pi}^{j}_{k}\\}$ (with guess $\omega$) denoted
by
$\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})\leftarrow\widehat{S\mathrm{Pa}}(X^{j}_{t})$.

8: Consider $X^{i}_{t-\tau}\in\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$.
Remove $X^{i}_{t-\tau}$ from $\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$ if
$X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\mid\left(\widehat{S\mathrm{Pa}}(X^{j}_{t})\cup\widehat{S\mathrm{Pa}}(X^{i}_{t-\tau})\right)\setminus
X^{i}_{t-\tau}$ using a CI Test with samples $t\in\\{t:t\geq
2\tau_{\text{ub}},t\in\widehat{\Pi}^{j}_{k}\\}$.

9: Store $\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$ for
$X^{j}_{t},t\in\\{t:t\geq 2\tau_{\text{ub}},t\in\widehat{\Pi}^{j}_{k}\\}$.

10: end for

11: end for

12:
$\widehat{\omega}_{j}\leftarrow\arg\min_{\omega\in[\omega_{\text{ub}}]}\max_{k\in[\omega]}|\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t\in\widehat{\Pi}^{j}_{k}})|$.

13: Set
$\widehat{\mathrm{Pa}}(X^{j}_{t})\leftarrow\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})$
for $X^{j}_{t},t\in\\{t:t\geq 2\tau_{\text{ub}}\\}$.

14:end for

15:return $\hat{\omega}_{j}$ and $\widehat{\mathrm{Pa}}(X^{j}_{t})\ \forall
j\in[n],t\geq 2\tau_{\text{ub}}$.

In this section, we propose an algorithm called $\text{PCMCI}_{\Omega}$, and
in section 3.1, we present a theorem demonstrating the soundness of
$\text{PCMCI}_{\Omega}$ and its ability to recover the causal graph. Our
algorithm $\text{PCMCI}_{\Omega}$ builds on the Algorithm PCMCI in Runge et
al. (2019). Additional details about PCMCI are provided in the supplementary
material.

Overview of Algorithm 1 $\text{PCMCI}_{\Omega}$. We assume that the
periodicity and time lag are upper bounded by $\omega_{\text{ub}}$ and
$\tau_{\text{ub}}$ respectively. Using PCMCI Runge et al. (2019), we obtain a
superset of parents for every variable $X^{j}_{t}$ denoted by
$\widehat{S\mathrm{Pa}}(X^{j}_{t})$ (line 2). Our goal is to identify the
correct set of parents along with its periodicity for every variable in $V$.
For a variable $X^{j}_{t}$, we guess its periodicity $\omega$ by iterating
over all possible values in $[\omega_{\text{ub}}]$. Next, we construct time
partition subsets $\widehat{\Pi}^{j}_{k},\;k\in[\omega]$ based on the guess of
periodicity $\omega$. In each time partition subset, we maintain a parent set,
denoted by $\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$, initializing it with
the superset $\widehat{S\mathrm{Pa}}(X^{j}_{t})$. Then we test the causal
relations between $X^{i}_{t-\tau}\in\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$
and $X^{j}_{t}$ using a CI test on the sample $t\in\widehat{\Pi}^{j}_{k}$
(lines 6-10).

For each guess $\omega$, every variable in $\mathbf{X}^{j}$ should have its
estimated parent set (line 9), and there are total $\omega$ potentially
different parent set index $\mathrm{pInd}^{j}_{k\in[\omega]}$ in
$\mathbf{X}^{j}$. We return an estimate $\hat{\omega}_{j}$ that maximizes the
sparsity of the causal graph (Lemma 3.4). Therefore, we select the value of
$\omega\in[\omega_{\text{ub}}]$ that minimizes the maximum value of
$|\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})|,\;t\in[T]$ as the estimator of
$\omega_{j}$ (line 12).

###  3.1 Theoretical Guarantees

Our main theorem shows that $\text{PCMCI}_{\Omega}$ 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 $\hat{\mathcal{G}}$ be the estimated graph using the Algorithm
$\text{PCMCI}_{\Omega}$. Under assumptions A1-A7 and with an oracle (infinite
sample size limit), we have that:

$\displaystyle\hat{\mathcal{G}}=\mathcal{G}$ (13)

almost surely.

Lemma 3.2 and Lemma 3.3 jointly state that if CI tests are conducted on
samples generated by different causal mechanisms, the obtained parent sets
$\widehat{S\mathrm{Pa}}(X^{j}_{t})$ should be the superset of the union of the
true and illusory parent sets $\cup_{k}\mathrm{Pa}_{k}(X^{j}_{t})$. That is,
the estimated graph should be denser than the correct graph. The true parent
set can then be obtained by directly testing the independent relations between
the target variable $X^{j}_{t}$ and the variables in
$\widehat{S\mathrm{Pa}}(X^{j}_{t})$, assuming a consistent CI test. Note that
the CI tests in our algorithm are assumed to be consistent given i.i.d.
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.i.d.
samples can be seamlessly integrated into our algorithm.

######  Lemma 3.2.

Denote that $\\{\mathrm{Pa}_{k}(X^{j}_{t})\\}_{k\in[\omega_{j}]}$ contain the
true and illusory parent sets, where $\omega_{j}$ is the true periodicity of
$\mathbf{X}^{j}$. For any random variable $X_{t}^{j}$ with large enough $t$,
under assumptions A1-A7 and with an oracle (infinite sample size limit), we
have:

$\displaystyle
p\bigg{(}p(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))\neq
p(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\setminus
y)\bigg{)}=1,\;\forall y\in\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})$
(14)

Here,
$p(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))=\lim_{T\rightarrow\infty}\hat{p}(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))$
where $\hat{p}(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))$
is an estimated conditional distribution using all samples
$t\in[\tau_{\max}+1,T]$:

$\displaystyle\hat{p}(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))$
$\displaystyle=\frac{\sum_{t}\mathbbm{1}(X^{j}_{t},\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))}{\sum_{t}\mathbbm{1}(\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))}.$
(15)

###### Proof sketch.

We argue that the estimated conditional distribution in Eq.(15) can be written
as a linear combination of $\hat{p}(X^{j}_{t}|\mathrm{Pa}(X^{j}_{t}))$ where
$t\in\Pi^{j}_{k},k\in[\omega_{j}]$, i.e., as a mixture of conditional
distributions. The coefficients in the linear function, say
$\alpha_{k},k\in[\omega_{j}]$, can be further decomposed based on a finer time
partition called the _homogenous time partition_ , which consists of subsets
constructed according to the Markov chains $\\{Z^{q}_{n}\\}_{q\in[\delta]}$
corresponding to the time series. Based on Assumption A7, the Markov chains
are stationary and ergodic. Therefore, after sufficiently large time steps,
the distribution of $\\{Z^{q}_{n}\\}_{q\in[\delta]}$ will be invariant across
$n$ as it achieves unique equilibrium. With this type of stationary sample, we
can express $\alpha_{k}$ 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.

∎

######  Lemma 3.3.

Let $\widehat{S\mathrm{Pa}}(\mathbf{X}^{j}_{t})$ denote the estimated superset
of parent set for $\mathbf{X}^{j}\in V$ obtained from the Algorithm 1 (line
2). $\\{\mathrm{Pa}_{k}(X^{j}_{t})\\}_{k\in[\omega_{j}]}$ contain the true and
illusory parent sets, where $\omega_{j}$ is the true periodicity of
$\mathbf{X}^{j}$. Under assumptions A1-A7 and with an oracle (infinite sample
size limit), we have that:

$\displaystyle\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\subseteq\widehat{S\mathrm{Pa}}(X^{j}_{t}),\;\forall
t\in[\tau_{\max}+1,T]$

almost surely.

###### Proof sketch.

Assume the contrary, i.e., there exists
$s\in\cup_{k}\mathrm{Pa}_{k}(X_{t}^{j})\setminus\widehat{S\mathrm{Pa}}(X^{j}_{t})$.
From Lemma 3.2, we have
$X^{j}_{t}\cancel{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}s\left|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\setminus
s\right.$. By the Definition 2.4, we have
$\mathrm{Pa}(X^{j}_{t})\subset\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})$.
If $s\not\in\mathrm{Pa}(X^{j}_{t})$, by the causal Markov property (Assumption
A2), the dependence relation can not be true because $s$ is a non-descendant
of $X^{j}_{t}$. If $s\in\mathrm{Pa}(X^{j}_{t})$, our Algorithm would have
concluded that
$X_{t}^{j}\cancel{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}s\left|\widehat{S\mathrm{Pa}}(X^{j}_{t})\right.$
(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 $\omega_{j}$ for
$\mathbf{X}^{j}$ through Lemma 3.4.

######  Lemma 3.4.

Let $\omega_{j}$ denote the true periodicity for $\mathbf{X}^{j}\in V$ and
$\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t\in\Pi^{j}_{k}})$ denote the estimated
parent set for $X^{j}_{t}$ obtained from Algorithm 1 line 9, where
$t\in\Pi^{j}_{k}$. Define:

$\displaystyle\widehat{\omega}_{j}=\arg\min_{\omega\in[\omega_{\text{ub}}]}\max_{k\in[\omega]}|\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t\in\Pi^{j}_{k}})|$
(16)

Under assumptions A1-A7 and with an oracle (infinite sample limit), we have
that $\hat{\omega}_{j}=\omega_{j},\;\forall j\in[n]$ almost surely.

###### Proof sketch.

Assume the contrary that $\hat{\omega}_{j}\neq\omega_{j}$, then in the
Algorithm 1, we have an incorrect time partition $\widehat{\Pi}^{j}$. Hence,
CI tests that are performed use samples with different causal mechanisms.
$\hat{p}(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))$ in
Eq.(15) is estimated from a mixture of two or more time partition subsets, say
$\Pi^{j}_{1}$ and $\Pi^{j}_{2}$. We can apply Lemma 3.2 with
$\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})$. With Lemma 3.3,
$\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})\subseteq\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})$.
Hence, for $\hat{\omega}_{j}$,
$|\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})|\geq|\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})|$
using mixture samples $t\in\cup_{k=1}^{2}\Pi^{j}_{k}$. For true $\omega_{j}$,
we have
$|\widehat{\mathrm{Pa}}_{\omega_{j}}(X^{j}_{t})|=|\mathrm{Pa}(X^{j}_{t})|$.
With Assumption A6 the Hard Mechanism Change,
$|\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})|>|\mathrm{Pa}(X^{j}_{t})|$ so that
$\omega_{j}$ always leads to a smaller size of estimated parent sets than
$\hat{\omega}_{j}$, contrary to the definition of $\hat{\omega}_{j}$. Hence,
$\hat{\omega}_{j}=\omega_{j}$.

∎

##  4 Experiments

###  4.1 Experiments on Continuous-valued Time Series

(a) Accuracy of $\hat{\omega}$

(b) Runtime (in sec.) when $\omega_{\max}$ varies

Figure 2:  $\text{PCMCI}_{\Omega}$ is tested on 5-variate time series with
$\tau_{\max}=5$. Set $\tau_{\text{ub}}=15,\omega_{\text{ub}}=15$ for all
variables. Every line corresponds to a different time series length. Every
marker corresponds to the average accuracy rate or average running time over
100 trials. a) The accuracy rate of $\hat{\omega}$ for different time series
lengths and different $\omega_{\max}$. b) Illustration of Runtime (in sec.)
when $\omega_{\max}$ varies.

To validate the correctness and effectiveness of our algorithm, we perform a
series of experiments. The Python code is provided at
https://github.com/CausalML-Lab/PCMCI-Omega. In this section, we test four
algorithms111We did not conduct experiments on JIT-LiNGAM because this is from
a very recent paper Fujiwara et al. (2023) and is considered concurrent per
NeurIPS policy., $\text{PCMCI}_{\Omega}$, 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-valued
time series with exponential noise and binary-valued time series are in the
supplementary material.

Following Runge et al. (2019), we generate the continuous-valued time series
in three steps:

  1. 1.

Construct an $n$-variate time series $V$ with length $T$ using independent and
identical (Standard Gaussian or Exponential) noise temporarily. Determine
$\tau_{\max}$ and $\omega_{\max}$ where
$\omega_{\max}=\max\\{\omega_{j}\\}_{j\in\\{n\\}}$. After making sure that one
univariate time series, say $\mathbf{X}^{j}$, has periodicity $\omega_{\max}$,
the periodicity of the remaining time series $\mathbf{X}^{i},i\neq j$ is
randomly selected from $\\{1,\cdots,\omega_{\max}\\}$ respectively.

  2. 2.

Randomly generate $\omega_{j}$ binary edge matrices with shape
$(n,\tau_{\max})$ for each time series $\mathbf{X}^{j},j\in[n]$. $1$ denotes
an edge and $0$ denotes no edge. Each binary matrix represent one parent set
index $\mathrm{pInd}^{j}_{k},k\in[\omega_{j}]$. Randomly generate $\omega_{j}$
coefficient matrices with shape $(n,\tau_{\max})$ for each time series
$\mathbf{X}^{j},j\in[n]$. One binary edge matrix and one coefficient matrix
jointly determine one causal mechanism. Hence, total $\omega_{j}$ causal
mechanisms are constructed. Here, make sure that $V$ satisfies Assumption A6.

  3. 3.

Starting from time point $t>\tau_{\max}$, generate vector $\mathbf{X}_{t}$
over time according to all the causal mechanisms of $V$, until $t$ achieves
$T$.

Following the previous work in Huang et al. (2020),

$F_{1}$ score, Adjacency Precision, and Adjacency Recall are used to measure
the performance of the algorithms. The details of calculating these metrics
are described in 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.

A correct estimator $\hat{\omega}$ is the prerequisite for obtaining the
correct causal graph.

(a) Gaussian noise

(b) Performance when $\tau_{\max}$ or $n$ varies

Figure 3: 4 algorithms are tested on 5-variate time series. Set
$\tau_{\text{ub}}=15,\omega_{ub}=15$ for all variables. Every line corresponds
to a different algorithm. Every marker corresponds to the average performance
over 100 trials.

Fig.2(a)

shows the accuracy rate of $\hat{\omega}$ for different time lengths $T$.
Here, elements in $\\{N\omega_{j}\\}$ where
$N\in[\frac{\omega_{\text{ub}}}{\omega_{j}}]$ are all treated as correct
estimations. By Definition 2.2 and 2.3, the multiple of $\omega_{j}$ is still
associated with a correct causal graph. However, it leads to a finer time
point partition $\Pi^{j}$, decreasing the sample size used in each CI test
from approximately $T/\omega_{j}$ to approximately $T/(N\omega_{j})$. The
accuracy rate is sensitive to $\omega_{\max}$ for small $T$. This result
verifies that algorithm $\text{PCMCI}_{\Omega}$ has the capacity to detect the
true periodicity of each $\mathbf{X}^{j}\in V$ with a large enough time
length.

We evaluate the performances of PCMCIΩ on continuous-valued time series with
Gaussian noise shown as Fig.3(a). As $T$ increases, it is natural to see a
continuous improvement in performance.

The sub-figures show that all three evaluation metrics decrease when
$\omega_{\max}$ gets larger. The precision of PCMCIΩ is always far better than
other algorithms when $\omega_{\max}$ is not equal to 1. Given the fact that
the parent sets $\widehat{\mathrm{Pa}}(X^{j}_{t})\ \forall j,t$ obtained from
PCMCIΩ are subsets of the parent set $\widehat{S\mathrm{Pa}}(X^{j}_{t})\
\forall j,t$ estimated from PCMCI, the recall rate of PCMCI should be the
upper bound of the recall rate of PCMCIΩ. This assertion has been verified as
the red recall line of PCMCIΩ is always below the blue recall line of PCMCI as
$T$ increases.

In Fig.3(a), the recall of PCMCIΩ is worse than PCMCI for $T=500$. In this
regime, the accuracy rate of $\hat{\omega}$ is low, shown as the dark blue
line in Fig.2(a). Small sample sizes in CI tests may result in a sparser
causal graph. Hence the number of true positive edges may decrease. This is a
common problem for many constraint-based algorithms, but it hurts PCMCIΩ the
most because in PCMCIΩ, the sample sizes in each CI test are approximate
$T/\hat{\omega}$ instead of $T$. As $T$ increases, the red recall line of
PCMCIΩ 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.

We also observe the performance of our algorithm as $\tau_{\max}$ and $N$
varies in Fig.3(b). As the performance of $\text{PCMCI}_{\Omega}$ is
consistent over $n$-variate time series with different $n$, large
$\tau_{\max}$ may lead to a smaller precision and recall rate.

###  4.2 Case Study

Here, we construct an experiment with a 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^{\text{wp}}_{t}$ denote the monthly mean of surface pressure in the West
Pacific, $X^{\text{cp}}_{t}$, $X^{\text{ep}}_{t}$ and $X^{\text{ta}}_{t}$
denote the monthly mean of air temperature in the Central Pacific, East
Pacific, and tropical Atlantic, respectively.

The parent sets for each variable obtained from PCMCIΩ algorithms are shown in
Table 1. Sets of true and illusory parents of a variable at time $t$ are
separated by curly braces. For instance, variable $X^{\text{wp}}_{t}$ with
$\hat{\omega}_{\text{wp}}=1$ means that the causal mechanism of the surface
pressure in the West Pacific remains invariant over time with the estimated
parent set
$\\{X^{\text{wp}}_{t-1},X^{\text{wp}}_{t-2},X^{\text{ep}}_{t-1},X^{\text{ta}}_{t-1}\\}$.
Only time series $\mathbf{X}^{\text{cp}}$ 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^{\text{cp}}_{t}$ imply
that the causal effect from the tropical Atlantic air temperature
$X^{\text{ta}}_{t-1}$ to the Central Pacific air temperature
$X^{\text{cp}}_{t}$ would disappear every quarter of a year. Note that we do
not have a ground truth in this case, and we do not possess the necessary
knowledge in this area, so the significance of these results is under-
explored. More discussion about this application can be found in the
supplementary materials.

Table 1: Climate application results estimated from PCMCIΩ. PCMCIω
---
$X$ | $\hat{\omega}$ |  $\\{\widehat{\mathrm{Pa}}_{k}\\}_{k\in[\hat{\omega}]}$: true and illusory parent sets
$X^{\text{wp}}_{t}$ | 1 | $\\{X^{\text{wp}}_{t-1},X^{\text{wp}}_{t-2},X^{\text{ep}}_{t-1},X^{\text{ta}}_{t-1}\\}$
$X^{\text{cp}}_{t}$ | 3 |  $\\{X^{\text{cp}}_{t-1}\\}$; |  $\\{X^{\text{cp}}_{t-1},X^{\text{cp}}_{t-2},X^{\text{ta}}_{t-1}\\}$; | $\\{X^{\text{cp}}_{t-1},X^{\text{ta}}_{t-1}\\}$
$X^{\text{ep}}_{t}$ | 1 | $\\{X^{\text{ep}}_{t-1},X^{\text{ta}}_{t-1},X^{\text{ta}}_{t-2},X^{\text{cp}}_{t-1}\\}$
$X^{\text{ta}}_{t}$ | 1 | $\\{X^{\text{ta}}_{t-1},X^{\text{wp}}_{t-1}\\}$

##  5 Conclusions

In this paper, we propose a non-parametric, constraint-based causal discovery
algorithm $\text{PCMCI}_{\Omega}$ 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 on
continuous-valued and discrete-valued time series data. The algorithm
$\text{PCMCI}_{\Omega}$ has the capacity to reveal the existence of
periodicity of causal mechanisms in real-world datasets.

##  6 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.

## References

  * Pearl (1980) Pearl, J. Causality: models, reasoning, and inference. 1980.
  * Peters et al. (2017) Peters, J.; Janzing, D.; Schölkopf, B. _Elements of causal inference: foundations and learning algorithms_ ; The MIT Press, 2017.
  * Spirtes et al. (2000) Spirtes, P.; Glymour, C. N.; Scheines, R.; Heckerman, D. _Causation, prediction, and search_ ; MIT press, 2000.
  * Bergsma (2004) Bergsma, W. P. _Testing conditional independence for continuous random variables_ ; Citeseer, 2004.
  * Zhang et al. (2012) Zhang, K.; Peters, J.; Janzing, D.; Schölkopf, B. Kernel-based conditional independence test and application in causal discovery. _arXiv preprint arXiv:1202.3775_ 2012,
  * Shah and Peters (2020) Shah, R. D.; Peters, J. The hardness of conditional independence testing and the generalised covariance measure. 2020,
  * Runge et al. (2019) Runge, J.; Nowack, P.; Kretschmer, M.; Flaxman, S.; Sejdinovic, D. Detecting and quantifying causal associations in large nonlinear time series datasets. _Science advances_ 2019, _5_ , eaau4996.
  * Gerhardus and Runge (2020) Gerhardus, A.; Runge, J. High-recall causal discovery for autocorrelated time series with latent confounders. _Advances in Neural Information Processing Systems_ 2020, _33_ , 12615–12625.
  * Chu and Glymour (2008) Chu, T.; Glymour, C. Search for Additive Nonlinear Time Series Causal Models. _Journal of Machine Learning Research_ 2008, _9_.
  * Entner and Hoyer (2010) Entner, D.; Hoyer, P. O. On causal discovery from time series data using FCI. _Probabilistic graphical models_ 2010, 121–128.
  * Malinsky and Spirtes (2018) Malinsky, D.; Spirtes, P. Causal structure learning from multivariate time series in settings with unmeasured confounding. Proceedings of 2018 ACM SIGKDD workshop on causal discovery. 2018; pp 23–47.
  * Ikram et al. (2022) Ikram, A.; Chakraborty, S.; Mitra, S.; Saini, S.; Bagchi, S.; Kocaoglu, M. Root Cause Analysis of Failures in Microservices through Causal Discovery. _Advances in Neural Information Processing Systems_ 2022, _35_ , 31158–31170.
  * Hyvärinen et al. (2010) Hyvärinen, A.; Zhang, K.; Shimizu, S.; Hoyer, P. O. Estimation of a structural vector autoregression model using non-gaussianity. _Journal of Machine Learning Research_ 2010, _11_.
  * Peters et al. (2013) Peters, J.; Janzing, D.; Schölkopf, B. Causal inference on time series using restricted structural equation models. _Advances in neural information processing systems_ 2013, _26_.
  * Pamfil et al. (2020) Pamfil, R.; Sriwattanaworachai, N.; Desai, S.; Pilgerstorfer, P.; Georgatzis, K.; Beaumont, P.; Aragam, B. Dynotears: Structure learning from time-series data. International Conference on Artificial Intelligence and Statistics. 2020; pp 1595–1605.
  * Assaad et al. (2022) Assaad, C. K.; Devijver, E.; Gaussier, E. Survey and evaluation of causal discovery methods for time series. _Journal of Artificial Intelligence Research_ 2022, _73_ , 767–819.
  * Han et al. (2002) Han, F.; Subramanian, S.; Price, E. R.; Nadeau, J.; Strohl, K. P. Periodic breathing in the mouse. _Journal of Applied Physiology_ 2002, _92_ , 1133–1140.
  * Nakamura et al. (2003) Nakamura, A.; Fukuda, Y.; Kuwaki, T. Sleep apnea and effect of chemostimulation on breathing instability in mice. _Journal of Applied Physiology_ 2003, _94_ , 525–532.
  * Carskadon et al. (2005) Carskadon, M. A.; Dement, W. C.; others Normal human sleep: an overview. _Principles and practice of sleep medicine_ 2005, _4_ , 13–23.
  * Komarzynski et al. (2018) Komarzynski, S.; Huang, Q.; Innominato, P. F.; Maurice, M.; Arbaud, A.; Beau, J.; Bouchahda, M.; Ulusakarya, A.; Beaumatin, N.; Breda, G.; others Relevance of a mobile internet platform for capturing inter-and intrasubject variabilities in circadian coordination during daily routine: pilot study. _Journal of Medical Internet Research_ 2018, _20_ , e204.
  * Gong et al. (2015) Gong, M.; Zhang, K.; Schoelkopf, B.; Tao, D.; Geiger, P. Discovering temporal causal relations from subsampled data. International Conference on Machine Learning. 2015; pp 1898–1906.
  * Huang et al. (2019) Huang, B.; Zhang, K.; Gong, M.; Glymour, C. Causal discovery and forecasting in nonstationary environments with state-space models. International conference on machine learning. 2019; pp 2901–2910.
  * Krich et al. (2020) Krich, C.; Runge, J.; Miralles, D. G.; Migliavacca, M.; Perez-Priego, O.; El-Madany, T.; Carrara, A.; Mahecha, M. D. Estimating causal networks in biosphere–atmosphere interaction with the PCMCI approach. _Biogeosciences_ 2020, _17_ , 1033–1061.
  * Qu et al. (2021) Qu, Y.; Montzka, C.; Vereecken, H. Causation discovery of weather and vegetation condition on global wildfire using the PCMCI Approach. 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS. 2021; pp 8644–8647.
  * Zou et al. (2022) Zou, L.; Zha, Y.; Diao, Y.; Tang, C.; Gu, W.; Shao, D. Coupling the Causal Inference and Informer Networks for Short-term Forecasting in Irrigation Water Usage. _Water Resources Management_ 2022, 1–23.
  * Menegozzo et al. (2020) Menegozzo, G.; Dall’Alba, D.; Fiorini, P. Causal interaction modeling on ultra-processed food manufacturing. 2020 IEEE 16th International Conference on Automation Science and Engineering (CASE). 2020; pp 200–205.
  * Peterson (2022) Peterson, S. Comparison of Lasso Granger and PCMCI for Causal Feature Selection in Multivariate Time Series. Ph.D. thesis, The University of Arizona, 2022.
  * Arvind et al. (2021) Arvind, D.; Maiya, S.; Sedeño, P. A. Identifying causal relationships in time-series data from a pair of wearable sensors. 2021 IEEE 17th International Conference on Wearable and Implantable Body Sensor Networks (BSN). 2021; pp 1–4.
  * Castri et al. (2022) Castri, L.; Mghames, S.; Hanheide, M.; Bellotto, N. Causal discovery of dynamic models for predicting human spatial interactions. International Conference on Social Robotics. 2022; pp 154–164.
  * Castri et al. (2023) Castri, L.; Mghames, S.; Hanheide, M.; Bellotto, N.; others Enhancing Causal Discovery from Robot Sensor Data in Dynamic Scenarios. 2023,
  * Malinsky and Spirtes (2019) Malinsky, D.; Spirtes, P. Learning the structure of a nonstationary vector autoregression. The 22nd International Conference on Artificial Intelligence and Statistics. 2019; pp 2986–2994.
  * Huang et al. (2020) Huang, B.; Zhang, K.; Zhang, J.; Ramsey, J. D.; Sanchez-Romero, R.; Glymour, C.; Schölkopf, B. Causal Discovery from Heterogeneous/Nonstationary Data. _J. Mach. Learn. Res._ 2020, _21_ , 1–53.
  * Fujiwara et al. (2023) Fujiwara, D.; Koyama, K.; Kiritoshi, K.; Okawachi, T.; Izumitani, T.; Shimizu, S. Causal Discovery for Non-stationary Non-linear Time Series Data Using Just-In-Time Modeling. 2nd Conference on Causal Learning and Reasoning. 2023.
  * Serfozo (2009) Serfozo, R. _Basics of applied stochastic processes_ ; Springer Science & Business Media, 2009.
  * Karlin (2014) Karlin, S. _A first course in stochastic processes_ ; Academic press, 2014.
  * Le et al. (2016) Le, T. D.; Hoang, T.; Li, J.; Liu, L.; Liu, H.; Hu, S. A fast PC algorithm for high dimensional causal discovery with multi-core PCs. _IEEE/ACM transactions on computational biology and bioinformatics_ 2016, _16_ , 1483–1495.
  * Bertsekas and Tsitsiklis (2008) Bertsekas, D.; Tsitsiklis, J. N. _Introduction to probability_ ; Athena Scientific, 2008; Vol. 1.

Appendix

##  Appendix A PCMCI Algorithm

The PCMCI algorithm is proposed by Runge et al. (2019), aiming to detect time-
lagged causal relations in a window causal graph. 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 where Assumption
A4-A5 should be satisfied. In the second stage, Momentary Conditional
Independence tests (MCI) are used to further remove the false positive edges
caused by autocorrelations in time series data. More specifically, these two
steps can be briefly formalized as follows:

  * •

PC1 in Algorithm A1: Condition selection stage. PC1 is a variant of the
skeleton-discovery part of the PC algorithm in a more robust version named
stable-PC Le et al. (2016). The goal in this stage is to obtain a superset of
the parents $\widehat{\mathrm{Pa}}(X^{j}_{t})$ for all variables
$X^{j\in[n]}_{t\in[\tau_{\max}+1,T]}\in\mathbf{V}$. Initialize
$\widehat{\mathrm{Pa}}(X^{j}_{t})=\\{X^{i}_{t-\tau}\\}_{i\in[n],\tau\in[\tau_{\max}]}$.
$\widehat{\mathrm{Pa}}(X^{j}_{t})$ will remove $X^{i}_{t-\tau}$ if

$\displaystyle X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|\widehat{\mathrm{Pa}}(X^{j}_{t})\backslash\\{X^{i}_{t-\tau}\\}\right.$
(17)

  * •

MCI in Algorithm A2: Causal discovery stage. In this stage, do MCI tests for
all variable pairs $(X^{i}_{t-\tau},X^{j}_{t})$ with $i,j\in[n]$ and time
delays $\tau\in[\tau_{\max}]$:

$\displaystyle
MCI(X^{i}_{t-\tau},X^{j}_{t}|\widehat{\mathrm{Pa}}(X^{j}_{t})\backslash\\{X^{i}_{t-\tau}\\},\widehat{\mathrm{Pa}}(X^{i}_{t-\tau}))$
(18)

where $\widehat{\mathrm{Pa}}(X^{j}_{t})$ and
$\widehat{\mathrm{Pa}}(X^{i}_{t-\tau})$ are estimated from the PC1 stage.

Note that $\tau_{\max}$ in this section is the same as $\tau_{\text{ub}}$ in
the main paper, serving as the upper bound for the time lag that exhibits
causal effects. On the other hand, $\tau_{\max}$ in the main paper denotes the
maximum time lag observed within the multivariate time series. Essentially, in
the main paper, $\tau_{\text{ub}}$ is a parameter that must be fed into the
algorithm, and $\tau_{\max}$ is observed from the true causal graph. As a
default, we assume $\tau_{\text{ub}}$ is configured with a value greater than
$\tau_{\max}$, ensuring that the algorithm uncovers the correct causal
relations. See Fig.4 for more detail.

Figure 4: Set $\tau_{\text{ub}}$ to be 5, then all parent candidates of
variables at $t=15$ are included in the large orange box, ranging from $t=10$
to $t=14$. Consequently, the algorithm will only examine causal effects with a
time lag not exceeding 5. In the causal graph, $\tau_{\max}$ 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
$\tau_{1}=2,\tau_{2}=3,\tau_{3}=1$, respectively, and $\tau_{\max}$ represents
the largest value among these three maximum lags. Algorithm A1 $PC_{q_{\max}}$

1:Input: A $n$-variate time series
$V=(\mathbf{X}^{1},\mathbf{X}^{2},\mathbf{X}^{3},...,\mathbf{X}^{n})$, target
time series $\mathbf{X}^{j}$, maximum time lag $\tau_{\max}$, significance
threshold $\alpha_{PC}$, maximum condition dimension $p_{\max}$ (default
$p_{\max}=n\tau_{\max}$), maximum number of combinations $q_{\max}$ (default
$q_{\max}$ = 1), conditional independence test function $CI$

2:function ${CI(X,Y,\mathbf{Z})}$

3: Test $X\perp\\!\\!\\!\perp Y|\mathbf{Z}$ using test statistic measure $I$

4: return $p\text{-value}$, test statistic value $I$

5:Initialize preliminary set of parents
$\widehat{Pa}(X^{j}_{t})=\\{X^{i}_{t-\tau}:i\in\\{1,...,n\\},\tau\in\\{1,...,\tau_{\max}\\}\\}$

6:Initialize dictionary of test statistic values
$I^{\min}(X^{i}_{t-\tau}\rightarrow X^{j}_{t})=\infty\;\forall
X^{i}_{t-\tau}\in\widehat{Pa}(X^{j}_{t})$

7:for $p=0,1,2,...,p_{\max}$ do

8: if $|\widehat{Pa}(X^{j}_{t})|-1<p$ then

9: Break for-loop

10: end if

11: for all $X^{i}_{t-\tau}$ in $\widehat{Pa}(X^{j}_{t})$ do

12: $q=-1$

13: for all lexicographically chosen subsets
$\mathcal{S}\subseteq\widehat{Pa}(X^{j}_{t})\setminus\\{X^{i}_{t-\tau}\\}$
with $|\mathcal{S}|=p$ do

14: $q=q+1$

15: if $q\geq q_{\max}$ then

16: Break from inner for-loop

17: end if

18: Run CI test to obtain $(p\text{-value},I)\leftarrow
CI(X^{i}_{t-\tau},X^{j}_{t},\mathcal{S})$

19: if $|I|<I^{\min}(X^{i}_{t-\tau}\rightarrow X^{j}_{t})$
then$\;\;\;\;\;\;\;\;\;\;\;\;\;\rhd\text{Store min. $I$ of parent among all
tests}$

20: $I^{\min}(X^{i}_{t-\tau}\rightarrow X^{j}_{t})=|I|$

21: end if

22: if $p\text{-value}>\alpha_{PC}$
then$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\rhd\text{Removed only after all
$X^{i}_{t-\tau}$ have been tested}$

23: Mark $X^{i}_{t-\tau}$ for removal from $\widehat{Pa}(X^{j}_{t})$

24: Break from inner for-loop

25: end if

26: end for

27: end for

28: Remove non-significant parents from $\widehat{Pa}(X^{j}_{t})$

29: Sort parents in $\widehat{Pa}(X^{j}_{t})$ by
$I^{\min}(X^{i}_{t-\tau}\rightarrow X^{j}_{t})$ from largest to smallest

30:end for

31:return $\widehat{Pa}(X^{j}_{t})$

Algorithm A2 $MCI$

1:Input: A $n$-variate time series
$V=(\mathbf{X}^{1},\mathbf{X}^{2},\mathbf{X}^{3},...,\mathbf{X}^{n})$, sorted
parents $\widehat{Pa}(X^{j}_{t})$ for all variables $X^{j}$ estimated with
Algorithm A1, maximum time lag $\tau_{\max}$, maximum number $p_{X}$ of
parents of variable $X^{i}$, and conditional independence test function $CI$

2:for all $(X^{i}_{t-\tau},X^{j}_{t})$ with
$i,j\in\\{1,...,n\\},\tau\in\\{0,...,\tau_{\max}\\}$, excluding
$(X^{j}_{t},X^{j}_{t})$ do

3: Remove $X^{i}_{t-\tau}$ from $\widehat{Pa}(X^{j}_{t})$ if necessary

4: Define $\widehat{Pa}_{p_{X}}(X^{i}_{t-\tau})$ as the first $p_{X}$ parents
from $\widehat{Pa}(X^{i}_{t})$, shifted by $\tau$

5: Run MCI test to obtain $(p\text{-value},I)\leftarrow
CI(X^{i}_{t-\tau},X^{j}_{t},\mathbf{Z}=\\{\widehat{Pa}(X^{j}_{t}),\widehat{Pa}_{p_{X}}(X^{i}_{t-\tau})\\})$

6:end for

7:Optionally adjust $p\text{-value}$ of all links by False Discovery Rate-
approach (FDR)

8:return $p\text{-value}$ and MCI test statistic values

##  Appendix B $\text{PCMCI}_{\Omega}$

For simplicity’s sake, define sets: $[b]:=\\{1,2,...,b\\}$ and
$[a,b]:=\\{a,a+1,...,b\\}$, where $a,b\in\mathbb{N}$.

Algorithm B1 $PCMCI_{\Omega}$

1:Input: A $n$-variate time series
$V=(\mathbf{X}^{1},\mathbf{X}^{2},\mathbf{X}^{3},...,\mathbf{X}^{n})$,
periodicity upper bound $\omega_{\text{ub}}$, time lag upper bound
$\tau_{\text{ub}}$. By default, we assume $\tau_{\text{ub}}$ and
$\omega_{\text{ub}}$ are larger than their true value.

2:A superset of parent set is obtained using PCMCI with $\tau_{\text{ub}}$ and
denote it by $\widehat{S\mathrm{Pa}}(X^{j}_{t})\ \forall j,t$.

3:for $\mathbf{X}^{j}$ where $j\in[n]$ do

4: for a guess $\omega\in[\omega_{ub}]$ of $\omega_{j}$ do

5: Let $\widehat{\Pi}^{j}:=\\{\widehat{\Pi}^{j}_{k}|k\in[\omega]\\}$ where
$\widehat{\Pi}^{j}_{k}=\\{2\tau_{\text{ub}}+k,2\tau_{\text{ub}}+\omega+k,2\tau_{\text{ub}}+2\omega+k,\cdots\\}$.

6: for $k\in[\omega]$ do

7: Initialize the parent set for $X^{j}_{t},t\in\\{t:t\geq
2\tau_{\text{ub}},t\in\widehat{\Pi}^{j}_{k}\\}$ (with guess $\omega$) denoted
by
$\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})\leftarrow\widehat{S\mathrm{Pa}}(X^{j}_{t})$.

8: Consider $X^{i}_{t-\tau}\in\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$.
Remove $X^{i}_{t-\tau}$ from $\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$ if
$X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\mid\left(\widehat{S\mathrm{Pa}}(X^{j}_{t})\cup\widehat{S\mathrm{Pa}}(X^{i}_{t-\tau})\right)\setminus
X^{i}_{t-\tau}$ using a CI Test with samples $t\in\\{t:t\geq
2\tau_{\text{ub}},t\in\widehat{\Pi}^{j}_{k}\\}$.

9: Store $\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t})$ for
$X^{j}_{t},t\in\\{t:t\geq 2\tau_{\text{ub}},t\in\widehat{\Pi}^{j}_{k}\\}$.

10: end for

11: end for

12:

13: if there exists turning points $S_{j},\;S_{j}\in[\omega_{\text{ub}}]$ then

14: $\widehat{\omega}_{j}\leftarrow\min S_{j}$

15: else

16:

17:
$\widehat{\omega}_{j}\leftarrow\arg\min_{\omega\in[\omega_{\text{ub}}]}\max_{k\in[\omega]}|\widehat{\mathrm{Pa}}_{\omega}(X^{j}_{t\in\widehat{\Pi}^{j}_{k}})|$.

18: end if

19: Set
$\widehat{\mathrm{Pa}}(X^{j}_{t})\leftarrow\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})$
for $X^{j}_{t},t\in\\{t:t\geq 2\tau_{\text{ub}}\\}$.

20:end for

21:return $\hat{\omega}_{j}$ and $\widehat{\mathrm{Pa}}(X^{j}_{t})\ \forall
j\in[n],t\geq 2\tau_{\text{ub}}$.

##  Appendix C Soundness of PCMCIΩ

Figure 5: Partial causal graph for 3-variate time series
$V=\\{\mathbf{X}^{1},\mathbf{X}^{2},\mathbf{X}^{3}\\}$ with a Semi-Stationary
SCM where $\tau_{\max}=3$, $\omega_{1}=3$, $\omega_{2}=2$, $\omega_{3}=1$,
$\Omega=6$ and $\delta=6$. The first 3(=$\tau_{\max}$) time slices
$\\{\mathbf{X}_{t}\\}_{1\leq t\leq 3}$ are the starting points. The same color
edges denote the same causal mechanism. E.g. for $\mathbf{X}^{1}$: there are 3
($=\omega_{j}$) time partition subsets $\\{\Pi^{1}_{k}\\}_{1\leq k\leq 3}$.
The time points $t$ of nodes $X^{1}_{t}$ sharing the same filling color are in
the same time partition subsets. The time points $t$ of nodes $X^{1}_{t}$
sharing both the same filling color and the same outline shape are in the same
homogenous time partition subsets. There are 6 ($=\delta$) different Markov
chains in this multivariate time series $V$, and the first element of these 6
Markov chains is shown as $\\{Z^{q}_{1}\\}_{1\leq q\leq 6}$ and are tinted
with a gradient of blue hues. $Z^{1}_{1}$ and $Z^{1}_{2}$ denote the first two
elements of the first Markov chain while $Z^{2}_{1}$ and $Z^{2}_{2}$ denote
the first two elements of the second Markov chain.

###  C.1 Stationary Markov Chain

Claim: Any discrete-valued time series $V$ with Semi-Stationary Structural
Causal Model (SCM) satisfying assumption A1, A2, A4, A5 can be written as a
Markov chain $\\{Z_{n}\\}$ as long as this Markov chain satisfies
$\mathrm{Pa}(Z_{n})\subset Z_{n}\cup Z_{n-1}$ for all $n$, where $Z_{n}$ is a
set of variables in $V$. This Markov chain has a finite number of states if
all time series in $V$ are discrete-valued time series.

Note that when the notation $n$ is related to a Markov chain $Z_{n}$, it means
the running index. In the context of $X^{j\in[n]}_{t}$, $n$ represents the
index of component time series within the $n$-variate time series.

To simplify, assume that one associated Markov chain of
$V=\\{\mathbf{X},\mathbf{Y}\\}$ has $Z_{n}=\\{X_{t},Y_{t},X_{t-1},Y_{t-1}\\}$
with $t\in\\{t\in\mathbb{N}^{+}:,t\leq T\\}$ satisfying
$\mathrm{Pa}(Z_{n})\subset Z_{n}\cup Z_{n-1}$. Here, the notation for the time
points of variables is simplified as $t$ and $t-1$, even though it should be a
function of $n$, the running index of the Markov chain. Note that
$Z_{n-1}=\\{X_{t-2},Y_{t-2},X_{t-3},Y_{t-3}\\}$ rather than
$\\{X_{t-1},Y_{t-1},X_{t-2},Y_{t-2}\\}$, as the simplified notation could
erroneously suggest the latter sequence. A simple proof is shown below through
Markov assumption (A2).

###### Proof.

$\displaystyle p(Z_{n}|Z_{n-1},Z_{n-2},...)$ (19)
$\displaystyle=p(X_{t},Y_{t},X_{t-1},Y_{t-1}|Z_{n-1},Z_{n-2},...)$ (20)
$\displaystyle=p(X_{t}|Z_{n}\cup Z_{n-1}\setminus
X_{t},Z_{n-2},\cdots)p\bigg{(}Y_{t}|Z_{n}\cup Z_{n-1}\setminus(X_{t}\cup
Y_{t}),Z_{n-2},\cdots\bigg{)}\cdots$ (21)
$\displaystyle=p\bigg{(}X_{t}|\mathrm{Pa}(X_{t}),Z_{n}\cup
Z_{n-1}\setminus(X_{t}\cup\mathrm{Pa}(X_{t}))\bigg{)}$ (22)
$\displaystyle\;\;\;\;\times p\bigg{(}Y_{t}|\mathrm{Pa}(Y_{t}),Z_{n}\cup
Z_{n-1}\setminus(X_{t}\cup Y_{t}\cup\mathrm{Pa}(Y_{t}))\bigg{)}$
$\displaystyle\;\;\;\;\times p\bigg{(}X_{t-1}|\mathrm{Pa}(X_{t-1}),Z_{n}\cup
Z_{n-1}\setminus(X_{t}\cup Y_{t}\cup
X_{t-1}\cup\mathrm{Pa}(X_{t-1}))\bigg{)}\cdots$
$\displaystyle=p(X_{t}|Z_{n}\cup Z_{n-1}\setminus
X_{t})p\bigg{(}Y_{t}|Z_{n}\cup Z_{n-1}\setminus(X_{t}\cup
Y_{t})\bigg{)}\cdots$ (23)
$\displaystyle=p(X_{t},Y_{t},X_{t-1},Y_{t-1}|Z_{n-1})$ (24)
$\displaystyle=p(Z_{n}|Z_{n-1})$ (25)

∎

Assume that the space of both $X_{t}$ and $Y_{t}$ with $t<T$ are $\\{1,2\\}$.
There are total $2^{4}=16$ states of Markov Chain
$\\{Z_{n}\\}=\\{\\{X_{t},Y_{t},X_{t-1},Y_{t-1}\\}\\}$. The transition
probability $\mathbf{P}$ for this Markov Chain is illustrated as a $16\times
16$ matrix:

$\mathbf{P}=\blockarray{ccccc}&\mbox{\scriptsize$(1,1,1,1)$
}\mbox{\scriptsize$(2,1,1,1)$ }\mbox{\scriptsize$\cdots$}\\\
\block{c[ccc]c}\mbox{\scriptsize$(1,1,1,1)$}\text{$p_{1,1}$}\text{$p_{1,2}$}\text{$\cdots$}\\\
\mbox{\scriptsize$(2,1,1,1)$}\text{$p_{2,1}$}\text{$p_{2,2}$}\text{$\cdots$}\\\
\mbox{\scriptsize$\cdots$}\text{$\cdots$}\text{$\cdots$}\text{$\cdots$}{\tiny\text{16$\times$16}}\\\
$

where $(1,1,1,1)$ means $X_{t}=1,Y_{t}=1,X_{t-1}=1,Y_{t-1}=1$. Each row in
this transition probability matrix is a conditional distribution of $Z_{n}$
given one realization of $Z_{n-1}$. Each entry is a probability of having one
specific realization of $Z_{n}$ given one realization of $Z_{n-1}$. This
probability can be decomposed by conditional distributions based on Markov
assumption (A2). Take $p_{1,1}$ as an example:

$\displaystyle p_{1,1}$
$\displaystyle=p(X_{t}=1,Y_{t}=1,X_{t-1}=1,Y_{t-1}=1|X_{t-2}=1,Y_{t-2}=1,X_{t-3}=1,Y_{t-3}=1)$
(26)
$\displaystyle=p(X_{t}=1|\mathrm{Pa}(X_{t}))p(Y_{t}=1|\mathrm{Pa}(Y_{t}))p(X_{t-1}=1|\mathrm{Pa}(X_{t-1}))p(Y_{t-1}=1|\mathrm{Pa}(Y_{t-1}))$
(27)

where $\mathrm{Pa}(.)$ here are realizations, not random variables.

For time series $V$ with Semi-Stationary SCM, there are (potentially) $\delta$
different Markov chains $\\{Z^{q}_{n}\\},q\in[\delta]$:

$\displaystyle Z^{q}_{n}$
$\displaystyle=\\{\mathbf{X}_{\tau_{\max}+q+(n-1)\delta},\mathbf{X}_{\tau_{\max}+q+1+(n-1)\delta},...,\mathbf{X}_{\tau_{\max}+q-1+n\delta}\\},$

where $n\in\\{n:n\in\mathbb{N}^{+},\tau_{\max}+q-1+n\delta\leq T\\}$,
$\delta=\lceil\frac{\tau_{\max}+1}{\Omega}\rceil\Omega$. As proved in the
claim, such a Markov chain exists as long as $\mathrm{Pa}(Z^{q}_{n})\subset
Z^{q}_{n}\cup Z^{q}_{n-1}$ for all $n$. The value of $\delta$ can guarantee
the existence of such Markov chain because $\delta$ is larger than
$\tau_{\max}+1$ and is a multiple of $\Omega$, that is, a multiple of all
$\\{\omega_{j}\\}_{j\in[n]}$. By doing so, $\mathrm{Pa}(Z^{q}_{n})\subset
Z^{q}_{n}\cup Z^{q}_{n-1}$ is satisfied; for any variable $X^{j}_{t}$, there
exists $q\in[\delta]$ and $n\in\mathbb{N}^{+}$ such that variable $X^{j}_{t}$
and its parent set $\mathrm{Pa}(X^{j}_{t})$ can be included in $Z^{q}_{n}$;
and the causal mechanism generating $Z^{q}_{n}$ is invariant for different
$n$. The state space of $\\{Z^{q}_{n}\\}$ is the set containing all possible
realizations of
$\\{\mathbf{X}_{\tau_{\max}+q+(i-1)+(n-1)\delta}\\}_{i\in[\delta],n\in\mathbb{N}}$.
The transition probabilities between the states are the product of associated
causal mechanisms based on Markov assumption (A2).

Determined by the starting slice $\mathbf{X}_{t}$ where
$\tau_{\max}<t\leq\tau_{\max}+\delta$, there should be $\delta$ potentially
different Markov chains $\\{Z^{q}_{n}\\}$ where $1\leq q\leq\delta$. To be
more specific, those Markov chains are:

$\displaystyle\text{Markov Chain 1: }Z^{1}_{n}$
$\displaystyle=\\{\mathbf{X}_{\tau_{\max}+1+(n-1)\delta},\mathbf{X}_{\tau_{\max}+2+(n-1)\delta},...,\mathbf{X}_{\tau_{\max}+n\delta}\\},$
(28)
$\displaystyle\;\;\;\;\;\text{where}\;n\in\\{n:n\in\mathbb{N}^{+},\tau_{\max}+n\delta\leq
T\\}.$ $\displaystyle\text{Markov Chain 2: }Z^{2}_{n}$
$\displaystyle=\\{\mathbf{X}_{\tau_{\max}+2+(n-1)\delta},\mathbf{X}_{\tau_{\max}+3+(n-1)\delta},...,\mathbf{X}_{\tau_{\max}+1+n\delta}\\},$
(29)
$\displaystyle\;\;\;\;\;\text{where}\;n\in\\{n:n\in\mathbb{N}^{+},\tau_{\max}+1+n\delta\leq
T\\}.$ $\displaystyle\vdots$ $\displaystyle\text{Markov Chain $\delta$:
}Z^{\delta}_{n}$
$\displaystyle=\\{\mathbf{X}_{\tau_{\max}+n\delta},\mathbf{X}_{\tau_{\max}+1+n\delta},...,\mathbf{X}_{\tau_{\max}-1+(n+1)\delta}\\},$
(30)
$\displaystyle\;\;\;\;\;\text{where}\;n\in\\{n:n\in\mathbb{N}^{+},\tau_{\max}-1+(n+1)\delta\leq
T\\}.$

Given Irreducible and Aperiodic Markov Chain assumption (A7), discrete-time
Markov chain $\\{Z^{q}_{n}\\}_{0<n},\;q\in[\delta]$ with finite states should
be a stationary and ergodic Markov chain, and there is a unique stationary
distribution $\pi_{q}$ (Bertsekas and Tsitsiklis (2008), Karlin (2014)).
Additionally, the large power of the associate transition matrix
$\mathbf{P}_{q}$ will eventually converge to a matrix in which each row is the
stationary distribution $\pi_{q}$. Equivalently,

$\displaystyle\lim_{n\rightarrow\infty}p(Z^{q}_{n}=a|Z^{q}_{1}=b)=p(Z^{q}_{n}=a),\forall
a,b\in S.$ (31)

where $S$ is the state space of $Z^{q}_{n}$.

In other words, after a sufficiently long time, equivalently, $n$ is large
enough, the distribution of $\\{Z^{q}_{n}\\}$ does not change with increasing
$n$. That is, for large enough $n$:

$\displaystyle p(Z^{q}_{n_{1}})=p(Z^{q}_{n_{2}}),\forall n_{1},n_{2}>n.$ (32)

Returning from the stationary and ergodic Markov chains
$\\{Z^{q}_{n}\\},\;q\in[\delta]$ back to the original data $V$ through Eq.(28)
to Eq.(30), the distribution of the original data $V$ must adhere to the
following condition:

$\displaystyle
p(\mathbf{X}_{\tau_{\max}+q+n_{1}\delta},\mathbf{X}_{\tau_{\max}+q+1+n_{1}\delta},...,\mathbf{X}_{\tau_{\max}+q+\delta-1+n_{1}\delta})$
$\displaystyle=p(\mathbf{X}_{\tau_{\max}+q+n_{2}\delta},\mathbf{X}_{\tau_{\max}+q+1+n_{2}\delta},...,\mathbf{X}_{\tau_{\max}+q+\delta-1+n_{2}\delta})$
(33)

for any $q\in[\delta]$ and $n_{1},n_{2}>n$.

Given these clarifications, we can naturally introduce a more refined time
partition that is based on, yet finer than, the time partition defined in
Definition 2.3 in the main paper.

######  Definition C.1 (Homogenous Time Partition).

For a univariate time series $\mathbf{X}^{j}$ in a Semi-Stationary SCM with
periodicity $\omega_{j}$, the time partition $\Pi^{j}_{k}$ of $\mathbf{X}^{j}$
can be further divided into a series of non-overlapping and non-empty subsets
$\\{\pi^{j}_{(k,s)}\\}_{1\leq s\leq\frac{\delta}{\omega_{j}}}$. For each
$t\in[\tau_{\max}+1,T]$, there exists $k\in[\omega_{j}]\text{ so that
}t\in\Pi^{j}_{k}$ and further there exists
$s\in[\frac{\delta}{\omega_{j}}]\text{ so that }t\in\pi^{j}_{(k,s)}$.
$\pi^{j}_{(k,s)}$ can be written as:

$\displaystyle\pi^{j}_{(k,s)}:=\\{t:\tau_{\max}+1\leq t\leq
T,(t\bmod\omega_{j})+1=k,(t\bmod\frac{\delta}{\omega_{j}})+1=s\\}.$ (34)

With this definition, we have
$\cup_{s=1}^{\frac{\delta}{\omega_{j}}}\pi^{j}_{(k,s)}=\Pi^{j}_{k}$. While
time partition $\Pi^{j}_{k}$ guarantees that all variables in
$\\{X^{j}_{t}\\}_{t\in\Pi^{j}_{k}}$ share the same causal mechanism,
homogenous time partition $\pi^{j}_{(k,s)}$ guarantees that all variables in
$\\{X^{j}_{t}\\}_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}$ share the same
distribution where $t^{\prime}$ represent the steps needed by the associated
Markov chain to achieve equilibrium.

Fig.5 shows a partial causal graph for a 3-variate time series with Semi-
Stationary SCM. $\tau_{\max}=3$ means that the causal mechanisms start from
$t=4$, and the random variables with $t\in\\{1,2,3\\}$ are random noises. For
the first time series $\mathbf{X}^{1}$, the periodicity $\omega_{1}$ is 3. And
the periodicity of the time series $\mathbf{X}^{2}$ and $\mathbf{X}^{3}$ is 2
and 1, respectively. The periodicity of the whole time series $V$ is obtained
by $\text{LCM}(3,2,1)=6$.
$\delta=\lceil\frac{\tau_{\max}+1}{\Omega}\rceil\Omega=\lceil\frac{3+1}{6}\rceil\times
6=1\times 6=6$.

In Fig.5, periodicity $\omega_{1}=3$ means that the causal mechanisms repeat
every three time points and hence there are three time partition subsets
$\Pi^{1}_{k},k\in[3]$. More specifically,
$\Pi^{1}_{1}=\\{4,7,10,13,...,4+3N,...\\},\Pi^{1}_{2}=\\{5,8,11,14,...,5+3N,...\\},\Pi^{1}_{3}=\\{6,9,12,15,...,6+3N,...\\}$
where $N\in\mathbb{N}^{+}$. Random variables $\\{X^{1}_{t}\\}$ with $t$ in the
same time partition subset share the same causal mechanism. However, they may
not share the same marginal distribution.

Still in Fig.5, based on the definition of homogenous time partition, time
partition subset $\Pi^{1}_{1}$ for $\mathbf{X}^{1}$ can be further decomposed
as $\pi^{1}_{(k=1,s=1)}=\\{4,10,...,4+\delta
N,...\\},\pi^{1}_{(k=1,s=2)}=\\{7,13,...,7+\delta N,...\\}.$ where
$s\in[\frac{\delta}{\omega_{1}}]$. After a long run $n$, $Z^{1}_{n}$ and
$Z^{1}_{n+1}$ will eventually share the same distribution, that is, all the
variables inside $Z^{q}_{n}$ will share the same joint or marginal
distribution as the corresponding variables inside $Z^{q}_{n+1}$. To
illustrate this, we assume that this Markov chain has already achieved its
equilibrium at time point $t=4$. Based on Eq.(28) and Eq.(33), we have:

$\displaystyle
p(\mathbf{X}_{4},\mathbf{X}_{5},...,\mathbf{X}_{9})=p(\mathbf{X}_{10},\mathbf{X}_{11},...,\mathbf{X}_{15})=p(\mathbf{X}_{16},\mathbf{X}_{17},...,\mathbf{X}_{21})=\cdots$
(35)

From the identical joint distribution, we can further have:

$\displaystyle p(X^{1}_{4})=p(X^{1}_{10})=p(X^{1}_{16})=\cdots$ (36)

as $X^{1}_{4}\in\mathbf{X}_{4}$, $X^{1}_{10}\in\mathbf{X}_{10}$ and
$X^{1}_{16}\in\mathbf{X}_{16}$.

Therefore, for sufficiently large values of $t$ ensuring that ${Z^{1}_{n}}$
has reached its stationary distribution, all variables within
$\\{X^{j}_{t}\\}_{t\in\pi^{j}{(k,s)}}$ will share the same distribution.

In Fig.5, there are $6(=\delta)$ potentially different Markov chains
$\\{Z^{q}_{n}\\},q\in[\delta]$ in $V$. For any time window with length
$\delta$, $\\{\mathbf{X}_{t},...,\mathbf{X}_{t+\delta-1}\\}$, there exists
$q\in[\delta],n\in\mathbb{N}^{+}$, so that this time window can be completely
included in $Z^{q}_{n}$. For instance, set
$\\{\mathbf{X}_{5},\mathbf{X}_{6},...,\mathbf{X}_{10}\\}$ is in $Z^{2}_{1}$,
which is the first element of Markov chain $\\{Z^{2}_{n}\\}$.

Constructing Markov chains and applying the Irreducible and Aperiodic Markov
Chain assumption (A7) enable us to obtain a consistent estimator for the
conditional and joint distributions of interest.

###  C.2 Consistent Estimator

The conditional distributions for variables in
$\\{X^{j}_{t}\\}_{t\in\Pi_{k}^{j}}$ are the same, that is,
$p(x^{j}_{t_{1}}|\mathrm{Pa}(x_{t_{1}}^{j}))=p(x^{j}_{t_{2}}|\mathrm{Pa}(x_{t_{2}}^{j})),\;\forall
t_{1},t_{2}\in\Pi_{k}^{j}$. For simplicity, denote

$\displaystyle
p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})):=p(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})),\;\forall
t\in\Pi_{k}^{j}$ (37)

Consider an indicator function such that
$\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))=1$ if configuration
$(x_{t}^{j},\mathrm{Pa}(x_{t}^{j}))$ has realized, otherwise
$\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))=0$.

Since every $t\in\pi^{j}_{(k,s)}$ is apart from each other with $N\delta$
steps where $N\in\mathbb{N}^{+}$, and there must exist $q\in[\omega_{j}]$ and
$n_{1}\in\mathbb{N}^{+}$ so that $\\{x^{j}_{t},\mathrm{Pa}(x_{t}^{j})\\}\in
Z^{q}_{n_{1}}$, then for the same $q$, there must exist another $n_{2}$ so
that $\\{x^{j}_{t+N\delta},\mathrm{Pa}(x_{t+N\delta}^{j})\\}\in
Z^{q}_{n_{2}}$. Hence, we have
$\\{\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\\}_{t\in\pi^{j}_{(k,s)}}=\\{f(Z^{q(t)}_{n_{1}(t)})\\}_{t\in\pi^{j}_{(k,s)}}$
with some function $f:\mathbb{R}^{n\times\delta}\rightarrow\mathbb{R}^{1}$
satisfying $E|f(Z^{q(t)}_{n_{1}(t)})|<\infty$. Since the value of $t$
determines $q$ and $n_{1}$, we use $q(t)$ and $n_{1}(t)$ to emphasize their
relations. For large enough $t>t^{\prime}$,
$\\{\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\\}_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}$
are identical samples where $t^{\prime}$ is the time point needed by the
associate Markov chain to achieve its equilibrium after $n_{1}(t^{\prime})$
steps.

Without loss of generality, we assume $T$ is a multiple of $\delta$ all the
time.

We can construct an estimator of $p(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$ with
large enough $t$ as:

$\displaystyle\hat{p}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))=\frac{\delta}{T}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$
(38)

where $k,s$ is determined by $t$ and there must exist one and only one $k,s$
satisfying $t\in\pi^{j}_{(k,s)}$. Now, we are going to show this estimator is
consistent.

We first decompose the estimator into two parts: time point $t\leq t^{\prime}$
and $t>t^{\prime}$, where $t^{\prime}$ represents the time point when the
equilibrium of the associated Markov chain is achieved.

$\displaystyle\frac{\delta}{T}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$
(39) $\displaystyle=\frac{\delta}{T}\bigg{(}\sum_{t\leq
t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))+\sum_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\bigg{)}$
(40) $\displaystyle=\frac{\delta}{T}\sum_{t\leq
t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))+\frac{\delta}{T}\sum_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$
(41) $\displaystyle=\frac{\delta}{T}\sum_{t\leq
t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))+\frac{\delta}{T-t^{\prime}}\frac{T-t^{\prime}}{\delta}\frac{\delta}{T}\sum_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$
(42) $\displaystyle=\frac{\delta}{T}\sum_{t\leq
t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))+\frac{\delta}{T-t^{\prime}}\frac{T-t^{\prime}}{\delta}\frac{\delta}{T}\sum_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$
(43) $\displaystyle=\frac{\delta}{T}\sum_{t\leq
t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))+\frac{T-t^{\prime}}{T}\bigg{(}\frac{\delta}{T-t^{\prime}}\sum_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\bigg{)}$
(44)

Take a limit of Eq.(39), we have:

$\displaystyle\lim_{T\rightarrow\infty}\frac{\delta}{T}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$
(46) $\displaystyle=\lim_{T\rightarrow\infty}\frac{\delta}{T}\sum_{t\leq
t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))+\lim_{T\rightarrow\infty}\frac{T-t^{\prime}}{T}\bigg{(}\frac{\delta}{T-t^{\prime}}\sum_{t>t^{\prime},t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\bigg{)}$
(47)
$\displaystyle=0+\lim_{T\rightarrow\infty}\frac{T-t^{\prime}}{T}\bigg{(}\frac{1}{n_{1}(T)-n_{1}(t^{\prime})}\sum^{n_{1}(T)}_{n_{1}(t)>n_{1}(t^{\prime})}f(Z^{q(t)}_{n_{1}(t)})\bigg{)},\text{where
}t>t^{\prime},t\in\pi^{j}_{(k,s)}$ (48)
$\displaystyle\xlongequal{\text{Birkhoff's Ergodic
Theorem}}0+E\bigg{(}f(Z^{q(t)}_{n_{1}(t)})\bigg{)}$ (49)
$\displaystyle=E\bigg{(}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\bigg{)},\text{where
}t>t^{\prime},t\in\pi^{j}_{(k,s)}$ (50)
$\displaystyle=p(x^{j}_{t},\mathrm{Pa}(x_{t}^{j})),\text{where
}t>t^{\prime},t\in\pi^{j}_{(k,s)}$ (51)

Denote

$\displaystyle
p_{t\in\pi^{j}_{(k,s)}}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))\coloneqq
p(x^{j}_{t},\mathrm{Pa}(x_{t}^{j})),\text{where
}t>t^{\prime},t\in\pi^{j}_{(k,s)}$ (52)

Based on the definition of homogenous time partition and time partition,
$p_{t\in\pi^{j}_{(k,s)}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))=p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})),\;\forall
s\in[\frac{\delta}{\omega_{j}}]$.

Similar to Eq.(38), one estimator of
$p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})),\;\forall
k=[\omega_{j}]$ is

$\displaystyle\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
$\displaystyle=\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))}{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\mathrm{Pa}(x_{t}^{j}))}$
(53)
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))}{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(\mathrm{Pa}(x_{t}^{j}))}$
(54)
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\frac{T}{\delta}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))}{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\frac{T}{\delta}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(\mathrm{Pa}(x_{t}^{j}))}$
(55)

Take a limit of Eq.(53), we have:

$\displaystyle\lim_{T\rightarrow\infty}\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
(56) $\displaystyle\xlongequal{\text{Eq.\eqref{limit of conditional
distribution}}}\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))}{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\mathrm{Pa}(x_{t}^{j}))}$
(57)
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))p_{t\in\pi^{j}_{(k,s)}}(\mathrm{Pa}(x_{t}^{j}))}{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\mathrm{Pa}(x_{t}^{j}))}$
(58)
$\displaystyle\xlongequal{p_{t\in\pi^{j}_{(k,s)}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))\text{
are same for all
}s}\frac{p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\mathrm{Pa}(x_{t}^{j}))}{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\mathrm{Pa}(x_{t}^{j}))}$
(59) $\displaystyle=p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
(60)

Hence, $\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$ is a
consistent estimator of
$p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$.

Similarly, we construct an estimator of
$p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$ where $t\in[T]$:

$\displaystyle\hat{p}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$
$\displaystyle=\sum_{t}\mathbbm{1}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$
(61)
$\displaystyle=\frac{\sum_{t}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}.$
(62)

We will prove that this estimator is converged as $T$ goes to infinity in
Lemma D.2. Hence, it is a consistent estimator.

In this section, we have proved that
$\hat{p}(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$ in Eq.(38) is a consistent
estimator of $p(x^{j}_{t},\mathrm{Pa}(x_{t}^{j}))$ using samples with $t$ in
the same homogenous time partition subset and
$\hat{p}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$ in Eq.(53) is a consistent
estimator of $p(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$ using samples with $t$ in
the same time partition subset.

##  Appendix D Theorem

######  Theorem D.1.

Let $\widehat{\mathcal{G}}$ be the estimated graph using the Algorithm
$\text{PCMCI}_{\Omega}$. Under assumptions A1-A7 and with an oracle (infinite
sample size limit), we have that:

$\displaystyle\widehat{\mathcal{G}}=\mathcal{G}$ (63)

almost surely.

######  Lemma D.2.

Denote that $\\{\mathrm{Pa}_{k}(X^{j}_{t})\\}_{k\in[\omega_{j}]}$ contain the
true and illusory parent sets, where $\omega_{j}$ is the true periodicity of
$\mathbf{X}^{j}$. For any random variable $X_{t}^{j}$ with large enough $t$,
under assumptions A1-A7 and with an oracle (infinite sample size limit), we
have:

$\displaystyle
p\bigg{(}p(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))\neq
p(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\setminus
y)\bigg{)}=1,\;\forall y\in\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})$
(64)

Here,
$p(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))=\lim_{T\rightarrow\infty}\hat{p}(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))$.

###### Proof.

We first prove that there exist a sequence of coefficients
$\\{\alpha_{k}\\}_{k\in[\omega_{j}]}$ satisfying
$\sum_{k=1}^{\omega_{j}}\alpha_{k}=1$ so that:
$\forall\text{ configuration }\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})$,

$\displaystyle\hat{p}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))=\sum_{k=1}^{\omega_{j}}\alpha_{k}\hat{p}_{k}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
(65)

If this is correct, then
$\hat{p}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$ would be a consistent
estimator of $p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$. Based on
Eq.(62), we have:

$\displaystyle\hat{p}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$ (66)
$\displaystyle=\frac{\sum_{t}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(67)
$\displaystyle=\frac{\sum_{t}\sum_{k}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))\mathbbm{1}(t\in\Pi_{k}^{j})}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(68)
$\displaystyle=\sum_{k}\frac{\sum_{t}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))\mathbbm{1}(t\in\Pi_{k}^{j})}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(69)
$\displaystyle=\sum_{k}\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(70)
$\displaystyle=\sum_{k}\big{(}\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\cancel{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}}\frac{\cancel{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}}{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\big{)}$
(71)
$\displaystyle=\sum_{k}\big{(}\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{k}\mathrm{Pa}_{t\in\Pi_{k}^{j}}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\big{)}$
(72)
$\displaystyle=\sum_{k}\big{(}\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\big{)}$
(73)
$\displaystyle=\sum_{k}\big{(}\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}),\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus\mathrm{Pa}(x_{t}^{j}))\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\big{)}$
(74)
$\displaystyle=\sum_{k}\big{(}\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}\big{)}$
(75)
$\displaystyle=\sum_{k}\alpha_{k}(T)\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})),$
(76)
$\displaystyle\text{where}\;\alpha_{k}(T)=\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}.$
(77)

Using the same logic in Eq.(46)-(51), we can decompose the numerator and
denominator of $\alpha_{k}$ with homogenous time partition until each
component converges to a stationary distribution.

$\displaystyle\alpha_{k}(T)$
$\displaystyle=\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{k}\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(78)
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\frac{T}{\delta}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\frac{T}{\delta}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(79) $\displaystyle\lim_{T\rightarrow\infty}\alpha_{k}(T)$
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}$
(80)

Without loss of generality, assume
$y\in\mathrm{Pa}(x_{t}^{j}),\text{where}\;t\in\Pi_{j}^{1}$ and
$y\notin\mathrm{Pa}(x_{t}^{j}),\text{where}\;t\notin\Pi_{j}^{1}$ . Then we
have

$\displaystyle\hat{p}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)$
$\displaystyle=\frac{\sum_{t}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus y)}$ (81)
$\displaystyle=\sum_{k=2}^{\omega_{j}}\big{(}\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}\big{)}$ (82)
$\displaystyle+\frac{\sum_{t\in\Pi_{1}^{j}}\mathbbm{1}(x^{j}_{t},\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{t\in\Pi_{1}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}\frac{\sum_{t\in\Pi_{1}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus y)}$
$\displaystyle=\sum_{k=2}^{\omega_{j}}\beta_{k}(T)\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))+\beta_{1}(T)\hat{p}_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus
y)$ (83) $\displaystyle\text{where}\;\beta_{k}(T)$
$\displaystyle=\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{t}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus y)}$ (84)

Similarly, we have:

$\displaystyle\beta_{k}(T)$
$\displaystyle=\frac{\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{k}\sum_{t\in\Pi_{k}^{j}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}$ (85)
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\frac{T}{\delta}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}\frac{T}{\delta}\sum_{t\in\pi^{j}_{(k,s)}}\mathbbm{1}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}$ (86) $\displaystyle\lim_{T\rightarrow\infty}\beta_{k}(T)$
$\displaystyle=\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}$ (87)

Proving $p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))\neq
p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus y)$ is equal to
proving:

$\displaystyle
p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))-p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)\neq 0$ (88)

Substitutes Eq.(76) and Eq.(83) in Eq.(88), we have the following derivation:

$\displaystyle
p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))-p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)$ (89)
$\displaystyle=\lim_{T\rightarrow\infty}\bigg{(}\sum_{k=1}^{\omega_{j}}\alpha_{k}(T)\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))\bigg{)}-$
(90)
$\displaystyle\lim_{T\rightarrow\infty}\bigg{(}\sum_{k=2}^{\omega_{j}}\beta_{k}(T)\hat{p}_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))+\beta_{1}(T)\hat{p}_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus
y)\bigg{)}$
$\displaystyle=\sum_{k=1}^{\omega_{j}}\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))}p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
(91)
$\displaystyle\;\;\;\;-\bigg{(}\sum_{k=2}^{\omega_{j}}\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
$\displaystyle\;\;\;\;+\frac{\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\Pi_{(1,s)}^{j}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}{\sum_{k}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)}p_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus y)\bigg{)}$

After equating the denominators, the numerator is:

$\displaystyle\bigg{(}\sum_{k=1}^{\omega_{j}}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))\bigg{)}\bigg{(}\sum_{k=1}^{\omega_{j}}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)\bigg{)}$
$\displaystyle-\bigg{(}\sum_{k=2}^{\omega_{j}}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$
$\displaystyle+\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\Pi_{(1,s)}^{j}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)p_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus y)\bigg{)}$
$\displaystyle\times\bigg{(}\sum_{k=1}^{\omega_{j}}\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))\bigg{)}$
(92)

For the sake of simplicity, denote

$\displaystyle a_{k}$
$\displaystyle\coloneqq\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$
(93) $\displaystyle b_{k}$
$\displaystyle\coloneqq\sum_{s=1}^{\frac{\delta}{\omega_{j}}}p_{t\in\pi^{j}_{(k,s)}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
y)$ (94) $\displaystyle c_{k}$ $\displaystyle\coloneqq
p_{t\in\Pi_{k}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j}))$ (95) $\displaystyle
c_{1}^{\prime}$ $\displaystyle\coloneqq
p_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus y)$ (96)

After substituting the simple notations in Eq.(92):

$\displaystyle(\sum_{k=1}^{\omega_{j}}a_{k}c_{k})(\sum_{k=1}^{\omega_{j}}b_{k})-(\sum_{k=2}^{\omega_{j}}b_{k}c_{k}+b_{1}c_{1}^{\prime})(\sum_{k=1}^{\omega_{j}}a_{k})$
(97)
$\displaystyle=\sum_{k=1}^{\omega_{j}}(c_{k}-c_{1^{\prime}})a_{k}b_{1}+\sum_{k=1}^{\omega_{j}}\sum_{i>1,i\neq
k}^{\omega_{j}}(c_{k}-c_{i})a_{k}b_{i}$ (98)
$\displaystyle=b_{1}\sum_{k=1}^{\omega_{j}}c_{k}a_{k}-c_{1^{\prime}}b_{1}\sum_{k=1}^{\omega_{j}}a_{k}+\sum_{k=1}^{\omega_{j}}c_{k}a_{k}\sum_{i>1,i\neq
k}^{\omega_{j}}b_{i}-\sum_{k=1}^{\omega_{j}}a_{k}\sum_{i>1,i\neq
k}^{\omega_{j}}c_{i}b_{i}$ (99)

Define

$\displaystyle V_{t}=\\{\mathbf{X}_{t^{\prime}}|0<t^{\prime}<t\\}$ (100)

That is, $V_{t}$ contains all the nodes before time point $t$.

Denote $\\{b_{t_{i}}\\}_{i\in[n]}=\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})$ and
assume $\\{b_{t_{i}}\\}_{{1\leq i\leq n_{1}<n}}=\mathrm{Pa}_{1}(x_{t}^{j})$,
where $t\in\Pi_{(k,s)}^{j}$

We express $p_{t\in\Pi_{(k,s)}^{j}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$ by
marginalizing all other random variables occurring before the latest variables
in $\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})$ and utilizing the Causal Markov
assumption (A2):

$\displaystyle p_{t\in\Pi_{(k,s)}^{j}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))$
(101)
$\displaystyle=p(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})|t\in\Pi_{(k,s)}^{j})$
(102)
$\displaystyle=\sum_{V_{h}\setminus\\{b_{t_{i}}\\}_{i\in[n]}}p(b_{t_{1}},b_{t_{2}},...b_{t_{n}},V_{h}\setminus\\{b_{t_{i}}\\}_{i\in[n]}|h=\max\\{t_{i},1\leq
i\leq n\\},t\in\Pi_{(k,s)}^{j})$ (103)
$\displaystyle=\sum_{\\{\mathrm{Pa}(b_{t_{i}})\\}_{i\in[n]}}p(b_{t_{i}}|\mathrm{Pa}(b_{t_{i}}))\sum_{V_{\tau_{\max}}}\sum_{\tau_{\max}<t^{\prime}\leq
h}\sum_{j\in[n]}\sum_{x^{j}_{t^{\prime}},\mathrm{Pa}(x^{j}_{t^{\prime}})}p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}p(V_{\tau_{\max}})$
(104)

Note that $x^{j}_{t^{\prime}}\in V_{h}\setminus\\{b_{t_{i}}\\}_{i\in[n]}$.

This joint distribution is now represented by conditional distributions of one
related variable given its parents.

Similarly, assume $y=b_{t_{1}}$, we have

$\displaystyle
p_{t\in\Pi_{(k,s)}^{j}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus y)$ (105)
$\displaystyle=p_{t\in\Pi_{(k,s)}^{j}}(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus
b_{t_{1}})$ (106) $\displaystyle=p(\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}\setminus
y)|t\in\Pi_{(k,s)}^{j})$ (107)
$\displaystyle=\sum_{V_{t_{n}}\setminus\\{b_{t_{i}}\\}_{i\neq
1}}p(b_{t_{2}},...b_{t_{n-1}},b_{t_{n}},V_{t_{n}}\setminus\\{b_{t_{i}}\\}_{i\neq
1}|h=\max\\{t_{i},2\leq i\leq n\\})$ (108)
$\displaystyle=\sum_{\\{\mathrm{Pa}(b_{t_{i}})\\}_{i\neq
1}}p(b_{t_{i}}|\mathrm{Pa}(b_{t_{i}}))\sum_{V_{\tau_{\max}}}\sum_{\tau_{\max}<t^{\prime}\leq
h}\sum_{j\in[n]}\sum_{x^{j}_{t^{\prime}},\mathrm{Pa}(x^{j}_{t^{\prime}})}p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}p(V_{\tau_{\max}})$
(109)

Note that $x^{j}_{t^{\prime}}\in V_{t_{n}}\setminus\\{b_{t_{i}}\\}_{i\neq 1}$.

$p_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus y)$ can also be
represented by those conditional distributions based on Bayes rule.

$\displaystyle p_{t\in\Pi_{1}^{j}}(x^{j}_{t}|\mathrm{Pa}(x_{t}^{j})\setminus
b_{t_{1}})$ (110)
$\displaystyle=\frac{p_{t\in\Pi_{1}^{j}}(x^{j}_{t},b_{t_{2}},..,b_{t_{n_{1}}})}{p_{t\in\Pi_{1}^{j}}(b_{t_{2}},...,b_{t_{n_{1}}})}$
(111)
$\displaystyle=\frac{\sum_{b_{t_{1}}}p_{t\in\Pi_{1}^{j}}(x^{j}_{t},b_{t_{1}},...,b_{t_{n_{1}}})}{\sum_{b_{t_{1}}}p_{t\in\Pi_{1}^{j}}(b_{t_{1}},...,b_{t_{n_{1}}})}$
(112)
$\displaystyle=\frac{\sum_{b_{t_{1}}}p_{t\in\Pi_{1}^{j}}(x^{j}_{t},\mathrm{Pa}_{1}(x_{t}^{j}))}{\sum_{b_{t_{1}}}p_{t\in\Pi_{1}^{j}}(\mathrm{Pa}_{1}(x_{t}^{j}))}$
(113)
$\displaystyle=\frac{\sum_{b_{t_{1}}}p(x^{j}_{t}|\mathrm{Pa}_{1}(x_{t}^{j}))p_{t\in\Pi_{1}^{j}}(b_{t_{1}},...,b_{t_{n_{1}}})}{{\sum_{b_{t_{1}}}p_{t\in\Pi_{1}^{j}}(b_{t_{1}},...,b_{t_{n_{1}}})}}$
(114)
$\displaystyle=\frac{\sum_{b_{t_{1}}}p(x^{j}_{t}|\mathrm{Pa}_{1}(x_{t}^{j}))\sum_{V_{h}\setminus\\{b_{t_{i\in[n_{1}]}}\\}}p_{t\in\Pi_{1}^{j}}(b_{t_{1}},...b_{t_{n_{1}}},V_{h}\setminus\\{b_{t_{i\in[n_{1}]}}\\}|h=\max\\{t_{i\in[n_{1}]}\\})}{{\sum_{b_{t_{1}}}\sum_{V_{h}\setminus\\{b_{t_{i\in[n_{1}]}}\\}}p_{t\in\Pi_{1}^{j}}(b_{t_{1}},...b_{t_{n_{1}}},V_{h}\setminus\\{b_{t_{i\in[n_{1}]}}\\}|h=\max\\{t_{i\in[n_{1}]}\\})}}$
(115) $\displaystyle=\frac{\sum_{b_{t_{1}}}AB}{\sum_{b_{t_{1}}}CD}$ (116)

where

$\displaystyle
A=p(x^{j}_{t}|\mathrm{Pa}_{1}(x_{t}^{j}))\sum_{\\{\mathrm{Pa}(b_{t_{i}})\\}_{i\in[n_{1}]}}p(b_{t_{i}}|\mathrm{Pa}(b_{t_{i}}))$
(117) $\displaystyle B=\sum_{V_{\tau_{\max}}}\sum_{\tau_{\max}<t^{\prime}\leq
h}\sum_{j\in[n]}\sum_{x^{j}_{t^{\prime}},\mathrm{Pa}(x^{j}_{t^{\prime}})}p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}p(V_{\tau_{\max}})$
(118) $\displaystyle
C=\sum_{\\{\mathrm{Pa}(b_{t_{i}})\\}_{i\in[n_{1}]}}p(b_{t_{i}}|\mathrm{Pa}(b_{t_{i}}))$
(119) $\displaystyle D=\sum_{V_{\tau_{\max}}}\sum_{\tau_{\max}<t^{\prime}\leq
h}\sum_{j\in[n]}\sum_{x^{j}_{t^{\prime}},\mathrm{Pa}(x^{j}_{t^{\prime}})}p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}p(V_{\tau_{\max}})$
(120)

Note that $t\in\Pi_{1}^{j}$ for distributions in above section from Eq.(110)
to Eq.(120) and that $x^{j}_{t^{\prime}}\in
V_{h}\setminus\\{b_{t_{i}}\\}_{i\in[n_{1}]}$.

Hence, every term in Eq.(99) can be expressed as a function of those
conditional distributions. Substituting Eq.(104), Eq.(109) and Eq.(116) in
Eq.(99), we have a polynomial equation only composed of conditional
distributions
$\\{p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}\\}_{j\in[n],t^{\prime}\leq
t}$ except the joint distribution of the starting points $p(V_{\tau_{\max}})$.
Note that the conditional distributions of variables in
$\\{X^{j}_{t}\\}_{t\in\Pi^{j}_{k}},j\in[n],k\in[\omega_{j}]$ are the same.
Since sets do not allow duplicate values, set
$\\{p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}\\}_{j\in[n],t^{\prime}\leq
t}$ contains only different conditional distributions. There should be
potentially total $\sum_{j=1}^{n}\omega_{j}$ 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 $\sum_{y}p(x|y)=1$, all components in the set
$\\{p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}\\}_{j\in[n],t^{\prime}\leq
t}$ are mutually independent, and $p(V_{\tau_{\max}})$ is also independent of
all the causal mechanisms because the first starting points are random noises.
That is, upon adjustments, all the terms in Eq.(99) should be rendered
independent of each other, without any imposed constraints across them.

After expanding all the summations in Eq.(99), the coefficients of this
polynomial equation are either $1$ or $-1$. Each coefficient is accompanied by
one unique monomial as index $(k,s)$ in the joint distribution
$p_{t\in\pi^{j}_{k,s}}$ determined a unique product of conditional
distributions, i.e., with a different pair of $(k,s)$, the product should be
different. Considering all random and independent conditional distributions in
$\\{p\big{(}x^{j}_{t^{\prime}}|\mathrm{Pa}(x^{j}_{t^{\prime}})\big{)}\\}_{j\in[n],t^{\prime}\leq
t}$, the polynomial is not identically zero, and the probability of choosing a
root of this polynomial is zero.

Denote the polynomial equation in Eq.(99) as A, we have:

$\displaystyle p(A=0)=0$ (121)

Back to the original Eq.(89), we finally have
$p\Big{(}p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j}))\neq
p(x^{j}_{t}|\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})\setminus y)\Big{)}=1,\;\forall
y\in\cup_{h}\mathrm{Pa}_{h}(x_{t}^{j})$.

∎

######  Lemma D.3.

Let $\widehat{S\mathrm{Pa}}(\mathbf{X}^{j}_{t})$ denote the estimated superset
of parent set for $\mathbf{X}^{j}\in V$ obtained from the Algorithm B1 (line
2). $\\{\mathrm{Pa}_{k}(X^{j}_{t})\\}_{k\in[\omega_{j}]}$ contain the true and
illusory parent sets, where $\omega_{j}$ is the true periodicity of
$\mathbf{X}^{j}$. Under assumptions A1-A7 and with an oracle (infinite sample
size limit), we have:

$\displaystyle\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\subseteq\widehat{S\mathrm{Pa}}(X^{j}_{t}),\;\forall
t\in[\tau_{\max}+1,T]$

almost surely.

###### Proof.

Assume the contrary, i.e., there exists
$s\in\cup_{k}\mathrm{Pa}_{k}(X_{t}^{j})\setminus\widehat{S\mathrm{Pa}}(X^{j}_{t})$.
From Lemma D.2, we have
$X^{j}_{t}\cancel{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}s\left|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\setminus
s\right.$. By the Definition 2.4, we have
$\mathrm{Pa}(X^{j}_{t})\subset\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})$.
If $s\not\in\mathrm{Pa}(X^{j}_{t})$, by the causal Markov property (A2), the
dependence relation can not be true, because $s$ is a non-descendant of
$X^{j}_{t}$. If $s\in\mathrm{Pa}(X^{j}_{t})$, our Algorithm would have
concluded that
$X_{t}^{j}\cancel{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}s\left|\widehat{S\mathrm{Pa}}(X^{j}_{t})\right.$
(line 2) with a consistent CI test, evident from the causal Markov property,
contradicting our assumption. Hence, the lemma. ∎

######  Lemma D.4.

Let $\widehat{Pa}(X^{j}_{t})$ denote the estimated parent set for
$\mathbf{X}^{j}\in V$ obtained from the Algorithm B1 (line 19) assuming that
true $\omega_{j}$ has obtained (line 17).
$\\{\mathrm{Pa}_{k}(X^{j}_{t})\\}_{k\in[\omega_{j}]}$ contain the true and
illusory parent sets. Under assumptions A1-A7 and with an oracle (infinite
sample size limit), we have:

$\displaystyle\widehat{Pa}(X^{j}_{t})=Pa(X^{j}_{t}),\;\forall
t\in[\tau_{\max}+1,T]$ (122)

almost surely.

###### Proof.

From Lemma D.3,

$\displaystyle
Pa(X^{j}_{t})\subset\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})\subseteq\widehat{S\mathrm{Pa}}(X^{j}_{t}),\;\forall
t\in[\tau_{\max}+1,T],j\in[n]$ (123)

In Runge et al. (2019), the author proved
$\widehat{Pa}(X^{j}_{t})=Pa(X^{j}_{t})$ if we run PCMCI on stationary time
series. Using the same logic, we have the following proof.

Suppose $X^{i}_{t-\tau}\notin\widehat{Pa}(X^{j}_{t})$ but $X^{i}_{t-\tau}\in
Pa(X^{j}_{t})$. With a consistent conditional independence test and correct
time partition, the $MCI$ test (line 8 in Algorithm B1) will remove
$X^{i}_{t-\tau}$ from $\widehat{Pa}_{\omega_{j}}(X^{j}_{t})$ if and only if:

$\displaystyle X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|\widehat{SPa}(X^{j}_{t})\setminus\\{X^{i}_{t-\tau}\\},\widehat{SPa}(X^{i}_{t-\tau})\right.$
(124)

Based on Eq.(123), the rule is equivalent to removing $X^{i}_{t-\tau}$ from
$\widehat{Pa}_{\omega_{j}}(X^{j}_{t})$ if and only if:

$\displaystyle X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|\Bigg{\\{}Pa(X^{j}_{t})\setminus\\{X^{i}_{t-\tau}\\},Pa(X^{i}_{t-\tau})\right.,$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\widehat{SPa}(X^{j}_{t})\setminus(Pa(X^{j}_{t})\cup\\{X^{i}_{t-\tau}\\}),\widehat{SPa}(X^{i}_{t-\tau})\setminus
Pa(X^{i}_{t-\tau})\Bigg{\\}}$ (125) $\displaystyle\Rightarrow$ $\displaystyle
X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|Pa(X^{j}_{t})\setminus\\{X^{i}_{t-\tau}\\},Pa(X^{i}_{t-\tau})\right.$
(126)

Based on Causal Markov Condition assumption (A2) and Faithfulness Condition
(A3), from Eq.(126) we have $X^{i}_{t-\tau}\notin Pa(X^{j}_{t})$. In other
words, if $X^{i}_{t-\tau}\notin\widehat{Pa}(X^{j}_{t})$ then
$X^{i}_{t-\tau}\notin Pa(X^{j}_{t})$. That is,
$Pa(X^{j}_{t})\subseteq\widehat{Pa}(X^{j}_{t})$

Suppose $X^{i}_{t-\tau}\in\widehat{Pa}(X^{j}_{t})$ but $X^{i}_{t-\tau}\notin
Pa(X^{j}_{t})$. By the contraposition of Faithfulness (A1), we know that
$X^{i}_{t-\tau}\cancel{\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}}X^{j}_{t}\left|\widehat{Pa}(X^{j}_{t})\setminus\\{X^{i}_{t-\tau}\\},\widehat{Pa}(X^{i}_{t-\tau})\right.$.
Denote
$W=\big{\\{}\widehat{Pa}(X^{j}_{t})\setminus\\{Pa(X^{j}_{t}),X^{i}_{t-\tau}\\}\big{\\}}\cup\big{\\{}\widehat{Pa}(X^{i}_{t-\tau})\setminus
Pa(X^{i}_{t-\tau})\big{\\}}$. Since $X^{i}_{t-\tau}\notin Pa(X^{j}_{t})$,
based on Causal Markov Condition assumption (A2),

$\displaystyle W\cup X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|Pa(X^{j}_{t})\right.$ $\displaystyle\xRightarrow{}W\cup
X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|Pa(X^{j}_{t}),Pa(X^{i}_{t-\tau})\right.$
$\displaystyle\xRightarrow{\text{Weak Union}}X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|\\{Pa(X^{j}_{t}),Pa(X^{i}_{t-\tau})\\}\cup W\right.$
$\displaystyle\xRightarrow{}X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}\left|\widehat{Pa}(X^{j}_{t})\setminus\\{X^{i}_{t-\tau}\\},\widehat{Pa}(X^{i}_{t-\tau})\right.$

This is contrary to the assumption so that there is no such $X^{i}_{t-\tau}$
satisfying $X^{i}_{t-\tau}\in\widehat{Pa}(X^{j}_{t})$ but
$X^{i}_{t-\tau}\notin Pa(X^{j}_{t})$. In other words, if
$X^{i}_{t-\tau}\in\widehat{Pa}(X^{j}_{t})$, then $X^{i}_{t-\tau}\in
Pa(X^{j}_{t})$. That is, $\widehat{Pa}(X^{j}_{t})\subseteq Pa(X^{j}_{t})$.
Combined with the previous conclusion that
$Pa(X^{j}_{t})\subseteq\widehat{Pa}(X^{j}_{t})$, we have
$\widehat{Pa}(X^{j}_{t})=Pa(X^{j}_{t})$.

∎

Based on Lemma D.2, Lemma D.3 and Lemma D.4, we can identify the true
$\omega_{j}$ for $\mathbf{X}^{j}$ through Lemma D.5.

######  Lemma D.5.

Let $\omega_{j}$ denote the true periodicity for $\mathbf{X}^{j}\in V$ and
$\widehat{\mathrm{Pa}}(X^{j}_{t\in\Pi^{j}_{k}})$ denote the estimated parent
set for $X^{j}_{t}$ obtained from Algorithm B1 where $t\in\Pi^{j}_{k}$.
Define:

$\displaystyle\widehat{\omega}_{j}=\arg\min_{\omega\in[\omega_{\text{ub}}]}\max_{k\in[\omega]}|\widehat{\mathrm{Pa}}(X^{j}_{t\in\Pi^{j}_{k}})|$
(127)

Under assumptions A1-A7 and with an oracle (infinite sample limit), we have
that $\hat{\omega}_{j}=\omega_{j},\;\forall j\in[n]$ almost surely.

###### Proof.

Assume the contrary that $\hat{\omega}_{j}\neq\omega_{j}$, then in the
Algorithm B1, we have an incorrect time partition $\widehat{\Pi}^{j}$. Hence,
CI tests that are performed use samples with different causal mechanisms.
$\hat{p}(X^{j}_{t}|\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t}))$ in
Eq.(66) is estimated from a mixture of two or more time partition subsets, say
$\Pi^{j}_{1}$ and $\Pi^{j}_{2}$. We can apply Lemma D.2 where
$\cup_{k=1}^{\omega_{j}}\mathrm{Pa}_{k}(X^{j}_{t})$ is replaced by
$\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})$ and then in Lemma D.3,
$\widehat{S\mathrm{Pa}}(X^{j}_{t})$ is replaced by
$\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})$ and hence
$\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})\subseteq\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})$
where $\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})$ is obtained from
samples with $t$ from the mixture of two different partition subsets (line 8).
Hence, with $\hat{\omega}_{j}$,
$|\widehat{\mathrm{Pa}}_{\hat{\omega}_{j}}(X^{j}_{t})|\geq|\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})|$
using mixture samples $t\in\cup_{k=1}^{2}\Pi^{j}_{k}$. However, with true
$\omega_{j}$, we have
$|\widehat{\mathrm{Pa}}_{\omega_{j}}(X^{j}_{t})|=|\mathrm{Pa}(X^{j}_{t})|$
based on Lemma D.4. With Assumption A6 the Hard Mechanism Change,
$|\cup_{k=1}^{2}\mathrm{Pa}_{k}(X^{j}_{t})|>|\mathrm{Pa}(X^{j}_{t})|$ so that
$\omega_{j}$ always leads to a smaller size of estimated parent sets than
$\hat{\omega}_{j}$, contrary to the definition of $\hat{\omega}_{j}$. Hence,
$\hat{\omega}_{j}=\omega_{j}$. ∎

With those lemmas, we can prove Theorem 1.

###### Proof.

Assuming that a correct $\omega_{j}$ has already been obtained, from Lemma D.4
we have

$\displaystyle\widehat{Pa}(X^{j}_{t})=Pa(X^{j}_{t}),\;\forall t\geq
2\tau_{\text{ub}},j\in[n]$

From Lemma D.5, we know that a correct $\omega_{j}$ must be obtained with
consistent CI tests, that is, $\hat{\omega}_{j}=\omega_{j},\forall j\in[n]$.
Therefore from Algorithm B1, we have

$\displaystyle\widehat{Pa}(X^{j}_{t})=Pa(X^{j}_{t}),\;\forall t\geq
2\tau_{\text{ub}},j\in[n]$

If the causal mechanism is fixed across time, i.e., $\omega_{j}=1,j\in[n]$,
the proof of PCMCI Runge et al. (2019) showed that for all $\mathbf{X}^{j}\in
V$,

$\displaystyle X^{i}_{t-\tau}\rightarrow
X^{j}_{t}\notin\mathcal{G}\xRightarrow{}X^{i}_{t-\tau}\rightarrow
X^{j}_{t}\notin\widehat{\mathcal{G}}$ $\displaystyle X^{i}_{t-\tau}\rightarrow
X^{j}_{t}\in\mathcal{G}\xRightarrow{}X^{i}_{t-\tau}\rightarrow
X^{j}_{t}\in\widehat{\mathcal{G}}$

Therefore $\widehat{\mathcal{G}}=\mathcal{G}$.

If $\exists\omega_{j}>1$, we can simply separate the whole graph $\mathcal{G}$
into sub graphs $\\{\mathcal{G}^{\omega_{j}}_{k}\\}_{k\in[\omega_{j}]}$
consisting of only target variable $X^{j}_{t}$ with corresponding
$t\in\\{\Pi^{j}_{k}\\}_{k\in[\omega_{j}]}$ and parent variables
$X^{i}_{t^{\prime}}\in\mathrm{Pa}(X^{j}_{t})$. Focusing only on one time
partition subset $\Pi^{j}_{k},k\in[\omega_{j}]$, we have

$\displaystyle\widehat{\mathcal{G}}^{\omega_{j}}_{k}=\mathcal{G}^{\omega_{j}}_{k}$
(128)

for any $k\in[\omega_{j}]$ and $j\in[n]$ based on the proof of Proposition 1
in the supplementary materials of Runge et al. (2019).
Each sub-graph $\mathcal{G}^{\omega_{j}}_{k}$ includes only variable
$X^{j}_{t}$, the edges entering $X^{j}_{t}$ for time points $t\in\Pi^{j}_{k}$
and the corresponding parent variables $X^{i}_{t^{\prime}}\in Pa(X^{j}_{t})$.
Given $\Pi^{j}=\underset{k\in[\omega_{j}]}{\cup}\Pi^{j}_{k}$ and
$V=\underset{j\in[n]}{\cup}\mathbf{X}^{j}$, we have:

$\displaystyle\widehat{\mathcal{G}}=\underset{j\in[n],\;k\in[\omega_{j}]}{\cup}\widehat{\mathcal{G}}^{\omega_{j}}_{k}$
(129)
$\displaystyle\mathcal{G}=\underset{j\in[n],\;k\in[\omega_{j}]}{\cup}\mathcal{G}^{\omega_{j}}_{k}$
(130)

On the basis of Eq.(128), we finally have:

$\displaystyle\widehat{\mathcal{G}}=\mathcal{G}$

∎

##  Appendix E Turning Points

Figure 6: In the above illustration of the "turning point," the sizes of
parent sets for different estimates $\hat{\omega}_{j}$ are depicted as
$|\widehat{\mathrm{Pa}}_{k}(X^{j}_{t})|,k\in[\hat{\omega}_{j}]$. It is worth
noting that $\widehat{\mathrm{Pa}}_{k}(X^{j}_{t})$ represents either the true
parent set or the illusory parent set of $X^{j}_{t}$. In this context, we are
interested in the sizes of these parent sets. The first occurrence of the
"turning point" happens at $\hat{\omega}_{j}=3$ since the sizes of parent sets
obtained when $\hat{\omega}_{j}=2$ and $\hat{\omega}_{j}=4$ are larger than
the corresponding size when $\hat{\omega}_{j}=3$, respectively. The term
"turning point" denotes that as $\hat{\omega}_{j}$ increases, the size of the
parent set initially decreases and then starts increasing once the local
minimum is reached. The corresponding relations exist because as long as
$\hat{\omega}_{j}$ is not a multiple of the true $\omega_{j}$, the estimated
time partition subsets with $\hat{\omega_{j}}$ must be a mixture of some
correct time partition subsets with $\omega_{j}$. Therefore, it is reasonable
to use this trick rather than looking at the maximum size of the parent sets
$\widehat{\mathrm{Pa}}_{k}(X^{j}_{t}),k\in[\hat{\omega}_{j}]$ (line 17 in
Algorithm B1).

Given infinite samples, our estimate $\hat{\omega}_{j}$ (line 17 in Algorithm
B1) is the exact value $\omega_{j}$ (see Lemma D.5). However, for finite
samples, estimating $\omega_{j}$ by the equation in line 17 in Algorithm B1
does not yield good performance when $T$ is small. While searching, larger
guesses $\omega$ lead to finer time partitions in $\Pi_{j}$, resulting in
smaller sizes for $\Pi^{k}_{j}$ (see Line 5 in Algorithm B1). Due to the power
limit of CI tests on a smaller sample given by $\Pi^{k}_{j}$, the number of
false negative edges increases. In order to solve this issue, we introduce
turning points. A turning point is a guess $\hat{\omega}$ satisfying:

$\displaystyle\max_{t}|\widehat{\mathrm{Pa}}_{\hat{\omega}}(X^{j}_{t})|<\min\\{\max_{t}|\widehat{\mathrm{Pa}}_{\hat{\omega}-1}(X^{j}_{t})|,\max_{t}|\widehat{\mathrm{Pa}}_{\hat{\omega}+1}(X^{j}_{t})|\\}$

where $|\widehat{\mathrm{Pa}}_{\hat{\omega}}(X^{j}_{t})|$ is the estimated
parent set for $X^{j}_{t}$ with periodicity guess $\hat{\omega}$. See line 19
in Algorithm B1.

We illustrate it with a special example in Fig.6. If there are several turning
points, then $\hat{\omega}_{j}$ is the first turning point. If there is no
turning point, then we obtain $\hat{\omega}_{j}$ using Line 17 of Algorithm
B1.

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 to a multiple of the true periodicity when $T$ 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. In these experiments, the accuracy rate is calculated by
considering
$\\{N\omega_{j}\\}_{N\in\lfloor\frac{\omega_{\text{ub}}}{\omega_{j}}\rfloor}$
as correct estimations.

##  Appendix F Computational Complexity

Executing the PCMCI algorithm on the entire time series constitutes the
initial phase of the proposed approach (Algorithm B1 line 2). The algorithm’s
worst-case overall computational complexity is
$O(n^{3}\tau_{\text{ub}}^{2})+O(n^{2}\tau_{\text{ub}})$, discussed in Runge et
al. (2019). Here, the symbol $n$ denotes $n$-variate time series and
$\tau_{\text{ub}}$ represents the upper boundary for time lags.

The subsequent computational load stemming from the remaining components of
our algorithm follows a complexity of
$O(\omega^{2}_{\text{ub}}n^{2}\tau_{\text{ub}})$ . This encompasses the
$O(n^{2}\tau_{\text{ub}})$ complexity associated with conducting Momentary
Conditional Independence (MCI) tests on all $n$ univariate time series. The
parameter $\omega^{2}_{\text{ub}}$ here arises due to the search procedure
involving $\omega$, iterating through values from 1 to $\omega_{\text{ub}}$
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).

##  Appendix G Experiments

All experiments, including those detailed in the main paper, are conducted on
a 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 ($F_{1}$ 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,\tau_{\max}+1]$. Within the three-dimensional binary
array, the value 1 signifies an edge pointing from one variable to another
with a specific time lag, while 0 indicates the absence of an edge. For
instance, if element $[i,j,k]$ in the matrix is 1, then there is an edge
pointing from $X^{i}_{t-k}$ to $X^{j}_{t}$. In semi-stationary time series,
due to the presence of multiple causal mechanisms, the binary edge matrix is a
four-dimensional array with dimensions $[n,\Omega,n,\tau_{\max}+1]$, where
$\Omega$ is defined as Eq.(7) in the main paper. This expanded binary matrix
is constructed based on the edge matrix of each variable $X^{j}_{t},j\in[n]$,
through repetition. For instance, if $\Omega=2\omega{j}$, setting the third
dimension of the large binary matrix to $j$ should yield $\omega_{j}$
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,\Omega,n,\tau_{\max}+1]$ and one obtained from the algorithm
with dimensions $[n,\widehat{\Omega},n,\tau_{\max}+1]$. If the estimator
$\widehat{\Omega}$ is incorrect, those two binary arrays will have different
sizes, so we can not directly compare them. To solve this problem, we do the
same operation and calculate the least common multiple of $\Omega$ and
$\widehat{\Omega}$. Denoting this least common multiple as
LCM($\Omega$,$\widehat{\Omega}$), we create two four-dimensional binary arrays
with dimensions of $[n,\text{LCM}(\Omega,\widehat{\Omega}),n,\tau_{\max}+1]$
based on the true edge array and the 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 a 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, if the
finding is not correct from an expert’s viewpoint, the following assumptions
may be violated:

  * •

Assumption A4 No Contemporaneous Causal Effects: There is a possibility of
potential causal effects from $X^{\text{ta}}_{t}$ to $X^{\text{cp}}_{t}$ that
the algorithm is unable to capture.

  * •

Assumption A6 Hard Mechanism Change combined with limited power of CI tests:
If there is a soft mechanism change in the variables, the reliability of the
CI test of two 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 $\hat{\omega}$.

We provide a sound and robust algorithm for experts in various fields who are
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 is a temporal extension of LiNGAM and LiNGAM is an
algorithm designed for non-Gaussian data, following the work in Pamfil et al.
(2020), we also construct experiments on continuous-valued time series data
with Exponential noise. Shown as Fig.7(a), the performance of
$\text{PCMCI}_{\Omega}$, 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 on Binary Time Series

Similar to the process of generating continuous-valued time series, the
generation of binary time series also involves three steps. However, the main
difference lies in the last two steps. In 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.

(a) Exponential noise

(b) Binary time series

Figure 7: a)$F_{1}$ Score, Adjacency Precision, and Adjacency Recall when
$\omega_{\max}$ varies for experiments on continuous-valued time series with
Exponential noise, length $T=\\{500,2000,8000\\}$, $\tau_{\max}=5$ and $n=5$.
b) $F_{1}$ Score, Adjacency Precision, and Adjacency Recall when
$\omega_{\max}$ varies for experiments on binary time series with length
$T=\\{4000,8000,12000\\}$, $\tau_{\max}=3$ and $n=3$.

For discrete-valued time series, a longer time length is required. To evaluate
performance, we conduct a series of experiments following the same methodology
as described in section 4.1. Fig.7(b) illustrates the variation in
comprehensive performance with respect to $\omega_{\max}$. PCMCIΩ demonstrates
a similar performance to PCMCI in terms of the $F{1}$ score, indicating a
well-balanced trade-off between precision and recall. This outcome is expected
since discrete-valued time series demand larger sample sizes, and the
increases in $\omega_{\max}$ negatively impact the power of MCI tests. This
observation is further supported by Fig.8(a), where an increase in time length
$T$ from 4000 to 12000 does not lead to a significant improvement in the
accuracy rate of $\hat{\omega}$, while the accuracy decreases rapidly with
higher values of $\omega_{\max}$.

Comparing these results to the experiments conducted on continuous-valued time
series, it becomes evident that the demand for efficient samples is even more
substantial for binary time series, and the influence of increasing
$\omega_{\max}$ on performance becomes more pronounced.

Fig.8(b) shows how the performance of the algorithm varies across
$\tau_{\max}$ and the same trade-off between recall and precision has been
shown.

(a) Accuracy of $\hat{\omega}$

(b) Performance when $\tau_{\max}$ varies

Figure 8: $\text{PCMCI}_{\Omega}$ is tested on 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 $\hat{\omega}$ for different time
series lengths and different $\omega_{\max}$. b) $F_{1}$ Score, Adjacency
Precision, and Adjacency Recall when $\tau_{\max}$ varies for experiments with
time series length $T=4000$, $\omega_{\max}=3$ and $n=3$.

###  G.4 More experiments on Continuous-valued time series

In this section, we conduct more experiments with continuous-valued time
series with Gaussian noises.

In Fig.9(a), we test our algorithm with and without utilizing the turning
point rule. See lines 13-14 in Algorithm B1 and section E for more information
about the turning point rule. Let PCMCIΩ TP denote the version of PCMCIΩ that
the turning point rule is utilized in choosing $\omega$. PCMCIΩ non-TP means
that the turning point rule is not applied and $\omega$ is chosen directly
according to Lemma D.5.

Fig.9(a) shows that the algorithm PCMCIΩ non-TP and PCMCIΩ TP have similar
performance with various $T$ and $\omega_{\max}$. With $T=500$, PCMCIΩ non-TP
yields slightly larger standard errors for those metrics, compared to PCMCIΩ
TP. As time length $T$ increases, the performance of the algorithm PCMCIΩ non-
TP has consistently increased and is even slightly better than PCMCIΩ TP. The
consistent performance of PCMCIΩ under different chosen rules of $\omega$
supports our theoretical result; that is, the correct periodicity leads to the
most sparse causal graph.

In Fig.9(b), non-stationary time series are produced instead of semi-
stationary ones. Consequently, the causal mechanisms for each univariate time
series no longer appear sequentially and periodically. The proposed method
performs slightly better in terms of F1 score and precision. However, the
recall rate is the worst compared to other baselines.

In Fig.9(c), we conduct experiments in the nonlinear setting. The proposed
algorithms PCMCIΩ TP and PCMCIΩ non-TP perform the best.

In Fig.9(d), with $\omega_{\text{ub}}<\omega_{\max}$, the performance of the
proposed algorithm is significantly worse compared to the scenario where
$\omega_{\text{ub}}>\omega_{\max}$. However, with
$\omega_{\text{ub}}<\omega_{\max}$, the proposed algorithm can still detect a
less dense graph in comparison to other baselines. Based on these outcomes, it
is essential to maintain a slightly higher $\omega_{\text{ub}}$ without
significantly impacting the number of efficient samples utilized in each CI
test.

(a) Performance with and without turning point

(b) Performance in non-stationary setting without periodicity

(c) Performance in non-linear SCM

(d) Performance when $\omega_{\text{ub}}>\omega_{max}$ and
$\omega_{\text{ub}}<\omega_{max}$

Figure 9: Multiple algorithms are tested on 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 $\omega$
supports our theoretical result; that is, the correct periodicity $\omega$
leads to the most sparse causal graph. In (b), data sets are in a non-
stationary setting without periodicity. In (c), the structural causal model
(SCM) is non-linear. In (d), algorithm PCMCIΩ are tested under conditions that
$\omega_{\text{ub}}>\omega_{\max}$ and $\omega_{\text{ub}}<\omega_{\max}$
respectively.