# ChronoEpilogi: Scalable Time-Series Variable Selection with Multiple Solutions

Etienne Vareille, Michele Linardi, Ioannis Tsamardinos, Vassilis Christophides

► To cite this version:

Etienne Vareille, Michele Linardi, Ioannis Tsamardinos, Vassilis Christophides. ChronoEpilogi: Scalable Time-Series Variable Selection with Multiple Solutions. 2024. hal-04990168

HAL Id: hal-04990168

https://hal.science/hal-04990168v1

Preprint submitted on 13 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ChronoEpilogi: Scalable Time-Series Variable Selection with Multiple Solutions

Etienne Vareille¹ Michele Linardi¹ Ioannis Tsamardinos² Vassilis Christophides¹

¹ETIS UMR-8051 Laboratory, CY Cergy Paris Universite, ENSEA, CNRS

²Computer Science Department, University of Crete, Heraklion, Greece

## Abstract

We consider the problem of selecting all the minimal-size subsets of multivariate time-series (TS) variables whose past leads to an optimal predictive model for the future (forecasting) of a given target variable (multiple feature selection problem for times-series). Identifying these subsets leads to gaining insights, domain intuition, and a better understanding of the data-generating mechanism; it is often the first step in causal modeling. While identifying a single solution to the feature selection problem suffices for forecasting purposes, identifying all such minimal-size, optimally predictive subsets is necessary for knowledge discovery and important to avoid misleading a practitioner. We develop the theory of multiple feature selection for time-series data, propose the ChronoEpilogi algorithm, and prove its soundness and completeness under two mild, broad, non-parametric distributional assumptions, namely *Compositionality* of the distribution and *Interchangeability* of time-series variables in solutions. Experiments on synthetic and real datasets demonstrate the scalability of ChronoEpilogi to hundreds of TS variables and its efficacy in identifying multiple solutions. In the real datasets, ChronoEpilogi is shown to reduce the number of TS variables by 96% (on average) by conserving or even improving forecasting performance. Furthermore, it is on par with Group Lasso performance, with the added benefit of providing multiple solutions.

## 1 Introduction

In analysis tasks involving TS with hundreds or even tens of thousands of variables (e.g., in manufacturing, environmental monitoring, energy grids, etc.), selecting appropriate series for TS forecasting not only has the potential to improve models' performance, but also to gain intuition, discover new knowledge, and understand the data-generating mechanism (causal structure). *Time-series variable selection (TVS)* is defined as the problem of discovering a *minimal-size subset of the available, measured TS variables whose past values optimally predict the future of a target TS variable T*. Such a series subset is called a **Markov Boundary (MB)** of *T* in the causal discovery literature [Pea88], denoted as **MB(T)**. Hence, **MB(T)** filters out the TS variables that are both *irrelevant* and *redundant* in forecasting *T*.

Identifying a small-size **MB** leads to a computationally more efficient model for *T* when real-time predictions are required. Moreover, it facilitates all subsequent analysis operations such as modeling, explanation calculations (e.g., SHAP values, feature importance), visualizations, and interpretations [VALC23]. Perhaps most importantly, identifying the **MB** of every series could be the first step in *causal modeling* of the data-generating mechanism. Indeed, under broad causal assumptions, the **MB(T)** contains the observable direct causes of *T*, i.e., the quantities that could causally influence and optimize the future values of *T* [LTB+16]. It is not an accident then that variable selection is the first step in several causal discovery algorithms for time-series or cohort data [AST+10b, AST+10a]. However, we note that the **MB(T)** may also contain variables confounded by latent ones (in general, variables connected to *T* by a collider path). Distilling the possible direct causes of *T* out of its
$\mathbf{MB}(T)$ requires further analysis with causal discovery algorithms, which is outside the scope of our work (readers are referred to [ADG22] for a recent survey).

Both in theory and in practice, there may be numerous Markov Boundaries of $T$. This is, for example, the case if a series subset {$X, Y, Z$} is 1-to-1 deterministically related with another subset {$A, B$}, in which case, either subset could substitute for the other in an optimal forecasting model. We call such subsets **informationally equivalent series**. Equivalent series define an equivalence class. The total number of Markov Boundaries of $T$ is exponential to the number of equivalence classes to which the members of $\mathbf{MB}(T)$ belong. This is because we can construct a new $\mathbf{MB}(T)$ by picking any of its subsets and substituting it with an equivalent one from its equivalence class. Determinism is not the only reason for multiple MBs in practice: when the sample size is finite it may be impossible to distinguish the true $\mathbf{MB}(T)$ from some other subset that leads to a forecasting model of statistically indistinguishable performance. The presence of an (often exponential) number of multiple MBs has been established in cohort, cross-sectional data [TCP$^{+22}$, BT21, SLLA13]. For time-series data there have not been any algorithms that return multiple solutions to the TVS problem [YLL21, SLL$^{+15}$].

Identifying all subsets $\mathbf{MB}(T)$ (referred to as the *multiple TVS (MTVS) problem*) is crucial. While finding any one $\mathbf{MB}(T)$ suffices for forecasting, this is insufficient for knowledge discovery and model interpretation. A practitioner may find it misleading to construct an optimal and minimal-size forecasting model by filtering out certain series although there exists another $\mathbf{MB}(T)$ that contains them. In essence, multiple $\mathbf{MB}(T)$s indicate the presence of multiple causal models and explanations that fit the data equally well. Moreover, if two series $X$ and $Y$ are informationally equivalent w.r.t. forecasting $T$, a variable selection algorithm may inconsistently select between them during cross-validation or bootstrapping - arbitrarily returning either $X$ or $Y$ at each step. Feature Importance XAI methods may suffer from similar instability, on top of already known misleading interpretations. For instance, Shapley values currently suffer from the inclusion of unrealistic data instances when features are correlated, even for linear models [AJL21].

In this paper, we define the novel problem of *multiple time-series variable selection* (MTVS), we introduce the concepts of TS informational equivalence and provide a taxonomy to characterize TS variables w.r.t. to their presence in multiple $\mathbf{MB}(T)$: *indispensable*, variables that belong in all $\mathbf{MB}(T)$, *replaceable*, variables that belong in some $\mathbf{MB}(T)$ but could be replaced by some other series in the same equivalence class, *redundant*, variables that are informative but do not belong in any $\mathbf{MB}(T)$, and *irrelevant*, variables that are completely uninformative. We describe the property of *Compositionality* of the data distribution that allows the construction of *sound* yet greedy MTVS algorithms. Additionally, we define the property of *Interchangeability*, which specifies that any two $\mathbf{MB}(T)$s can be decomposed into pairs of equivalent single variables. Under the assumption of Interchangeability, one variable can be substituted with another variable to create a new MB. This property enables the design of algorithms that are *complete* in identifying all MBs without requiring exhaustive enumeration. This is particularly valuable in data distributions with an exponential number of MBs.

We proceed with designing a MTVS algorithm, named **ChronoEpilogi**¹² that is *sound* and *complete* in terms of identifying and returning all MBs of $T$, under Composition, Interchangeability, and other broad, non-parametric assumptions. ChronoEpilogi extends previous variable selection algorithms designed for cohort data [BT21] to a TS forecasting setting and the identification of multiple solutions (MBs). It returns a reference $\mathbf{MB}(T)$ and a list of equivalence classes. Each of the classes contains TS that are informationally equivalent and could substitute one for the other to construct another MB and an equally-performing forecasting model. Experiments on synthetic data demonstrate that ChronoEpilogi has *near perfect average causal f1-score, with stable performance as the MTS dimensionality increases* to thousands of available TS. Moreover, the greedy approximation of the ideal ChronoEpilogi achieves a 70% speedup with only a 0.05 decrease in f1 score for causal discovery (*Claim 1*). Furthermore, ChronoEpilogi variants *are on par with Group Lasso performance* [SFHT13] (*Claim 4*) -arguably the most scalable algorithm for the single solution TVS problem-, with the added benefit of providing multiple solutions in real datasets from the Monash archive [GBW$^{+21}$] like Electricity and Traffic (*Claim 3*). On both real and synthetic datasets the multiple solutions produced by ChronoEpilogi have similar forecasting performance to the unique solution produced by GroupLasso when averaging across all selected targets (*Claim 5*) while when causal ground truth is

<footnote n="1">¹Chronos greek and Epilogi greek</footnote>

<footnote n="2">²https://github.com/ev07/ChronoEpilogi</footnote>
available ChronoEpilogi outperforms GroupLasso in terms of causal f1-score (Claim 2). Finally, both variable selection algorithms actually conserve or even improve the models' performance compared to models trained on the original MTS (Claim 6). We also investigate how SHAP explanations of regression models might misrepresent the role of variables belonging to some but not all Markov boundaries of the modeled target. We claim that the SHAP importance of each equivalent set of variables is distributed among equivalent variables, hence leading to underestimations of the importance of equivalent sets when considered individually (Claim 7). This attribution is unstable: on different data splits, different variables among an equivalence set obtain high importance (Claim 8).

## 2 Preliminaries

We represent a multivariate time series by a set of univariate TS $\mathbf{X} = \{X^1, ..., X^N\}$, where each $X^i \in \mathbb{R}^M$ ($1 \le i \le N$) are regularly sampled observations of a time series variable. We denote with $X_t^i \in \mathbb{R}$ the $t^{th}$ occurrence of $X^i$, and by $X_{t-L,t-1}^i \in \mathbb{R}^L$ the last $L$ occurrences of $X^i$ before the instant $t$. Note that, in the following part, bold letters (e.g., $\mathbf{X}$) refer to sets of TS variables. The TS forecasting task we consider is to predict the values of a target TS variable $T_t$ ($T = X^i$ for some $i$), where $L \le t \le M$, with a model $f$ parameterized by $\Theta$ that uses up until $L$ past timestamps (for each $t$) of the multivariate TS $\mathbf{X}$. Hence, $T_t == f_{\Theta}(X_{t-L,t-1})$.

Conditional independence of two TS variables $X_t$ and $Y_{t'}$ given a third one $Z_{t''}$, written as $X_{t'} \perp Y_{t'} | Z_{t''}$, is defined as $P(X_t, Y_{t'} | Z_{t''}) = P(X_t | Z_{t''}) \cdot P(Y_{t'} | Z_{t''})$. In our work, we assume stationarity, namely the conditional distribution of $T_t$ as a function of $L$ previous lags of $\mathbf{X}$ does not depend on $t$. Stationarity implies temporal consistency, where the conditional independence relations do not rely on $t$. In this respect, we consider the principle of temporal precedence according to which causes (i.e., independent variables) occur before an effect (i.e., outcome). We exclusively test conditional independence relations between a target $T_t$ and variables with timestamps ranging between $t-L$ and $t-1$. Due to stationarity, we can remove the index and denote $X_{t-L,t-1}$ as $\mathbf{X}$.

**Definition 2.1 (Information Equivalence (IEQ))**. Two TS variables $X_a$ and $Y_b$ are information equivalent (shorthand: equivalent) with respect to a target $T_t$ given conditioning TS variable set $Z$, noted $X_a \equiv_T Y_b | Z$, iff they are made independent of $T_t$ by conditioning on the other and on $Z$ itself. Formally: $X_a \equiv_T Y_b | Z \iff X_a \perp T_t | Y_b, Z$ and $Y_b \perp T_t | X_a, Z$. We write shortly $X_a \equiv_T Y_b$ when the conditioning set is $Z = \{\}$.

**Definition 2.2 (Markov Boundary MB($T_t$))**. Given a set of TS variables $\mathbf{X}$ and a target variable $T$ in this set, $M(T_t)$ is a Markov blanket of $T_t$ iff: $T_t \perp (\mathbf{V} \setminus M) | M$. The set $M(T_t)$ is a Markov boundary of $T_t$ iff $\forall \mathbf{M}' \subsetneq \mathbf{M}, T_t \perp (\mathbf{V} \setminus \mathbf{M}') | \mathbf{M}'$.

**Definition 2.3 (Multiple Time-series variable selection (MTVS))**. Let MB($T_t$) a reference MB of a target $T_t$ TS variable. The solution of the MTVS problem consists into finding all MBs, denoted by $\mathcal{M}$. All MB($T_t$) are information equivalent. For all $\mathcal{MB}_i(T_t) \in \mathcal{M}$ the forecasting models $T_t = f_{\Theta}(\mathcal{MB}_i(T_t))$ are equally-performing according to a metric (e.g., $R^2$).

**Example 2.1.** We consider a hypothetical water flow monitoring system involving two rivers A, B and a small hydroelectric station on river B. Let MTS $X, Z \in \mathbb{R}$ be the inflows of the two confluent rivers A and B. The dam is controlled such that if $Z_t \le z_{th}$, it does not produce energy. Otherwise, it diverts a flow $z_{th}$ to power production. The power production flow mixes with river A first for total flow $Y_t = f(Z_t) + X_t$, then with the rest of river B for total flow $T_t = X_t + Z_t$, with $f(Z_t) \in \{0, z_{th}\}$. In this situation, X and Y are deterministically related only when taking account Z, and information equivalent for T given Z only. Both X, Z and Y, Z are MB of T. In this case, solving the MTVS problem relates to discovering both subsets without misidentifying singletons and X, Y as MB.

TS variables selected by a MTVS algorithm can be characterized according to whether they belong to all, at least one, or none of the Markov boundaries of a target $T_t$ [BT21]. A variable $X_a$ is said irreplaceable iff it is part of all information equivalent Markov boundaries of $T_t$, i.e., $X_a \in \mathcal{MB}(T_t) \forall \mathcal{MB}(T_t) \in \mathcal{M}$. A variable $X_a$ is said replaceable iff it is part of at least one $\mathcal{MB}(T_t)$ but not indispensable, i.e., $\exists \mathcal{MB}(T_t), \mathcal{MB}'(T_t) \in \mathcal{M}, X_a \in \mathcal{MB}(T_t) \land X_a \notin \mathcal{MB}'(T_t)$. A variable $X_a$ is redundant iff $\exists Z \subseteq \mathbf{X} \setminus \{X_a\}, X_a \perp T_t | Z$ and $\forall \mathcal{MB}_i(T_t) \in \mathcal{M}, X_a \notin \mathcal{MB}_i(T_t)$ while it is called irrelevant iff $\forall Z \subseteq \mathbf{X} \setminus \{X_a\}, X_a \perp T_t | Z$.

If Causal Markov Condition and Faithfulness [CMJSV$^+$23] hold, each target $T_t$ has the guarantee to have a unique Markov Boundary, which is also the unique solution of the MTVS problem. However, common probability distributions in real settings might violate faithfulness [URBY13], hence the
MB($T_t$) might not be unique. Jointly Gaussian random variables (RV) with a singular covariance matrix or with deterministic relations do no respect faithfulness<fnref n="3" /> [Mee95, Lem07]. While faithfulness does not hold, the weaker *Composition* property relaxes the structure of faithful data. As a matter of fact, faithfulness implies composition, while the opposite relation does not necessarily hold [Pea88].

**Assumption 2.1 (Composition).** For any subset of RVs $X$, $Y$, $T$ and conditioning set $Z$:

$$X \perp T | Z \text{ and } Y \perp T | Z \Rightarrow X \cup Y \perp T | Z \quad (1)$$

Composition is a general property of the joint probability distribution of a set of RVs regardless of their temporal context, hence we dropped index $t$ from formal definitions. The reciprocal property is called *Decomposition*, and it is always true in any probability distribution. Many common probability distributions that violate faithfulness actually satisfy composition.

**Example 2.2.** Consider a TS $X$ of $n$ covariates for which composition holds. Any deterministic transformation of $X$, namely $Y = f(X)$ of size $m$, where $Y_i = f_i(X_{\sigma(i)})$, with $f_i$ invertible and $\sigma_i$ a mapping from $[[1, m]]$ to $[[1, n]]$, also satisfies composition.

**Example 2.3.** Jointly Gaussian distributions [JW07] also satisfy composition, as they are a special case where pairwise independence is equivalent to mutual independence. If independence relations $X \perp T | Z$ and $Y \perp T | Z$ hold, the union $X \cup Y$ contains only RVs that are pairwise independent from each variable in $T$ given $Z$. This implies $X \cup Y \perp T | Z$.

For probability distributions where all information equivalences are caused by invertible deterministic transformations between individual RVs (singletons), the MBs of any target are *interchangeable*: the variables are equivalent regardless of the conditioning set. More generally, we define interchangeability for two Markov Boundaries, where each variable in a MB(T) is equivalent to a variable in the other MB(T) conditioned on the remaining MB(T) variables. We assume that all MB are interchangeable.

**Assumption 2.2 (Interchangeability).** Two MBs of a target $T_t$, $MB(T_t)$ and $MB'(T_t)$ are interchangeable iff:

$$\forall MB(T_t), MB'(T_t) \in \mathcal{M}, \forall X \in MB(T_t), \exists Y \in MB'(T_t), (\{X\} \cup MB'(T_t) \setminus \{Y\}) \in \mathcal{M}$$

**Example 2.4.** In the water flow example 2.1, $Z$ is irreplaceable, and $X$ and $Y$ are replaceable by each other. The markov boundary structure $\mathcal{M} = (\{Z, X\}, \{Z, Y\})$ satisfies Interchangeability, as $\mathcal{M} = \{Z\} \times \{X, Y\}$ is a cartesian product of equivalence classes.

# 3 The ChronoEpilogi Algorithm

In this section, we present the details of our MTVS algorithm, named ChronoEpilogi. We propose two versions of the algorithm: (1) Forward Backward Equivalent (FBE - Algorithm 1) and (2) a computationally optimized version (approximate) named Forward Equivalent (FE - Algorithm 2).

**Algorithm 1 ChronoEpilogi-FBE**

**Require:** TS $X$, target $T$, max lag $L$, threshold params $[\alpha, \gamma, \delta]$

1: set $S \leftarrow \text{FORWARD}(X, T, L, \alpha)$

2: $S \leftarrow \text{BACKWARD}(X, T, L, S, \gamma)$

3: set $\mathcal{M} \leftarrow \text{EQUIV}(X, T, L, S, \delta)$

4: return $\mathcal{M} \triangleright$ set of eq. Markov bound.

**Algorithm 2 ChronoEpilogi-FE**

**Require:** TS $X$, target TS variable $T$, max lag $L$, threshold params $[\alpha, \delta]$

1: $\mathcal{M} \leftarrow \text{FORW-EQUIV}(X, T, L, \alpha)$

2: return $\mathcal{M} \triangleright$ set of eq. Markov bound.

3:

4:

In FBE, we first select informative TS random variables (RV) for predicting the values of a target $T$, adopting a greedy heuristic (FORWARD routine, line 1). Then, a backward phase iteratively removes redundant RVs from the selection (BACKWARD routine, line 2). Lastly, equivalent Markov boundaries are discovered (EQUIV routine, line 3) by checking if any of the selected RVs can be replaced by another one. In ChronoEpilogi-FE, we select equivalent Markov boundaries during informative RV selection (FORW-EQUIV routine, line 1). Such choice permits us to compute conditional independences on smaller TS variable sets. To further reduce time complexity, we also

<footnote n="3">³Conditions leading to faithfulness violations have been studied from the standpoint of probability distributions and graphs [Sad17, ZZM17] as well as applications in homeostatic systems [ZS08], evolutionary systems [And13], gene expression [SLLA13] or spatio-temporal records of meteorological phenomena [YWD+17].</footnote>
**Algorithm 3 FORWARD**

**Require:** TS $\mathbf{X} \in \mathbb{R}^{N \times M}$, target $T$, max lag $L$, stopping threshold $\alpha$

1: set $\mathbf{S} \leftarrow \{T_{1,M-1}\}$

2: model $m \leftarrow$ fit and save $f_{\Theta}$ such that $T_{L+1} = f_{\Theta}(\mathbf{S}_{1,L}), \dots, T_M = f_{\Theta}(\mathbf{S}_{M-L,M-1})$

3: repeat

4: $R \leftarrow m.residuals \in \mathbb{R}^{M-L}$

5: $S_{new} \leftarrow \arg \min_{X' \in \mathbf{X} \setminus \mathbf{S}} \text{Lag-Pearson-pval}(R, X', L) \triangleright \text{Select } S_{new} \text{ s.t. } S_{new} \not\perp R$

6: $\mathbf{S} \leftarrow \mathbf{S} \cup S_{new}$

7: model $m' \leftarrow m$

8: model $m \leftarrow$ fit and save $f_{\Theta}$ on $\mathbf{S}$

9: until pvalue of Likelihood ratio test of (H0: $m' = m$) $\ge \alpha$ $\triangleright$ Verify $S_{new} \not\perp T_{t+1} | \mathbf{S} \setminus S_{new}$

10: return $\mathbf{S} \setminus S_{new}$

**Algorithm 4 BACKWARD**

**Require:** TS $\mathbf{X}$, target $T$, time pred. $t$, max lag $L$, sel. variables $\mathbf{S}$, stopping threshold $\gamma$

1: repeat

2: for $S_{del} \in \mathbf{S}$ do

3: if pvalue of Likelihood ratio test between $f_{\Theta}(\mathbf{S} \setminus \{S_{del}\})$ and $f_{\Theta}(\mathbf{S}) \ge \gamma$ then

4: $\mathbf{S} \leftarrow \mathbf{S} \setminus \{S_{del}\}$

5: until no change in $\mathbf{S}$

6: return $\mathbf{S}$

approximate the search of equivalent Markov boundaries using forecasting model residuals. In Example 2.1, the forward phase might select an upstream redundant variable $U$ first, then $Z$ and $X$ before terminating. The backward phase would test $U \perp T | X, Z$, removing $U$ from the selected set. Finally, the equivalence phase would test $Y \equiv_T X | Z$ and produce the solution space $\mathcal{M}$.

We evaluate ChronoEpilogi considering AutoRegressive Distributed Lags (ARDL) (linear) forecasting model [HPS84, PSS01]. An ARDL model of orders $p, q$ uses lags of both the target ($T$) and other TS RVs ($\mathbf{X}$) as predictors, with $T_t = a_0 + \sum_{i=1}^p a_i T_{t-i} + \sum_{i=1}^k \sum_{j=1}^q b_j X_{t-j}^i + \epsilon_t$. The model includes autoregressive terms ($a_j$) and other explanatory variable terms ($b_j$). It is generally required that $\epsilon_t$ is an iid centered normal noise, but this assumption can be relaxed as we can build estimators in the presence of autocorrelated noise and noise with heteroscedastic components [Whi80, HPS84]. Under sufficient assumptions, an ARDL model can be estimated using Ordinary Least Squares (OLS) procedures, and correspond to a Maximum Likelihood estimation model. Consequently, this estimation has the advantages of the OLS method, with guaranteed convergence and fast computation.

Hereafter, we detail sub-routines of ChronoEpilogi. The *forward* phase (Algorithm 3) iteratively builds a first selected set starting from the past $L$ lags of $T$ (line 1), then incrementally selects new RVs of $\mathbf{X}$. At each iteration, a selected variable ($S_{new}$) is the one maximizing the statistical importance of Pearson correlation between all windows of length $L$ (the predictors) and the residuals of model $m$ (line 5). Algorithm 7 (see Appendix) contains the correlation computation pseudo-code, which iterates all the windows in a candidate variable $X'$ (line 2) to compute p-values of Pearson correlation with the forecasting model residuals. The selection in the FORWARD routine terminates when the last forecasting model (built over $\mathbf{S}$) is not statistically better than the previous one (line 9). We use Likelihood ratio test [Kin98], as they suit ARDL models. The *backward* phase (Algorithm 4) iteratively removes redundant RVs from the selected RVs in the *forward* phase. The main loop of the algorithm tests each one of these RVs. It stops as soon as removing any variable degrades the predictive performance of the *best-so-far* forecasting models (line 3), ensuring that the produced variable set is minimal and equally predictive as the *forward* phase solution. We prove that the *forward* and *backward* phases provide an exact solution of the MTVS problem (See appendix D.1, D.2). The *equivalent* search phase (Algorithm 5) tests the equivalence between each selected TS variable $\mathbf{S}$ (line 2), and any non-selected variable $S^d$ (line 4). In line 5 we compare the equivalence replacing the two Rvs in a new forecasting model. If this latter is comparable with the baseline, we store the equivalent variable sets in dictionary $\mathcal{Q}$ (line 6).
**Algorithm 5 EQUIV**

**Require:** TS $\mathbf{X}$, target $T$, time pred. $t$, max lag $L$, sel. variables $\mathbf{S}$, equivalence threshold $\delta$

1: dictionary $Q \leftarrow \{\cdot\}$

2: **for** $\mathbf{S} \in \mathbf{S}$ **do**

3: $Q[\mathbf{S}] \leftarrow \{\cdot\}$

4: **for** $\mathbf{S}^d \in (\mathbf{X} \setminus \mathbf{S})$ **do**

5: **if** pvalue of LR test between $f_{\Theta}(\{\mathbf{S}^d\} \cup \mathbf{S} \setminus \{\mathbf{S}\})$ and $f_{\Theta}(\{\mathbf{S}^d\} \cup \mathbf{S}) \ge \delta$ **then**

6: $Q[\mathbf{S}] \leftarrow Q[\mathbf{S}] \cup \{\mathbf{S}^d\}$

7: $\mathcal{M} \leftarrow \{S_1\} \cup Q[S_1] \times \dots \times \{S_n\} \cup Q[S_n]$

8: **return** $\mathcal{M}$

**Algorithm 6 FORW_EQUIV**

**Require:** MTS $\mathbf{X}$, target $T$, maximal lag $L$, stopping threshold $\alpha$, equivalence threshold $\delta$

1: set $\mathbf{S} \leftarrow \{T_{1,M-1}\}$

2: set $\mathbf{R} \leftarrow \mathbf{X} \setminus \mathbf{S}$

3: dictionary $Q \leftarrow \{\cdot\}$

4: **model** $m \leftarrow$ fit and save $f_{\Theta}$ such that $T_{L+1} = f_{\Theta}(\mathbf{S}_{1,L}), \dots, T_M = f_{\Theta}(\mathbf{S}_{M-L,M-1})$

5: **repeat**

6: $\textit{Res} \leftarrow \textit{m.residuals} \in \mathbb{R}^{M-L}$

7: $S_{new} \leftarrow \arg \min_{X' \in \mathbf{R}} \textit{Lag-Pearson-pval}(\textit{Res}, X', L)$ ▷ Select $S_{new}$ s.t. $S_{new} \not\perp \textit{Res}$

8: $\mathbf{S} \leftarrow \mathbf{S} \cup \{S_{new}\}$

9: $\mathbf{R} \leftarrow \mathbf{R} \setminus \{S_{new}\}$

10: $Q[S_{new}] \leftarrow \textit{FIND-EQUIVALENCES}(\textit{Res}, S_{new}, \mathbf{R}, \delta)$ ▷ Test $C \equiv_{\textit{Res}} S_{new}$ for $C \in \mathbf{R}$

11: $\mathbf{R} \leftarrow \mathbf{R} \setminus Q[S_{new}]$

12: $m' \leftarrow m$

13: **model** $m \leftarrow$ fit and save $f_{\Theta}$ on $\mathbf{S}$

14: **until** pvalue of Likelihood ratio test of (H0: $m' = m$) $\ge \alpha$ ▷ Verify $S_{new} \not\perp T_{t+1}|\mathbf{S} \setminus S_{new}$

15: **return** $\mathbf{S} \setminus S_{new}, Q$

Once the equivalent Markov Boundaries $\mathcal{M}$ are obtained, we can characterize irreplaceable variables and replaceable variables, as irreplaceable variables are the unique members of their equivalence class. Conversely, non-unique variables in their equivalence class are replaceable. In Section D (see Appendix), we provide soundness and completeness proofs of the multiple solutions computed by Algorithm 5. In Algorithm 6 we report the pseudocode of FORW_EQUIV routine, which is used in FE to select informative variables and to estimate equivalent sets over forecasting model residuals. Algorithm 8 in Appendix details the residual-based equivalence search (FIND-EQUIVALENCES) employed during the forward phase (Algorithm 6, line 10). Such a routine tests statistical redundancy by modeling residuals with model $g_{\Theta}(\cdot)$. In this sense, if a non-selected variable $X$ is equivalent to the tested variable, $S_{new}$ is added to the equivalent set (line 7).

**Complexity analysis** Given a TS $\mathbf{X} \in \mathbb{R}^{N \times M}$, and denoted by $\mathbf{S}$ the selected TS variables in the forward phase, both ChronoEpilogi versions (FBE and FE) run $\mathcal{O}(N|\mathbf{S}|)$ conditional tests. In the FORWARD phase, the sub-routine Lag-Pearson-pval (Algorithm 7) runs in $\mathcal{O}(L n \log(n))$ if Fast Fourier Transform is adopted [MZZ+22], where $L$ is the number of lags to predict a target and $n = M - L$ the size of residuals (number of predicted targets). Fitting an ARDL model through matrix inversion requires $\mathcal{O}(|\mathbf{S}|^3 L^3 M)$. Overall, the forward, backward, equivalent phases take $\mathcal{O}(|\mathbf{S}| N L M + |\mathbf{S}|^4 L^3 M)$, $\mathcal{O}(|\mathbf{S}|^5 L^3 M)$, $\mathcal{O}(|\mathbf{S}|^4 L^3 N M)$ time respectively.

The two variants of ChronoEpilogi, satisfy different properties. FBE (Alg. 1) is an ideal version of our algorithm with provable soundness and completeness under mild assumptions. FE (Alg. 2) is a greedy approximation of FBE. To simplify the proofs of FBE theoretical guarantees we rely on a practical commonly made in causal discovery algorithms that the conditional independence tests are correct <fnref n="4" /> [PNBT07, BT21]. We should stress this requires considering different statistical tests for different distributions [SP20]. We prove that FBE is **sound** as any set of TS variables returned as a solution in $\mathcal{M}$ is a Markov blanket. In other words, our algorithm does not return false positives. We

<footnote n="4">⁴Alternatively, algorithms might rely on an oracle producing variable sets with given properties [SLLA13].</footnote>
then examine the conditions under which FBE is **complete**, i.e. all Markov boundaries are solutions in $\mathcal{M}$ discovered by the algorithm. This corresponds to the lack of false negatives in the equivalent Markov boundary discovery. Theorems proofs are given in Appendix D.

**Theorem 3.1 (Soundness and Completeness of FBE).** *Assuming composition, interchangeable Markov boundaries, and perfect tests for the correlation, termination, elimination and equivalences, FBE computes all Markov boundaries $\mathcal{M}$ of a target $T_t$ and only Markov boundaries.*

The BACKWARD, EQUIV and FORWARD phases rely on model-based independence tests that might be replaced by any conditional independence (CI) test. In other words, CI tests performed in line 5 of the FORWARD algorithm (Algorithm 3) compute correlations over forecasting model residuals. We argue that in the context of linear regression with joint normal variables, measuring the residual correlation is equivalent to a model-based conditional independence test. This is the case of linear models of joint normal variables for which $R_\parallel Z$, and $R$ is normal jointly with all other variables. We notice that residual-based tests are not necessarily equivalent to full conditional tests. Yet no conditional independence remains undetected by a residual test. Formally:

**Theorem 3.2.** *Given $T$ and $X$, two TS variables and $Z$ a conditioning set. Let $R$ be the residual variable corresponding to the regression of $T$ on $Z$. Assume that the modeling achieves independence of the residuals: $R_\parallel Z$, then: $T_\parallel X|Z \Rightarrow R_\parallel X$. Additionally, if for all candidate variable $X$, $R_\parallel X \Rightarrow R_\parallel X|Z$, then $R_\parallel X \Rightarrow T_\parallel X|Z$.*

# 4 Experimental Framework

**Synthetic dataset:** To build synthetic datasets with multiple Markov boundaries that respect composition and interchangeability, we build MTS starting from faithful distributions. Each MTS starts from a Vector AutoRegressive process with 20 variables, with different maximal lags and Markov boundary sizes. We then make copies of some of the MB variables to obtain replaceable variables. Those copies can be randomly shifted forward to change the lag of the relation. The same process is used to obtain redundant variables, copying correlated features not in the original Markov boundary. Irrelevant variables are then added to the dataset by sampling other VAR processes. For each MB size [2, 5, 10], total number of variables [10,100,1000], and maximal lag [1,5,10], we sample 10 datasets per configuration. We verify that the noise intensity varies, as the $R^2$ of a model trained with a MB ranges from 0.02 to 0.9 over the different data instances.

**Real datasets** We evaluate our approach on five forecasting datasets covering different domains: Electricity (consumption), Solar (production), S.F. Traffic [GBW+21], METR-LA, and PEMS-BAY [LYSL18] (transport). They are commonly used in recently proposed deep forecasting models [JZL+23, ZY22, WCWW22]. Electricity has 321 TS and 26304 observations, Traffic 862 TS and 17544 observations, Solar 137 TS and 52560 observations, METR-LA 207 TS and 34272 observations, PEMS-BAY 325 TS and 52116 observations. We evaluate 10 randomly chosen targets per dataset.

**Cross Validation protocol** In TS data, observations are generally not independent, so data splitting for forecasting tasks must ensure that the train split precedes the test split (Forward Chaining Cross Validation). Therefore, we split the dataset along time into a Tuning and Holdout set, respectively for training and evaluating feature selection algorithms and forecasting models. The Tuning set is itself separated into five folds along the time axis, to conduct hyperparameter optimization. We optimize each considered pair of TS selection algorithm and forecasting algorithm together, for maximal average predictive performance $R^2$ over all folds. Due to the consequent training time of deep learning models, it is impractical to tune TSS (Time Series variables Selection) algorithms together. We first tune each TSS algorithm with Support Vector Regression (SVR) models as proxies, then tune deep forecasters for the tuned TSS parameters.

**Baselines and Forecasting models** We compare our algorithm with the only linear scalable baseline, GroupLasso [NST13], and with no selection in the case of Electricity, Traffic, and Solar. For forecasting, we use the ARDL model [PSS01], as it is a standard linear model for MTS data. Real datasets are forecasted using nonlinear models: Support Vector Regression (SVR) with nonlinear kernel, and deep forecasters like DeepAR [SFGJ20] and Temporal Fusion Transformer (TFT) [LALP21]. The input window size is 10 for synthetic MTS and 96 for the five real datasets, similarly to a recent benchmark [SWX+23]. The tuned parameters are described in Appendix (Table 3). We report the performance of the best forecasting model for each target.
<figure><img src="image_2.png" /><figcaption>Figure 1: Multiple solutions of ChronoEpilogi on target T264 of Traffic. MB are build by combining irreplaceable variables (in blue) with members of the replaceable equivalence class (in yellow).</figcaption></figure>

**Evaluation metrics** We evaluate the quality of Markov boundaries based on three criteria: computation time, predictive performance, size, and on synthetic datasets, f1-score using causal ground truth. Predictive performance is measured by $R^2$, RMSE and MAPE. When more than 20 multiple solutions are obtained by ChronoEpilogi, we sample 20 solutions among the solution space, evaluate each of them, and report the mean. When the causal graph is known, we measure the f1-score of identifying respectively irreplaceable and replaceable TS variables. Since GroupLasso returns a unique solution, we also compute the f1-score of identifying all causal variables (irreplaceable and replaceable). Formally, denoting TP the number of correctly detected variables, FN for undetected correct variables, FP for wrongly detected variables, $f1 = \frac{2TP}{2TP+FP+FN}$.

**SHAP explanations** To assess the impact of replaceable variables on SHAP explanations, we train Support Vector Regressors (SVR) and XGBoost Regressors (XGB) for each of the 90 synthetic MTS with 100 covariates (denoted as full MTS), using a window of 10 timesteps. We then compute the SHAP values of random holdout samples using the KernelSHAP approximation [LL17], which we sum along the time axis to obtain the overall importance of variables (per additivity property of SHAP values [KVSF20]). We consider the average of those values as global explanations [BSC$^+$21]. We obtain average SHAP values $s_1^i...s_{100}^i$ corresponding to the global importance of each variable in the MTS number $i$. The total importance of each equivalent variable set $S$ becomes $s_{tot} = \sum_{i \in S} |s_i|$. Then, we compute the contribution percentage of the most important variables (until top-14) in the equivalent set $|s_i/s_{tot}|$. For comparison with SHAP explanations over a unique MB, we create reduced MTS where we remove all but one replaceable variable per equivalent set. For instance, if in MTS $i$ variables $\{1, 4, 16, 78\}$ are equivalent, we remove $\{4, 16, 78\}$ and compute the reduced importance $sr_1^i$ of the representative 1. We obtain instability results by repeating model training and explanation over different (time-contiguous) data splits of each MTS. Stability is measured as proposed by Nogeira et al. [NSB18], applied to the top $|MB|$ variables of each explanation, with $|MB|$ the theoretical size of a Markov Boundary. The metric accounts for covariate size, hence results are comparable between full MTS and reduced MTS.

# 5 Experimental Results

Variant FBE has near perfect average causal f1-score (Table 1), with stable performance as the number of TS increases 2b). **Claim 1:** Variant FE has a 70% speedup compared to FBE with only a decrease of causal f1 metric of 0.05 (Table 1), due to a less precise equivalence discovery phase leading to lower replaceable f1-score (Fig.2b). As the number of TS grows, FE runtime grows slower than FBE (Fig.2a). Given its high comparative effectiveness and efficiency, we decided to use FE on real datasets. **Claim 2:** FBE and FE achieve better results than GroupLasso (GL) w.r.t. causal f1-score, despite similar $R^2$ scores (Table 1). We attribute this difference to GL sensitivity in tuning regularization parameters, compared to ChronoEpilogi thresholds (see Appendix Fig.5). **Claim 3:** FE uncovers multiple MBs in Electricity (for 3 targets) and Traffic (for 7 targets), PEMS-BAY (3 targets), METR-LA (7 targets). In Solar only one MB is found for all targets. The average size of selected TS variables (Table 2) is less than 4% of each MTS. The maximal selected size is below 5% all except Electricity (11%). The solutions have close to identical $R^2$ in Electricity, and spread over 0.02 around the median for Traffic (Fig.2c). Fig.1 illustrates a compact representation of the solution space for one target of Traffic. **Claim 4:** On averaging across all targets the multiple solutions produced by ChronoEpilogi-FE have similar predictive performance to the unique solution produced by GroupLasso (GL). We use the total 50 targets to conduct a Friedman-Nemenyi statistical test on method ranking [Her20], and conclude that GL and FE are not statistically different, but **Claim 6:** variable selection do improves models' performance compared to models trained on the original MTS (full analysis Appendix E.2). **Claim 5:** In Electricity and Solar, GL produces slightly smaller solutions than FE while in the rest datasets, FE produces 8 times smaller solutions than GL (Table 2).
<figure><img src="image_3.png" /><figcaption>Figure 2: Performance of ChronoEpilogi versions on (a)(b) synthetic MTS and (c) real MTS.</figcaption></figure>

<table><caption>Table 1: Computation, predictive and causal performance of tuned ChronoEpilogi variants (FBE, FE) and GroupLasso (GL) on the synthetic dataset, over the 270 synthetic MTS. FE has comparable execution time and predictive power to GL, with a 30% increase in causal f1-score.</caption><thead><tr><th></th><th>Time(s)</th><th>R<sup>2</sup></th><th>MB size</th><th>Number of MB</th><th>causal f1-score</th></tr></thead><tbody><tr><td>FBE</td><td>375±772</td><td>0.373±0.236</td><td>5.69±3.33</td><td>8.87e+05±3.7e+06</td><td>0.99±0.064</td></tr><tr><td>FE</td><td>118±178</td><td>0.373±0.236</td><td>5.74±3.31</td><td>6.87e+07±7.2e+08</td><td>0.94±0.188</td></tr><tr><td>GL</td><td>120±149</td><td>0.337±0.239</td><td>7.50±10.2</td><td>NA</td><td>0.57±0.349</td></tr></tbody></table>

Note that computing one solution with FE forward phase is up to two orders of magnitude faster than GL. Table 2 reports the best forecaster after careful hyperparameter tuning per target.

**Shapley values instability over replaceable variables** **Claim 7:** Models tend to select one or a few variables as important predictors. For XGB models, we report in Fig. 9 (see Appendix) the average contribution of the top important variables for each equivalent set of variables. Clearly, **importance is increasingly shared as the number of equivalent variables grows. Claim 8:** The most important variable of an equivalent set is highly unstable to data resampling. The stability [NSB18] of the top ranked variable is respectively 0.05±0.11 and 0.02±0.09 for XGB and SVR, where -0.05 means total randomness and 1 deterministic selection. Additionally, the presence of replaceable variables impacts the stability of the important variables of the entire explanation. The top $|MB(T)|$ ranked variables overall over full MTS are significantly more unstable compared to reduced MTS. We applied paired Wilcoxon signed rank test and concluded that top variables in full MTS are less stable than when a unique MB(T) is left, with p-value 4.2e-13.

## 6 Related work

Most of the state-of-the-art variable selection algorithms for MTS have focused on selecting a unique concise solution. Scalable methods include using bivariate linear VAR models to select causally related pairs of TS, without considering multivariate interactions [SLL+15]. There are mRMR approaches for MTS [HRL15] with Dynamic Time Wrapping distances [RGFO17], but with a quadratic complexity in the number of TS due to distance computations, which are proven to be impractical for large dataset sizes. We selected GroupLasso [NST13] as baseline in our experiments, as it can specify groups of TS [LNZ15][HMNK10] while applying Lasso-type optimization with linear modeling. It is worth noticing that Lasso-type algorithms have been extended to identify multiple solutions in i.i.d data [PLCT18] without, however, proving any formal property.

As a matter of fact, the problem of computing multiple solutions for the variable selection problem is still in its infancy. In particular, we are not aware of any algorithm applied to MTS data. Theoretically sound methods focus on different distributional properties and assumptions. TIE* [SLLA13] assumes that a single Markov boundary discovery algorithm for non-faithful data can be called as an oracle. All MB(T) are discovered using exhaustive combinatorial conditional tests ($O(s.n^s)$) where s is the maximal size of a considered MB and n the number of variables [BT21]). TMFBS [BT21] proposes forward-backward MB discovery algorithms that aims to decrease the number of redundant
<table><caption>Table 2: Forecasting performance of ChronoEpilogi (FE), GroupLasso (GL) and No Selection (NS) over real datasets. We report the number of times each model was selected (TFT/DeepAR/SVR), and the time spent in the forward/equivalence phases (time F/E). FE multiple solutions are on par to the performance of the unique GL solution while their size in Traffic is 8 time smaller than GL.</caption><thead><tr><th>Dataset</th><th>TSS</th><th>R<sup>2</sup> ↑</th><th>rmse ↓</th><th>mape ↓</th><th>size</th><th>time F/E</th><th>#model</th><th>#MB</th></tr></thead><tbody><tr><td>Electricity</td><td>FE</td><td><strong>0.940</strong></td><td><strong>0.226</strong></td><td><strong>2.388</strong></td><td>10.9</td><td><strong>31</strong> /2411</td><td>6/0/4</td><td>1e+14</td></tr><tr><td>Electricity</td><td>GL</td><td>0.934</td><td>0.236</td><td>2.649</td><td><strong>6.5</strong></td><td>127</td><td>3/0/7</td><td>NA</td></tr><tr><td>Electricity</td><td>NS</td><td>0.863</td><td>0.349</td><td>4.935</td><td>321.0</td><td>NA</td><td>0/0/10</td><td>NA</td></tr><tr><td>Solar</td><td>FE</td><td><strong>0.985</strong></td><td><strong>0.109</strong></td><td>0.664</td><td>5.1</td><td><strong>3</strong> /261</td><td>0/0/10</td><td>1</td></tr><tr><td>Solar</td><td>GL</td><td><strong>0.984</strong></td><td><strong>0.111</strong></td><td><strong>0.623</strong></td><td><strong>4.7</strong></td><td>119</td><td>0/0/10</td><td>NA</td></tr><tr><td>Solar</td><td>NS</td><td>0.968</td><td>0.159</td><td>1.607</td><td>138.0</td><td>NA</td><td>0/0/10</td><td>NA</td></tr><tr><td>Traffic</td><td>FE</td><td>0.783</td><td>0.442</td><td>39.745</td><td><strong>11.7</strong></td><td><strong>12</strong> / 1248</td><td>6/0/4</td><td>7e+24</td></tr><tr><td>Traffic</td><td>GL</td><td><strong>0.797</strong></td><td><strong>0.431</strong></td><td><strong>28.130</strong></td><td>79.1</td><td>211</td><td>3/0/7</td><td>NA</td></tr><tr><td>Traffic</td><td>NS</td><td>0.740</td><td>0.491</td><td>42.371</td><td>863.0</td><td>NA</td><td>1/0/9</td><td>NA</td></tr><tr><td>PEMS-BAY</td><td>FE</td><td>0.860</td><td>0.358</td><td><strong>1.208</strong></td><td><strong>4.1</strong></td><td><strong>143</strong> /6250</td><td>8/0/2</td><td>6e+4</td></tr><tr><td>PEMS-BAY</td><td>GL</td><td><strong>0.867</strong></td><td><strong>0.355</strong></td><td>1.217</td><td>37.7</td><td>765</td><td>10/0/0</td><td>NA</td></tr><tr><td>PEMS-BAY</td><td>NS</td><td>0.820</td><td>0.957</td><td>1.343</td><td>325</td><td>NA</td><td>9/0/1</td><td>NA</td></tr><tr><td>METR-LA</td><td>FE</td><td>0.886</td><td>0.374</td><td>1.121</td><td><strong>3.8</strong></td><td><strong>100</strong> /2246</td><td>8/0/2</td><td>2e+4</td></tr><tr><td>METR-LA</td><td>GL</td><td>0.864</td><td>0.401</td><td>1.230</td><td>41.9</td><td>306</td><td>7/0/3</td><td>NA</td></tr><tr><td>METR-LA</td><td>NS</td><td><strong>0.896</strong></td><td><strong>0.363</strong></td><td><strong>1.054</strong></td><td>207</td><td>NA</td><td>10/0/0</td><td>NA</td></tr></tbody></table>

conditional independence tests operated by TIE* to $O(n^s)$. KIAMB is an iterative forward-backward algorithm [PNBT07] that requires composition and heuristic selections to find one MB of a target, with every MB having a non-zero probability of selection. Identifying all MBs might require an exponential number of runs of KIAMB [SLLA13]. In construct, ChronoEpilogi computes provable equivalent subsets of variables in $O(ns)$ conditional tests under composition and interchangeability assumptions, while the forward-backward phase is based on a heuristic. We should mention that there also exist related works that rely on pure associational criteria where heuristics are used without any causal guarantee, especially in the field of gene expression [LLZ10a] [LLZ10b] (see [SLLA13] for an in-depth review). Finally, several works highlight the importance of a structured representation of multiple MBs for interpretability concerns, especially in the presence of an exponential number of MBs of a target of interest [TLT18][LPL$^+$23]. ChronoEpilogi is the first algorithm that provides a compact representation of mutually equivalent variables for MTS. Redundant and irreplaceable variables can be easily distinguished at first glance (see Figure 1).

# 7 Conclusions

Overall, the paper's contributions are (a) the adaptation to time series and further development of the theory of multiple variable selection, (b) the design of the first variable selection algorithm for time-series data, called ChronoEpilogi, that scales to thousands of available series, (c) the conditions of soundness and completeness of ChronoEpilogi, and (d) the empirical evaluations of ChronoEpilogi demonstrating the presence of multiple solutions in real data, achieving on par performance against Group Lasso and no selection, while reducing the number of TS required to build the model by 96% (on average) by conserving or even improving forecasting performance. The reduced model could be employed as the final model for production, or as a surrogate to a model using all available time series (assuming it is better performing) to facilitate interpretation, visualization, and explanations. Finally, we leave as future work the construction of ensemble models trained on several or all MBs as a means to create forecasters more robust to noise, faulty sensors, or other systematic errors.

We note that the conclusions are limited to the scope of the experimental study. The latter could benefit from a larger scope of experiments with more real and synthetic datasets of varying size, statistical properties, and data types (e.g., discrete time-series). In addition, different variants of the main algorithm could employ non-linear models to compute residuals and non-linear correlation methods for heuristically selecting time series. The implementation has not been optimized at the low level to reach higher computational gains. Other distributional properties that could lead to greedy yet sound and complete multiple time series selection algorithms could be explored. We also hope to relax the assumption of stationarity in future works, as most practical applications have trends/seasonality, different train-test distributions, or change points.
Acknowledgements

This work has received funding from the Horizon Europe Framework Programme under Grant agreement No 101135775 (PANDORA) and ANR under the Grant agreement 24-CE23-6509 (AIDA).

References

[ADG22] Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. *J. Artif. Int. Res.*, 73, may 2022.

[AJL21] Kjersti Aas, Martin Jullum, and Anders Løland. Explaining individual predictions when features are dependent: More accurate approximations to shapley values. *Artificial Intelligence*, 298:103502, 2021.

[And13] Holly Andersen. When to expect violations of causal faithfulness and why it matters. *Philosophy of Science*, 80(5):672–683, 2013.

[AST+10a] C.F. Aliferis, A.R. Statnikov, I. Tsamardinos, S. Mani, and X.D. Koutsoukos. Local causal and markov blanket induction for causal discovery and feature selection for classification part i: Analysis and extensions. *The Journal of Machine Learning Research*, 2010.

[AST+10b] C.F. Aliferis, A.R. Statnikov, I. Tsamardinos, S. Mani, and X.D. Koutsoukos. Local causal and markov blanket induction for causal discovery and feature selection for classification part ii: Analysis and extensions. *The Journal of Machine Learning Research*, 2010.

[BSC+21] João Bento, Pedro Saleiro, André F. Cruz, Mário A.T. Figueiredo, and Pedro Bizarro. Timeshap: Explaining recurrent models through sequence perturbations. In *Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining*, KDD '21, page 2565–2573, New York, NY, USA, 2021. Association for Computing Machinery.

[BT21] Giorgos Borboudakis and Ioannis Tsamardinos. Extending greedy feature selection algorithms to multiple solutions. *Data Mining and Knowledge Discovery*, 35(4):1393–1434, Jul 2021.

[CMJSV+23] Manuel Castro, Pedro Ribeiro Mendes Júnior, Aurea Soriano-Vargas, Rafael de Oliveira Werneck, Maiara Moreira Gonçalves, Leopoldo Lusquino Filho, Renato Moura, Marcelo Zampieri, Oscar Linares, Vitor Ferreira, Alexandre Ferreira, Alessandra Davólio, Denis Schiozer, and Anderson Rocha. Time series causal relationships discovery through feature importance and ensemble models. *Scientific Reports*, 13(1):11402, Jul 2023.

[GBW+21] Rakshitha Godahewa, Christoph Bergmeir, Geoffrey I. Webb, Rob J. Hyndman, and Pablo Montero-Manso. Monash time series forecasting archive, 2021.

[Her20] Steffen Herbold. Autorank: A python package for automated ranking of classifiers. *Journal of Open Source Software*, 5(48):2173, 2020.

[HMNK10] Stefan Haufe, Klaus-Robert Müller, Guido Nolte, and Nicole Krämer. Sparse causal discovery in multivariate time series. In *causality: objectives and assessment*, pages 97–106. PMLR, 2010.

[HPS84] David F. Hendry, Adrian R. Pagan, and J.Denis Sargan. Chapter 18 dynamic specification. volume 2 of *Handbook of Econometrics*, pages 1023–1100. Elsevier, 1984.

[HRL15] Min Han, Weijie Ren, and Xiaoxin Liu. Joint mutual information-based input variable selection for multivariate time series modeling. *Engineering Applications of Artificial Intelligence*, 37:250–257, 2015.

[JW07] Richard A Johnson and Dean W. Wichern. *Applied Multivariate Statistical Analysis*. Pearson, 6th edition, 2007.
[JZL<sup>+</sup>23] Ming Jin, Yu Zheng, Yuan-Fang Li, Siheng Chen, Bin Yang, and Shirui Pan. Multivariate Time Series Forecasting With Dynamic Graph Neural ODEs. *IEEE Transactions on Knowledge and Data Engineering*, 35(9):9168–9180, September 2023. Conference Name: IEEE Transactions on Knowledge and Data Engineering.

[Kin98] Gary King. *Unifying Political Methodology: The Likelihood Theory of Statistical Inference*. University of Michigan Press, 1998.

[KVSF20] I. Elizabeth Kumar, Suresh Venkatasubramanian, Carlos Scheidegger, and Sorelle A. Friedler. Problems with shapley-value-based explanations as feature importance measures. In *Proceedings of the 37th International Conference on Machine Learning*, ICML’20, 2020.

[LALP21] Bryan Lim, Sercan Ö. Arık, Nicolas Loeff, and Tomas Pfister. Temporal fusion transformers for interpretable multi-horizon time series forecasting. *International Journal of Forecasting*, 37(4):1748–1764, 2021.

[Lem07] Jan. Lemeire. *Learning causal models of multivariate systems: and the value of it for the performance modeling of computer programs*. Vrije Universiteit Brussel, Brussels, Belgium, 2007.

[LL17] Scott M. Lundberg and Su-In Lee. A unified approach to interpreting model predictions. In *Proceedings of the 31st International Conference on Neural Information Processing Systems*, NIPS’17, page 4768–4777, Red Hook, NY, USA, 2017. Curran Associates Inc.

[LLZ10a] Huawen Liu, Lei Liu, and Huijie Zhang. Ensemble gene selection by grouping for microarray data classification. *Journal of Biomedical Informatics*, 43(1):81–87, feb 2010.

[LLZ10b] Huawen Liu, Lei Liu, and Huijie Zhang. Ensemble gene selection for cancer classification. *Pattern Recognition*, 43(8):2763–2772, aug 2010.

[LNZ15] Yanming Li, Bin Nan, and Ji Zhu. Multivariate sparse group lasso for the multivariate multiple linear regression with an arbitrary group structure. *Biometrics*, 71(2):354–363, 2015.

[LPL<sup>+</sup>23] Kleanthi Lakiotaki, Zaharias Papadovasilakis, Vincenzo Lagani, Stefanos Fafalios, Paulos Charonyktakis, Michail Tsagris, and Ioannis Tsamardinos. Automated machine learning for genome wide association studies. *Bioinformatics*, 39(9):545, sep 2023.

[LTB<sup>+</sup>16] Vincenzo Lagani, Sofia Triantafillou, Gordon Ball, Jesper Tegner, and Ioannis Tsamardinos. Probabilistic computational causal discovery for systems biology. *Uncertainty in Biology*, 2016.

[LYSL18] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. In *International Conference on Learning Representations (ICLR '18)*, 2018.

[Mee95] Christopher Meek. Strong completeness and faithfulness in bayesian networks. In *Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence*, UAI’95, page 411–418, San Francisco, CA, USA, 1995. Morgan Kaufmann Publishers Inc.

[MZZ<sup>+</sup>22] Abdullah Mueen, Sheng Zhong, Yan Zhu, Michael Yeh, Kaveh Kamgar, Krishnamurthy Viswanathan, Chetan Gupta, and Eamonn Keogh. The fastest similarity search algorithm for time series subsequences under euclidean distance, Aug 2022.

[NSB18] Sarah Nogueira, Konstantinos Sechidis, and Gavin Brown. On the stability of feature selection algorithms. *Journal of Machine Learning Research*, 18(174):1–54, 2018.

[NST13] Trevor Hastie Noah Simon, Jerome Friedman and Robert Tibshirani. A sparse-group lasso. *Journal of Computational and Graphical Statistics*, 22(2):231–245, 2013.
[Pea88] Judea Pearl. *Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference*. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988.

[PLCT18] Yannis Pantazis, Vincenzo Lagani, Paulos Charonyktakis, and Ioannis Tsamardinos. Multiple equivalent solutions for the lasso, 2018.

[PNBT07] Jose M. Peña, Roland Nilsson, Johan Björkegren, and Jesper Tegnér. Towards scalable and data efficient learning of Markov boundaries. *International Journal of Approximate Reasoning*, 45(2):211–232, July 2007.

[PSS01] M. Hashem Pesaran, Yongcheol Shin, and Richard J. Smith. Bounds testing approaches to the analysis of level relationships. *Journal of Applied Econometrics*, 16(3), 2001.

[RGFO17] Milos Radovic, Mohamed Ghalwash, Nenad Filipovic, and Zoran Obradovic. Minimum redundancy maximum relevance feature selection approach for temporal gene expression data. *BMC Bioinformatics*, 18(1):9, January 2017.

[Sad17] Kayvan Sadeghi. Faithfulness of probability distributions and graphs. *J. Mach. Learn. Res.*, 18(1):5429–5457, jan 2017.

[SFGJ20] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. *International Journal of Forecasting*, 36(3):1181–1191, 2020.

[SFHT13] Noah Simon, Jerome Friedman, Trevor Hastie, and Robert Tibshirani. A sparse-group lasso. *Journal of computational and graphical statistics*, 22(2):231–245, 2013.

[SLL+15] Youqiang Sun, Jiuyong Li, Jixue Liu, Christopher W. K. Chow, Bing-Yu Sun, and Rujing Wang. Using causal discovery for feature selection in multivariate numerical time series. *Mach. Learn.*, 101(1-3):377–395, 2015.

[SLLA13] Alexander Statnikov, Nikita I. Lytkin, Jan Lemeire, and Constantin F. Aliferis. Algorithms for Discovery of Multiple Markov Boundaries. *Journal of machine learning research : JMLR*, 14:499–566, feb 2013.

[SP20] Rajen D. Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. *The Annals of Statistics*, 48(3):1514–1538, jun 2020. Publisher: Institute of Mathematical Statistics.

[SWX+23] Zezhi Shao, Fei Wang, Yongjun Xu, Wei Wei, Chengqing Yu, Zhao Zhang, Di Yao, Guangyin Jin, Xin Cao, Gao Cong, et al. Exploring progress in multivariate time series forecasting: Comprehensive benchmarking and heterogeneity analysis. *arXiv preprint arXiv:2310.06119*, 2023.

[TCP+22] Ioannis Tsamardinos, Pavlos Charonyktakis, Georgios Papoutsoglou, Giorgos Borboudakis, Kleanthi Lakiotaki, Jean-Claude Zenklusen, Hartmut Juhl, Ekaterini Chatzaki, and Vincenzo Lagani. Just add data: automated predictive modeling for knowledge discovery and feature selection. *NPJ Precis Oncol*, 6, jun 2022.

[TLT18] Michail Tsagris, Vincenzo Lagani, and Ioannis Tsamardinos. Feature selection for high-dimensional temporal data. *BMC Bioinformatics*, 19(1):17, jan 2018.

[URBY13] Caroline Uhler, Garvesh Raskutti, Peter Bühlmann, and Bin Yu. Geometry of the faithfulness assumption in causal inference. *The Annals of Statistics*, 41(2), April 2013. arXiv:1207.0547 [math, stat].

[VALC23] Etienne Vareille, Adel Abbas, Michele Linardi, and Vassillis Christophides. Evaluating explanation methods of multivariate time series classification through causal lenses. In *10th IEEE International Conference on Data Science and Advanced Analytics, DSAA 2023, Thessaloniki, Greece, October 9-13, 2023*, pages 1–10. IEEE, 2023.

[WCWW22] Chunnan Wang, Xingyu Chen, Chengyue Wu, and Hongzhi Wang. AutoTS: Automatic Time Series Forecasting Model Design Based on Two-Stage Pruning, March 2022. arXiv:2203.14169 [cs].
[Whi80] Halbert White. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. *Econometrica*, 48(4):817–838, 1980.

[YLL21] Kui Yu, Lin Liu, and Jiuyong Li. A unified view of causal and non-causal feature selection. *ACM Trans. Knowl. Discov. Data*, 15(4), apr 2021.

[YWD$^{+}$17] Kui Yu, Xindong Wu, Wei Ding, Yang Mu, and Hao Wang. Markov Blanket Feature Selection Using Representative Sets. *IEEE Transactions on Neural Networks and Learning Systems*, 28(11):2775–2788, November 2017. Conference Name: IEEE Transactions on Neural Networks and Learning Systems.

[ZS08] Jiji Zhang and Peter Spirtes. Detection of unfaithfulness and robust causal inference. *Minds and Machines*, 18:239–271, 2008.

[ZY22] Yunhao Zhang and Junchi Yan. Crossformer: Transformer Utilizing Cross-Dimension Dependency for Multivariate Time Series Forecasting. sep 2022.

[ZZM17] Zhalama, Jiji Zhang, and Wolfgang Mayer. Weakening faithfulness: some heuristic causal discovery algorithms. *International journal of data science and analytics*, 3:93–104, 2017.
# A Additional information on preliminaries

**Information Equivalence** We remark that the equivalence between variables properly defines an equivalence relation. From its definition, it is straightforwardly reflexive ($X \equiv_T X | Z$), symmetric ($X \equiv_T Y | Z \implies Y \equiv_T X | Z$) and transitive.

## A.1 Properties

Independence relationships validate certain properties in any distribution. Contraction and weak union [Pea88] link independences of sets of variables to conditional independences:

**Lemma A.1 (Weak Union).** Let $X, Y, Z, W$ be sets of variables, then $X \perp Y \cup W | Z \implies X \perp Y | Z \cup W$.

**Lemma A.2 (Contraction).** Let $X, Y, Z, W$ be sets of variables, then $X \perp Y | Z$ and $X \perp W | Z \cup Y \implies X \perp Y \cup W | Z$.

Equivalent variables conditioned on a certain set have the same conditional information on the target. They also contain the same information on the target when considering only information that they have in common with another variable $B$. Precisely, the following property holds for any probability distribution:

**Lemma A.3 (Lemiere [Lem07]).** Let $X, Y, T, B$ be random variables and $Z$ a conditioning set. If $X \equiv_T Y | Z$, then:

$$B \perp T | Z, X \iff B \perp T | Z, Y \qquad (2)$$

Joint normal distributions have a specific independence structure: pairwise independence implies mutual independence in such distributions, which is not the case generally.

**Lemma A.4.** [JW07] Let $U, V$ be two sets of random variables such that $U \cup V$ is joint normal. Then:

$$U \perp V \iff \forall U \in U, \forall V \in V, U \perp V$$

When the covariance matrix of a joint gaussian distribution is invertible, the distribution is faithful (factorising into a bayesian network). In this circumstance, there is a unique Markov boundary. For equivalence relations to exist in jointly gaussian distributions, the covariance matrix must be singular.

## A.2 Assumptions

**Composition** Equivalent search algorithms [SLLA13], [PNBT07] can alternatively rely on a local version of composition holding only for the target variable.

**Assumption A.1 (Local Composition [SLLA13]).** The local composition assumption with respect to a target $T$ holds when for all $X, Y, Z$, when Eq.1 holds with $T = T$.

**Limitations of Composition** Composition is essentially a weaker version of "pairwise independence". Mutual independence would require that the joint probability $P[X, Y, T]$ can be factored into $P[X]P[Y]P[T]$ when $X, Y, T$ are mutually independent. Instead, composition requires that when $X, T$ and $Y, T$ are pairwise independent, the joint probability $P[X, Y, T]$ factorizes as $P[X, Y]P[T]$.

The typical example of a XOR operation on two Bernoulli variables $X, Y$ with probability 0.5, and $T = X \text{xor} Y$ invalidates composition. Indeed, $X$ and $Y$ individually bring no information on $T$, while together define entirely $T$.

Another way of seeing composition is regarding interaction among variables. In clinical trials, it is often the case that the causal effect of two context variables $X, Y$ is greater (or lesser) than the sum of the individual effects of $X$ and $Y$. This is called interaction. Composition asks that two variables $X$ and $Y$ that are individually independent from $T$, cannot interact such that $T$ becomes dependent of them. Note that interaction is possible in the case of any of $X$ and $Y$ already being dependent on $T$.
**Interchangeability** Interchangeability asks for the set of multiple solution $\mathcal{M}$ to factorize as a cartesian product of equivalence class of individual variables. This allows the solution space to be quickly discovered by a linear time algorithm.

In general, equivalences can occur at the level of group of variables, and not just individual variables. For instance, let $Z$ be a standard normal random variable. Define $X = \mathbb{I}[Z > 0].Z$ and $Y = \mathbb{I}[Z < 0].Z$. Since $X$ has the value of $Z$ if $Z > 0$ and $Y$ the converse, then $Z = X + Y$. In this situation, for any target, $Z \equiv_T X, Y$. Interchangeability requires that $Z$ cannot be part of a markov boundary of $T$.

It is possible to relax Interchangeability to a looser version by introducing additional complexity to the algorithm. If instead of testing singletons for equivalences, the algorithm also tests pairs of variables, more complex equivalence relations can be integrated. To take into account all possible equivalence relations among group of variables, exponential search is required [SLLA13].

## B Pseudocode of Lag-Pearson-pval routine

**Algorithm 7** Lag-Pearson-pval

**Require:** residuals $R$, candidate variable $X' \in \mathbb{R}^M$, maximum lag $L$
1: $p \leftarrow 1$
2: **for** $1 \le l \le L$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp; $p \leftarrow \min(p, \text{Pearson-corr}(R, X'_{l,M-L+l-1}))$
&nbsp;&nbsp;&nbsp;&nbsp; **return** $p$

**Algorithm 7** Lag-Pearson-pval

**Require:** residuals $R$, candidate variable $X' \in \mathbb{R}^M$, maximum lag $L$

1: $p \leftarrow 1$

2: **for** $1 \le l \le L$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp; $p \leftarrow \min(p, \text{Pearson-corr}(R, X'_{l,M-L+l-1}))$

**return** $p$

## C Pseudocode of equivalence computation of version FE

The equivalence search (Algorithm 8) of the approximate FE version builds surrogate models $g_\Theta$ to predict the residuals instead of the target $T$. Its structure correspond to the definition of Information Equivalence, in that it tests proxies for the two assertions $X \perp Res|S_{new}$, and $S_{new} \perp Res|X$. Establishing equivalence for residuals is not always guaranteed to ascertain equivalence for $T$ (see Theorems), but we empirically find it a good proxy for our synthetic datasets.

**Algorithm 8** FIND-EQUIVALENCES

**Require:** residuals $Res$, selected variable $S_{new}$, remaining variables $R$, maximum lag $L$, equivalence threshold $\delta$

1: set $E \leftarrow \{\}$

2: **for** $X \in R$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp; $p \leftarrow \text{pvalue of Likelihood ratio test between } Res = g_\Theta(X \cup S_{new}) \text{ and } Res = g_\Theta(S_{new})$

4: &nbsp;&nbsp;&nbsp;&nbsp; **if** $p \ge \delta$ **then** &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\triangleright X$ is redundant given $S_{new}$ to predict $Res$

5: &nbsp;&nbsp;&nbsp;&nbsp; $p \leftarrow \text{pvalue of Likelihood ratio test between } Res = g_\Theta(X \cup S_{new}) \text{ and } Res = g_\Theta(X)$

6: &nbsp;&nbsp;&nbsp;&nbsp; **if** $p \ge \delta$ **then** &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\triangleright S_{new}$ is redundant given $X$ to predict $Res$

7: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{E} \leftarrow E \cup \{X\}$

**return** $E$

## D Theorems and proofs

**Theorem D.1 (Forward Soundness).** Assuming composition (Assumption 2.1) and perfect statistical tests for the correlation and termination, the set $S$ returned by the FORWARD and by the FORWARD-EQUIVALENCES procedures is a Markov blanket of $T_t$.
**Forward Soundness.** We reuse the proof of KIAMB [PNBT07]. At each iteration, the association heuristic selects $S_{new}$ such that $S_{new} \perp T_t|S$ whenever such a variable exists. The learning model used will improve in this case and only in this case. When the forward phase ends, for all remaining variables $C \in R$, then $C \perp T_t|S$. By composition, this ensures that $R \perp T_t|S$, hence $S$ is a Markov Blanket of the target $T_t$. □

**Theorem D.2 (Backward Soundness).** Assuming composition (Assumption 2.1) and perfect tests for the correlation, termination, and elimination, the set $S$ returned by the succession of the FORWARD and BACKWARD procedures is a Markov boundary of $T_t$.

**Backward Soundness.** We reuse the proof of KIAMB [PNBT07]. We assume that the set $S$ produced by the forward phase is a Markov Blanket of $T_t$. Assuming that every conditional independence test of $T_t \perp S_{del}|S \setminus \{S_{del}\}$ is correct. At each deletion step, assuming $V \setminus S \perp T_t|S$ and $T_t \perp S_{del}|S \setminus \{S_{del}\}$. In this case, by contraction (A.2), $T_t \perp S_{del} \cup (V \setminus S)|S \setminus \{S_{del}\}$. This guarantees that $S \setminus \{S_{del}\}$ is still a Markov blanket.

By the end of the backward phase, for all $S \in S$, $T_t \perp S|S \setminus \{S\}$. By contradiction, suppose that there is $M \subset S$ a Markov blanket of $T_t$. Let $X \in S \setminus M$, let $Y = S \setminus M \setminus \{X\}$. Since $M$ is a Markov blanket of $T_t$, $X \perp T_t|M$ and $Y \perp T_t|M$. Applying composition, $X, Y \perp T_t|M$, then applying weak union (A.1), $X \perp T_t|M \cup Y$. This contradicts the end condition of the backward phase, so the $S$ has to be minimal for inclusion among Markov blankets. □

**Theorem D.3 (Soundness of FBE).** Assuming composition and perfect tests for the correlation, termination, elimination and equivalences, any $MB_1...MB_n$ in the set $M$ computed by FBE is a Markov blanket of $T_t$.

**FBE Soundness (Thm.D.3).** In the following, we define a useful notation $S_k$ where $S$ is an indexed set from 1 to $n$ and $1 \le k \le n$, as the set of elements of $S$ indexed from 1 to $k$.

First, we remark that it is enough to prove that $MB_n \equiv_T S_n$, where $MB_n$ is one of the produced solutions $MB_1, ..., MB_n \in M$. In that case, since $S_n$ is a Markov blanket, for any variable of the MTS $B \in X$, then $B \perp T_t|S_n$ holds. By Lemma A.3 applied to the equivalence $MB_n \equiv_T X_n$, we have $B \perp T_t|MB_n$. Hence, $MB_n$ is a Markov blanket of $T_t$.

To prove $MB_n \equiv_T S_n$, we proceed by recursion. Suppose that $MB_{k-1} \equiv_T S_{k-1}|S_k...S_n$. We seek to prove that $MB_k \equiv_T S_k|S_{k+1}...S_n$. For simplicity of notation, in the following proof, denote $Z = S_{k+1}...S_n$. By definition of equivalence, we need to prove that

$$ \begin{cases} MB_k \perp T_t | S_k, Z \\ S_k \perp T_t | MB_k, Z \end{cases} $$

The first condition is true since $S_n$ is a Markov Blanket of $T_t$. To prove the second condition, by the property of decomposition, and by the assumption of composition (Assumption 2.1), it is necessary and sufficient to prove that

$$ S_k \perp T_t | MB_k, Z \qquad (3) $$

$$ S_{k-1} \perp T_t | MB_k, Z \qquad (4) $$

By using Lemma A.3 with equivalence $MB_{k-1} \equiv_T S_{k-1}|S_k, Z$ (recursion supposition), equation (3) is equivalent to $S_k \perp T_t|S_n^{-k} \cup \{MB_k\}$ which is the condition to include $MB_k$ in the equivalence of $S_k$ in our algorithm.

By using Lemma A.3 with equivalence $MB_k \equiv_T S_k|S_n^{-k}$ (implied by the test in the algorithm), equation (4) is equivalent to $S_{k-1} \perp T_t|MB_{k-1}, S_k, Z$. This last relation is true due to the recursion supposition.

The initial condition holds since $S_1 \equiv_T MB_1|S_2...S_n$. Hence, $MB_n \equiv_T MB_n$. □

**Corollary D.3.1 (of Thm.D.3).** Assuming composition, assuming that all Markov boundaries have the same size and perfect tests for the correlation, termination, elimination and equivalences, any equivalent solution is a Markov boundary.
**Corollary.** All returned solutions are Markov Blankets (Thm.D.3), and at least one is a Markov Boundary (Thm.D.2). Since by definition of Markov Boundaries, any Markov Blanket is the superset of at least one Markov Boundary, then all returned solutions must be of the same size as the Markov Boundary contained. Hence, all solutions are Markov Boundaries. □

**Completeness (Thm.3.1).** Direction $\implies$: suppose that the equivalence search starting from a selected set $S$ find out all Markov boundaries of $T_t$. From the proof of Thm D.3, we can deduce that for whichever indexing of the reference set $S$, denoting $S_1...S_n$, and for whichever equivalent boundary $M$ produced by the algorithm, there is an ordering $M_1...M_n$ such that

$$ \forall i, M_i \equiv_T S_i | S \setminus S_i \qquad (5) $$

This is the case since the proof of Thm D.3 does not rely on a particular ordering of $S$.

Let us consider the reference boundary $S$ and any other boundary $M$ of $T_t$. By supposition, $M$ is included in the solutions of the algorithm. Let us consider a reference indexing $\mathcal{I}$ of $S$ and $M$ as is provided by the successive equivalent sets of our algorithm. Define any function $f: [|1, n|] \to 0, 1$, resulting in a intermediary set $A$.

There is a permutation $\sigma: [|1, n|] \to [|1, n|]$ such that $f \circ \sigma$ is monotonically decreasing. We reindex $S'_i = S_{\sigma^{-1}(i)}$ and similarly for $M'_i$. Since Eq.(5) holds for any common reindexing of $S$ and $M$, it holds for $S'$ and $M'$. Also, there is an integer $k$ for which the set $M'_k \cup (S' \setminus S'_k)$ is identical to $A$. By Thm.D.3, $M'_k \cup (S' \setminus S'_k)$ is a Markov boundary.

Direction $\Leftarrow$: suppose that for any correctly indexed $M$, any function $f: [|1, n|] \to 0, 1$ induces a set $A$ and $A$ is a Markov boundary of $T_t$. In particular, for any $k \in [|1, n|]$, the function $f_k: i \mapsto 1$ if $k = i$ else 0 induces a Markov boundary $S^{-k} \cup M_k$. Therefore, we deduce that $S_k \perp T_t | S^{-k} \cup M_k$ for all $k$, which is exactly the conditions verified by the algorithm. Hence, each $M_k$ is included in an equivalence set, and $M$ is an equivalent Markov boundary produced by the algorithm. □

**Residual tests (Thm.3.2).** Since $T_t \perp X | Z$, then $T_t, Z \perp X | Z$. Any deterministic transformation of $T_t, Z$ cannot add information so $f(T_t, Z) \perp X | Z$. Using $Res \perp Z$, by contraction (Lemma A.2), $Res \perp X, Z$. Thus (decomposition) $Res \perp X$. □

**Residual tests with composition (Thm.3.2).** Starting from $Res \perp X, Z$ and $Res \perp Z$, by weak union (LemmaA.1), $Res \perp X | Z$. Naturally, $Res, Z \perp X | Z$. Any deterministic transformation of $Res, Z$ cannot add information and $T_t = f(Res, Z)$. So $T_t \perp X | Z$. □

**Residual tests for normal variables (Thm.3.2).** Since the distribution is joint normal, then $Res = T_t - (a_0 + \sum_i a_i Z_i)$ is also a normal variable and $V \cup \{Res\}$ is also joint normal. By lemma A.4, we can verify that this distribution verifies composition.

The residuals must be independent from $Z$ as joint normal variables are independent if and only if they are uncorrelated. Thus, pairing $Res \perp X$ and $Res \perp Z$ with composition (Assump.2.1), $Res \perp X, Z$. □

# E Additional information on experiments

## E.1 Tuning study

We optimize the hyperparameters for each algorithm and forecaster we consider (see Table 3). We use the python library optuna with Grid Search optimization. ChronoEpiloghi thresholds cover several orders of value between $10^{-20}$ and 0.05. As we observed that for GroupLasso regularization parameter $10^{-20}$, all TS where select no matter the dataset, and at 0.1, no TS was selected other than the target past, GroupLasso group regularization parameter ranges 25 values within this range on a logarithmic scale.

We describe the pipeline used for tuning and evaluation in Fig.3. Cross-validation on time series is done in a fixed start window growing along time. The holdout set size is consistently 1000 for synthetic MTS, and 30 percent of the MTS for real dataset. The smallest training window size must
<figure><img src="image_4.png" /><figcaption>Figure 3: General experimental pipeline. We use K-fold cross validation suited to time series, for the tuning of hyperparameters of both TSS algorithms and forecasters.</figcaption></figure>

<table><caption>Table 3: Algorithms hyperparameters</caption><thead><tr><td>Model</td><td>Tuned hyperparameters</td></tr></thead><tbody><tr><td>GroupLasso [NST13]</td><td>group regularization coefficient</td></tr><tr><td>ARDL [PSS01]</td><td>Trend (no trend, constant, linear)</td></tr><tr><td>SVR</td><td>Kernel (rbf, sigmoid), regularization weight</td></tr><tr><td>DeepAR [SFGJ20]</td><td>Layer size, #layers, dropout rate</td></tr><tr><td>TFT [LALP21]</td><td>Layer size, #layers, dropout rate, #attention heads</td></tr></tbody></table>

cover at least 40 percents of the Tuning split, so to obtain meaningful results and retain enough samples.

Over Fig.4, we examine how metrics vary when the thresholds are chosen to maximize the $R^2$ of the forward-backward solution $S$. We observe that the forward and backward threshold depend on the $R^2$ value, hence the proportion useful information variation compared to the total variation of the data. We confirm that different thresholds do not result in vastly different performances.

The standard deviation of the $R^2$ of the top 5 configurations for each MTS on average is 0.001 for FBE and FE compared to 0.02 for GroupLasso. See 5 for the boxplot of the distribution. This indicate a much higher impact of hyperparameter tuning on GroupLasso than ChronoEpilogi, and participates to the difficulty of tuning GroupLasso. Practictioners should expect spending more time tuning regularization-based solution than p-value-based solutions.

## E.2 Statistical significance of results on real datasets

The statistical analysis was conducted for 3 populations with 50 paired samples. The family-wise significance level of the tests is $\alpha=0.050$. We rejected the null hypothesis that the population is normal for the populations ChronoEpilogi ($p=0.000$), GroupLasso ($p=0.000$), and NoSelection ($p=0.000$). Therefore, we assume that not all populations are normal. Because we have more than two populations and some of them are not normal, we use the non-parametric Friedman test as omnibus test to determine if there are any significant differences between the median values of the populations.
<figure><img src="image_5.png" /><figcaption>Figure 4: Optimal hyperparameters choice influence on performance metrics for version FBE. Lines correspond to the forward phase threshold and backward phase threshold. Columns correspond to metrics (R2, f1score on the detection of irreplaceable variables, f1-score on the detection of replaceable variables.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Standard deviation of test R2 of the top 5 configurations, over the 270 MTS of the synthetic dataset, for version FBE, FE and GroupLasso.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Evaluation of FBE (solid lines) and FE (dashed lines) over the synthetic dataset. Column 1: computing time of each method. Columns 2-5: Value of each causal metric for varying data characteristics. For columns 1-4, we report the performance of maximal sample size (8000).</figcaption></figure>
<table><caption>Table 4: FBE predictive performance on Traffic and METR-LA. We observe that version FBE increases predictive performance compared to FE while diminishing selected size, at the cost of higher computational time</caption><thead><tr><th>dataset</th><th>R2</th><th>rmse</th><th>mape</th><th>size</th><th>time</th><th>#MB</th></tr></thead><tbody><tr><td>METR-LA</td><td>0.922622</td><td>0.316716</td><td>1.135145</td><td>2.1</td><td>1629.063206</td><td>2.8</td></tr><tr><td>Traffic</td><td>0.817231</td><td>0.413731</td><td>30.916130</td><td>5.3</td><td>14777.317261</td><td>36288.3</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Multiple solutions predictive performance over METR-LA and Traffic. When more than 20 solutions are found, we randomly select 20 of them for evaluation. We observe that both the number of targets for which there are multiple solutions decreases, with a diminution of the number of solutions overall.</figcaption></figure>

We use the post-hoc Nemenyi test to infer which differences are significant. We report the median (MD), the median absolute deviation (MAD) and the mean rank (MR) among all populations over the samples. Differences between populations are significant, if the difference of the mean rank is greater than the critical distance CD=0.469 of the Nemenyi test. We reject the null hypothesis (p=0.000) of the Friedman test that there is no difference in the central tendency of the populations ChronoEpilogi (MD=0.289+-0.112, MAD=0.105, MR=2.580), GroupLasso (MD=0.287+-0.106, MAD=0.106, MR=2.320), and NoSelection (MD=0.442+-0.310, MAD=0.257, MR=1.100). Therefore, we assume that there is a statistically significant difference between the median values of the populations. Based on the post-hoc Nemenyi test, we assume that there are no significant differences within the following groups: ChronoEpilogi and GroupLasso. All other differences are significant.

## E.3 ChronoEpilogi-FBE on Traffic and METR-LA

For comparison purpose, we ran the full version of the algorithm on two of the real datasets. The performance is reported in Table 4. The number of target for which multiple solutions are found is reduced to 2 in METR-LA and 6 on Traffic (Figure 8). We observe that overall, results from synthetic datasets are verified: FBE decreases selected size and diminishes the total number of solutions, with higher computational time, resulting in an increase in predictive performance. We confirm that version FE tend to select false positives, but keep the constatation that multiple solution exist in real datasets.

## E.4 SVR models outperforming TFT models on Solar

A surprising experimental result was that at the end of parameter tuning, SVRModels were outperforming TFT models during cross validation. As we suspected that cross validation could influence negatively the TFT performance (earlier split having less data, and Transformer-based models notoriously requiring consequent datasets), we tried a different pipeline where tuning and testing was done on a single training/validation/test split with 70,20,10 repartition. The result are in Table 5.
<table><caption>Table 5: Forecasting performance on Solar with modified tuning protocol. TFT models win on most targets. The higher computation time is due to running the experiment on a 16 cores machine (instead of 36) due to server availability.</caption><thead><tr><th>TSS</th><th>R<sup>2</sup></th><th>rmse</th><th>mape</th><th>size</th><th>time F/E</th><th>#model</th><th>#MB</th></tr></thead><tbody><tr><td>ChronoEpilogi-FE</td><td>0.990</td><td>0.087</td><td>0.730</td><td>5.1</td><td>3/473</td><td>10/0/0</td><td>1</td></tr><tr><td>GroupLasso</td><td>0.989</td><td>0.089</td><td>0.600</td><td>4.7</td><td>296</td><td>8/0/2</td><td>NA</td></tr><tr><td>NoSelection</td><td>0.985</td><td>0.100</td><td>0.612</td><td>138.0</td><td>NA</td><td>10/0/0</td><td>NA</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 8: PEMS-BAY dataset visualisation. The sensors (dots) are overlaid with the highways and primary road networks of San Jose. Red sensors are either redundant or irrelevant. The green square figures the target variable. Blue sensors are irreplaceable and the two orange ones are replaceable sensors in a same equivalence class. We observe, confirming the observations of the original paper [LYSL18], that important predictors are close spatially.</figcaption></figure>

We note that no matter the predictive model, ChronoEpilogi does not return multiple solutions on this dataset. Additionally, the results using the modified pipeline do not change the method ranking of ChronoEpilogi, GroupLasso and NoSelection.

## E.5 Feature Importance under the existence of multiple Markov Boundaries

It is known that Shapley values currently suffer from inclusion of unrealistic data instances when features are correlated, the explanations even of linear models can be misleading [AJL21]. In this section, we investigate how SHAP explanations of regression models might misrepresent the role of variables belonging to some but not all Markov boundaries of the modeled target.

To assess the impact of replaceable variables on SHAP explanations, we tune and train Support Vector Regressors (SVR) and XGBoost Regressors (XGB) for each of the 90 synthetic MTS with 100 covariates (denoted as *full MTS*), using a window of 10 timesteps. We then compute the SHAP values of random holdout samples using the KernelSHAP approximation [LL17], that we sum along the time axis to obtain the overall importance of variables (per additivity property of SHAP values [KVSF20]). We consider the average of those values as global explanations [BSC$^{+}$21]. We obtain average SHAP values $s_i^1...s_{100}^i$ corresponding to the global importance of each variable in the MTS number *i*.
To evaluate the total importance for each set of equivalent variables, without the disturbance brought by the presence of replaceable variables, we apply the same average SHAP computations as previously to the 90 MTS, where we remove all but one replaceable variable per equivalent set (denoted as reduced MTS), to compare with SHAP explanation of unique Markov Boundary processes. For instance, if in MTS *i* variables {1, 4, 16, 78} are equivalent, we remove {4, 16, 78} and compute the reduced importance *sr*<sub>1</sub><sup>*i*</sup> of the representative 1. Note that variable 1 is irreplaceable in the reduced MTS.

For each equivalent set {$j_1, ..., j_c$} in MTS *i*, we compute the contribution of each variable *j* ∈ {$j_1, ..., j_c$} to the total absolute importance of the set, as $\frac{|s_j^i|}{\sum_k |s_{j_k}^i|}$. We plot the contribution of the most important variable, second most important, etc (Fig. 9). We choose to use absolute importance values as averages are not easily comparable with relative importance. We observe that **importance is not equally shared among equivalent variables**. Models tend to select one or a few variables as important predictors. Additionally, **the importance of the top variable relative to the other diminishes as the number of equivalent variables grows.**

<figure><img src="image_10.png" /><figcaption>Figure 9: Average proportion of importance of top1, top2, ... variables in an equivalent set among the total importance of the equivalent set for a) SVR model b) XGB model. Models tend to focus predominantly on one variable, but with diminishing importance as the number of equivalent variables grows.</figcaption></figure>

<table><caption>**Table 6: The most important variable of an equivalent set is highly unstable to data resampling.**
For each synthetic MTS with 100 variables, the stability (against data variability) of the most important (top1) variable in each equivalent set is measured with the stability [NSB18] metric. Bootstraps are computed by subsampling contiguous intervals of timestamps, for 21 evaluations. The metric adjusts for equivalent set size. We observe that the top1 variable of an equivalent set is highly unstable. Relevance estimation based on feature stability might fail when replaceable variables exist.</caption><thead><tr><td>Model</td><td>Average stability (higher is more stable)</td><td>Standard deviation of stability</td><td>Theoretical min (most unstable): -1/(M-1) with M number of bootstraps</td><td>Theoretical max (most stable)</td></tr></thead><tbody><tr><td>XGB</td><td>0.052</td><td>0.110</td><td>-0.05</td><td>1</td></tr><tr><td>SVR</td><td>0.020</td><td>0.089</td><td>-0.05</td><td>1</td></tr></tbody></table>

## E.6 Hardware and Codebase

We make available our anonymized codebase at https://github.com/ev07/ChronoEpilogi.

We run our experiments on servers running Ubuntu 22.04.4 LTS, with 36 cores, 1TB RAM, GPU Quadro RTX 8000 with driver version 550.54.14 and CUDA version 12.4.

We include in the repository a requirements file `requirements.txt` listing the necessary dependencies.
<table><caption>Table 7: Presence of replaceable variable increases the instability of top important variables.
For each synthetic MTS with 100 variables, all variables but one in each equivalent sets are removed to make a reduced MTS with a unique Markov Boundary. XGB models are computed and SHAP explanations produced. We compute the stability [NSB18] of the $|MB(T)|$ most important variables in the entire explanation, for the original MTS (with multiple MB(T)) and the reduced MTS (with a unique MB(T)). Note that the metric adjusts for total number of variables in the dataset [NSB18]. We observe that explanations are significantly more unstable when replaceable variables are present.</caption><thead><tr><td>Model</td><td>Average Stability</td><td>Standard deviation of Stability</td><td>P-value of a paired wilcoxon signed rank test (lower indicates that XGB (reduced) is more stable)</td></tr></thead><tbody><tr><td>XGB (original)</td><td>0.380</td><td>0.140</td><td rowspan="2">4.2 e-13</td></tr><tr><td>XGB (reduced)</td><td>0.516</td><td>0.195</td></tr></tbody></table>