Reconstructing regime-dependent causal relationships from observational time series

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 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 with a regime learning optimisation approach. PCMCI allows for linear and nonlinear, high-dimensional time series causal discovery. Our method, Regime-PCMCI, is evaluated on a number of numerical experiments demonstrating that it can distinguish regimes with different causal directions, time lags, effects and sign of causal links, as well as changes in the variables' autocorrelation. Further, Regime-PCMCI is employed to observations of El Ni\~no Southern Oscillation and Indian rainfall, demonstrating skill also in real-world datasets.


Introduction
Understanding causal relationships [25,36] among different processes is an ubiquitous task in many scientific disciplines as well as engineering (e.g., in the context of climate research [10,33,20,19,14,31], econometrics [12,22], or molecular dynamics [11]).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 [31].Based on that graph the processes that generated the data can then be modelled in the framework of structural causal models (SCMs) [25] 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 data-driven causal discovery methods, but many challenges emerging from the dynamic nature of such datasets have not yet been met.Further, causal knowledge cannot be gained from data alone and each method comes with its particular set of assumptions [36] about properties of the underlying processes and the observed data.See [31] for an overview of methodological frameworks, challenges, and application scenarios.
A particular and wide-spread challenge is regime-dependence, a common property of nonlinear dynamical systems that can also be described as one form of non-stationary behaviour.Regime-dependence means that the causal relationships between the considered processes vary depending on some prevailing background regime that may be modelled as switching between different states.Further, often such regimes have strong persistence, that is, they operate and affect causal relations on much longer time scales 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 [39,34].It is, however, generally assumed that ENSO does only marginally affect Indian rainfall in winter [24].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 time scale (several months) than the causal relations among ENSO and Indian rainfall (several weeks).

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 [31].In general, causal (network) discovery methods can be classified into classical Granger causality approaches [12,3], constraint-based causal network learning algorithms [36], score-based Bayesian network learning methods [18,6], structural causal models [27], and state-space reconstruction methods [37,2].
Here we focus on the constraint-based framework which has the advantage that it can flexibly account for nonlinear causal relations and different data-types (continuous and categorical, univariate and multivariate).PCMCI adapts this framework to the time series case yielding high detection power also in high-dimensional and strongly autocorrelated time series settings.However, one of the general assumptions of PCMCI (as well as of other causal discovery algorithms) is stationarity, i.e., that at least the existence or absence of a causal link does not change over the considered time series segment [28].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 which constitute a particular case of latent confounding.
Some recent work addresses causal discovery in the presence of non-stationarity.[21] model non-stationarity in the form of (continuous) stochastic trends in a linear autoregressive framework.[41] 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 modelled as latent drivers.In [26,7] assumed known non-stationary regimes are exploited to estimate causal relations also in the presence of general latent confounders.
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 (RAM) were introduced already in 1990 [42].These can yield physically well interpretable results that, however, require well-chosen ancillary variables and a seasonal index which are not learned from data.Thus, RAM requires a priori knowledge of the regimes, which one often aims to learn rather than enforce.Furthermore, the autoregressive framework only permits linear relationships.In the context of discrete state spaces regime dependent causal discovery has been considered in [11].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 [22].Specifically, the regime assignments are computed by sampling from a Markov chain.
A more general framework to handle discrete regimes is the Markov-switching ansatz of [9], which flexibly models regime-dependence utilizing the assumption of a finite number of regimes and a level of persistency in the transitions between different regimes.This ansatz has been successfully realised in combination with many different model assumptions (e.g., see [13]) here we want to explore it for causal networks and combine it with PCMCI [29], a constraint-based time series causal discovery method [36].We call our method Regime-PCMCI.
The remainder of the paper is structured as follows: In section 2 the underlying mathematical problem, concepts, and key assumptions are formalised, and a motivating example is discussed to provide some intuition.Our novel method Regime-PCMCI is then presented in section 3.These theoretical and algorithmic parts are complemented by a thorough numerical investigation of the proposed method in various artificial settings in section 4. Finally, in section 5, Regime-PCMCI is applied to a real-world dataset from climate science, addressing the changing relationships of ENSO and rainfall over India.

Problem setting
Let {X t } t∈Z be a sequence of real-valued N X dimensional random variables X t ∈ R N X where t is associated with time.A realisation over the time interval [0, T ] of this stochastic process is denoted {x t } t∈[0,T ] and we assume that it is possible to obtain observations of these realisations.We assume that the underlying process is modelled by a regime-stationary discrete-time structural causal model (SCM) Here the measurable functions g j t depend non-trivially on all their arguments, the noise variables η j t are jointly independent and are assumed to be stationary, i.e., η j t ∼ D j for all t for some distribution D, and the sets P j t ⊂ (X t−1 , X t−2 , . ..) define the causal parents of X j t .Here we assume lagged relationships, but this is not a necessity.In contrast to approaches assuming stationarity, both g j t and P j t are allowed to depend on regimes in time as further formalized in Assumption 1 (section 2.2).Then the problem setting considered in this manuscript is of the nature of the following inverse problem with ] where g j t belong to an appropriate functions space for each t and j. τ max is the maximum considered time lag.In other words, the aim is to fit a set of unknown parameters Θ t on the basis of an observed time series {x t } t∈[0,T ] .In the next section we will discuss the particular structure of the parameters Θ t we are interested in.

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 [25,36].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 X j t with j = 1, . . ., N X at each time t ∈ Z. Variables X i t−τ and X j t for a time lag τ > 0 and a given t are connected by a lag-specific directed link "X i t−τ → X j t " if X i t−τ ∈ P j t for a particular t.We denote the maximum ground truth time lag of any parent as τ P max .
For a more detailed introduction the reader is referred to [29].In the following we will use graphs and networks interchangeably.
The collection of parent sets for all components at time t is denoted P t = {P 1 t , . . ., P N X t }.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 g j t (P j t , η j t ) in SCM (1) corresponding to these links which we here restrict to an appropriate function class as modelled in Eq. ( 2).If we assume linear functions with coefficients Φ j t , then the inverse problem Eq. (2) simplifies to Thus for a given time series x t ∈ R N X and with t ∈ [0, T ] and functional G t the aim is to find the unknown parameters

Persistence
As mentioned above, in many application areas non-stationarity may be modelled not in form of abrupt or continuous changes, but via piece-wise constant regimes [17,40,11].These regimes will further exhibit a certain persistent behaviour.
In order to capture non-stationary systems with these properties we will restrict our inference to regime-dependent persistent dynamics.Assumption 1: Denote the parents and functional dependency of a given variable j for a regime k as P j t = P j k and g j t (P j t , η j t ) = g j k (P j k , η j t ).We call a regime (N M , N K )-persistent if the parents and functional dependencies are stationary for an average of N M consecutive time steps t.Further, we assume that there is a finite number of regimes on the whole time domain, i.e., k ∈ {1, . . ., N K }.
Note that persistence enters here via a regime average persistence N M , which naturally implies a finite number of regimes N K ≤ T /N M .Under Assumption 1 the considered linear inverse problem (3) reduces to finding a set of parameters and the change points between the regimes given by the regime assigning process with Γ(t) ∈ [0, 1] N K ×T .For example, component k of the regime assigning process can be of the form γ k = (0, 1, 1, 0, 0, 0, 1, . ..) ∈ [0, 1] T .Regime k is active for all time steps for which γ k (t) = 1.

Motivating example
Before we introduce our novel regime detecting 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 X 1 and X 2 and the associated causal graphs as shown in Figure 1,a.Variable X 1 linearly influences X 2 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.The cross-correlation of X 1 and X 2 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 X 1 on X 2 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 (Figure 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 as well as the regime-dependent causal relations.

Method
Our approach is designed to alternate between learning the regimes and 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 [29] as a method that adapts the constraint-based causal discovery framework to the time series case.

Causal discovery
The constraint-based framework has the advantage that it can flexibly account for nonlinear causal relations and different data-types (continuous and categorical, univariate and multivariate) since it is based on conditional independence defined as follows.Two variables X and Y are conditionally independent given a (potentially multivariate ) variable Z, denoted where p denotes associated probability density functions.
There exist a large variety of conditional independence tests, see [28,29] for a discussion.If relationships are assumed linear, as is the case in practical examples of the present work, a partial correlation can be used.
As is explained in detail in [29], PCMCI is based on a variant of the PC algorithm (names after its inventors Peter Spirtes and Clark Glymour [36]) combined with the momentary conditional independence (MCI) test.It consists of two stages: (i) PC 1 condition selection to identify relevant conditions B j t for all time series variables X j t and (ii) the MCI test to test whether X i t−τ → X j t with MCI: Thus, MCI conditions on both the parents of X j t and the time-shifted parents of X i t−τ .These two stages serve the following purposes: PC 1 is a time-lagged causal discovery algorithm based on the PC-stable algorithm [8] that removes irrelevant lagged conditions (up to some τ max ) for each variable by iterative conditional independence testing.A liberal significance level α PC in the tests lets PC 1 adaptively converge to typically only few relevant conditions that include the causal parents with high probability, but might also include some false positives.The MCI test then addresses false positive control for the highly-interdependent time series case, which is why we chose it here.A causal interpretation of the relationships estimated with PCMCI comes from the standard assumptions in the constraint-based framework [36,28,29], namely causal sufficiency, the Causal Markov condition, Faithfulness, non-contemporaneous effects, and stationarity within the regimes as further discussed below.As demonstrated in [29], 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 τ max , and the significance levels α in MCI and α PC in PC 1 .We discuss the selection of these parameters in section 3.6.
PCMCI is applied to sample subsets of the time series pertaining to an estimated regime k in an iterative step of our method, initialised with some random regime assignment.Given a significance level α, the output of PCMCI is the set of parents P k = {P 1 k , . . ., P N X k } for all time series variables for that regime, Based on these parents with associated causal links, causal effects that quantify the strength of a link can be estimated.

Regime learning
Given an estimated set of parents P k for each regime, the regime variable for the next iteration is updated assuming a particular non-stationary setting of finite metastable regimes as defined in Assumption 1 (section 2.2).This learning approach is based on ideas first proposed in [13] and later extended to many different models [9].
In the following we focus on the linear setting, the nonlinear extension is discussed in section 3.5.To learn the regime parameters for the inverse problem (3) introduced in section 2, given P k (the output from PCMCI), we define a cost functional subject to constraints and where d is a distance measure such as the squared euclidean distance • 2 2 and γ is a regime assigning process describing the weight of the individual networks at each time t.
The format of L(Γ, P) 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).Note that the persistence enters the functional in form of a regularization (see Constraint 10).An alternative option is to add a regularisation term that enforces some form of smoothness of Γ (e.g., Tikhonov regularisation [38]).
The free tuning parameter N C is related to the average regime duration of N M time steps as follows: an average regime duration of N M in all N K regimes is implemented by choosing N C ≈ T /(N M N K ).Note that in practice, the average regime switching time N M might not be exactly known.However, we expect in many application areas that prior domain knowledge on reasonable time scales of regime switching is available.The choice of parameters, including choices of value N K , will be discussed in section 3.6.

Algorithm
The Regime-PCMCI algorithm iterates over two major estimation steps: (Step 1) causal discovery to obtain P k and fit the coefficients Φ k (or a nonlinear function) and (Step 2) regime learning to update the regime variable Γ.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.

Step 1: Causal discovery and model estimation
The first step is to estimate a set of parents {P k } (q) and coefficients {Φ k } (q) with k ∈ {1, . . ., N K } on the basis of a fixed {Γ(t)} (q) obtained in step 2 of the previous iteration (see lines of Algorithm 1 and section 3.3.2).In the first iteration, the regimes are assigned randomly.{P k } (q) and {Φ k } (q) are estimated on the basis of a subset of the time series x t with t ∈ {ϒ k } (q) := t : {γ k (t)} (q) ≥ 0.5 (11) for each regime k.The regime-dependent parents {P k } (q) are estimated via PCMCI.
To solve Eq. ( 3), we assume a functional relationship G that relates each variable to its parents P k (for each regime).Here we assume linear functions g j k implying that the coefficients Φ k can be estimated from the following regression model for each fixed k: x for t ∈ {ϒ k } (q) .In other words for every k ∈ {1, . . ., N K } the following optimisation has to be solved for t ∈ {ϒ k } (q) .Note that the coefficients not indicated as relevant via the parent set are defined to be zero, i.e., Φ j k (i, τ) := 0 for X i t−τ / ∈ P j,(q) k .

Step 2: Regime learning
Step 2 is to determine an optimal regime assigning process for the parents and {Φ k } (q) coefficients (see second bullet in q-loop of Algorithm 1).For this the following optimisation problem needs to be solved subject to the constraints ( 9) and (10), and where for each k ∈ {1, . . ., Since the first τ max time steps cannot be predicted, we choose to set those to x j k,t = x j k,t and to not consider this portion of the time series in the algorithm evaluation.
In order to search for the global minimum, the algorithm is run for a number N A of different initializations of {Γ} (0) (annealing).The annealing run with the lowest cost functional objective is chosen as optimal fit.Note that the individual annealing steps are embarrassingly parallelizable.

Reconstruction of time series
A single prediction from Eq. ( 15) can be derived as the weighted sum over k But note this is never used in the code (only (15) via its presence in ( 14) is used).

Nonlinear functions
It is important to mention that the choice of functions g j k in the learning problem (2) should be determined according to the considered applications and on assumptions on the data.Further, the conditional independence test used in PCMCI should cover at least an equally expressive functional dependency class.For example, if Gaussian processes are used to estimate g j k , then the Gaussian Process Distance Correlation (GPDC) test (see [29]) can be used.Consequently, a nonlinear version of the presented Regime-PCMCI would require a different cost functional.The complexity of the assumed model would increase significantly due to the two-fold presence of non-linearity (one through the regime-dependence and the other one via nonlinear causal relations).Therefore, we here restricted the focus to linear functions g j k .Addressing nonlinearity in combination with the considered non-stationarity will be explored in subsequent research.

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 τ max , and the significance levels α in MCI and α PC in PC 1 .α PC should be regarded as a hyper-parameter and can be chosen based on model-selection criteria such as the Akaike Information Criterion (AIC) [1] or cross-validation.τ max could be incorporated into this model selection.But since PCMCI is not very sensitive to this parameter [29] (as opposed to, e.g., Granger causality), its choice can be based on lagged correlation functions, see [29] for a discussion.The choice of conditional independence test is a modelling assumption guided by the assumed nonlinearity of the underlying process and also finite sample considerations.Finally, α is chosen based on the desired level of false positives.
Determining a suitable choice of the unknown number of regimes N K 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) [5] which has been derived in the context of regression models and since been adapted to various other model scenarios including graphs [35].
An IC is designed to capture the goodness of fit penalised 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. ( 17) relates to the number of parameters required to describe Γ which can be fully determined via the change points.The second term in Eq. ( 17) counts the number of relevant parents, that is, the non-zero coefficients Φ j k (i, τ).Here we use the corrected Akaike Information criterion (AICc) first proposed in [15] to estimate N K .Note that we use the corrected version of the original AIC [1] 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 [23] for a more detailed discussion).The choice of N K is numerically investigated in section 4.3.
The number of iteration steps N Q should be chosen to ensure that the optimisation process converges.In our experiments we found with exploratory testings that N Q shows convergence after about 10-20 iterations for all examples investigated.The number of annealing steps N A should be chosen to ensure we can span a large number of local solutions to this non-convex optimisation problem (Eq.( 8)).However, computational time will set a limit to a too high parameter.Note, however, that this part is embarrassingly parallelisable.

Numerical investigation
In the following we investigate the performance of Regime-PCMCI by means of several toy examples.The artificial data is designed to test the methods robustness and accuracy with respect to various potential scenarios that could occur in real applications.At first low dimensional (N X = 2) causal relations are studied as the results can be interpreted more easily.Next, we also consider higher dimensional settings (N X = 10).The reference time series are generated with the following SCM time series model: with predefined {Γ(t)} ref , {Φ k } ref , and {σ 2 } ref .Note that the reference set of parents is specified by the non-zero coefficients {Φ j k (i, τ)} ref .

Low dimensional data with two underlying regimes
First we focus on a simple setting of two regimes, i.e. {N K } ref = 2, and a two dimensional underlying process X t ∈ R 2 (i.e., N X = 2).Our aim is to test the performance of Regime-PCMCI for different elemental features that can change between regimes.For brevity, links X i t−τ → X j t will be called auto-links or auto-dependencies for i = j and cross links for i = j.We consider the following scenarios as summarised in Table 1: 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, which is a realistic yet challenging assumption for many algorithms.

Experiment settings
We design five toy models, in network terms, corresponding to different sets of parents defined via the references parameters {Φ j k (i, τ)} ref given in columns 4 to 5 of Table 1.Further, synthetic regime assigning processes {Γ(t)} ref are generated for all examples.More specifically, {γ 1 (t)} ref is designed to consist of 41 alternating windows, i.e., {N C } ref = 40 regime transitions.The length of these windows is randomly selected to be between 70 and 100 and the constraint (9) imposes The final length of the time series is capped at T = 3, 000 to ensure equally-long regime assignment time series.
Then an artificial time series x t via (19) with {σ 2 } ref = 1 is generated.Note that the stochastic process ( 19) can be exactly reconstructed via the coefficients {Φ j k (i, τ)} ref , their activation {Γ(t)} ref and a specific realisation of the innovation term ε j t .The PCMCI parameters are chosen as follows: partial correlation as a conditional independence test, α = 0.01, α PC = 0.2 as recommended in [32], τ max = 3, and masking type 'y' (see the documentation of tigramite for the definition of masking types).The number of regimes was set to N K = 2 and the maximum number of regime transitions is N C = 40, i.e., correct guess on number of regimes and switches (model selection for N K is investigated in section 4.3).The number of iterations is N Q = 20 and the number of annealings is N A = 50.A summary of the parameters is shown in Table 2.We generate N R = 100 time series realisations for each example.

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 1, focusing on one of the N R synthetic datasets.Table 3 shows summary statistics over all N R runs.
The case sign X 1 X 2 is discussed in detail.The ground-truth regime evolution and networks are shown in the top part of panels a and b in Figure 2; in the middle part of both panels their Regime-PCMCI reconstruction is shown; in the bottom part the difference between reconstructed and true regimes are presented to visually inspect the accuracy.The reconstructed regime assigning process for each regime matches the truth in 99.6% of time steps (97% average value over N R , see Table 3).The corresponding networks have all and only the correct links (TPR = 0.99 and FPR = 0.01 average value over N R ); their linear causal effect is also well estimated with each link correct up to ±0.02 (N R -averaged error per link is 0.028 (9%)).
The other four cases are presented in Figures 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 optimising, which we mitigate by the annealing steps as mentioned in section 3.6.A similar challenge is found for some high dimensional runs for which we refer to section 4.4.
Table 3 shows the summary results over the N R realisations.The estimation errors are presented in terms of the regime assigning process (second column), the network structure (third to sixth column), the causal effects of links (seventh to tenth columns) and the overall reconstructed time series (last column).The second column, ∆γ%, 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).In terms of networks, the link detection performance is evaluated via the true positive (TPR) and false positive rates (FPR).Further, 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 ∆Φ, the average difference between the reconstructed linear coefficient and the reference values of the ground truth links.∆Φ is also expressed as percentage, i.e. each difference is weighted by the absolute value of the ground truth coefficient.The last column, ε, is the expected prediction error per variable and per time step and is computed as ε = L/(N X T ) with L defined in Eq (8) and with N X 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 A.  1.

Example
In summary, Table 3 shows that : • ∆γ%: on average, the regime assigning process is reconstructed correctly in ∼ 94% of the time steps for all cases except causal effect.The causal effect and lag examples are the hardest to infer, with 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, there are 7 instances where ∆γ < 15% (one presented in Figure 6) and those, as expected from PCMCI, give very good network fit.We notice     that these runs do not correspond to the lowest objective values of the N R 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, dynamic wise, are strong enough to be detectable.
• FPR: Ideally the false positive rate should be upper-bounded by α = 0.01.This is also the case if we assume the correct regimes (see column FPR ref ).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.
• ∆Φ%: Errors in parents' detection (either due to false positives (FPR) or to false negatives (missed links, FNR = 1-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 −2 implying a relative error of about 10%.

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 N K > 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).

Experiment settings
The artificial time series is generated via a regime dependent causal graph that is designed by combining two of the regimes settings presented in section 4.1, namely sign X 1 X 2 change and arrow inversion (for details see Table 4).The regime assigning reference process {Γ} 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 5.
Table 4: Artificial model configuration for a low dimensional example with N K = 3 underlying regimes.

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

CI test
Table 5: Method parameters for a low dimensional example with N K = 3 underlying regimes.6: Results for N K = 3 experiments averaged over N R = 100 realisations generated for each example described in Table 4.

Regime parameter selection
We investigate how parameter selection of the number of regimes affects the results by means of the AICc scores defined in (18).We investigate two test scenarios of for The reference value for the number of switches is on average (due to randomisation of We note that the lowest N K 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 towards the number of parameters.Thus a higher number of regimes N K does not necessarily result in an increase of the total number of parameters.In other words the penalisation is not becoming stronger with higher values of N K .Concluding, it is clearly visible that no significant improvement is gained by increasing the number of N K 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.

High dimensional linear network
In this section the algorithm is evaluated on high-dimensional datasets, with each dataset consisting of N X = 10 interacting variables.The background regimes are generated with two regular alternating regimes of 300 time steps each, for a total length T = 15, 000.The network structures are randomly generated from a family of linear networks defined via the parameters shown in Table 7, 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 x t ∈ R 10 are generated with model (19) and for N R = 70 realisations.Regime-PCMCI is then run with the settings shown in Table 8.The results are shown in Table 9, which is structured like Table 3. Regime-PCMCI performs very well even in this challenging setting.Notably, individual runs can perform extremely well, with ∆γ reaching as low as 0.02%, and a total of 53 runs below total average of ∆γ = 11.7% (second row in Table 9).The other 7 runs are responsible for most of the deviation of the average statistics from the reference values (first row).
As in the causal effect case, there is a mismatch between runs with the lowest prediction errors ε and the lowest error on regime-assigning process ∆γ, meaning that we cannot use a filtering on ε to find the best performing runs.This behaviour 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).

Computational complexity
Table 10 shows some indicators of performance of the method: the fraction of N R 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 N R set of runs (the code run parallel over the N A annealings and using 4 to 6 CPUs per job).Most of the examples reach local minima in more than 50% of the N R runs, while the percentage is very low for causal effect (second column).We note that examples with high percentage of local minima correspond also to quick convergence in terms of iterations steps (third column).They are also associated with better regimes reconstruction (see

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 the introduction.We are interested in if, for given time series of ENSO and AIR, our method is able to distinguish 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 [39].Further, climate data is typically very noisy with causal relationships being diluted by other, often unknown processes given a complex coupled climate system [40].
Our input data consist of monthly observations of ENSO and AIR, for the years 1871 to 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) [4] 1 .Data for AIR anomalies (with the climatology subtracted) are provided by the Indian Institute of Tropical Meteorology (IITM) [16] 2 .We choose the following parameters of Regime-PCMCI: For the regime part, we set N K = 2 and N C = 292, which is equivalent to assuming two seasons per year.For the PCMCI settings, we use a significance level α = 0.01 (α PC = 0.2).Further, we use a maximum time-lag of two months, i.e. τ max = 2.The optimisation is run N A = 100 annealing times, to span many local minima, with each annealing allowed for up to N Q = 100 iteration steps to converge.
Among the annealing steps, which correspond to different random initial guesses on the regime-assigning process Γ, some clearly performed better in terms of fitting the data.We estimate the average prediction error associated with each annealing, ε (A.4), and Figure 9,a shows it for all annealings (ranked according to ε).A red box highlights the top performing cluster (13 runs).
All of the top 13 annealing 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 −0.4,meaning that a one standard deviation increase in ENSO results in a reduction of 0.4 standard deviations in AIR (Figure 9,c) .This negative dependence is well documented in the literature [39].During regime 2, in contrast, ENSO and AIR are, on average, almost independent, with only a very weak link (−0.05, not shown) detected from AIR to ENSO.More importantly, our results indicate a clear seasonal dependence.Figure 9,d shows the number of months assigned to each regime (normalised by the number one would expect on the hypothesis of no seasonality, see figure caption).A clear peak in summer months is found for regime 1.More precisely, most of the months between June to 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.
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 (N C and N K ) and an interpretation of the results.This is yet a common caveat to many data-driven approaches, which we nevertheless want to stress strongly.

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 [33].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 analysed for many different artificially generated causal scenarios and for varying regimes showing that the method covers a wide range of settings (see Figures 2-5 and Table 3).The performance of the algorithm is maintained also for high-dimensional settings with 10 variables (see Table 9) as well as for more than two regimes (see Figure 7 and Table 6).We found limitations of the method for the case where only the causal effect strength of a link changes between regimes (see Figure 6) which seems to be hard to detect with our optimisation scheme and requires further investigation.Further, the capability of Regime-PCMCI was verified by means of a well documented climate example using real data of ENSO and Indian rainfall (see Figure 9).Overall, the proposed method presents itself as a promising approach in the context of nonlinear 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) do not require this assumption [28] and can be interpreted as an absence of a causal relation already under the weaker Faithfulness assumption [29].While for PCMCI asymptotic consistency was shown [29], 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, as already discussed in section 3.5, the PCMCI algorithm allows for nonlinear causal links [29] and thus a nonlinear extension of the Regime-PCMCI is a logical next step.Recent extensions of PCMCI to the case of not only lagged, but also contemporaneous causal relations can also be integrated [30].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, it would be interesting to utilise 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.7 Author's contribution E.S., J.dW., M.K. and J.R. designed the research, E.S. mainly performed the research, E.S., J.dW., M.K. and J.R. analyzed the results and wrote the manuscript.

Figure 1 :
Figure 1: Motivating example.(a) Regime dependent ground truth: regime-assigning process and regime-dependent networks.The links are labelled with the associated linear coefficient Φ j k (i, τ) and lag τ.The sign of the coefficient is highlighted by the color (red for positive, blue for negative).(b) Network reconstruction with PCMCI estimated from the whole time series, i.e. if links are wrongly assumed to be stationary.

Figure 2 :
Figure 2: Example case Sign X 1 X 2 .(a) The ground-truth regime-assigning process, {γ} re f (top), the Regime-PCMCI reconstructed process, {γ} reco.(middle) and the difference between the two, ∆γ (bottom).(b) The ground-truth networks for each regime (top), the Regime-PCMCI reconstructed networks (middle) and the difference between the two (bottom).The links are labelled with the associated linear coefficient Φ j k (i, τ) and the lag τ.The sign of the coefficient is highlighted by the color (red for positive, blue for negative).

Figure 4 :
Figure 4: Example case Arrow direction.See description in Figure 2.

Figure 5 :
Figure 5: Example case Lag.See description in Figure 2.

Figure 6 :
Figure 6: Example case Causal effect.See description in Figure 2.

Figure 7 :
Figure 7: Example with N k = 3 regimes for the case Sign X 1 X 2 and arrow direction.See description in Figure 2 but with three regimes.
3 for a selection of the examples defined in the previous sections 4.1 and 4.2.The PCMCI parameters are as in those sections while N R = 29, N Q = 20 and N A = 20.The resulting AICc values are displayed in Figure 8.The N C value is changed adaptively for each N K to ensure a similar N M value for the different number of regimes, i.e., N C

Figure 8 :
Figure 8: Numerical investigation of AICc values for runs with different N K and (a) {N K } ref = 2 for three networks examples (sign X 1 X 2 , arrow and lag change) and (b) {N K } ref = 3 for the sign X 1 X 2 and arrow change example.In each example, individual dots represent the value attained by the N R = 29 runs, and the dashed line goes through the mean values of each set.The vertical grey bar highlights the ground-truth number of regimes {N K } ref .

Figure 9 :
Figure 9: Climate example.(a) Prediction error for each annealing step in ascending order, lowest 13 annealings highlighted in red box.All the other panels refer to this selection.(b) Regime learning: regime-assigning process corresponding to the best annealing (rank 0) (top) and departure from this estimate of the remaining best 12 annealings (in percentage difference).(c) Network learning: mean networks per regime, each causal effect is the mean of of the corresponding coefficient in the individual 13 annealings.(d) Seasonality of the regimes: Number of years per month m assigned to each regime (N k m ), normalised by N * m , which refers to the expected number of months assigned to a given regime if one assumes equal probability 1/N K of assigning a month to one of the two regimes.Thus, here, N * m = 13 • T /(12 * N K )).
Number of assumed regimes N K -Maximum number of transitions within a single regime N C -Maximum time lag τ max -Functional model G -Conditional independence test according to G (e.g.partial correlation for linear G) -Significance level α (and α PC for PC 1 step)

Table 1 :
Artificial model configurations for different low dimensional experiments with N K = 2 underlying regimes.

Table 2 :
Method parameters for low dimensional examples with N K = 2 underlying regimes.

Table 3 :
Results for N K = 2 experiments averaged over N R = 100 realisations generated for each example described in Table

Table 7 :
High dimensional networks parameters

Table 8 :
Method parameters for high dimensional experiments with two underlying regimes.

Table 9 :
Results for high-dimensional experiments over N R = 70 realisations generated for each example described in

Table 3 ,
6,9), confirming that a clear cost functional minimum (as shown from the second and third column) is linked to better detection.Finally, the runtime is quite fast: the low dimensional examples take between 10 and 20 minutes for N K = 2 and 45 minutes for N K = 3 to complete 100 runs.The high dimensional example takes just below 3 hours for 70 runs.

Table 10 :
Summary performance statistics of all examples.The third column is the average value over the respective N R .

Table 11 :
Abbreviations used throughout the manuscript.

Table 12 :
Notation used throughout the manuscript.