On Selecting Stable Predictors in Time Series Models

TL;DR

This paper introduces a new predictor selection method for time series that accounts for data dependence, improving over traditional methods like lasso, with demonstrated benefits on simulated and real datasets.

Contribution

It extends feature selection to dependent data using mixing processes, providing finite-sample guarantees and practical applicability for time series analysis.

Findings

Improved predictor stability in dependent data settings

Finite-sample performance guarantees

Effective on both simulated and real datasets

Abstract

We extend the feature selection methodology to dependent data and propose a novel time series predictor selection scheme that accommodates statistical dependence in a more typical i.i.d sub-sampling based framework. Furthermore, the machinery of mixing stationary processes allows us to quantify the improvements of our approach over any base predictor selection method (such as lasso) even in a finite sample setting. Using the lasso as a base procedure we demonstrate the applicability of our methods to simulated and several real time series datasets.

Tables

Table 1. Table 1: RMSE/MAPE on the test period compared to competing methods. Here q 1 = ⌊ 0.2 ​ p ⌋ subscript 𝑞 1 0.2 𝑝 q_{1}=\lfloor 0.2p\rfloor and q 2 = ⌊ 0.6 ​ p ⌋ subscript 𝑞 2 0.6 𝑝 q_{2}=\lfloor 0.6p\rfloor .

Method	BeijingPM2.5	ShanghaiPM2.5	IndoorTemp	AirQuality
$q_1$-BPA	$1.946/294.159$	$2.052/44.844$	$1.197/4.345$	$1.875/4.487$
$q_1$-Lasso	$1.994/292.555$	$2.267/49.415$	$1.209/4.385$	$2.247/5.608$
$q_2$-BPA	$2.368/333.700$	$2.616/56.98$	$1.324/4.837$	$3.744/9.621$
$q_2$-Lasso	$2.424/343.126$	$2.873/62.607$	$1.446/5.403$	$4.116/10.455$
Lasso	$2.584/359.337$	$3.490/77.106$	$1.507/5.657$	$5.728/15.086$
ElasticNet	$2.561/357.331$	$3.487/76.937$	$1.508/5.657$	$5.716/15.01$
AdapLasso	$2.597/360.919$	$3.497/77.139$	$1.212/4.396$	$5.527/14.366$
AIC	$1.868/300.788$	$3.209/70.978$	$1.271/4.739$	$5.665/14.804$

Equations

$${\mathbf{y}}_t=\sum _{i=1}^{\tilde{p}}{\mathbf{A}}_i{\mathbf{y}}_{t-i}+\sum _{j=0}^s{\mathbf{C}}_j{\mathbf{x}}_{t-j}+{\mathbf{u}}_t$$

$$\underset{\mathbf{B}\in d\times (k\tilde{d}+ms)}{\min }\sum _{t=\max (\tilde{p},s)}^{T}L({\mathbf{y}}_{t+1},B{\mathbf{Z}}_t)+\lambda {\Vert \mathbf{B}\Vert }_1$$

$${\mathbf{Z}}_t={\left[\begin{array}{c}{\mathbf{y}}_t\dots {\mathbf{y}}_{t-\tilde{p}+1}{\mathbf{x}}_t\dots {\mathbf{x}}_{t-s+1}\end{array}\right]}^{\top }$$

$$\mathbf{B}=\left[\begin{array}{c}{\mathbf{A}}_1\dots {\mathbf{A}}_{\tilde{p}}{\mathbf{C}}_1\dots {\mathbf{C}}_s\end{array}\right]$$

$$p_{k,l_T}=\mathcal{P}(k\in {\hat{S}}_{l_T})=\mathbb{E}\left[1_{k\in {\hat{S}}_{l_T}}\right]$$

$$\mathbf{O}={\cup }_{j=1}^{{\mu }_T}{O}_j\text{,}$$

$$O_j=\{i:2(j-1)a_T+1\le i\le (2j-1)a_T\}$$

$$\mathbf{E}={\cup }_{j=1}^{{\mu }_T}{E}_j\text{,}$$

$$E_j=\{i:2(j-1)a_T+1\le i\le 2j{a}_T\}$$

$${\Pi }_B^{\text{av}}(k)=\frac{1}{2B}\sum _{j=1}^{2B}\text{mathbbm}{1}_{k\in {\hat{S}}_{|{\mathbf{O}}_j|}}$$

$${\Pi }_B^{\text{sim}}(k)=\frac{1}{B}\sum _{j=1}^B\text{mathbbm}{1}_{k\in {\hat{S}}_{|{\mathbf{O}}_{2j-1}|}}\text{mathbbm}{1}_{k\in {\hat{S}}_{|{\mathbf{O}}_{2j}|}}$$

$$\frac{T{c}_{\beta }}{C(\varphi ,B)}\max \left(\frac{1}{{\left(p_0/p\right)}^2},\frac{1}{{\left(1-p_0/p\right)}^2}\right)\le a_T^{r+1},r>0$$

$$C(\varphi ,B)=$$

$$\displaystyle\left\{\begin{array}[]{ll}\frac{1}{2\left(2\phi-1-1/(2B)\right)}\mbox{ \text{if}}\\ \phi\in\left(\min\left(\frac{1}{2}+\theta^{2},\frac{1}{2}+1/(2B)+\frac{3}{4}\theta^{2}\right),\frac{3}{4}\right]\\ \\ \frac{4\left(1-\phi+1/(2B)\right)}{1+1/B}\mbox{ \text{if} }\phi\in\left(\frac{3}{4},1\right]\end{array}\right.$$

$$\mathbb{E}\left[|{\hat{S}}^{\text{av}}\cap L_{\theta }|\right]$$

$$\mathbb{E}\left[|{\hat{N}}^{\text{av}}\cap H_{\theta }|\right]$$

$$2\theta C(\varphi ,B)\mathbb{E}\left[|{\hat{S}}_{l_T}\cap L_{\theta }|\right]$$

$$=\frac{2q^{2}C(\varphi ,B)}{p}=lp$$

$$Y_t=\sum _{i=1}^3{\theta }_i{Y}_{t-i}+\sum _{j=1}^{50}{\beta }_j{X}_t^j+{ϵ}_t\sim \mathcal{N}(0,{\sigma }^2)$$

$$\mathbf{Z}({\mathbf{O}}_j)=\{{\mathbf{Z}}_i,i\in {\mathbf{O}}_j\},\mathbf{Z}({\mathbf{E}}_j)=\{{\mathbf{Z}}_i,i\in {\mathbf{E}}_j\}$$

$$\mathbb{E}\left[h(\mathbf{O})\right]$$

$$\le \mathbb{E}\left[\left|h(\mathbf{O})-h(\tilde{\mathbf{O}})\right|\right]+\mathbb{E}\left[h(\tilde{\mathbf{O}})\right]$$

$$\mathbb{E}\left[h(\mathbf{O})\right]$$

$$0$$

$$=1-2{\Pi }_B^{\text{av}}(k)+{\Pi }_B^{\text{sim}}(k)$$

$$\mathcal{P}(k\in {\hat{S}}^{\text{av}})$$

$$\le \mathcal{P}\left({\Pi }_B^{\text{sim}}(k)\ge 2\varphi -1\right)$$

$$=\mathbb{E}\left[\text{mathbbm}{1}_{{\Pi }_B^{\text{sim}}(k)\ge 2\varphi -1}\right]$$

$$\le ({\mu }_T-1){\beta }_{a_T}+\mathbb{E}\left[h(\tilde{\mathbf{O}})\right]$$

$$=({\mu }_T-1){\beta }_{a_T}+{\mathcal{P}}_{\tilde{O}}\left({\Pi }_B^{\text{sim}}(k)\ge 2\varphi -1\right)$$

$$\le ({\mu }_T-1){\beta }_{a_T}+C(\varphi ,B){\mathbb{E}}_{\tilde{\mathcal{O}}}\left[{\Pi }_B^{\text{sim}}(k)\right]$$

$$=({\mu }_T-1){\beta }_{a_T}+C(\varphi ,B)p_{k,l_T}^2$$

$$C(\varphi ,B)=$$

$$\displaystyle\left\{\begin{array}[]{ll}\frac{1}{2\left(2\phi-1-1/(2B)\right)}\mbox{ \text{if}}\\ \phi\in\left(\min\left(\frac{1}{2}+\theta^{2},\frac{1}{2}+1/(2B)+\frac{3}{4}\theta^{2}\right),\frac{3}{4}\right]\\ \\ \frac{4\left(1-\phi+1/(2B)\right)}{1+1/B}\mbox{ \text{if} }\phi\in\left(\frac{3}{4},1\right]\end{array}\right.$$

$$\frac{T{c}_{\beta }{p}^2}{2C(\varphi ,B)p_0^2}\le a_T^{r+1}$$

$$\mathcal{P}(k\in {\hat{S}}^{\text{av}})\le 2C(\varphi ,B)p_{k,l_T}^2$$

$$\mathbb{E}\left[|{\hat{S}}^{\text{av}}\cap L_{\theta }|\right]$$

$$=\sum _{k=1}^p\mathcal{P}(k\in {\hat{S}}^{\text{av}})\text{mathbbm}{1}_{p_{k,l_T}\le \theta }$$

$$\le 2C(\varphi ,B)\sum _{k=1}^p{p}_{k,l_T}^2\text{mathbbm}{1}_{p_{k,l_T}\le \theta }$$

Full text

1 Introduction

Variable or Predictor selection is an important problem in machine learning and statistics and has received a lot of attention in the literature. These methods are valuable in various applications in biology and finance for both the very large $n$ (samples) or the high dimensional ($p>>n$) setting. In many cases of interest the objective of variable selection is to uncover an underlying causal mechanism and hence the method in question must exhibit robustness to noise in the data.

Most predictor selection methods discussed in the literature hold only for i.i.d data. The topic is relatively unexplored in the case of time series problems where variables correspond to the choice of exogenous predictors or even the endogenous lags to include in the regression model111We use predictors to mean either a (possibly lagged) exogenous predictor or a endogenous lag.. As time series data becomes pervasive with the advent of sensor enabled devices and networks of these devices (for e.g. Yu et al. (2017)), such predictor selection is poised to become widely applicable. Some recent applications already employ external time series of web searches to predict disease prevalence such as Influenza (GFT of Ginsberg et al. (2009), ARGO of Yang et al. (2015)) and twitter based sentiment time series for stock prediction Si et al. (2013). These approaches attempt to mitigate the problem of searching through millions of candidate exogenous time series that may be predictive of the underlying metric. In such applications selecting the correct time series is particularly crucial as with incorrect predictors, disease or other forecasts can deviate significantly from reality (see Lazer et al. (2014) for a discussion). Another important application lies in the domain of causal inference (Brodersen et al. (2015)), wherein the impact of an intervention on a metric is evaluated using a Bayesian structural time series model. This is achieved by forecasting, using metric history and exogenous predictors (pseudo-controls), beyond the time of intervention. Finally, comparing the difference between the forecasts and what really happened, the event’s impact is estimated. It is evident that the quality of the forecast is crucial in this task and that the choice of the pseudo-controls or the correct lags is very important and so we would like to have sufficient confidence in this choice.

1.1 Contribution

In this paper we propose and analyze novel and efficient stable procedures for predictor selection in time series inspired by the framework developed in Shah and Samworth (2013). To that end we first describe block sampling for time series and then propose stability measures called block pair average (BPA) and simultaneous block selection (SBS) with corresponding error control bounds. We empirically validate our procedures on several real time series and show both qualitative and quantitative improvements over competing methods. To the best of our knowledge these are the first predictor selection methods with finite sample guarantees that have applications to interpretable forecasting, classification and other time series domains.

1.2 Related Work

The foundation of our work is based on extensions to sub-sampling based methods for i.i.d data. In this setting Meinshausen and Bühlmann (2010) proposed Stability Selection (SS) as a repeated sub-sampling based methodology to improve the performance of any variable selection technique and also provided bounds on the number of selected false positives. This method is an improvement over standard variable selection techniques as it is usually non-trivial to provide error control bounds for methods run on the complete data. More recently Shah and Samworth (2013) proposed Complementary Pairs Stability Selection (CPSS) procedures that provide significant improvements over SS. Their performance bounds do not explicitly depend on signal and noise variables (that are usually unknown) but instead depends on the number of low selection probability variables that are included and on the number of high selection probability variables that are excluded by the CPSS procedures. This approach doesn’t require any dependence on restrictive and unverifiable assumptions such as exchangeability which is required for the analysis in Meinshausen and Bühlmann (2010). The central idea of stability selection using both SS and CPSS procedures is repeated execution of a base procedure (e.g. lasso - Tibshirani (1996) ) on subsamples of $\lfloor n/2\rfloor$ data points to identify variables that show up often in the selected set. In time series applications, the error control yielded by these stability procedures does not hold as the sub-sampling does not account for the underlying dependence.

Most existing predictor selection methods in time series are largely based on heuristics, Ng et al. (2013) or simply use plain lasso Yang et al. (2015); Buncic and Tischhauser (2017) on the entire data and it is non-trivial to provide guarantees for such methods. For the specific case of vector auto regression (VAR) models, Song and Bickel (2011) propose a grouped penalty based approach that provably identifies relevant lags and predictors in the asymptotic $d$ (number of time series) and $\tilde{p}$ (number of lags) regime. Our method is of a fundamentally distinct flavor in that we provide quantifiable improvement over any base predictor selection method, including the method in Song and Bickel (2011), even in the finite data (sample or dimension) setting. Moreover our approach also works for the more general VAR-X (VAR with exogenous) model and in general is independent of the base predictor selection mechanism.

2 Preliminaries

In this section we establish notation and preliminary details for the models we work with. Let ${\bf y}_{t}\in\mathbb{R}^{d}$ and ${\bf x}_{t}\in\mathbb{R}^{m}$ denote strictly stationary sequences. Before we present an example of variable selection procedure in the time series domain first consider the general VAR-X model Ltkepohl (2007)

[TABLE]

where each $\mathbf{A}_{i}\in\mathbb{R}^{d\times d}$ and $\mathbf{C}_{j}\in\mathbb{R}^{d\times m}$. Model (1) captures the effect of endogenous and exogenous variables at different lags on the response variable and has been a staple in econometrics with more recent applications in genetics Larvie et al. (2016) and renewable energy forecasting Cavalcante et al. (2016).

This model can be $\ell_{1}$-regularized and thus employed as a variable selection procedure, for example to determine stable cross sectional relations in ${\bf y}_{t}$ or to identify predictive exogenous series in ${\bf x}_{t}$.

[TABLE]

where $L$ is a convex loss and each ${\bf Z}_{t}\in\mathbb{R}^{d\tilde{p}+ms}$ and $\mathbf{B}\in\mathbb{R}^{d\times(dp+ms)}$ are stacked such that

[TABLE]

Since ${\bf y}_{t}$ and ${\bf x}_{t}$ is a stationary mixing sequence it follows that ${\bf Z}_{t}$ is also stationary (See Bradley (2005)). Unless stated we are not assuming Model 2 as the underlying data generating process (DGP). The model and its sparse estimation only serves as a means for predictor selection as described in the later sections.

3 Stable Predictor Selection

Consider the stationary sequence ${\bf Z}_{1},\ldots,{\bf Z}_{T}$ in $\mathbb{R}^{p}$ . In the case of Model 2, $p=d\tilde{p}+ms$. Let us define a variable selection procedure $\hat{S}_{l_{T}}=\hat{S}_{l_{T}}({\bf Z}_{i_{1}},\ldots,{\bf Z}_{i_{l_{T}}})$, as an estimator of the set of signal variables $S\subset\{1,...,p\}$, that takes as input a dependent random sample of length $l_{T}$ 222Corresponds to the proportion of the entire time series length ($T$) used in sampling. Specified later. and takes values in all subsets of $\{1,\ldots,p\}$. Define the selection probability of a variable index $k\in\{1,\ldots,p\}$ as

[TABLE]

The improved stability framework as presented in Shah and Samworth (2013) executes the base selection procedure $\hat{S}_{\lfloor n/2\rfloor}$ on $B$ i.i.d random samples of size $\lfloor n/2\rfloor$ to then combine the variable selection estimates. In our setting because the stationary sequence is dependent standard stability sub-sampling doesn’t apply, so instead using the independent block technique from Yu (1994) we create “almost" independent blocks and then transfer the stability error control to i.i.d blocks that have the same distribution as the original blocks.

Divide the sequence ${\bf Z}_{1},\ldots,{\bf Z}_{T}$ into $2\mu_{T}$ blocks of length $a_{T}$ and assume that $T=2\mu_{T}a_{T}$ without any loss of generality. Let $O$ and $E$ be the sets that denote the indices in the odd and even blocks respectively such that

[TABLE]

The intuition behind this approach is that for an appropriately mixing sequence the odd blocks are roughly independent, provided $a_{T}$ is large enough (See Figure 1). And as we show later in the analysis, if we create “independent" blocks with each new block having the same distribution as its dependent counterpart we can use this duality to work with the original dependent time series.

The block construction described leads to measures of stability for predictors that we define next.

3.1 Definitions

3.1.1 Block Pair Average (BPA)

Consider a randomly selected (without replacement) subset of the sequence of odd blocks $\mathbf{O}_{j}=\left(O_{j_{1}},\ldots,O_{j_{\lfloor\mu_{T}/2\rfloor}}\right)$.

Definition 1**.**

Block Pair Average: Let $\left\{\left(\mathbf{O}_{2j-1},\mathbf{O}_{2j}\right):j=1,\ldots\,B\right\}$ be the set of randomly selected pairs of sequence of blocks such that $\mathbf{O}_{2j-1}\cap\mathbf{O}_{2j}=\emptyset$. For some $\phi\in(0,1]$ the block average selection estimator is then defined as $\hat{S}^{\text{av}}=\{k:\Pi^{\text{av}}_{B}(k)\geq\phi\}$ where

[TABLE]

Algorithm 1 shows the use of the BPA measure in the context of lasso base predictor.

3.1.2 Simultaneous Block Selection (SBS)

Definition 2**.**

Simultaneous Block Selection: Let $\left\{\left(\mathbf{O}_{2j-1},\mathbf{O}_{2j}\right):j=1,\ldots\,B\right\}$ be the set of randomly selected pairs of sequence of blocks such that $\mathbf{O}_{2j-1}\cap\mathbf{O}_{2j}=\emptyset$. For some $\phi\in(0,1]$ the simultaneous selection estimator is then defined as $\hat{S}^{\text{sim}}=\{k:\Pi^{\text{sim}}_{B}(k)\geq\phi\}$ where

[TABLE]

Note that measures (6) and (7) are dependent data analogues of the i.i.d measures presented in Shah and Samworth (2013) and serve the same purpose of reducing variance (via averaging) as estimators of $p_{k,l_{T}}$. For these we have $l_{T}=\lfloor\mu_{T}/2\rfloor a_{T}$.

Finally to analyze these measures we need the following definitions

Definition 3**.**

For some $\theta\in(0,1]$, define $L_{\theta}=\{k:p_{k,l_{T}}\leq\theta\}$ be the set of low selection probability features under $\hat{S}_{l_{T}}$ and $H_{\theta}=\{k:p_{k,l_{T}}>\theta\}$ denote the set of features that have high selection probability and $\hat{N}_{l_{T}}=\{1,\ldots,p\}-\hat{S}_{l_{T}}$.

The quantities $\mathbb{E}\left[|\hat{S}^{\text{av}}\cap L_{\theta}|\right]$ and $\mathbb{E}\left[|\hat{N}^{\text{av}}\cap H_{\theta}|\right]$ denote the expected number of low (noise) and high (signal) probability predictors that are included and excluded by our procedure. We analyze these measures against a base predictor to quantify the improvements that our procedures yield.

3.2 Assumptions

We will refer to the following assumptions when necessary

Assumption 1: Probability of selection under $\hat{S}_{l_{T}}$ is better than random. Thus $p_{k,l_{T}}\geq p_{0}/p$ where $p_{0}$ is the number of signal variables. 2. 2.

Assumption 2: The process ${\bf Z}_{1},\ldots,{\bf Z}_{T}$ is algebraically $\beta$-mixing (See Mohri and Rostamizadeh (2009)) implying there exists some constant $c_{\beta}>0$ such that $\beta_{t}\leq c_{\beta}/t^{r}$ for $r>0$. 3. 3.

Assumption 3: The process ${\bf Z}_{1},\ldots,{\bf Z}_{T}$ is geometrically $\beta$-mixing (See Bradley (2005)) implying there exists some constant $c_{\beta}>0$ such that $\beta_{t}\leq\exp(-c_{\beta}t)$. 4. 4.

Assumption 4: The distribution of the sequence $\Pi^{\text{sim}}_{B}(k)$ is unimodal for each $k\in L_{\theta}$. This assumption is discussed in Section $3.2$ of Shah and Samworth (2013) and reflects the empirical observation that for many different DGPs the distribution of $\Pi^{\text{sim}}_{B}(k)$ is consistently unimodal and leads to a sharper version of Markov inequality employed in the analysis.

There are several equivalent definitions of $\beta$-mixing in the literature but we refer to Definition-$2.2$ in Yu (1994).

4 Block Pair Averaging - Analysis

For the BPA measure we have the following

Theorem 1**.**

If $\phi\in\{\frac{1}{2}+\frac{1}{B},\frac{1}{2}+\frac{3}{2B},...,1\}$ and $\theta<\frac{1}{\sqrt{3}}$, then under Assumptions 1, 2 and 4 and when

[TABLE]

where

[TABLE]

then we have for the block average procedure (6) and with $l_{T}=\lfloor\mu_{T}/2\rfloor,a_{T}=\lfloor T/4\rfloor$

[TABLE]

Remark 1- $C(\phi,B)$ is the constant in the improved Markov inequality that was developed for a finite grid ($\phi$) in Shah and Samworth (2013). 2. 2.

Remark 2 When $\phi=0.8$ and $\theta=0.2$ then we get $\mathbb{E}\left[|\hat{S}^{\text{av}}\cap L_{\theta}|\right]\leq 0.28\ \mathbb{E}\left[|\hat{S}_{l_{T}}\cap L_{\theta}|\right]$, showing that the BPA measure selects only at most $28\%$ of the low selection probability predictors selected by the base procedure. 3. 3.

Remark 3- For geometric mixing sequences (Assumption 3) we require $a_{T}=\mathcal{O}(\log T)$ for the Theorem 1 to hold. 4. 4.

Remark 4- The improved empirical performance obtained using the $r$-concavity assumption also carries over to our case but for simplicity and generality we retain only the unimodal assumption. See Shah and Samworth (2013) for more details. 5. 5.

Remark 5 We can quite easily extend the BPA measure to the case of $d$-BPA wherein we sample $d$ blocks without replacement to get more conservative measures. This becomes obvious in the steps leading upto (19) in the proof of Theorem 1. However, for simplicity we retain only the discussion for $d=2$. The i.i.d analogue of this approach seems non-trivial. 6. 6.

Remark 6 The SBS measure is used in the proof of Theorem 1 and may also be used in practice besides BPA and its analysis is almost exactly similar.

5 Empirical Results

To convince the reader of the utility of our methods we evaluate them on simulated and real time series data using lasso as a base predictor. For the real time series we evaluate forecast errors on a hold-out test period using the selected predictors from Algorithm 1 and compare them to other predictor selection methods such as standard lasso, elastic net Zou and Hastie (2005), adaptive lasso Zou (2006) and AIC based maximum lag selection for both $\tilde{p}$ and $s$ (See Model 1). For these methods we use the standard method of cross-validation. Note that the standard method of CV works for auto-regressive models as long as the errors are assumed to be uncorrelated. See Bergmeir et al. (2018) for a discussion.

The post selection training is restricted to least squares regression. The train-test split is $67\%/33\%$ and the predictions are made in a rolling fashion with each test data point added to the training set and re-trained to forecast the next step. On the simulated data we test the robustness of our method by adding varying degree of noise to the simulated data and reporting true positive and false positive rates (TPR/FPR). We de-trend the real data in case of any obvious violations of stationarity.

Before we get into the details of the empirical results we first discuss the choice of parameters and how to use Theorem 1 for the case of lasso as a base selection method.

5.1 Theory in Practice

Since $\mathbb{E}\left[|\hat{S}_{l_{T}}\cap L_{\theta}|\right]$ is not known in practice we approximate it by $q=\mathbb{E}\left[|\hat{S}_{l_{T}}|\right]$. Setting $q$ is equivalent to choosing the regularization parameter $\lambda$, we can vary it until lasso selects $q$ predictors. Once $q$ is set we use $\theta=q/p$ to denote the irrelevant variables as those having less than average selection probability.

The appropriate threshold can then be determined based on Theorem 1 by solving for $\phi$

[TABLE]

where $l\in(0,1]$ is pre-specified. In general we can specify any two of $q$, $\phi$ and $l$ to get the stable predictor set. We will present results with different values of $q$ for some threshold $\phi\in(0.5,0.9]$. $B$ is set to $50$ and $a_{T}$ was set to some multiple of seasonality in the data for all experiments. This reflects the empirical assumption that across seasons the data points have a weak dependence. In the absence of seasonality we didn’t see a significant difference against other choices of $a_{T}$ such as $\sqrt{T}$ or $\log T$. See Algorithm 1 for a complete recipe of our method using the lasso base predictor. In the algorithm $\mathbf{Y}(\mathbf{O}_{l})$ and $\mathbf{Z}(\mathbf{O}_{l})$ correspond to the data points from the corresponding block sequence (see also the proof of Theorem 1). The regularizer sequence is the commonly used in the glmnet package Friedman et al. (2010). Note that in the proofs we require $\lfloor T/4\rfloor$ of the data to be used for the comparison against the base procedure, however in Algorithm 1 and empirically we use all the data for determining the $q$-estimated active set.

5.2 Simulated Data

We simulate a synthetic time series with $T=10000$ from the following auto-regressive with exogenous variables (AR(3)-X) model where the exogenous time series are simulated from a simple AR($2$) model with Gaussian innovations.

[TABLE]

We also simulate $100$ time series (AR($3$)) for the dataset along with $2$ i.i.d noise time series to create a more realistic setting of several candidates with varying possibility of being included in the stable set. Our objective is to recover a stable set (lags and exogenous) using Algorithm 1 for a AR(3)-X model with $104$ exogenous time series. The parameters are normally distributed, $\beta_{i}\sim\mathcal{N}(0,0.05)$ for all $i=1,...,50$.

We vary the amount of noise $(\mathcal{N}(0,\sigma^{2}))$ added to the predictors and the response and plot the TPR/FPR with $1$-standard deviation error bars. To estimate the stable predictors using the BPA measure we set $q=40$ and $\phi=0.8$. Figure 2a) shows that the stable set is more robust to growing noise levels compared to the competing approaches. The TPR for lasso and Elastic net is higher since it selects many more predictors and consequently also has a much higher FPR. In contrast adaptive lasso is too conservative and Algorithm 1 has an FPR comparable to it but much better TPR. Note that for $\sigma^{2}=0$ only less than half of the original predictors in Model (14) show up in the selected set, since a pure least squares estimation of the model reveals that only a fraction of the $50$ predictors are statistically significant.

6 Real Data

6.1 Chinese Cities PM2.5

The Chinese cities PM2.5 dataset, PM (2) is an hourly sampled measurement of particulate matter concentration alongside several other exogenous meteorological time series. Our goal is to forecast the PM2.5 concentration of Beijing and Shanghai using exogenous and endogenous signals. The data specific parameters settings are $T=52854$, $a_{T}=3*24$ (thrice the daily seasonality), $\tilde{p}_{\text{max}}=24$ and $s_{\text{max}}=3$ (max endogenous and exogenous lags to select from).

6.2 Indoor Temperature

The indoor temperature dataset SML is a one minute sampled ($15$ minute smoothed) set of measurements by a monitoring system that captures various attributes such as CO2 concentration, humidity, lighting, outdoor conditions etc. The idea behind this dataset is to aid in accurate forecasting of the indoor temperature to better regulate the energy consumption of the HVAC system. For this analysis we only consider the impact of the exogenous signal on the temperature so that $\tilde{p}_{\text{max}}=0$ and $s_{\text{max}}=2$. Other parameters are $T=4137$, $a_{T}=96$.

6.3 Air Quality - Temperature Prediction

The air quality dataset air is an hourly sampled set of measurements by chemical sensors that capture various attributes such as CO, NO2 and O3 concentration. While the original purpose of this dataset is to estimate sensor estimation quality, we repurpose it to model temperature as a function of these concentrations. For this analysis we set $\tilde{p}_{\text{max}}=24$ and $s_{\text{max}}=3$. Other parameters are $T=9358$, $a_{T}=24$.

6.4 Results and Discussion

Table 1 reveals that the BPA measure performs better on all but one datasets. For the BeijingPM2.5 dataset it appears that all variables contribute to the improved RMSE. Note that the choice of $q$ does make a difference to the results, as too small a value would exclude signal variables (appears to be the case for BeijingPM2.5) and a large $q$ would add noisy variables to the selected set. Some prior knowledge of the domain can be useful here, but in its absence a largish $q$ (for example $\lfloor 0.5p\rfloor$) seems to be a safe choice (due to the guaranteed error control for the BPA measure) as the results indicate that for most datasets RMSE only degrades a little. From a causality perspective Figure 2b) reveals that dew point (temperature at which dew droplets form) and precipitation are important variables. This is corroborated by meteorological studies such as Liang et al. (2015).

7 Conclusion

Filling the gap in the feature selection literature for dependent data, we proposed novel stable predictor selection techniques in the time series setting. For robustness scores that can be used with any selection method, we provided theoretical results guaranteeing error control. The stable predictors selected by our method were shown to have superior predictive performance on several real datasets. In the future we intend to explore such a scheme and measures for nonstationary and heteroscedastic time series.

8 Appendix

We gather all the proofs in this section

8.1 Proof of Theorem 1

Proof.

Denote the sequence of random variables that correspond to the $\mathbf{O}_{j}$ and $\mathbf{E}_{j}$ indices as

[TABLE]

Consider a sequence of i.i.d blocks $\{{\bf\tilde{Z}}(\mathbf{O}_{j}):j=1,\ldots,\mu_{T}\}$ where ${\bf\tilde{Z}}(\mathbf{O}_{j})=\{{\bf\tilde{Z}}_{i}:i\in\mathbf{O}_{j}\}$ such that the sequence is independent of ${\bf Z}_{1},\ldots,{\bf Z}_{T}$ and each block has the same distribution as a block from the original sequence (${\bf Z}(\mathbf{O}_{j})$).

Let $h(\mathbf{O})$ be a bounded measurable function on the set of selected blocks and $\tilde{\mathbf{O}}$ be the corresponding i.i.d sequence of blocks from ${\bf\tilde{Z}}(\mathbf{O}_{j})$ then

[TABLE]

where the expectation is w.r.t to the distribution of the original and constructed sequences. For a measurable and bounded function $h$ such that $|h|\leq M$ we have from Lemma 4.1 of Yu (1994)

[TABLE]

Next we have

[TABLE]

Using (17) for $h(\mathbf{O})=\mathbbm{1}_{\Pi^{\text{sim}}_{B}(k)\geq 2\phi-1}$, using the unimodal Markov inequality for the simultaneous selector on the constructed independent blocks $\{\tilde{\mathbf{O}}\}$, from Theorem 3 of (Shah and Samworth, 2013) and finally using the fact that the sequence is stationary (all blocks have the same distribution)

[TABLE]

where

[TABLE]

Using the assumption $\left(p_{0}/p\right)\leq p_{k,a_{T}}$ (Assumption 1) it is easy to see that when $\mu_{T}(c_{\beta}/a^{r}_{T})\leq C(\phi,B)p^{2}_{k,l_{T}}$

[TABLE]

This implies from (19) that when (25) holds

[TABLE]

Following the proof of Theorem $1$ in Shah and Samworth (2013) and using (26) we have

[TABLE]

For the second part replace $\hat{S}^{\text{av}}$ by $\hat{N}^{\text{av}}$ as an estimator of $N$, the set of noise variables. Let $\Pi^{\text{av}}_{B,\hat{N}}$ and $\Pi^{\text{sim}}_{B,\hat{N}}$ correspond to the noise variable estimators. Let $h(\tilde{\mathbf{O}})=\mathbbm{1}_{\Pi^{\text{sim}}_{B,\hat{N}}}$

[TABLE]

Just like in the proof of the first part when $Tc_{\beta}/\left(2C(\phi,B)\left(1-p_{0}/p\right)^{2}\right)\leq a^{r+1}_{T}$

[TABLE]

[TABLE]

∎