ChronoEpilogi: Scalable Time-Series Variable Selection with Multiple Solutions

ChronoEpilogi: Scalable Time-Series Variable Selection with Multiple Solutions Etienne Vareille1Michele Linardi1Ioannis Tsamardinos2Vassilis Christophides1 1ETIS UMR-8051 Laboratory, CY Cergy Paris Universite, ENSEA, CNRS 2Computer 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 38th Conference on Neural Information Processing Systems (NeurIPS 2024). 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 MB(T) belong. This is because we can construct a new 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 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 MB(T) (referred to as the multiple TVS (MTVS) problem) is crucial. While finding any one 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 MB(T) that contains them. In essence, multiple 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 crossvalidation 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 MB(T) :indispensable, variables that belong in all MB(T) ,replaceable, variables that belong in some 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 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 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 12 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 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 1Chronos greek and Epilogi greek 2https://github.com/ev07/ChronoEpilogi 2 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 X={X1, ..., XN} , where each Xi∈RM ( 1≤i≤N ) are regularly sampled observations of a time series variable. We denote with Xi t∈R the tth occurrence of Xi , and by Xi t−L,t−1∈RL the last L occurrences of Xi before the instant t . Note that, in the following part, bold letters (e.g., X ) refer to sets of TS variables. The TS forecasting task we consider is to predict the values of a target TS variable Tt ( T = Xi for some i ), where L≤t≤M , with a model f parameterized by Θ that uses up until L past timestamps (for each t) of the multivariate TS X. Hence, Tt== fΘ(Xt−L,t−1). Conditional independence of two TS variables Xt and Yt′ given a third one Zt′′ , written as Xt′ |= Yt′|Zt′′ , is defined as P(Xt, Yt′|Zt′′ ) = P(Xt|Zt′′ ).P(Yt′|Zt′′ ) . In our work, we assume stationarity, namely the conditional distribution of Tt as a function of L previous lags of 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 Tt and variables with timestamps ranging between t−Land t−1. Due to stationarity, we can remove the index and denote Xt−L,t−1as X. Definition 2.1 (Information Equivalence (IEQ) ).Two TS variables Xa and Yb are information equivalent (shorthand: equivalent) with respect to a target Tt given conditioning TS variable set Z , noted Xa≡TYb|Z , iff they are made independent of Tt by conditioning on the other and on Z itself. Formally: Xa≡TYb|Z⇐⇒ Xa |= Tt|Yb,Zand Yb |= Tt|Xa,Z . We write shortly Xa≡TYb when the conditioning set is Z={}. Definition 2.2 (Markov Boundary MB(Tt) ).Given a set of TS variables X and a target variable T in this set, M(Tt) is a Markov blanket of Tt iff: Tt |= (V\M)|M . The set M(Tt) is a Markov boundary of Ttiff ∀M′⫋M, Tt⊥⊥ (V\M′)|M′. Definition 2.3 (Multiple Time-series variable selection (MTVS)).Let MB(Tt) a reference MB of a target TtTS variable. The solution of the MTVS problem consists into finding all MBs, denoted by M . All MB(Tt) are information equivalent. For all MBi(Tt)∈ M the forecasting models Tt=fΘ(MBi(Tt)) are equally-performing according to a metric (e.g., R2). 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 ∈R be the inflows of the two confluent rivers A and B . The dam is controlled such that if Zt≤zth , it does not produce energy. Otherwise, it diverts a flow zth to power production. The power production flow mixes with river A first for total flow Yt=f(Zt) + Xt , then with the rest of river B for total flow Tt=Xt+Zt , with f(Zt)∈ {0, zth} . 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 Tt [ BT21 ]. A variable Xa is said irreplaceable iff it is part of all information equivalent Markov boundaries of Tt , i.e., Xa∈ MB(Tt)∀MB(Tt)∈ M . A variable Xa is said replaceable iff it is part of at least one MB(Tt) but not indispensible, i.e., ∃MB(Tt),MB′(Tt)∈ M, Xa∈MB(Tt)∧Xa/∈MB′(Tt) . A variable Xa is redundant iff ∃Z⊆X\ {Xa}, Xa⊥⊥ Tt|Z and ∀MBi(Tt)∈ M, Xa/∈MBi(Tt) while it is called irrelevant iff ∀Z⊆X\ {Xa}, Xa |= Tt|Z. If Causal Markov Condition and Faithfulness [ CMJSV+23 ] hold, each target Tt 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 3 MB(Tt) might not be unique. Jointly Gaussian random variables (RV) with a singular covariance matrix or with deterministic relations do no respect faithfulness 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,Tand conditioning set Z: X |= T|Zand Y |= T|Z=⇒X∪Y |= T|Z(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 Yi=fi(Xσ(i)) , with fi invertible and σ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 |= T|Z and Y |= T|Z hold, the union X∪Y contains only RVs that are pairwise independent from each variable in Tgiven Z. This implies X∪Y |= 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 Tt , M B(Tt) and M B′(Tt) are interchangeable iff: ∀M B(Tt),M B′(Tt)∈ M,∀X∈M B(Tt),∃Y∈M B′(Tt),({X}∪M B′(Tt)\{Y})∈ 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 M= ({Z, X},{Z, Y }) satisfies Interchangeability, as M={Z}×{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 [α, γ, δ] 1: set S←FORWARD(X, T, L,α) 2: S←BACKWARD(X, T, L,S,γ) 3: set M ← EQUIV(X, T, L,S,δ) 4: return M▷set of eq. Markov bound. Algorithm 2 ChronoEpilogi-FE Require: TS X , target TS variable T , max lag L, threshold params [α, δ] 1: M ← FORW-EQUIV(X, T, L,α) 2: return M▷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 3Conditions 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 phenomenons [YWD+17]. 4 Algorithm 3 FORWARD Require: TS X∈RN×M, target T, max lag L, stopping threshold α 1: set S← {T1,M−1} 2: model m←fit and save fΘsuch that TL+1 =fΘ(S1,L), ..., TM=fΘ(SM−L,M−1) 3: repeat 4: R←m.residuals ∈RM−L 5: Snew ←arg minX′∈X\SLag-Pearson-pval(R, X′, L)▷Select Snew s.t. Snew ⊥⊥ R 6: S←S∪Snew 7: model m′←m 8: model m←fit and save fΘon S 9: until pvalue of Likelihood ratio test of (H0: m′=m)≥α ▷ Verify Snew ⊥⊥ Tt+1|S\Snew 10: return S\Snew Algorithm 4 BACKWARD Require: TS X, target T, time pred. t, max lag L, sel. variables S, stopping threshold γ 1: repeat 2: for Sdel ∈Sdo 3: if pvalue of Likelihood ratio test between fΘ(S\ {Sdel})and fΘ(S)≥γthen 4: S←S\ {Sdel} 5: until no change in S 6: return 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 |= T|X, Z , removing U from the selected set. Finally, the equivalence phase would test Y≡TX|Zand produce the solution space 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 ( X ) as predictors, with Tt=a0+Pp i=1 ai.Tt−i+Pk i=1 Pq j=1 bj.Xi t−j+ϵt . The model includes autoregressive terms (aj) and other explanatory variable terms ( bj ). It is generally required that ϵ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 X . At each iteration, a selected variable ( Snew ) 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). Algo. 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 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 (Algortihm 5) tests the equivalence between each selected TS variable S (line 2), and any non-selected variable Sd (line 4). In line 5we 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 Q(line 6). 5 Algorithm 5 EQUIV Require: TS X, target T, time pred. t, max lag L, sel. variables S, equivalence threshold δ 1: dictionary Q← {:} 2: for S∈Sdo 3: Q[S]← {} 4: for Sd∈(X\S)do 5: if pvalue of LR test between fΘ({Sd} ∪ S\ {S})and fΘ({Sd} ∪ S)≥δthen 6: Q[S]←Q[S]∪ {Sd} 7: M←{S1} ∪ Q[S1]×... × {Sn} ∪ Q[Sn] 8: return M Algorithm 6 FORW_EQUIV Require: MTS X, target T, maximal lag L, stopping threshold α, equivalence threshold δ 1: set S← {T1,M−1} 2: set R←X\S 3: dictionary Q← {:} 4: model m←fit and save fΘsuch that TL+1 =fΘ(S1,L), ..., TM=fΘ(SM−L,M−1) 5: repeat 6: Res ←m.residuals ∈RM−L 7: Snew ←arg min{X′∈R}Lag-Pearson-pval(Res, X′, L)▷Select Snew s.t. Snew ⊥⊥ Res 8: S←S∪ {Snew} 9: R←R\ {Snew} 10: Q[Snew]←FIND-EQUIVALENCES(Res, Snew,R, δ)▷Test C≡Res Snew for C∈R 11: R←R\Q[Snew] 12: m′←m 13: model m←fit and save fΘon S 14: until pvalue of Likelihood ratio test of (H0: m′=m)≥α ▷ Verify Snew ⊥⊥ Tt+1|S\Snew 15: return S\Snew , Q Once the equivalent Markov Boundaries 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 Algortihm 5. In Algorithm 6we 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 8in 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Θ(.) . In this sense, if a non-selected variable X is equivalent to the tested variable, Snew is added to the equivalent set (line 7). Complexity analysis Given a TS X∈RN×M , and denoted by S the selected TS variables in the forward phase, both ChronoEpilogi versions ( FBE and FE ) run O(N|S|) conditional tests. In the FORWARD phase, the sub-routine Lag-Pearson-pval (Algorithm 7) runs in O(L n log(n)) if Fast Fourier Transform is adopted [MZZ+22], where Lis 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 O(|S|3L3M) . Overall, the forward, backward, equivalent phases take O(|S|N L M +|S|4L3M),O(|S|5L3M),O(|S|4L3N 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 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 M is a Markov blanket. In other words, our algorithm does not return false positives. We 4Alternatively, algorithms might rely on an oracle producing variable sets with given properties [SLLA13]. 6 then examine the conditions under which FBE is complete, i.e. all Markov boundaries are solutions in 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 Mof a target Ttand 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 5of 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 |= 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 |= Z , then: T |= X|Z=⇒R |= X . Additionally, if for all candidate variable X , R |= X=⇒R |= X, Z, then R |= X=⇒T |= 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 R2 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 R2 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. 7 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). 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 R2 , 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 = 2T P 2T P +F P +F N . 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 si 1...si 100 corresponding to the global importance of each variable in the MTS number i . The total importance of each equivalent variable set S becomes stot =Pi∈S|si| . Then, we compute the contribution percentage of the most important variables (until top-14) in the equivalent set |si/stot| . 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 sri 1 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 R2 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 R2in Electricity, and spread over 0.02 around the median for Traffic (Fig.2c). Fig.1illustrates 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). 8 0 1000 Number of TS (a) 0 500 1000 Computing time (s) 101102103 Number of TS (b) 0.7 0.8 0.9 1.0 Causal metrics metric iF1 rF1 (a) computing time (b) f1 score for replaceable (rF1) and irreplaceable (iF1) variables of FBE (solid lines) and FE (dashed lines) over the synthetic dataset. FE has up to 70% speedup while losing at most 0.1 in rF1. pems-bay metr-la traffic electricity 0.4 0.6 0.8 1.0 R2 2e+00 6e+01 6e+05 1e+04 3e+00 4e+01 3e+03 2e+05 1e+01 1e+04 4e+02 2e+19 2e+01 3e+00 2e+09 7e+12 1e+15 2e+07 2e+09 (c) Boxplot of the R2 of the 20 sampled MB for each target where ChronoEpilogi finds multiple MBs. Red dots figure the performance of GroupLasso. The number of MB detected is reported for each target. Figure 2: Performance of ChronoEpilogi versions on (a)(b) synthetic MTS and (c) real MTS. 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. Time(s) R2MB size Number of MB causal f1-score FBE 375±772 0.373±0.236 5.69±3.33 8.87e+05±3.7e+06 0.99±0.064 FE 118±178 0.373±0.236 5.74±3.31 6.87e+07±7.2e+08 0.94±0.188 GL 120±149 0.337±0.239 7.50±10.2 NA 0.57±0.349 Note that computing one solution with FE forward phase is up to two orders of magnitude faster than GL. Table 2reports 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 |M B(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.ns) 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 9 Interchangeability Interchangeability asks for the set of multiple solution 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=I[Z > 0].Z and Y=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≡TX, 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 7details the Pearson correlation based test. We note that in this version, a p-value is obtained for each of the lags, the minimum of which is returned. This procedure returns p-values that are likely lower than the actual p-value, due to using the minimum as an estimator. Hence, testing for previous and current model statistical difference is kept separate from the correlation routine, ensuring that Forward phase termination relies on a single p-value. Algorithm 7 Lag-Pearson-pval Require: residuals R, candidate variable X′∈RM, maximum lag L 1: p←1 2: for 1≤l≤Ldo 3: p←min(p, 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Θ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 |= Res|Snew , and Snew |= 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 Snew , remaining variables R , maximum lag L , equivalence threshold δ 1: set E← {} 2: for X∈Rdo 3: p←pvalue of Likelihood ratio test between Res =gΘ(X∪Snew)and Res =gΘ(Snew) 4: if p≥δthen ▷ X is redundant given Snew to predict Res 5: p←pvalue of Likelihood ratio test between Res =gΘ(X∪Snew)and Res =gΘ(X) 6: if p≥δthen ▷ Snew is redundant given Xto predict Res 7: E←E∪ {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 FORWARDEQUIVALENCES procedures is a Markov blanket of Tt. 16 Forward Soundness. We reuse the proof of KIAMB [ PNBT07 ]. At each iteration, the association heuristic selects Snew such that Snew ⊥⊥ Tt|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∈R , then C |= Tt|S . By composition, this ensures that R |= Tt|S , hence S is a Markov Blanket of the target Tt. 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 Tt. 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 Tt . Assuming that every conditional independence test of Tt |= Sdel|S\{Sdel} is correct. At each deletion step, assuming V\S |= Tt|S and Tt |= Sdel|S\{Sdel} . In this case, by contraction (A.2), Tt |= Sdel ∪(V\S)|S\ {Sdel}. This guarantees that S\ {Sdel} is still a Markov blanket. By the end of the backward phase, for all S∈S , Tt⊥⊥ S|S\ {S} . By contradiction, suppose that there is M⊂S a Markov blanket of Tt . Let X∈S\M , let Y=S\M\ {X} . Since M is a Markov blanket of Tt , X |= Tt|M and Y |= Tt|M . Applying composition, X, Y |= Tt|M , then applying weak union (A.1), X |= Tt|M∪Y . This contradicts the end condition of the backward phase, so the Shas 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 MB1...MBn in the set M computed by FBE is a Markov blanket of Tt. FBE Soundness (Thm.D.3). In the following, we define a useful notation Sk where S is an indexed set from 1 to nand 1≤k≤n, as the set of elements of Sindexed from 1to k. First, we remark that it is enough to prove that M Bn≡TSn , where M Bn is one of the produced solutions MB1, ..., MBn∈ M . In that case, since Sn is a Markov blanket, for any variable of the MTS B∈X , then B |= Tt|Sn holds. By Lemma A.3 applied to the equivalence M Bn≡TXn , we have B |= Tt|M Bn. Hence, M Bnis a Markov blanket of Tt. To prove M Bn≡TSn , we proceed by recursion. Suppose that M Bk−1≡TSk−1|Sk...Sn . We seek to prove that M Bk≡TSk|Sk+1...Sn . For simplicity of notation, in the following proof, denote Z=Sk+1...Sn. By definition of equivalence, we need to prove: M Bk |= Tt|Sk,Z Sk |= Tt|M Bk,Z The first condition is true since Sn is a Markov Blanket of Tt . 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 Sk |= Tt|M Bk,Z(3) Sk−1 |= Tt|M Bk,Z(4) By using Lemma A.3 with equivalence M Bk−1≡TSk−1|Sk,Z (recursion supposition), equation (3) is equivalent to Sk |= Tt|Sn−k∪{MBk} which is the condition to include MBk in the equivalence of Skin our algorithm. By using Lemma A.3 with equivalence MBk≡TSk|Sn−k (implied by the test in the algorithm), equation (4) is equivalent to Sk−1 |= Tt|M Bk−1, Sk,Z . This last relation is true due to the recursion supposition. The initial condition holds since S1≡TMB1|S2...Sn. Hence, M Bn≡TM Bn. 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. 17 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 =⇒ : suppose that the equivalence search starting from a selected set S find out all Markov boundaries of Tt . From the proof of Thm D.3, we can deduce that for whichever indexing of the reference set S , denoting S1....Sn , and for whichever equivalent boundary Mproduced by the algorithm, there is an ordering M1...Mnsuch that ∀i, Mi≡TSi|S\Si(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 Tt . By supposition, M is included in the solutions of the algorithm. Let us consider a reference indexing I of S and M as is provided by the successive equivalent sets of our algorithm. Define any function f: [|1, n|]→0,1 , resulting in a intermediary set A. There is a permutation σ: [|1, n|]→[|1, n|] such that f◦σ is monotonically decreasing. We reindex S′ i=Sσ−1(i)and similarily for M′ i. Since Eq.(5) holds for any common reindexing of Sand M, it holds for S′ and M′ . Also, there is an integer k for which the set M′ k∪(S′\S′ k) is identical to A . By Thm.D.3,M′ k∪(S′\S′ k)is a Markov boundary. Direction ⇐= : suppose that for any correctly indexed M , any function f: [|1, n|]→0,1 induces a set A and A is a Markov boundary of Tt . In particular, for any k∈[|1, n|] , the function fk:i7→ 1if k=ielse 0 induces a Markov boundary S−k∪Mk . Therefore, we deduce that Sk |= Tt|S−k∪Mk for all k , which is exactly the conditions verified by the algorithm. Hence, each Mk is included in an equivalence set, and M is an equivalent Markov boundary produced by the algorithm. Residual tests (Thm.3.2). Since Tt |= X|Z , then Tt,Z |= X|Z . Any deterministic transformation of Tt,Z cannot add information so f(Tt,Z) |= X|Z . Using Res |= Z , by contraction (Lemma A.2), Res |= X, Z. Thus (decomposition) Res |= X. Residual tests with composition (Thm.3.2). Starting from Res |= X, Z and Res |= Z , by weak union (LemmaA.1), Res |= X|Z . Naturally, Res, Z |= X|Z . Any deterministic transformation of Res, Z cannot add information and Tt=f(Res, Z). So Tt |= X|Z. Residual tests for normal variables (Thm.3.2). Since the distribution is joint normal, then Res = Tt−(a0+PiaiZi)is also a normal variable and V∪ {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 |= X and Res |= Z with composition (Assump.2.1), Res |= 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. ChronoEpilogi 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 18 For each fold: - apply TSS, get solution S - fit Forecaster (MSE) on S - fit Forecaster on validation (R2) over selected series S Time Series Selection HP Optimization Evaluation Compute R2 on the Holdout set TSS precision, recall (if synthetic data) TSS size, time spent Fold 1 Fold K Train Validation TIME Holdout Tuning Hyperparam Sampler Select parameters with best average R2 TSS variant Forecaster Train best configuration on Tuning (TSS algorithm and Forecaster) Find the equivalent markov blankets EVALUATION PIPELINE 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. Table 3: Algorithms hyperparameters Model Tuned hyperparameters GroupLasso [NST13] group regularization coefficient ARDL [PSS01] Trend (no trend, constant, linear) SVR Kernel (rbf, sigmoid), regularization weight DeepAR [SFGJ20] Layer size, #layers, dropout rate TFT [LALP21] Layer size, #layers, dropout rate, #attention heads 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 R2 of the forward-backward solution S . We observe that the forward and backward threshold depend on the R2 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 R2 of the top 5 configurations for each MTS on average is 0.001 for FBE and FE compared to 0.02 for GroupLasso. See 5for 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. 19 Figure 4: Optimal hyperparameters choice influence on performance metrics for version FBE. Lines corresponds 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. ChronoEpilogi-FBE ChronoEpilogi-FE GroupLasso 0.00 0.05 0.10 0.15 0.20 0.25 Std. dev. of the R2 of the top 5 configurations for each MTS 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. 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). 20 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 dataset R2 rmse mape size time #MB METR-LA 0.922622 0.316716 1.135145 2.1 1629.063206 2.8 Traffic 0.817231 0.413731 30.916130 5.3 14777.317261 36288.3 metr-la traffic 0.7 0.8 0.9 1.0 R2 Figure 7: Multiple solutions predicitive 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. 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 mutliple 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. 21 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. TSS R2rmse mape size time F/E #model #MB ChronoEpilogi-FE 0.990 0.087 0.730 5.1 3/473 10/0/0 1 GroupLasso 0.989 0.089 0.600 4.7 296 8/0/2 NA NoSelection 0.985 0.100 0.612 138.0 NA 10/0/0 NA 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. 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 si 1...si 100 corresponding to the global importance of each variable in the MTS number i. 22 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 sri 1 of the representative 1 . Note that variable 1 is irreplaceable in the reduced MTS. For each equivalent set {j1, ..., jc} in MTS i , we compute the contribution of each variable j∈ {j1, ..., jc} to the total absolute importance of the set, as |si j| Pk|si jk| . 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. 2 4 6 8 10 12 14 Number of variables in the equivalent set 0.0 0.2 0.4 0.6 0.8 1.0 Contribution |Si/Stot| top1 top2 top3 top4 top5 top6 top7 top8 top9 top10 top11 top12 top13 top14 2 4 6 8 10 12 14 Number of variables in the equivalent set 0.0 0.2 0.4 0.6 0.8 1.0 Contribution |Si/Stot| top1 top2 top3 top4 top5 top6 top7 top8 top9 top10 top11 top12 top13 top14 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. 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. Model Average stability (higher is more stable) Standard deviation of stability Theoretical min (most unstable): -1/(M-1) with M number of bootstraps Theoretical max (most stable) XGB 0.052 0.110 -0.05 1 SVR 0.020 0.089 -0.05 1 E.6 Hardware and Codebase We make available our anonymized codebase at https://github.com/ev07/ChronoEpilogi. We run our experiments on servers runing 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. 23 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. Model Average Stability Standard deviation of Stability P-value of a paired wilcoxon signed rank test (lower indicates that XGB (reduced) is more stable) XGB (original) 0.380 0.140 4.2 e-13 XGB (reduced) 0.516 0.195 NeurIPS Paper Checklist 1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? Answer: [Yes] Justification: At the beginning of Section 5, we added a detailed list of the claims made in the abstract and introduced in the last part of the introduction. Here, we support each claim by providing all the empirical and statistical evidences we have observed. Guidelines: • The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 2. Limitations Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: Limitations related to the scope of the experimental study are discussed in Conclusions Guidelines: • The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper. • The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. 24 • The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren’t acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 3. Theory Assumptions and Proofs Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: In Section 2, we state the main assumptions in our algorithm, namely Composition and Interchangeability, under which we prove soundness and completeness of our algorithm in Section D(Appendix). Guidelines: • The answer NA means that the paper does not include theoretical results. • All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems. • The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 4. Experimental Result Reproducibility Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: ChronoEpilogi is detailed in Section 3, where we provide the pseudocode of all core functions along with the critical details that could help a third-party implementation of our algorithm. Reproducibility of our experimental results is facilitated by the description of the executed pipelines as well as by the availability of the code and data used. Guidelines: • The answer NA means that the paper does not include experiments. • If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same 25