ChronoEpilogi: Scalable Time Series Selection with Multiple Solutions

Authors: etienne  vareille
michele  linardi
ioannis  tsamardinos
vassilis  christophides
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 11:32:27

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.

Main paper: 
 introduction: 
 In analysis tasks involving TS with hundreds or even tens of thousands of variables (e.g., in manufacturing, environmental monitoring, energy grids, etc.), selecting appropriate series for TS forecasting not only has the potential to improve models' performance, but also to gain intuition, discover new knowledge, and understand the data-generating mechanism (causal structure). Time-series variable selection (TVS) is defined as the problem of discovering a minimal-size subset of the available, measured TS variables whose past values optimally predict the future of a target TS variable T . Such a series subset is called a Markov Boundary (MB) of T in the causal discovery literature [Pea88] , denoted as MB(T). Hence, MB(T ) filters out the TS variables that are both irrelevant and redundant in forecasting T .Identifying a small-size MB leads to a computationally more efficient model for T when real-time predictions are required. Moreover, it facilitates all subsequent analysis operations such as modeling, explanation calculations (e.g., SHAP values, feature importance), visualizations, and interpretations [VALC23] . Perhaps most importantly, identifying the MB of every series could be the first step in causal modeling of the data-generating mechanism. Indeed, under broad causal assumptions, the MB(T ) contains the observable direct causes of T , i.e., the quantities that could causally influence and optimize the future values of T [LTB + 16]. It is not an accident then that variable selection is the first step in several causal discovery algorithms for time-series or cohort data [AST + 10b, AST + 10a]. However, we note that the MB(T ) may also contain variables confounded by latent ones (in general, variables connected to T by a collider path). Distilling the possible direct causes of T out of its MB(T ) requires further analysis with causal discovery algorithms, which is outside the scope of our work (readers are referred to [ADG22] for a recent survey).Both in theory and in practice, there may be numerous Markov Boundaries of T . This is, for example, the case if a series subset {X, Y, Z} is 1-to-1 deterministically related with another subset {A, B}, in which case, either subset could substitute for the other in an optimal forecasting model. We call such subsets informationally equivalent series. Equivalent series define an equivalence class. The total number of Markov Boundaries of T is exponential to the number of equivalence classes to which the members of MB(T ) belong. This is because we can construct a new MB(T ) by picking any of its subsets and substituting it with an equivalent one from its equivalence class. Determinism is not the only reason for multiple MBs in practice: when the sample size is finite it may be impossible to distinguish the true MB(T ) from some other subset that leads to a forecasting model of statistically indistinguishable performance. The presence of an (often exponential) number of multiple MBs has been established in cohort, cross-sectional data [TCP + 22, BT21, SLLA13]. For time-series data there have not been any algorithms that return multiple solutions to the TVS problem [YLL21, SLL + 15]. Identifying all subsets MB(T ) (referred to as the multiple TVS (MTVS) problem) is crucial. While finding any one MB(T ) suffices for forecasting, this is insufficient for knowledge discovery and model interpretation. A practitioner may find it misleading to construct an optimal and minimal-size forecasting model by filtering out certain series although there exists another MB(T ) that contains them. In essence, multiple MB(T )s indicate the presence of multiple causal models and explanations that fit the data equally well. Moreover, if two series X and Y are informationally equivalent w.r.t. forecasting T , a variable selection algorithm may inconsistently select between them during crossvalidation or bootstrapping -arbitrarily returning either X or Y at each step. Feature Importance XAI methods may suffer from similar instability, on top of already known misleading interpretations. For instance, Shapley values currently suffer from the inclusion of unrealistic data instances when features are correlated, even for linear models [AJL21] .In this paper, we define the novel problem of multiple time-series variable selection (MTVS), we introduce the concepts of TS informational equivalence and provide a taxonomy to characterize TS variables w.r.t. to their presence in multiple MB(T ): indispensable, variables that belong in all MB(T ), replaceable, variables that belong in some MB(T ) but could be replaced by some other series in the same equivalence class, redundant, variables that are informative but do not belong in any MB(T ), and irrelevant, variables that are completely uninformative. We describe the property of Compositionality of the data distribution that allows the construction of sound yet greedy MTVS algorithms. Additionally, we define the property of Interchangeability, which specifies that any two MB(T )s can be decomposed into pairs of equivalent single variables. Under the assumption of Interchangeability, one variable can be substituted with another variable to create a new MB. This property enables the design of algorithms that are complete in identifying all MBs without requiring exhaustive enumeration. This is particularly valuable in data distributions with an exponential number of MBs.We proceed with designing a MTVS algorithm, named ChronoEpilogi 12 that is sound and complete in terms of identifying and returning all MBs of T , under Composition, Interchangeability, and other broad, non-parametric assumptions. ChronoEpilogi extends previous variable selection algorithms designed for cohort data [BT21] to a TS forecasting setting and the identification of multiple solutions (MBs). It returns a reference MB(T ) and a list of equivalence classes. Each of the classes contains TS that are informationally equivalent and could substitute one for the other to construct another MB and an equally-performing forecasting model. Experiments on synthetic data demonstrate that ChronoEpilogi has near perfect average causal f1-score, with stable performance as the MTS dimensionality increases to thousands of available TS. Moreover, the greedy approximation of the ideal ChronoEpilogi achieves a 70% speedup with only a 0.05 decrease in f1 score for causal discovery (Claim 1). Furthermore, ChronoEpilogi variants are on par with Group Lasso performance [SFHT13] (Claim 4) -arguably the most scalable algorithm for the single solution TVS problem-, with the added benefit of providing multiple solutions in real datasets from the Monash archive [GBW + 21] like Electricity and Traffic (Claim 3). On both real and synthetic datasets the multiple solutions produced by ChronoEpilogi have similar forecasting performance to the unique solution produced by GroupLasso when averaging across all selected targets (Claim 5) while when causal ground truth is available ChronoEpilogi outperforms GroupLasso in terms of causal f1-score (Claim 2). Finally, both variable selection algorithms actually conserve or even improve the models' performance compared to models trained on the original MTS (Claim 6). We also investigate how SHAP explanations of regression models might misrepresent the role of variables belonging to some but not all Markov boundaries of the modeled target. We claim that the SHAP importance of each equivalent set of variables is distributed among equivalent variables, hence leading to underestimations of the importance of equivalent sets when considered individually (Claim 7). This attribution is unstable: on different data splits, different variables among an equivalence set obtain high importance (Claim 8). 
preliminaries: 
 We represent a multivariate time series by a set of univariate TS X = {X 1 , ..., X N }, where eachX i ∈ R M (1 ≤ i ≤ N )are regularly sampled observations of a time series variable. We denote with X i t ∈ R the t th occurrence of X i , and by X i t-L,t-1 ∈ R L the last L occurrences of X i before the instant t. Note that, in the following part, bold letters (e.g., X) refer to sets of TS variables. The TS forecasting task we consider is to predict the values of a target TS variable T t (T = X i for some i), where L ≤ t ≤ M , with a model f parameterized by Θ that uses up until L past timestamps (for each t) of the multivariate TS X. Hence, T t == f Θ (X t-L,t-1 ).Conditional independence of two TS variables X t and Y t ′ given a third one Z t ′′ , written asX t ′ |= Y t ′ |Z t ′′ , is defined as P (X t , Y t ′ |Z t ′′ ) = P (X t |Z t ′′ ).P (Y t ′ |Z t ′′ ).In our work, we assume stationarity, namely the conditional distribution of T t as a function of L previous lags of X does not depend on t. Stationarity implies temporal consistency, where the conditional independence relations do not rely on t. In this respect, we consider the principle of temporal precedence according to which causes (i.e., independent variables) occur before an effect (i.e., outcome). We exclusively test conditional independence relations between a target T t and variables with timestamps ranging between t -L and t -1. Due to stationarity, we can remove the index and denote X t-L,t-1 as X. Definition 2.1 (Information Equivalence (IEQ) ). Two TS variables X a and Y b are information equivalent (shorthand: equivalent) with respect to a target T t given conditioning TS variable set Z, noted X a ≡ T Y b |Z, iff they are made independent of T t by conditioning on the other and on Z itself. Formally:X a ≡ T Y b |Z ⇐⇒ X a |= T t |Y b , Z and Y b |= T t |X a , Z.We write shortly X a ≡ T Y b when the conditioning set is Z = {}. Definition 2.2 (Markov Boundary MB(T t )). Given a set of TS variables X and a target variableT in this set, M (T t ) is a Markov blanket of T t iff: T t |= (V \ M )|M . The set M (T t ) is a Markov boundary of T t iff ∀M ′ ⫋ M , T t ̸⊥ ⊥ (V \ M ′ )|M ′ .Definition 2.3 (Multiple Time-series variable selection (MTVS)). Let MB(T t ) a reference MB of a target T t TS variable. The solution of the MTVS problem consists into finding all MBs, denoted by M. All M B(T t ) are information equivalent. For all MB i (T t ) ∈ M the forecasting models T t = f Θ (MB i (T t )) are equally-performing according to a metric (e.g., R 2 ). Example 2.1. We consider a hypothetical water flow monitoring system involving two rivers A, B and a small hydroelectric station on river B. Let MTS X, Z ∈ R be the inflows of the two confluent rivers A and B. The dam is controlled such that if Z t ≤ z th , it does not produce energy. Otherwise, it diverts a flow z th to power production. The power production flow mixes with river A first for total flow Y t = f (Z t ) + X t , then with the rest of river B for total flow T t = X t + Z t , with f (Z t ) ∈ {0, z th }. In this situation, X and Y are deterministically related only when taking account Z, and information equivalent for T given Z only. Both X, Z and Y, Z are MB of T . In this case, solving the MTVS problem relates to discovering both subsets without misidentifying singletons and X, Y as MB.TS variables selected by a MTVS algorithm can be characterized according to whether they belong to all, at least one, or none of the Markov boundaries of a target T t [BT21]. A variable X a is said irreplaceable iff it is part of all information equivalent Markov boundaries of T t , i.e., X a ∈ MB(T t )∀MB(T t ) ∈ M. A variable X a is said replaceable iff it is part of at least one MB(T t ) but not indispensible, i.e., ∃MB(T t ), MB ′ (T t ) ∈ M, X a ∈ MB(T t ) ∧ X a / ∈ MB ′ (T t ). A variable X a is redundant iff ∃Z ⊆ X \ {X a }, X a ̸⊥ ⊥ T t |Z and ∀MB i (T t ) ∈ M, X a / ∈ MB i (T t ) while it is called irrelevant iff ∀Z ⊆ X \ {X a }, X a |= T t |Z.∀M B(T t ), M B ′ (T t ) ∈ M, ∀X ∈ M B(T t ), ∃Y ∈ M B ′ (T t ), ({X}∪M B ′ (T t )\{Y }) ∈ MExample 2.4. In the water flow example 2.1, Z is irreplaceable, and X and Y are replaceable by each other. The markov boundary structure M = ({Z, X}, {Z, Y }) satisfies Interchangeability, as M = {Z} × {X, Y } is a cartesian product of equivalence classes.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Limitations related to the scope of the experimental study are discussed in Conclusions• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
the chronoepilogi algorithm: 
 In this section, we present the details of our MTVS algorithm, named ChronoEpilogi. We propose two versions of the algorithm: (1) Forward Backward Equivalent (FBE -Algorithm 1 ) and (2) a computationally optimized version (approximate) named Forward Equivalent (FE -Algorithm 2).Algorithm 1 ChronoEpilogi-FBE Require: TS X, target T , max lag L, threshold params [α, γ, δ] 1: set S ← FORWARD(X, T, L,α) 2: S ← BACKWARD(X, T, L,S,γ) 3: set M ← EQUIV(X, T, L,S,δ) 4: return M ▷ set of eq. Markov bound.Algorithm 2 ChronoEpilogi-FE Require: TS X, target TS variable T , max lag L, threshold params [α, δ] 1: M ← FORW-EQUIV(X, T, L,α) 2: return M ▷ set of eq. Markov bound.In FBE, we first select informative TS random variables (RV) for predicting the values of a target T , adopting a greedy heuristic (FORWARD routine, line 1). Then, a backward phase iteratively removes redundant RVs from the selection (BACKWARD routine, line 2). Lastly, equivalent Markov boundaries are discovered (EQUIV routine, line 3) by checking if any of the selected RVs can be replaced by another one. In ChronoEpilogi-FE, we select equivalent Markov boundaries during informative RV selection (FORW-EQUIV routine, line 1). Such choice permits us to compute conditional independences on smaller TS variable sets. To further reduce time complexity, we also S ← S \ {S del } 5: until no change in S 6: return S approximate the search of equivalent Markov boundaries using forecasting model residuals. In Example 2.1, the forward phase might select an upstream redundant variable U first, then Z and X before terminating. The backward phase would test U |= T |X, Z, removing U from the selected set. Finally, the equivalence phase would test Y ≡ T X|Z and produce the solution space M.Algorithm 3 FORWARD Require: TS X ∈ R N ×M , target T , max lag L, stopping threshold α 1: set S ← {T 1,M -1 } 2: model m ← fit and save f Θ such that T L+1 = f Θ (S 1,L ), ..., T M = f Θ (S M -L,M -1 ) 3: repeat 4: R ← m.residuals ∈ R M -L 5: S new ← arg min X ′ ∈X\S Lag-Pearson-pval(R, X ′ , L) ▷ SelectWe evaluate ChronoEpilogi considering AutoRegressive Distributed Lags (ARDL) (linear) forecasting model [HPS84, PSS01] . An ARDL model of orders p, q uses lags of both the target (T ) and other TS RVs (X) as predictors, with T t = a 0 + p i=1 a i .T t-i + k i=1 q j=1 b j .X i t-j + ϵ t . The model includes autoregressive terms (a j ) and other explanatory variable terms (b j ). It is generally required that ϵ t is an iid centered normal noise, but this assumption can be relaxed as we can build estimators in the presence of autocorrelated noise and noise with heteroscedastic components [Whi80, HPS84] . Under sufficient assumptions, an ARDL model can be estimated using Ordinary Least Squares (OLS) procedures, and correspond to a Maximum Likelihood estimation model. Consequently, this estimation has the advantages of the OLS method, with guaranteed convergence and fast computation.Hereafter, we detail sub-routines of ChronoEpilogi. The forward phase (Algorithm 3) iteratively builds a first selected set starting from the past L lags of T (line 1), then incrementally selects new RVs of X. At each iteration, a selected variable (S new ) is the one maximizing the statistical importance of Pearson correlation between all windows of length L (the predictors) and the residuals of model m (line 5). Algo. 7 (see Appendix) contains the correlation computation pseudo-code, which iterates all the windows in a candidate variable X ′ (line 2) to compute p-values of Pearson correlation with the forecasting model residuals. The selection in the FORWARD routine terminates when the last forecasting model (built over S) is not statistically better than the previous one (line 9). We use Likelihood ratio test [Kin98] , as they suit ARDL models. The backward phase (Algorithm 4) iteratively removes redundant RVs from the selected RVs in the forward phase. The main loop of the algorithm tests each one of these RVs. It stops as soon as removing any variable degrades the predictive performance of the best-so-far forecasting models (line 3), ensuring that the produced variable set is minimal and equally predictive as the forward phase solution. We prove that the forward and backward phases provide an exact solution of the MTVS problem (See appendix D.1,D.2). The equivalent search phase (Algortihm 5) tests the equivalence between each selected TS variable S (line 2), and any non-selected variable S d (line 4). In line 5 we compare the equivalence replacing the two Rvs in a new forecasting model. If this latter is comparable with the baseline, we store the equivalent variable sets in dictionary Q (line 6).Require: TS X, target T , time pred. t, max lag L, sel. variables S, equivalence threshold δ 1: dictionary Q ← {:} 2: for S ∈ S do 3:Q[S] ← {} 4: for S d ∈ (X \ S) do 5: if pvalue of LR test between f Θ ({S d } ∪ S \ {S}) and f Θ ({S d } ∪ S) ≥ δ then 6: Q[S] ← Q[S] ∪ {S d } 7: M ← {S 1 } ∪ Q[S 1 ] × ... × {S n } ∪ Q[S n ] 8: return M Algorithm 6 FORW_EQUIV Require: MTS X, target T , maximal lag L, stopping threshold α, equivalence threshold δ 1: set S ← {T 1,M -1 } 2: set R ← X \ S 3: dictionary Q ← {:} 4: model m ← fit and save f Θ such that T L+1 = f Θ (S 1,L ), ..., T M = f Θ (S M -L,M -1 ) 5: repeat 6: Res ← m.residuals ∈ R M -L 7: S new ← arg min {X ′ ∈R} Lag-Pearson-pval(Res, X ′ , L) ▷ Select S new s.t. S new ̸⊥ ⊥ Res 8: S ← S ∪ {S new } 9: R ← R \ {S new } 10: Q[S new ] ← FIND-EQUIVALENCES(Res, S new , R, δ) ▷ Test C ≡ Res S new for C ∈ R 11: R ← R \ Q[S new ] 12: m ′ ← m 13:model m ← fit and save f Θ on S 14: until pvalue of Likelihood ratio test of (H0:m ′ = m) ≥ α ▷ Verify Snew ̸⊥ ⊥ Tt+1|S \ Snew 15: return S \ Snew, QOnce the equivalent Markov Boundaries M are obtained, we can characterize irreplaceable variables and replaceable variables, as irreplaceable variables are the unique members of their equivalence class. Conversely, non-unique variables in their equivalence class are replaceable. In Section D (see Appendix), we provide soundness and completeness proofs of the multiple solutions computed by Algortihm 5. In Algorithm 6 we report the pseudocode of FORW_EQUIV routine, which is used in FE to select informative variables and to estimate equivalent sets over forecasting model residuals. Algorithm 8 in Appendix details the residual-based equivalence search (FIND-EQUIVALENCES) employed during the forward phase (Algorithm 6, line 10). Such a routine tests statistical redundancy by modeling residuals with model g Θ (.) . In this sense, if a non-selected variable X is equivalent to the tested variable, S new is added to the equivalent set (line 7).Complexity analysis Given a TS X ∈ R N ×M , and denoted by S the selected TS variables in the forward phase, both ChronoEpilogi versions (F BE andF E) run O(N |S|) conditional tests. In the FORWARD phase, the sub-routine Lag-Pearson-pval (Algorithm 7) runs in O(L n log(n)) if Fast Fourier Transform is adopted [MZZ + 22],where L is the number of lags to predict a target and n = M -L the size of residuals (number of predicted targets). Fitting an ARDL model through matrix inversion requires O(|S| 3 L 3 M ). Overall, the forward, backward, equivalent phases takeO(|S| N L M +|S| 4 L 3 M ), O(|S| 5 L 3 M ), O(|S| 4 L 3 N M ) time respectively.The two variants of ChronoEpilogi, satisfy different properties. FBE (Alg. 1) is an ideal version of our algorithm with provable soundness and completeness under mild assumptions. FE (Alg. 2) is a greedy approximation of FBE. To simplify the proofs of FBE theoretical guarantees we rely on a practical commonly made in causal discovery algorithms that the conditional independence tests are correct4 [PNBT07, BT21] . We should stress this requires considering different statistical tests for different distributions [SP20] . We prove that FBE is sound as any set of TS variables returned as a solution in M is a Markov blanket. In other words, our algorithm does not return false positives. We then examine the conditions under which FBE is complete, i.e. all Markov boundaries are solutions in M discovered by the algorithm. This corresponds to the lack of false negatives in the equivalent Markov boundary discovery. Theorems proofs are given in Appendix D.Theorem 3.1 (Soundness and Completeness of FBE). Assuming composition, interchangeable Markov boundaries, and perfect tests for the correlation, termination, elimination and equivalences, FBE computes all Markov boundaries M of a target T t and only Markov boundaries.The BACKWARD, EQUIV and FORWARD phases rely on model-based independence tests that might be replaced by any conditional independence (CI) test. In other words, CI tests performed in line 5 of the FORWARD algorithm (Algorithm 3) compute correlations over forecasting model residuals. We argue that in the context of linear regression with joint normal variables, measuring the residual correlation is equivalent to a model-based conditional independence test. This is the case of linear models of joint normal variables for which R |= Z, and R is normal jointly with all other variables. We notice that residual-based tests are not necessarily equivalent to full conditional tests. Yet no conditional independence remains undetected by a residual test. Formally: Theorem 3.2. Given T and X, two TS variables and Z a conditioning set. Let R be the residual variable corresponding to the regression of T on Z. Assume that the modeling achieves independence of the residuals: R |= Z, then:T |= X|Z =⇒ R |= X. Additionally, if for all candidate variable X, R |= X =⇒ R |= X, Z, then R |= X =⇒ T |= X|Z.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: In Section 2, we state the main assumptions in our algorithm, namely Composition and Interchangeability, under which we prove soundness and completeness of our algorithm in Section D (Appendix).Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
experimental framework: 
 Synthetic dataset: To build synthetic datasets with multiple Markov boundaries that respect composition and interchangeability, we build MTS starting from faithful distributions. Each MTS starts from a Vector AutoRegressive process with 20 variables, with different maximal lags and Markov boundary sizes. We then make copies of some of the MB variables to obtain replaceable variables. Those copies can be randomly shifted forward to change the lag of the relation. The same process is used to obtain redundant variables, copying correlated features not in the original Markov boundary. Irrelevant variables are then added to the dataset by sampling other VAR processes. For each MB size [2, 5, 10], total number of variables [10,100,1000], and maximal lag [1,5,10], we sample 10 datasets per configuration. We verify that the noise intensity varies, as the R 2 of a model trained with a MB ranges from 0.02 to 0.9 over the different data instances.Real datasets We evaluate our approach on five forecasting datasets covering different domains: Electricity (consumption), Solar (production), S.F. In TS data, observations are generally not independent, so data splitting for forecasting tasks must ensure that the train split precedes the test split (Forward Chaining Cross Validation). Therefore, we split the dataset along time into a Tuning and Holdout set, respectively for training and evaluating feature selection algorithms and forecasting models. The Tuning set is itself separated into five folds along the time axis, to conduct hyperparameter optimization. We optimize each considered pair of TS selection algorithm and forecasting algorithm together, for maximal average predictive performance R 2 over all folds. Due to the consequent training time of deep learning models, it is impractical to tune TSS (Time Series variables Selection) algorithms together. We first tune each TSS algorithm with Support Vector Regression (SVR) models as proxies, then tune deep forecasters for the tuned TSS parameters.Baselines and Forecasting models We compare our algorithm with the only linear scalable baseline, GroupLasso [NST13] , and with no selection in the case of Electricity, Traffic, and Solar. For forecasting, we use the ARDL model [PSS01] , as it is a standard linear model for MTS data. Real datasets are forecasted using nonlinear models: Support Vector Regression (SVR) with nonlinear kernel, and deep forecasters like DeepAR [SFGJ20] and Temporal Fusion Transformer (TFT) [LALP21] . The input window size is 10 for synthetic MTS and 96 for the five real datasets, similarly to a recent benchmark [SWX + 23]. The tuned parameters are described in Appendix (Table 3 ). We report the performance of the best forecasting model for each target. Evaluation metrics We evaluate the quality of Markov boundaries based on three criteria: computation time, predictive performance, size, and on synthetic datasets, f1-score using causal ground truth.Predictive performance is measured by R 2 , RM SE and M AP E. When more than 20 multiple solutions are obtained by ChronoEpilogi, we sample 20 solutions among the solution space, evaluate each of them, and report the mean. When the causal graph is known, we measure the f1-score of identifying respectively irreplaceable and replaceable TS variables. Since GroupLasso returns a unique solution, we also compute the f1-score of identifying all causal variables (irreplaceable and replaceable). Formally, denoting TP the number of correctly detected variables, FN for undetected correct variables, FP for wrongly detected variables, f 1 = 2T P 2T P +F P +F N . SHAP explanations To assess the impact of replaceable variables on SHAP explanations, we train Support Vector Regressors (SVR) and XGBoost Regressors (XGB) for each of the 90 synthetic MTS with 100 covariates (denoted as full MTS), using a window of 10 timesteps. We then compute the SHAP values of random holdout samples using the KernelSHAP approximation [LL17] , which we sum along the time axis to obtain the overall importance of variables (per additivity property of SHAP values [KVSF20] ). We consider the average of those values as global explanations [BSC + 21]. We obtain average SHAP values s i 1 ...s i 100 corresponding to the global importance of each variable in the MTS number i. The total importance of each equivalent variable set S becomes s tot = i∈S |s i |. Then, we compute the contribution percentage of the most important variables (until top-14) in the equivalent set |s i /s tot |. For comparison with SHAP explanations over a unique MB, we create reduced MTS where we remove all but one replaceable variable per equivalent set. For instance, if in MTS i variables {1, 4, 16, 78} are equivalent, we remove {4, 16, 78} and compute the reduced importance sr i 1 of the representative 1. We obtain instability results by repeating model training and explanation over different (time-contiguous) data splits of each MTS. Stability is measured as proposed by Nogeira et al. [NSB18] , applied to the top |M 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.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: ChronoEpilogi is detailed in Section 3, where we provide the pseudocode of all core functions along with the critical details that could help a third-party implementation of our algorithm. Reproducibility of our experimental results is facilitated by the description of the executed pipelines as well as by the availability of the code and data used. 
experimental results: 
 Variant FBE has near perfect average causal f1-score (Table 1 ), with stable performance as the number of TS increases 2b). Claim 1: Variant FE has a 70% speedup compared to FBE with only a decrease of causal f1 metric of 0.05 (Table 1 ), due to a less precise equivalence discovery phase leading to lower replaceable f1-score (Fig. 2b ). As the number of TS grows, FE runtime grows slower than FBE (Fig. 2a ). Given its high comparative effectiveness and efficiency, we decided to use FE on real datasets. Claim 2: FBE and FE achieve better results than GroupLasso (GL) w.r.t. causal f1-score, despite similar R 2 scores (Table 1 ). We attribute this difference to GL sensitivity in tuning regularization parameters, compared to ChronoEpilogi thresholds (see Appendix Fig. 5 ). Claim 3: FE uncovers multiple MBs in Electricity (for 3 targets) and Traffic (for 7 targets), PEMS-BAY (3 targets), METR-LA (7 targets). In Solar only one MB is found for all targets. The average size of selected TS variables (Table 2 ) is less than 4% of each MTS. The maximal selected size is below 5% all except Electricity (11%). The solutions have close to identical R 2 in Electricity, and spread over 0.02 around the median for Traffic (Fig. 2c ). Fig. 1 illustrates a compact representation of the solution space for one target of Traffic. Claim 4: On averaging across all targets the multiple solutions produced by ChronoEpilogi-FE have similar predictive performance to the unique solution produced by GroupLasso (GL). We use the total 50 targets to conduct a Friedman-Nemenyi statistical test on method ranking [Her20] , and conclude that GL and FE are not statistically different, but Claim 6: variable selection do improves models' performance compared to models trained on the original MTS (full analysis Appendix E.2). Claim 5: In Electricity and Solar, GL produces slightly smaller solutions than FE while in the rest datasets, FE produces 8 times smaller solutions than GL (Table 2 ). [HRL15] with Dynamic Time Wrapping distances [RGFO17] , but with a quadratic complexity in the number of TS due to distance computations, which are proven to be impractical for large dataset sizes. We selected GroupLasso [NST13] as baseline in our experiments, as it can specify groups of TS [LNZ15] [HMNK10] while applying Lasso-type optimization with linear modeling. It is worth noticing that Lasso-type algorithms have been extended to identify multiple solutions in i.i.d data [PLCT18] without, however, proving any formal property.As a matter of fact, the problem of computing multiple solutions for the variable selection problem is still in its infancy. In particular, we are not aware of any algorithm applied to MTS data. Theoretically sound methods focus on different distributional properties and assumptions. TIE* [SLLA13] assumes that a single Markov boundary discovery algorithm for non-faithful data can be called as an oracle.All MB(T ) are discovered using exhaustive combinatorial conditional tests (O(s.n s ) where s is the maximal size of a considered MB and n the number of variables [BT21] ). TMFBS [BT21] proposes forward-backward MB discovery algorithms that aims to decrease the number of redundant [SLLA13] . In construct, ChronoEpilogi computes provable equivalent subsets of variables in O(ns) conditional tests under composition and interchangeability assumptions, while the forward-backward phase is based on a heuristic. We should mention that there also exist related works that rely on pure associational criteria where heuristics are used without any causal guarantee, especially in the field of gene expression [LLZ10a] [LLZ10b] (see [SLLA13] for an in-depth review). Finally, several works highlight the importance of a structured representation of multiple MBs for interpretability concerns, especially in the presence of an exponential number of MBs of a target of interest [TLT18] [LPL + 23]. ChronoEpilogi is the first algorithm that provides a compact representation of mutually equivalent variables for MTS. Redundant and irreplaceable variables can be easily distinguished at first glance (see Figure 1 ).Question: Does the paper provide open access to the data and code, with sufficient instructions 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 available at: https://anonymous.4open.science/r/ChonoEpilogi/. We include guidelines to run the scripts of our experiments.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] 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. Justification: Anonymity of our submission is guaranteed in all sections. The research proposed in this paper has no impact and potential harmful consequences for the society, nor does it use personal data.Guidelines:• The answer NA means that the authors have not 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., if there is a special consideration due to laws or regulations in their jurisdiction). 
conclusions: 
 Overall, the paper's contributions are (a) the adaptation to time series and further development of the theory of multiple variable selection, (b) the design of the first variable selection algorithm for time-series data, called ChronoEpilogi, that scales to thousands of available series, (c) the conditions of soundness and completeness of ChronoEpilogi, and (d) the empirical evaluations of ChronoEpilogi demonstrating the presence of multiple solutions in real data, achieving on par performance against Group Lasso and no selection, while reducing the number of TS required to build the model by 96% (on average) by conserving or even improving forecasting performance. The reduced model could be employed as the final model for production, or as a surrogate to a model using all available time series (assuming it is better performing) to facilitate interpretation, visualization, and explanations. Finally, we leave as future work the construction of ensemble models trained on several or all MBs as a means to create forecasters more robust to noise, faulty sensors, or other systematic errors.We note that the conclusions are limited to the scope of the experimental study. The latter could benefit from a larger scope of experiments with more real and synthetic datasets of varying size, statistical properties, and data types (e.g., discrete time-series). In addition, different variants of the main algorithm could employ non-linear models to compute residuals and non-linear correlation methods for heuristically selecting time series. The implementation has not been optimized at the low level to reach higher computational gains. Other distributional properties that could lead to greedy yet sound and complete multiple time series selection algorithms could be explored. We also hope to relax the assumption of stationarity in future works, as most practical applications have trends/seasonality, different train-test distributions, or change points. 
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 not been able to detail these important issues.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or 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 not 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, it 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 as 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 or unintentional) misuse of the technology. • If 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). 
safeguards: 
 Question: 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 (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We have not identified any risk in this sense.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: 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?Answer : [Yes] :Justification: We cite all models and third-part implementation we used in the implementation of our solution. Our own assets including code and produced datasets under free and open source license Apache V. 2.0, as we include a modification of a 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 a 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 be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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 not available online, the authors are encouraged to reach out to the asset's creators. 
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: //anonymous.4open.science/r/ChonoEpilogi/, along with comprehensive documentation following the Code and Data Submission Guidelines provided by the NeurIPS conference committee.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of 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 is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
 
 Appendix: Information Equivalence We remark that the equivalence between variables properly defines an equivalence relation. From its definition, it is straightforwardly reflexive (X ≡ T X|Z), symmetric (X ≡ T Y |Z =⇒ Y ≡ T X|Z) and transitive.Independence relationships validate certain properties in any distribution. Contraction and weak union [Pea88] link independences of sets of variables to conditional independences:Lemma A.1 (Weak Union). Let X, Y , Z, W be sets of variables, then X |= Y ∪ W |Z =⇒ X |= Y |Z ∪ W . Lemma A.2 (Contraction). Let X, Y , Z, W be sets of variables, then X |= Y |Z and X |= W |Z ∪ Y =⇒ X |= Y ∪ W |Z.Equivalent variables conditioned on a certain set have the same conditional information on the target.They also contain the same information on the target when considering only information that they have in common with another variable B. Precisely, the following property holds for any probability distribution:Lemma A.3 (Lemiere [Lem07]). Let X, Y, T, B be random variables and Z a conditioning set. If X ≡ T Y |Z, then:EQUATIONJoint normal distributions have a specific independence structure: pairwise independence implies mutual independence in such distributions, which is not the case generally.Lemma A.4. [JW07] Let U , V be two sets of random variables such that U ∪ V is joint normal. Then:U |= V ⇐⇒ ∀U ∈ U , ∀V ∈ V , U |= VWhen the covariance matrix of a joint gaussian distribution is invertible, the distribution is faithful (factorising into a bayesian network). In this circumstance, there is a unique Markov boundary. For equivalence relations to exist in jointly gaussian distributions, the covariance matrix must be singular.Composition Equivalent search algorithms [SLLA13] , [PNBT07] can alternatively rely on a local version of composition holding only for the target variable. Assumption A.1 (Local Composition [SLLA13] ). The local composition assumption with respect to a target T holds when for all X, Y , Z, when Eq.1 holds with T = T . The typical example of a XOR operation on two Bernouilli variables X, Y with probability 0.5, and T = XxorY invalidates composition. Indeed, X and Y individually bring no information on T , while together define entirely T .Another way of seeing composition is regarding interaction among variables. In clinical trials, it is often the case that the causal effect of two context variables X, Y is greater (or lesser) than the sum of the individual effects of X and Y . This is called interaction. Composition asks that two variables X and Y that are individually independent from T , cannot interact such that T becomes dependent of them. Note that interaction is possible in the case of any of X and Y already being dependent on T .Interchangeability Interchangeability asks for the set of multiple solution M to factorize as a cartesian product of equivalence class of individual variables. This allows the solution space to be quickly discovered by a linear time algorithm.In general, equivalences can occur at the level of group of variables, and not just individual variables. For instance, let Z be a standard normal random variable. Define X = I[Z > 0].Z and Y = I[Z < 0].Z. Since X has the value of Z if Z > 0 and Y the converse, then Z = X + Y . In this situation, for any target, Z ≡ T X, Y . Interchangeability requires that Z cannot be part of a markov boundary of T .It is possible to relax Interchangeability to a looser version by introducing additional complexity to the algorithm. If instead of testing singletons for equivalences, the algorithm also tests pairs of variables, more complex equivalence relations can be integrated. To take into account all possible equivalence relations among group of variables, exponential search is required [SLLA13] .Algorithm 7 details the Pearson correlation based test. We note that in this version, a p-value is obtained for each of the lags, the minimum of which is returned. This procedure returns p-values that are likely lower than the actual p-value, due to using the minimum as an estimator. Hence, testing for previous and current model statistical difference is kept separate from the correlation routine, ensuring that Forward phase termination relies on a single p-value.Algorithm 7 Lag-Pearson-pvalRequire: residuals R, candidate variable X ′ ∈ R M , maximum lag L 1: p ← 1 2: for 1 ≤ l ≤ L do 3: p ← min(p, Pearson-corr(R, X ′ l,M -L+l-1 )) return pThe equivalence search (Algorithm 8) of the approximate FE version builds surrogate models g Θ to predict the residuals instead of the target T . Its structure correspond to the definition of Information Equivalence, in that it tests proxies for the two assertions X |= Res|S new , and S new |= Res|X. Establishing equivalence for residuals is not always guaranteed to ascertain equivalence for T (see Theorems), but we empirically find it a good proxy for our synthetic datasets.Require: residuals Res, selected variable S new , remaining variables R, maximum lag L, equivalence threshold δ FBE Soundness (Thm.D.3). In the following, we define a useful notation S k where S is an indexed set from 1 to n and 1 ≤ k ≤ n, as the set of elements of S indexed from 1 to k.1: set E ← {} 2: for X ∈ R do 3: p ← pvalue of Likelihood ratio test between Res = g Θ (X ∪ S new ) and Res = g Θ (S new ) 4: if p ≥ δ then ▷ XFirst, we remark that it is enough to prove that M B n ≡ T S n , where M B n is one of the produced solutions M B 1 , ..., M B n ∈ M. In that case, since S n is a Markov blanket, for any variable of the M B k |= T t |S k , Z S k |= T t |M B k , ZThe first condition is true since S n is a Markov Blanket of T t . To prove the second condition, by the property of decomposition, and by the assumption of composition (Assumption 2.1), it is necessary and sufficient to prove thatEQUATIONBy using Lemma A.3 with equivalence Corollary. All returned solutions are Markov Blankets (Thm.D.3), and at least one is a Markov Boundary (Thm.D.2). Since by definition of Markov Boundaries, any Markov Blanket is the superset of at least one Markov Boundary, then all returned solutions must be of the same size as the Markov Boundary contained. Hence, all solutions are Markov Boundaries.M B k-1 ≡ T S k-1 |S k , Z (recursion supposition), equation (3) is equivalent to S k |= T t |S n -k ∪{M B k }Completeness (Thm.3.1). Direction =⇒ : suppose that the equivalence search starting from a selected set S find out all Markov boundaries of T t . From the proof of Thm D.3, we can deduce that for whichever indexing of the reference set S, denoting S 1 ....S n , and for whichever equivalent boundary M produced by the algorithm, there is an ordering M There is a permutation σ : [|1, n|] → [|1, n|] such that f • σ is monotonically decreasing. We reindex S ′ i = S σ -1 (i) and similarily for M ′ i . Since Eq.( 5) holds for any common reindexing of S and M , it holds for S ′ and M ′ . Also, there is an integer k for which the setEQUATIONM ′ k ∪ (S ′ \ S ′ k ) is identical to A. By Thm.D.3, M ′ k ∪ (S ′ \ S ′ k) is a Markov boundary. Direction ⇐=: suppose that for any correctly indexed M , any function f : [|1, n|] → 0, 1 induces a set A and A is a Markov boundary of T t . In particular, for any k ∈ [|1, n|], the function f k : i → 1 if k = i else 0 induces a Markov boundary S -k ∪ M k . Therefore, we deduce that S k |= T t |S -k ∪ M k for all k,We optimize the hyperparameters for each algorithm and forecaster we consider (see Table 3 ). We use the python library optuna with Grid Search optimization. ChronoEpilogi thresholds cover several orders of value between 10 -20 and 0.05. As we observed that for GroupLasso regularization parameter 10 -20 , all TS where select no matter the dataset, and at 0.1, no TS was selected other than the target past, GroupLasso group regularization parameter ranges 25 values within this range on a logarithmic scale.We describe the pipeline used for tuning and evaluation in Fig. 3 . Cross-validation on time series is done in a fixed start window growing along time. The holdout set size is consistently 1000 for synthetic MTS, and 30 percent of the MTS for real dataset. Train best configuration on Tuning (TSS algorithm and Forecaster)Find the equivalent markov blanketsFigure 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.Tuned hyperparameters GroupLasso [NST13] group regularization coefficient ARDL [PSS01] Trend (no trend, constant, linear) SVR Kernel (rbf, sigmoid), regularization weight DeepAR [SFGJ20] Layer size, #layers, dropout rate TFT [LALP21] Layer size, #layers, dropout rate, #attention heads cover at least 40 percents of the Tuning split, so to obtain meaningful results and retain enough samples.Over Fig. 4 , we examine how metrics vary when the thresholds are chosen to maximize the R 2 of the forward-backward solution S. We observe that the forward and backward threshold depend on the R 2 value, hence the proportion useful information variation compared to the total variation of the data. We confirm that different thresholds do not result in vastly different performances.The standard deviation of the R 2 of the top 5 configurations for each MTS on average is 0.001 for FBE and FE compared to 0.02 for GroupLasso. See 5 for the boxplot of the distribution. This indicate a much higher impact of hyperparameter tuning on GroupLasso than ChronoEpilogi, and participates to the difficulty of tuning GroupLasso. Practictioners should expect spending more time tuning regularization-based solution than p-value-based solutions.The statistical analysis was conducted for 3 populations with 50 paired samples. The family-wise significance level of the tests is alpha=0.050. We rejected the null hypothesis that the population is normal for the populations ChronoEpilogi (p=0.000), GroupLasso (p=0.000), and NoSelection (p=0.000). Therefore, we assume that not all populations are normal. Because we have more than two populations and some of them are not normal, we use the non-parametric Friedman test as omnibus test to determine if there are any significant differences between the median values of the populations. When more than 20 solutions are found, we randomly select 20 of them for evaluation. We observe that both the number of targets for which there are multiple solutions decreases, with a diminution of the number of solutions overall.We use the post-hoc Nemenyi test to infer which differences are significant. We report the median (MD), the median absolute deviation (MAD) and the mean rank (MR) among all populations over the samples. Differences between populations are significant, if the difference of the mean rank is greater than the critical distance CD=0.469 of the Nemenyi test. We reject the null hypothesis (p=0.000) of the Friedman test that there is no difference in the central tendency of the populations ChronoEpilogi (MD=0.289+-0.112, MAD=0.105, MR=2.580), GroupLasso (MD=0.287+-0.106, MAD=0.106, MR=2.320), and NoSelection (MD=0.442+-0.310, MAD=0.257, MR=1.100). Therefore, we assume that there is a statistically significant difference between the median values of the populations. Based on the post-hoc Nemenyi test, we assume that there are no significant differences within the following groups: ChronoEpilogi and GroupLasso. All other differences are significant.For comparison purpose, we ran the full version of the algorithm on two of the real datasets. The performance is reported in Table 4 . The number of target for which multiple solutions are found is reduced to 2 in METR-LA and 6 on Traffic (Figure 8 ). We observe that overall, results from synthetic datasets are verified: FBE decreases selected size and diminishes the total number of solutions, with higher computational time, resulting in an increase in predictive performance. We confirm that version FE tend to select false positives, but keep the constatation that mutliple solution exist in real datasets.A surprising experimental result was that at the end of parameter tuning, SVRModels were outperforming TFT models during cross validation. As we suspected that cross validation could influence negatively the TFT performance (earlier split having less data, and Transformer-based models notoriously requiring consequent datasets), we tried a different pipeline where tuning and testing was done on a single training/validation/test split with 70,20,10 repartition. The result are in Table 5 . We note that no matter the predictive model, ChronoEpilogi does not return multiple solutions on this dataset. Additionally, the results using the modified pipeline do not change the method ranking of ChronoEpilogi, GroupLasso and NoSelection.It is known that Shapley values currently suffer from inclusion of unrealistic data instances when features are correlated, the explanations even of linear models can be misleading [AJL21] . In this section, we investigate how SHAP explanations of regression models might misrepresent the role of variables belonging to some but not all Markov boundaries of the modeled target.To assess the impact of replaceable variables on SHAP explanations, we tune and train Support Vector Regressors (SVR) and XGBoost Regressors (XGB) for each of the 90 synthetic MTS with 100 covariates (denoted as full MTS), using a window of 10 timesteps. We then compute the SHAP values of random holdout samples using the KernelSHAP approximation [LL17] , that we sum along the time axis to obtain the overall importance of variables (per additivity property of SHAP values [KVSF20] ). We consider the average of those values as global explanations [BSC + 21]. We obtain average SHAP values s i 1 ...s i 100 corresponding to the global importance of each variable in the MTS number i.To evaluate the total importance for each set of equivalent variables, without the disturbance brought by the presence of replaceable variables, we apply the same average SHAP computations as previously to the 90 MTS, where we remove all but one replaceable variable per equivalent set (denoted as reduced MTS), to compare with SHAP explanation of unique Markov processes. For instance, if in MTS i variables {1, 4, 16, 78} are equivalent, we remove {4, 16, 78} and compute the reduced importance sr i 1 of the representative 1. Note that variable 1 is irreplaceable in the reduced MTS.For each equivalent set {j 1 , ..., j c } in MTS i, we compute the contribution of each variable j ∈ {j 1 , ..., j c } to the total absolute importance of the set, as|s i j | k |s i j k| . We plot the contribution of the most important variable, second most important, etc (Fig. 9 ). We choose to use absolute importance values as averages are not easily comparable with relative importance. We observe that importance is not equally shared among equivalent variables. Models tend to select one or a few variables as important predictors. Additionally, the importance of the top variable relative to the other diminishes as the number of equivalent variables grows. We include in the repository a requirements file requirements.txt listing the necessary dependencies.Table 7 : Presence of replaceable variable increases the instability of top important variables.For each synthetic MTS with 100 variables, all variables but one in each equivalent sets are removed to make a reduced MTS with a unique Markov Boundary. XGB models are computed and SHAP explanations produced. We compute the stability [NSB18] of the |M B(T )| most important variables in the entire explanation, for the original MTS (with multiple MB(T)) and the reduced MTS (with a unique MB(T)). Note that the metric adjusts for total number of variables in the dataset [NSB18] . We observe that explanations are significantly more unstable when replaceable variables are present. Justification: At the beginning of Section 5, we added a detailed list of the claims made in the abstract and introduced in the last part of the introduction. Here, we support each claim by providing all the empirical and statistical evidences we have observed.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to 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 a 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 submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe 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., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.Chronos greek and Epilogi greekhttps://github.com/ev07/ChronoEpilogiConditions leading to faithfulness violations have been studied from the standpoint of probability distributions and graphs[Sad17,ZZM17] as well as applications in homeostatic systems[ZS08], evolutionary systems[And13], gene expression[SLLA13] or spatio-temporal records of meteorological phenomenons [YWD + 17].Alternatively, algorithms might rely on an oracle producing variable sets with given properties[SLLA13]. 


Decision: 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.

Reviews: review: paper_summary: paper_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. 
The 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. 
The 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. 
The 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. 
The 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. 
This work would add some valuable knowledge and insights to the community with interest in causal modeling and interpretable learning in time series data.
strength_weakness: strength_weakness: The performance of the proposed methods against the baseline of Group Lasso on real world data sets is not exactly compelling. 
More clarify on the relative advantage of the proposed method(s) would be valuable. 
The 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.
questions: 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: limitations: 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.
paper_summary: paper_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."
strength_weakness: strength_weakness: 1. The proposed algorithm was only compared against GroupLasso and not against any other among the related work mentioned in the paper.
questions: questions: 1. Line 176 should say "forward" rather than "backward."
2. I suspect that algorithm 3 step 9 and algorithm 4 step 3 should have $\geq$ in place of $\leq$.
3. In line 260, how is TSS defined?
limitations: limitations: Yes
paper_summary: paper_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:

1. **Theory Development**: Introduces the problem of multiple time-series variable selection (MTVS) and the concepts of informational equivalence and interchangeability in TS data.
2. **Algorithm Design**: Proposes ChronoEpilogi, which identifies all Markov Boundaries under broad, non-parametric assumptions.
3. **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.
strength_weakness: strength_weakness: While being a very interesting paper, there are some avenues for improvement: 
1. While the authors discuss the scalability, the algorithm's computational complexity seems to be high for very large datasets, potentially limiting its practicality.
2. 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
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
questions: questions: It would be great if the authors can discuss about the computational complexity and the limitations stemming from their assumptions
limitations: limitations: See weakness above
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.;

- 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;

- 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.
questions: questions: - Why did you choose ARDL over other models like RNNs which are also autoregressive models and are certainly more efficient (Equation 2)?

- 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)?

- What can explain the outperformance of SVR over TFT and DeepAR with the Solar Dataset?

- How do the authors expect their algorithm to perform when faced with sparse multivariate time series, given that missing value can affect causal inference?

- 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: limitations: - The authors have identified the limitations of their work and plan to address them in their future work;

- On line 176 it should be forward instead of backward;

-  The reviewer suggests improving definition 2, which is the core definition of the paper.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

2. 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. 

3. Prior work discussion is incomplete: for instance, signature paths and feature selection should be discussed and compared to the method. Some specific examples include
    - Cross-correlation analysis,
    - Signature transforms https://arxiv.org/abs/1603.03788
questions: 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?
2. 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.
3. Please comment on weakness #1, which I have listed above.


## Additional comments:

Additional comments
I am really surprised that the primary area is interpretability and explainability, even though the paper almost does not discuss these aspects of the solution.

What is V in tvs?

“under Composition, Interchangeability, and other broad, non-parametric assumptions” – could you list all the assumptions here?

L102 : conditional independence is not written well

L205 RVS is not properly capitalized

Eq 2: what does the dot mean? Is it a misprint?

Table 2:
What do you mean by size?
No selection yields worse results. I wonder, what happens if you use a stronger base model for forecasting, e.g., LSTM?

Re: Statistical significance:
The paper does not report statistical significance in Table 2.
limitations: 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.

Comments: title: Additional Contributions
Acknowledging author response
Response regarding additional evaluation
comment: Dear reviewer, thank you a lot for the time you have spent on our paper and for the constructive feedback. We believe this work will be valuable for researchers of our domain.
Thank you for responding! Your clarifications helped me much better understand the contribution. I will raise my score to 5.
Dear authors, thank you for your detailed answers. My concerns were clarified. I keep the score unchanged.
We thank you for your positive feedback on our answer during the rebuttal. 

We will definitively further detail domain discovery under the light of the Data Explainability in the paper introduction, add new claims regarding the limitation of Feature Importance methods in the presence of redundant features and include in the main paper and the appendix the experiments that support them.

We should stress that that this a fourth contribution of our work that could not been made without the previous three : 

(a) ChronoEpilogi algorithm variants 

(b) Theoretical guarantees of the algorithm 

(c) quantitative evaluation of ChronoEpilogi results with synthetic and real datasets. 

More precisely, we will make the following new claims, based on the SHAP value experiments we have included in the rebuttal:

Claim 7: The importance score of variable with equivalences decreases by half in our experiments, as the number of equivalent variables grows. 

Claim 8: the top |MB(T)| predictors in explanations are significantly more unstable to data changes in presence of multiple Markov Boundaries compared to distributions with only one MB(T) (pvalue 4e-13).

The above claims can be validated from the one pdf  page we added in our global answer. The experimentation protocol (models libraries, SHAP explainer implementation, tuning, instability metric) we used will be also detailed. 

We believe that, on the additional page of the camera-ready version we will be able to include the main results of our new experiments in the paper beside the appendix.
Thank you for the response!

The response clarifies my concern. I think this is a very interesting aspect of your work about data explainability. In my opinion, viewing this work from this perspective would make it a really valuable contribution.

Currently, my concern is that adding such a massive claim would make the paper significantly different from the version under review. Could you briefly list all the changes you plan to incorporate?
Thanks for the detailed response
Thank you for elaborating on your on-going efforts on additional empirical evaluation. It would be good (and strengthen the paper)  to include the results of those (regarding both Group Lasso and SHAP) in an updated version of the paper.