Variable Screening for High Dimensional Time Series

TL;DR

This paper extends variable screening methods to high-dimensional linear models with dependent and heavy-tailed data, introducing a new GLSS procedure that leverages serial correlation for improved performance.

Contribution

It analyzes the theoretical properties of SIS under dependence and heavy tails, and proposes GLSS, a novel screening method utilizing serial correlation, with proven sure screening properties.

Findings

GLSS outperforms SIS in many dependent data scenarios.

Both procedures achieve sure screening under certain dependence and moment conditions.

Simulations and real data application demonstrate effectiveness of the proposed methods.

Abstract

Variable selection is a widely studied problem in high dimensional statistics, primarily since estimating the precise relationship between the covariates and the response is of great importance in many scientific disciplines. However, most of theory and methods developed towards this goal for the linear model invoke the assumption of iid sub-Gaussian covariates and errors. This paper analyzes the theoretical properties of Sure Independence Screening (SIS) (Fan and Lv [J. R. Stat. Soc. Ser. B Stat. Methodol. 70 (2008) 849-911]) for high dimensional linear models with dependent and/or heavy tailed covariates and errors. We also introduce a generalized least squares screening (GLSS) procedure which utilizes the serial correlation present in the data. By utilizing this serial correlation when estimating our marginal effects, GLSS is shown to outperform SIS in many cases. For both procedures we prove sure screening properties, which depend on the moment conditions, and the strength of dependence in the error and covariate processes, amongst other factors. Additionally, combining these screening procedures with the adaptive Lasso is analyzed. Dependence is quantified by functional dependence measures (Wu [Proc. Natl. Acad. Sci. USA 102 (2005) 14150-14154]), and the results rely on the use of Nagaev-type and exponential inequalities for dependent random variables. We also conduct simulations to demonstrate the finite sample performance of these procedures, and include a real data application of forecasting the US inflation rate.

Figures

Tables

Table 1. Table 1: Case 1

	SIS			GLSS
$(\gamma ,\alpha )$	(.4,.6)	(.5,.8)	(.6,.9)	(.4,.6)	(.5,.8)	(.6,.9)
Gaussian
$p_n=1000$	.95	.63	.15	.99	.99	.98
$p_n=5000$	.62	.11	.01	.95	.95	.97
$t_5$
$p_n=1000$	.58	.26	.06	.83	.84	.83
$p_n=5000$	.21	.01	0	.55	.49	.50

Table 2. Table 2: Case 2: Scenario A

	SIS			GLSS
$\alpha$	.4	.6	.8	.4	.6	.8
Gaussian
$p_n=1000$	.83	.73	.55	.95	.90	.90
$p_n=5000$	.38	.30	.07	.63	.63	.57
$t_5$
$p_n=1000$	.44	.42	.21	.56	.56	.53
$p_n=5000$	.01	.04	0	.16	.14	.16

Table 3. Table 3: Case 2: Scenario B

	SIS			GLSS
$\alpha$	.4	.6	.8	.4	.6	.8
Gaussian
$p_n=1000$	.90	.82	.68	.99	1.00	1.00
$p_n=5000$	.71	.64	.26	.95	.97	.98
$t_5$
$p_n=1000$	.76	.63	.40	.92	.90	.92
$p_n=5000$	.37	.26	.06	.76	.74	.75

Table 4. Table 4: Case 3: Scenario A

	GLSS-AdaLasso			AdaLasso
$\alpha$	.4	.5	.6	.4	.5	.6
Gaussian
$p_n=1000$	.79	.65	.49	.60	.49	.35
$p_n=5000$	.84	.65	.46	.66	.43	.29
$t_5$
$p_n=1000$	.45	.37	.23	.32	.22	.14
$p_n=5000$	.36	.32	.18	.24	.18	.10

Table 5. Table 5: Case 3: Scenario B

	GLSS-AdaLasso			AdaLasso
$\alpha$	.4	.5	.6	.4	.5	.6
Gaussian
$p_n=1000$	.86	.72	.59	.57	.49	.34
$p_n=5000$	.69	.59	.43	.60	.44	.25
$t_5$
$p_n=1000$	.48	.41	.22	.30	.19	.10
$p_n=5000$	.35	.25	.19	.25	.16	.11

Table 6. Table 6: Inflation Forecasts: 12 month horizon

	CPI-MSE	CPI-MAE	PPI-MSE	PPI-MAE
AR(4)	1.00	1.00	1.00	1.00
Lasso	.94	.99	.69	.89
Adaptive Lasso	1.08	1.05	.80	.99
SIS-Lasso	.96	.97	.76	.95
SIS-Adaptive Lasso	1.03	1.00	.82	1.00
GLSS-Lasso	.84	.98	.65	.87
GLSS-Adaptive Lasso	.94	1.00	.70	.92
AR(4) + 4 Factors	1.18	.99	1.08	1.09

Equations

$$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$$

$${ϵ}_i=g\left(\dots ,e_{i-1},e_i\right)$$

$${\delta }_q({ϵ}_i)=||{ϵ}_i-g\left({\mathcal{F}}_i^{\ast }\right)|{|}_q=(E|{ϵ}_i-g\left({\mathcal{F}}_i^{\ast }\right){|}^q{)}^{1/q}$$

$${\Delta }_{0,q}(\mathbf{ϵ})=\sum _{i=0}^{\infty }{\delta }_q({ϵ}_i)<\infty$$

$${\theta }_q({ϵ}_l)=||\text{E}\left({ϵ}_l|{\mathcal{F}}_0\right)-\text{E}\left({ϵ}_l|{\mathcal{F}}_{-1}\right)|{|}_q=||{\mathcal{P}}_0{ϵ}_l|{|}_q$$

$${{\mathbf{x}}_{\mathbf{i}}}^{(n)}=\mathbf{h}\left(\dots ,{\mathbf{\eta }}_{i-1}^{(n)},{\mathbf{\eta }}_i^{(n)}\right)$$

$${\Phi }_{m,q}(\mathbf{x})=\sum _{i=m}^{\infty }\underset{j\le p_n}{\max }{\delta }_q(X_{ij})<\infty$$

$$\hat{\mathbf{\rho }}=({\hat{\rho }}_1,\dots ,{\hat{\rho }}_{p_n})\text{, where }{\hat{\rho }}_j=(\sum _{t=1}^n{X}_{tj}^2{)}^{-1}(\sum _{t=1}^n{X}_{tj}{Y}_t)$$

$${\mathcal{M}}_{\ast }=\left\{1\le i\le p_n:{\beta }_i\ne 0\right\}$$

$${\hat{\mathcal{M}}}_{{\gamma }_n}=\left\{1\le i\le p_n:|{\hat{\rho }}_i|\ge {\gamma }_n\right\}$$

$${\delta }_q(X_{ij})=||A_{i,j}{\mathbf{\eta }}_0-A_{i,j}{\mathbf{\eta }}_0^{\ast }|{|}_q\le 2|A_{i,j}|||{\eta }_{0,1}|{|}_q$$

$${\tilde{\mathbf{x}}}_t$$

$${\vartheta }_n=\frac{s_n{n}^{\omega }{K}_{x,r}^r}{(n/s_n{)}^{r/2-r\kappa /2}}+\frac{n^{\iota }{K}_{x,r}^{\tau }{K}_{ϵ,q}^{\tau }}{n^{\tau -\tau \kappa }}+\exp \left(-\frac{n^{1-2\kappa }}{s_n^2{K}_{x,r}^4}\right)+\exp \left(-\frac{n^{1-2\kappa }}{K_{x,r}^2{K}_{ϵ,q}^2}\right)$$

$$P\left(\underset{j\le p_n}{\max }|{\hat{\rho }}_j-{\rho }_j|>c_2{n}^{-\kappa }\right)$$

$$P\left({\mathcal{M}}_{\ast }\subset {\hat{\mathcal{M}}}_{{\gamma }_n}\right)$$

$$P\left(|{\hat{\mathcal{M}}}_{{\gamma }_n}|\le O(n^{2\kappa }{\lambda }_{max}(\Sigma ))\right)$$

$$P\left(\underset{j\le p_n}{\max }|{\hat{\rho }}_j-{\rho }_j|>c_2{n}^{-\kappa }\right)\le$$

$$+O\left(p_n\exp {\left(-\frac{n^{1/2-\kappa }}{{\upsilon }_x{\upsilon }_{ϵ}}\right)}^{{\tilde{\alpha }}^{\prime }}\right)$$

$$P\left({\mathcal{M}}_{\ast }\subset {\hat{\mathcal{M}}}_{{\gamma }_n}\right)\ge 1$$

$$-O\left(s_n\exp {\left(-\frac{n^{1/2-\kappa }}{{\upsilon }_x{\upsilon }_{ϵ}}\right)}^{{\tilde{\alpha }}^{\prime }}\right)$$

$$P\left(|{\hat{\mathcal{M}}}_{{\gamma }_n}|\le O(n^{2\kappa }{\lambda }_{max}(\Sigma ))\right)\ge 1$$

$$-O\left(p_n\exp {\left(-\frac{n^{1/2-\kappa }}{{\upsilon }_x{\upsilon }_{ϵ}}\right)}^{{\tilde{\alpha }}^{\prime }}\right)$$

$$Y_t=X_{tk}{\rho }_k+{ϵ}_{t,k}$$

$${\tilde{\beta }}_k^M=({\mathbf{X}}_k^T{\Sigma }_k^{-1}{\mathbf{X}}_k{)}^{-1}{\mathbf{X}}_k^T{\Sigma }_k^{-1}\mathbf{y}$$

$${\hat{\Sigma }}_{k,l_n}={\left({\hat{\gamma }}_{i-j,k}\text{mathbbm}{1}_{|i-j|\le l_n}\right)}_{1\le i,j\le n}$$

$${\hat{\beta }}_k^M=({\mathbf{X}}_k^T{\hat{\Sigma }}_{k,l_n}^{-1}{\mathbf{X}}_k{)}^{-1}{\mathbf{X}}_k^T{\hat{\Sigma }}_{k,l_n}^{-1}\mathbf{y}$$

$${\beta }_k^M=E(y_t-\sum _{i=1}^{L_k}{\alpha }_i{y}_{t-i})(X_{t,k}-\sum _{i=1}^{L_k}{\alpha }_i{X}_{t-i,k})/(E(X_{t,k}-\sum _{i=1}^{L_k}{\alpha }_i{X}_{t-i,k}{)}^2)$$

$${\delta }_{q^{\prime }}({ϵ}_{t,k})=$$

$$\le \sum _{i\in M_{\ast }}|{\beta }_i|{\delta }_{q^{\prime }}(X_{ti})+{\delta }_{q^{\prime }}({ϵ}_t)+|{\rho }_k|{\delta }_{q^{\prime }}(X_{tk})$$

$$a_n=l_n\left[\frac{n^{\iota }{l}_n^{{\tau }^{\prime }}{K}_{x,r}^{{\tau }^{\prime }}{\tilde{K}}_{ϵ,q^{\prime }}^{{\tau }^{\prime }}}{n^{{\tau }^{\prime }-{\tau }^{\prime }\kappa }}+\frac{n^{\zeta }{l}_n^{q^{\prime }/2}{\tilde{K}}_{ϵ,q^{\prime }}^{q^{\prime }}}{n^{q^{\prime }/2-q^{\prime }\kappa /2}}+\frac{n^{\omega }{l}_n^{r/2}{K}_{x,r}^r}{n^{r/2}}\right]$$

$$b_n=l_n\left[\exp {\left(-\frac{n^{1/2}}{l_n{\upsilon }_x^2}\right)}^{\tilde{\alpha }}+\exp {\left(-\frac{n^{1/2-\kappa }}{l_n{\upsilon }_x\varphi }\right)}^{{\tilde{\phi }}^{\prime }}+\exp {\left(-\frac{n^{1/2-\kappa }}{l_n{\varphi }^2}\right)}^{\tilde{\phi }}\right]$$

$$P(||{\hat{\Sigma }}_{k,l_n}-{\Sigma }_k|{|}_2>c{n}^{-\kappa })\le O(a_n)$$

$$P(||{\hat{\Sigma }}_{k,l_n}-{\Sigma }_k|{|}_2>c{n}^{-\kappa })\le O(b_n)$$

Full text

62F07,

62J07,

High-dimensional Statistics,

sparsity,

Lasso,

Time Series,

functional dependence measure,

Variable Selection,

Nagaev inequality,

Sure Independence Screening,

journalname

,

1 Introduction

With the advancement of data acquisition technology, high dimensionality is a characteristic of data being collected in fields as diverse as health sciences, genomics, neuroscience, astronomy, finance, and macroeconomics. Applications where we have a large number of predictors for a relatively small number of observations are becoming increasingly common. For example, in disease classification we usually have thousands of variables, such as expression of genes, which are collected, while the sample size is usually in the tens. Other examples include fMRI data, where the number of voxels can number in the thousands and far outnumber the observations. For an overview of high dimensionality in economics and finance, see [22]. For the biological sciences, see [24, 5] and references therein. The main goals in these situations according to [4] are:

• •

To construct as effective a method as possible to predict future observations.

• •

To gain insight into the relationship between features and response for scientific purposes, as well as hopefully, to construct an improved prediction method.

More formally we are dealing with the case where

[TABLE]

with $\bm{y}=\left(Y_{1},\ldots,Y_{n}\right)^{T}$ being an $n$-vector of responses, $\bm{X}=\left(\bm{x_{1}},\ldots,\bm{x_{n}}\right)^{T}$ being an $n\times p_{n}$ random design matrix, and $\bm{\epsilon}=\left(\epsilon_{1},\ldots,\epsilon_{n}\right)^{T}$ is a random vector of errors. In addition, when the dimensionality of the predictors ($p_{n}$) is large we usually make the assumption that the underlying coefficient vector $(\bm{\beta})$ is sparse. Sparsity is a characteristic that is frequently found in many scientific applications [20],[31]. For example, in disease classification it is usually the case that only a small amount of genes are relevant to predicting the outcome.

Indeed, there are a wealth of theoretical results and methods that are devoted to this issue. Our primary focus is on screening procedures. Sure Independence Screening (SIS) as originally introduced in [20], was applicable to the linear model, and is based on a ranking of the absolute values of the marginal correlations of the predictors with the response. This method allows one to deal with situations in which the number of predictors is of an exponential order of the number of observations, which they termed as ultrahigh dimensionality. Further work on the topic has expanded the procedure to cover the case of generalized linear models [25], non-parametric additive models [19], Cox proportional hazards model [18], single index hazard rate models [27], and varying coefficient models [23]. Model-free screening methods have also been developed. For example; screening using distance correlation was analyzed in [35], a martingale difference correlation approach was introduced in [42], additional works include [57], [30] among others. For an overview of works related to screening procedures, one can consult [37]. The main result introduced with these methods is that, under appropriate conditions, we can reduce the predictor dimension from size $p_{n}=O\left(\exp\left(n^{\alpha}\right)\right)$, for some $\alpha<1$, to a size $d_{n}$, while retaining all the relevant predictors with probability approaching 1.

Another widely used class of methods is based on the penalized least squares approach. An overview of these methods is provided in [21] and [5]. Examples of methods in this class are the Lasso [45], and the adaptive Lasso [58]. Various theoretical results have been discovered for these class of methods. They broadly fall into analyzing the prediction error $|X(\hat{\beta}-\beta)|_{2}^{2}$, parameter estimation error $|\hat{\beta}-\beta|_{1}$, model selection consistency, as well as limiting distributions of the estimated parameters (see [9] for a comprehensive summary). Using screening procedures in conjunction with penalized least squares methods, such as the adaptive Lasso, presents a powerful tool for variable selection. Variable screening can allow us to quickly reduce the parameter dimension $p_{n}$ significantly, which weakens the assumptions needed for model selection consistency of the adaptive Lasso [29, 39].

A key limitation of the results obtained for screening methods, is the assumption of independent observations. In addition, it is usually assumed that the covariates and the errors are sub-Gaussian (or sub-exponential). However, there are many examples of real world data where these assumptions are violated. Data which is observed over time and/or space such as meteorological data, longitudinal data, economic and financial time series frequently exhibit covariates and/or errors which are serially correlated. One specific example is the case of fMRI time series, where there can exist a complicated spatial-temporal dependence structure in the errors and the covariates (see [48]). Another example is in forecasting macroeconomic indicators such as GDP or inflation rate, where we have large number of macroeconomic and financial time series, along with their lags, as possible covariates. Examples of heavy tailed and dependent errors and covariates can be found most prominently in financial, insurance and macroeconomic data.

These examples stress why it is extremely important for variable selection methods to be capable of handling scenarios where the assumption of independent sub-Gaussian (or sub-exponential) observations is violated. Some works related to this goal for the Lasso include [46], which extended the Lasso to jointly model the autoregressive structure of the errors as well as the covariates. However, their method is applicable only to the case where $p_{n}<n$, and they assume an autoregressive structure where the order of the process is known. Whereas [53] studied the theoretical properties of the Lasso assuming a fixed design in the case of heavy tailed and dependent errors. Additionally [3], and [33] investigated theoretical properties of the Lasso for high-dimensional Gaussian processes. Most recently [39] analyzed the adaptive Lasso for high dimensional time series while allowing for both heavy tailed covariate and errors processes, with the additional assumption that the error process is a martingale difference sequence.

Some works related to this goal for screening methods include [36], which allows for heavy tailed errors and covariates. Additionally [10], [54], and [57] also relax the Gaussian assumption, with the first two requiring the tails of the covariates and the response to be exponentially light, while the latter allows for heavy tailed errors provided the covariates are sub-exponential. Although these works relax the moment and distributional assumptions on the covariates and the response, they still remain in the framework of independent observations. A few works have dealt with correlated observations in the context of longitudinal data (see [13],[56]). However, the dependence structure of longitudinal data is too restrictive to cover the type of dependence present in most time series. Most recently [12] proposed a non-parametric kernel smoothing screening method applicable to time series data. In their work they assume a sub-exponential response, covariates that are bounded and have a density, as well as assuming the sequence $\{(Y_{i},\bm{x}_{i})\}$ is strong mixing, with the additional assumption that the strong mixing coefficients decay geometrically. These assumptions can be quite restrictive; they exclude, for example, heavy tailed time series, and discrete valued time series which are common in fields such as macroeconomics, finance, neuroscience, amongst others [16].

In this work, we study the theoretical properties of SIS for the linear model with dependent and/or heavy tailed covariates and errors. This allows us to substantially increase the number of situations in which SIS can be applied. However, one of the drawbacks to using SIS in a time series setting is that the temporal dependence structure between observations is ignored. In an attempt to correct this, we introduce a generalized least squares screening (GLSS) procedure, which utilizes this additional information when estimating the marginal effect of each covariate. By using GLS to estimate the marginal regression coefficient for each covariate, as opposed to OLS used in SIS, we correct for the effects of serial correlation. Our simulation results show the effectiveness of GLSS over SIS, is most pronounced when we have strong levels of serial correlation and weak signals. Using the adaptive Lasso as a second stage estimator after applying the above screening procedures is also analyzed. Probability bounds for our combined two stage estimator being sign consistent are provided, along with comparisons between our two stage estimator and the adaptive Lasso as a stand alone procedure.

Compared to previous work, we place no restrictions on the distribution of the covariate and error processes besides existence of a certain number of finite moments. In order to quantify dependence, we rely on the functional dependence measure framework introduced by [49], rather than the usual strong mixing coefficients. Comparisons between functional dependence measures and strong mixing assumptions are discussed in section 2. For both GLSS and SIS, we present the sure screening properties and show the range of $p_{n}$ can vary from the high dimensional case, where $p_{n}$ is a power of $n$, to the ultrahigh dimensional case discussed in [20]. We detail how the range of $p_{n}$ and the sure screening properties are affected by the strength of dependence and the moment conditions of the errors and covariates, the strength of the underlying signal, and the sparsity level, amongst other factors.

The rest of the paper is organized as follows: Section 2 reviews the functional and predictive dependence measures which will allow us to characterize the dependence in the covariate ($\bm{x_{i}},i=1,...,n$) and error processes. We also discuss the assumptions placed on structure of the covariate and error processes; these assumptions are very mild, allowing us to represent a wide variety of stochastic processes which arise in practice. Section 3 presents the sure screening properties of SIS under a range of settings. Section 4 introduces the GLSS procedure and presents its sure screening properties. Combining these screening procedures with the adaptive Lasso will discussed in Section 5. Section 6 covers simulation results, while section 7 discusses an application to forecasting the US inflation rate. Lastly, concluding remarks are in Section 8, and the proofs for all the results follow in the appendix.

2 Preliminaries

We shall assume the error sequence is a strictly stationary, ergodic process with the following form:

[TABLE]

Where $g(\cdot)$ is a real valued measurable function, and $e_{i}$ are iid random variables. This representation includes a very wide range of stochastic processes such as linear processes, their non-linear transforms, Volterra processes, Markov chain models, non-linear autoregressive models such as threshold auto-regressive (TAR), bilinear, GARCH models, among others (for more details see [50],[49]). This representation allows us to use the functional and predictive dependence measures introduced in [49]. The functional dependence measure for the error process is defined as the following:

[TABLE]

where $\mathcal{F}_{i}^{*}=\left(\ldots,e_{-1},e_{0}^{*},e_{1},\ldots,e_{i}\right)$ with $e_{0}^{*},e_{j},j\in\mathbb{Z}$ being iid. Since we are replacing $e_{0}$ by $e_{0}^{*}$, we can think of this as measuring the dependency of $\epsilon_{i}$ on $e_{0}$ as we are keeping all other inputs the same. The cumulative functional dependence measure is defined as $\Delta_{m,q}(\bm{\epsilon})=\sum_{i=m}^{\infty}\delta_{q}(\epsilon_{i})$. We assume weak dependence of the form:

[TABLE]

The predictive dependence measure is related to the functional dependence measure, and is defined as the following:

[TABLE]

where $\mathcal{F}_{i}=\left(\ldots,e_{-1},e_{0},e_{1},\ldots,e_{i}\right)$ with $e_{i},i\in\mathbb{Z}$ being iid. The cumulative predictive dependence measure is defined as $\Theta_{q}(\bm{\epsilon})=\sum_{l=0}^{\infty}\theta_{q}(\epsilon_{l})$, and by Theorem 1 in [49] we obtain $\Theta_{q}(\bm{\epsilon})\leq\Delta_{0,q}(\bm{\epsilon})$.

Similarly the covariate process is of the form:

[TABLE]

Where $\bm{\eta}_{i}^{(n)}\in\mathcal{R}^{p_{n}},i\in\mathbb{Z}$, are iid random vectors, $\bm{h(\cdot)}=(h_{1}(\cdot)\ldots,h_{p_{n}}(\cdot))$, $\bm{x_{i}}^{(n)}=(X_{i1},...,X_{ip_{n}})$ and $X_{ij}=h_{j}(...,\bm{\eta}_{i-1}^{(n)},\bm{\eta}_{i}^{(n)})$. The superscript $(n)$ denotes that the dimension of vectors is a function of $n$, however for presentational clarity we suppress the superscript $(n)$ from here on and use $\bm{x}_{i}$ and $\bm{\eta}_{i}$ instead. Let $\mathcal{H}_{i}^{*}=(\ldots,\bm{\eta}_{-1},\bm{\eta}^{*}_{0},\bm{\eta}_{1},\ldots,\bm{\eta}_{i})$. As before the functional dependence measure is $\delta_{q}(X_{ij})=||X_{ij}-h_{j}\left(\mathcal{H}_{i}^{*}\right)||_{q}$ and the cumulative dependence measure for the covariate process is defined as:

[TABLE]

The representations (2), and (6), along with the functional and predictive dependence measures have been used in various works including [52],[55], and [53] amongst others. Compared to strong mixing conditions, which are often difficult to verify, the above dependence measures are easier to interpret and compute since they are related to the data generating mechanism of the underlying process [50]. In many cases using the functional dependence measure also requires less stringent assumptions. For example, consider the case of a linear process, $\epsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j}$, with $e_{i}$ iid. Sufficient conditions for a linear process to be strong mixing involve: the density function of the innovations ($e_{i}$) being of bounded variation, restrictive assumptions on the decay rate of the coefficients ($f_{j}$), and invertibility of the process (see Theorem 14.9 in [14] for details). Additional conditions are needed to ensure strong mixing if the innovations for the linear process are dependent [17].

As a result many simple processes can be shown to be non-strong mixing. A prominent example involves an AR(1) model with iid Bernoulli (1/2) innovations: $\epsilon_{i}=\rho\epsilon_{i-1}+e_{i}$ is non-strong mixing if $\rho\in(0,1/2]$ [2]. These cases can be handled quite easily in our framework, since we are not placing distributional assumptions on the innovations, $e_{i}$, such as the existence of a density. For linear processes with iid innovations, representation (2) clearly holds and (4) is satisfied if $\sum_{j=0}^{\infty}|f_{j}|<\infty$. For dependent innovations, suppose we have: $e_{i}=h(\ldots,a_{i-1},a_{i})$, where $h(\cdot)$ is a real valued measurable function and $a_{i},i\in\mathbb{Z}$, are iid. Then $\epsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j}$, has a causal representation, and satisfies (4) if: $\sum_{i=0}^{\infty}\delta_{q}(e_{i})<\infty$, and $\sum_{j=0}^{\infty}|f_{j}|<\infty$ (see [51]).

3 SIS with Dependent Observations

Sure Independence Screening, as introduced by Fan and Lv [20], is a method of variable screening based on ranking the magnitudes of the $p_{n}$ marginal regression estimates. Under appropriate conditions, this simple procedure is shown to possess the sure screening property. The method is as follows, let:

[TABLE]

Therefore, $\hat{\rho}_{j}$ is the OLS estimate of the linear projection of $Y_{t}$ onto $X_{tj}$. Now let

[TABLE]

and let $|\mathcal{M}_{*}|=s_{n}<<n$ be the size of the true sparse model. We then sort the elements of $\bm{\hat{\rho}}$ by their magnitudes. For any given $\gamma_{n}$, define a sub-model

[TABLE]

and let $|\hat{\mathcal{M}}_{\gamma_{n}}|=d_{n}$ be the size of the selected model. The sure screening property states that for an appropriate choice of $\gamma_{n}$, we have $P\left(\mathcal{M}_{*}\subset\hat{\mathcal{M}}_{\gamma_{n}}\right)\rightarrow 1$.

Throughout this paper let: $Y_{t}=\sum_{i=1}^{p_{n}}X_{ti}\beta_{i}+\epsilon_{t}$, $\bm{x_{t}}=(X_{t1},...,X_{tp_{n}})$, $\Sigma=cov(\bm{x_{t}})$, and $\bm{X}_{k}$ be $k^{th}$ column of $\bm{X}$. In addition, we assume $Var(Y_{t}),Var(X_{tj})=O(1)$, $\forall j\leq p_{n}$. Note that $\bm{x}_{t}$ can contain lagged values of $Y_{t}$. Additionally, let $\rho_{j}=(E(X_{tj}^{2}))^{-1}E(X_{tj}Y_{t})$, and $\mathcal{M}_{\gamma_{n}}=\left\{1\leq i\leq p:|\rho_{i}|\geq\gamma_{n}\right\}$. For a vector $\bm{a}=(a_{1},...,a_{n})$, $sgn(\bm{a})$ denotes its sign vector, with the convention that $sgn(0)=0$, and $|\bm{a}|_{p}^{p}=\sum_{i=1}^{n}|a_{i}|^{p}$. For a square matrix $\bm{A}$, let $\lambda_{\min}(\bm{A})$ and $\lambda_{\max}(\bm{A})$, denote the minimum eigenvalue, and maximum eigenvalue respectively. For any matrix $\bm{A}$, let $||\bm{A}||_{\infty}$, and $||\bm{A}||_{2}$ denote the maximum absolute row sum of $\bm{A}$, and the spectral norm of $\bm{A}$ respectively. Lastly we will use $C,c$ to denote generic positive constants which can change between instances.

3.1 SIS with dependent, heavy tailed covariates and errors

To establish sure screening properties, we need the following conditions:

Condition A: $|\rho_{k}|\geq c_{1}n^{-\kappa}$ for $k\in M_{*}\textrm{, }\kappa<1/2$

Condition B: $E(\epsilon_{t}),E(X_{tj}),E(X_{tj}\epsilon_{t})=0\textrm{ }\forall j,t$.

Condition C: Assume the error and the covariate processes have representations (2), and (6) respectively. Additionally, we assume the following decay rates $\Phi_{m,r}(\bm{x})=O(m^{-\alpha_{x}}),\Delta_{m,q}(\bm{\epsilon})=O(m^{-\alpha_{\epsilon}})$, for some $\alpha_{x},\alpha_{\epsilon}>0$, $q>2,r>4$ and $\tau=\frac{qr}{q+r}>2$.

Condition A is standard in screening procedures, and it assumes the marginal signals of the active predictors cannot be too small. Condition B assumes the covariates and the errors are contemporaneously uncorrelated. This is significantly weaker than independence between the error sequence and the covariates usually assumed. Condition C presents the structure, dependence and moment conditions on the covariate and error processes. Notice that higher values of $\alpha_{x},\alpha_{\epsilon}$ indicate weaker temporal dependence.

Examples of error and covariate processes which satisfy Condition C are: If $\epsilon_{i}$ is a linear process, $\epsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j}$ with $e_{i}$ iid and $\sum_{j=0}^{\infty}|f_{j}|<\infty$ then $\delta_{q}(\epsilon_{i})=|f_{i}|||e_{0}-e_{0}^{*}||_{q}$. If $f_{i}=O(i^{-\beta})$ for $\beta>1$ we have $\Delta_{m,q}=O(m^{-\beta+1})$ and $\alpha_{\epsilon}=\beta-1$. We have a geometric decay rate in the cumulative functional dependence measure, if $\epsilon_{i}$ satisfies the geometric moment contraction (GMC) condition, see [41]. Conditions needed for a process to satisfy the GMC condition are given in Theorem 5.1 of [41]. Examples of processes satisfying the GMC condition include stationary, causal finite order ARMA, GARCH, ARMA-GARCH, bilinear, and threshold autoregressive processes, amongst others (see [50] for details).

For the covariate process, if we assume $\bm{x}_{i}$ is a vector linear process: $\bm{x}_{i}=\sum_{l=0}^{\infty}A_{l}\bm{\eta}_{i-l}$. Where $A_{l}$ are $p_{n}\times p_{n}$ coefficient matrices and $\bm{\eta}_{i}=(\eta_{i1},\ldots,\eta_{ip_{n}})$ are iid random vectors with $cov(\bm{\eta_{i}})=\Sigma_{\eta}$. For simplicity, assume $\eta_{i,j}(j=1,\ldots,p_{n})$ are identically distributed, then

[TABLE]

where $A_{i,j}$ is the $j^{th}$ column of $A_{i}$. If $||A_{i}||_{\infty}=O(i^{-\beta})$ for $\beta>1$, then $\Phi_{m,q}=O(m^{-\beta+1})$.

In particular for stable VAR(1) processes, $\bm{x}_{t}=B_{1}\bm{x}_{t-1}+\bm{\eta}_{t}$, $\Phi_{m,q}(\bm{x})=O(||\tilde{B}_{1}||_{2}^{m})$ [11]. For stable VAR($k$) processes, $\bm{x}_{t}=\sum_{i=1}^{k}B_{i}\bm{x}_{t-i}+\bm{\eta}_{t}$, we can rewrite this as a VAR(1) process, $\tilde{\bm{x}}_{t}=\tilde{B}_{1}\tilde{\bm{x}}_{t-1}+\tilde{\bm{\eta}}_{t}$, with:

[TABLE]

And by section 11.3.2 in [38], the process $\tilde{\bm{x}}_{t}$ is stable if and only if $\bm{x}_{t}$ is stable. Therefore if $\tilde{B}_{1}$ is diagonalizable, we have $O(a^{m})$, where $a$ represents the largest eigenvalue in magnitude of $\tilde{B}_{1}$. And by the stability of $\bm{x}_{t}$, $a\in(0,1)$. Additional examples of error and covariate processes which satisfy Condition C are given in [52] and [53] respectively.

Define $\alpha=min(\alpha_{x},\alpha_{\epsilon})$ and let $\omega=1$ if $\alpha_{x}>1/2-2/r$, otherwise $\omega=r/4-\alpha_{x}r/2$. Let $\iota=1$ if $\alpha>1/2-1/\tau$, otherwise $\iota=\tau/2-\tau\alpha$. Additionally, let $K_{\epsilon,q}=\sup_{m\geq 0}(m+1)^{\alpha_{\epsilon}}\Delta_{m,q}(\bm{\epsilon})$ and $K_{x,r}=\max_{j\leq p_{n}}\sup_{m\geq 0}(m+1)^{\alpha_{x}}\sum_{i=m}^{\infty}\delta_{r}(X_{ij})$. Given Condition C, it follows that $K_{\epsilon,q},K_{x,r}<\infty$. For ease of presentation we let:

[TABLE]

The following theorem gives the sure screening properties, and provides a bound on the size of the selected model:

Theorem 1**.**

Suppose Conditions A,B,C hold.

1. (i)

For any $c_{2}>0$, we have:

[TABLE] 2. (ii)

For $\gamma_{n}=c_{3}n^{-\kappa}$ with $c_{3}\leq c_{1}/2$, we have:

[TABLE] 3. (iii)

For $\gamma_{n}=c_{3}n^{-\kappa}$ with $c_{3}\leq c_{1}/2$, we have:

[TABLE]

In Theorem 1 we have two types of bounds, for large $n$ the polynomial terms dominate, whereas for small values of $n$ the exponential terms dominate. The covariate dimension ($p_{n}$) can be as large as $o(min(\frac{s_{n}(n/s_{n})^{r/2-r\kappa/2}}{n^{\omega}},\frac{n^{\tau-\tau\kappa}}{n^{\iota}}))$. The range of $p_{n}$ depends on the dependence in both the covariate and the error processes, the strength of the signal ($\kappa$), the moment condition, and the sparsity level ($s_{n}$). If we assume $s_{n}=O(1),r=q$, and $\alpha\geq 1/2-2/r$ then $p_{n}=o(n^{r/2-r\kappa/2-1})$. For the case of iid errors and covariates, we would replace $K_{x,r},K_{\epsilon,q}$ in Theorem 1 with $\max_{j\leq p_{n}}||X_{ij}||_{r/2}$ and $||\epsilon_{i}||_{q}$ respectively. Therefore for the case of weaker dependence in the covariate and error processes (i.e. $\alpha_{x}>1/2-2/r$ and $\alpha>1/2-1/\epsilon$), our range for $p_{n}$ is reduced only by a constant factor. However, our range for $p_{n}$ is significantly reduced in the case of stronger dependence in the error or covariate processes (i.e. either $\alpha_{x}<1/2-2/r$ or $\alpha_{\epsilon}<1/2-2/q$). For instance if $\alpha_{x}=\alpha_{\epsilon}$ and $q=r$, our range for $p_{n}$ is reduced by a factor of $n^{r/4-\alpha r/2}$ in the case of stronger dependence.

In the iid setting, to achieve sure screening in the ultrahigh dimensional case, [20] assumed the covariates and errors are jointly normally distributed. Future works applicable to the linear model, such as [25],[19] among others, relaxed this Gaussian assumption, but generally assumed the tails of the covariates and errors are exponentially light. Compared to the existing results for iid observations, our moment conditions preclude us from dealing with the ultrahigh dimensional case. However, our setting is far more general in that it allows for dependent and heavy tailed covariates and errors. In addition, we allow for the covariates and error processes to be dependent on each other, with the mild restriction that $E(X_{tj}\epsilon_{t})=0$, $\forall j\leq p_{n}$.

3.2 Ultrahigh Dimensionality under dependence

It is possible to achieve the sure screening property in the ultrahigh dimensional setting with dependent errors and covariates. However, we need to make stronger assumptions on the moments of both the error and covariate processes. Until now we have assumed the existence of a finite $qth$ moment, which restricted the range of $p$ to a power of $n$. If the error and covariate processes are assumed to follow a stronger moment condition, such as $\Delta_{0,q}(\bm{\epsilon})<\infty$ and $\Phi_{0,q}(\bm{x})<\infty$ for arbitrary $q>0$, we can achieve a much larger range of $p_{n}$ which will cover the ultrahigh dimensional case discussed in [20]. More formally, we have:

Condition D: Assume the error and the covariate processes have representations (2), and (6) respectively. Additionally assume $\upsilon_{x}=\sup_{q\geq 2}q^{-\tilde{\alpha}_{x}}\Phi_{0,q}(\bm{x})<\infty$ and $\upsilon_{\epsilon}=\sup_{q\geq 2}q^{-\tilde{\alpha}_{\epsilon}}\Delta_{0,q}(\bm{\epsilon})<\infty$, for some $\tilde{\alpha}_{x},\tilde{\alpha}_{\epsilon}\geq 0.$

By Theorem 3 in [53], Condition D implies the tails of the covariate and error processes are exponentially light. There are a wide range of processes which satisfy the above condition. For example, if $\epsilon_{i}$ is a linear process: $\epsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j}$ with $e_{i}$ iid and $\sum_{l=0}^{\infty}|f_{l}|<\infty$ then $\Delta_{0,q}(\epsilon_{l})=||e_{0}-e_{0}^{*}||_{q}\sum_{l=0}^{\infty}|f_{l}|$. If we assume $e_{0}$ is sub-Gaussian, then $\tilde{\alpha}_{\epsilon}=1/2$, since $||e_{0}||_{q}=O(\sqrt{q})$. Similarly if $e_{i}$ is sub-exponential we have $\tilde{\alpha}_{\epsilon}=1$. More generally, for $\epsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j}^{p}$, if $e_{i}$ is sub-exponential, we have $\tilde{\alpha}_{\epsilon}=p$. Similar results hold for vector linear processes discussed previously.

Condition D is primarily a restriction on the rate at which $||\epsilon_{i}||_{q},\max_{j\leq p_{n}}||X_{ij}||_{q}$ increase as $q\rightarrow\infty$. We remark that, for any fixed $q$, we are not placing additional assumptions on the temporal decay rate of the covariate and error processes besides requiring $\Delta_{0,q}(\bm{\epsilon})$,$\Phi_{0,q}(\bm{x})<\infty$. In comparison, in the ultrahigh dimensional setting, [12] requires geometrically decaying strong mixing coefficients, in addition to requiring sub-exponential tails for the response. As an example, if we assume $\epsilon_{i}=\sum_{j=0}^{\infty}f_{j}e_{i-j}$, geometrically decaying strong mixing coefficients would require the coefficients, $f_{j}$, to decay geometrically. Whereas in Condition D, the only restrictions we place on the coefficients, $f_{j}$, is absolute summability.

Theorem 2**.**

Suppose Conditions A,B,D hold. Define $\tilde{\alpha}^{\prime}=\frac{2}{1+2\tilde{\alpha}_{x}+2\tilde{\alpha}_{\epsilon}}$, and $\tilde{\alpha}=\frac{2}{1+4\tilde{\alpha}_{x}}$.

1. (i)

For any $c_{2}>0$ we have:

[TABLE] 2. (ii)

For $\gamma_{n}=c_{3}n^{-\kappa}$ with $c_{3}\leq c_{1}/2$, we have:

[TABLE] 3. (iii)

For $\gamma_{n}=c_{3}n^{-\kappa}$ with $c_{3}\leq c_{1}/2$, we have:

[TABLE]

From Theorem 2, we infer the covariate dimension ($p_{n}$) can be as large as

$o(\min[\exp\left(\frac{Cn^{1/2-\kappa}}{s_{n}}\right)^{\tilde{\alpha}}/s_{n},\exp(Cn^{1/2-\kappa})^{\tilde{\alpha}^{\prime}}])$. As in Theorem 1, the range of $p_{n}$ depends on the dependence in both the covariate and the error processes, the strength of the signal ($\kappa$), the moment condition, and the sparsity level ($s_{n}$). For the case of iid covariates and errors, we would replace $\upsilon_{x}$ and $\upsilon_{\epsilon}$ with $\mu_{r/2}=\max_{j\leq p_{n}}||X_{ij}||_{r/2}$ and $||\epsilon_{i}||_{q}$ respectively. In contrast to Theorem 1, temporal dependence affects our range of $p_{n}$ only by a constant factor.

If we assume $s_{n}=O(1)$, and both the covariate and error processes are sub-Gaussian we obtain $p_{n}=o(\exp(n^{\frac{1-2\kappa}{3}}))$, while for sub-exponential distributions we obtain $p_{n}=o(\exp(n^{\frac{1-2\kappa}{5}}))$. In contrast, Fan and Lv [20], assuming independent observations, allow for a larger range $p_{n}=o(\exp(n^{1-2\kappa})$. However, their work relied critically on the Gaussian assumption. Fan and Song [25], relax the Gaussian assumption by allowing for sub-exponential covariates and errors, and our rates are similar to theirs up to a constant factor. Additionally, in our work we relax the sub-exponential assumption, provided the tails of the covariates and errors are exponentially light.

4 Generalized Least Squares Screening (GLSS)

Consider the marginal model:

[TABLE]

where $\rho_{k}$ is the linear projection of $y_{t}$ onto $X_{tk}$. In SIS, we rank the magnitudes of the OLS estimates of this projection. In a time series setting, if we are considering the marginal model (14) it is likely the case that the marginal errors ($\epsilon_{t,k}$) will be serially correlated. This holds even if we assume that the errors ($\epsilon_{t}$) in the full model (1) are serially uncorrelated. A procedure which accounts for this serial correlation, such as Generalized Least Squares (GLS), will provide a more efficient estimate of $\rho_{k}$.

We first motivate our method by considering a simple univariate model. Assume $Y_{t}=\beta X_{t}+\epsilon_{t}$ and the errors follow an AR(1) process, $\epsilon_{t}=\rho\epsilon_{t-1}+\theta_{t}$, where $\theta_{t}$, and $X_{t}$ are iid standard Gaussian. We set $\beta=.5$, $n=200$, and estimate the model using both OLS and GLS for values of $\rho$ ranging from .5 to .95. The mean absolute errors for both procedures is plotted in figure 1. We observe that the performance of OLS steadily deteriorates for increasing values of $\rho$, while the performance of GLS stays constant. This suggests that a screening procedure based on GLS estimates will be most useful in situations where we have weak signals and high levels of serial correlation.

The infeasible GLS estimate for $\rho_{k}$ is:

[TABLE]

Where $\bm{X}_{k}$ is the $k^{th}$ column of $\bm{X}$, and $\Sigma_{k}=(\gamma_{i-j,k})_{1\leq i,j\leq n}$ is the auto-covariance matrix of $\bm{\epsilon}^{k}=(\epsilon_{t,k},t=1,...,n)$. Given that $\Sigma_{k}$ needs to be estimated to form our GLS estimates, we use the banded autocovariance matrix estimator introduced in [52], which is defined as:

[TABLE]

Where $l_{n}$ is our band length, $\hat{\gamma}_{r,k}=\frac{1}{n}\sum_{t=1}^{n-|r|}\hat{\epsilon}_{t,k}\hat{\epsilon}_{t+|r|,k}$, with $\hat{\epsilon}_{t,k}=y_{t}-X_{tk}\hat{\rho}_{k}$, and $\hat{\rho}_{k}$ is the OLS estimate of $\rho_{k}$. Our GLS estimator is now:

[TABLE]

When $E(\bm{\epsilon}^{k}|\bm{X}_{k})=0$, by the Gauss-Markov theorem it is clear that $\tilde{\beta}_{k}^{M}$ is efficient relative to the OLS estimator. [1] showed that under non-stochastic regressors and appropriate conditions on the error process, a two stage sieve type GLS estimator has the same limiting distribution as the infeasible GLS estimator $\tilde{\beta}_{k}^{M}$. In the appendix, we provide the appropriate conditions under which our GLS estimator, $\hat{\beta}_{k}^{M}$, and the infeasible GLS estimate, $\tilde{\beta}_{k}^{M}$, have the same asymptotic distribution.

For positive definite $\Sigma_{k}$, the banded estimate for $\Sigma_{k}$ is not guaranteed to be positive definite, however it is asymptotically positive definite (see Lemma 1). For small samples, we can preserve positive definiteness by using the tapered estimate: $\hat{\Sigma}_{k}*R_{l_{n}}$, where $R_{l_{n}}$ is a positive definite kernel matrix, and $*$ denotes coordinate-wise multiplication. For example, we can choose $R_{l_{n}}=(\max(1-\frac{|i-j|}{l_{n}},0))_{1\leq i,j\leq n}$. We need the following conditions for the sure screening property to hold:

Condition E: Assume the marginal error process, $\epsilon_{t,k}$, is a stationary AR($L_{k}$) process, $\epsilon_{t,k}=\sum_{i=1}^{L_{k}}\alpha_{i}\epsilon_{t-i,k}+e_{t}$. Where $L_{k}<K<\infty$, $\forall k\leq p_{n}$.

Condition F: For $k\in M_{*}\textrm{, }\kappa<1/2$: $\beta_{k}^{M}=E(y_{t}-\sum_{i=1}^{L_{k}}\alpha_{i}y_{t-i})(X_{t,k}-\sum_{i=1}^{L_{k}}\alpha_{i}X_{t-i,k})/(E(X_{t,k}-\sum_{i=1}^{L_{k}}\alpha_{i}X_{t-i,k})^{2})\geq c_{6}n^{-\kappa}$.

Condition G: Assume $E(X_{tk}),E(\epsilon_{t}),E(\bm{X}_{k}^{T}\Sigma_{k}^{-1}\bm{\epsilon})=0$

Condition H: Assume $\epsilon_{t,k}$, $\epsilon_{t}$ are of the form (2), and the covariate process is of the form (6). Additionally we assume the following decay rates $\Delta_{m,q}(\bm{\epsilon})=O(m^{-\alpha_{\epsilon}})$, $\Phi_{m,r}(\bm{x})=O(m^{-\alpha_{x}}),\chi_{m,q^{\prime}}=\sum_{i=m}^{\infty}\max_{k\leq p_{n}}\delta_{q}(\epsilon_{i,k})=O(m^{-\alpha})$, for some $\alpha_{x},\alpha_{\epsilon}>0$, $\alpha=min(\alpha_{x},\alpha_{\epsilon})$, and $q^{\prime}=\min(q,r)\geq 4$.

Condition I: Assume $\epsilon_{t,k}$, $\epsilon_{t}$ are of the form (2), and the covariate process is of the form (6). Additionally assume $\upsilon_{x}=\sup_{q\geq 4}q^{-\tilde{\alpha}_{x}}\Phi_{0,q}(\bm{x})<\infty$

,$\upsilon_{\epsilon}=\sup_{q\geq 4}q^{-\tilde{\alpha}_{\epsilon}}\Delta_{0,q}(\bm{\epsilon})<\infty,\phi=\sup_{q\geq 4}q^{-\varphi}\chi_{0,q}<\infty$ for some $\tilde{\alpha}_{x},\tilde{\alpha}_{\epsilon}\geq 0$, and $\varphi=\max(\tilde{\alpha}_{\epsilon},\tilde{\alpha}_{x})$.

In Condition E, we can let the band length $K$ diverge to infinity at a slow rate, e.g $O(\log(n))$, for simplicity we set $K$ to be a constant. Assuming a finite order AR model for the marginal error process is reasonable in most practical situations, since any stationary process with a continuous spectral density function can be approximated arbitrarily closely by a finite order linear AR process (see corollary 4.4.2 in [7]). For further details on linear AR approximations to stationary processes, see [1] and [8]. We remark that compared to previous works [1, 34], knowledge about the structure of the marginal errors is not necessary in estimating $\beta_{k}^{M}$, since we use a non-parametric estimate of $\Sigma_{k}$. Therefore Condition E is assumed strictly for technical reasons.

For Condition F, from (14), we have $\beta_{k}^{M}=\rho_{k}$, iff $E(\epsilon_{t,k}-\sum_{i=1}^{L_{k}}\alpha_{i}\epsilon_{t-i,k})(X_{t,k}-\sum_{i=1}^{L_{k}}\alpha_{i}X_{t-i,k})=0$. When $\beta_{k}^{M}\neq\rho_{k}$, recall that:

[TABLE]

If we assume the cross covariance, $\gamma_{\bm{X}_{k},Y}(h)$, is proportional to $E(X_{tk}Y_{t})$, i.e. $\gamma_{\bm{X}_{k},Y}(h)\propto E(X_{tk}Y_{t})$, for $h\in\{-L_{k},\ldots,-1,1,\ldots,L_{k}\}$, then $\beta_{k}^{M}\propto\rho_{k}$ whenever $|\beta_{k}^{M}|>0$. And for $|\rho_{k}|>0$, it is likely the case that $\beta_{k}^{M}\propto\rho_{k}$ if we assume $\gamma_{\bm{X}_{k},Y}(h)\propto E(X_{tk}Y_{t}),\textrm{ for }h\in\{-L_{k},\ldots,-1,1,\ldots,L_{k}\}$. When $\beta_{k}^{M}\neq\rho_{k}$, we believe the advantage in using GLSS is due to the GLS estimator being robust to serial correlation in the marginal error process (see the appendix for details).

For Condition H, since $\epsilon_{t,k}=Y_{t}-X_{tk}\rho_{k}$, we have $\epsilon_{t,k}=r_{k}(\ldots,\bm{\theta}_{t-1},\bm{\theta}_{t})$, where $r_{k}(\cdot)$ is a measurable function and $\bm{\theta}_{t}=(\bm{\eta}_{t},e_{t})$. If we assume $e_{t}$, and $\bm{\eta}_{i}$ are independent for $i\neq t$, then $\bm{\theta}_{i}$ are iid. We then have:

[TABLE]

Therefore, $\chi_{m,q^{\prime}}=O(m^{-\alpha})$, if we assume $\sum_{i\in M_{*}}|\beta_{i}|=O(1)$, $\Delta_{m,q}(\bm{\epsilon})=O(m^{-\alpha_{\epsilon}})$, and $\Phi_{m,r}(\bm{x})=O(m^{-\alpha_{x}})$.

For GLSS; define $\hat{\mathcal{M}}_{\gamma_{n}}=\left\{1\leq i\leq p_{n}:|\hat{\beta}_{k}^{M}|\geq\gamma_{n}\right\}$, $\alpha=min(\alpha_{x},\alpha_{\epsilon})$,

$\tau=\frac{qr}{q+r},\tau^{\prime}=\frac{qq^{\prime}}{q+q^{\prime}}=\min(q/2,\tau)$. Let $\iota=1$ if $\alpha>1/2-1/\tau^{\prime}$, otherwise $\iota=\tau^{\prime}/2-\tau^{\prime}\alpha$. Let $\zeta=1$, if $\alpha>1/2-2/q^{\prime}$, otherwise $\zeta=q^{\prime}/4-\alpha q^{\prime}/2$ and let $\omega=1$, if $\alpha_{x}>1/2-2/r$, otherwise $\omega=r/4-\alpha_{x}r/2$. Additionally, let $K_{x,r}=\max_{j\leq p_{n}}\sup_{m\geq 0}(m+1)^{\alpha_{x}}\sum_{i=m}^{\infty}\delta_{r}(X_{ij})$, $\tilde{K}_{\epsilon,q^{\prime}}=\max_{k\leq p_{n}}\sup_{m\geq 0}(m+1)^{\alpha}\sum_{i=m}^{\infty}\delta_{q^{\prime}}(\epsilon_{i,k})$. Given Condition H, it follows that $K_{x,r},\tilde{K}_{\epsilon,q^{\prime}}<\infty$. For the case of exponentially light tails, we define $\tilde{\varphi}^{\prime}=\frac{2}{1+2\tilde{\alpha}_{x}+2\varphi},\tilde{\varphi}=\frac{2}{1+4\varphi}$, and $\tilde{\alpha}=\frac{2}{1+4\tilde{\alpha}_{x}}$. Lastly, for ease of presentation let:

[TABLE]

We first present the following lemma, which provides deviation bounds on $||\hat{\Sigma}_{k,l_{n}}-\Sigma_{k}||_{2}$. This lemma, which is of independent interest, will allow us to obtain deviation bounds on our GLSS estimates.

Lemma 1**.**

Assume the band length, $l_{n}=c\log(n)$ for sufficiently large $c>0$.

1. (i)

Assume Condition H holds. For $\kappa\in[0,1/2)$ we have the following:

[TABLE] 2. (ii)

Assume Condition I holds. For $\kappa\in[0,1/2)$ we have the following:

[TABLE]

The following theorem gives the sure screening properties of GLSS:

Theorem 3**.**

Assume the band length, $l_{n}=c\log(n)$ for sufficiently large $c>0$.

1. (i)

Assume Conditions E,F,G,H hold, for any $c_{2}>0$ we have:

[TABLE] 2. (ii)

Assume Conditions E,F,G,H hold, then for $\gamma_{n}=c_{5}n^{-\kappa}$ with $c_{5}\leq c_{6}/2$:

[TABLE] 3. (iii)

Assume Conditions E,F,G,I hold, for any $c_{2}>0$ we have:

[TABLE] 4. (iv)

Assume Conditions E,F,G,I hold, then for $\gamma_{n}=c_{5}n^{-\kappa}$ with $c_{5}\leq c_{6}/2$:

[TABLE]

In Lemma 1, the rate of decay also depends on the band length ($l_{n}$). The band length primarily depends on the decay rate of the autocovariances of the process $\epsilon_{t,k}$. Since we are assuming an exponential decay rate, we can set $l_{n}=O(\log(n))$. If $\gamma_{i,k}=O(i^{-\beta})$ for $\beta>1$, then we require $l_{n}^{-\beta+1}=o(n^{-\kappa})$. We omit the exponential terms in the bounds for part part (i) of Lemma 1, and parts (i), and (ii) of Theorem 3 to conserve space and provide a cleaner result. For GLSS, the range for $p_{n}$ also depends on the band length ($l_{n}$), in addition to the moment conditions and the strength of dependence in the covariate and error processes. For example, if we assume $r=q$, and $\alpha\geq 1/2-2/r$ then $p_{n}=o(n^{r/2-r\kappa/2-1}/l_{n}^{r/2+1})$. Compared to SIS, we have a lower range of $p_{n}$ by a factor of $l_{n}^{r/2+1}$. We conjecture that this is due to our proof strategy, which relies on using a deviation bound on $||\hat{\Sigma}_{k,l_{n}}-\Sigma_{k}||_{2}$, and uses the functional dependence measure, rather than autocorrelation, to quantify dependence. In practice, we believe using GLSS, which corrects for serial correlation, and uses an estimator with lower asymptotic variance will achieve better performance. We illustrate this in more detail in our simulations section, and in the appendix (section 9.2).

Similar to SIS, we can control the size of the model selected by GLSS. For the case when $\beta_{k}^{M}=\rho_{k}\textrm{ }\forall k$, the bound on the selected model size is the same as in SIS. However, we need to place an additional assumption when $\beta_{k}^{M}\neq\rho_{k}$: If the cross covariance, $\gamma_{\bm{X}_{k},Y}(h)\propto E(X_{tk}Y_{t})$, for $h\in\{-L_{k},\ldots,-1,1,\ldots,L_{k}\}$, we can bound the selected model size by the model size selected by SIS. More formally we have:

Corollary 4**.**

Assume the cross covariance, $\gamma_{\bm{X}_{k},Y}(h)\propto E(X_{k,t}Y_{t})$, for $h\in\{-L_{k},\ldots,-1,1,\ldots,L_{k}\}$

1. (i)

Assume Conditions E,F,G,H hold, then for $\gamma_{n}=c_{5}n^{-\kappa}$ with $c_{5}\leq c_{6}/2$:

[TABLE] 2. (ii)

Assume Conditions E,F,G,I hold, then for $\gamma_{n}=c_{5}n^{-\kappa}$ with $c_{5}\leq c_{6}/2$:

[TABLE]

5 Second Stage Selection with Adaptive Lasso

The adaptive Lasso, as introduced by [58], is the solution to the following:

[TABLE]

and $\hat{\beta}_{I,j}$ is our initial estimate. For sign consistency; when $p_{n}>>n$, the initial estimates can be the marginal regression coefficients provided the design matrix satisfies the partial orthogonality condition as stated in [29], or we can use the Lasso as our initial estimator provided the restricted eigenvalue condition holds (see [39]). Both of these conditions can be stringent when $p_{n}>>n$. This makes the adaptive Lasso a very attractive option as a second stage variable selection method, after using screening to significantly reduce the dimension of the feature space. We have the following estimator:

[TABLE]

Where $\bm{X}_{\hat{\mathcal{M}}_{\gamma_{n}}}$ denotes the $n\times d_{n}$ submatrix of $\bm{X}$ that is obtained by extracting its columns corresponding to the indices in $\hat{\mathcal{M}}_{\gamma_{n}}$. We additionally define $\bm{X}_{\mathcal{M}_{\gamma_{n}}}$ accordingly. Our initial estimator $\hat{\bm{\beta}_{I}}=(\hat{\beta}_{I,1},\ldots,\hat{\beta}_{I,d_{n}})$ is obtained using the Lasso. Let $\hat{\Sigma}_{\mathcal{M}_{\gamma_{n}}}=\bm{X}_{\mathcal{M}_{\gamma_{n}}}^{T}\bm{X}_{\mathcal{M}_{\gamma_{n}}}/n$, and let $\Sigma_{\mathcal{M}_{\gamma_{n}}}$ be its population counterpart. Our two stage estimator, $\hat{\bm{\beta}}_{\hat{\mathcal{M}}_{\gamma_{n}}}$, is then formed by inserting zeroes corresponding to the covariates which were excluded in the screening step, and inserting the adaptive Lasso estimates, $\tilde{\bm{\beta}}_{\hat{\mathcal{M}}_{\gamma_{n}}}$, for covariates which were selected by the screening step. We need the following conditions for the combined two stage estimator to achieve sign consistency:

Condition J: The matrix $\Sigma_{\mathcal{M}_{\frac{\gamma_{n}}{2}}}$ satisfies the restricted eigenvalue condition,

RE($s_{n}$,3)(see [6] for details):

[TABLE]

where $\bm{v}=(v_{1},\ldots,v_{d_{n}^{\prime}})$ and $\bm{v}_{S}=(v_{i},i\in S),\bm{v}_{S^{c}}=(v_{i},i\in S^{c}).$

Condition K: Let $\lambda_{n}$ and $\lambda_{I,n}$ be the regularization parameters of the adaptive lasso and the initial lasso estimator respectively. For some $\psi\in(0,1)$, we assume:

[TABLE]

Condition L: Let $\beta_{min}=\min_{i\leq s_{n}}|\beta_{i}|$, and $w_{\max}=\max_{i\leq s_{n}}w_{i}>0$. Assume $\beta_{min}>\frac{2}{w_{\max}}$ and $\beta_{min}>2c\frac{\lambda_{I,n}s_{n}}{\phi_{0}n}$.

Condition J allows us to use the Lasso as our initial estimator. Notice that we placed the RE($s_{n}$,3) assumption on the matrix $\Sigma_{\mathcal{M}_{\frac{\gamma_{n}}{2}}}$, rather than the matrix $\hat{\Sigma}_{\hat{\mathcal{M}}_{\gamma_{n}}}$, given the indices in $\hat{M}_{\gamma_{n}}$ are random as a result of our screening procedure. Recall that for SIS, $\mathcal{M}_{\frac{\gamma_{n}}{2}}=\left\{1\leq i\leq p:|\rho_{i}|\geq\gamma_{n}/2\right\}$, and $|\mathcal{M}_{\frac{\gamma_{n}}{2}}|=d_{n}^{\prime}=O(d_{n})$, and for GLSS we have a similar definition. Therefore, we are placing the RE($s_{n}$,3) assumption on the population covariance matrix of a fixed set of $d_{n}^{\prime}$ predictors. Conditions K and L are standard assumptions, and are similar to the ones used in [39]. Condition K primarily places restrictions on the rate of increase of $\lambda_{n}$, and $\lambda_{I,n}$. Condition L places a lower bound on the magnitude of the non-zero parameters which decays with the sample size.

The next theorem deals with the two stage SIS-Adaptive Lasso estimator. A very similar result applies to the two stage GLSS-Adaptive Lasso estimator, if we replace Conditions A,B,C (resp. D) with Conditions E,F,G,H (resp. I), to avoid repetition we omit the result. For the following theorem, the terms $\iota,\omega,K_{x,r}\textrm{, and }K_{\epsilon,q}$ have been defined in the paragraph preceding Theorem 1, and $\tilde{\alpha}^{\prime},\tilde{\alpha}$ have been defined in Theorem 2.

Theorem 5**.**

1. (i)

Assume Conditions A,B,C,J,K,L hold, then for $\gamma_{n}=c_{3}n^{-\kappa}$ with $c_{3}\leq c_{1}/2$ we have:

[TABLE] 2. (ii)

Assume Conditions A,B,C,J,K,L hold, then for $\gamma_{n}=c_{3}n^{-\kappa}$ with $c_{3}\leq c_{1}/2$ we have:

[TABLE]

To achieve sign consistency for the case of finite polynomial moments we require:

Condition M: Assume $\lambda_{n}n^{\psi/2-1/2}\rightarrow\infty$ and $p_{n}=o(\min(\frac{s_{n}(n/s_{n})^{r/2-r\kappa/2}}{n^{\omega}},\frac{n^{\tau-\tau\kappa}}{n^{\iota}}))$, $d_{n}^{\prime}=o(min((n/s_{n})^{r/4-\omega/2},\lambda_{n}^{\tau}n^{\tau\psi/2-\iota}))$

For the case of exponential moments, we require:

Condition N: Assume $\lambda_{n}n^{\psi/2-1/2}\rightarrow\infty$,

$p_{n}=o(\min(\exp\left(\frac{Cn^{1/2-\kappa}}{s_{n}}\right)^{\tilde{\alpha}}/s_{n},\exp(Cn^{1/2-\kappa})^{\tilde{\alpha}^{\prime}}))$,

and $d_{n}^{\prime}=o(\min(\exp\left(\frac{n^{1/2}}{s_{n}}\right)^{\tilde{\alpha}/2},\exp\left(\lambda_{n}n^{\psi/2-1/2}\right)^{\tilde{\alpha}^{\prime}}))$

From Conditions M, N, and Theorem 3, we see an additional benefit of using the two stage selection procedure as opposed to using the adaptive Lasso as a stand alone procedure. For example, if we assume $d_{n}\leq n^{2\kappa}\lambda_{\max}(\Sigma)=O(n)$, and that both the error and covariate processes are sub-Gaussian, we obtain $p_{n}=o(\exp(n^{\frac{1-2\kappa}{3}}))$ for the two stage estimator. By setting $d_{n}^{\prime}=p_{n}$, we obtain the result when using the adaptive Lasso as a stand alone procedure, with the Lasso as its initial estimator. Under the scenario detailed above, the dimension of the feature space, which depends on $\lambda_{n}$ and $\psi$, for the stand alone adaptive Lasso can be at most $p_{n}=o(\exp(n^{\frac{1}{6}}))$. Therefore for $\kappa<1/4$, we obtain a larger range for $p_{n}$ and a faster rate of decay using the two stage estimator. For $\kappa\geq 1/4$ it is not clear whether the two stage estimator has a larger range for $p_{n}$, compared to using the adaptive Lasso alone.

The sign consistency of the stand alone adaptive Lasso estimator in the time series setting was established in [39]. Their result was obtained under strong mixing assumptions on the covariate and error processes, with the additional assumption that the error process is a martingale difference sequence. Additionally, in the ultrahigh dimensional setting they require a geometric decay rate on the strong mixing coefficients. In contrast, we obtain results for both the two stage and stand alone adaptive lasso estimator, and our results are obtained using the functional dependence measure framework. Besides assuming moment conditions, we are not placing any additional assumptions on the temporal decay of the covariate and error processes other than $\Delta_{0,q}(\bm{\epsilon})$,$\Phi_{0,q}(\bm{x})<\infty$. Furthermore, we weaken the martingale difference assumption they place on the error process, thereby allowing for serial correlation in the error process. Finally, by using Nagaev type inequalities introduced in [53], our results are easier to interpret and also allow us obtain a higher range for $p_{n}$.

6 Simulations

In this section, we evaluate the performance of SIS, GLSS, and the two stage selection procedure using the adaptive Lasso. For GLSS instead of using the banded estimate for $\Sigma_{k}$ we use a tapered estimate: $\hat{\Sigma}_{k}*R_{l_{n}}$, where $\hat{\Sigma}_{k}=(\hat{\gamma}_{i-j,k})_{1\leq i,j\leq n}$ and $R_{l_{n}}=(\max(1-\frac{|i-j|}{l_{n}},0))_{1\leq i,j\leq n}$ is the triangular kernel. We fix $l_{n}=15$, and we observed the results were fairly robust to the choice of $l_{n}$. In our simulated examples, we fix $n=200,s_{n}=6$ and $d_{n}=n-1$, while we vary $p_{n}$ from 1000 to 5000. We repeat each experiment 200 times. For screening procedures, we report the proportion of times the true model is contained in our selected model. For the two stage procedure using the adaptive Lasso, we report the proportion of times there was a $\lambda_{n}$ on the solution path which selected the true model.

Case 1: Uncorrelated Features

Consider the model (1), for the covariate process we have:

[TABLE]

Where $A_{1}=diag(\gamma)$, and we vary $\gamma$ from .4 to .6. We set $\bm{\eta}_{t}\sim N(0,\Sigma_{\eta})$, or $\bm{\eta}_{t}\sim t_{5}(0,V)$ in which case the covariance matrix is $\Sigma_{\eta}=(5/3)*V$. For this scenario we will be dealing with uncorrelated predictors, we set $\Sigma_{\eta}=I_{p_{n}}$. For the error process, we have an AR(1) process: $\epsilon_{i}=\alpha\epsilon_{i-1}+e_{i}$. We let $\alpha$ vary from .6 to .9, and let $e_{i}\sim t_{5}$ or $e_{i}\sim N(0,1)$. We set $\bm{\beta}=(\bm{\beta}_{1},\bm{\beta}_{2})$, where $\bm{\beta}_{1}=(.5,.5,.5,.5,.5,.5)$ and $\bm{\beta}_{2}=\bm{0}$. Even though the features are uncorrelated, this is still a challenging setting, given the low signal to noise ratio along with heavy tails and serial dependence being present.

The results are displayed in table 1. The entries below “Gaussian” correspond to the setting where both $e_{i}$ and $\bm{\eta}_{i}$ are drawn from a Gaussian distribution. Accordingly the entries under “$t_{5}$” correspond to the case where $e_{i}$ and $\bm{\eta}_{i}$ are drawn from a $t_{5}$ distribution. We see from the results that the performance of SIS, and GLSS are comparable when $p_{n}=1000$, with moderate levels of temporal dependence, along with Gaussian covariates and errors. Interestingly, in this same setting, switching to heavy tails seems to have a much larger effect on the performance of SIS vs GLSS. In all cases, the performance of GLSS appears to be robust to the effects of serial correlation in the covariate and the error processes. Whereas, for SIS the performance severely deteriorates as we increase the level of serial correlation. For example, for our highest levels of serial correlation, SIS nearly always fails to contain the true model.

Case 2: Correlated Features

We now compare the performance of SIS and GLSS for the case of correlated predictors. We have two scenarios:

Scenario A: The covariate process is generated from (25), with $A_{1}=diag(.4)$. $\bm{\eta}_{t}\sim N(0,\Sigma_{\eta})$, or $\bm{\eta}_{t}\sim t_{5}(0,V)$, with $\Sigma_{\eta}=\{.3^{|i-j|}\}_{i,j\leq p_{n}}$ for both cases. Therefore $\Sigma=\sum_{i=0}^{\infty}.4^{2i}\Sigma_{\eta}$. We set $\bm{\beta}_{1}=(1,-1,1,-1,1,-1)$ and $\bm{\beta}_{2}=\bm{0}$. We have an AR(1) process for the errors: $\epsilon_{i}=\alpha\epsilon_{i-1}+e_{i}$, we vary $\alpha$ from .4 to .8, and set $e_{i}\sim t_{5}$ or $e_{i}\sim N(0,1)$

Scenario B: The covariate process is generated from (25), with $A_{1}=\{.4^{|i-j|+1}\}_{i,j\leq p_{n}}$. And $\bm{\eta}_{t}\sim N(0,\Sigma_{\eta})$, or $\bm{\eta}_{t}\sim t_{5}(0,V)$, with $\Sigma_{\eta}=I_{p_{n}}$ for both cases. Therefore $\Sigma=\sum_{i=0}^{\infty}(A_{1}^{T})^{i}A_{1}^{i}$. We set $\bm{\beta}_{1}=(1,-1,1,-1,1,-1)$ and $\bm{\beta}_{2}=\bm{0}$. We have an AR(1) process for the errors: $\epsilon_{i}=\alpha\epsilon_{i-1}+e_{i}$, and we vary $\alpha$ from .4 to .8. The errors are generated in the same manner as in scenario A above.

The results are displayed in tables 2, and 3 respectively. In scenario A, we have a Toeplitz covariance matrix for the predictors, and moderate levels of serial dependence in the predictors. The trends are similar to the ones we observed in case 1. The performance of SIS is sensitive to the effects of increasing the serial correlation in the errors, with the effect of serial dependence being more pronounced as we encounter heavy tail distributions. In contrast, increasing the level of serial dependence has a negligible impact on the performance of GLSS. For scenario B, we observe similar trends as in scenario A.

Case 3: Two Stage Selection

We test the performance of the two stage GLSS-AdaLasso procedure. We also compare its performance with using the adaptive Lasso on its own. We use the Lasso as our initial estimator and select $\lambda_{I,n}$ using the modified BIC introduced in [47]. [26] extended the theory of the modified BIC to the case where $p>n$, $p=o(n^{a}),\textrm{ }a>1$, and independent observations. We conjecture that the same properties hold in a time series setting. We have two scenarios:

Scenario A: The covariate process is generated from (25), with $A_{1}=diag(.4)$. And $\bm{\eta}_{t}\sim N(0,\Sigma_{\eta})$, or $\bm{\eta}_{t}\sim t_{5}(0,V)$, with $(\Sigma_{\eta})_{i,j}=\{.8^{|i-j|}\}_{i,j\leq p_{n}}$. We set $\bm{\beta}_{1}=(.5,.5,.5,.5,.5,.5)$ and $\bm{\beta}_{2}=\bm{0}$. We have an AR(1) process for the errors: $\epsilon_{i}=\alpha\epsilon_{i-1}+e_{i}$, we vary $\alpha$ from .4 to .6, and set $e_{i}\sim t_{5}$ or $e_{i}\sim N(0,1)$

Scenario B: The covariate process is generated from (25), with $A_{1}=\{.4^{|i-j|+1}\}_{i,j\leq p_{n}}$. And $\bm{\eta}_{t}\sim N(0,\Sigma_{\eta})$, or $\bm{\eta}_{t}\sim t_{5}(0,V)$, with $(\Sigma_{\eta})_{i,j}=.8$ for $i\neq j$ and 1 otherwise. We set $\bm{\beta}_{1}=(.75,.75,.75,.75,.75,.75)$ and $\bm{\beta}_{2}=\bm{0}$. The errors are generated the same as in scenario A above.

In both scenarios we have a high degree of correlation between the predictors, low signal to noise ratio, along with mild to moderate levels of serial correlation in the covariate and error processes. The results are displayed in tables 4 and 5 for scenarios A and B respectively. We observe that the two stage estimator outperforms the standalone adaptive Lasso for both scenarios, with the difference being more pronounced in scenario B. For both scenarios, going from mild to moderate levels of serial correlation in the errors appears to significantly deteriorate the performance of the adaptive Lasso. This affects our results for the two stage estimator primarily at the second stage of selection. This sensitivity to serial correlation appears to increase as we encounter heavy tailed distributions.

7 Real Data Example: Forecasting Inflation Rate

In this section we focus on forecasting the 12 month ahead inflation rate. We use two major monthly price indexes as measures of inflation: the consumer price index (CPI), and the producer price index less finished goods (PPI). Specifically we are forecasting:

[TABLE]

Therefore the above quantities are approximately the percentage change in CPI or PPI over 12 months. Our data was obtained from the supplement to [32], and it consists of 132 monthly macroeconomic variables from January 1960 to December 2011, for a total of 624 observations. Apart from $\log(CPI)$ and $\log(PPI)$ which we are treating as $I(1)$, the remaining 130 macroeconomic time series have been transformed to achieve stationarity according to [32]. Treating $\log(CPI)$, and $\log(PPI)$ as $I(1)$, has been found to provide an adequate description of the data according to [44],[43],[39].

We consider forecasts from 8 different models. Similar to [39, 44] our benchmark model is an AR(4) model: $\hat{y}_{t+12}^{12}=\hat{\alpha}_{0}+\sum_{i=0}^{3}\hat{\alpha}_{i}y_{t-i}$ , where $y_{t}=1200\times log(CPI_{t}/CPI_{t-1})$ when forecasting CPI, and $y_{t}=1200\times log(PPI_{t}/PPI_{t-1})$ when forecasting PPI. For comparison, we also consider an AR(4) model augmented with 4 factors. Specifically we have:

[TABLE]

Where $\bm{\hat{F}}_{t}$ are four factors which are estimated by taking the first four principal components of the 131 predictors along with three of their lags. We also consider forecasts estimated by the Lasso and the adaptive Lasso. And lastly we include forecasts estimated by the following two stage procedures: GLSS-Lasso, GLSS-adaptive Lasso, SIS-Lasso, and SIS-Adaptive Lasso. Our forecasting equation for the penalized regression and two stage forecasts is:

[TABLE]

Where $\bm{x}_{t}$ consists of $y_{t}$ and three of its lags along with the other 131 predictors and three of their lags, additionally we also include the first four estimated factors $\hat{F}_{t}$. Therefore $\bm{x}_{t}$ consists of 532 covariates in total. For each of the two stage methods, we set $d_{n}=\lceil n/\log(n)\rceil=73$ for the first stage screening procedure. For the second stage selection, and the standalone lasso/adaptive lasso models, we select the tuning parameters and initial estimators using the approach described in section 6.

We utilize a rolling window scheme, where the first simulated out of sample forecast was for January 2000 (2000:1). To construct this forecast, we use the observations between 1960:6 to 1999:1 (the first five observations are used in forming lagged covariates and differencing) to estimate the factors, and the coefficients. Therefore for the models described above, $t$=1960:6 to 1998:1. We then use the regressor values at $t$=1999:1 to form our forecast for 2000:1. Then the next window uses observations from 1960:7 to 1999:2 to forecast 2000:2. Using this scheme, in total we have 144 out of sample forecasts, and for each window we use $n=451$ observations for each regression model. The set-up described above allows us to simulate real-time forecasting.

Table 6 shows the mean squared error (MSE), and the mean absolute error (MAE) of the resulting forecasts relative to the MSE and MAE of the baseline AR(4) forecasts. We observe that the two stage GLSS methods clearly outperform the benchmark AR(4) model, and appear to have the best forecasting performance overall for both CPI and PPI, with the difference being more substantial when comparing by MSE. Furthermore GLSS-lasso and GLSS-adaptive Lasso do noticeably better than their SIS based counterparts with the differences being greater when forecasting PPI. We also note that the widely used factor augmented autoregressions do worse than the benchmark model AR(4) model.

8 Discussion

In this paper we have analyzed the sure screening properties of SIS in the presence of dependence and heavy tails in the covariate and error processes. In addition, we have proposed a generalized least squares screening (GLSS) procedure, which utilizes the serial correlation present in the data when estimating our marginal effects. Lastly, we analyzed the theoretical properties of the two stage screening and adaptive Lasso estimator using the Lasso as our initial estimator. These results will allow practitioners to apply these techniques to many real world applications where the assumption of light tails and independent observations fails.

There are plenty of avenues for further research, for example extending the theory of model-free screening methods such as distance correlation, or robust measures of dependence such as rank correlation to the setting where we have heavy tails and dependent observations. Other possibilities include extending the theory in this work, or to develop new methodology for long range dependent processes, or certain classes of non-stationary processes. Long range dependence, is a property which is prominent in a number of fields such as physics, telecommunications, econometrics, and finance (see [40] and references therein). If we assume the error process ($\epsilon_{i}$) is long range dependent, then by the proof of Theorem 1 in [52] we have $\Delta_{0,q}(\bm{\epsilon})=\infty$. A similar result holds for the covariate process, therefore we may need to use a new dependence framework when dealing with long range dependent processes. Lastly, developing new methodology which aims to utilize the unique qualities of time series data such as serial dependence, and the presence of lagged covariates, would be a particularly fruitful area of future research.

9 Appendix

9.1 Proofs of Results

Proof of Theorem 1.

We first prove part (i), we start by obtaining a bound on:

[TABLE]

Let $T_{1}=\sum_{t=1}^{n}X_{tj}^{2}/n$, $T_{2}=\sum_{t=1}^{n}X_{tj}Y_{t}/n$. Then $|\hat{\rho}_{j}-\rho_{j}|=|T_{2}/T_{1}-E(T_{2})/E(T_{1})|=\\ |(T_{1}^{-1}-E(T_{1})^{-1})(T_{2}-E(T_{2}))+(T_{2}-E(T_{2}))/E(T_{1})+(T_{1}^{-1}-E(T_{1})^{-1})E(T_{2})|$

Therefore:

[TABLE]

For the RHS of (30), we obtain:

[TABLE]

Therefore it suffices to focus on terms (31), (32). For (31), recall that

Recall that $T_{2}=\sum_{t=1}^{n}X_{tj}(\bm{x}_{t}\bm{\beta}+\epsilon_{t})/n=\sum_{t=1}^{n}X_{tj}(\sum_{k=1}^{p_{n}}X_{tk}\beta_{k}+\epsilon_{t})/n$. Now we let:

[TABLE]

By Condition B, $E(X_{tj}\epsilon_{t})=0$, therefore

[TABLE]

Recall that $\sum_{k=1}^{p_{n}}\mathbbm{1}_{|\beta_{k}|>0}=s_{n}$, thus:

[TABLE]

From section 2 in [53]: $||X_{ij}||_{r}\leq\Delta_{0,r}(\bm{X}_{j})\leq\Phi_{0,r}(\bm{x})$. Using this we compute the cumulative functional dependence measure of $X_{tk}X_{tj}$ as:

[TABLE]

Therefore we obtain: $\sup_{m}(m+1)^{\alpha_{x}}\sum_{t=m}^{\infty}||X_{tj}X_{tk}-X_{tj}^{*}X_{tk}^{*}||_{r/2}\leq 2K_{x,r}^{2}$. Combining this with (36), and Theorem 2 in [53], yields:

[TABLE]

Similarly for $X_{tj}\epsilon_{t}$, by using Holder’s inequality we obtain:

[TABLE]

Therefore $\sup_{m}(m+1)^{\alpha}\sum_{t=m}^{\infty}||X_{tj}\epsilon_{t}-X_{tj}^{*}\epsilon_{t}^{*}||_{\tau}\leq 2K_{x,r}K_{\epsilon,q}$. Using Theorem 2 in [53], we obtain:

[TABLE]

For (32), assuming $E(X_{ij}^{2})=O(1)\textrm{ }\forall j\leq p_{n}$, and $\max_{j\leq p_{n}}E(X_{tj}Y_{t})<L<\infty$ we obtain:

[TABLE]

We set $M<\min_{j\leq p_{n}}E(X_{ij}^{2})-\epsilon$, for $\epsilon>0$. We then have:

[TABLE]

We can then bound the above two equations similar to (38). By combining (33)(35),(38),(40),(41), along with union bound we obtain:

[TABLE]

To prove part (ii), we follow the steps in the proof of Theorem 2 in [36]. Let $\mathcal{A}_{n}=\{\max_{k\in M_{*}}|\hat{\rho}_{k}-\rho_{k}|\leq\frac{c_{1}n^{-\kappa}}{2}\}$. On the set $\mathcal{A}_{n}$, by Condition A, we have:

[TABLE]

Hence by our choice of $\gamma_{n}$, we obtain $P\left(\mathcal{M}_{*}\subset\hat{\mathcal{M}}_{\gamma_{n}}\right)>P(\mathcal{A}_{n})$. By applying part (i), the result follows.

For part (iii) we follow the steps in the proof of Theorem 3 in [36]. Using $Var(Y_{t}),Var(X_{tj})=O(1)\textrm{ for }j\leq p_{n}$, along with Condition B, we obtain $\sum_{k=1}^{p_{n}}\rho_{k}^{2}=O(\lambda_{\max}(\Sigma))$. Then on the set $\mathcal{B}_{n}=\{\max_{k\leq p_{n}}|\hat{\rho}_{k}-\rho_{k}|\leq c_{4}n^{-\kappa}\}$, the number of $\{k:|\hat{\rho}_{k}|>2c_{4}n^{-\kappa}\}$ cannot exceed the number of $\{k:|\rho_{k}|>c_{4}n^{-\kappa}\}$ which is bounded by $O(n^{2\kappa}\lambda_{\max}(\Sigma))$. Therefore, by setting $c_{4}=c_{3}/2$ we obtain:

[TABLE]

The result then follows from part (i).

∎

Proof of Theorem 2.

We follow the steps from the proof of Theorem 1. Let $\bm{T}=(T_{1},\ldots,T_{n})$ where $T_{i}=X_{ij}X_{ik}$, and let $\bm{R}=(R_{1},\ldots,R_{n})$ where $R_{i}=X_{ij}\epsilon_{i}$. We need to bound the sums: $\sum_{i=1}^{n}(T_{i}-E(T_{i}))/n$ and $\sum_{i=1}^{n}R_{i}/n$.

By Theorem 1 in [49], $\Theta_{q}(\bm{T})\leq\Delta_{0,q}(\bm{T})$, and from Section 2 in [53]: $||X_{ij}||_{q}\leq\Delta_{0,q}(\bm{X}_{j})\leq\Phi_{0,q}(\bm{x})$. Additionally, by Holders inequality we have

[TABLE]

Using these, along with Condition D we obtain:

[TABLE]

Combining the above and using Theorem 3 in [53], we obtain:

[TABLE]

Similarly, using the same procedure we obtain:

[TABLE]

Now using the above bounds and following the steps in the proof of Theorem 1 we obtain the results.

∎

Proof of Lemma 1.

By the proof of Theorem 2 in [52], we have:

[TABLE]

Recall that $\hat{\rho}_{k}$ is the OLS estimate of the marginal projection, by (14) we have $\hat{\epsilon}_{t,k}=\epsilon_{t,k}-X_{tk}(\hat{\rho}_{k}-\rho_{k})=\epsilon_{t,k}-X_{tk}(\frac{\sum_{j=1}^{n}X_{jk}\epsilon_{j,k}/n}{\sum_{j=1}^{n}X_{jk}^{2}/n})$. Which gives us:

[TABLE]

By Condition E and $l_{n}=c\log(n)$, for sufficiently large $c$, we have: $\sum_{i=l_{n}+1}^{\infty}|\gamma_{i,k}|=o(n^{-\kappa})$, so we focus on the term $\sum_{i=1}^{l_{n}}|\hat{\gamma}_{i,k}-\gamma_{i,k}|$ in (49). We then have:

[TABLE]

And

[TABLE]

For (54), the bias $|E(\sum_{t=1}^{n-|i|}\frac{\epsilon_{t,k}\epsilon_{t+|i|,k}}{n}-\gamma_{i,k}|\leq\frac{i\gamma_{i,k}}{n}$. Using the techniques in the proof of Theorem 1 we can then bound (53). For (55) we have:

[TABLE]

And $P\left(\left|\frac{1}{n}\sum_{t=1}^{n-|i|}\epsilon_{t,k}X_{t+|i|,k}\right|>M\right)$

[TABLE]

And we set $M>\max_{k\leq p_{n}}\max_{i\leq l_{n}}2|E(\epsilon_{t,k}X_{t+|i|,k})|+\epsilon$, for some $\epsilon>0$. Similarly we have $P(\frac{|\sum_{j=1}^{n}X_{jk}\epsilon_{j,k}/n|}{\sum_{j=1}^{n}X_{jk}^{2}/n}>cn^{-\kappa}/Ml_{n})$

[TABLE]

Where we set $M_{1}<\min_{j\leq p_{n}}E(X_{ij}^{2})-\epsilon$, for $\epsilon>0$. The same method we used for (53) can be applied to (56), (57). Using the techniques in the proof of Theorem 1, and (52), we obtain the result. For (ii), we follow the same procedure as in (i), and apply the methods seen in the proof of Theorem 2.

∎

Proof of Theorem 3.

For (i), as before we start with a bound on: $P(|\hat{\beta}_{k}^{M}-\beta_{k}^{M}|>c_{2}n^{-\kappa})$. Using Condition E, we can write:

[TABLE]

After combining this with (14), it suffices to obtain a bound for:

[TABLE]

Similar to the proof of Theorem 1 we let $T_{1}=\bm{X}_{k}^{T}\hat{\Sigma}_{k,l_{n}}^{-1}\bm{X}_{k}/n$,

$T_{2}=\bm{X}_{k}^{T}\hat{\Sigma}_{k,l_{n}}^{-1}\bm{\epsilon}^{k}/n$, $T_{3}=E(\bm{X}_{k}^{T}\Sigma_{k}^{-1}\bm{X}_{k})$, and $T_{4}=E(\bm{X}_{k}^{T}\Sigma_{k}^{-1}\bm{\epsilon}^{k})$. Then:

[TABLE]

Following the steps in the proof of Theorem 1, it suffices to focus on the terms:

[TABLE]

We then have:

[TABLE]

We first deal with the term $\bm{X}_{k}^{T}\Sigma_{k}^{-1}\bm{\epsilon}^{k}/n$. We can rewrite this term as $\tilde{\bm{X}}_{k}^{T}\tilde{\bm{\epsilon}}^{k}/n$, where $\tilde{\bm{X}}_{k}=V_{k}\bm{X}_{k},\tilde{\bm{\epsilon}}^{k}=V_{k}\bm{\epsilon}^{k}/n$, $V_{k}$ is a lower triangle matrix and the square root of $\Sigma_{k}^{-1}$. Ignoring the first $L_{k}$ observations, we can express:

[TABLE]

, where $(\alpha_{1,k},\ldots,\alpha_{L_{k},k})$ are the autoregressive coefficients of the process $\epsilon_{t,k}$.

We compute the cumulative functional dependence measure of $\tilde{X}_{t,k}\tilde{\epsilon}_{t,k}$ as:

[TABLE]

We have: $||\tilde{X}_{l,k}-\tilde{X}_{l,k}^{*}||_{r}\leq||X_{l,k}-X_{l,k}^{*}||_{r}+\sum_{i=1}^{L_{k}}|\alpha_{i}|||X_{k,l-i}-X_{k,l-i}^{*}||_{r}$. And by our assumptions $||\tilde{\epsilon}_{l,k}-\tilde{\epsilon}_{l,k}^{*}||_{q^{\prime}}=0$, for $l>0$. From which we obtain:

[TABLE]

Using Theorem 2 in [53]:

[TABLE]

For the term $|\bm{X}_{k}^{T}(\hat{\Sigma}_{k,l_{n}}^{-1}-\Sigma_{k}^{-1})\bm{\epsilon}^{k}/n|$, using Cauchy-Schwarz inequality:

[TABLE]

Using (69) we obtain:

[TABLE]

Where the right hand side of (70) is:

[TABLE]

Let $M=M_{1}M_{2}$, where $M_{1}\geq\max_{k\leq p_{n}}E(X_{i,k}^{2})+\epsilon$, and $M_{2}=\max_{k\leq p_{n}}E(\epsilon_{i,k}^{2})+\epsilon$, for some $\epsilon>0$. The second term of (71) is:

[TABLE]

We can bound the above using the same techniques as in the previous proofs.

By Condition E, the spectral density of the process $\epsilon_{t,k},\forall k\leq p_{n}$ is bounded away from zero and infinity. Therefore, $0<C_{1}\leq\lambda_{min}(\Sigma_{k})\leq\lambda_{max}(\Sigma_{k})\leq C_{2}<\infty,\forall k\leq p_{n}$ [52]. We then use:

[TABLE]

Let $a_{1}\geq a_{2}\geq\ldots\geq a_{n}$ be the ordered eigenvalues of $\Sigma_{k}^{-\frac{1}{2}}\hat{\Sigma}_{k,l_{n}}\Sigma_{k}^{-\frac{1}{2}}$, therefore $||\Sigma_{k}^{\frac{1}{2}}\hat{\Sigma}_{k,l_{n}}^{-1}\Sigma_{k}^{\frac{1}{2}}-I_{n}||_{2}=\max_{i}|\frac{1}{a_{i}}-1|=\max_{i}|\frac{a_{i}-1}{a_{i}}|$. We then have

[TABLE]

Let $a_{j}=argmin_{a_{i}}|a_{i}^{-1}|$, using this and (73),(74) we obtain:

[TABLE]

Where $M_{3}\in(0,1-\epsilon)$ for $\epsilon>0$. We then have

[TABLE]

Combining the above with (75) and Lemma 1, we obtain:

[TABLE]

By (64),(68),(70),(72),(76) we obtain a bound for $P(|T_{2}-E(T_{2})|>Cn^{-\kappa})$. For the term $P(|T_{1}-E(T_{1})|>Cn^{-\kappa})$, we proceed in a similar fashion:

[TABLE]

We can then obtain a bound on the above terms by following a similar procedure as before. Combining these gives us the result for (i). For (ii), using the result from (i) we follow a similar procedure to the proof of Theorem 1. For (iii) and (iv) we follow the same procedure as (i) and (ii), and apply the methods seen in the proof of Theorem 2; we omit the details.

∎

Proof of Corollary 4.

Recall that:

[TABLE]

Therefore by our assumption, we have that $\beta_{k}^{M}\propto\rho_{k}$ whenever $\beta_{k}^{M}>0$. Using this we obtain $\sum_{k=1}^{p_{n}}(\beta_{k}^{M})^{2}=O(\sum_{k=1}^{p_{n}}\rho_{k}^{2})=O(\lambda_{\max}(\Sigma))$. We obtain the result, by following the procedure in the proof of Theorem 1 and using the results from Theorem 3.

∎

Proof of Theorem 5.

For simplicity we only prove part (i), the proof for part (ii) follows similarly. We will work on the following set $\mathcal{D}_{n}=\mathcal{A}_{n}\cap\mathcal{B}_{n}\cap\mathcal{C}_{n}$, where

[TABLE]

On the set $\mathcal{A}_{n}$, if we apply screening as a first stage procedure, by our choice of $\gamma_{n}$, we obtain:

[TABLE]

Next we need to use Lemma 7 and 8 in [39], specifically we need to show our reduced model satisfies conditions DGP 3,DESIGN, and WEIGHTS in [39]. On the set $\mathcal{B}_{n}$, by Lemma 1 in [39], we have $\phi_{\Sigma_{\hat{\mathcal{M}}_{\gamma_{n}/2}}}=\phi_{\Sigma_{\mathcal{M}_{\gamma_{n}/2}}}=\phi_{0}$. Therefore, we have:

[TABLE]

Using this along with Lemma 1 in [39] and Condition J, we have that DESIGN 3a is satisfied with $\phi_{\min}=\phi_{0}/16$, where $\inf_{v^{T}v=1}v^{T}\Sigma_{11}v>2\phi_{\min}>0$, and $\Sigma_{11}$ is the covariance matrix of the relevant predictors. On the set $\mathcal{D}_{n}$, by Conditions K and L in our work, and Lemma 2 and proposition 1 in [39], assumption WEIGHTS is satisfied. On the set $\mathcal{A}_{n}\cap\mathcal{B}_{n}$, DGP 3 and DESIGN 3b are satisfied, while DESIGN 2 is satisfied by Condition L.

Now by proposition 2, Lemmas 7 and 8 in [39] we obtain:

[TABLE]

$P(\mathcal{A}_{n}^{\complement})$ is given in Theorem 1 part i. For $P(\mathcal{B}_{n}^{\complement})$ using the method in the proof for Theorem 1, we obtain:

[TABLE]

And for $P(\mathcal{C}_{n}^{\complement})$:

[TABLE]

To prove part ii) we follow the same steps from part i). We obtain $P(\mathcal{A}_{n}^{\complement})$, $P(\mathcal{B}_{n}^{\complement}),P(\mathcal{C}_{n}^{\complement})$ by following the method in the proof of Theorem 2, and using Theorem 3 in [53]. ∎

9.2 Asymptotic Distribution of GLS estimator

Lemma 2**.**

Assume conditions E,F,G,H hold, then $\sqrt{n}(\hat{\beta}_{k}^{M}-\beta_{k}^{M})$ and $\sqrt{n}(\tilde{\beta}_{k}^{M}-\beta_{k}^{M})$ have the same asymptotic distribution.

Proof of Lemma 2.

It is clear that sufficient conditions for the feasible GLS estimator $\hat{\beta}_{k}^{M}$, and $\tilde{\beta}_{k}^{M}$ to have the same asymptotic distribution are [15]:

[TABLE]

By the proof of theorem 3, both these conditions are satisfied, therefore $\hat{\beta}_{k}^{M}$, and $\tilde{\beta}_{k}^{M}$ have the same asymptotic distribution. ∎

We use the above lemma, and rely on the asymptotic distribution of $\tilde{\beta}_{k}^{M}$ to provide an explanation for the superior performance of GLSS, and its robustness to increasing levels of serial correlation in $\epsilon_{t,k}$. We deal with three cases, and we assume an AR(1) process for the errors for simplicity and ease of presentation. The results can be generalized to AR$(p)$ processes, by using the moving average representation of $\epsilon_{t,k}$:

Case 1:

We start with the setting used in figure 1, assume $x_{t,k}$ is iid and $\epsilon_{t,k}=\alpha\epsilon_{t-1,k}+e_{t}$, with $x_{t,k}$, and $\epsilon_{t,k}$ being independent $\forall t$. Using Gordin’s central limit theorem [28], we calculate the asymptotic distribution of $\sqrt{n}(\tilde{\beta}_{k}^{M}-\beta_{k}^{M})\rightarrow N(0,J)$, where $J=\frac{\sigma_{e}^{2}}{\sigma_{x_{k}}^{2}(1+\alpha^{2})}$, $\sigma_{e}^{2}=var(e_{t})$, and $\sigma_{x_{k}}^{2}=var(x_{t,k})$ . Using the same methods we calculate the asymptotic distribution of the marginal OLS estimator as $\sqrt{n}(\hat{\rho}_{k}-\rho_{k})\rightarrow N(0,V)$, where $V=\frac{\sigma_{e}^{2}}{\sigma_{x_{k}}^{2}(1-\alpha^{2})}$. Therefore the variance of the OLS estimator increases without bound as $\alpha$ increases towards 1. Whereas the variance of the GLS estimator actually decreases as $\alpha$ increases.

Case 2:

We expand this to the case when $x_{t,k}$ is temporally dependent, for simplicity we let $x_{t,k}=\phi x_{t-1,k}+\eta_{t}$. We still assume $x_{j,k}$ and $\epsilon_{t}$ are independent $\forall j,t$, and $\epsilon_{t,k}=\alpha\epsilon_{t-1,k}+e_{t}$. This is the setting for the first model in the simulations section. Using Gordin’s central limit theorem, and elementary calculations: $\sqrt{n}(\tilde{\beta}_{k}^{M}-\beta_{k}^{M})\rightarrow N(0,J)$, where $J=\frac{(1-\phi^{2})\sigma_{e}^{2}}{(1+\alpha^{2}-2\phi\alpha)\sigma_{\eta}^{2}}$. And for the marginal OLS estimator $\sqrt{n}(\hat{\rho}_{k}-\rho_{k})\rightarrow N(0,V)$, where $V=\frac{(1+\phi^{2})\sigma_{e}^{2}}{(1-\alpha^{2})\sigma_{\eta}^{2}}$. We clearly see that for fixed $\phi$, the GLS estimate is robust to increasing $\alpha$, whereas the variance of the OLS estimator increases without bound as $\alpha$ increases towards 1. This sensitivity to $\alpha$ provides an explanation for the results seen in case 1 of the simulations, which show the performance of SIS severely deteriorates for high levels of serial correlation in $\epsilon_{t,k}$

Case 3:

In both the previous cases, it is easy to see the GLS estimator is asymptotically efficient to the OLS estimator. For the case where $\bm{X}_{k}=(x_{t,k},t=1,\ldots,n)$ and $\bm{\epsilon}^{k}=(\epsilon_{t,k},t=1,\ldots,n)$ are dependent on each other, it is more complicated. In this setting, it is likely the case that $\rho_{k}\neq\beta_{k}^{M}$. Assume $\epsilon_{t,k}=\alpha\epsilon_{t-1,k}+e_{t}$, and let $x_{t,k}-\alpha x_{t-1,k}=\tilde{x}_{t,k}$, and $W_{1}=\sum_{i=-\infty}^{\infty}\gamma(i)$, where $\gamma(i)=cov(\tilde{x}_{t,k}e_{t},\tilde{x}_{t-i,k}e_{t-i})$. We start by examining the asymptotic distribution of $\sqrt{n}(\tilde{\beta}_{k}^{M}-\beta_{k}^{M})\rightarrow N(0,J)$, where $J=W_{1}/(var(\tilde{x}_{t,k}))^{2}$. By the proof of theorem 1 in [52], $W_{1}\leq(\sum_{t=0}^{\infty}\delta_{2}(\tilde{x}_{t,k}e_{t}))^{2}$, which gives us:

[TABLE]

Where the last inequality follows from: $\delta_{2}(\tilde{x}_{t,k}e_{t})=||e_{0}||_{4}||\tilde{x}_{t,k}-\tilde{x}_{t,k}^{*}||_{4}+||\tilde{x}_{0,k}||_{4}||e_{t}-e_{t}{*}||_{4}$. Since $e_{t}$ is iid $||e_{t}-e_{t}{*}||_{4}=0,\forall t>0$. If we assume, $x_{t,k}=\phi x_{t-1,k}+\eta_{t}$, by writing $\tilde{x}_{t,k}=\eta_{t}+(\phi-\alpha)x_{t-1,k}$, we have:

[TABLE]

From these results we see that the asymptotic variance of the GLS estimator is bounded when $\alpha$ increases towards 1, and is largely robust to increasing levels of serial correlation in $\epsilon_{t,k}$. This result seems to provide an explanation for GLSS being robust to increasing levels of serial correlation in our simulations.

For the OLS estimator we obtain, $(\hat{\rho}_{k}-\rho_{k})\rightarrow N(0,V)$, where $V=W_{2}/(var(x_{t,k}))^{2}$ and $W_{2}=\sum_{i=-\infty}^{\infty}cov(x_{t,k}\epsilon_{t},x_{t-i}\epsilon_{t-i})$. As before, we can bound:

[TABLE]

We see the above bound is very sensitive to increasing serial correlation in $\epsilon_{t,k}$. Although this is an upper bound to the asymptotic variance, it seems to explain the deterioration in performance of SIS when increasing the serial correlation of $\epsilon_{t,k}$ in our simulations.

Acknowledgments

The author would like to thank the two anonymous referees, the Associate Editor, and the Editor for their helpful comments that have improved the paper.