Inferring causal relations from observational time series data is a key problem across science and engineering whenever experimental interventions are infeasible or unethical. Increasing data availability over the past few decades has spurred the development of a plethora of causal discovery methods, each addressing particular challenges of this difficult task. In this paper, we focus on an important challenge that is at the core of time series causal discovery: regime-dependent causal relations. Often dynamical systems feature transitions depending on some, often persistent, unobserved background regime, and different regimes may exhibit different causal relations. Here, we assume a persistent and discrete regime variable leading to a finite number of regimes within which we may assume stationary causal relations. To detect regime-dependent causal relations, we combine the conditional independence-based PCMCI method [based on a condition-selection step (PC) followed by the momentary conditional independence (MCI) test] with a regime learning optimization approach. PCMCI allows for causal discovery from high-dimensional and highly correlated time series. Our method, Regime-PCMCI, is evaluated on a number of numerical experiments demonstrating that it can distinguish regimes with different causal directions, time lags, and sign of causal links, as well as changes in the variables’ autocorrelation. Furthermore, Regime-PCMCI is employed to observations of El Niño Southern Oscillation and Indian rainfall, demonstrating skill also in real-world datasets.

Regime-dependent non-stationarity is a ubiquitous feature of physical systems, especially prominent in atmospheric sciences. This dependence can be looked at as an intermittent change in relationships defining the dynamics of a multivariate system, each of which can be described as a time series causal network. In this work, we develop a novel algorithm to detect regime-dependent causal relations that combines the constrained-based causal discovery algorithm PCMCI with a regime assigning linear optimization algorithm. Our method, Regime-PCMCI, is evaluated on a number of numerical experiments and demonstrates high performance in detecting a variety of regime-dependent features. Finally, Regime-PCMCI is applied to observations of El Niño Southern Oscillation and Indian rainfall, demonstrating skill in detecting well-known seasonal regimes in a real-world dataset.

## I. INTRODUCTION

Understanding causal relationships^{1,2} among different processes is a ubiquitous task in many scientific disciplines as well as engineering (e.g., in the context of climate research,^{3–8} econometrics,^{9,10} molecular,^{11} and animal group^{12} dynamics). Yet, the common approach to gaining causal knowledge by conducting experiments is often infeasible or unethical, for example, in Earth sciences. All that is often given is a set of time series describing these processes with no specific knowledge about the direction and form of their causal relationships available. The challenge, termed causal discovery, is then to reconstruct the underlying graph of causal relationships from time series data.^{8} Based on that graph, the processes that generated the data can be modeled in the framework of structural causal models (SCMs)^{1} to further understand causal relations, predict the effect of interventions, and for forecasting.

Today’s ever-growing abundance of time series datasets promises many application scenarios for the now numerous data-driven causal discovery methods. But many challenges emerging from the dynamic nature of such datasets have not yet been met. Furthermore, causal knowledge cannot be gained from data alone and each method comes with its particular set of assumptions^{2} about properties of the underlying processes and the observed data. Runge *et al.*^{8} recently provided an overview of current methodological frameworks, their application scenarios, and open challenges.

One such issue is posed by time-varying causal relationships, a frequent feature of both natural and artificial systems. In Earth sciences, for example, the dominant causal relationship between soil moisture and air temperature periodically reverses due to land-atmosphere feedbacks;^{13} in animal group dynamics, the leader–follower role of an individual often mutates in time;^{14} and in econometrics, the direction of influence between stock markets and macroeconomic variables is often dynamical.^{15}

A feature commonly observed in non-stationary dynamical systems is regime dependence. Regime dependence means that the causal relationships between the considered processes vary depending on some prevailing background regime that may be modeled as switching between different states. Furthermore, often such regimes have strong persistence, that is, they operate and affect causal relations on much longer timescales than the causal relations among the individual processes. In the climate system, for instance, several cases of such regime dependencies exist. For example, rainfall in India in summer is known to be influenced by the so-called El Niño Southern Oscillation (ENSO), an important mode of variability in the tropical Pacific affecting the large-scale atmospheric circulation and thereby weather patterns around the globe.^{16,17} It is, however, generally assumed that ENSO does only marginally affect Indian rainfall in winter.^{18} Thus, the causal relationships between ENSO and rainfall over India change dependent on the season that here defines the background regime and operates on a longer timescale (several months) than the causal relations among ENSO and Indian rainfall (several weeks).

### A. Existing work

Causal discovery has seen a steep rise with a plethora of novel approaches and methods in recent years. Each approach has different underlying assumptions and targets different real-world challenges as discussed in Runge *et al.*^{8} In general, causal (graph) discovery methods can be classified into classical Granger causality approaches,^{9,19} constraint-based causal network learning algorithms,^{2} score-based Bayesian network learning methods,^{20,21} structural causal models,^{22,23} and state-space reconstruction methods.^{24,25}

Here, we focus on the constraint-based framework, which has the advantage that it can flexibly account for linear and nonlinear causal relations and different data types (continuous and categorical, univariate and multivariate). In particular, we adopt the PCMCI algorithm^{26} [based on a condition-selection step (PC) followed by the momentary conditional independence (MCI) test] to reconstruct time series causal graphs. PCMCI is an adaptation of the constraint-based PC algorithm (named after its inventors Peter Spirtes and Clark Glymour^{2}) that addresses autocorrelation of time series via the use of a momentary conditional independence (MCI) test. PCMCI yields high detection power also in high-dimensional and strongly autocorrelated time series settings (see Sec. III A and Runge *et al.*^{26} for more details). However, one of the general assumptions of PCMCI (as well as of other causal discovery algorithms) is stationarity, i.e., that the existence or absence of a causal link does not change over the considered time series segment.^{27} While known changes in the background signal can be accounted for by restricting the time series to the stationary regimes, PCMCI cannot handle unknown background regimes that constitute a particular case of latent confounding.

Some recent work addresses causal discovery in the presence of non-stationarity. Malinsky and Spirtes^{28} model non-stationarity in the form of (continuous) stochastic trends in a linear autoregressive framework. Zhang *et al.*^{29} account for non-stationarity in the more general constraint-based framework. However, both address the case of a (smoothly) varying continuous background variable that continuously changes causal relations among the observed variables. This means that these methods will not output regime-dependent causal graphs, but a “summary” graph that accounts for regimes modeled as latent drivers. Peters, Bühlmann, and Meinshausen^{30} and Christiansen and Peters^{31} assumed that known non-stationary regimes are exploited to estimate causal relations also in the presence of general latent confounders. Furthermore, in the context of continuously varying causality, methods based on information transfer metrics have been proposed. In the field of animal group dynamics, for instance, detection of time-varying leader–follower relationships is achieved with the use of a time dependent transfer entropy.^{23,32,33} Applied to non-stationary climate systems, Hagan *et al.*^{34} proposed a Kalman filter estimate of the time-varying parameters for the Liang–Kleeman information flow. Benefits of these methods are the treatment of non-linearity^{32} and the identification of both timing and frequency of interactions.^{34} However, in these approaches, only bivariate influences are modeled, i.e., the effect of a third variable $Z$ on the estimated effect of $X$ onto $Y$ cannot be accounted for. The practical extension to high-dimensional systems and short time series also remains hard to address.

Currently, few methods exist that address the case of a discrete regime variable leading to distinct causal regimes that may be physically interpreted. For example, in the climate science context, regime-dependent autoregressive models (RAMs) were introduced already in 1990.^{35} These can yield physically well interpretable results that, however, require well-chosen ancillary variables and a seasonal index that are not learned from data. Thus, RAM requires *a priori* knowledge of the regimes, which one often aims to learn rather than enforce. In the context of discrete state spaces, regime-dependent causal discovery has been considered in Gerber and Horenko^{11} for Boolean variables. Non-stationary Boolean network models have also been considered in Porfiri and Marin^{23}—specifically, the approach is to fit an appropriate parameterization of associate transition probabilities. Another approach that has been proposed to model time dependent Granger (non-) causality is based on a Markov Switching VAR ansatz with an economics application in mind.^{10} Specifically, the regime assignments are computed by sampling from a Markov chain. Further methods have been proposed to obtain time step specific bivariate Granger Causality from partitioning the time-series into regular time windows.^{36,37}

A more general framework to handle discrete regimes is the Markov-switching ansatz of de Wiljes *et al.*,^{38} which flexibly models regime dependence utilizing the assumption of a finite number of regimes and a level of persistence in the transitions between different regimes. This ansatz has been successfully realized in combination with many different model assumptions (e.g., see Refs. 39 and 40). Here, we want to explore it for causal networks by combining it with PCMCI,^{26} a constraint-based time series causal discovery method.^{2} We call our method Regime-PCMCI.

The remainder of the paper is structured as follows: In Sec. II, the underlying mathematical problem, concepts, and key assumptions are formalized, and a motivating example is discussed to provide some intuition. Our novel method Regime-PCMCI is then presented in Sec. III. These theoretical and algorithmic parts are complemented by a thorough numerical investigation of the proposed method in various artificial settings in Sec. IV. Finally, in Sec. V, Regime-PCMCI is applied to a real-world dataset from climate science, addressing the changing relationships of ENSO and rainfall over India.

## II. PROBLEM SETTING

Let ${Xt}t\u2208Z$ be a sequence of real-valued $NX$ dimensional random variables $Xt\u2208RNX$, where $t$ is associated with time. A realization over the time interval $[0,T]$ of this stochastic process is denoted ${xt}t\u2208[0,T]$, and we assume that it is possible to obtain observations of these realizations. We assume that the underlying process is modeled by a regime-stationary discrete-time structural causal model (SCM),

Here, the measurable functions $gtj$ depend non-trivially on all their arguments, the noise variables $\eta tj$ are jointly independent and are assumed to be stationary, i.e., $\eta tj\u223cDj$ for all $t$ for some distribution $Dj$, and the sets $Ptj\u2282(Xt\u22121,Xt\u22122,\u2026)$ define the causal parents of $Xtj$. Note that we assume lagged causal relationships, but this is not a necessity since there exist causal discovery algorithms that can deal with contemporaneous causal links^{41} and also hidden confounders. In contrast to approaches assuming stationarity, both $gtj$ and $Ptj$ are allowed to depend on regimes in time as further formalized in Assumption 1 (Sec. II B).

List of notations . | Model . | Parameters . | . |
---|---|---|---|

{X_{t}}_{t∈[0,T]} | Stochastic process | N_{X} | Dimension of stochastic process |

$\eta tj$ | Noise variable of component $Xtj$ | T | Time length of stochastic process |

$Dj$ | Stationary noise distribution | N_{K} | Number of regimes |

$gtj$ | Structural causal model function | N_{C} | Max switches for each regime |

$Ptj$ | Causal parents of X^{j}, time dependent | N_{M} | Regime average persistence |

x_{t} | Realization of X_{t} | τ_{max} | Maximum causal time lag |

$G^t$ | Operator in inverse problem | CI test | Conditional independence test |

$g^tj$ | Components of operator $G^t$ | α_{PC} | Significance level for PC_{1} step |

Θ _{t} | Unknown parameter in inverse problem | α | Link significance level |

$L(\Gamma ,P,\Phi )$ | Cost functional | N_{Q} | Number of optimization iterations |

Γ(t) | Regime-assigning process | N_{A} | Number of annealings |

Φ_{t} | Linear link coefficients, time dependent | N_{R} | Number of realizations for a toy example |

γ_{k}(t) | Regime-assigning process for regime k | $N(0,{\sigma 2}ref)$ | Ground-truth Gaussian noise distribution |

$\Phi kj(i,\tau )$ | Linear link coefficient in regime k | ${\Phi kj(i,\tau )}ref$ | Ground-truth linear link coefficient |

$Pkj$ | Causal parents of X^{j} in regime k | N_{para} | Number of model parameters |

$\Upsilon k$ | Collection of time steps associated with regime k | $x^k,t$ | Reconstructed time-series for regime k |

List of notations . | Model . | Parameters . | . |
---|---|---|---|

{X_{t}}_{t∈[0,T]} | Stochastic process | N_{X} | Dimension of stochastic process |

$\eta tj$ | Noise variable of component $Xtj$ | T | Time length of stochastic process |

$Dj$ | Stationary noise distribution | N_{K} | Number of regimes |

$gtj$ | Structural causal model function | N_{C} | Max switches for each regime |

$Ptj$ | Causal parents of X^{j}, time dependent | N_{M} | Regime average persistence |

x_{t} | Realization of X_{t} | τ_{max} | Maximum causal time lag |

$G^t$ | Operator in inverse problem | CI test | Conditional independence test |

$g^tj$ | Components of operator $G^t$ | α_{PC} | Significance level for PC_{1} step |

Θ _{t} | Unknown parameter in inverse problem | α | Link significance level |

$L(\Gamma ,P,\Phi )$ | Cost functional | N_{Q} | Number of optimization iterations |

Γ(t) | Regime-assigning process | N_{A} | Number of annealings |

Φ_{t} | Linear link coefficients, time dependent | N_{R} | Number of realizations for a toy example |

γ_{k}(t) | Regime-assigning process for regime k | $N(0,{\sigma 2}ref)$ | Ground-truth Gaussian noise distribution |

$\Phi kj(i,\tau )$ | Linear link coefficient in regime k | ${\Phi kj(i,\tau )}ref$ | Ground-truth linear link coefficient |

$Pkj$ | Causal parents of X^{j} in regime k | N_{para} | Number of model parameters |

$\Upsilon k$ | Collection of time steps associated with regime k | $x^k,t$ | Reconstructed time-series for regime k |

The problem setting considered in this manuscript is of the nature of the following inverse problem:

with $G^t=[g^t1,\u2026,g^tNX]$, where $g^tj$ belongs to an appropriate function space for each $t$ and $j$. $\tau max$ is the maximum considered time lag. In other words, the aim is to fit a set of unknown parameters $\Theta t$ on the basis of an observed time series ${xt}t\u2208[0,T]$. In Sec. II A, we will discuss the particular structure of the parameters $\Theta t$ we are interested in. Please refer to Table I for a summary of the notation used throughout the text.

### A. Causal graphs

Representing causal relations between different processes as graphs (also referred to as networks) is common practice in the context of causal inference and causal discovery.^{1,2} For time series, we use the concept of time series graphs. The nodes in the time series graph associated with the SCM (1) are the individual time-dependent variables $Xtj$ with $j=1,\u2026,NX$ at each time $t\u2208Z$. Variables $Xt\u2212\tau i$ and $Xtj$ for a time lag $\tau >0$ and a given $t$ are connected by a lag-specific directed link, denoted $Xt\u2212\tau i\u2192Xtj$, when $Xt\u2212\tau i\u2208Ptj$ for a particular $t$. If a SCM is not given, another way to define links is that $Xt\u2212\tau i$ is not conditionally independent of $Xtj$ given the past of all variables, defined by , with denoting the absence of a (conditional) independence.^{27} Widely used in constraint-based methods, note that this definition cannot be used to define directed contemporaneous links. We denote the maximum ground-truth time lag of any parent as $\tau maxP$. For a more detailed introduction, the reader is referred to Runge *et al.*^{26} In the following, we will use graphs and networks interchangeably.

The collection of parent sets for all components at time $t$ is denoted $Pt={Pt1,\u2026,PtNX}$. This set of parents is part of the unknown parameters we want to infer. Note that their dimensionality is assumed finite, but not known *a priori*. The other quantity of interest is the functional form of the causal relations $gtj(Ptj,\eta tj)$ in SCM (1) corresponding to these links. We will assume a known function class $G^t(\u2026;\Phi t)$ of unknown coefficients$\Phi t={\Phi t1,\u2026,\Phi tNX}$ that are going to be inferred via

In other words, for a given time series, $xt\u2208RNX$ with $t\u2208[0,T]$ and known function class $G^t$ the aim is to find the unknown parameters $\Theta t=[Pt,\Phi t]$.

Note that we will specifically focus on linear function classes, as discussed in Sec. III.

### B. Persistence

As mentioned above, in many application areas, non-stationarity may be modeled not in the form of abrupt or continuous changes, but via piece-wise stationary regimes.^{11,42,43} These regimes will further exhibit a certain persistent behavior. In order to capture non-stationary systems with these properties, we will restrict our inference to regime-dependent persistent dynamics.

Denote the causal parents and functional dependency of a given variable $j$ for a regime $k$ as $Ptj=Pkj$ and $gtj(Ptj,\eta tj)=gkj(Pkj,\eta tj)$. We call a regime $(NM,NK)$-persistent if the parents and functional dependencies are stationary for an average of $NM$ consecutive time steps $t$ and that there is a finite number of regimes on the whole time domain, i.e., $k\u2208{1,\u2026,NK}$.

The persistency enters here via the regime average persistence $NM$, which also naturally implies a finite number of regimes $NK\u2264T/NM$.

Under Assumption 1, the considered inverse problem (3) reduces to finding the unknown parameters $\Theta t=[\Gamma (t),P,\Phi ]$ comprising (1) a set of regimes’ network parameters

and, to encode their time dependence, (2) the change points between the regimes given by the regime-assigning process

with $\Gamma (t)\u2208[0,1]NK\xd7T$. For example, component $k$ of the regime-assigning process can be of the form $\gamma k=(0,1,1,\u20260,1)\u2208[0,1]T$, indicating that regime $k$ is active for all time steps for which $\gamma k(t)=1$.

### C. Optimization problem

Finally, in order to solve the inverse problem (3) under the persistency Assumption 1, we can define a cost functional

subject to constraints

and

where $d$ is a distance measure such as the squared Euclidean distance $\u2225\u22c5\u222522$ and $\gamma k(t)$ can be regarded as the weight of the $k$ regime-specific network at each time $t$.

This learning approach is based on the ideas first proposed in Horenko^{39} and later extended to many different models.^{38} The format of $L(\Gamma ,P,\Phi )$ relies on the assumption that the system associated with the considered data exhibits metastability in time (see Assumption 1 that translates in the summation over $k$, controlled by the regime number $NK$). The desired level of persistence enters the functional in the form of a regularization [see Constraint (6), controlled by parameter $NC$]. An alternative option is to add a regularization term that enforces some form of smoothness of $\Gamma $ (e.g., Tikhonov regularization^{44}).

The tuning parameter $NC$ is related to the *average* regime duration of $NM$ time steps of Assumption 1 as follows: an average regime duration of $NM$ in all $NK$ regimes is implemented by choosing $NC\u2248T/(NMNK)$. Importantly, note that the regularization (6) ensures an average persistence on each regime without imposing that *individual* regime durations are a constant; in fact, they can be fully irregular within the bound of performing at maximum $NC$ switches. This regime learning method thus provides a simple, flexible, and computationally tractable strategy to go beyond the assumption of fixed length for each regime duration often employed in previous methods (e.g., Refs. 37 and 33).

Note that, in practice, it is reasonable to assume that an estimate of the average regime switching time $NM$ is available, consistent with a typical timescale of the application domain. The choice of parameters $NK$ and $NC$($NM$) will be discussed in Sec. III D.

### D. Motivating example

Before we introduce our novel regime detecting the causal discovery algorithm, we illustrate the underlying challenges of causal discovery in the face of regime dependence by giving a simple example. Consider the case of two background regimes and two time series $X1$ and $X2$ and the associated causal graphs as shown in Fig. 1(a). Variable $X1$ linearly influences $X2$ but the sign changes in time, alternating between a positive (during regime 1) and a negative (during regime 2) influence. Here, the two regimes alternate equidistantly in time. The cross correlation of $X1$ and $X2$ over the whole time period is zero because the opposite sign effects cancel each other out in the linear regression. Thus, any linear causal discovery method would fail in detecting the influence of $X1$ on $X2$ when no *a priori* knowledge on the two background regimes exists. For example, applying a linear version of PCMCI on the whole time sample would give a network of disconnected variables [Fig. 1(b)].

In contrast, if the regimes are known and PCMCI is applied to samples from both regimes separately, the positive and negative links are correctly detected (not shown). To deal with such problems automatically, our algorithm needs to learn both the regimes and the regime-dependent causal relations.

## III. METHOD

The proposed approach, Regime-PCMCI, is designed to solve the optimization problem (4) by alternating between learning the regimes and learning the causal graphs for each regime in an iterative fashion. In principle, any causal discovery method that yields a causal graph can be used. Here, we chose PCMCI^{26} as a well-tested method that adapts the constraint-based causal discovery framework to the time series case.

### A. Causal discovery

The constraint-based framework has the advantage that it can flexibly account for various functional causal relations and different data types (continuous and categorical, univariate and multivariate) since it is based on discovering causal links by means of conditional independence (see link definition in Sec. II A). Two variables $X$ and $Y$ are conditionally independent given a (potentially multivariate) variable $Z$, denoted $X\u22a5\u22a5Y|Z$, if

where $p$ denotes the associated probability density functions.

To practically test this relationship, there exist a large variety of conditional independence tests; see Runge *et al.*^{26,27} for a discussion. If relationships are assumed linear, as is the case of the present work, partial correlation can be used such that it can be shown that if the partial correlation between $X$ and $Y$ conditioned on $Z$ is significantly different from $0$.

As mentioned in Sec. I A, PCMCI is based on the constraint-based PC algorithm^{2} combined with the momentary conditional independence (MCI) test. It consists of two stages. (1) At first, the PC$1$ condition pre-selection method is run to identify relevant conditions $B^tj$ for all time series variables $Xtj$. More specifically, PC$1$ is a Markov set discovery algorithm based on the PC-stable algorithm^{47} that removes irrelevant conditions for each of the $NX$ variables by iterative independence testing. (2) Then, the MCI test is performed, to confirm whether $Xt\u2212\tau i\u2192Xtj$ by means of testing

Thus, MCI conditions on both the (potential) parents of

A causal interpretation of the relationships estimated with PCMCI comes from the standard assumptions in the constraint-based framework,^{2,26,27} namely, causal sufficiency, the causal Markov condition, faithfulness, non-contemporaneous effects, and stationarity within the regimes as further discussed below. As demonstrated in Runge *et al.*,^{26} PCMCI has high detection power and controlled false positives also in high-dimensional and strongly autocorrelated time series settings.

The main free parameters of PCMCI are the chosen conditional independence test, the maximum time lag $\tau max$, and the significance levels $\alpha $ in MCI and $\alpha PC$ in PC$1$. We discuss the selection of these parameters in Sec. III D. In terms of conditional independence test, note that PCMCI can be used in combination with linear or nonlinear tests and can therefore extract also non-linear causal relationships. In this work, we focus on linear systems and thus use PCMCI in conjunction with partial correlation.

In each iterative step of our approach, PCMCI is applied to the sample subset of the time series pertaining to the estimated $k$ regime. Given a significance level $\alpha $, the output of PCMCI is the set of parents $Pk={Pk1,\u2026,PkNX}$ for all time series variables for that regime,

### B. Regime-dependent causal discovery

The Regime-PCMCI algorithm iterates over two major estimation steps: (step 1) causal discovery to obtain $Pk$ and fit the coefficients $\Phi k$ and (step 2) regime learning to update the regime variable $\Gamma $.

To find good estimates of the parameters and the regime variable optimizing the cost functional (4), this two-step procedure is necessary. In fact, there are generally no analytic solutions to the problem available due to the complexity of the cost functional. Fixing one variable to estimate the other allows in both cases to solve the individual optimization step via linear programming. In Theorem 2.1 of Ref. 39 and Sec. II B in Ref. 48 it is shown that these types of algorithms monotonously decrease the value of $L$. It is important to note, however, that due to the non-convexity of the underlying problem the algorithm can be caught in regions of local minima. This issue is addressed via additional simulated annealing steps as discussed in more detail in Sec. III B 2.

In the following, $q$ indicates the current iteration. The superscript $(q)$ is added combined with brackets to the variables updated in each loop. The details of the consecutive subroutines are laid out below.

#### 1. Step 1: Causal discovery and model estimation

The first step is to estimate a set of parents ${Pk}(q)$ and coefficients ${\Phi k}(q)$ with $k\u2208{1,\u2026,NK}$ on the basis of a fixed ${\Gamma (t)}(q)$ obtained in step 2 of the previous iteration (first and second bullets in *q*-loop of of Algorithm 1). In the first iteration, the regimes are assigned randomly. ${Pk}(q)$ and ${\Phi k}(q)$ are estimated on the basis of a subset of the time series $xt$ with

for each regime $k$. The regime-dependent parents ${Pk}(q)$ are estimated via PCMCI.

As stated at the beginning of Sec. III, to solve Eq. (3), we assume a linear functional relationship that relates each variable to its parents $Pk$. It implies that coefficients $\Phi k$ can be estimated from the following regression model for each fixed $k$:

for $t\u2208{\Upsilon k}(q)$. In other words for every $k\u2208{1,\u2026,NK}$, the following optimization has to be solved:

for $t\u2208{\Upsilon k}(q)$. Note that the coefficients not indicated as relevant via the parent set are defined to be zero, i.e., $\Phi kj(i,\tau ):=0$ for $Xt\u2212\tau i\u2209Pkj,(q)$.

#### 2. Step 2: Regime learning

Step 2 is to determine an optimal regime-assigning process ${\Gamma t}(q+1)\u2208[0,1]NK\xd7T$ given the current estimates ${Pk}(q)$ for the parents and ${\Phi k}(q)$ coefficients (see third bullet in $q$-loop of Algorithm 1). In agreement with the cost functional (4), the following optimization problem needs to be solved: find

Since the first $\tau max$ time steps cannot be predicted, we choose to set those to $x^k,tj=xk,tj$ and to not consider this portion of the time series in the algorithm evaluation. This step can be solved with standard optimization linear programming routines.

In order to search for the global minimum of this non-convex problem, the algorithm is run for a number $NA$ of different initializations of ${\Gamma}(0)$ (annealing). The annealing run with the lowest cost functional objective is chosen as an optimal fit. Note that the individual annealing steps are *embarrassingly parallelizable*.

### C. Reconstruction of time series

A single prediction from Eq. (14) can be derived as the weighted sum over $k$

### D. Parameter selection

Regime-PCMCI involves a number of parameters that need to be chosen. They can be separated into parameters of the causal discovery method PCMCI and those of the regime learning part.

The main free parameters of PCMCI are the chosen conditional independence test, the maximum time lag $\tau max$, and the significance levels $\alpha $ in MCI and $\alpha PC$ in PC$1$. $\alpha PC$ should be regarded as a hyper-parameter and can be chosen based on model-selection criteria such as the Akaike information criterion (AIC)^{49} or cross-validation. $\tau max$ could be incorporated into this model selection. But since PCMCI is not very sensitive to this parameter^{26} (as opposed to, e.g., Granger causality), its choice can be based on lagged correlation functions, see Runge *et al.*^{26} for a discussion. The choice of conditional independence test is a modeling assumption guided by the assumed nonlinearity of the underlying process and also finite sample considerations. Finally, $\alpha $ is chosen based on the desired level of false positives.

The two free parameters in the regime learning step are the bound on the number of switches $NC$ and the number of regimes $NK$. Usually, $NC$ can be reasonably inferred from the application and given the number of regimes, as explained in Sec. II C. Here, we assume that $NC$ is known. Note that since $NC$ bounds the maximum number of switches between regimes, i.e., the optimal reconstructed number can be lower, and does not constrain the individual regime durations, its choice can account for a degree of error. Yet determining a suitable choice of the unknown number of regimes $NK$ is a difficult task. In particular, it is hard to find the right balance between avoiding to overfit and to choose appropriately complex models to describe a specific dataset and thus the underlying dynamics well. One way to assess this balance heuristically is to employ an information criterion (IC),^{50} which has been derived in the context of regression models and since been adapted to various other model scenarios including graphs.^{51}

An IC is designed to capture the goodness of fit penalized by the number of parameters in order to prefer models with as few parameters as possibles to avoid overfitting (parsimony). Here, the number of parameters is defined as

The first term in Eq. (16) relates to the number of parameters required to describe $\Gamma $, which can be fully determined via the change points. The second term in Eq. (16) counts the number of relevant parents, or equivalently the non-zero coefficients $\Phi kj(i,\tau )$. Here, we use the corrected Akaike information criterion (AICc) first proposed in Hurvich and Tsai^{52} to estimate $NK$, assuming known $NC$. Note that we use the corrected version of the original AIC^{49} to correct for small samples sizes relative to the number of parameters

where $L$ is the maximum value of the likelihood function for the model one assumes for the residuals (see Metzner *et al.*^{48} for a more detailed discussion). Note that the AICc also depends on $NC$ (as it enters the number of parameters $Npara$) and it is in general possible to simultaneously estimate $NK$ and $NC$.^{38,40} The choice of $NK$ is numerically investigated in Sec. IV C.

The number of iteration steps $NQ$ should be chosen to ensure that the optimization process converges. In our experiments, we found with exploratory testings that $NQ$ shows convergence after about 10–20 iterations for all examples investigated. The number of annealing steps $NA$ should be chosen to ensure spanning a large number of local solutions to this non-convex optimization problem [Eq. (4)]. Computational time will set a limit to a too high parameter. Note, however, that the annealing part is embarrassingly parallelizable.

## IV. NUMERICAL INVESTIGATION

In the following, we investigate the performance of Regime-PCMCI by means of several toy examples. The artificial data are designed to test the methods robustness and accuracy with respect to various potential scenarios that could occur in real applications. At first, low dimensional ($NX=2$) causal relations are studied as the results can be interpreted more easily. Next, we also consider higher dimensional settings ($NX=10$). The reference time series are generated with the following linear SCM time series model:

with predefined ${\Gamma (t)}ref$, ${\Phi k}ref$, and ${\sigma 2}ref$. Note that here we numerically investigate equally distributed noises for all variables ($\eta tj\u223cD\u2200j$), but we refer the interested reader to Appendix A for a treatment of heterogeneous noise distributions. Note that the reference set of parents is specified by the non-zero coefficients ${\Phi kj(i,\tau )}ref$.

### A. Low dimensional data with two underlying regimes

First, we focus on a simple setting of two regimes, i.e., ${NK}ref=2$, and a two dimensional underlying process $Xt\u2208R2$ (i.e., $NX=2$). Our aim is to test the performance of Regime-PCMCI for different elemental features that can change between regimes. For brevity, links $Xt\u2212\tau i\u2192Xtj$ will be called auto-links or auto-dependencies for $i=j$ and cross-links for $i\u2260j$. We consider the following scenarios as summarized in Table II: sign change of coefficient (in auto-link and cross-variables link), lag change (in cross-link), coefficient change (in auto-link), and child–parent inversion defined via an assortment of linear functions and associated coefficients. In all examples, each variable is also auto-linked at lag 1 (linear coefficient $0.2$), which is a realistic yet challenging assumption for many algorithms.

Example . | k = 1
. | k = 2
. | ${\Phi 1j(i,\tau )}ref$ . | ${\Phi 2j(i,\tau )}ref$ . |
---|---|---|---|---|

Arrow direction | X^{1} → X^{2} | X^{1} ← X^{2} | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 21(2,1)}ref=0.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | |||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | |||

Causal effect | $X1\u2192|a|X1$ | $X1\u2192|b|X1$ | ${\Phi 11(1,1)}ref=0.8$ | ${\Phi 21(1,1)}ref=0.1$ |

${\Phi 12(2,1)}ref=0.4$ | ${\Phi 22(2,1)}ref=0.4$ | |||

Lag | $X1\u2192\tau =1X2$ | $X1\u2192\tau =2X2$ | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 22(1,2)}ref=0.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | |||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | |||

Sign X^{1} | $X1\u2192|a|X1$ | $X1\u2192\u2212|a|X1$ | ${\Phi 11(1,1)}ref=0.8$ | ${\Phi 21(1,1)}ref=\u22120.8$ |

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | |||

Sign X^{1}X^{2} | $X1\u2192|a|X2$ | $X1\u2192\u2212|a|X2$ | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 22(1,1)}ref=\u22120.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | |||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ |

Example . | k = 1
. | k = 2
. | ${\Phi 1j(i,\tau )}ref$ . | ${\Phi 2j(i,\tau )}ref$ . |
---|---|---|---|---|

Arrow direction | X^{1} → X^{2} | X^{1} ← X^{2} | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 21(2,1)}ref=0.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | |||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | |||

Causal effect | $X1\u2192|a|X1$ | $X1\u2192|b|X1$ | ${\Phi 11(1,1)}ref=0.8$ | ${\Phi 21(1,1)}ref=0.1$ |

${\Phi 12(2,1)}ref=0.4$ | ${\Phi 22(2,1)}ref=0.4$ | |||

Lag | $X1\u2192\tau =1X2$ | $X1\u2192\tau =2X2$ | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 22(1,2)}ref=0.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | |||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | |||

Sign X^{1} | $X1\u2192|a|X1$ | $X1\u2192\u2212|a|X1$ | ${\Phi 11(1,1)}ref=0.8$ | ${\Phi 21(1,1)}ref=\u22120.8$ |

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | |||

Sign X^{1}X^{2} | $X1\u2192|a|X2$ | $X1\u2192\u2212|a|X2$ | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 22(1,1)}ref=\u22120.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | |||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ |

#### 1. Experiment settings

We design five toy models, in network terms, corresponding to different sets of parents defined via the reference parameters ${\Phi kj(i,\tau )}ref$ given in columns $4$–$5$ of Table II. Furthermore, synthetic regime-assigning processes ${\Gamma (t)}ref$ are generated for all examples. More specifically, ${\gamma 1(t)}ref$ is designed to consist of $41$ alternating windows, i.e., ${NC}ref=40$ regime transitions. The length of these windows is randomly selected to be between $70$ and $100$ and the constraint (5) imposes ${\gamma 2(t)}ref=1\u2212{\gamma 1(t)}ref$. The final length of the time series is capped at $T=3000$ to ensure equally long regime assignment time series.

Then, an artificial time series $xt$ via (18) with ${\sigma 2}ref=1$ is generated. Note that the stochastic process (18) can be exactly reconstructed via the coefficients ${\Phi kj(i,\tau )}ref$, their activation ${\Gamma (t)}ref$, and a specific realization of the innovation term $\eta tj$.

The PCMCI parameters are chosen as follows: partial correlation as a conditional independence test, $\alpha =0.01$, $\alpha PC=0.2$ as recommended in Runge *et al.*^{53} $\tau max=3$, and masking type “y” (see the documentation of tigramite for the definition of masking types). The number of regimes was set to $NK=2$, and the maximum number of regime transitions is $NC=40$, i.e., correct guess on number of regimes and switches (model selection for $NK$ is investigated in Sec. IV C). The number of iterations is $NQ=20$, and the number of annealings is $NA=50$. A summary of the parameters is shown in Table III. We generate $NR=100$ time series realizations for each example.

#### 2. Results

The ability of the proposed method to recover the networks and the regimes on the basis of the artificially designed time series are presented in the following. Figures 2–6 present results for each case in Table II, focusing on one of the $NR$ synthetic datasets. Table IV and Fig. 7 show summary statistics over all $NR$ runs.

Example . | Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. | $\u03f5^$ . |
---|---|---|---|---|---|---|---|---|---|---|

Arrow direction | 3.0 | 1.0 | 1.0 | 0.02 | 0.01 | 0.021 | 0.020 | 7.0 | 7.0 | 0.76 |

Causal effect | 43.0 | 0.81 | 0.98 | 0.11 | 0.01 | 0.286 | 0.020 | 120.0 | 10.0 | 0.68 |

Lag | 6.0 | 0.98 | 1.0 | 0.04 | 0.01 | 0.027 | 0.018 | 11.0 | 8.0 | 0.68 |

Sign X^{1} | 4.0 | 0.98 | 1.0 | 0.03 | 0.01 | 0.033 | 0.016 | 10.0 | 6.0 | 0.65 |

Sign X^{1}X^{2} | 3.0 | 0.99 | 1.0 | 0.01 | 0.01 | 0.028 | 0.019 | 9.0 | 7.0 | 0.75 |

Example . | Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. | $\u03f5^$ . |
---|---|---|---|---|---|---|---|---|---|---|

Arrow direction | 3.0 | 1.0 | 1.0 | 0.02 | 0.01 | 0.021 | 0.020 | 7.0 | 7.0 | 0.76 |

Causal effect | 43.0 | 0.81 | 0.98 | 0.11 | 0.01 | 0.286 | 0.020 | 120.0 | 10.0 | 0.68 |

Lag | 6.0 | 0.98 | 1.0 | 0.04 | 0.01 | 0.027 | 0.018 | 11.0 | 8.0 | 0.68 |

Sign X^{1} | 4.0 | 0.98 | 1.0 | 0.03 | 0.01 | 0.033 | 0.016 | 10.0 | 6.0 | 0.65 |

Sign X^{1}X^{2} | 3.0 | 0.99 | 1.0 | 0.01 | 0.01 | 0.028 | 0.019 | 9.0 | 7.0 | 0.75 |

The case *sign $X1X2$* is discussed in detail. The ground-truth regime evolution and networks are shown in the top part of panels a and b in Fig. 2; in the middle part of both panels, their Regime-PCMCI reconstruction is shown; and in the bottom part, the difference between reconstructed and true regimes is presented to visually inspect the accuracy. The reconstructed regime-assigning process for each regime matches the truth in 98.6% of time steps (97% average value over $NR$, see Table IV). The corresponding networks have all and only the correct links ($TPR=0.99$ and $FPR=0.01$ average value over $NR$); their linear causal effect is also well estimated with each link correct up to $\xb10.02$ [$NR$-averaged error per link is $0.028$ (9%)].

The other four cases are presented in Figs. 3–6. The case *causal effect*, and to a lesser extent *lag change*, are hardest to detect. This is because the difference between the individual regimes and a mixed state of the two is not very large and thus the detection is more challenging. This adds to the general challenge of non-convexity of the functional we are optimizing, which we mitigate by the annealing steps as mentioned in Sec. III D. A similar challenge is found for some high-dimensional runs for which we refer to Sec. IV D.

The average accuracy of Regime-PCMCI is estimated from $NR$=100 synthetic datasets per each example and is presented in Fig. 7 and Table IV.

For a compact overview of the results and to facilitate the comparison between examples, Fig. 7 focuses on two key statistics: the precision of the reconstructed regime-assigning process (lightblue box plot, $\Delta \gamma %$) and the precision of the reconstructed links’ causal effects (pink box blot, $\Delta \Phi %$). $\Delta \gamma %$ is the average percentage of wrongly estimated time steps per regime (the lower the better, note that this value is the same for $k=1,2$, by construction). $\Delta \Phi %$ is the average difference between the reconstructed linear coefficient and the reference values of the ground-truth links expressed as percentage, i.e., each difference is weighted by the absolute value of the ground-truth coefficient. The precise definition of the statistics can be found in Appendix B. These quantities provide a summary of the reconstructed regimes’ accuracy, previously shown in the bottom parts of Figs. 2–6, and are presented for all $NR$ runs. The examples are ranked in the order of decreasing performance.

As already seen from inspection of single runs, change of *arrow direction* and of *causal effect* are the hardest to detect (the highest error in Fig. 7). It is also clear that a worsening performance in regime detection results in a higher network error. Except for *causal effect*, regime error $\Delta \gamma %$ is between 1% and 7% and network error $\Delta \Phi %$ is between 1% and 10%. Note that the average $\Delta \Phi %$ (dark pink cross) is very close to $\Delta \Phi ref%$ (black dot), error for PCMCI runs with the ground-truth regime variable known (but causal structures unknown). This reference value sets the optimal baseline to which compare Regime-PCMCI performance and is met by our algorithm in all but one example.

Table IV shows a more detailed summary of the results over the $NR$ realizations. The estimation errors are presented in terms of the regime-assigning process (the second column), the network structure (third to sixth columns), the causal effects of links (seventh to tenth columns) and the overall reconstructed time series (the last column). The second column is $\Delta \gamma %$, the average percentage of wrongly estimated time steps per regime, introduced above. In terms of networks, the link detection performance is evaluated via the true positive (TPR) and false positive rates (FPR). Furthermore, we compare these with the reference FPR and TPR (superscript *ref *) if PCMCI is run with the ground-truth regime variable known (but causal structure unknown). The accuracy in links’ causal effects is assessed via $\Delta \Phi $, the average difference between the reconstructed linear coefficient and the reference values of the ground truth links, and in the percentage version ($\Delta \Phi %$, see above). The last column, $\epsilon ^$, is the expected prediction error per variable and per time step and is computed as $\epsilon ^=L/(NXT)$ with $L$ defined in Eq. (4) and $NX$ and $T$ referring to the number of variables, here two, and the length of the time series, respectively. The precise definition of all the above statistics can be found in Appendix B.

In summary, Table IV shows that

$\Delta \gamma %$: on average, the regime-assigning process is reconstructed correctly in $\u223c94%$ of the time steps for all cases except

*causal effect*. The*causal effect*and*lag*examples are the hardest to infer, with the causal effect being particularly deficient. In these examples, a mixed-regime state (e.g., arising from assigning a considerable fraction of wrong time steps to a regime) is still quite close to any of the true regimes. Therefore, the algorithm struggles to decide which time steps belong to which regime since they could fit both to some degree. Yet, also for*causal effect*there are seven instances where $\Delta \gamma <15%$ (one presented in Fig. 6) and those, as expected from PCMCI, give a very good network fit. We notice that these runs do not correspond to the lowest objective values of the $NR$ set (i.e., better fit), which shows that runs that end up in mixed states can still fit the data quite well. Also, we notice that the causal effect setup reaches local minima in 16% of the 100 runs, thus in 84% of the runs the algorithm cannot easily find a stable solution which points at a weaker confidence in the output.TPR: despite some errors in reconstructing the regime-assigning process, the TPR is always very close to 1. This can indicate that the true signals, dynamicwise, are strong enough to be detectable.

FPR: Ideally, the false positive rate should be upper-bounded by $\alpha =0.01$. This is also the case if we assume the correct regimes (see column $FPRref$). However, if the regimes are learned, in most of the examples the FPR value is higher due to errors in learning the regimes. If a wrong regime is learned, then both false positives and false negatives can occur. False negatives, i.e., missing links in the PC$1$ step of PCMCI, can lead to false positives in the MCI step.

$\Delta \Phi %$: Errors in parents’ detection [either due to false positives (FPR) or to false negatives (missed links, FNR = 1$\u2212$TPR)] surely impact the estimation of link effects. Since the TPR and FPR are good, except for the causal effects case, we expect to obtain also good results for the linear coefficients. This is indeed the case, as the difference is of order $10\u22122$, implying a relative error of about $10%$. Also, this matches very closely the optimal baseline for the network fit $\Delta \Phi ref%$.

### B. Low dimensional data with three underlying regimes

To illustrate how Regime-PCMCI deals with more than two regimes, we also considered a toy time series based on three different causal regimes. It is, of course, possible to consider the case $NK>3$, yet in applications it is often desirable to infer a few prominent and relevant regimes rather than having too many that are not interpretable anymore. In other words, the aim is to avoid overfitting and to increase the information gain by reducing the complexity of the assumed model (parsimony).

The artificial time series is generated via a regime-dependent causal graph that is designed by combining two of the regimes settings presented in Sec. IV A, namely, *sign $X1X2$* change and *arrow inversion* (for details, see Table V). The regime-assigning reference process ${\Gamma}ref$ is generated by randomly choosing between different persistence lengths of $60$, $70$, and $80$ time steps and iterating for $20$ times. The algorithm is run with free parameters in Table VI.

Example . | k = 1
. | k = 2
. | k = 3
. | $\Phi 1j(i,\tau )ref$ . | ${\Phi 2j(i,\tau )}ref$ . | ${\Phi 3j(i,\tau )}ref$ . |
---|---|---|---|---|---|---|

Sign and X^{1}X^{2}arrowdirection | $X1\u2192|a|X2$ | $X1\u2192\u2212|a|X2$ | $X2\u2192|a|X1$ | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 22(1,1)}ref=\u22120.8$ | ${\Phi 31(2,1)}ref=0.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | ${\Phi 31(1,1)}ref=0.2$ | ||||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | ${\Phi 32(2,1)}ref=0.2$ |

Example . | k = 1
. | k = 2
. | k = 3
. | $\Phi 1j(i,\tau )ref$ . | ${\Phi 2j(i,\tau )}ref$ . | ${\Phi 3j(i,\tau )}ref$ . |
---|---|---|---|---|---|---|

Sign and X^{1}X^{2}arrowdirection | $X1\u2192|a|X2$ | $X1\u2192\u2212|a|X2$ | $X2\u2192|a|X1$ | ${\Phi 12(1,1)}ref=0.8$ | ${\Phi 22(1,1)}ref=\u22120.8$ | ${\Phi 31(2,1)}ref=0.8$ |

${\Phi 11(1,1)}ref=0.2$ | ${\Phi 21(1,1)}ref=0.2$ | ${\Phi 31(1,1)}ref=0.2$ | ||||

${\Phi 12(2,1)}ref=0.2$ | ${\Phi 22(2,1)}ref=0.2$ | ${\Phi 32(2,1)}ref=0.2$ |

CI test . | τ_{max}
. | α
. | α_{PC}
. | N_{K}
. | N_{C}
. | N_{Q}
. | N_{A}
. |
---|---|---|---|---|---|---|---|

ParCorr | 3 | 0.01 | 0.2 | 3 | 40 | 20 | 50 |

CI test . | τ_{max}
. | α
. | α_{PC}
. | N_{K}
. | N_{C}
. | N_{Q}
. | N_{A}
. |
---|---|---|---|---|---|---|---|

ParCorr | 3 | 0.01 | 0.2 | 3 | 40 | 20 | 50 |

Figure 8 shows the results. There are only minimal deviations from the true reference values, which confirms that the proposed method is capable to deal with $NK>2$. This also holds for the summary results over $NR=100$ runs presented in Table VII. Yet, it is important to note that we chose a combination of causal graphs that performed well for $NK=2$, i.e., causal effect changes would also be difficult to detect for $NK=3$.

Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. | $\u03f5^$ . |
---|---|---|---|---|---|---|---|---|---|

4.0 | 0.98 | 1.0 | 0.05 | 0.01 | 0.033 | 0.020 | 10.0 | 7.0 | 0.5 |

Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. | $\u03f5^$ . |
---|---|---|---|---|---|---|---|---|---|

4.0 | 0.98 | 1.0 | 0.05 | 0.01 | 0.033 | 0.020 | 10.0 | 7.0 | 0.5 |

### C. Regime parameter selection

We investigate how parameter selection of the number of regimes affects the results by means of the AICc scores defined in (17). We investigate two test scenarios of ${NK}ref=2,3$ for a selection of the examples defined in Secs. IV A and IV B. The PCMCI parameters are as in Secs. IV A and IV B, while $NR=29$, $NQ=20$ and $NA=20$. The resulting AICc values are displayed in Fig. 9. The $NC$ value is changed adaptively for each $NK$ to ensure a similar $NM$ value for the different number of regimes, i.e.,

for $NK>{NK}ref$. The choice of $NC$ is based on $NM$ which in real life application can be chosen according to the considered processes and data as a good estimate of the timescale of regime changes is often available. Nevertheless, it is also possible to chose $NC$ via an information criterion simultaneously with $NK$ (e.g., in the context of regime-dependent clustering Falkena *et al.*^{40} or regime-dependent Markov regression de Wiljes *et al.*^{38}) The reference value for the number of switches is on average (due to randomization of ${\Gamma}ref$) ${NCref}=40$ for both ${NK}ref=2,3$.

We note that the lowest $NK$ at which the AICc plateaus is the ground-truth one. The plateau itself occurs due to the fact that only the links with non-zero causal effect values are counted toward the number of parameters. Thus, a higher number of regimes $NK$ does not necessarily result in an increase of the total number of parameters. In other words, the penalization is not becoming stronger with higher values of $NK$. Concluding, it is clearly visible that no significant improvement is gained by increasing the number of $NK$ beyond the reference number of regimes. Since the entry point to the plateau reveals the reference number of regimes, it seems possible to face scenarios where the true number of regimes is unknown.

### D. High-dimensional linear network

In this section, the algorithm is evaluated on high-dimensional datasets, with each dataset consisting of $NX=10$ interacting variables. The background regimes are generated with two regular alternating regimes of 300 time steps each, for a total length $T=15000$. The network structures are randomly generated from a family of linear networks defined via the parameters shown in Table VIII, where $L$ is the number of randomly drawn cross-variable links with random coefficients from the third column. Note that each variable is also auto-linked at lag 1 with coefficient randomly drawn from the fourth column. The time series $xt\u2208R10$ are generated with model (18) and for $NR=70$ realizations. Regime-PCMCI is then run with the settings shown in Table IX.

N_{X}
. | L . | $\Phi kj(i,\tau )$ . | $\Phi ki(i,\tau )$ . | max lag . |
---|---|---|---|---|

10 | 30 | [−0.4, 0.4] | [ 0.2, 0.5, 0.9] | 3 |

N_{X}
. | L . | $\Phi kj(i,\tau )$ . | $\Phi ki(i,\tau )$ . | max lag . |
---|---|---|---|---|

10 | 30 | [−0.4, 0.4] | [ 0.2, 0.5, 0.9] | 3 |

CI test . | τ_{max}
. | α
. | α_{PC}
. | N_{K}
. | N_{C}
. | N_{Q}
. | N_{A}
. |
---|---|---|---|---|---|---|---|

ParCorr | 4 | 0.05 | 0.2 | 2 | 49 | 30 | 50 |

CI test . | τ_{max}
. | α
. | α_{PC}
. | N_{K}
. | N_{C}
. | N_{Q}
. | N_{A}
. |
---|---|---|---|---|---|---|---|

ParCorr | 4 | 0.05 | 0.2 | 2 | 49 | 30 | 50 |

The results are shown in Table X, which is structured like Table IV except for TPR and FPR being estimated for the cross-variables links thus focusing on the connections between variables. All links are considered in $\Delta \Phi $. Regime-PCMCI performs very well even in this challenging setting. Notably, individual runs can perform extremely well, with $\Delta \gamma $ reaching as low as $0.02%$, and a total of 53 runs below total average of $\Delta \gamma =11.7%$ (the second row in Table X). The other seven runs are responsible for most of the deviation of the average statistics from the reference values (the first row).

Selection . | Δγ%
. | TPR_{cros}
. | $TPRcrosref$ . | FPR_{cros}
. | $FPRcrosref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. | $\u03f5^$ . | no. of runs . |
---|---|---|---|---|---|---|---|---|---|---|---|

All | 11.7 | 0.94 | 1.0 | 0.18 | 0.08 | 0.059 | 0.005 | 16.0 | 1.5 | 0.85 | 70 |

Δγ < 11.7% | 0.19 | 1.0 | 1.0 | 0.08 | 0.07 | 0.006 | 0.005 | 1.8 | 1.5 | 0.70 | 53 |

Selection . | Δγ%
. | TPR_{cros}
. | $TPRcrosref$ . | FPR_{cros}
. | $FPRcrosref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. | $\u03f5^$ . | no. of runs . |
---|---|---|---|---|---|---|---|---|---|---|---|

All | 11.7 | 0.94 | 1.0 | 0.18 | 0.08 | 0.059 | 0.005 | 16.0 | 1.5 | 0.85 | 70 |

Δγ < 11.7% | 0.19 | 1.0 | 1.0 | 0.08 | 0.07 | 0.006 | 0.005 | 1.8 | 1.5 | 0.70 | 53 |

As in the *causal effect* case, there is a mismatch between runs with the lowest prediction errors $\epsilon ^$ and the lowest error on the regime-assigning process $\Delta \gamma $, meaning that we cannot use a filtering on $\epsilon ^$ to find the best performing runs. This behavior can be explained from the tendency of the algorithm to still over-fit when too many degrees of freedom are available, as well as from the complexity of distinguishing different causal effects (a challenge already manifested in the *causal effect* case).

### E. Computational complexity

Table XI shows some indicators of the performance of the method: the fraction of $NR$ runs that correspond to a (local) minima, the average number of q-iterations needed to reach a local minima and the runtime for the whole $NR$ set of runs (the code run parallel over the $NA$ annealings and using 4–6 CPUs per job).

Example . | No. of local minima/N_{R} (%)
. | Iterations to minima . | Runtime (s)
. |
---|---|---|---|

Arrow direction | 92 | 7 | 600 |

Causal effect | 16 | 13 | 970 |

Lag | 60 | 11 | 1130 |

Sign X^{1} | 52 | 12 | 970 |

Sign X^{1}X^{2} | 70 | 9 | 700 |

Sign X^{1}X^{2} and arrow | 56 | 10 | 2670 |

High dimensional | 92 | 6 | 10 780 |

Example . | No. of local minima/N_{R} (%)
. | Iterations to minima . | Runtime (s)
. |
---|---|---|---|

Arrow direction | 92 | 7 | 600 |

Causal effect | 16 | 13 | 970 |

Lag | 60 | 11 | 1130 |

Sign X^{1} | 52 | 12 | 970 |

Sign X^{1}X^{2} | 70 | 9 | 700 |

Sign X^{1}X^{2} and arrow | 56 | 10 | 2670 |

High dimensional | 92 | 6 | 10 780 |

Most of the examples reach local minima in more than 50% of the $NR$ runs, while the percentage is very low for the *causal effect* (the second column). We note that examples with a high percentage of local minima correspond also to quick convergence in terms of iteration steps (the third column). They are also associated with better regimes reconstruction (see Tables IV, VII, and X), confirming that a clear cost functional minimum (as shown from the second and third columns) is linked to better detection. Finally, the runtime is quite fast: the low dimensional examples take between $10$ and $20$ min for $NK=2$ and 45 min for $NK=3$ to complete $100$ runs. The high-dimensional example takes just below $3$ h for $70$ runs.

## V. A REAL-WORLD EXAMPLE: THE EFFECT OF EL NIÑO SOUTHERN OSCILLATION ON INDIAN RAINFALL

We finally test the performance of Regime-PCMCI on real-world data and apply it to address the non-stationary relationship of El Niño Southern Oscillation (ENSO) and all-India rainfall (AIR) mentioned in Sec. I between the winter and summer months, i.e., the background regimes, and to detect a reported link from ENSO to AIR during summer.

This example can be considered a difficult case since the expected signal from ENSO to AIR is likely small compared to natural variability.^{16} Furthermore, climate data are typically very noisy with causal relationships being diluted by other, often unknown processes given a complex coupled climate system.^{43}

Our input data consist of monthly observations of ENSO and AIR, for the years 1871–2016, resulting in two time series consisting of 1740 monthly values each. More precisely, ENSO is represented by the so-called relative Nino3.4 index provided by the National Oceanic and Atmospheric Administration (NOAA).^{54,55} Data for AIR anomalies (with the climatology subtracted) are provided by the Indian Institute of Tropical Meteorology (IITM).^{56,57}

We choose the following parameters of Regime-PCMCI: for the regime part, we set $NK=2$ and $NC=292$, which is equivalent to assuming two seasons per year. For the PCMCI settings, we use a significance level $\alpha =0.01$ ($\alpha PC=0.2$). Furthermore, we use a maximum time lag of 2 months, i.e., $\tau max=2$. The optimization is run $NA=100$ annealing times, to span many local minima, with each annealing allowed for up to $NQ=100$ iteration steps to converge.

Among the annealing steps, which correspond to different random initial guesses on the regime-assigning process $\Gamma $, some clearly performed better in terms of fitting the data. We estimate the average prediction error associated with each annealing, $\epsilon ^$ (B 4), and Fig. 10(a) shows it for all annealings (ranked according to $\epsilon ^$). A red box highlights the top performing cluster (13 runs).

All of the top 13 annealings find a link from ENSO to AIR during one of their two regimes only (for simplicity hereafter called regime 1). In the following, we present results averaged over these annealings and plot links that surpass a strength of $0.1$.

The causal link from ENSO to AIR in regime 1 has an average standardized linear effect of $\u22120.4$, meaning that a one standard deviation increase in ENSO results in a reduction of $0.4$ standard deviations in AIR [Fig. 10(c)]. This negative dependence is well documented in the literature.^{16} During regime 2, in contrast, ENSO and AIR are, on average, almost independent, with only a very weak link ($\u22120.05$, not shown) detected from AIR to ENSO. More importantly, our results indicate a clear seasonal dependence. Figure 10(d) shows the number of months assigned to each regime (normalized by the number one would expect on the hypothesis of no seasonality, see the figure caption). A clear peak in summer months is found for regime 1. More precisely, most of the months between June and September are assigned to regime 1 ($70%$). These are the months in which the Indian summer Monsoon is active and for which a robust influence from ENSO has been shown. In contrast, months assigned to regime 2 are predominantly winter months ($60%$ of all December to March months). Thus, despite the relatively weak mean causal effect of ENSO on AIR during summer, and the large inter-annual variability, our algorithm successfully reconstructed this well-documented relationship given all-year time series of ENSO and AIR.

A method to detect long-term changes of this summer teleconnection has recently been proposed by Bódai *et al.*^{58} using ensemble climate models. The additional dimension provided by the ensemble members’ allows to compute year-dependent correlations to infer inter-annual changes. In contrast, our method uses a single realization of the dynamics (e.g., observations) to still obtain time-dependent statistics (networks), although in a finite number ($NK$). Note that long-term changes may still be detectable with Regime-PCMCI either as long-term changes to the persistence/start of a regime each year or with the emergence of a regime in a specific time period.

Overall, these results are promising and show the potential of Regime-PCMCI to detect regime-dependent causal structures in a system as complex as the climate system. On the other hand, it also shows that domain knowledge is required to assure a suitable choice of parameters ($NC$ and $NK$) and an interpretation of the results. This is yet a common caveat to many data-driven approaches, which we nevertheless want to stress strongly.

## VI. DISCUSSION AND CONCLUSIONS

Causal discovery is emerging as an important framework across many disciplines in science and engineering, but each discipline has particular challenges that novel methods need to address.^{4} We introduced a novel method, Regime-PCMCI, to learn regime-dependent causal relations, overcoming one of the key drawbacks of current causal recovery methods. The performance of Regime-PCMCI was analyzed for many different artificially generated causal scenarios and for varying regimes showing that the method covers a wide range of settings (see Figs. 2–5, 7 and Table IV). The performance of the algorithm is maintained also for high-dimensional settings with 10 variables (see Table X) as well as for more than two regimes (see Fig. 8 and Table VII). We found limitations of the method for the case where only the causal effect strength of a link changes between regimes (see Fig. 6), which seems to be hard to detect with our optimization scheme and requires further investigation. Furthermore, the capability of Regime-PCMCI was verified by means of a well-documented climate example using real data of ENSO and Indian rainfall (see Fig. 10). Overall, the proposed method presents itself as a promising approach in the context of non-stationary causal links manifested in regime changes in time.

Note that a causal interpretation of estimated links in our observational causal discovery framework still assumes causal sufficiency, that is, no unobserved common causes. However, estimated non-causality (zero coefficients) does not require this assumption^{27} and can be interpreted as an absence of a causal relation already under the weaker faithfulness assumption.^{26} While for PCMCI asymptotic consistency was shown,^{26} this is a more difficult task for Regime-PCMCI and deferred to further research.

There are several interesting aspects that could be explored in the future, building on the present work. These extensions can build on other causal discovery algorithms or extensions of PCMCI in the causal discovery step of our method. For example, the PCMCI algorithm allows for nonlinear causal links^{26} and thus a nonlinear extension of the Regime-PCMCI is a logical next step, e.g., Gaussian processes are used to estimate $gkj$, then the Gaussian Process Distance Correlation (GPDC) test (see Runge *et al.*^{26}) could potentially be used. In that, yet the Regime-PCMCI version would require a different cost functional and optimization approach. Recent extensions of PCMCI to the case of not only lagged, but also contemporaneous causal relations can also be integrated.^{41} Moreover, potentially it is also possible to better capture the *causal effect* case and it might be possible to learn a regime dependence of the noise term.

With respect to applications in climate science, it would be interesting to utilize the proposed method to study other links in the climate system that are likely regime dependent, but less understood than the presented El Niño-Indian rainfall example. Since Regime-PCMCI is formulated in general terms that are not only specific to climate datasets, problems of causal non-stationarity in other application areas could be explored.

## AUTHORS’ CONTRIBUTIONS

E.S., J.d.W., M.K., and J.R. designed the research; E.S. mainly performed the research; and E.S., J.d.W., M.K., and J.R. analyzed the results and wrote the manuscript.

## ACKNOWLEDGMENTS

E.S. was supported by the Centre for Doctoral Training in Mathematics of Planet Earth, UK EPSRC funded (Grant No. EP/L016613/1). M.K. and J.d.W. have been partially funded by the ERC Advanced Grant ACRCC (Grant No. 339390). The research of J.d.W. has been partially funded by Deutsche Forschungsgemeinschaft (DFG) (Grant No. SFB1294/1-318763901) and by the Simons CRM Scholar-in-Residence Program. M.K. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement (No. 841902). The optimization software by Gurobi$TM$ was used for this work. The authors would like to thank Giorgia di Capua for discussions on the climate example.

## DATA AVAILABILITY

The data that support the analysis of the first four sections of this study have been synthetically generated by the authors and can be fully reproduced using the equations and parameters described in the article. The data that support the findings of the last section of this study are openly available in the KNMI Climate Explorer at https://climexp.knmi.nl/, Refs. 55 and 57. PCMCI is part of the open-source Python package tigramite available at https://github.com/jakobrunge/tigramite, Ref. 59.

## NOMENCLATURE

- AIC
Akaike information criterion)

- AICc
Corrected Akaike information criterion

- ENSO
El Niño Southern Oscillation

- FPR
False positive rate

- MCI
Momentary conditional independence

- PCMCI
Causal discovery method

^{26}- RAM
Regime-dependent autoregressive model

- SCM
Structural causal model

- TPR
True positive rate

### APPENDIX A: HETEROGENEOUS NOISE

In the general framework laid out in Eq. (1), the noise variables $\eta tj$ are only assumed to be jointly independent and stationary, each distributed according to a distribution $Dj$. Given that the primary focus of this work is to detect regime-dependent causal structures rather than noise structures, the effective choice for noise distributions used to generate the data is a Gaussian with unit variance $\eta tj\u223cN(0,1)$ for all variables $j$ (Sec. IV).

Yet, this simplification does not necessarily represent the variability of processes in real-world scenarios. Here, the performance of the proposed Regime-PCMCI is exemplified for Gaussian noises with variable-specific variances $\eta tj\u223cN(0,\sigma j2)$ (see Appendix A 1) and noises from two different distributions, Gaussian and uniform (see Appendix A 2).

#### 1. Gaussian noises with variable-specific variances

The data are generated from model (18) with example *sign change $X1X2$* coefficients. The noise terms $\eta j$ are Gaussian distributed with a fixed variance for variable $X1$, $D1=N(0,1)$, and three different cases for variable $X2$, $D2=N(0,\sigma 22)$ with $\sigma 2=0.25$, $\sigma 2=0.5$, and $\sigma 2=2.0$. The Regime-PCMCI results, averaged over 100 different realizations of the regime-assigning processes, are presented in Table XII. The algorithm performs very well in the first two cases (average regime detection error $\Delta \gamma \u22641%$). This is to be expected since a smaller noise in $X2$ allows for a better fit of the data. The latter case is harder to infer since the noise on $X2$ is very large compared to the deterministic signal (here $\Delta \gamma \u224325%$).

Case . | Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. |
---|---|---|---|---|---|---|---|---|---|

σ_{2} = 0.25 | 0.3 | 1.0 | 1.0 | 0.01 | 0.01 | 0.010 | 0.010 | 5.2 | 5.2 |

σ_{2} = 0.5 | 0.8 | 1.0 | 1.0 | 0.02 | 0.01 | 0.013 | 0.013 | 5.3 | 5.2 |

σ_{2} = 2.0 | 24.0 | 0.85 | 1.0 | 0.02 | 0.01 | 0.16 | 0.03 | 36.0 | 9.0 |

Case . | Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. |
---|---|---|---|---|---|---|---|---|---|

σ_{2} = 0.25 | 0.3 | 1.0 | 1.0 | 0.01 | 0.01 | 0.010 | 0.010 | 5.2 | 5.2 |

σ_{2} = 0.5 | 0.8 | 1.0 | 1.0 | 0.02 | 0.01 | 0.013 | 0.013 | 5.3 | 5.2 |

σ_{2} = 2.0 | 24.0 | 0.85 | 1.0 | 0.02 | 0.01 | 0.16 | 0.03 | 36.0 | 9.0 |

#### 2. Different noise distributions

The data are generated from model (18) with example *sign change $X1X2$* coefficients. The noise terms $\eta j$ are set to follow completely different distributions: variable $X1$ is associated with a unit variance Gaussian noise, $D1=N(0,1)$, and variable $X2$ with uniformly distributed noise between $\xb11.5$, $D2=U(\u22121.5,1.5)$. The Regime-PCMCI results, averaged over 100 different realizations of the regime-assigning processes, are presented in Table XIII. This scenario gives results comparable to the ones presented in the paper for the same example, i.e., $\Delta \gamma \u22433%$.

Case . | Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. |
---|---|---|---|---|---|---|---|---|---|

Gauss, Unif | 3 | 1.0 | 1.0 | 0.01 | 0.01 | 0.025 | 0.019 | 8.0 | 7.0 |

Case . | Δγ%
. | TPR_{all}
. | $TPRallref$ . | FPR_{all}
. | $FPRallref$ . | ΔΦ . | ΔΦ^{ref}
. | ΔΦ % . | ΔΦ^{ref} %
. |
---|---|---|---|---|---|---|---|---|---|

Gauss, Unif | 3 | 1.0 | 1.0 | 0.01 | 0.01 | 0.025 | 0.019 | 8.0 | 7.0 |

To summarize, the results show that Regime-PCMCI can deal with specific heterogeneous noise distributions, even belonging to different families of distributions. Since the optimization method acts on regression residuals, we can speculate that we expect good performance as long as the noise terms are not too large in their magnitude and are not too skewed. An elaborate study of these conclusions and an investigation of the potential for generalization of the method to more extreme noise distributions is an interesting research aspect for the future.

### APPENDIX B: DEFINITION OF RESULT STATISTICS

#### 1. Regime-assigning process

#### 2. Link detection

*TPR*

Over the cross-variables links (in Table X):

*FPR*

Over the cross-variables links (in Table X):

#### 3. Link coefficients

can also be computed as average *percentage* error per regime,

#### 4. Prediction error

with **L** defined in Eq. (4).

## REFERENCES

*Nonlinear and Stochastic Climate Dynamics*(Cambridge University Press, 2017), pp. 184–208.

*Proceedings of the 26th International Joint Conference on Artificial Intelligence*(ACM, 2017), pp. 1347–1353.

*Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence*, edited by D. Sontag and J. Peters (AUAI Press, 2020).

*Second International Symposium on Information Theory, Tsahkadsor, Armenia, 2–8 September 1971*, edited by B. N. Petrov and F. Csáki (Akadémiai Kiadó, Budapest, 1973), pp. 267–281.