ChronoEpilogi: Scalable Time Series Selection with Multiple Solutions
Authors: Etienne Vareille
Michele Linardi
Ioannis Tsamardinos
Vassilis Christophides
Venue: NeurIPS 2024 poster (2024)
Keywords: Multivariate Time Series Causal Discovery
Forecasting
Explanations
Multiple Markov Boundaries
Primary area: interpretability_and_explainability

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 GroupLasso performance, with the added benefit of providing multiple solutions.

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 Commuterr 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 better understanding of the data-generating mechanism; it is often the first step in causal modeling. While identifying a single solution to the feature selec- tion problem suffices for forecasting purposes, identifying all such minimal-size, optimally predictive subsets is necessary for knowledge discovery and important avoid misleading 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 Interchangeabinng 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 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, is on par with Group Lasso performance, with the added benefit of providing multiple solutions.
Introduction
In analysis tasks involving TS with hundreds or even tens thousands of variables (c.g.. in manufac- turing, 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- seattaggnnrtinggggg 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 small-size MB leads to 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- getta-generatig 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 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 (NcurIPS 2024).

---

(T) requires further analysis with causal discovery algorithms, which is outside the scope of our work (readers are referred to [ADG22] for recent survey).
Both in theory and in practice, there may be numerous Markov Boundaries Tf This is, for example, the case if series subset (X, Y, } is I-to-1 deterministically related with another subset (4, 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 (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 (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 2TCP+ 1T21, SLLA13 13]. For time-series data there have not been any algorithms that return multiple solutions to the TVS problem [YLL221, SLL 15].
Identifying all subsets (T) (referred to as the multiple TVS (MTVS) problem) is crucial. While finding any one MB( suffices for forecasting, this 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( 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, 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 [AJL2I].
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. their presence in multiple (T): indispensable, variables that belong in all (T',), replaceable, variables that belong in some (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 (T), and irrelevant, variables that are completely uninformative. We describe the property of Compositionallty 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 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- non-amamettric assumptions. ChronoEpilogi extends previous variable selection algorithms designed for cohort data [BT21] ]] to 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-performe that ChronoEpilogi has near perfect average causal fl-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 fl score for causal discovery (Claim /). Furthermore, ChronoEpilogi variants are on par with Group Lasso performance [SFHT13] (Claim 4) aaruubly 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
forecasting model. Experiments on synthetic data demonstrate
I CChronos greek and Epilogi greek htttps:// //ggt/u///ttttt /oroo//po//r
2

---

available ChronoEpilogi outperforms GroupLasso in terms of causal fs-se (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 = [x'..
where each X E RM (1 < i S N) are regularly sampled observations of time series variable. We denote with X E the thh occurrence of X and by X E R² the last occurrences of X before the instant t. Note that, in the following part, bold letters (e.g.. X) refer to sets TS variables. The TS forecasting task we consider is to predict the values of target TS variable T; (I'= = X for some i), where S S M, with model ff parameterized by θ that uses up until L past timestamps (for each t) of the multivariate TS X. Hence, T₁ = fe(X.-L,t-1)
Conditional independence of two TS variables X, and Yr given third one Zer, written as X |x| Zer,, is defined as P(X,, Yz = P(X, Zzr)) Ze"). In our work, we assume stationarity, namely the conditional distribution of Tr as a function of L previous lags of X does not depend on L. Stationarity implies temporal consistency, where the conditional independence relations do not rely on In this respect, we consider the principle of temporal precedence according to which causes (i.c., independent variables) occur before an effect (i.e., outcome). We exclusively test conditional independence relations between target T3 and variables with timestamps ranging between t - L and t - 1. Due to stationarity, we can remove the index and denote Xt-L,t-1 as X.
Definition 2.1 (Information Equivalence (IEQ) ). Two TS variables X, and Y, are information equivalent (shorthand: equivalent) with respect to target Tú given conditioning TS variable set Z, noted X. F Yo Z. iff they are made independent of T, by conditioning on the other and on Z itself. Formally: X. =T Yo Z 2 IL Yo, Z and Y, | Z. We write shortly X. =T Yo when the conditioning set 8 ().
Definition 2.2 (Markov Boundary MT))..) Given a set of TS variables X and target variable T in this set, (T,) is a Markov blanket of Tr iff: IL |M)M.. The set (T,)) is a Markov boundary of T, iff VM' C &, (V M')
Definition 2.3 (Multiple Time-series variable selection (MTVS)). Let MB(T;) reference MB of a target T. TS variable. The solution of the MTVS problem consists into finding all MBs, denoted by M. All M are information equivalent. For all MB, (T.) E .M the forecasting models T. = fe (T)) are pereollyppooroingggg... according to a metric (e.g., R²).
Example 2.1. We consider hypothetical water flow monitoring system involving two rivers A, B and a small hydroelectric station on river B. Let MTS X, Z E R be the inflows of the two confluent rivers A and B. The dam is controlled such that if Z, S zth. it does not produce energy. Otherwise, diverts a flow Zeh. to power production. The power production flow mixes with river A first for total flow Y f(Z,)+ + X,, then with the rest of river B for total flow T, X, Ze. with f(Z)) E {0, ). 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, none of the Markov boundaries of a target T, [BT21]. A variable X. is said irreplaceable iff it is part of all information equivalent Markov boundaries of T., i.c., X, e E M. A variable X, is said replaceable iff it is part of at least one MB(T;) but not indispensible, i.e. HMB(T), MB' (Tt) E M, X, E MB(T)) A X, & MB' (T)). A variable X, is redundant iff 3Z C X (X. X, KL T Z and VMB (T)) M, & (T),T)) while it is called irrelevant iff 1 X
If Causal Markov Condition and Faithfulness [CMISV-23] 23] hold, each target 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
3

---

(B) might not be unique. Jointly Gaussian random variables (RV) with singular covariance matrix or with deterministic relations do no respect faithfulness [Mee95, Lem07]. While faithfulness does not hold, the weaker Composition property relaxes the structure of faithful data. As 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, Y.Y T and conditioning set Z:
X Z and Y 2 = ZUYLT
Composition is a general property of the joint probability distribution of a set of RVs regardless of their temporal context, hence we dropped index from formal definitions. The reciprocal property is called Decompesition, 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 = fXX) of size m, where Y; = fi fx(o(()). with fa invertible and 0g a mapping from [[1,m]] to [1, ra]]. 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 XI Z and Y Z hold, the union X U Y contains only RVs that are pairwise independent from each variable in T given Z. This implies X
For probability distributions where all information equivalences are caused by invertible deterministic transformations between individual RVs (singletons), the MBs of any target interchangeable: the variables are equivalent regardless of the conditioning set. More generally, we define interchangeabil- ity for two Markov Boundaries, where each variable in MB( equivalent a variable in the other MB(T) conditioned on the remaining MB(T) variables. We assume that all MB are interchangeable.
Assumption 2.2 (Interchaneeabilit)). Two MBs of target T,, MB) and M MB'(TT) (T;) are inter- changeable iff:
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 = (fZ,X), Y)) satisfies Interchangeability, as . {2} (X,Y) Y) is cartesian product of equivalence classes.
AM (T)), MBB( E M, E MY E MB) (IX)MM
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) computationally optimized version (approximate) named Forward Equivalent (FE- Algorithm 2).
Algorithm 1 ChronoEpiilgii-FEE
Algorithm ChronoEpilogi-FE
Require: TS X. target T,, max lag L, thresh-
Require: TS X, target TS variable T, max
old params [a,, 7, o]
lag L, threshold params [a, 6]
set S +- FORWARD(X, T, L,a)
1: M + FORW-EQUIV(X, T, I.,o)
2: S + BACKWARD(X, T. L.S.7)
return M set of eq. Markov bound.
3: set M + EQUIV(X, T. L,S,8)
3:
4: return M set of eq. Markov bound.
In FBE, we first select informative TS random variables (RV) for predicting the values of target T, adopting a greedy heuristic (FORWARD routine, line 1). Then, 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-FEE 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
*Conditions leading to faithfulness violations have been studied from the standpoint of probability distribu- tions 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].

---

Algorithm 3 FORWARD
Require: TS X RNNMM, target T, max lag L,, stopping threshold
1: set S <- -7)
2: model m +- fit and save fe such that TE+1 ff((S,,)) TM fe(SM-L.M-)
3: repeat
4: R m.r.eiduuuls RR--
5: Snew f arg einx' (ea
>> Select Snew S.t. Snew V R
6: S S
7: model m' +
8:
model m t and save fe S
9:
until pvalue of Likelihood ratio test of (HO: m' = m) 2 a
b Verify Swew H Te41 Snew
10: return S Snew
Algorithm 4 BAAKWARD
Require: TS X, target T, time pred. t, max lag L, sel. variables S, stopping threshold y
1: repeat
2: for Saet E S do
3: if pvalue of Likelihood ratio test between Se(S) (Saa)) and fe(S) 2 then
4:
S 5 S {Sddt]
5: until no change
6: return
approximate the search of equivalent Markov boundaries using forecasting model residuals. In Example 2.1, the forward phase might select upstream redundant variable U first, then Z and X before terminating. The backward phase would test U J_L,IX, Z, removing from the selected set. Finally, the equivalence phase would test Y = X Z and produce the solution space M.
We evaluate ChronoEpilogi considering AutoRegressive Distributed Lags (ARDL) (linear) forecasting model [HPS84, PSSO1]. An ARDL model of orders 9 uses lags of both the target (T) and other TS RVs (X) predictors, with T. ao Σ.= a,TTi- + Σ Σ byxX; + Et. The model includes autoregressive terms (a,) and other explanatory variable terms (b,). It is generally required that e, 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 Maximum Likelihood estimation model. Consequently, this estimation has the advantages of the OLS method, with guaranteed convergence and fast computation.
Hereafter, detail sub-routines of ChronoEpilogi. The forward phase (Algorithm 3) iteratively builds a first selected set starting from the past L lags of (line 1), then incrementally selects new RVs of X. At each iteration, selected variable (Snew) is the maximizing the statistical importance of Pearson correlation between all windows of length L (the predictors) and the residuals of model m (line 5). Algo. (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 (line 2), and any non-selected variable Sd (line 4). In line 5 we compare the equivalence replacing the two Rys in a new forecasting model. If this latter is comparable with the baseline, we store the equivalent variable sets in dictionary (line 6).
5

---

Algorithm 5 EQUIV
Require: TS X, target T, time pred. I, max lag L, sel. variables S, equivalence threshold à
1: dictionary Q < {:}
2: for S E E
3: Q[S] f
4: for 5d (X (.S) do
5:
if pvalue of LR test between U (8)) and fe((S)) U 8 then
6:
Q[s] + Q[ ( 7: M [
(
8: return M
Algorithm 6 FORW EQQIV
Require: MTS X, target T, maximal lag L, stopping threshold a, equivalence threshold
1: set S <- M
2: set R X S
3: dictionary Q 4 {:}
4: model m + fit and save fe such that T2+1 fe(Si,,))................))))))) TM fe(SM-L.M-)
5: repeat
6: Res m.residuals M-L
7: Snew - arg (X'ER) Lag-Pearson Pearson- Res, X', L) Select Snew s.t. Snew K Res
8:
S + SU S (Sew
9:
R R+ R (S new
S
10:
[[S f FIND-EQUIVALENCES/ Res, Snew, R)
Test C' =Ree Snew for C R
R + R Q[Snew]
12:
m' <- m
model m <- and save fe S
14:
until pvalue of Likelihood ratio test of (HO: m' = m) 2 a
D Verify Snew H.. T Saew
15:
return S S
Once the equivalent Markov Boundaries M are obtained, we can characterize irreplaceable variables and replaceable variables, 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 we report the pseudocode of FORW EQUIV routine, which is used in FE select informative variables and to estimate equivalent sets over forecasting model residuals. Algorithm in Appendix details the residual- based equivalence search REED--EIIVVLENEEEE employed during the forward phase (Algorithm 6, line 10). Such routine tests statistical redundancy by modeling residuals with model ge(.). In this sense, if non-selected variable X is equivalent to the tested variable, Snew is added to the equivalent (line 7).
Complexity analysis Given TS X NxM, and denoted by S the selected TS variables in the forward phase, both ChronoEpilogi versions (FBE and FE) run (N]SS) conditional tests. In the FORWARD phase, the sub-routine Lag-Pcarson-pva (Algorithm 7) runs in O(L log(n)) if Fast Fourier Transform is adopted [MZZ+222] where L is the number of lags to predict a target and M L the size of residuals (number of predicted targets). Fitting ARDL model through matrix inversion requires O 3 M). Overall, the forward, backward, equivalent phases take O(SS N L M +|s 3 M). O(SS) I M), O(IS) 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 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 correct" [PNBT07, BT21]. We should stress this requires considering different statistical tests for different distributions [SP20]. We prove that FBE is sound any of TS variables returned solution in M a Markov blanket. In other words, our algorithm does not return false positives. We
Alternatively, algorithms might rely oracle producing variable with given properties [SLLA13].

---

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 proof' 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 M of a target Tr 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.l. Z, and normal jointly with all other variables. We notice that residual-basee 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 conditioning set. Let R. be the residual variable corresponding to the regression Tf on Z. Assume that the modeling achieves independence of the residuals: R.I. Z, then: T T = R.L.X. Additionally, if all candidate variable X, R = R Z. then R 7 T
Experimental Framework
To build synthetic datasets with multiple Markov boundaries that respect com- position and interchangeability, we build MTS starting from faithful distributions. Each MTS starts from 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, 1000, and maximal lag [1,5,10], we sample 10 datasets per configuration. We verify that the noise intensity varies, as the R² of a model trained with a MB ranges from 0.02 to 0.9 over the different data instances.
Synthetic dataset:
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 [LYSLI8] (transport). They are commonly used in recently proposed deep forecasting models 2J, ZY22, WCWW22J. 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, GrsupLasso [NST13], and with no selection in the case of Electricity, Traffic, and Solar. For forecast- ing, we use the ARDL model [PSSO1], as it is 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 recent benchmark [SWX+233. The tuned parameters are described in Appendix (Table 3). We report the performance of the best forecasting model for each target.
7

---

T846, T525, T572, T316, T712, T801, T810, T577, T791
T426, T868. T56. T808, T34, T756
T264
T585
T646
T12
T767
T437
T368
Targel T264
T457. T315, T299, T333, T547. T62, T196, T138, T295, T488, T148, T309, T293, T628, TS37
Figure 1: Multiple solutions of ChronoEpilogi on target T264 of Traffic. MB are build by combining irrellaceable variables (in blue) with members of the replaceable equivalence class (in yellow).
Evaluation metrics We evaluate the quality of Markov boundaries based three criteria: computa- tion time, predictive performance, size, and on synthetic datasets, fl-score using causal ground truth. Predictive performance is measured by R², RMSE and M MPE.. 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 fl-score of identifying respectively irreplaceable and replaceable TS variables. Since GroupLasso returns a unique solution, we also compute the fl-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, fl 2TP4FP+FN
2TP
SHAP explanations To assess the impact of replaceable variables 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 window of 10 timesteps. We then compute the SHAP values of random holdout samples using the KerneISHAP approximation [LL17]. which we sum along the time axis to obtain the overall importance of variables (per additivity property of SHAP values [KV2)))]). We consider the average of those values as global explanations [BSC+21]. We obtain average SHAP values $100 corresponding the global importance of each variable in the MTS number i. The total importance of each equivalent variable set S becomes Stot = Eies | Then, we compute the contribution percentage of the most important variables (until top-14) in the equivalent set |s:/ Stor| |. For comparison with SHAP explanations over a unique MB, we create reduced MTS where remove all but one replaceable variable per equivalent set. For instance, in MTS i variables [1, 4, 16, 78} are equivalent, we remove {4, 16, 78} and compute the reduced importance sti of the representative 1. We obtain instability results by repeating model training and explanation over different túme-contiguous) data splits of each MTS. Stability is measured as proposed by Nogeira et al. [NSB18], applied to the top JM B) variables of each explanation, with M B] 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 fl-score (Table 1), with stable performance as the number TS increases 2b). Claim 1: Variant FE has a 70% speedup compared to FBE with only decrease of causal fl metric of 0.05 (Table 1), due to less precise equivalence discovery phase leading to lower replaceable fl fl-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 fl-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 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) less than 4% of each MTS. The maximal selected size is below 5% all except Electricity (11%). The solutions have close to identical R² in Electricity, and spread over 0.02 around the median for Traffic (Fi..2c). Fig.1 illustrates compact representation of the solution space for one target of Traffic. Claim 4: On averaging across all targets the multiple solutions produced by ChronoEpilogi-FF have similar predictive performance to the unique solution produced by GroupLasso (GL). We use the total 50 targets to conduct a Friedman- Friedman-Nemenyi Nemenyi statistical test on method ranking [Her20], and conclude that GL and FE are not statistically different, but Claim variable selection do improves models' performance compared to models trained on the original MTS (full analysis Appendix E.2). Claim 5: 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

---

Computing time (5)
1.0
Causal metrics
1.0
1000
0.8
metric
500
0.6
"Ge as 4c+01 1e-01
19+19 2e+09 20 30-00 2c+05 4e-02
1000
0.7
10a
10a
0.4
2e=00 1c104 30-03 1c404 26+01 7e 26+09
Be
30 000 1et15
Num ber of TS (a)
Number 1S
poms-bay
metria
traffic
electricity
(a) computing time (b) fl score for replaceable (rF1) and irreplaceable (iF1) variables of FBE (solid lines) and FF (dashed lines) over the synthetic dataset. FE has up 70% speedup while losing at most 0.1 in rF1.
(c) Boxplot of the R of the 20 sampled MB for each target where ChronoFpilogi finds multiple MBs. Red dots figure the performance of Groupl asso. The number of MB detected reported for each target.
Figure 2: Performance of ChronoEpilogi versions on (a)()) 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 30% increase in causal fl-score.
Time(s)
R
MB size
Number of MB
causal fl sc-re
FBE
375
5.69:13.33 0.373
5,87c+05133.7++06 6.87c+07:1.222+08
0.99 0.064
FE
118:178
0.373::.2366
5.74:33.1
0.94 0. 188
GL
120+:19
0.337+0.239
7.50
0.57 0.349 NA
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.
few Shapley values instability over replaceable variables Claim 7: Models tend to select 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 increasingly shared as the number of equivalent variables grows. Claim 8: The most important variable of equivalent set is highly unstable to data resampling. The stability [NSB18 18] 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 MT'. 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.2c-13.
Related work
Most of the state-of th- art variable selection algorithms for MTS have focused on selecting 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 the number of TS due to distance computations, which are proven to be impractical for large dataset sizes. We selected Groupl Lasso [NST13] as baseline in our experiments, as it can specify groups of TS [LNZ]5][HMKII] 10] while applying Lasso-type optimization with linear modeling. It is worth noticing that Lasso-type algorithms have been extended to identify multiple solutions i.i.d data [PLCTI8] 88] without, however, proving any formal property.
As 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 MTS data. Theoretically sound methods focus on different distributional properties and assumptions. TIE* (SLLA13] assumes that a single Markov boundary discovery algorithm for non-faithuul data can be called as an oracle. All MB(T) are discovered using exhaustive combinatorial conditional tests (O(s.n") where is the maximal size of a considered MB and n the number of variables [BT21]). TMFBS [BT21] proposes forward-backe MB discovery algorithms that aims to decrease the number of redundant

---

Table 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 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.
Dataset
TSS
R.2
imse |
mape
size
time F/E
#model
#MB
Electricity
FE
0.940
0.226
2.388
10.9
31
/2411
6/0/4
lc+14
Electricity
GL
0.934
0.236
2.649
6.5
127
3/0/7
NA
Electricity
NS
0.863
0.349
4.935
321.0
NA
0/0/10
NA
Solar
FE 0.985
0.109
0.664
5.1
3/261
0/0/10
1
Solar
GL
0.984
0.111
0.623
4.7
119
0/0/10
NA
Solar
NS
0.968
0.159
1.607
138.0
NA
0/0/10
NA
Traffic
FE
0.783
0.442
39.745
11.7
12
1248
6/0/4
7e+24
Traffic
GL
0.797
0,431
28.130
79.1
211
3/0/7
NA
Traffic
NS
0.740
0.491
42.371
863.0
NA
1/0/9
NA
PEMS-BAY
FE 0.860
0.358
1.208
4.1
143/62.0
8/0/2
6e+4
PEMS-BAN BAY
GL
0.867
0.355
1.217
37.7
765
10/0/0
NA
PEMS-BAY
NS
0.820
0.957
1.343
325
NA
9/0/1
NA
METR-LA
FE
0.886
0.374
1.121
3.8
100/2246
8/0/2
2e+4
METR-LA
GL
0.864
0.401
1.230
41.9
306
7/0/3
NA
METR-LA
NS
0.896
0.363
1.054
207
NA
10/0/0
NA
conditional independence tests operated by TIE* to O(n"). KIAMB an iterative heward-backwarr algorithm [PNBT07] that requires composition and heuristic selections to find one MB of a target, with every MB having a non-zero zero probability of selection. Identifying all MBs might require an exponential number of runs of KIAMB [SLLA13] 1.3]. In construct, ChronoEpilogg computes provable equivalent subsets of variables in O (ns)) conditional tests under composition and interchangeability assumptions, while the forward- barrd--bbwwrrdddd phase is based 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 [LLZI0a] [LLZI0b] (see [SLLA13] for an in-depth review). Finally, several works highlight the importance of structured representation of multiple MBs for interpretability concerns, especially in the presence of an exponential number of MBs of target of interest [TLT18][LPL [ ChronoEpilogi is the first algorithm that provides compact representation of mutually equivalent variables for MTS. Redundant and irreplaceable variables can be easily distinguished at first glance (see Figure 1).
Conclusions
Overall, the paper's contributions are (a) the adaptation time series and further development of the theory of multiple variable selection, (b) the design of the first variable selection algorithm for time-scries 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 surrogate to a model using all available time series (assuming it 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 larger scope of experiments with more real and synthetic datasets of varying size, statistical properties, and data types (e.g., discrete time time-serie)). In addition, different variants of the main algorithm could employ non-linear models to compute residuals and non-line 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/scasonality. different train-test distributions, or change points.
10

---

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 dependent: More accurate approximations to shapley values. Artift- cial Intelligence, 298: 103502, 2021.
[And13] Holly Andersen. When to expect violations of causal faithfulness and why it matters. Philosophy of Science, 80(5)6722668, 2013.
[AST+10M C.F. Aliferis, A.R. Statnikov, I. Tsamardinos, 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 Muchine Learning Research, 2010.
[AST+ 10b] C.F. Aliferis, A.R. Statnikov, I. Tsamardinos, 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 Bizarto. Timeshap: Explaining recurrent models through sequence perturbations. In Proceed- ings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, KDD 221, page 2563-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- 1333- 1434, jul 2021.
[CMISV+ 23] Manuel Castro, Pedro Ribeiro Mendes Júnior, Aurea Soriano- Vargas, Rafael de Oliveira Werneck, Maiara Moreira Gonçalves. Leopoldo Lusquino Filho, Re- nato Moura, Marcelo Zampieri, Oscar Linares, Vitor Ferreira, Alexandre Ferreira, Alessandra Davólio, Denis Schiozer, and Anderson Rocha. Time series causal re- lationships discovery through feature importance and ensemble models. Scientific Reports, 13(1): ((1):11400 1402, Jul 2023.
[GBW+21] Rakshitha Godahewa, Christoph Bergmeir, Geoffrey I. Webb, Rob J. Hyndman, and Pablo Montero-Manse 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 specifi- cation. 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.

---

[3] Ming Jin, Yu Zheng, Yuan-Fang Li, Siheng Chen, Bin Yang, and Shirui Pan. Multivari- ate Time Series Forecasting With Dynamic Graph Neural ODEs. IEEE Transactions Knowledge and Data Engineering. 35(9):9168 1111969998880 9180,, September 2023. Conference Name: TEEE 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 Venkatasuubammanian, 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, ICMI120. 2020.
[LALP21] Bryan Lim, Sercan Ô. Arı,, Nicolas Loeff, and Tomas Pfister. Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting, 37(4):1748-1766 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 predic- tions. 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 Hujjie Zhang. Ensemble gene selection for cancer classifi- cation. Pattern Recognition. 43(8):2763-2772, aug 2010.
[LZZ15] 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.
[L] 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()):545, sep 2023.
[LTB+16 16] Vincenzo Lagani, Sofia Triantafilouu, 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 Uncertainty in Artificial Intelligence, UA1'95, page 41 1-14 San Francisco, CA, USA, 1995. Morgan Kaufmann Publishers Inc.
[MZZ+22] Abdullah Mueen, Sheng Zhong, Yan Zhu, Michael Yeh, Kaveh Kamgar, Krishna- murthy Viswanathan, Chetan Gupta, and Eamonn Keogh. The fastest similarity search algorithm for time series subsequences under euclidean distance, Aug 2022.
[NSB18] 18] 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.
12

---

[Pea88] Judea Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988.
[PLCTI8] 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 1--22, July 2007.
[PSSO1] M. Hashem Pesaran, Yongcheol Shin, and Richard Smith. Bounds testing ap- proaches to the analysis of level relationships. Journal of Applied Econometrics, 16(3), 2001.
[RGFO17] Milos Radovic, Mohamed Ghalwash, Nenad Filipovic, and Zoran Obradovic. Mini- mum 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-5477, 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): 181-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-355 2015.
[SLLA13] Alexander Statnikov, Nikita I. Lytkin, Jan Lemeire, and Constantin F. Aliferis. Algo- rithms for Discovery of Multiple Markov Boundaries. Journal of machine learning research JMLR, 44:49-5666, 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-1558, jun 2020. Publisher: Institute of Mathematical Statistics.
[SWX+23] Zezhi Shao, Fei Wang, Yongjun Xu, Wei Wei, Chengqing Yu, Zhao Zhang, 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.0611 2023.
[TCP+22) loannis Tsamardinos, Pavlos Charonyktakis, Georgios Papoutsoglou, Giorgos Bor- boudakis, Kleanthi Lakiotaki, Jean-Claude Zenklusen, Hartmut Juhi, Ekaterini Chatzaki, and Vincenzo Lagani. Just add data: automated predictive modeling for knowledge discovery and feature selection. NPJ Precis Oncol, 6, jun 2022.
[TLTI8] Michail Tsagris, Vincenzo Lagani, and Ioannis Tsamardinos. Feature selection for sigh--imensional 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: .0477 [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.1 141.1990.. [cs].
13

---

[Whi80] Halbert White. A heteroskedaasicityycconnie consiseent covariance matrix estimator and direct test for heteroskedasticity. Econometrica, 48 1980.
[YLL2I] 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(1 (1:177777-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.
14

---

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 = X [), symmetric (X F Z + Y =r 2) and transitive.
A.1 Properties
Independence relationships validate certain properties in any distribution. Contraction and weak union [Pea88] link independences of sets variables to conditional independences:
Z. W be sets of variables, then XIL W Lemma A.1 (Weak Union). Let X,Y, K,Y,Z, X_LYZZ Z W.
Lemma A.2 (Contraction). Let X, Y, Z, W be sets of variables, then ll KLLYZ and X XI ZU Y
ILYUWIZ
Equivalent variables conditioned on 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, T be random variables and Z a conditioning set. If
X = Y Z, then:
B -
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 sets of random variables such that U U V is joint normal. Then:
ULV = YU W V.ULV
When the covariance matrix of a joint gaussian distribution is invertible, the distribution is faithful (factorising into 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 target T holds when for all X, Y, Z, when Eq.1 holds with T T.
Limitations of Composition Composition is essentially weaker version of "pairwise indepen- dence implies mutual independence". Mutual independence would require that the joint probability P[X, Y, T) can be factored into P X [PP()]P(TT) when X, Y, T are mutually independent. Instead, composition requires that when X. T and Y, T.T are pairwise independent, the joint probability F[X, Y. T] factorizes as FX,, (IT).)
The typical example of XOR operation on two Bernouilli variables X, with probability 0.5, and T = X xorY invalidates composition. Indeed, X and Y individually bring 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, 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 becomes dependent of them. Note that interaction is possible in the case of any of X and Y already being dependent on T.
15

---

Interchangeability Interchangeability asks for the set of multiple solution M to factorize as cartesian product of equivalence class of individual variables. This allows the solution space be quickly discovered by linear time algorithm.
at the level of group of variables, and not just individual variables. For instance, let Z be standard normal random variable. Define X = 1|2 > 0].Z and Y = 1|2 < 0 J. Since X has the value of >> 0 and Y the converse, then Z X Y. In this situation, for any larget, Z =T X, Interchangeability requires that Z cannot be part of a markov boundary T.
In general, equivalences can
possible to relax Interchangeabiilty 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 [SILAA13].
B Pseudocode of
routine
Algorithm 7 details the Pearson correlation based test. We note that in this version, p-value is obtained for each the lags, the minimum of which returned. This procedure returns p-values that likely lower than the actual p-value, due to using the minimum for previous and current model statistical difference is kept separate from the correlation routine, ensuring that Forward phase termination relies single p-value.
estimator. Hence, testing
Algorithm 7 Sog-Pearson-pval
Require: residuals R, candidate variable X' RR,, maximum lag L
1: 4-1
2: for 2/ L do
3: P min(p, Pearson-corr(R,22 XY,M -M--H-)) 1+1-1) return p
C Pseudocode of equivalence computation of version FE
The equivalence search (Algorithm 8) of the approximate FE version builds surrogate models Je to predict the residuals instead of the target T. Its structure correspond to the definition of Informa- tion Equivalence, in that it tests proxies for the two assertions X (ILRes] R and Snew IL Res X. Establishing equivalence for residuals not always guaranteed to ascertain equivalence for T (see Theorems), but empirically find good proxy for synthetic datasets.
Algorithm 8 FIND-EQUIVALENCES
Require: residuals Res, selected variable Snew. remaining variables R, maximum lag L, equiva- lence threshold 8
1: set E +- 0
2: for X R do
3:
< pvalue of Likelihood ratio test between Res ge U Unew)) and Res ge (Snew)
4:
if p 2 8 then
redundant given Snew to predict Res pvalue of Likelihood ratio test between Res = ge(X X U Unew) and Res ge (X))
5:
6:
p 2 8 then
> Snew is redundant given to predict Res
7:
E U [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 reurned by the FORWARD and by the FORWARD- EQUIVALENCES procedures is a Markov blanket Te.
16

---

Forward Soundness. We reuse the proof of KIAMB [PNBT07]. At each iteration, the association heuristic selects new such that Skew-Kk 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' E R. then C. By composition, this ensures that R, hence S is a Markov Blanket of the target T₁-
Theorem D.2 (Backward Soundness). Assuming composition (Assumption 2.1) and perfect tests for the correlation, lermination, and elimination, the set returned by the succession of the FORWARD and BACKWARD procedures is Markov boundary of 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, Assuming that every conditional independence test of T2 aaa
| S {Sae) } is correct. At each deletion step, assuming V V\ S_I.T. and IL (S) (Saa). In this case, by contraction (A.2), Sac Sea U (V SS)SS) \ (Saa): This guarantees that S \ {Sart) still a Markov blanket.
By the end of the backward phase, for all S E S, T S {S}. By contradiction, suppose that there is M C. Markov blanket of T.. Let X E S M, let Y S M \ {XX. Since M is a Markov blanket of T,, X JLT, and Y.LT; M. Applying composition, X, Y MM, then applying weak union (A.1), X X_ILT, M U 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 FBE). Assuming composition and perfect tests for the correlation, termination, elimination and equivalences any M MB,...d B,, in the set M computed by FBE is Markov blanket T
FBE Soundness (Thm.D.3) D.3). In the following, define useful notation S& where S is an indexed set from I to n and & S n,, as the set of elements of S indexed from 1 to k.
First, we remark that it is enough to prove that M B,, =r Sn. where B,, is one of the produced solutions M Bi..., M,, E EM. that case, since Sn a Markov blanket, for any variable of the MTS B E X, then B I S,, holds. By Lemma A.3 applied to the equivalence M MB, =r X, we have B B.LT, Bn, Hence, M B,, is a Markov blanket of T.
To prove MB, =T Sn, we proceed by recursion. Suppose that MB Fy 1k-1 We seek prove that M Bx =T Sk S..1....... For simplicity of notation, in the following proof, denote Z = Sk. By definition of equivalence, we need to prove:
MB LT, Sk Z M..
The first condition true since S, is Markov Blanket of 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
Sk. Il TV MB Z
Sk Ta MB. Z
By using Lemma A.3 with equivalence M M ΣT Sk-1/Sk, Z (recursion supposition), equation (3) is equivalent to Sy Ty S,, *U(MB&) which is the condition to include M BN in the equivalence of St in our algorithm.
By using Lemma A.3 with equivalence M =y S, - (implied by the test in the algorithm), equation (4) is equivalent to - B Sk, Z. This last relation true due to the recursion supposition.
The initial condition holds since S, =T MB.,S Hence, MB, =T MB..
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 Markov boundary.
17

---

Corollary. All returned solutions are Markov Blankets (Thm. and at least one is Markov Boundary (Thm.D.2). Since by definition of Markov Boundaries, any Markov Blanket 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 selected set S find out all Markov boundaries of Ty. From the proof of Thm D.3, we can deduce that for whichever indexing of the reference set S, denoting S,....a. and for whichever equivalent boundary M produced by the algorithm, there is an ordering M,..M,, such that
Vi, M4 T S S
(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 By supposition, M is included in the solutions of the algorithm. Let us consider reference indexing I of S and M as is provided by the successive equivalent sets of our algorithm. Define any function f [[1, - 0, 1, resulting in intermediary set A.
There a permutation a [[1, ra]] + ||1, r2 such that f a a is monotonically decreasing. We reindex S (1) and similarily for Mj. Since (5)(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 M2 U (S') (S) is identical to A. By D.3. M2 \S)\ (Ss)) is a Markov boundary.
Direction s: suppose that for any correctly indexed M, any function f [|], + 0, induces set A and A is Markov boundary of Tr. particular, for any k [|,, the function fü if
else induces a Markov boundary S U Mg. Therefore, we deduce that SJLTTS U Mr. for all k which is exactly the conditions verified by the algorithm. Hence, each M2 is included in an equivalence set, and M is an equivalent Markov boundary produced by the algorithm.
Residual tests (T2m.3.2). Since Ti Z. then T. Z. Any deterministic transformation of T,, Z cannot add information so f(T,, Z) [L Z. Using Res Iz by contraction (Lemma A.2),
Res XX Z. Thus (decomposition) Res ILX.
Residual tests with composition (Thm.3.2). Starting from ResILX, Z and Resll by weak union (LemmaA. 1), Res JLX Z. Naturally, Res, Z JLL Z. Any deterministic transformation of Res, Z cannot add information and T; = (Res, Z). So T
O
Residual tests for normal variables (Thm.3.2). Since the distribution is joint normal, then Res T. (uD + E, a; Z.Z is also normal variable and V U (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 and Res JL with composition (Assump.2.1) Res 11x 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 and 0.05. As 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 logarithmic scale.
We describe the pipeline used for tuning and evaluation in Fig.3. cross-validation time series done in a fixed start window growing along time. The holdout set size consistently 1000 for synthetic MTS, and 30 percent of the MTS for real dataset. The smallest training window size must
18

---

Fold 1
Fold K
EVALUATION PIPELINE
Time Series Selection HP Optimization
Tuning
TIME
Train
TSS variant
Hyperparam Sampler
For each fold:
aappyy get solution
fit Forecaster (MSE) on
validation (R2)
Holdout
Validation
over selected series S
Select parameters with best
average
Evaluation
Train best configuration on Tuning (TSS algorithm and Forecaster)
TSS precision, recall (ff synthetic datal TSS size, time spent
Find the equivalent markov blankets
Compute R2 the Holdout
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
GroupLasse [NST13]
group regularization coefficient
ARDL [PSSO1]
Trend (no trend, constant, linear)
SVR
Kerne (rbf, sigmoid), regularization weight
DeepAR [SFGJ20]
Layer size, #layers, dropout rate Layer size, #layers, dropout rate, #attention heads
TFT [LALP21]
cover at least 40 percents of the Tuning split, obtain meaningful results and retain enough samples.
Over Fig.4, we examine how metrics vary when the thresholds are chosen to maximize the R of the forward- backwrd solution S. We observe that the forward and backward threshold depend on the R² 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² 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 tuning GroupLasss Lasso. Practictioners should expect spending more time tuning regularization-base solution than p-value-based solutions.
E.2 Statistical significance of results on real datasets
The statistical analysis was conducted for 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 normal for the populations ChronoEpiilogi (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

---

implaresible Eiscore
1.0
80
1.0
1.0
0.9
0.0
0.7
1e-20
1e:10
1e-05
0.001
0.02
1e:20
1e-10
1e:05
0,002
0.01
1e:20 1e:10
1e:05 0.001
0.01
forward_ threshold
thwerd_thresholdd
forward threshold
1.0
traplacaatee.n Sixcor
1.0
replaccable riscore
1.0
0.9
2
0.5
0.5
0.8
0.0
16/20
3.00
0..0002..00
0 oo? 0.003
0.01 0.05
16:20
2.00 0.01
0.03
backward
threshold
backward thresheld
Figure 4: Optimal hyperparameters choice influence on performance metrics for version FBE. Lines corresponds the forward phase threshold and backward phase threshold. Columns correspond to metrics (R2, flscore on the detection of irreplaceable variables, fl fl-score on the detection of replaceable variables.
R2 the contiquratins for coch MTS
0.25
0.15
0.10
0.00
ChrnoFpllogi-FFE
ChronoFpilogi-FF
Groupt
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.
1.00
1.00
1.00
1.0
Computinn time
Cousal metris
0.8
500
metric
500
1000
10a
10a
2500 5000 7500 Size training MTS
Number TS
Number of TS
Theoretical MB size
Theoretical maximal lag
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, 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
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
1.0
0.9
R2
0.8
0.7
metr-la
traffic
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 the central tendency of the populations ChronoEpi- logi (MD=002899+0)) 12, MAD=O. 1005, MR=2.580), GroupLasso (MD-0.2877-0.1 106, MAD=0.100 MR=2.320), and NoSelection (MD=0.442+-0.3U MAD=0.257, MR=1.100), Therefore, we assume that there is 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 significant.
E.3 ChronoEpilogi-FBE on Traffic and METR-LA
For comparison purpose, we ran the full version of the algorithm two of the real datasets. The performance is reported 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 increase in predictive performance. We confirm that version FE tend to select false positives, but keep the constatation that mutliple solution exist 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 outper- forming TFT models during cross validation. As we suspected that cross validation could influence negatively the TFT performance (earlier split having less data, and Transformer-base models notori- ously requiring consequent datasets), we tried a different pipeline where tuning and testing was done on single training/validation/ees split with 70,20, 10 repartition. The result are in Table 5.
21

---

Table 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
R²
imse
mape
size
time F/E
#model
#MB
ChronoEpilooii-FE
0.990
0.087
0.730
3/473
10/0/0
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
Target
1
2
1
Figure 8: PEMS-BAY dataset visualisation. The sensors (dots) 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 irreplaceable and the two orange ones 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 matter the predictive model, ChronoFppilogi does return multiple solutions this dataset. Additionally, the results using the modified pipeline do not change the method ranking of ChronoFpilogi, Groupl aasso and NoSelection.
E.5
Feature Importance under the existence of multiple Markov Boundaries
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 [AJI121]. 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, Vector Regressors (SVR) and XGBoost Regressors (XGB) for each of the 90 synthetic MTS with 100 covariates (denoted as full MTS), using window of 10 timesteps. We then compute the SHAP values of random holdout samples using the KemelSHAP approximation [LL17]. that we sum along the time axis to obtain the overall importance of variables (per additivity property of SHAP values [KV)).0]). We consider the average of those values as global explanations [BSC+21]. We obtain average SHAP values s1...f100 corresponding the global importance of each variable in the MTS number i.
tune and train Support
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 variables [1,4, 16, 78) are equivalent, we remove [4, 16, 78} and compute the reduced importance sri of the representative 1. Note that variable irreplaceable in the reduced MTS.
For each equivalent ffi.....ee in MTS i, we compute the contribution of each variable j [e)l
(o.) to the total absolute importance of the set, as EieT We plot the contribution of most important variable, second most important, etc (Fig. 9). We choose to use absolute importance values averages are not easily comparable with relative importance. We observe that importance not equally shared among equivalent variables. Models tend to select one or a few variables important predictors. Additionally, the importance of the top variable relative to the other diminishes as the number of equivalent variables grows.
10
top2
txp3
top?
lops
topt2
lap12
lap21
Contriboziin [SitSto1)
088
top3
top4
topa
topo
tepl
toppe
Contribution [S1/5a0]]
6a
Top
10p4
fnoil
Ino13
Dop14
top5
teplo
1065
Dop.00
0.-
64
00
12
6
10
Number of variables in the equivalent sat
Number of variables the equivalent set
Figure 9: Average proportion of importance of topl, top2. variables in 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 equivalent set highly unstable to data resampling. For each synthetic MTS with 100 variables, the stability (against data variability) of the most important (topl) 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 topl 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))) with M number of bootstraps
Theoretical max (most stable)
XGB
0.052
0.110
-0.05
SVR
0.020
0.089
-0.05
E.6
Hardware and Codebase
We make available anonymized codebase at https://githu:
wot.///o///hoooooppiiioo/// We run our experiments on servers runing Ubuntu 22.04.4 LTS, with cores, 1TB RAM, GPU Quadro RTX 8000 with driver version 550.54. 14 and CUDA version 12.4.
We include in the repository requirements file requirements. txt listing the necessary dependen- cies.
23

---

Table 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 reduced MTS with unique Markov Boundary. XGB models are computed and SHAP explanations produced. We compute the stability [NSB 18] of the IMB(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 [NSB 18]. 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 the last part of the introduction. Here, 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 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 other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Limitations related to the scope of the experimental study discussed in Conclusions
Justification:
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 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 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 few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.

---

The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails 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, worse outcome might be that reviewers discover limitations that aren't It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant 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, state the main assumptions in algorithm, namely Com- position 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 cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide 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 ex- perimental results of the paper to the extent that affects the main claims and/or conclusions the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: ChronoEpilogi is detailed in Section 3, where provide the pseudocode of core functions along with the critical details that could help a third-party implementation of algorithm. Reproducibility of experimental results is facilitated by the description of the executed pipelines as well as by the availability 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 perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same
25

---

dataset, provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c)
If the contribution is new model (e.g.. a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset instructions for how construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The source code and the data to reproduce the experiments are publicly avail- able at: https: / / /anonymous. 4open. sciince/rr/honoEpplnic We include guide- lines to run the scripts of our experiments.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.. cc/ public/guides/
cetbbbmsssoooooiioo PiPiiyyy for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / //nips. wwtppcccccc//cccccccccccccccctlt wppppc/gguudds/////ppppppopppapp gwtggguiee/////
gewidee/ooeeeuamiiiiii/ Cot///oooeeeieieeeeeeeeeeeeee iettmmmsssinn wwwssnoollli onPolicy) for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question:
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
26

---

Justification: In Section 4, we describe all the details about synthetic and real data used in the validation, along with cross-validation protocol, forecasting model adopted and evaluation metrics. In Appendix E we detail hyperparameter tuning including optimizer, and data splits for cross validation.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We report the standard deviation of metrics in Table and Figure 2a and Figure 2c. We use statistical testing to compare ChronoFpilogi, Groupl and NoSelection on real datasets, and measure Shapley value instability.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: describe the characteristics of the experimentation server in appendix E.6.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make into the paper).
27

---

9. Code Of Ethics
Question:
Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: hctps://neurrpe: httpeerr.......upiiiieeeeee Sptublic/Eticssuiddeee.
Answer: [Yes]
Justification: Anonymity of our submission guaranteed in all sections. The research proposed in this paper has no impact and potential harmful consequences for the society, nor does use personal data.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e-g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: We do not identify any negative societal impacts of our work. On the contrary producing reduced dimensionality forecasters opens interesting perspectives for Sustainable or Green AI. Due to space limitations we have been able detail these important issues.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (c.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We have not identified any risk in this sense.
Guidelines:
The answer NA means that the paper poses such risks.
28

---

Released models that have high risk for misuse or dual use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]:
Justification: We cite all models and third-part implementation used in the implementa- tion of our solution. Our own assets including code and produced datasets under free and open source license Apache V. 2.0, include modification of code component under this license.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswi thcode. www.coo//aaaaaaa com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors encouraged to reach out the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: dataset/code/model are available in an anonymized repository here: https: / /aanonymous 4open. science/s wcwincc/rr/honooBeeeeeeeeeeeee /ChonoEpilogi// along with comprehensive documen- tation following the Code and Data Submission Guidelines provided by the NeurlPS confer- ence committee.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/code/moodl as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects

---

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
We have not made use of crowdsourcing our evaluation.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: Our research does not concern experiments with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
30


Reviews:
[{"review_id": "fKrIh6VE5d", "reviewer": "Reviewer_hgHX", "summary": "This paper considers the problem of finding all minimal subsets of variables for optimal prediction of time series data, coining the term \"Markov Boundaries\" for those minimal subsets constituting Markov Blankets for the target time series variables in question. \nThe paper then proposes novel algorithms for this problem,FBE and FE, and prove the soundness and completeness of FBE (FE is a faster approximate algorithm) and empirically evaluate their performance. \nThe experiments are conducted using both synthetic data (with ground truth causal structure) and real world data, and compare the performance of the proposed algorithms against baselines of Group Lasso (GL) and No variable selection, with respect metrics including predictive accuracy, accuracy of causal structure learning (for synthetic data), computation time and solution size. \nThe results presented validate a number of claims about the proposed methods, the notable ones being that they are more accurate at uncovering the ground truth causal structure on synthetic data and FE roughly comparable to GL on real world data sets in terms of accuracy and computation time, sometimes significantly out-performing it in terms of solution size. \nThe problem formulating is apparently novel and interesting, and the proposed methods are also novel and theoretically sound (and complete). The empirical results show that they are at least competitive to the standard baselines. \nThis work would add some valuable knowledge and insights to the community with interest in causal modeling and interpretable learning in time series data.", "questions": "One wonders if there are ways to use Group Lasso to obtain multiple solutions of the type obtained by the proposed methods, for example, by performing multiple randomized runs with perturbation and aggregating the outputs. A comparison with such heuristics would be of interest.", "limitations": "The authors do mention that the claims are limited to the scope of their empirical evaluation which could be enhanced. \nIt would be additionally desirable to address the question of quantifying the statistical confidence of the outputs of the algorithms.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem formulation is novel and interesting and well motivated practically. \nThe proposed solution is novel and sound and complete.\nThe empirical evaluation is reasonable.", "weaknesses": "The performance of the proposed methods against the baseline of Group Lasso on real world data sets is not exactly compelling. \nMore clarify on the relative advantage of the proposed method(s) would be valuable. \nThe optimal algorithm, FBE, is not evaluated on real world data sets, which I assume is due to computational complexity. It would be beneficial to know if any evaluation (even if partial) could be performed on FBE on the real world data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you very much for all the constructive comments, remarks, and questions that helped us to improve the quality of our manuscript. Hereafter, we provide answers to all your questions and comments pointed out in the weaknesses and limitations sections of your review.  \n\n**Question 1) One wonders if there are ways to use Group Lasso to obtain multiple solutions of the type obtained by the proposed methods, for example, by performing multiple randomized runs with perturbation and aggregating the outputs. A comparison with such heuristics would be of interest.**\n\nOne could extend Lasso-type algorithms based on global optimization to identify multiple solutions. An example of preliminary work in this direction or cross-sectional i.i.d. data is at https://arxiv.org/abs/1710.04995. However, how to extend this or similar types of algorithms for the Group Lasso and time-series data is still an open problem. In addition, proving soundness and completeness for algorithms based on Lasso is also challenging.\n\n**Limitation 1) The authors do mention that the claims are limited to the scope of their empirical evaluation which could be enhanced. It would be additionally desirable to address the question of quantifying the statistical confidence of the outputs of the algorithms.**\n\n\nIn the current version of the paper, we evaluate ChonoEpilogi and its competitors with experiments conducted over three real TS datasets. For each one, we test 10 different targets.  This dataset number does not allow us to apply a significant statistical test on the algorithms’ predictive performance. We are currently experimenting with additional real TS datasets fitting multiple Markov Boundaries (MB) that allow us to compute a critical difference diagram of the statistically significant ranking of GroupLasso and ChronoEpilogi performance using the post-hoc Nemenyi test. \nMoreover, to assess the stability of Feature Importance methods like SHAP in the presence of replaceable variables, we are currently experimenting with synthetic datasets of varying redundancy degrees using the Wilcoxon signed rank tests. In the new datasets, results are similar to the ones shown so far and strengthen our claims.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ulcFsLRQAa", "reviewer": "Reviewer_PVJ5", "summary": "This paper handles the problem of selecting all the minimal-size subsets of multivariate time series variables such that the past leads to an optimal predictive model for the forecast of a given target variable, which is essentially a time series feature selection problem. Past algorithms have worked to select a single such subset. The proposed algorithm is relatively efficient, in that it does not take as much longer than finding a single subset as one would think, but leading to more insight and better \"Markov blankets.\"", "questions": "1. Line 176 should say \"forward\" rather than \"backward.\"\n2. I suspect that algorithm 3 step 9 and algorithm 4 step 3 should have $\\geq$ in place of $\\leq$.\n3. In line 260, how is TSS defined?", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The paper handles an important problem in a clever way and is explained quite clearly.\n2. The experimental results are convincing and actually include running time, which is often omitted.\n3. The theoretical results look correct, although admittedly I did not comb through the proofs in great detail.", "weaknesses": "1. The proposed algorithm was only compared against GroupLasso and not against any other among the related work mentioned in the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you very much for all the constructive comments, remarks, and questions that helped us to improve the quality of our manuscript. Hereafter, we provide answers to all your questions and comments pointed out in the weaknesses and limitations sections of your review.  \n\n**Question 1) Line 176 should say \"forward\" rather than \"backward.\"**\n\nThank you for finding this typo, we have fixed it in the paper.\n\n**Question 2) I suspect that algorithm 3 step 9 and algorithm 4 step 3 should have ≥ in place of ≤.**\n\nThis one is another typo since we stop when we fail to reject the hypothesis stating the models are different. The smaller the p-value, the more different the tested models. We corrected it. \n\n**Question 3) In line 260, how is TSS defined?**\n\nTSS stands for Time Series variable Selection and will repeat the definition of the acronym at its first appearance in the introduction.\n\n**Weakness 1) The proposed algorithm was only compared against GroupLasso and not against any other among the related work mentioned in the paper.**\n\nThis is due to the fact that competitor solutions for multiple MBs discovery assume for iid data hence do not provide open-source code for MTS data, and thus their implementation is challenging and a research work in itself. We should note that time series require specific precautions in terms of estimators since data is not iid. Moreover, the only algorithms with theoretical guarantees for multiple MB discovery are TIE* and KIAMB which require an exponential runtime (see related work) and could not run for the high-dimensional TS datasets of our experimental evaluation.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ynEwtfT4ew", "reviewer": "Reviewer_UZ5R", "summary": "The authors presents **ChronoEpilogi**, an algorithm for multiple feature selection in multivariate time-series (TS) forecasting. This approach aims to identify all minimal-size subsets of TS variables (Markov Boundaries) that optimally predict a given target variable's future. The key contributions are:\n\n1. **Theory Development**: Introduces the problem of multiple time-series variable selection (MTVS) and the concepts of informational equivalence and interchangeability in TS data.\n2. **Algorithm Design**: Proposes ChronoEpilogi, which identifies all Markov Boundaries under broad, non-parametric assumptions.\n3. **Experiments and Results**: Demonstrates ChronoEpilogi's scalability to hundreds of TS variables and its effectiveness in reducing the number of variables while maintaining or improving forecasting performance.", "questions": "It would be great if the authors can discuss about the computational complexity and the limitations stemming from their assumptions", "limitations": "See weakness above", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Some of the key strengths of the paper are:\n\n1. The paper proposes a novel theoretical foundation for multiple feature selection in TS data including concepts like informational equivalence and interchangeability. Combining these concepts the authors have been able to propose an empirical method to detect Markov Boundaries\n2. Furthermore, the proposed algorithm ChronoEpilogi is shown to handle large datasets effectively, making it suitable for real-world applications with numerous TS variables. This scalability is crucial to real world usage of the algorithm\n3. Another key contribution for the paper is that ChronoEpilogi aims to identify all minimal-size subsets, offering multiple valid forecasting models and insights.", "weaknesses": "While being a very interesting paper, there are some avenues for improvement: \n1. While the authors discuss the scalability, the algorithm's computational complexity seems to be high for very large datasets, potentially limiting its practicality.\n2. Some of at the assumptions that the algorithm relies on such as Compositionality and Interchangeability may not hold in all real-world scenarios, potentially affecting its generalizability. The authors should consider discussing the limitations in their papers and how well the assumptions hold in practice\n3.The authors have provided thorough experimentations. However, to justify the practicality of the algorithm, it would be interesting to report additional validation on more diverse and complex real-world datasets", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you very much for all the constructive comments, remarks, and questions that helped us to improve the quality of our manuscript. Hereafter, we provide answers to all your questions and comments pointed out in the weaknesses and limitations sections of your review.  \n\n**Question 1) It would be great if the authors can discuss about the computational complexity and the limitations stemming from their assumptions.**\n**Weakness 1) While the authors discuss the scalability, the algorithm's computational complexity seems to be high for very large datasets, potentially limiting its practicality.**\n\nChronoEpilogi complexity essentially depends on the number and the complexity of tests run by each phase of the algorithm. Note that a model is constructed only at the end of each iteration step and not for each candidate subset of TS variables. The exact version of ChronoEpilogi, termed FBE, is indeed expensive (see experiments with synthetic data in Figure 2) but allows us to prove soundness and completeness theorems. For this reason, we have also proposed a greedy approximated version termed FE that avoids running the backward phase and computes the equivalence directly in the forward phase, thus pruning new model generations for all the selected TS variables discovered in the forward phase. As we can see in Table 1, in synthetic datasets fitting multiple MBs, FE runs in 30% of FBE time on average, with only a decrease of causal f1 metric of 0.05 while in real TS datasets, FE is able to find the first MB(T) up to 2 orders of magnitude faster than GroupLasso. Such results suggest the practicality of the proposed algorithm.\n\n**Weakness 2) Some of the assumptions that the algorithm relies on such as Compositionality and Interchangeability may not hold in all real-world scenarios, potentially affecting its generalizability. The authors should consider discussing the limitations in their papers and how well the assumptions hold in practice**\n\nPlease note that without any assumption regarding the joint probability distribution, causal variable selection is NP-hard even for linear regression problems (https://www.tandfonline.com/doi/abs/10.1080/00949658208810560). Rather than simply dropping the strong Faithfulness assumption we have considered weaker assumptions like Compositionality and Interchangeability allowing us to prove ChronoEpilogi soundness and completeness.  Several practical and useful data distributions satisfy Composability.  To that extent, in example 2.1, we show that any deterministic transformation satisfies composition.  In example 2.2, we show that Joint Gaussian Distributions also satisfy Composability.  We will also include in the appendix examples of distributions that do not respect Compositionality, starting with the typical example of non-faithful XOR operation between Bernouilli variables with p=0.5.  ChronoEpilogi (and its variants) is an effective tool that allows us to experimentally verify the Compositionality and Interchangeability properties in real TS datasets. To the best of our knowledge, our submission is the first work that reports multiple MBs on popular datasets used for MTS forecasting problems related to these two weak assumptions of our algorithm.\n\n**Weakness 3.The authors have provided thorough experimentations. However, to justify the practicality of the algorithm, it would be interesting to report additional validation on more diverse and complex real-world datasets**\n\nWe are currently extending our experimental study to two new TS datasets: METR-LA and PEMS-BAY (speed reading in road networks) that we found to exhibit multiple MBs (https://openreview.net/forum?id=SJiHXGWAZ). Current MTS datasets commonly used for forecasting do suffer from lack of domain diversity (https://arxiv.org/pdf/2310.06119) (https://arxiv.org/html/2403.20150v1). Please note that the challenge of the experimental evaluation of our work is not to measure forecasting performance with different models but the quantitative evaluation of ChronoEpilogi's multiple solutions against our main competitor namely GroupLasso. The newly produced results also confirm ChronoEpilogi's excellence and they will be included in Table 2 and the appendix of our paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "G5OG1DKMR2", "reviewer": "Reviewer_fxYh", "summary": "The authors propose a scalable algorithm called ChronoEpilogi that aims to select multiple subsets (Markov Boundaries) of time series (TS) features in order to better understand the underlying data generation process and to provide better explanations of downstream forecasting tasks. Through extensive experiments, the authors show that time series forecasting models perform better when fed with these subsets of TS features (individually) than when fed with all TS features.", "questions": "- Why did you choose ARDL over other models like RNNs which are also autoregressive models and are certainly more efficient (Equation 2)?\n\n- Is the model redesigned at each iteration of the Forward phase, or does it remain unchanged? The reviewer asks this question because the size  of the inputs $\\textbf{S}$ may vary at each iteration (see Algorithm 3 line 8)?\n\n- What can explain the outperformance of SVR over TFT and DeepAR with the Solar Dataset?\n\n- How do the authors expect their algorithm to perform when faced with sparse multivariate time series, given that missing value can affect causal inference?\n\n- The reviewer does not understand the relevance of reporting the standard deviation in columns MB size and Number of MB (Table 1). Could you please elaborate on this?", "limitations": "- The authors have identified the limitations of their work and plan to address them in their future work;\n\n- On line 176 it should be forward instead of backward;\n\n-  The reviewer suggests improving definition 2, which is the core definition of the paper.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "-  Originality: Although the problem addressed in this paper is not new and the proposed solution is based on the combination of existing methods/models, the originality lies in the fact that, unlike previous work, the proposed algorithms provide a compact representation of mutually equivalent variables for multivariate time series (MTS). In addition, the redundant and irreplaceable variables in MTS can be relatively easy to identify.\n\n- Quality: The document is well structured and the assertions are fairly well supported.\n\n- Clarity: Overall, the document is well written and pleasant to read. However, the reviewer suggests improving definition 2. What does V stand for?\n\n- Significance: The reviewer believes the proposed algorithm can be used as an alternative solution to provide explainability in diverse and sensitive fields such as medicine and autonomous vehicles.  Indeed, in these fields, when it comes to forecasting tasks, identifying the  the time series  features that influence the decision is as important as model accuracy.", "weaknesses": "- The experiment is not conducted with multivariate time series that have a high rate of missing values. This is a crucial aspect that the study should have taken into account, as missing values are inherent in time series and may affect causal inference.;\n\n- The authors simply identify subsets of time series variables without providing concrete explanations that could have strengthened their claim. For example, the authors should have taken a few subsets (Markov Boundaries)  in any task and shown how they are actually relevant to the target;\n\n- The reviewer understands the usefulness of greedy heuristics to speed up the algorithm. However, this heuristic has the disadvantage of providing suboptimal results. Although the authors demonstrate its effectiveness, it would be interesting in future work to test it on additional datasets covering different domains.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you very much for all the constructive comments, remarks, and questions that helped us to improve the quality of our manuscript. Hereafter, we provide answers to all your questions and comments pointed out in the weaknesses and limitations sections of your review.\n\n**Question 1)**\n\nIn a preliminary phase of experimentation, we indeed planned to include LSTM models to tackle nonlinear tasks. Our observation was that applying the forward backward phases to find an MB(T) in nonlinear synthetic datasets (https://doi.org/10.1613/jair.1.13428) had poor causal recovery when using LSTMs compared to linear methods. We explained this by the variability of DNN training which has higher fluctuations in R2 compared to the increase in R2 expected due to adding a correct MB(T) variable. Obtaining sufficient stability might require more data, finer human analysis of model training routines, and extensive bootstrapping to remove fluctuations. Therefore, we did not include LSTM in ChronoEpilogi.\n\n**Question 2)**\n\nAt every iteration, we build a new model with the variables selected at that step. It is important to note that ChronoEpilogi avoids the creation of a new model for every possible variable set by examining the correlation of the model residuals with the variable set selected at the previous step. \n\n**Question 3)**\n\nThe discrepancy is likely to be influenced by cross-validation, where the training set in the earliest split is 50% smaller than the one of the last split. SVR results are stable across folds, while TFT models have more unstable results. We plan to add, in the appendix, new results for this experiment with a different cross-validation strategy to confirm our hypothesis.\n\n\n**Question 4)**\n\nChronoEpilogi requires two main operators: (a) a modeler, and (b) an estimator of pairwise correlation between the residuals and a Time Series variable. One could impute missing values before constructing a model or computing correlation (https://folia.unifr.ch/unifr/documents/309429). Alternatively, one could employ models that inherently fit models in the presence of missing values (https://doi.org/10.1038/s41598-018-24271-9). The absolute performance of such a version of our algorithm, as well as its relative performance against similar types of algorithms on datasets with missing values, is an open research question for causal discovery methods (https://doi.org/10.1613/jair.1.13428). In this paper, we propose the first algorithm for Time Series Selection with Multiple Solutions, so that we can extend it in the future to handle datasets with missing values.  \n\n**Question 5)** \n\nIn Table 1, we observe that ChronoEpilogi discovers (on average) smaller solutions (with smaller variance) than GL. Standard deviation and mean are computed across the 270 synthetic MTS, where for each MTS ChronoEpilogi and GL are applied once. For each MTS, we obtain the size of the solutions found by GroupLasso and ChronoEpilogi (all solutions produced by ChronoEpilogi are the same size) by computing their mean and standard deviation. \nThe standard deviation of  MB size is relevant to our experiment, as it allows us to understand that GroupLasso computes feature sets with varying lengths compared to those of ChronoEpilogi. GL tends to include too many variables generating overfitting in linear models.\nWe also observe that GL tends to dismiss some irreplaceable variables that should have a low impact on predictive performance.\nEvaluating the standard deviation of solution size is also relevant to compare FBE and FE (approximate version). We observe that FE tends to create more solutions than FBE, with higher fluctuations. As reported in Table 2 (real datasets), both methods can provide parsimonious solutions.  To better appreciate the quality of the approximation offered by FE we will include in the appendix also the characteristics of the solutions computed by FBE in real datasets. \n\n**Weakness 1)**\n \nWe consider this as an important issue, but it is orthogonal to our current submission also as in recent empirical studies of causal discovery methods (https://doi.org/10.1613/jair.1.13428).\n\n**Weakness 2)**\n \nWe started investigating selected subsets in real TS datasets as Traffic (see Figure 1).  This study is interesting to reveal the reasons why multiple replaceable TS variables may exist, as for instance the spatial proximity of traffic sensors included in the dataset.  However, such root cause analysis could not be completed without the involvement of experts from different domains, a task that is not possible with the public TS datasets used in our experimental evaluation.  We are planning in the future to conduct a deeper analysis of proprietary datasets in the context of an ongoing research project. As noted in the paper conclusions, the discovery of multiple MBs aims to facilitate a better understanding (via appropriate visualization tools) of the underlying data generation mechanism, as well as 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.\n\n**Weakness 3)**\n\nCurrent MTS datasets commonly used for forecasting do suffer from lack of domain diversity (https://arxiv.org/pdf/2310.06119) (https://arxiv.org/html/2403.20150v1). We are considering datasets from other domains (such as maintenance and water flow forecasting), but those are proprietary and outside of the scope of this paper due to each dataset's specific challenges.\nWe are currently experimenting with two more additional datasets: METR-LA and PEMS-BAY (speed readings in road networks) that we found to exhibit multiple MBs. The new results will be included in Table 2 and the appendix of the camera-ready version.\n\n\n**Limitation 2)**\nFixed\n\n**Limitation 3)**\n\nWe have improved definitions 2.1 and 2.3, which we include in the Camera-Ready version of the paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "untddUnYmW", "reviewer": "Reviewer_5YNE", "summary": "The authors consider the problem of feature selection when forecasting multivariate time series. They propose a novel algorithm called ChronoEpilogi based on identifying a Markov boundary of the time series variables. They experimentally and theoretically validate the findings.", "questions": "1. The abstract states, “Identifying these subsets leads to gaining insights, domain intuition, and a better understanding of the data-generating mechanism.” Is this claim supported by the experiments or otherwise in the main text?\n2. I have the same question about “identifying all such minimal-size, optimally predictive subsets is necessary for knowledge discovery and important to avoid misleading a practitioner.” I guess feature selection helps with interpretability, but I think the main text does not discuss it in detail.\n3. Please comment on weakness #1, which I have listed above.\n\n\n## Additional comments:\n\nAdditional comments\nI am really surprised that the primary area is interpretability and explainability, even though the paper almost does not discuss these aspects of the solution.\n\nWhat is V in tvs?\n\n“under Composition, Interchangeability, and other broad, non-parametric assumptions” – could you list all the assumptions here?\n\nL102 : conditional independence is not written well\n\nL205 RVS is not properly capitalized\n\nEq 2: what does the dot mean? Is it a misprint?\n\nTable 2:\nWhat do you mean by size?\nNo selection yields worse results. I wonder, what happens if you use a stronger base model for forecasting, e.g., LSTM?\n\nRe: Statistical significance:\nThe paper does not report statistical significance in Table 2.", "limitations": "In weaknesses, I’ve already commented that some crucial assumptions have not been discussed. These assumptions are critical to the framework, in my opinion.", "rating": 5, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. A significant problem to tackle, \n2. Good formalization of the problem,\n3. The paper is generally well-written.", "weaknesses": "1. Not all limitations are discussed: for instance, the model assumes that the selected set of variables remains fixed over time. When we deploy time series models, it is important that the method should work with different train and test time segments. However, since the set of variables is selected, generally, it might not be applicable to other train/test splits and thus might cause issues in real-world use.\n\n2. In my view, the experiments could be more comprehensive. It would be beneficial to consider other forecasting models, baselines, and datasets to provide a more robust evaluation of the model's performance. \n\n3. Prior work discussion is incomplete: for instance, signature paths and feature selection should be discussed and compared to the method. Some specific examples include\n    - Cross-correlation analysis,\n    - Signature transforms https://arxiv.org/abs/1603.03788", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you very much for all the constructive comments, remarks, and questions that helped us to improve the quality of our manuscript. We hope this response addresses your concerns effectively!\n \n**Q1**\nWe argue that finding one MB primarily serves to build a surrogate model with improved performance. On the other hand, finding multiple MBs is essential for domain discovery and helps practitioners understand the underlying data generation mechanism used to train a model. This is about Data Explainability rather than Model Explainability and our experiments with real datasets (see Figure 1 and Table 2) provide evidence that it is an issue of practical concern. To the best of our knowledge, our submission is the first work that reports multiple MBs on popular datasets used for MTS forecasting problems.\n \n**Q2**\nThe main text contains a quantitative evaluation of explanations using Markov boundaries (MBs). We show that the ChronoEpilogi algorithm reduces the number of time series to 5% of the original size, with increased performance to no selection (by 0.02 to 0.08 R2) (see Table 2 in the paper). Multiple signatures (a.k.a explanations) provide a complete set of explanations. The MBs contain the direct causes of the outcome and facilitate the interpretation of the data generation mechanism. A qualitative and subjective evaluation of the benefits for interpretation to a human user would require a user study and is out of the scope of this paper. \nA qualitative and subjective evaluation of the benefits for interpretation to a human user would require a user study and is out of the scope of this paper. We are highlighting in the paper conclusions potential usages of the multiple solutions computed by ChronoEpilogi e.g., constructing an ensemble model trained on several or all MBs as a means to create forecasters more robust to noise, faulty sensors, or other systematic errors.\n \n**W1** \nWe agree that data stationarity (with a unique distribution over the entire dataset) is a limitation of our current work, as in most cause discovery algorithms (https://doi.org/10.1613/jair.1.13428). There exist indeed many scenarios where such an assumption is broken, including datasets with trends, datasets where train and tests are both stationary but do not come from the same distribution, or datasets with change points. This is an open problem in the causal discovery literature that deserves further research to identify in a first phase realistic evolution patterns of Markov Boundaries, unlike the strong hypothesis of a few previous works. Additionally, all experiments follow the Forward Chaining Cross Validation protocol. We will add these clarifications to the Camera-Ready version of the paper.\n \n**W2** \nRegarding competitors, algorithms for multiple MBs discovery assume iid data and do not provide open source code for MTS that we could include in our experiments. Their implementation is challenging as time series require specific precautions in terms of estimators since data is not iid. Moreover, the only algorithms with theoretical guarantees for multiple MB discovery are TIE* and KIAMB which require an exponential runtime and could not run for the high-dimensional TS datasets of our experimental evaluation. See also our global rebuttal answer.\n\n**W3** \nWe thank you for pointing out the Signature transforms work. However, we cannot understand how we could leverage Signature Transforms  (a Feature Extraction method) to perform Feature Selection with multiple solutions in Multivariate Time Series. Getting more suggestions and pointers will be valuable and could allow us to add more discussion and material in the new paper version. We will be very willing to understand how this ML technique can be relevant to the discovery of multiple equivalent solutions (subsets of features) discovery for time series forecasting (MBo).\n\n**AC1**\nChronoEpilogi is the first algorithm for Time-series Variable Selection with Multiple Solutions that opens new perspectives in Data Explainability compared to Model Explainability heavily studied in the xAI literature. We argue that existing feature importance methods like SHAP are challenged by the presence of multiple MB in the datasets. This is a new claim that we are planning to support with additional experiments that will appear in the camera-ready version. It is known that Shapley values currently suffer from the inclusion of unrealistic data instances when features are correlated, even for linear models. We are currently investigating 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 score of each equivalent set of variables is distributed among equivalent variables, hence leading to underestimations of an equivalent set importance, when considered individually. This attribution is unstable: on different data splits, different variables among an equivalence set obtain high importance. See global review. To the best of our knowledge, these results have not been reported in the existing xAI literature.\n\n**AC7**\nFor each MTS, we obtain the size of the solutions found by GroupLasso and ChronoEpilogi (all solutions produced by ChronoEpilogi are the same size) by computing their mean and standard deviation. \n \n**AC8**\nWe evaluated ChonoEpilogi and its competitors with experiments conducted over three real TS datasets. For each one, we test 10 different targets. We are currently experimenting with additional real datasets fitting multiple Markov Boundaries (MB) that allow us to reliably compute a critical difference diagram of the statistically significant ranking of GroupLasso and ChronoEpilogi performance using the post-hoc Nemenyi test. In the new TS datasets, results are similar to the ones shown so far and strengthen our claims.\n\nThe typos have been corrected.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all the reviewers for their constructive comments, remarks, and questions that helped us to improve the quality of our manuscript. Hereafter, we provide a summary of the main modifications and additions we will bring to the final version of the paper. 

**1) Extension of the empirical datasets and consideration of other forecasting methods and baselines.**

In terms of additional datasets, we are currently extending our experimental study to two new TS datasets (speed reading in road networks), namely, METR-LA and PEMS-BAY (https://arxiv.org/abs/1707.01926) that we found to exhibit multiple MBs. With 5 datasets in total and 10 targets modeled each, we can compute critical difference diagrams of the statistically significant ranking of GroupLasso, ChronoEpilogi, and NoSelection performance. The newly produced results also confirm ChronoEpilogi's excellence and will be included in Table 2 of the camera-ready version of the paper. Note that current MTS datasets commonly used for forecasting do suffer from a lack of domain diversity (https://arxiv.org/pdf/2310.06119) (https://arxiv.org/html/2403.20150v1).
Regarding additional forecasting models, we could include Nhits and DecoderMLP available in the dedicated library we are using (https://pytorch-forecasting.readthedocs.io/en/stable/models.html). Please note that the challenge of the experimental evaluation of our work is not to measure forecasting performance with different models but the quantitative evaluation of ChronoEpilogi's multiple solutions against our main competitor namely GroupLasso.
Unfortunately, comparison with additional baselines is not possible as competitor solutions for multiple MBs discovery assume iid data and how to adapt the methods for MTS data would require a paper in itself.
We will also provide FBE results on real datasets where possible.

**2) Further evaluation of ChronoEpilogi relevance to xAi.**

ChronoEpilogi is the first algorithm for Time Series Selection with Multiple Solutions and opens new perspectives in Data Explainability compared to Model Explainability which has been heavily studied in the xAI literature.  We argue that existing feature importance methods like SHAP are challenged by the presence of multiple MB in the TS datasets.  This is a new claim that we are planning to support in the camera-ready version of the paper and we submitted as a pdf file.  It is known that Shapley values currently suffer from the inclusion of unrealistic data instances when features are correlated, even for linear models (https://www.sciencedirect.com/science/article/pii/S0004370221000539). We are currently investigating 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 score of each equivalent set of variables is distributed among equivalent variables, hence leading to underestimations of an equivalent set importance, when considered individually. This attribution is unstable: on different data splits, different variables among an equivalence set obtain high importance. The importance of the top variable diminishes as the number of equivalent variables grows, sharing more and more importance with the remaining variables. We report the results of these experiments on the page attached to the global rebuttal. To the best of our knowledge, these results have not been reported in the existing xAI literature.

Meta-review: Five reviewers all recommend acceptance. The reviews are detailed, and the authors have responded carefully. The research is thorough and novel.

This paper is somewhat interdisciplinary, because it is about feature selection in the unusual context where the features are entire time series, and with explainability as a main motivation, and therefore it likely deserves acceptance. This paper was considered for a spotlight, but due to the extremely limited space for spotlights, and when considering all of the potential spotlights across the entire NeurIPS accepted paper set, and ranking them, this paper will be instead be considered as a poster presentation. We do hope the authors will present this original interdisciplinarity as such to the audience, as making the motivation and contribution of the work easy to understand will be important, so please consider this to the extent possible.
Decision: Accept (poster)