for each
S[FOR]ForEach[1] #1

Machine learning algorithms are designed to capture complex relationships between features. In this context, the high dimensionality of data often results in poor model performance, with the risk of overfitting. Feature selection, the process of selecting a subset of relevant and non-redundant features, is, therefore, an essential step to mitigate these issues. However, classical feature selection approaches do not inspect the causal relationship between selected features and target, which can lead to misleading results in real-world applications. Causal discovery, instead, aims to identify causal relationships between features with observational data. In this paper, we propose a novel methodology at the intersection between feature selection and causal discovery, focusing on time series. We introduce a new causal feature selection approach that relies on the forward and backward feature selection procedures and leverages transfer entropy to estimate the causal flow of information from the features to the target in time series. Our approach enables the selection of features not only in terms of mere model performance but also captures the causal information flow.

In this context, we provide theoretical guarantees on the regression and classification errors for both the exact and the finite-sample cases. Finally, we present numerical validations on synthetic and real-world regression problems, showing results competitive w.r.t. the considered baselines.

§ INTRODUCTION

Machine learning <cit.> techniques are designed for predictive modelling, data analysis, and decision-making and are applied in numerous fields such as medicine, finance, and social sciences. One of the most significant challenges in ML is dealing with observational data, which is collected without controlling for the variables of interest or the confounding variables. This type of data can lead to biased results, hindering the generalizability and interpretability of the ML models. To reduce the number of dimensions and avoid memory and modelling issues like the curse of dimensionality, ML models are trained over a set of variables identified by feature selection [Li et al., 2017], which usually chooses the features with the largest association with the target, without inspecting the causal relationships  [Peters, Janzing, and Schölkopf, 2017].

Causal discovery is the process of identifying causal relationships among variables from observational data without external interventions [Spirtes et al., 2000, Eberhardt, 2017]. Its primary purpose is to undercover the direction of causality and estimate the magnitude and uncertainty of causal effects. This approach not only helps in developing more accurate and robust models but also provides a deeper understanding of the underlying mechanisms governing the data [Pearl, 2009].

The importance of observational causal discovery for time series data in ML lies on the fact that a model trained on causal features enables researchers to predict the outcomes of an intervention or a policy change in a robust and interpretable way [Schölkopf, 2022]. This is particularly useful in situations where it is not feasible or ethical to intervene, such as in healthcare, social sciences, or economics [Imbens and Rubin, 2015]. Furthermore, considering a physical system such as climate science, the causal knowledge of direct and indirect effects among variables can help to understand the physical system, filtering the effect of spurious correlations. Causal discovery can also guide the choice of the features by performing a causal feature selection that identifies those with the most relevant flow of information to the target. However, causal discovery for time series is challenging on complex systems due to unobserved confounders, high dimensionality, autocorrelated variables, and time lags [Runge et al., 2019]. Moreover, causal discovery algorithms usually focus on theoretical guarantees on the (asymptotic) convergence to the real causal graph, or to an equivalence class of graphs, without inspecting the effect of the selection of the candidate causes of a variable on its prediction.

In this paper, we propose a novel methodology at the intersection between feature selection and causal discovery. We start in Section <ref> by reviewing the main approaches and assumptions of causal discovery for time series. Then, in Section <ref>, we introduce the notation and formulation of the problem. In Section <ref>, we introduce the novel causal feature selection methodology that relies on the forward and backward feature selection procedures and leverages the transfer entropy (TE) [Schreiber, 2000] as building block to guide the selection of the features. TE is a causal quantity specifically designed to estimate the flow of information between features and the target in time series.

Unlike traditional feature selection and similar to Granger causality, our causal approach exploits the property of TE to filter out the autoregressive component of the target in order to measure the asymmetric flow of information from a feature to it without the confounding association due to the target itself. Furthermore, conditioning on the other features also makes it possible to filter out the information due to the presence of a causal variable that a feature and the target have in common.
In this setting, we provide theoretical guarantees about the regression and classification error due to the reduction of features. In addition, in Section <ref>, we analyse the finite-sample scenario, exploiting a concentration bound of a specific kernel-based TE estimator available in the literature [Singh and Póczos, 2014].  Finally, in Section <ref>, we provide numerical simulations in regression settings, showing the capability of the proposed approaches to identify causal features in synthetic experiments and to achieve competitive performance on real-world datasets.

§ RELATED WORKS

In this section, we review state-of-the-art causal discovery approaches, focusing on their application to time series. [Regarding feature selection methods, we refer the interested reader to the survey Li et al., 2017 for an overview.]

Classical Causal DiscoveryClassical causal discovery relies on four main assumptions [Pearl et al., 2000].

Causal Markov and causal faithfulness ensure a correspondence between the (conditional) probability independence between variables and the absence of a causal relationship in the causal graph. Then, acyclicity avoids the presence of cyclical causal relationships, and causal sufficiency guarantees the absence of latent confounders. The benchmark non-parametric algorithm of causal discovery that assumes all these four assumptions to hold is the PC algorithm [Spirtes et al., 2000]. Non-parametric approaches that relax parts of the assumptions or parametric approaches that assume a linear or non-linear structural causal model are revised in [Eberhardt, 2017, Glymour, Zhang, and Spirtes, 2019].

Causal Discovery for Time Series
In  [Runge et al., 2019], the difficulties of extending the causal discovery framework to high-dimensional time series are extensively analysed. The classic extension of causal discovery to time series is the notion of Granger causality [Granger, 1969], which evaluates the impact of the past value of a feature on the prediction of the target as its causal strength [Ding, Chen, and Bressler, 2006]. This approach may suffer from the curse of dimensionality and autocorrelation effects, as highlighted in [Runge et al., 2019], where the PCMCI algorithm has been proposed to deal with these effects. Variations of Granger causality-based methods and of the PCMCI algorithm that deals with contemporaneous causal effects and latent confounders can be found in [Moneta et al., 2011, Malinsky and Spirtes, 2018, Runge, 2020, Gerhardus and Runge, 2020]. Other existing approaches assume specific structures of the causal models or apply neural networks to extend the methods to highly non-linear contexts. A broader overview can be found in Moraffah et al., 2021.

Similar to the approaches proposed in this paper, some recent works  [Mastakouri and Schölkopf, 2020, Mastakouri, Schölkopf, and Janzing, 2021, Dong and Kluger, 2022] focus on identifying sets of causally relevant features rather than the complete causal graph, although they are not directly comparable because of different settings and assumptions. In particular,  [Mastakouri and Schölkopf, 2020, Mastakouri, Schölkopf, and Janzing, 2021] introduce and apply the SyPI algorithm, which is guaranteed to asymptotically identify (some) ancestors of a target variable, also in the presence of latent confounders, under several assumptions on the relationships between variables. [Dong and Kluger, 2022] focus on the identification of a set of causal features, exploiting a neural network to optimize a specifically designed loss to asymptotically guarantee the identification of a set of features causally interconnected.

Causal Information Measures
Several indexes from information theory are available to quantify the causal information that flows from one variable to another.
Transfer entropy <cit.> evaluates the information flow from one variable to another, computing the information that the past values of the first have on the second. Similarly, directed information [Massey et al., 1990] computes the cumulative sum of information shared between all the past values of the first variable and the actual value of the second one. These two quantities are strictly related [Liu and Aviyente, 2012], and they can be applied to test for Granger causality between variables [Amblard and Michel, 2012].

§ PRELIMINARIES

Time Series
We consider a supervised learning problem with $D$ features and a scalar target. At the time $t$, we denote with $Y_t$ the target's current random value and with $X_i^t$ the current random value of the $i$-th feature. Features and target are multivariate time series $X_i=(X_i^{t},X_i^{t-1},\dots)$ for $i\in [D] \coloneqq \{1,\dots,D\}$ and $Y=(Y_{t},Y_{t-1},\dots)$ that form a nonlinear, discrete-time, stationary vector autoregressive process. We denote with $X_A^{\xtpast}$ the last $L$ values of the set of features having indices in $A\subseteq [D]$ and with $X_{A}^{\xtpast}$ the last $L$ of the set of features having indices in set $A$ (for simplicity, $X^{\xtpast} = X_{[D]}^{\xtpast}$ and $X_i^{\xtpast} = X_{\{i\}}^{\xtpast}$). We also denote $X_A\cup X_i := X_{A\cup \{i\}}$ and $X_A\setminus X_i := X_{A\setminus \{i\}}$. Finally, $Y_{\ytpast}$ is the set of the last $M$ values of the target. We assume that the absolute value of the target is bounded by a constant $B$ uniformly over time. We denote with $\mathbb{E}_Z[f(Z)]$ and $\mathbb{V}\mathrm{ar}_Z[f(Z)]$ the expected value and the variance of a function $f(\cdot)$ of the random variable $Z$  its distribution.

Causal Model
Together with the four classical assumptions of causal discovery introduced in the previous section (causal sufficiency, acyclicity, causal Markov assumption, faithfulness), we assume that the actual target may only depend on its last $M$ values and the last $L$ values of the features.

The (stationary) autoregressive processes are described by the following structural causal model (SCM):

\begin{align*}
&X_i^t = f_i(\text{pa}(X_i^t),\epsilon_i), \quad i\in [D], \qquad
Y_t = f(\text{pa}(Y_t),\epsilon),
\end{align*}

where $f_i, f$ control the relationship between each variable and its causes, $\epsilon_i,\epsilon$ are independent unobserved noise terms and $\text{pa}(X_i^t)$ (resp. $\text{pa}(Y_t)$) is the set of unknown causes of feature $X_i$ (resp. the target $Y_t$), assumed to depend only on some past values of the considered features (resp. target).

Information Theory
Mutual information <cit.> is an index considered in information theory to evaluate the amount of information that a feature $X$ shares with the target $Y$, defined as:

\begin{align*}
I(X;Y)&\coloneqq \mathbb{E}_{X,Y}\left[\log \frac{p(X,Y)}{p(X)p(Y)}\right],  %\\&\ = \int\int p(x,y) \log \frac{p(x,y)}{p(x)p(y)}\, \mathrm{d}x \mathrm{d}y,
\end{align*}

where $p(\cdot,\cdot)$ is the joint p.d.f. of the random variables $(X,Y)$. Conditional mutual information (CMI) is obtained conditioning w.r.t. a third random variable $Z$, $I(X;Y|Z)\coloneqq\E_Z[I(X|Z;Y|Z)]$.

Transfer entropy <cit.> evaluates the information flow from one variable to another, computing the information that the past values of the first have on the second. This is the directional counterpart of mutual information, which is symmetric and, therefore, cannot identify the causal direction of the flow of information.
Similarly to CMI, conditional transfer entropy [Shahsavari Baboukani et al., 2020] extends the definition of TE to a conditional set that only considers the additional amount of information flow from $X$ to $Y$ given the flow from $Z$ to $Y$. We denote the TE between $X$ and $Y$ as $TE_{X\rightarrow Y}:=I(Y_t;X^{\xtpast}|Y_{\ytpast})$ and the conditional transfer entropy between $X$ and $Y$ given $Z$ as $TE_{X\rightarrow Y| Z}:=I(Y_t;X^{\xtpast}|Y_{\ytpast},Z^{\xtpast})$. $M$ and $L$ determine the maximum number of previous values of the features and target that impact the current value of the target.

§ TRANSFER ENTROPY FOR CAUSAL FEATURE SELECTION

In this section, we introduce two theoretical results that relate the effect of selecting features based on the TE to the regression and classification errors (Section <ref>). Then, we apply these findings to devise novel forward and backward causal feature selection algorithms (Section <ref>). The results we present extend the findings of [Beraha et al., 2019] in two directions: (i) we move from an i.i.d. supervised learning setting to a time series setting, (ii) we provide a result that is expressed in terms of TE, which represents a natural expression of the loss of causal information for time series.

§.§ Regression and Classification Error Bounds

We first provide an upper bound to the Mean Squared Error (MSE) of the regression performed to estimate the value of the target $Y_t$ at time $t$, considering the last $M$ values of the target and the last $L$ values of the selected features.[Note that the bound holds for a generic subset of features, not necessarily the result of a feature selection.]  The proofs of the results are reported in Appendix <ref>.

Let $ \Bar{A} \subseteq [D]$ be a set of indices and let $A = [D] \setminus \Bar{A}$ be its complementary.

Then, considering a lag of $L,M \in \mathbb{N}_{\ge 0}$ for the features and the target respectively, the regression error suffered by considering only the set of features $X_{\Bar{A}}$ is bounded by:

\begin{equation}\label{eq:regrErr}
\begin{aligned}
\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}[(Y_t-g(X_{\abar}^{\xtpast},Y_{\ytpast}))^2] \\ \leq \sigma^2 + 2B^2\cdot TE_{X_{A}\rightarrow Y| X_{\abar}},
\end{aligned}
\end{equation}

where $\mathcal{G}_{\abar}$ is the set of all functions that map $X_{\abar}$ to $\mathbb{R}$, 8.59$\sigma^2 = \E_{X,Y}[(Y_t-\E[Y_t|X^{\xtpast},Y_{\ytpast})^2]$ is the irreducible error, and $B$ is the maximum absolute value of the target.

The result of the theorem complies with the intuition that, considering a subset of features, the expected MSE is bounded by the irreducible expected MSE (unavoidable even when considering the entire set of features), plus an index of the information flow from the discarded features to the target $TE_{X_{A}\rightarrow Y| X_{\abar}}$.
A similar result holds for the Bayes error in classification settings.

Let $A$ and $\abar$ be defined as in Theorem <ref>. Then, considering a lag of $L,M \in \mathbb{N}_{\ge 0}$ for the features and the target respectively, the Bayes error suffered by considering only the set of features $X_{\Bar{A}}$ is bounded by:

\begin{equation}\label{eq:classErr}
\begin{aligned}
\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}\left[\mathbb{1}_{Y_t\neq g\left(X_{\abar}^{\xtpast},Y_{\ytpast}\right)}\right] \leq \epsilon + \sqrt{2\cdot TE_{X_{A}\rightarrow Y| X_{\abar}}},
\end{aligned}
\end{equation}

where 8.59$\epsilon = \E_{X,Y}\left[\mathbb{1}_{Y_t\neq \arg\max_{y_t\in\mathcal{Y}}p(y_t|X^{\xtpast},Y_{\ytpast})}\right]$ is the irreducible Bayes error and $\mathbb{1}$ is the indicator function.

The presented results enable the derivation of algorithms based on conditional TE to perform causal feature selection in a sound way. Indeed, we can easily control the loss in prediction performance via the conditional transfer entropy of the discarded features, for both regression and classification.

§.§ Algorithms

Building on the presented results, in this section, we derive Backward and Forward Transfer Entropy Feature Selection (TEFS) for causal feature selection via conditional TE. For both of them, we provide guarantees on the regression/classification error at the end of the selection process.

Backward Approach

Algorithm <ref> reports the pseudo-code of the Backward TEFS algorithm, defined in terms of the feature and target lags $M$ and $L$ and in terms of the maximum information loss $\delta$.

The algorithm begins by assuming that all features are significant to the target variable (line <ref>). Then, it iteratively selects the feature with the smallest information flow, given the remaining features (line <ref>), and removes it if the loss of information flow keeps the error below $\delta$ (lines <ref>-<ref>). When removing a feature leads to a cumulative information loss greater than the maximum desired (line <ref>), Backward TEFS stops and returns the remaining features (line <ref>). These are the candidate causal drivers of the current value of the target in terms of conditional TE.

The two hyperparameters $M$ and $L$ determine the number of past time steps that the algorithm is considering to evaluate the importance of the features on the current target and filter out the autoregressive component. They balance the positivity-unconfoundedness tradeoff (and, more generally, the bias-variance tradeoff in ML). Indeed, considering a larger lag for a covariate may lead to a better understanding of its impact on the actual target and increases the adherence to the causal sufficiency assumption.

On the other hand, given a time series of length $T$, increasing the number of past values of the features (and of the target)
reduces the number of available samples (at the limit case, considering $L=M=T$ we only have one sample with the last observed target $Y_T$ and all past features and target to condition on). This may lead to overfitting and the curse of dimensionality [D'Amour et al., 2020]. Finally, the hyperparameter $\delta$ is the maximum information loss, in terms of TE, that the user is willing to lose in order to reduce the number of features, as theoretically justified by the next theorem.

Backward TEFS: Backward Transfer Entropy Feature Selection
[1]
$D$ features with $L$ past values, $\{X_1^{\xtpast},\dots,X_D^{\xtpast}\}$; target $Y_t$ and its $M$ past values $Y_{\ytpast}$; $N$ samples; maximum information loss $\delta$
set of selected features $X_{\abar}$

$X_{{A}} \leftarrow \{\}$Initialization: no removed features
$X_{{\abar}} \leftarrow \{X_1^{\xtpast},\dots,X_D^{\xtpast}\}$

$TE\_loss\leq \frac{\delta}{2B^2}$

$TE\_loss\leq \frac{\delta^2}{2}$ for classification
$X_{\hat{i}}\leftarrow\arg\min_{X_i\in X_{\abar}} TE_{X_{i}\rightarrow Y|X_{\abar}\setminus X_{i}} $
$TE\_loss \leftarrow TE\_loss+TE_{X_{\hat{i}}\rightarrow Y|X_{\abar}\setminus X_{\hat{i}}}$
$X_{{A}} \leftarrow X_{{A}} \cup X_{\hat{i}}$
$X_{{\abar}} \leftarrow X_{{\abar}} \setminus X_{\hat{i}}$

By definition of TE, the two algorithms evaluate the contribution of the last $L$ values $X_i^{\xtpast}$ of any feature $i$ as a whole. One may consider each value $X_i^{t-\tau},\tau\in\{1,\dots,L\}$ of each feature at a different temporal lag as a different feature and if needed remove each one separately. Although this would be a forcing of the concept of TE, it would lead to a more accurate identification of the causal interaction between the features and the target, since each of the remaining features would identify a specific feature and time lag, similar to the PCMCI approach. However, this would lead to a larger feature space, an increased computational complexity and a larger risk of the curse of dimensionality. For this reason, in line with the definition of conditional TE and Granger causality, the proposed approach considers the contribution of a specific feature at different time lags as a whole.

Forward TEFS: Forward Transfer Entropy Feature Selection
[1]
$D$ features with $L$ past values, $\{X_1^{\xtpast},\dots,X_D^{\xtpast}\}$; target $Y_t$ and its $M$ past values $Y_{\ytpast}$; $N$ samples; minimum information gain $\Delta$
set of selected features $X_{A}$

$X_{{A}} \leftarrow \{\}$Initialization: no selected features
$X_{{\abar}} \leftarrow \{X_1^{\xtpast},\dots,X_D^{\xtpast}\}$

$TE\_cumulated\leq \frac{\Delta}{2B^2}$

$TE\_cumulated\leq \frac{\Delta^2}{2}$ for classification
$X_{\hat{i}}\leftarrow\arg\max_{X_i\in X_{\abar}} TE_{X_{i}\rightarrow Y|X_{A}} $
$TE\_cumulated \leftarrow TE\_cumulated+TE_{X_{\hat{i}}\rightarrow Y|X_{A}}$
$X_{{A}} \leftarrow X_{{A}} \cup X_{\hat{i}}$
$X_{{\abar}} \leftarrow X_{{\abar}} \setminus X_{\hat{i}}$

The following lemma proves that the overall information loss by removing the set of features $X_A$ can be computed as the cumulative transfer entropy loss evaluated at each step.

Let $T(X_A)$ be the conditional TE between the target and the discarded features $X_A$, given its complementary $X_{\abar}$

$(T(X_A) \coloneqq TE_{X_{A}\rightarrow Y|X_{\abar}})$.

Let also $T_{\abar}=\min_{X_i\in X_{\abar}} TE_{X_{i}\rightarrow Y|X_{\abar}\setminus X_{i}}$ be the minimum conditional TE of a single feature among the selected ones and the target, given the other selected features, and $X_{\hat{i}}=\arg\min_{X_i\in X_{\abar}} TE_{X_{i}\rightarrow Y|X_{\abar}\setminus X_{i}}$ its argument. Then it holds:

\begin{equation}\label{eq:recursiveBackward}
T(X_A\cup X_{\hat{i}}) = T(X_A) + T_{\abar}.
\end{equation}

Lemma <ref> can be exploited to justify the choice of the stopping condition shown in Algorithm <ref>. Indeed the following theorem shows that, by considering the proposed stopping criterion, the increase of the ideal regression (or classification) error is bounded by the maximum information loss $\delta$.

The regression error suffered considering the variables $X_{\abar}$ selected by Backward TEFS is bounded by:

\begin{equation*}
\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}[(Y_t-g(X_{\abar}^{\xtpast},Y_{\ytpast}))^2] \leq \sigma^2+\delta.
\end{equation*}

Similarly, the classification error is bounded by:

\begin{equation*}
\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}\left[\mathbb{1}_{Y_t\neq g(X_{\abar}^{\xtpast},Y_{\ytpast})}\right] \leq \epsilon + \delta.
\end{equation*}

At the first iteration of the algorithm $T(X_A)=T(\emptyset)=0$. Iteratively applying Equation (<ref>), at the end of the algorithm we have 8.59$TE_{X_{A}\rightarrow Y| X_{\abar}}=\sum_{k=1}^K TE_{X_{k}\rightarrow Y|X_{\abar_k}\setminus X_{k}}$, where $K$ is the number of iterations and $X_k,X_{\abar_k}$ are the removed feature and the set of candidate selected features during the $k$-th iteration. Since the algorithm stops when this quantity reaches the threshold (line <ref>), at the end we have $TE_{X_{A}\rightarrow Y| X_{\abar}}\leq \frac{\delta}{2B^2}$($\frac{\delta^2}{2}$ in classification). Finally, starting from Equation (<ref>) (Equation (<ref>) in classification), substituting this upper bound on $TE_{X_{A}\rightarrow Y| X_{\abar}}$ the result follows.

Forward Approach

The pseudo-code of Forward TEFS can be found in Algorithm <ref>,  defined in terms of the feature and target lags $M$ and $L$ and in terms of the minimum information gain $\Delta$. As in the backward case, $X_{A}$ is initialised to the empty set, and all features are in the complementary set $X_{\abar}$ (lines <ref>-<ref>). In this forward case, however, the selected features are represented by the set $X_A$. At each iteration, the algorithm selects the feature, among those not selected, that maximises the TE with the target, conditioned to the already selected features (line <ref>).
Finally, the stopping condition ensures that the procedure ends when the amount of information flowing from the features to the target is sufficiently large, depending on the hyperparameter $\Delta$ (line <ref>). A small value of the minimum information gain $\Delta$ leads to the selection of fewer features and a greater loss, a large value induces the selection of many features and a smaller loss of information. Furthermore, the hyperparameters $L,M$ regulate the bias-variance tradeoff as discussed in Remark <ref>.

The following lemma proves that the overall information obtained with the Forward TEFS, which selects the set $X_A$, is the cumulative conditional TE maximised at each step.

Let $X_i\in X_{\abar}$ be a feature in the complementary set of $X_A$. Then, the amount of information that flows from $X_A$ and the feature $X_i$ to the target $Y$ in terms of TE is:

\begin{equation}\label{eq:TEaddingAFeature}
TE_{X_{A}\cup X_i \rightarrow Y} = TE_{X_{A}\rightarrow Y} + TE_{X_{i}\rightarrow Y|X_{A}}.
\end{equation}

Intuitively, Lemma <ref> shows that the information flowing from $X_A$ and the feature $X_i$ to the target $Y$ is the information flowing from $X_A$ to the target plus the additional information flowing from $X_i$ to the target, given the information already flowing from the features in $X_A$.

Finally, in the following theorem we derive an upper bound on the ideal error, which identifies in the hyperparameter $\Delta$ the way to regulate the minimum information gain desired in the selection procedure of Algorithm <ref>.

The regression error suffered considering the variables $X_{A}$ selected by Forward TEFS is bounded by:
\begin{equation*}
\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}[(Y_t-g(X_{A}^{\xtpast},Y_{\ytpast}))^2] \leq \sigma^2+2B^2 \cdot TE_{X\rightarrow Y}-\Delta.
\end{equation*}

Similarly, the classification error is bounded by:
\begin{equation*}
\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}\left[\mathbb{1}_{Y_t\neq g(X_{A}^{\xtpast},Y_{\ytpast})}\right] \leq \epsilon + \sqrt{2\cdot TE_{X\rightarrow Y}-\Delta^2}.
\end{equation*}

Similarly to Theorem <ref>, the proof exploits the recursive relationship of Lemma <ref> to bound the final amount of information with the hyperparameter $\Delta$. Then, reformulating the upper bound on the ideal error of Theorem <ref> and <ref>, the results follow. Details can be found in Appendix <ref>.

§ TRANSFER ENTROPY ESTIMATION AND FINITE-SAMPLE ANALYSIS

In this section, we extend the theoretical analysis by discussing the effect of estimating the TE with a finite number of samples. We provide an upper bound on the regression and classification error that holds in high probability with a specific kernel estimator of the TE.

TE Estimation
As discussed, TE is a specific instance of a CMI, whose finite-sample estimation is a challenging task. Existing approaches rely on kernel density or $k$-nearest neighbour estimators. A traditional $k$-nearest neighbour estimator for continuous variables has been proposed in [Kraskov, Stögbauer, and Grassberger, 2004], and it has been extended in more recent works [Gao et al., 2017, Runge, 2018]. The usual theoretical guarantee of these approaches is asymptotic consistency, while a finite-sample analysis is little studied and remains largely an open problem. Considering continuous features and a continuous target,  [Singh and Póczos, 2014, Singh and Póczos, 2014] show a finite-sample exponential concentration inequality for density functional estimators and extend this theoretical result to conditional density functionals, respectively. In our analysis, we consider the mirrored kernel density estimator proposed in Section 4 of [Singh and Póczos, 2014]. It first estimates probability density functions through a modified kernel density estimator.
Then, the final estimator is obtained by plugging the estimated density into the functional of interest (i.e., CMI for our purposes).

The resulting CMI estimate 8.59$\hat{I}(X;Y|Z)$, made with $n$ samples, enjoys the following concentration bound with probability at least $1-\eta$ (Section 6.2 of  Singh and Póczos, 2014):

\begin{equation}\label{eq:CMIconcentration}
|I(X;Y|Z)-\hat{I}(X;Y|Z)| \leq \underbrace{\sqrt{\frac{4C_V^2}{n}\log\left(\frac{2}{\eta}\right)}}_{\text{(i)}}+\underbrace{C_{B}n^{\frac{-\beta}{\beta+d}}}_{\text{(ii)}},
\end{equation}

where $d=d_x+d_y+d_z$ is the sum of the dimensions of the three vectors $X$, $Y$, and $Z$. The terms $C_V,C_B>0$ are constant the number of sample $n$ but exponentially dependent on the dimension $d$ of the variables.[Intuitively, high dimensional vectors lead to larger uncertainty on the estimator, which may suffer from curse of dimensionality.]
$\beta\in [0,1]$ is the order of the Hölder space to which the joint density $p(x,y,z)$ belongs (e.g., $\beta=1$ implies Lipschitz continuity).[We refer the interested reader to Section 2 of [Singh and Póczos, 2014] for a detailed presentation.]
Intuitively, the term (i) in  Equation (<ref>) accounts for the variance of the density estimate and converges to zero with the optimal rate $\mathcal{O}(n^{-\frac{1}{2}})$. Instead, the term (ii) refers to the bias of the density estimator and displays a convergence rate of order  $\mathcal{O}(n^{\frac{-\beta}{\beta+d}})$, which is therefore the overall convergence rate.

Finite-Sample Bound
Recalling Theorem <ref> and <ref> in the ideal case (i.e., when the TE is assumed to be exact), the following result rephrases them in the finite-sample case (i.e., when the TE is estimated from the kernel density estimator discussed above). The proof can be found in Appendix <ref>.

Three variables.

Five variables.

Ten variables.

Causal graphs considered for synthetic experiments. Each edge reports the time lag and the interaction coefficient.

The regression error suffered considering the variables $X_{\abar}$ selected by Backward TEFS with the estimator of CMI proposed in  [Singh and Póczos, 2014], with $n$ samples, is bounded, with probability at least $1-\eta$, by:

\begin{align*}
&\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}[(Y_t-g(X_{\abar}^{\xtpast},Y_{\ytpast}))^2] \\ &\leq \sigma^2+\delta+2B^2K\left( \sqrt{\frac{4C_V^2}{n}\log\left(\frac{2^{d+1}}{\eta}\right)}+C_{B}n^{\frac{-\beta}{\beta+d}}\right) ,
\end{align*}

where $K$ is the number of iterations performed by the algorithm and $d$ is the dimension of the full set of variables ($d=LD+M+1$).

The regression error suffered considering the selected variables $X_{A}$ by running Forward TEFS, under the same conditions,  is bounded, with probability at least $1-\eta$, by:

\begin{align*}
\begin{split}
&\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}[(Y_t-g(X_{A}^{\xtpast},Y_{\ytpast}))^2] \\&\leq \sigma^2+2B^2 \cdot TE_{X\rightarrow Y}-\Delta\\& +2B^2K\left( \sqrt{\frac{4C_V^2}{n}\log\left(\frac{2^{d+1}}{\eta}\right)}+C_{B}n^{\frac{-\beta}{\beta+d}} \right).
\end{split}
\end{align*}

Analogous results hold in classification, extending the ideal Bayes error of Theorem <ref>, <ref> in the finite-sample case.

§ EXPERIMENTAL VALIDATION

In this section, we conduct an experimental validation of the proposed algorithms on both synthetic (Section <ref>) and real-world (Section <ref>) domains.[Code and datasets can be found at: <https://www.dropbox.com/scl/fo/ctvxqh8q131f8sim4s22o/h?dl=0   rlkey=5hsbpold82bi8vilu44ubb508>.]

§.§ Synthetic Experiments

We analyse the behaviour of the forward and backward TEFS the forward and backward feature selection algorithms based on CMI proposed in [Beraha et al., 2019] and with the PCMCI algorithm [Runge et al., 2019], a state-of-the-art causal discovery approach. Details and additional results can be found in Appendix <ref>.

Figure <ref> shows the graph built to perform the experiments in three, five, and ten dimensions, with the time lag of the interaction between each pair of variables and the coefficient of their linear interaction. Indeed, the structural causal model of each variable is assumed to be linear, with additive Gaussian noise. Each of the three experiments is repeated ten times with different seeds, considering the coefficients reported in the figure, $N=300$ samples, maximum time lags $M=L=2$ for the two and five-dimensional data, and $M=L=3$ for the ten-dimensional data, and additive Gaussian noise with standard deviation $s=0.1$. Then, the individual effect of each of these parameters is studied, fixing the other parameters. To further inspect the impact of high dimensionality, we also consider $D\in\{15,20,40,60,80,100\}$, preserving the relationships and coefficients as in Figure <ref> and adding independent variables that can be autocorrelated or interact with each other.

Each model is evaluated with the true positive rate (TPR) and false positive rate (FPR) on the number of causal links detected, leading to some considerations on the algorithms. Both forward and backward CMI feature selection do not consider the target in the conditioning. Hence they tend to select features that are caused by the target as well. On the contrary, forward and backward TEFS and the PCMCI algorithm tend to identify the features causally relevant to the target, along with its autoregressive component. Additionally, PCMCI tends to consider more features to avoid false negatives at the cost of a higher false positive rate. In contrast, the two versions of the TEFS algorithm are more focused on identifying the causal features that are most significant in terms of regression performance. This translates into better control of the FPR, at the cost of possibly not considering some less relevant causal links, worsening the TPR. Table <ref> shows, for the three-dimensional setting, the complete results in terms of TPR and FPR,

In Appendix <ref>, it is possible to find the same results for a larger number of features $D\in\{5,10,15,20,40,60,80,100\}$.

Three-dimensional experiments on synthetic datasets. Each experiment has been repeated ten times, with different seeds, varying one parameter at a time.

Forward CMI     Backward CMI     PCMCI     lightgray!25 Forward TEFS     lightgray!25 Backward TEFS
Benchmark                 lightgray!25     lightgray!25

TPR     $0.5$     $0.5$     $1.0$     lightgray!25 $1.0$     lightgray!25 $1.0$

FPR     $1$     $1$     $0.1$     lightgray!25 $0$     lightgray!25 $0$

Increasing Noise                 lightgray!25     lightgray!25

TPR     $[0.5, 0.5, 0.5, 0.5, 0.45]$     $[0.5, 0.5, 0.5, 0.5, 0.45]$      $[1, 1, 1, 1, 1]$     lightgray!25 $[1, 1, 1, 0.95, 0.9]$     lightgray!25 $[1, 1, 1, 1, 1]$

FPR     $[1,1,1,1,1]$     $[1,1,1,1,1]$     $[0.1, 0, 0.1, 0.3, 0.4]$     lightgray!25 $[0, 0, 0, 0, 0]$     lightgray!25 $[0, 0, 0, 0, 0]$

Increasing $\tau$                 lightgray!25     lightgray!25

TPR     $[0.5, 0.5, 0.5, 0.5]$     $[0.5, 0.5, 0.5, 0.5]$      $[1, 1, 1, 1]$     lightgray!25 $[1, 1, 1, 1]$     lightgray!25 $[1, 1, 1, 1]$

FPR     $[1,1,1,1]$     $[1,1,1,1]$     $[0.1, 0, 0.3, 0.3]$     lightgray!25 $[0, 0, 0, 0]$     lightgray!25 $[0, 0, 0, 0]$

Increasing N                 lightgray!25     lightgray!25

TPR     $[0.5, 0.5, 0.5, 0.5, 0.5]$     $[0.5, 0.5, 0.5, 0.5, 0.5]$     $[0.9, 1, 1, 1, 1]$     lightgray!25 $[1, 1, 1, 1, 1]$     lightgray!25 $[1, 1, 1, 1, 1]$

FPR     $[1,1,1,1,1]$     $[1,1,1,1,1]$     $[0.5, 0.2, 0.1, 0.1, 0.1]$     lightgray!25 $[0, 0, 0, 0, 0]$     lightgray!25 $[0, 0, 0, 0, 0]$

Random Coefficients                 lightgray!25     lightgray!25

TPR     $[0.5,0.5,0.5,0.5,0.3]$     $[0.5,0.5,0.5,0.5,0.3]$     $[1,1,1,0.75,1]$     lightgray!25 $[1,1,1,1,1]$     lightgray!25 $[1,1,1,1,1]$

FPR     $[1, 1, 1, 1, 1]$     $[1, 1, 1, 1, 1]$     $[0.2,0.2,0.1,0.1,0.2]$     lightgray!25 $[0,0,0,0,0.1]$     lightgray!25 $[0,0,0,0,0.1]$

8.58.5Real-world experiments for three different datasets. Each experiment has been repeated three times, considering different time-depths $M=L\in\{1,2,3\}$. The table reports the set of IDs of the selected features and the related $R^2$ test score.

Best_CMI     Best_CI     lightgray!25 Forward TEFS      lightgray!25 Backward TEFS
Climate, 5 features             lightgray!25     lightgray!25

$M=L=1$     $\{0,3,4\}\rightarrow0.34$     $\{0,2,4\}\rightarrow0.45$     lightgray!25 $\{0,4\}\rightarrow0.44$     lightgray!25 $\{0,4\}\rightarrow0.44$

$M=L=2$     $\{0,1\}\rightarrow0.36$     $\{0,1,3,4\}\rightarrow0.47$     lightgray!25 $\{0,4\}\rightarrow0.42$     lightgray!25 $\{0,4\}\rightarrow0.42$

$M=L=3$     $\{0,1\}\rightarrow0.35$     $\{0,3,4\}\rightarrow0.48$     lightgray!25 $\{0,4\}\rightarrow0.47$     lightgray!25 $\{0,4\}\rightarrow0.47$

Climate, 15 features             lightgray!25     lightgray!25

$M=L=1$     $\{0,2,3,11\}\rightarrow0.29$     $\{0,2,7\}\rightarrow0.45$     lightgray!25 $\{0,7\}\rightarrow0.44$     lightgray!25 $\{0,11\}\rightarrow0.44$

$M=L=2$     $\{0,2,11\}\rightarrow0.30$     $\{5,7,11\}\rightarrow0.47$     lightgray!25 $\{0,7\}\rightarrow0.43$     lightgray!25 $\{0,11\}\rightarrow0.42$

$M=L=3$     $\{0,5\}\rightarrow0.28$     $\{2,7\}\rightarrow0.48$     lightgray!25 $\{7\}\rightarrow0.47$     lightgray!25 $\{11\}\rightarrow0.43$

Benchmark, 5 features             lightgray!25     lightgray!25

$M=L=1$     $\{3,4\}\rightarrow0.09$     $\{3,4\}\rightarrow0.09$     lightgray!25 $\{3,4\}\rightarrow0.09$     lightgray!25 $\{3,4\}\rightarrow0.09$

$M=L=2$     $\{3\}\rightarrow0.07$     $\{3,4\}\rightarrow0.13$     lightgray!25 $\{4\}\rightarrow0.09$     lightgray!25 $\{4\}\rightarrow0.09$

$M=L=3$     $\{0,3\}\rightarrow0.07$     $\{0,3,4\}\rightarrow0.13$     lightgray!25 $\{3\}\rightarrow0.11$     lightgray!25 $\{3\}\rightarrow0.11$

§.§ Real-World Experiments

We conclude the validation with three applications to real-world data. In particular, we firstly analyse a climatic datasets where the target represents the vegetation state of a basin of the Po River, and the five features are meteorological variables related to temperature and precipitation. To increase the dimensionality of the data, we also consider a similar dataset with fifteen meteorological features. Finally, to further test the algorithm on a benchmark dataset, we consider a dataset composed of five variables available in an online causal repository (<https://causeme.uv.es/>), where we consider in turn each variable to be the target and the other four variables the candidate causal features.
Ideally, we would know the real underlying causal relationships to validate the algorithms, which is unrealistic in complex problems. For this reason, motivated by the feature selection inspiration of this work, we propose the regression performance (in terms of coefficient of determination $R^2$) to draw some conclusions on the effectiveness of the methods. To further inspect the performance of the proposed approaches, along with forward and backward feature selection based on CMI, we consider as causal discovery baselines PCMCI and FullCI, which are implemented in the same library (<https://github.com/jakobrunge/tigramite.git>). These two methods are applied considering both linear and non-linear independence tests. Table <ref> shows the selected set of features and the associated $R^2$ test score, considering increasing past depth $L=M\in\{1,2,3\}$. For clarity, only the best-performing result among the linear and non-linear PCMCI and FullCI applications (best_CI) is reported. The same applies to the forward and backward CMI approaches (best_CMI). Furthermore, the three datasets in the table are one climatic target with five features, the climatic target with fifteen features, and the benchmark dataset with five variables, where the third variable is taken as the target. The complete results for all the considered datasets can be found in Appendix <ref>. As discussed for the synthetic experiments, from the results, it is possible to observe the tendency of the Forward and Backward TEFS, PCMCI and FullCI, to consider a reduced set of the most related causal features in terms of TE. In the empirical results, this leads to an empirical performance that preserves most of the information while significantly reducing dimensionality. Therefore, the $R^2$ scores of the proposed algorithms are similar to or slightly worse than the best causal discovery counterparts, which consider more features as they aim to identify all possible causal links.
More extensive applications of Forward and Backward TEFS on complex real data are the subject of future research.

§ CONCLUSIONS AND FUTURE DEVELOPMENTS

In this paper, we introduced two novel causal feature selection approaches for time series, Forward and Backward TEFS, relying on conditional TE to guide the selection. We provided an analysis of the information loss due to the selection procedure in the exact and in the finite-sample cases. Our theoretical results validate the soundness of the two approaches, formally highlighting the role of the hyperparameters in controlling the final selection.

Then, we conducted synthetic experiments in comparison with standard CMI feature selection and the state-of-the-art PCMCI, highlighting the abilities of our method to successfully uncover causal relationships.

Finally, we performed experiments on real-world datasets (where no ground truth of causal relationships is present), proving that our methods display competitive results in terms of regression scores.
A promising future development is related to the invariance framework, which relies on the intuition that a variable causing the target maintains the same relationship with the target itself, even if its
own distribution changes [Arjovsky et al., 2019, Bühlmann, 2020]. The concentration bound of Equation (<ref>) could be therefore exploited to test the capability of the selected features to generalize across heterogeneous environments.

[Amblard and Michel, 2012]
Amblard, P.-O.; and Michel, O. J. 2012.
The relation between Granger causality and directed information
theory: A review.
Entropy, 15(1): 113–143.

[Arjovsky et al., 2019]
Arjovsky, M.; Bottou, L.; Gulrajani, I.; and Lopez-Paz, D. 2019.
Invariant risk minimization.
arXiv preprint arXiv:1907.02893.

[Beraha et al., 2019]
Beraha, M.; Metelli, A. M.; Papini, M.; Tirinzoni, A.; and Restelli, M. 2019.
Feature selection via mutual information: New theoretical insights.
In 2019 international joint conference on neural networks
(IJCNN), 1–9. IEEE.

[Bishop, 2006]
Bishop, C. M. 2006.
Pattern Recognition and Machine Learning (Information Science
and Statistics).
Berlin, Heidelberg: Springer-Verlag.
ISBN 0387310738.

[Bühlmann, 2020]
Bühlmann, P. 2020.
Invariance, causality and robustness.

[Cover and Thomas, 2006]
Cover, T. M.; and Thomas, J. A. 2006.
Elements of Information Theory (Wiley Series in
Telecommunications and Signal Processing).
USA: Wiley-Interscience.
ISBN 0471241954.

[D'Amour et al., 2020]
D'Amour, A.; Ding, P.; Feller, A.; Lei, L.; and Sekhon, J. 2020.
Overlap in Observational Studies with High-Dimensional Covariates.

[Ding, Chen, and Bressler, 2006]
Ding, M.; Chen, Y.; and Bressler, S. L. 2006.
Granger causality: basic theory and application to neuroscience.
Handbook of time series analysis: recent theoretical
developments and applications, 437–460.

[Dong and Kluger, 2022]
Dong, M.; and Kluger, Y. 2022.
GEASS: Neural causal feature selection for high-dimensional
biological data.
In The Eleventh International Conference on Learning

[Eberhardt, 2017]
Eberhardt, F. 2017.
Introduction to the foundations of causal discovery.
International Journal of Data Science and Analytics, 3:

[Gao et al., 2017]
Gao, W.; Kannan, S.; Oh, S.; and Viswanath, P. 2017.
Estimating mutual information for discrete-continuous mixtures.
Advances in neural information processing systems, 30.

[Gerhardus and Runge, 2020]
Gerhardus, A.; and Runge, J. 2020.
High-recall causal discovery for autocorrelated time series with
latent confounders.
Advances in Neural Information Processing Systems, 33:

[Glymour, Zhang, and Spirtes, 2019]
Glymour, C.; Zhang, K.; and Spirtes, P. 2019.
Review of causal discovery methods based on graphical models.
Frontiers in genetics, 10: 524.

[Granger, 1969]
Granger, C. W. 1969.
Investigating causal relations by econometric models and
cross-spectral methods.
Econometrica: journal of the Econometric Society, 424–438.

[Imbens and Rubin, 2015]
Imbens, G. W.; and Rubin, D. B. 2015.
Causal inference in statistics, social, and biomedical
Cambridge University Press.

[Kraskov, Stögbauer, and Grassberger, 2004]
Kraskov, A.; Stögbauer, H.; and Grassberger, P. 2004.
Estimating mutual information.
Physical review E, 69(6): 066138.

[Li et al., 2017]
Li, J.; Cheng, K.; Wang, S.; Morstatter, F.; Trevino, R. P.; Tang, J.; and Liu,
H. 2017.
Feature selection: A data perspective.
ACM computing surveys (CSUR), 50(6): 1–45.

[Liu and Aviyente, 2012]
Liu, Y.; and Aviyente, S. 2012.
The relationship between transfer entropy and directed information.
In 2012 IEEE Statistical Signal Processing Workshop (SSP),
73–76. IEEE.

[Malinsky and Spirtes, 2018]
Malinsky, D.; and Spirtes, P. 2018.
Causal structure learning from multivariate time series in settings
with unmeasured confounding.
In Proceedings of 2018 ACM SIGKDD workshop on causal
discovery, 23–47. PMLR.

[Massey et al., 1990]
Massey, J.; et al. 1990.
Causality, feedback and directed information.
In Proc. Int. Symp. Inf. Theory Applic.(ISITA-90), 303–305.

[Mastakouri and Schölkopf, 2020]
Mastakouri, A.; and Schölkopf, B. 2020.
Causal analysis of Covid-19 spread in Germany.
Advances in Neural Information Processing Systems, 33:

[Mastakouri, Schölkopf, and Janzing, 2021]
Mastakouri, A. A.; Schölkopf, B.; and Janzing, D. 2021.
Necessary and sufficient conditions for causal feature selection in
time series with latent common causes.
In International Conference on Machine Learning, 7502–7511.

[Moneta et al., 2011]
Moneta, A.; Chlaß, N.; Entner, D.; and Hoyer, P. 2011.
Causal search in structural vector autoregressive models.
In NIPS Mini-Symposium on Causality in Time Series, 95–114.

[Moraffah et al., 2021]
Moraffah, R.; Sheth, P.; Karami, M.; Bhattacharya, A.; Wang, Q.; Tahir, A.;
Raglin, A.; and Liu, H. 2021.
Causal inference for time series analysis: Problems, methods and
Knowledge and Information Systems, 63: 3041–3085.

[Pearl, 2009]
Pearl, J. 2009.
Cambridge university press.

[Pearl et al., 2000]
Pearl, J.; et al. 2000.
Models, reasoning and inference.
Cambridge, UK: CambridgeUniversityPress, 19(2).

[Peters, Janzing, and Schölkopf, 2017]
Peters, J.; Janzing, D.; and Schölkopf, B. 2017.
Elements of causal inference: foundations and learning
The MIT Press.

[Runge, 2018]
Runge, J. 2018.
Conditional independence testing based on a nearest-neighbor
estimator of conditional mutual information.
In International Conference on Artificial Intelligence and
Statistics, 938–947. PMLR.

[Runge, 2020]
Runge, J. 2020.
Discovering contemporaneous and lagged causal relations in
autocorrelated nonlinear time series datasets.
In Conference on Uncertainty in Artificial Intelligence,
1388–1397. PMLR.

[Runge et al., 2019]
Runge, J.; Bathiany, S.; Bollt, E.; Camps-Valls, G.; Coumou, D.; Deyle, E.;
Glymour, C.; Kretschmer, M.; Mahecha, M. D.; Muñoz-Marí, J.; et al.
Inferring causation from time series in Earth system sciences.
Nature communications, 10(1): 2553.

[Runge et al., 2019]
Runge, J.; Nowack, P.; Kretschmer, M.; Flaxman, S.; and Sejdinovic, D.
Detecting and quantifying causal associations in large nonlinear time
series datasets.
Science advances, 5(11): eaau4996.

[Schölkopf, 2022]
Schölkopf, B. 2022.
Causality for machine learning.
In Probabilistic and Causal Inference: The Works of Judea
Pearl, 765–804.

[Schreiber, 2000]
Schreiber, T. 2000.
Measuring information transfer.
Physical review letters, 85(2): 461.

[Shahsavari Baboukani et al., 2020]
Shahsavari Baboukani, P.; Graversen, C.; Alickovic, E.; and Østergaard, J.
Estimating Conditional Transfer Entropy in Time Series Using Mutual
Information and Nonlinear Prediction.
Entropy, 22(10).

[Singh and Póczos, 2014]
Singh, S.; and Póczos, B. 2014a.
Exponential concentration of a density functional estimator.
Advances in Neural Information Processing Systems, 27.

[Singh and Póczos, 2014]
Singh, S.; and Póczos, B. 2014b.
Generalized exponential concentration inequality for Rényi
divergence estimation.
In International Conference on Machine Learning, 333–341.

[Spirtes et al., 2000]
Spirtes, P.; Glymour, C. N.; Scheines, R.; and Heckerman, D. 2000.
Causation, prediction, and search.
MIT press.

§ TRANSFER ENTROPY FOR CAUSAL FEATURE SELECTION: ADDITIONAL PROOFS AND RESULTS

This section contains proofs and additional results of Section <ref> and <ref>. To simplify the notation, $X_A^{-t},X_{\abar}^{-t}$ will denote the values of the set of features with indices in $A$ or its complementary $\abar$ in the timesteps $\xtpast$. Moreover, $Y_{-t}$ will denote the last values of the target $Y_{\ytpast}$ and $Y_t$ will be the current value of the target itself.

Since the best regression model, given the inputs, is the conditional expected value of the target, the following equality holds:

\begin{equation*}
\begin{aligned}
\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}[(Y_t-g(X_{\abar}^{-t},Y_{-t}))^2] \\= \E_{X,Y}[(Y_t-\E[Y_t|X_{\abar}^{-t},Y_{-t}])^2].
\end{aligned}
\end{equation*}

Adding and subtracting $\E[Y_t|X^{-t},Y_{-t}]$ to the inner term and writing the outer expected value as integral, the term above is equal to:

\begin{equation*}
\begin{aligned}
&\int p(x^{-t},y_{-t}) \int p(y_t|x^{-t},y_{-t}) \big( y_t-\E[Y_t|X_{\abar}^{-t},Y_{-t}] \\ &\pm \E[Y_t|X^{-t},Y_{-t}] \big)^2 dy_tdx^{-t}dy_{-t}\\
=& \int p(x^{-t},y_{-t}) \int p(y_t|x^{-t},y_{-t}) \big( y_t-\E[Y_t|X^{-t},Y_{-t}] \big)^2 dy_tdx^{-t}dy_{-t}\\
+&\int p(x^{-t},y_{-t}) \int p(y_t|x^{-t},y_{-t}) \left(\E[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X_{\abar}^{-t},Y_{-t}] \right)^2\\ %dy_tdx^{-t}dy_{-t} \\
+&2\E_{X,Y}[( y_t-\E[Y_t|X^{-t},Y_{-t}] )\cdot (\E[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X_{\abar}^{-t},Y_{-t}] )].
\end{aligned}
\end{equation*}

We now analyse each term of the right-hand side of this equation.

The first term is the irreducible mean squared error:

\begin{equation*}
\begin{aligned}
\int p(x^{-t},y_{-t}) & \int p(y_t|x^{-t},y_{-t}) \left( y_t-\E[Y_t|X^{-t},Y_{-t}] \right)^2 dy_tdx^{-t}dy_{-t} \\
&=\E_{X,Y}[(Y_t-\E[Y_t|X^{-t},Y_{-t})^2] = \sigma^2.
\end{aligned}
\end{equation*}

The third term is equal to zero:

\begin{equation*}
\begin{aligned}
\E_{X,Y}[ &( Y_t-\E[Y_t|X^{-t},Y_{-t}] )\cdot (\E[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X_{\abar}^{-t},Y_{-t}] )]\\
=&\E_{X,Y_{-t}}\big[\E_{Y_t}[( Y_t-\E[Y_t|X^{-t},Y_{-t}] ) \\ &\cdot (\E[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X_{\abar}^{-t},Y_{-t}] )]|X^{-t},Y_{-t}\big]\\
=&\E_{X,Y_{-t}}[( \E_{Y_t}[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X^{-t},Y_{-t}] )\\&\cdot (\E[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X_{\abar}^{-t},Y_{-t}] )]=0,
\end{aligned}
\end{equation*}

where the second equality holds since only the first term ($Y_t$) depends on $Y_t$ itself.

Finally, the second term can be bounded as:

\begin{align*}
& \int p(x^{-t},y_{-t}) \int p(y_t|x^{-t},y_{-t}) \\&\left(\E[Y_t|X^{-t},Y_{-t}]-\E[Y_t|X_{\abar}^{-t},Y_{-t}] \right)^2 dy_tdx^{-t}dy_{-t}\\
= &  \int p(x^{-t},y_{-t}) \int p(y_t|x^{-t},y_{-t}) \\&\left(\int y_t\cdot (p(y_t|x^{-t},y_{-t})-p(y_t|x_{\abar}^{-t},y_{-t}))dy_t\right)^2 dy_tdx^{-t}dy_{-t}\\
\leq & B^2 \int p(x^{-t},y_{-t}) \\&\left(\int |p(y_t|x^{-t},y_{-t})-p(y_t|x_{\abar}^{-t},y_{-t})|dy_t\right)^2 dx^{-t}dy_{-t}\\
= & B^2 \int p(x^{-t},y_{-t}) \\& \left(2\cdot D_{TV}(p(y_t|x^{-t},y_{-t})||p(y_t|x_{\abar}^{-t},y_{-t}))\right)^2 dx^{-t}dy_{-t}\\
\leq & B^2 \int p(x^{-t},y_{-t}) \\&\left(2\cdot \sqrt{\frac{1}{2}D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_t|x_{\abar}^{-t},y_{-t}))}\right)^2 dx^{-t}dy_{-t}\\
= & 2B^2 \int p(x^{-t},y_{-t}) D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_t|x_{\abar}^{-t},y_{-t})) dx^{-t}dy_{-t}\\
= & 2B^2 \E_{X,Y_{-t}}[D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_t|x_{\abar}^{-t},y_{-t}))]\\
= & 2B^2I(Y_t,X_A^{-t}|X_{\abar}^{-t},Y_{-t}) = 2B^2\cdot TE_{X_{A}\rightarrow Y| X_{\abar}}.
\end{align*}

The first equality rewrites the inner expected values as integrals. Then, the inequality holds applying the absolute value in the inner term and recalling that $|Y_t|\leq B$ by assumption. The second equality follows since the total variation distance is half the $L_1$ norm. Finally, the second inequality follows from Pinsker's inequality and the last three equations simply rearrange the terms and apply the definitions of conditional mutual information and transfer entropy.

Combining the three terms analysed, the result of the theorem follows.

The optimal class that can be predicted, given the inputs, is the one with the largest probability. Therefore, it is equal to 7.58$y_t^* = \arg\max_{y_t\in\mathcal{Y}}p(y_t|X^{-t},Y_{-t})$ or 7.58$y^*_{t,\abar} = \arg\max_{y_t\in\mathcal{Y}}p(y_t|X_{\abar}^{-t},Y_{-t})$ respectively if the full set of features or only the selected ones is considered.

The expected prediction error considering the selected features, for the best model, is therefore equal to:

\begin{equation*}
\begin{aligned}
&\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}\left[\mathbb{1}_{\{Y_t\neq g(X_{\abar}^{-t},Y_{-t})\}}\right] = \E_{X,Y}\left[ \mathbb{1}_{\{Y_t\neq Y^*_{t,\abar}\}} \right] \\&= \epsilon + \E_{X,Y}\left[ \mathbb{1}_{\{Y_t\neq Y^*_{t,\abar}\}}-\mathbb{1}_{\{Y_t\neq Y^*_{t}\}} \right],
\end{aligned}
\end{equation*}

where the second equality follows summing and subtracting the term $\epsilon=\E_{X,Y}\left[ \mathbb{1}_{\{Y_t\neq Y^*_{t}\}}\right]$, that is equal to the irreducible Bayes error suffered considering the entire set of features.

We now focus on the expected value obtained, rewriting it as an integral:

\begin{align*}
& \E_{X,Y} \left[ \mathbb{1}_{\{Y_t\neq Y^*_{t,\abar}\}}-\mathbb{1}_{\{Y_t\neq Y^*_{t}\}} \right] \\
& = \int p(x^{-t},y_{-t}) ( \E_{y_t}[ \mathbb{1}_{\{Y_t\neq Y^*_{t,\abar}\}}  |X^{-t},Y_{-t}] \\& - \E_{y_t}[ \mathbb{1}_{\{Y_t\neq Y^*_{t}\}} |X^{-t},Y_{-t}]) dx^{-t}dy_{-t}.
\end{align*}

Focusing on the inner expected values:

\begin{align*}
&\E_{y_t}\left[ \mathbb{1}_{\{Y_t\neq Y^*_{t}\}}|X^{-t},Y_{-t} \right] = 1-\E_{y_t}\left[ \mathbb{1}_{\{Y_t= Y^*_{t}\}}|X^{-t},Y_{-t} \right] \\ & = 1-p(Y^*_{t}|X^{-t},Y_{-t}),
\end{align*}

where the equalities follow from the definition of the expected value of the complement of a set and of the expected value of the indicator function. Intuitively, the equality states that the expected error when the target is estimated with its maximum probability given the features and the past of the target is equal to $1$ minus the probability of correctly estimating the target with its maximum conditional expected probability given the features and the past values of the target. A similar result holds conditioning only on the selected features, allowing rewriting the equation above as:

\begin{align*}
\E_{X,Y} & \left[ \mathbb{1}_{\{Y_t\neq Y^*_{t,\abar}\}}-\mathbb{1}_{\{Y_t\neq Y^*_{t}\}} \right] \\
& = \int p(x^{-t},y_{-t}) (p(y^*_{t}|x^{-t},y_{-t}) \\
& -p(y^*_{t,\abar}|x^{-t},y_{-t})) dx^{-t}dy_{-t}\\
& = \int p(x^{-t},y_{-t}) (p(y^*_{t}|x^{-t},y_{-t}) -p(y^*_{t,\abar}|x^{-t},y_{-t}) \\& \pm p(y^*_{t,\abar}|x_{\abar}^{-t},y_{-t}) ) dx^{-t}dy_{-t}.\\
\end{align*}

Focusing on each inner term individually:

\begin{align*}
&p(y^*_{t}|x^{-t},y_{-t}) -p(y^*_{t,\abar}|x_{\abar}^{-t},y_{-t}) \\& = \max_{y_t\in \mathcal{Y}} p(y_t|x^{-t},y_{-t})-\max_{y_t\in \mathcal{Y}} p(y_{t}|x_{\abar}^{-t},y_{-t})\\
& \leq \max_{y_t\in \mathcal{Y}} | p(y_t|x^{-t},y_{-t})-p(y_{t}|x_{\abar}^{-t},y_{-t}) | \\
&\leq D_{TV}(p(y_t|x^{-t},y_{-t})||p(y_{t}|x_{\abar}^{-t},y_{-t}))\\
& \leq \sqrt{\frac{1}{2}D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_{t}|x_{\abar}^{-t},y_{-t}))},
\end{align*}

where the first equality follows by definition of $y_t^*,y^*_{t,\abar}$, the first inequality holds since the difference of maxima is smaller than or equal to the maximum absolute value of the difference and the second inequality holds by definition of total variation distance. Finally, applying Pinsker's inequality, the third inequality holds. Similarly, the second term inside the integral can be bounded as follows:

\begin{align*}
&p(y^*_{t,\abar}|x_{\abar}^{-t},y_{-t})-p(y^*_{t,\abar}|x^{-t},y_{-t}) \\
&\leq \sqrt{\frac{1}{2}D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_{t}|x_{\abar}^{-t},y_{-t}))}.
\end{align*}

Therefore, the integral becomes:

\begin{align*}
&\int p(x^{-t},y_{-t})  (p(y^*_{t}|x^{-t},y_{-t}) -p(y^*_{t,\abar}|x^{-t},y_{-t}) \\&\pm p(y^*_{t,\abar}|x_{\abar}^{-t},y_{-t}) ) dx^{-t}dy_{-t}\\
& \leq \int p(x^{-t},y_{-t}) \sqrt{2D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_{t}|x_{\abar}^{-t},y_{-t}))} dx^{-t}dy_{-t}\\
& \leq \sqrt{2\int p(x^{-t},y_{-t}) D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_{t}|x_{\abar}^{-t},y_{-t}))} dx^{-t}dy_{-t}\\
& = \sqrt{2\E_{X,Y_{-t}} [D_{KL}(p(y_t|x^{-t},y_{-t})||p(y_{t}|x_{\abar}^{-t},y_{-t}))] } \\
& = \sqrt{2I(Y_t;X_A^{-t}|X_{\abar}^{-t},Y_{-t})} = \sqrt{2\cdot TE_{X_{A}\rightarrow Y| X_{\abar}}}.
\end{align*}

Substituting this result in the first equation of this proof we obtain Equation <ref> of the main paper, proving the theorem.

Let 7.58$X_{\hat{i}}=\arg\min_{X_i\in X_{\abar}} TE_{X_{i}\rightarrow Y|X_{\abar}\setminus X_{i}}$, then $T(X_A\cup X_{\hat{i}})$ is defined as:

\begin{align*}
&T(X_{A}\cup X_{\hat{i}}) = TE_{X_{A}\cup X_{\hat{i}}\rightarrow Y|X_{\abar}\setminus X_{\hat{i}}}\\&=I(Y_t;\{X_{A}\cup X_{\hat{i}}\}^{\xtpast}|\{X_{\abar}\setminus X_{\hat{i}}\}^{\xtpast},Y_{\ytpast}).
\end{align*}

Then, applying multiple times the chain rule of the conditional mutual information and the definition of 7.58$T(X_A)=TE_{X_{A}\rightarrow Y|X_{\abar}}$ and 7.58$T_{\abar}=\min_{X_i\in X_{\abar}} TE_{X_{i}\rightarrow Y|X_{\abar}\setminus X_{i}}$, the following equalities hold, proving the theorem:

\begin{align*}
&T(X_{A}\cup X_{\hat{i}}) = TE_{X_{A}\cup X_{\hat{i}}\rightarrow Y|X_{\abar}\setminus X_{\hat{i}}}\\
& = I(Y_t;X^{\xtpast},Y_{\ytpast}) - I(Y_t;\{{X_{\abar}\setminus X_{\hat{i}}}\}^{\xtpast},Y_{\ytpast})\\
& = I(Y_t;X_{\abar}^{\xtpast},Y_{\ytpast}) \\&+ I(Y_t;X_{A}^{\xtpast}|X_{\abar}^{\xtpast},Y_{\ytpast})\\
& \quad -I(Y_t;\{{X_{\abar}\setminus X_{\hat{i}}}\}^{\xtpast},Y_{\ytpast})\\
& = I(Y_t;\{X_{\abar}\setminus X_{\hat{i}}\}^{\xtpast},X_{\hat{i}}^{\xtpast},Y_{\ytpast}) + T(X_A)\\
& \quad - I(Y_t;\{{X_{\abar}\setminus X_{\hat{i}}}\}^{\xtpast},Y_{\ytpast})\\
& = T(X_A) + I(Y_t;X_{\hat{i}}^{\xtpast}|\{X_{\abar}\setminus X_{\hat{i}}\}^{\xtpast},Y_{\ytpast})\\
& \quad \pm I(Y_t;\{{X_{\abar}\setminus X_{\hat{i}}}\}^{\xtpast},Y_{\ytpast})\\
& = T(X_A)\cup T_{\abar}.
\end{align*}

Let us consider $TE_{X_{A}\cup X_i \rightarrow Y}$, which is defined in terms of conditional mutual information as:

$$ TE_{X_{A}\cup X_i \rightarrow Y} = I(Y_t;\{X_{A}\cup X_i\}^{\xtpast}|Y_{\ytpast}). $$

Therefore, applying the chain rule the result follows:

\begin{align*}
TE_{X_{A}\cup X_i \rightarrow Y} & = I(Y_t;\{X_{A}\cup X_i\}^{\xtpast}|Y_{\ytpast})\\
&  =   I(Y_t;X_i^{\xtpast}|X_{A}^{\xtpast},Y_{\ytpast}) \\
&+ I(Y_t;X_{A}^{\xtpast}|Y_{\ytpast})\\
& =
TE_{X_{i}\rightarrow Y|X_{A}} + TE_{X_{A}\rightarrow Y}.
\end{align*}

In the forward approach of Algorithm <ref>, which selects the set of features $X_A$, the ideal regression and classification error are respectively equal to 7.58$\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}[(Y_t-g(X_{A}^{\xtpast},Y_{\ytpast}))^2]$ and 7.58$\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}\left[\mathbb{1}_{Y_t\neq g(X_{A}^{\xtpast},Y_{\ytpast})}\right]$. To prove the theorem, we first need the following remark, which allows us to rewrite in a forward fashion the upper bound to the ideal error.

From Theorem <ref> and <ref> of the main paper, applying the chain rule of conditional mutual information, the ideal regression mean squared error and classification Bayes error can be bounded as:

\begin{align*}
&\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}[(Y_t-g(X_{A}^{\xtpast},Y_{\ytpast}))^2] \\
&\leq \sigma^2 + 2B^2\cdot TE_{X_{\abar}\rightarrow Y| X_{A}} = \sigma^2 + 2B^2\cdot (TE_{X\rightarrow Y} - TE_{X_A \rightarrow Y}),\\
&\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}\left[\mathbb{1}_{Y_t\neq g(X_{A}^{\xtpast},Y_{\ytpast})}\right] \leq \epsilon + \sqrt{2\cdot TE_{X_{\abar}\rightarrow Y| X_{A}}}\\
&= \epsilon + \sqrt{2\cdot (TE_{X\rightarrow Y} - TE_{X_A \rightarrow Y})}.
\end{align*}

This rephrasing of the statements shows that the maximisation performed in the feature selection approach of Algorithm <ref> is consistent with the theory, since the maximisation of the term $TE_{X_A \rightarrow Y}$ minimizes the error.

We are now ready to prove the relationship between the regression and classification error and the chosen amount of information desired, $\Delta$.

Recalling that at the first iteration of Algorithm <ref> we initialize the set of selected features to be empty, $X_A=\emptyset$, it follows that at the beginning $TE_{X_A\rightarrow Y}=0$. Recursively applying Equation <ref>, after the final iteration we have $TE_{X_A\rightarrow Y} = \sum_{k=1}^K TE_{X_k\rightarrow Y | X_{A_k} },$ where $K$ is the number of iterations of the while cycle and $X_k,X_{A_k}$ are respectively the selected feature and the set of selected features during the k-th iteration. Therefore, when the algorithm stops we have $TE_{X_A\rightarrow Y} \geq \frac{\Delta}{2B^2}$ (or $\frac{\Delta^2}{2}$ in classification). Therefore, the theorem follows by substituting these values in the reformulation of the ideal regression and classification error introduced in Remark <ref>

In the backward case, the following inequality holds from Theorem <ref>:

\begin{align*}
&\inf_{g\in \mathcal{G}_{\abar}}\E_{X,Y}[(Y_t-g(X_{\abar}^{\xtpast},Y_{\ytpast}))^2] \\&\leq \sigma^2+2B^2\cdot \sum_{k=1}^K TE_{X_{k}\rightarrow Y|X_{\abar_k}\setminus X_{k}},
\end{align*}

where $K$ is the number of iterations of the while cycle and $X_k,X_{\abar_k}$ are respectively the removed feature and the set of candidate selected features during the $k-th$ iteration. Substituting each transfer entropy term with the confidence interval of its estimator from Equation <ref>, and recalling that at the end of the algorithm 7.58$\sum_{k=1}^K \hat{TE}_{X_{k}\rightarrow Y|X_{\abar_k}\setminus X_{k}}\leq \frac{\delta}{2B^2}$, the result follows. Notice that, to guarantee that any possible subset of features contemporaneously satisfies the concentration inequality of Equation <ref> with probability $\eta$, a term $2^d$ is added in the concentration exploiting the union bound over all possible subsets of features.

Similarly, in the forward case, the following inequality holds from Theorem <ref>:

\begin{equation*}
\inf_{g\in \mathcal{G}_{A}}\E_{X,Y}[(Y_t-g(X_{A}^{\xtpast},Y_{\ytpast}))^2] \leq
\sum_{k=1}^K TE_{X_k\rightarrow Y | X_{A_k} },
\end{equation*}

where $K$ is the number of iterations of the while cycle and $X_k,X_{A_k}$ are respectively the selected feature and the set of selected features during the k-th iteration. Again, substituting each transfer entropy term with the high-probability confidence interval of its estimator from Equation <ref>, and recalling that at the end of the algorithm the cumulative selected information 7.58$\sum_{k=1}^K \hat{TE}_{X_k\rightarrow Y | X_{A_k} }\geq \frac{\delta}{2B^2}$, the result follows.

§ EXPERIMENTAL ANALYSIS: ADDITIONAL DETAILS AND RESULTS

In this section, the experiments performed both with synthetic data and with real-world datasets are discussed in detail. We ran the experiments on a server with 88 Intel(R) Xeon(R) CPU E7-8880 v4 @ 2.20GHz cpus and 94 GB of RAM. We used 8 cpus of the server and the overall computational time was less than 24 hours. As reported in the main paper, code and datasets can be found at <https://www.dropbox.com/scl/fo/ctvxqh8q131f8sim4s22o/h?dl=0   rlkey=5hsbpold82bi8vilu44ubb508> and they will be available on GitHub after the double-blind review process.

§.§ Synthetic Experiments

In the main paper, Figure <ref> shows the three causal graphs considered for the three, five and ten-dimensional synthetic experiments, together with the time lag and the coefficient selected to build the benchmark dataset. The graphs are designed to consider some classical challenges of time-series causal discovery: autocorrelation, time lags, noise, and spurious association due to common causal confounders. As a first experiment, in all three cases, $N=300$ samples have been drawn (corresponding to considering time series of length equal to $300$). In particular, each node without causal parents has been associated with a feature of the real-world experiment on droughts that considers fifteen features discussed in the next subsection ($X_0$ for the three-dimensional dataset, $X_1,X_3$ for the five-dimensional dataset, $X_0,X_5,X_7$ for the ten-dimensional one). In this way, a realistic time series is considered for each of these nodes. Then, following the causal dependencies, the other features and the target have been designed, considering a linear contribution of each of their (lagged) parent variables through the coefficient reported in Figure <ref>. Moreover, Gaussian noise with variance $\sigma=0.1$ has been added to each feature.
Considering a varying seed from $0$ to $9$, the algorithms have been run ten times to produce the results in terms of true positive rate (TPR) and false positive rate (TPR). In particular, forward and backward CMI feature selection have been considered as baselines, applying the algorithms provided in [Beraha et al., 2019]. Moreover, the PCMCI algorithm [Runge et al., 2019] has been considered as state-of-the-art causal discovery method, applying the implementation provided by the authors of the method in the Pyhton library Tigramite (<https://github.com/jakobrunge/tigramite.git>). Given the linear nature of the problem and the much faster computational time, PCMCI algorithm has been run considering linear independence tests based on partial correlations, leaving the other parameters with default values. Moreover, to detect the majority of causal links, the proposed algorithms Forward and Backward TEFS, have been run considering respectively a small value of the maximum information loss $\delta=10^{-6}$ and a large value of the required minimum information gain $\Delta=100$.

To further analyse the impact of the coefficients, the number of samples, the noise and the time lags, different values of each of the parameters have been considered, keeping the others fixed as discussed above and repeating the experiments for each configuration of value ten times, with seeds from $0$ to $9$. In particular, the coefficients have been randomly sampled five additional times from a uniform distribution in the interval $[-1,-0.5]\cup [0.5,1]$ (not considering the interval around zero $[-0.5,0.5]$ to avoid too weak causal links). The variance of the Gaussian noise has been increased five times, considering the interval $[0.1,0.3,0.5,0.7,0.9]$ for the three-dimensional setting and $[0.01,0.05,0.1,0.15,0.2]$ for the five and ten-dimensional ones (to avoid too noisy experiments). Then, the number of samples has been varied, considering five increasing values $[100,200,300,400,500]$. Finally, the maximum time lag that the algorithms should consider has been tested, considering the exact value $2$ together with $3,4,5$. As discussed in the main paper, in most cases the proposed methods are able to correctly detect the correct causal links, with TPR better or equivalent to the PCMCI algorithm, controlling the FPR better. Moreover, the proposed methods are sufficiently robust increase of noise, a decrease of the number of samples and with different causal links (identified by different random coefficients). Additionally, they have a satisfactory performance also considering different time depths, with a worsening of the TPR only with large depth.

In the main paper, Table <ref> shows the TPR and FPR for all the considered algorithms for the three-dimensional setting. In this section, Table <ref> and <ref> respectively show the same results for the five and ten-dimensional synthetic experiments. The TPR and FPR have been obtained considering the true causal links between the features and the target, for all ten repetitions of each experiment with different seeds. As an example, in the three-dimensional experiment, there are two correct and one false causal link between the features and the target (also considering its autoregressive component). Therefore, if the links are correctly detected during all ten repetitions, the TPR is equal to $\frac{20}{20}=1$. On the contrary, if an algorithm detects one incorrect causal link out of the repetitions, the FPR is equal to $\frac{1}{10}=0.1$. Finally, we consider the PCMCI to detect a link between a feature and the target if the algorithm detects one link between them, independently from the time lag that it selects. Indeed, PCMCI algorithm tries to identify not only the causal features but also the time lags that have causal effects on the actual target (see Remark <ref> in the main paper for a broader discussion on the choice to consider all the contribution of a feature as a whole for the forward and backward TEFS).

Five-dimensional experiments on synthetic datasets. Each experiment has been repeated ten times, with different seeds, varying one parameter at a time.

Forward CMI     Backward CMI     PCMCI     lightgray!25 Forward TEFS (ours)     lightgray!25 Backward TEFS (ours)
Benchmark                 lightgray!25     lightgray!25

TPR     $0.47$     $0.47$     $0.67$     lightgray!25 $0.93$     lightgray!25 $0.93$

FPR     $0.50$     $0.50$     $0.40$     lightgray!25 $0$     lightgray!25 $0$

Increasing Noise                 lightgray!25     lightgray!25

TPR     $[0.67,0.60,0.47,0.47,0.40]$     $[0.67,0.60,0.47,0.47,0.40]$      $[0.37, 0.43, 0.67, 0.80, 0.83]$     lightgray!25 $[1,1,0.93,0.87,0.77]$     lightgray!25 $[1,1,0.93,0.87,0.77]$

FPR     $[0.50,0.50,0.50,0.50,0.50]$     $[0.50,0.50,0.50,0.50,0.50]$     $[0.05, 0.05, 0.40, 0.65, 0.50]$     lightgray!25 $[0, 0, 0, 0, 0]$     lightgray!25 $[0, 0, 0, 0, 0]$

Increasing $\tau$                 lightgray!25     lightgray!25

TPR     $[0.67,0.33,0.30,0.27]$     $[0.67,0.33,0.30,0.27]$      $[0.37,0.03,0.63,0.73]$     lightgray!25 $[1,0.67,0.53,0.47]$     lightgray!25 $[1,0.67,0.67,0.67]$

FPR     $[0.50,0.50,0.50,0.50]$     $[0.50,0.50,0.50,0.50]$     $[0.05,0.05,0.60,0.60]$     lightgray!25 $[0,0.5,0,0]$     lightgray!25 $[0,0.5,0,0]$

Increasing N                 lightgray!25     lightgray!25

TPR     $[0.43, 0.57, 0.47, 0.63, 0.60]$     $[0.43, 0.57, 0.47, 0.63, 0.60]$     $[0.30,0.50,0.67,0.77,0.57]$     lightgray!25 $[0.70,0.90,0.93,1,0.97]$     lightgray!25 $[0.73,0.90,0.93,1,0.97]$

FPR     $[0.50,0,0.50,0.50,0.50]$     $[0.50,0,0.50,0.50,0.50]$     $[0,0.5,0.4,0.4,0.4]$     lightgray!25 $[0.2,0.05,0,0,0]$     lightgray!25 $[0.15,0.05,0,0,0]$

Random Coefficients                 lightgray!25     lightgray!25

TPR     $[0.43,0.67,0.67,0.43,0.67]$     $[0.43,0.67,0.67,0.43,0.67]$     $[0.63,0.77,0.77,0.90,0.77]$     lightgray!25 $[0.90,1,1,0.67,1]$     lightgray!25 $[0.90,1,1,0.67,1]$

FPR     $[0.50,0,0.50,0.50,0.55]$     $[0.50,0,0.50,0.50,0.50]$     $[0.15,0.30,0.45,0.30,0.45]$     lightgray!25 $[0,0.05,0,0,0]$     lightgray!25 $[0,0,0,0,0]$

Ten-dimensional experiments on synthetic datasets. Each experiment has been repeated ten times, with different seeds, varying one parameter at a time.

Forward CMI     Backward CMI     PCMCI     lightgray!25 Forward TEFS (ours)     lightgray!25 Backward TEFS (ours)
Benchmark                 lightgray!25     lightgray!25

TPR     $0.67$     $0.67$     $0.97$     lightgray!25 $1.0$     lightgray!25 $1.0$

FPR     $0.01$     $0.01$     $0.66$     lightgray!25 $0$     lightgray!25 $0$

Increasing Noise                 lightgray!25     lightgray!25

TPR     $[0.67, 0.67, 0.67, 0.67, 0.67]$     $[0.67, 0.67, 0.67, 0.67, 0.67]$     $[1, 0.97, 0.97, 0.77, 0.87]$     lightgray!25 $[1, 1, 1, 1, 1]$     lightgray!25 $[1, 1, 1, 1, 1]$

FPR     $[0,0,0.01, 0.09, 0.14]$     $[0,0,0.01,0.09,0.13]$     $[0.71, 0.57, 0.66, 0.57, 0.6]$     lightgray!25 $[0, 0, 0, 0, 0]$     lightgray!25 $[0, 0, 0, 0, 0]$

Increasing $\tau$                 lightgray!25     lightgray!25

TPR     $[0.67, 0.67, 0.67,0.43]$     $[0.67, 0.67, 0.67,0.43]$     $[0.53, 0.97, 0.97, 1]$     lightgray!25 $[1, 1, 1, 0.7]$     lightgray!25 $[1, 1, 1, 0.7]$

FPR     $[0.03, 0.01, 0, 0]$     $[0.06, 0.01, 0, 0]$     $[0.44, 0.66, 0.73, 0.83]$     lightgray!25 $[0.04, 0, 0, 0.01]$     lightgray!25 $[0.04, 0, 0, 0.01]$

Increasing N                 lightgray!25     lightgray!25

TPR     $[0.67,0.67,0.67,0.67,0.67]$     $[0.67,0.67,0.67,0.67,0.67]$     $[0.67, 0.9, 0.97, 0.9, 0.9]$     lightgray!25 $[0.6, 0.97, 1,1,1]$     lightgray!25 $[0.93, 1, 1,1,1]$

FPR     $[0.1,0.04,0.01,0.01,0]$     $[0.2, 0.04, 0.01, 0.01, 0]$     $[0.49,0.37,0.66,0.5,0.56]$     lightgray!25 $[0.03, 0.01, 0, 0, 0]$     lightgray!25 $[0.11, 0, 0, 0, 0]$

Random Coefficients                 lightgray!25     lightgray!25

TPR     $[0.67,0.63,0.67,0.36,0.67]$     $[0.67,0.63,0.67,0.36,0.67]$     $[0.97,0.77,0.87,0.97,0.93]$     lightgray!25 $[1,0.47,0,0.67,1]$     lightgray!25 $[1,0.6,0.97,0.67,1]$

FPR     $[0.01,0.14,0.23,0.19,0]$     $[0.01,0.14,0.16,0.16,0]$     $[0.66,0.36,0.54,0.16,0.63]$     lightgray!25 $[0,0,0,0,0]$     lightgray!25 $[0,0.03,0,0,0]$

To inspect also the capability of the proposed methods to identify causal relationships with a large number of features, we also repeated the synthetic analysis with $D\in\{15,20,40,60,80,100\}$, preserving the relationships and coefficients of the first nine features and the target as in Figure <ref> and adding independent variables that can only be autocorrelated or interact with each other. In particular, we iteratively add features in groups of three, where the first one is randomly sampled from a standard normal distribution, and the second one linearly depends (with coefficient randomly sampled from a uniform distribution in $[-1,1]$) on the value of the first of one timestep before plus Gaussian noise, and the third depends on the second of two timesteps before and on its value of three timesteps before (again, with random coefficients sampled from a uniform distribution in $[-1,1]$) plus Gaussian noise. This group of three variables increases the number of features without impacting the causal relationships between the features and the target. In this context, the experiments have been repeated ten times with fixed benchmark parameters as described for the ten-dimensional experiment, varying ten different seeds. The associated TPR and FPR are reported in Table <ref> and show the capability of the proposed methods to control the FPR the PCMCI method, which is less accurate in discarding features, selecting more than ten features in high dimensions. Moreover, from the result, it is also possible to notice that the forward version of the TEFS algorithm is more stable in high-dimensional contexts, due to the huge number of features in the conditioning phase of the first iterations of the backward approach, which makes more difficult the estimate of the contribution of each feature on the target.

Larger dimensional experiments on synthetic datasets. Each experiment has been repeated ten times, with different seeds.

Forward CMI     Backward CMI     PCMCI     lightgray!25 Forward TEFS (ours)     lightgray!25 Backward TEFS (ours)
$D=15$                 lightgray!25     lightgray!25

TPR     $0.67$     $0.67$     $0.73$     lightgray!25 $1.0$     lightgray!25 $1.0$

FPR     $0.01$     $0.01$     $0.52$     lightgray!25 $0$     lightgray!25 $0$

$D=20$                 lightgray!25     lightgray!25

TPR     $0.67$     $0.67$     $0.73$     lightgray!25 $1.0$     lightgray!25 $1.0$

FPR     $0.01$     $0.01$     $0.29$     lightgray!25 $0$     lightgray!25 $0$

$D=40$                 lightgray!25     lightgray!25

TPR     $0.67$     $0.67$     $0.73$     lightgray!25 $1.0$     lightgray!25 $1.0$

FPR     $0.01$     $0.01$     $0.25$     lightgray!25 $0$     lightgray!25 $0$

$D=60$                 lightgray!25     lightgray!25

TPR     $0.67$     $0.67$     $0.73$     lightgray!25 $1$     lightgray!25 $0.93$

FPR     $0.01$     $0.01$     $0.22$     lightgray!25 $0$     lightgray!25 $0.01$

$D=80$                 lightgray!25     lightgray!25

TPR     $0.67$     $0.57$     $0.73$     lightgray!25 $1.0$     lightgray!25 $0.90$

FPR     $0.01$     $0.20$     $0.20$     lightgray!25 $0$     lightgray!25 $0.18$

$D=100$                 lightgray!25     lightgray!25

TPR     $0.67$     $0.53$     $0.73$     lightgray!25 $1.0$     lightgray!25 $0.90$

FPR     $0.01$     $0.28$     $0.18$     lightgray!25 $0$     lightgray!25 $0.29$

§.§ Real-World Experiments

This subsection provides additional details on the real-world datasets and experiments performed in order to analyse the behaviour of the proposed algorithms in real-world scenarios. In particular, three datasets with five features have been first considered. These climate datasets are composed by climatological features and a scalar target. The targets represent the state of vegetation of three sub-basins of the Po River (Vegetation Health Index)[They have been extracted from "Peter Zellner. Vegetation health index - 231 m 8 days (version 1.0) [data set]. Eurac Research, 2022".]. The features are observational data of temperature and precipitation[Features have been retrieved from  "Richard C. Cornes, Gerard van der Schrier, Else J. M. van den Besselaar, and Philip D. Jones. An ensemble version of the e-obs temperature and precipitation data sets. Journal of Geophysical Research: Atmospheres, 123(17):9391–9409, 2018.", "K. Didan. Myd13q1 modis/aqua vegetation indices 16-day l3 global 250m sin grid v006 [data set]. NASA EOSDIS Land Processes DAAC., 2015." and averaged over the basin they refer to.]. Moreover, to consider a larger dataset, a similar dataset has been considered, with the same index as the target and fifteen features retrieved from the same sources and representing the precipitation and temperature over the sub-basin considered and its two neighbouring sub-basins. These four datasets have been designed by the authors and they are available in the repository of this work. It is expected that precipitation and temperature have an impact on the state of the vegetation, with some timesteps of lag. The timestep considered for each of the four datasets is a week and the four datasets contain $N=639$ train samples (from 2001 to 2014) and $N=228$ test samples (from 2015 to 2019).
In addition, to further analyse the behaviour of the proposed algorithms on a benchmark dataset, we consider a dataset composed of five variables available in an online causal repository (<https://causeme.uv.es/>), where we consider in turn each variable to be the target and the other four variables the candidate causal features (for a total of five experiments). In particular, in the repository mentioned above, we consider the dataset TESTWEATHNOISE with $N=1000$ samples ($600$ samples for training and the following $400$ to test) and five variables, that contains classical causal discovery challenges (autocorrelation, time delays, non-linearity, observational noise).

Together with the two proposed methods (Forward and Backward TEFS), as done for the synthetic experiments, we consider the forward and backward CMI feature selection described in [Beraha et al., 2019], together with the PCMCI algorithm [Runge et al., 2019], considering both linear and non-linear (CMI-based) independence tests. Additionally, to also consider a baseline method related to classical Granger causality, we also consider the FullCI method implemented in the same library (<https://github.com/jakobrunge/tigramite.git>), considering again both the linear and non-linear versions. Since there is no ground truth on the causal relationships between the features and the target, we consider the $R^2$ test score performing linear regression as an index of the performance of the algorithms (considering as inputs also the last $M$ values of the target for our approaches and when selected by the PCMCI and Full_CI algorithms). This is motivated by the nature of the proposed methodology, which is designed to select features through a causal quantity but finally aimed to select the most relevant features in terms of regression performance. Since the conditioning on the last values of the target itself is the basic concept to measure causality with transfer entropy, we also provide the $R^2$ test score considering only the autoregressive component. Moreover, we consider $L=M\in\{1,2,3\}$ as past timesteps considered by the algorithms, since in real-world problems there is no knowledge on the number of past values that have an impact on the current value of the target.

Real-world experiments for three climate datasets with five features and one climate dataset with fifteen features. Each experiment has been repeated three times, considering different time-depths $M=L\in\{1,2,3\}$. The table reports the set of IDs of the selected features and the corresponding $R^2$ test score.

3 c Climate, 5 features, first dataset

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.28$     $\{\}\rightarrow0.28$     $\{\}\rightarrow0.27$

Forward CMI     $\{0,3,4\}\rightarrow0.34$     $\{0,1\}\rightarrow0.36$     $\{0,1\}\rightarrow0.35$

Backward CMI     $\{0,1,2,3,4\}\rightarrow0.34$     $\{0,3,4\}\rightarrow0.26$     $\{0,1,3\}\rightarrow0.35$

FullCI_linear     $\{0,2,4\}\rightarrow0.45$     $\{0,1,3,4\}\rightarrow0.47$      $\{0,3,4\}\rightarrow0.48$

FullCI_nonLinear     $\{0,4\}\rightarrow0.44$     $\{0,4\}\rightarrow0.42$      $\{0,3,4\}\rightarrow0.48$

PCMCI_linear     $\{0,2,4\}\rightarrow0.45$     $\{0,2,4\}\rightarrow0.43$      $\{0,2,3,4\}\rightarrow0.47$

PCMCI_nonLinear     $\{0,1,2,4\}\rightarrow0.43$     $\{0,2,3,4\}\rightarrow0.42$      $\{0,1,3,4\}\rightarrow0.47$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{0,4\}\rightarrow0.44$       lightgray!25 $\{0,4\}\rightarrow0.42$     lightgray!25 $\{0,4\}\rightarrow0.47$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{0,4\}\rightarrow0.44$       lightgray!25 $\{0,4\}\rightarrow0.42$     lightgray!25 $\{0,4\}\rightarrow0.47$

3 c Climate, 5 features, second dataset

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.29$     $\{\}\rightarrow0.30$     $\{\}\rightarrow0.29$

Forward CMI     $\{0,1,2,3,4\}\rightarrow0.24$     $\{1,2,3,4\}\rightarrow0.23$     $\{2,3\}\rightarrow0.25$

Backward CMI     $\{0,1,2,3,4\}\rightarrow0.24$     $\{1,2,3,4\}\rightarrow0.23$     $\{2,3\}\rightarrow0.25$

FullCI_linear     $\{0,1,2\}\rightarrow0.41$     $\{0,2,4\}\rightarrow0.45$     $\{2,3,4\}\rightarrow0.46$

FullCI_nonLinear     $\{0,1,2\}\rightarrow0.41$     $\{2\}\rightarrow0.42$     $\{2\}\rightarrow0.47$

PCMCI_linear     $\{1,2,3\}\rightarrow0.43$     $\{0,1,2,3\}\rightarrow0.47$     $\{1,2,3\}\rightarrow0.47$

PCMCI_nonLinear     $\{1,2,3,4\}\rightarrow0.43$     $\{1,2,3,4\}\rightarrow0.40$     $\{1,2,3,4\}\rightarrow0.46$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{2\}\rightarrow0.43$      lightgray!25 $\{2\}\rightarrow0.42$     lightgray!25 $\{2\}\rightarrow0.47$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{2\}\rightarrow0.43$      lightgray!25 $\{2\}\rightarrow0.42$      lightgray!25 $\{2\}\rightarrow0.47$

3 c Climate, 5 features, third dataset

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.29$     $\{\}\rightarrow0.30$     $\{\}\rightarrow0.29$

Forward CMI     $\{0,1,2,3\}\rightarrow0.17$     $\{0,1,2,3\}\rightarrow0.19$     $\{0,1,2\}\rightarrow0.19$

Backward CMI     $\{0,1,2,3\}\rightarrow0.17$     $\{0,1,2,3\}\rightarrow0.19$     $\{2,3,4\}\rightarrow0.15$

FullCI_linear     $\{1,2,3,4\}\rightarrow0.44$     $\{0,1,2\}\rightarrow0.40$     $\{0,1,2\}\rightarrow0.45$

FullCI_nonLinear     $\{1,2,3\}\rightarrow0.40$     $\{1\}\rightarrow0.38$     $\{1\}\rightarrow0.43$

PCMCI_linear     $\{1,2\}\rightarrow0.40$     $\{1,2,3\}\rightarrow0.40$     $\{1,2\}\rightarrow0.44$

PCMCI_nonLinear     $\{1,2\}\rightarrow0.40$     $\{1,2,4\}\rightarrow0.41$     $\{1,2\}\rightarrow0.44$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{1,2\}\rightarrow0.40$     lightgray!25 $\{1\}\rightarrow0.38$     lightgray!25 $\{1\}\rightarrow0.43$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{1,2\}\rightarrow0.40$     lightgray!25 $\{1\}\rightarrow0.38$     lightgray!25 $\{1\}\rightarrow0.43$

3 c Climate, 15 features

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.28$     $\{\}\rightarrow0.28$     $\{\}\rightarrow0.27$

Forward CMI     $\{0,3,4,11,12\}\rightarrow0.29$     $\{0,2,11\}\rightarrow0.30$     $\{0,5\}\rightarrow0.28$

Backward CMI     $\{0,2,3,11\}\rightarrow0.29$     $\{0,3,4\}\rightarrow0.26$     $\{1,2,3,12\}\rightarrow0.26$

FullCI_linear     $\{2,5,7,12,14\}\rightarrow0.45$     $\{5,7,11\}\rightarrow0.47$     $\{12\}\rightarrow0.32$

FullCI_nonLinear     $\{12\}\rightarrow0.32$     $\{\}\rightarrow0.28$     $\{13\}\rightarrow0.29$

PCMCI_linear     $\{0,2,7\}\rightarrow0.45$     $\{0,2,7,9,13\}\rightarrow0.42$     $\{2,7\}\rightarrow0.48$

PCMCI_nonLinear     $\{2,7\}\rightarrow0.44$     $\{0,2,7,9,10\}\rightarrow0.41$     $\{0,7,12,13\}\rightarrow0.48$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{0,7\}\rightarrow0.44$      lightgray!25 $\{0,7\}\rightarrow0.43$     lightgray!25 $\{7\}\rightarrow0.47$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{0,11\}\rightarrow0.44$     lightgray!25 $\{0,11\}\rightarrow0.42$     lightgray!25 $\{11\}\rightarrow0.43$

Real-world experiments for benchmark dataset, with five variables, considering one variable at time as target. Each experiment has been repeated three times, considering different time-depths $M=L\in\{1,2,3\}$. The table reports the set of IDs of the selected features and the corresponding $R^2$ test score.

3 c Benchmark, variable 0 as target

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.10$     $\{\}\rightarrow0.12$     $\{\}\rightarrow0.12$

Forward CMI     $\{2\}\rightarrow0.01$     $\{3\}\rightarrow0.01$     $\{1\}\rightarrow0.01$

Backward CMI     $\{1,3\}\rightarrow0.01$     $\{3\}\rightarrow0.01$     $\{1\}\rightarrow0.01$

FullCI_linear     $\{\}\rightarrow0.10$     $\{\}\rightarrow0.12$      $\{\}\rightarrow0.12$

FullCI_nonLinear     $\{\}\rightarrow0.10$     $\{\}\rightarrow0.12$      $\{\}\rightarrow0.12$

PCMCI_linear     $\{\}\rightarrow0.10$     $\{\}\rightarrow0.12$      $\{\}\rightarrow0.12$

PCMCI_nonLinear     $\{\}\rightarrow0.10$     $\{\}\rightarrow0.12$      $\{\}\rightarrow0.12$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{\}\rightarrow0.10$      lightgray!25 $\{\}\rightarrow0.12$     lightgray!25 $\{\}\rightarrow0.12$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{\}\rightarrow0.10$       lightgray!25 $\{\}\rightarrow0.12$     lightgray!25 $\{\}\rightarrow0.12$

3 c Benchmark, variable 1 as target

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0$     $\{\}\rightarrow0$     $\{\}\rightarrow0$

Forward CMI     $\{0\}\rightarrow0$     $\{2,3\}\rightarrow0$     $\{3\}\rightarrow0$

Backward CMI     $\{2,3\}\rightarrow0$     $\{2,3\}\rightarrow0$     $\{3\}\rightarrow0$

FullCI_linear     $\{\}\rightarrow0$     $\{2,4\}\rightarrow0$     $\{2,4\}\rightarrow0$

FullCI_nonLinear     $\{\}\rightarrow0$     $\{\}\rightarrow0$     $\{4\}\rightarrow0$

PCMCI_linear     $\{\}\rightarrow0$     $\{2,4\}\rightarrow0$     $\{2,4\}\rightarrow0$

PCMCI_nonLinear     $\{\}\rightarrow0$     $\{2\}\rightarrow0$     $\{2\}\rightarrow0$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{\}\rightarrow0$      lightgray!25 $\{\}\rightarrow0$     lightgray!25 $\{\}\rightarrow0$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{\}\rightarrow0$      lightgray!25 $\{\}\rightarrow0$      lightgray!25 $\{\}\rightarrow0$

3 c Benchmark, variable 2 as target

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.01$     $\{\}\rightarrow0.01$     $\{\}\rightarrow0.03$

Forward CMI     $\{3,4\}\rightarrow0.09$     $\{3\}\rightarrow0.07$     $\{0,3\}\rightarrow0.07$

Backward CMI     $\{3,4\}\rightarrow0.09$     $\{3\}\rightarrow0.07$     $\{0,3\}\rightarrow0.07$

FullCI_linear      $\{3,4\}\rightarrow0.09$     $\{3,4\}\rightarrow0.13$     $\{0,3,4\}\rightarrow0.13$

FullCI_nonLinear      $\{3,4\}\rightarrow0.09$     $\{3,4\}\rightarrow0.13$     $\{3,4\}\rightarrow0.12$

PCMCI_linear      $\{3,4\}\rightarrow0.09$     $\{3,4\}\rightarrow0.13$     $\{3,4\}\rightarrow0.12$

PCMCI_nonLinear      $\{3,4\}\rightarrow0.09$     $\{3,4\}\rightarrow0.13$     $\{3,4\}\rightarrow0.12$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{3,4\}\rightarrow0.09$     lightgray!25 $\{4\}\rightarrow0.09$     lightgray!25 $\{3\}\rightarrow0.11$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{3,4\}\rightarrow0.09$     lightgray!25 $\{4\}\rightarrow0.09$     lightgray!25 $\{3\}\rightarrow0.11$

3 c Benchmark, variable 3 as target

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0.01$     $\{\}\rightarrow0.01$     $\{\}\rightarrow0.01$

Forward CMI     $\{0,4\}\rightarrow0.01$     $\{1,2\}\rightarrow0$     $\{2\}\rightarrow0$

Backward CMI     $\{2\}\rightarrow0$     $\{1,2\}\rightarrow0$     $\{2\}\rightarrow0$

FullCI_linear     $\{2,4\}\rightarrow0.01$     $\{2,4\}\rightarrow0.03$     $\{2,4\}\rightarrow0.03$

FullCI_nonLinear     $\{2\}\rightarrow0.01$     $\{\}\rightarrow0.01$     $\{1,2\}\rightarrow0.01$

PCMCI_linear     $\{2,4\}\rightarrow0.01$     $\{2,4\}\rightarrow0.03$     $\{2,4\}\rightarrow0.03$

PCMCI_nonLinear     $\{2\}\rightarrow0.01$     $\{2,4\}\rightarrow0.03$     $\{2,4\}\rightarrow0.03$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{4\}\rightarrow0.02$     lightgray!25 $\{\}\rightarrow0.01$     lightgray!25 $\{\}\rightarrow0.01$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{4\}\rightarrow0.02$     lightgray!25 $\{\}\rightarrow0.01$     lightgray!25 $\{\}\rightarrow0.01$

3 c Benchmark, variable 4 as target

$M=L=1$     $M=L=2$     $M=L=3$

Autoregression     $\{\}\rightarrow0$     $\{\}\rightarrow0$     $\{\}\rightarrow0$

Forward CMI     $\{2,3\}\rightarrow0$     $\{2,3\}\rightarrow0.01$     $\{2\}\rightarrow0.03$

Backward CMI     $\{2,3\}\rightarrow0$     $\{2,3\}\rightarrow0.01$     $\{2\}\rightarrow0.03$

FullCI_linear     $\{2,3\}\rightarrow0$     $\{2,3\}\rightarrow0.01$     $\{2,3\}\rightarrow0.03$

FullCI_nonLinear     $\{3\}\rightarrow0$     $\{2\}\rightarrow0.03$     $\{2,3\}\rightarrow0.03$

PCMCI_linear     $\{2,3\}\rightarrow0$     $\{2,3\}\rightarrow0.01$     $\{2,3\}\rightarrow0.03$

PCMCI_nonLinear     $\{\}\rightarrow0$     $\{2,3\}\rightarrow0.01$     $\{2,3\}\rightarrow0.03$

lightgray!25 Forward TEFS (ours)     lightgray!25 $\{2\}\rightarrow0$      lightgray!25 $\{2\}\rightarrow0.03$     lightgray!25 $\{2\}\rightarrow0.03$

lightgray!25 Backward TEFS (ours)     lightgray!25 $\{2\}\rightarrow0$     lightgray!25 $\{2\}\rightarrow0.03$     lightgray!25 $\{2\}\rightarrow0.03$

Table <ref> contains the subset of selected features and the $R^2$ test score with linear regression of the three climatological datasets with five features and the one with fifteen features, considering the three values of $L,M$ and the feature selection and causal discovery algorithms discussed above. From the results it is possible to see that the proposed algorithms almost always select the same set of features, showing a stable behaviour. Moreover, as already discussed, the two proposed methods tend to identify the subset of the most (causally) relevant features, preserving the majority of information, that is highlighted by the competitive results. PCMCI and FullCI algorithms have a less stable behaviour and tend to identify larget subsets of features, sometimes slightly improving the performance in terms of $R^2$ score, at the cost of considering more features. Finally, Table <ref> reports the subset of selected features and the associated $R^2$ test score with linear regression of the benchmark dataset, considering, in turn, each variable as a target and the others as candidate causal features. From the results, the first variable (variable 0) seems to be an autoregressive process. The second variable (variable 1) seems to be a random variable with no dependencies on its previous value or the previous values of any other feature. For the third variable (variable 2) all methods obtain similar performances, showing that the autoregressive component is negligible for this variable and variable 3 and 4 are the relevant ones for its prediction. Regarding the fourth variable (variable 3), it is mostly a noise signal, with the proposed approaches that do not identify significant features in the majority of cases, while PCMCI and its variants identify slight importance of variable 2,4. Finally, for the fifth variable (variable 4), considering one timestep before, the other variables and their previous value do not provide information for its prediction, while a little information seems to be provided by variable 2 considering two and three timesteps before.