On LASSO for Predictive Regression

On LASSO for Predictive Regression

\affil[a]{Department of Economics, University of Illinois. 1407 W. Gregory Dr., 214 David Kinley Hall, Urbana, IL 61801, USA.} \affil[b]{Department of Economics, the Chinese University of Hong Kong, Sha Tin, New Territories, Hong Kong SAR, China.} \affil[c]{Department of Economics, University of Southern California. Kaprielian Hall, 3620 South Vermont Avenue, Los Angeles, CA 90089, USA.}

Key words: cointegration, nonstationary time series, machine learning, shrinkage estimation, variable selection

JEL code: C22, C53, C61

\thispagestyle{empty}

{Introduction}

Predictive regressions are used extensively in empirical finance. A leading example is the stock return regression, for which predictability has long been a primary goal. We focus on this type of predictive regression in this paper. The first central econometric issue in these models is severe test size distortion in the presence of highly persistent predictors coupled with regression endogeneity. When persistence and endogeneity are high, the conventional inferential apparatus designed for stationary data can be misleading. Another major challenge in predictive regressions is the low signal-to-noise ratio (SNR). One explanation is that the noise in short-term predictive models may obscure the signal given by slow-moving (persistent) predictors.

Even though certain predictors have been suggested based on economic theory, there is no consensus regarding the predictive ability of any predictor. A researcher may use ex ante a pool of candidate regressors, hoping to identify a few important regressors. However, the more variables the researcher includes, the greater is the need for a data-driven routine for variable selection as many of these variables ex post demonstrate little or no predictive ability because of the competitive nature of the market. A low SNR naturally calls for variable selection.

Advancements in machine learning techniques, driven by an unprecedented abundance of data sources from many disciplines, offer valuable opportunities for economic data analysis. In the era of big data, shrinkage methods are becoming increasingly popular in econometric inference and prediction because of their variable selection and regularization properties. In particular, the least absolute shrinkage and selection operator tibshirani1996regression has received intensive study in the past two decades.

This paper investigates LASSO methods in predictive regressions. The LASSO methods are attractive in predictive regressions because they enable researchers to identify pertinent and exclude irrelevant predictors. However, time series predictors in predictive regressions carry heterogeneous degrees of persistence, which we call mixed roots. Some may exhibit short memory (e.g., Treasury bills), whereas others may be highly persistent (e.g., most financial/macro predictors). Moreover, a multitude of persistent predictors can be cointegrated. For example, the dividend price (DP) ratio is essentially a cointegrating residual between dividend and price, and the so-called cay data lettau2001consumption is another cointegrating residual between consumption, asset holdings, and labor income.

Given the difficulty of classifying time series predictors, we take an agnostic view of the types of the regressors and examine whether the LASSO methods can cope with heterogeneous regressors. We consider a predictive regression with three categories of predictors: short memory (I(0)) regressors, cointegrated local-to-unity persistent regressors, and non-cointegrated persistent regressors. Our focus is on variable selection. It is well known that in cross-sectional regressions, the adaptive LASSO zou2006adaptive enjoys variable selection consistency, which means that the estimated zeros exactly coincide with the true zero coefficients when the sample size is sufficiently large. The heterogeneous time series regressors challenge the conventional wisdom.

The main contribution of this paper is its study and refinement of Alasso in the time series context. As we consider a fixed number of regressors, we use OLS as Alasso's initial estimator. Alasso, with a proper choice of tuning parameter, consistently selects the non-cointegrating stationary and persistent variables, but it may over-select inactive cointegrating variables, where “inactive” means that these variables have zero regression coefficients, so they are excluded in the true data-generating process (DGP). The convergence rate of the OLS estimator is $\sqrt{n}$ for the cointegrated variables, and the resulting Alasso penalty weight is too small to eliminate the true inactive cointegrated groups because each cointegrating time series behaves as a local unit root process. Nevertheless, as is formally stated in Theorem (ref), with probability approaching one (wpa1), variables forming an inactive cointegrating group cannot all survive Alasso's selection.

This finding suggests a straightforward method of reclaiming the desirable oracle property: simply run a second-round Alasso on the variables selected by the first-round Alasso. As the first-round Alasso has broken the cointegration relationship between inactive cointegrating variables, the over-selected variables become non-cointegrated local-to-unity regressors in the second-round regression. In the post-selection OLS, the speed of convergence of these variables is boosted from the slow $\sqrt{n}$ rate to the fast $n$ rate. Then, the second-round Alasso can successfully suppress all of the inactive coefficients wpa1. We call this post-selection Alasso procedure twin adaptive LASSO and abbreviate it as TAlasso.

TAlasso achieves the oracle property fan2001variable, which implies that the rate of convergence is the same as the OLS when the relevant variables are known in advance and that the variable selection is consistent in the presence of a mixture of heterogeneous regressors. To the best of our knowledge, this paper is the first to establish these desirable properties in a nonstationary time series context. Note that the name TAlasso distinguishes it from the post-selection double LASSO belloni2014inference. Double LASSO is named for the double inclusion of selected variables in cross-sectional data to correct for LASSO's shrinkage bias in order to achieve uniform statistical inference. In contrast, our TAlasso features double exclusion in predictive regressions with mixed roots.

In developing asymptotic theory, we consider a linear process of time series innovations that encompasses a general ARMA structure arising in many practical applications, such as long-horizon return prediction. To focus on the distinctive feature of nonstationary time series, we adopt a simple asymptotic framework in which the number of regressors $p$ is fixed and the number of time periods $n$ passes to infinity. Our exploration in this paper is a stepping stone toward automated variable selection in high-dimensional predictive regressions.

We also explore the plain LASSO tibshirani1996regression, the standardized LASSO (Slasso, in which the $l_{1}$-penalty is multiplied by the sample standard deviation of each regressor) in the predictive regression with heterogeneous regressors. We call the phenomenon variable screening effect whereby an estimator coerces some coefficients to exactly zero fanlv2008jrssb. We find that neither Plasso nor Slasso maintains variable screening and consistent coefficient estimation in all of the components simultaneously, as Plasso imposes equal penalty weight regardless of the nature of the regressor, whereas Slasso's penalty only changes according to the scale variation of individual regressors and ignores their connection via cointegration.

In Monte Carlo simulations, we use various DGPs to examine the finite sample performance of Alasso/TAlasso compared with Plasso/Slasso, assessing their mean squared prediction errors and variable selection success rates. TAlasso is much more capable in selecting the correct model, thus improving prediction accuracy. These LASSO methods are further evaluated in a real data application. We use the 12 predictors from welch2008comprehensive to predict the S&P 500 returns. These 12 predictors exemplify the necessity of considering the three types of heterogeneous regressors. Please see Section (ref) for the details. Alasso and TAlasso attain better forecasting performance than their competitors and are robust in various estimation windows and prediction horizons.

\paragraph{Literature Review}

Since the original LASSO paper by tibshirani1996regression and the basis pursuit of \citet*{chen2001atomic}, a variety of important extensions of LASSO have been proposed, such as Alasso zou2006adaptive and the elastic net zou2005regularization. In econometrics, caner2009lasso, caner2014adaptive, shi2016estimation, \citet*{su2016identifying}, and kock2019inference use LASSO-type procedures in cross-sectional and panel data models. \citet*{belloni2011L1}, \citet*{belloni2012sparse}, \citet*{belloni2014inference}, and \citet*{Belloni2015uniforminf} develop methodologies and inferential theories in a variety of microeconometric settings.

ng2013variable surveys the variable selection methods in predictive regressions. Unlike the extensive literature on cross-sectional environments, the theoretical properties of shrinkage methods have received relatively little attention in time series models despite great empirical interest in macro and financial applications giannone2017economic,gu2018empirical,chinco2019sparse,rapach2019industry,freyberger2020dissecting. medeiros2016L1 study Alasso in high-dimensional stationary time series, and kock2015oracle discuss LASSO in a vector autoregression (VAR) system. In time series forecasting, inoue2008useful apply various model selection and model averaging methods to forecast U.S. consumer price inflation. sun2020time study model averaging methods under structural changes, and hirano2017forecasting develop a local asymptotic framework with independently and identically distributed (iid) orthonormalized predictors to study the risk properties of several machine learning estimators. There are even fewer papers on LASSO with nonstationary data. caner2013alternative discuss the bridge estimator, a generalization of LASSO, for the augmented Dickey-Fuller test in autoregressions, and kock2016consistent studies Alasso under the same setup. Earlier works on variable selection in the I(1) case utilize the different stochastic limit of penalty functions with stochastically trending regressors, such as phillips1994posterior and ng2001lag. For example, the stochastic limit of penalty factors is different with integrated regressors compared with the stationary case ng2001lag.

In predictive regressions, kostakis2014robust, lee2016predictive, and phillips2013predictive,phillips2016robust provide inferential procedures in the presence of multiple predictors with various degrees of persistence. xu2018testing studies variable selection and inference with possible cointegration between the I(1) predictors. In a vector error correction model (VECM), liao2015automated use Alasso for cointegration rank selection. Under a similar setting with high-dimensional I(0) and I(1) regressors and one cointegration group, koo2016high investigate Plasso's variable estimation consistency and non-standard asymptotic distribution. We differ from koo2016high in two aspects. (i) In an asymptotic framework that accommodates high-dimensional regressors, they focus on the rate of convergence of Plasso. In contrast, our focus is the oracle property in the mixed root model with three types of regressors. (ii) The unit root variables enter their model via the cointegrating groups only, so they do not deal with the balance of the predictive regression, whereas in our analysis, we allow local-to-unity persistent variables and local-to-zero small coefficients. Finally, smeekes2018macroeconomic demonstrate LASSO's numerical performance via simulations and empirical examples, while we are the first to systematically explore the theory concerning variable selection under mixed regressor persistence.

\paragraph{Notation}

We use standard notation. We define $\left\Vert \cdot\right\Vert _{1}$ and $\left\Vert \cdot\right\Vert $ as the usual vector $l_{1}$- and $l_{2}$-norms, respectively. The arrows $\Longrightarrow$ and $\overset{p}{\to}$ represent weak convergence and convergence in probability, respectively. $\asymp$ indicates “of the same asymptotic order, ” and $\sim$ signifies “being distributed as” either exactly or asymptotically, depending on the context. $\left\lfloor \cdot\right\rfloor $ is the floor function. $a \wedge b = \min\{a,b\}$. The symbols $O\left(1\right)$ and $o(1)$ ($O_{p}\left(1\right)$ and $o_{p}(1)$) denote (stochastically) asymptotically bounded or negligible quantities. For a generic set $M$, let $\vert M\vert$ be its cardinality. For a generic vector $\theta=(\theta_{j})_{j=1}^{p}$ with $p\geq\vert M\vert$, let $\theta_{M}=(\theta_{j})_{j\in M}$ be the subvector of $\theta$ associated with the index set $M$. $I_{p}$ is the $p\times p$ identity matrix, $\boldsymbol{1}_{p}$ is a column of $p$ ones, and $I\left(\cdot\right)$ is the indicator function.

The rest of this paper is organized as follows. Section (ref) introduces unit root regressors into a simple LASSO framework to fix ideas of Alasso. This model is substantially generalized in Section (ref) to include I(0), local-to-unity, and cointegrated regressors, and the asymptotic properties of Alasso and TAlasso are established. The mixed root model is carried over into Section (ref) to explore the behaviors of Plasso and Slasso. The theoretical results are confirmed through a set of empirically motivated simulation designs in Section (ref). Finally, we examine stock return predictions via these LASSO methods in Section (ref).

{Adaptive LASSO with Unit Roots}

In this section, we study LASSO with $p$ unit root regressors. To fix ideas, we investigate the asymptotic behavior of Alasso under a simple nonstationary regression model to better understand the technical issues arising from nonstationary predictors. Section (ref) generalizes the model to include the I(0), (local) unit roots, and cointegrated predictors.

Model

Assume that the dependent variable $y_{i}$ is generated from a linear model

where $n$ is the sample size. The $p\times1$ true coefficient is $\beta_{n}^{\ast}=(\beta_{jn}^{\ast}=\beta_{j}^{0\ast}/n^{\delta_{j}})_{j=1}^{p}$, where $\beta_{j}^{0\ast}\in\mathbb{R}$ is a fixed constant independent of the sample size and $\delta_{j}\in\left(0,1\right)$ is the local-to-zero rate, following phillips2013predictive and timmermann2017monitoring. If $\beta_{j}^{0\ast}=0$, then $\beta_{jn}^{\ast}$ remains zero regardless of the sample size; if $\beta_{j}^{0\ast}\neq0$, then it varies with $n$. This type of local-to-zero coefficient is designed to balance the I(0)--I(1) relation between the stock return and unit root predictors, and to model the weak SNR in predictive regressions phillips2015halbert. The $1\times p$ regressor vector $x_{i\cdot}=(x_{i1},\ldots,x_{ip})$ follows a pure unit root process

where $e_{k\cdot}=(e_{k1},...,e_{kp})$ is the innovation. For simplicity, we assume the initial value $e_{0\cdot}=0$. We use the following iid assumption on the innovations.

The regression equation ((ref)) can be equivalently written as

where $y=(y_{1},\ldots,y_{n})^{\prime}$ is the $n\times1$ response vector, $u=(u_{1},\ldots,u_{n})^{\prime}$, $x_{j}=(x_{1j},\ldots,x_{nj})^{\prime}$, and $X=(x_{1},\ldots,x_{p})$ is the $n\times p$ predictor matrix. This pure I(1) regressor model in ((ref)) is a direct extension of the common predictive regression application with a single unit root predictor (e.g., DP ratio). The mixed roots case in Section (ref) has multiple predictors and is more realistic in practice.

The literature focuses on the non-standard statistical inference caused by persistent regressors and weak signals. The asymptotic theory is usually confined to a small number of candidate predictors. Following the literature on predictive regressions, we consider the asymptotic framework in which $p$ is fixed and the sample size $n\rightarrow\infty$. This simple asymptotic framework allows us to concentrate on the contrast between the standard iid setting and the predictive regression involving nonstationary regressors.

In this model, the unknown coefficients $\beta_{n}^{\ast}$ can be obtained from the data by running OLS \[ \widehat{\beta}^{\mathrm{ols}}=\arg\min_{\beta}\Vert y-X\beta\Vert^{2}, \] whose asymptotic behavior is well understood phillips1987time. Assumption (ref) implies the following functional central limit theorem:

To represent the asymptotic distribution of the OLS estimator, define $u_{i}^{+}=u_{i}-\Sigma_{eu}^{\prime}\Sigma_{ee}^{-1}e_{i\cdot}^{\prime}$, and then $n^{-1/2}\sum_{i=1}^{\left\lfloor nr\right\rfloor }u_{i}^{+}\Longrightarrow B_{u^{+}}(r)$. By definition, $\mathrm{cov}\left(e_{ij},u_{i}^{+}\right)=0$ for all $j$ so that \[ \frac{X^{\prime}u}{n}\Longrightarrow\zeta:=\int_{0}^{1}B_{e}(r)dB_{u^{+}}(r)+\int_{0}^{1}B_{e}(r)\Sigma_{eu}^{\prime}\Sigma_{ee}^{-1}dB_{e}(r)^{\prime}, \] which is the sum of a (mixed) normal random vector and a non-standard random vector. The OLS limit distribution is

where $\Omega:=\int_{0}^{1}B_{e}(r)B_{e}(r)^{\prime}dr$. This result implies that when we inflate $\widehat{\beta}_{j}^{\mathrm{ols}}$ by the factor $n^{\delta_{j}}$ so that its magnitude is comparable to the constant $\beta_{j}^{0\ast}$, we attain consistency in that \[ n^{\delta_{j}}\left(\widehat{\beta_{j}}^{\mathrm{ols}}-\beta_{jn}^{\ast}\right)=n^{\delta_{j}}\widehat{\beta_{j}}^{\mathrm{ols}}-\beta_{j}^{0\ast}=O_{p} (n^{\delta_{j}-1})=o_{p}\left(1\right)\ \mbox{for all }j\leq p. \]

Even if $\beta_{j}^{0*}\neq0$, when $\delta_{j}$ is close to 1, the signal of $x_{j}$ is weak, and therefore the rate of convergence is slow. However, some true coefficients $\beta_{j}^{0\ast}$ in ((ref)) could be exactly zero, where the associated predictors would be redundant (inactive) in the regression. Let $M^{\ast}=\{j:\beta_{j}^{0\ast}\neq0\}$ be the index set of regressors relevant to the regression, $p^{\ast}=\left\vert M^{\ast}\right\vert $, and $M^{\ast c}=\left\{ 1,\ldots,p\right\} \backslash M^{\ast}$ be the set of redundant regressors. For simplicity, we refer to $M^{*}$ as the active set, meaning that it plays an active role in the regression, and we call $M^{*c}$ the inactive set. If we had knowledge about $M^{\ast}$, ideally, we would estimate the unknown parameters by OLS in the active set $M^{\ast}$ only. Define $\widehat{\beta}^{\mathrm{ora}}=\left(\widehat{\beta}_{M^{*}}^{\mathrm{ora}\prime},\widehat{\beta}_{M^{*c}}^{\mathrm{ora}\prime}\right)^{\prime}$, where \[ \widehat{\beta}_{M^{*}}^{\mathrm{ora}}=\arg\min_{\beta}\Vert y-\sum_{j\in M^{\ast}}x_{j}\beta_{j}\Vert^{2} \] and $\widehat{\beta}_{M^{*c}}^{\mathrm{ora}}=\boldsymbol{0}$. We call $\widehat{\beta}^{\mathrm{ora}}$ the oracle estimator, which is based on the infeasible “oracle” information of $M^{*}$. Eq.((ref)) implies that its asymptotic distribution is \[ n\left(\widehat{\beta}^{\mathrm{ora}}-\beta_{n}^{\ast}\right)_{M^{\ast}}\Longrightarrow\Omega_{M^{\ast}}^{-1}\zeta_{M^{\ast}}, \] where $\Omega_{M^{\ast}}$ is the $p^{\ast}\times p^{\ast}$ submatrix $\left(\Omega_{jj^{\prime}}\right)_{j,j^{\prime}\in M^{\ast}}$ and $\zeta_{M^{\ast}}$ is the $p^{\ast}\times1$ subvector $\left(\zeta_{j}\right)_{j\in M^{\ast}}$.

{Adaptive LASSO }

We study the asymptotic behavior of Alasso in predictive regressions with these pure unit root regressors. Alasso for ((ref)) is defined as

where the weight $\widehat{\tau}_{j}=|\widehat{\beta}_{j}^{\mathrm{init}}|^{-\gamma}$ for some initial estimator $\widehat{\beta}_{j}^{\mathrm{init}}$, and $\lambda_{n}$ and $\gamma$ are the two tuning parameters. In practice, $\gamma$ is often fixed at either 1 or 2, and $\lambda_{n}$ is selected as the primary tuning parameter. Thus, we discuss the case of a fixed $\gamma\geq1$ and the initial estimator $\widehat{\beta}^{\mathrm{init}}=\widehat{\beta}^{\mathrm{ols}}$ because of the fixed $p$ setting.

Alasso enjoys the oracle property in regressions with weakly dependent regressors medeiros2016L1. The following Theorem (ref) confirms that Alasso maintains the oracle property in regressions with unit root regressors. Let $\widehat{M}^{\mathrm{A}}=\{j:\widehat{\beta}_{j}^{\mathrm{A}}\neq0\}$ be Alasso's estimated active set, and $\bar{\delta}=\max\{(\delta_{j})_{j=1}^{p}\}$ be the fastest speed such that $\beta_{jn}^{*}$ shrinks to 0.

As shown in Theorem (ref) (a), the estimated active set $\widehat{M}^{\mathrm{A}}$ coincides with the true active set $M^{*}$ wpa1; in other words, $\widehat{\beta}_{j}^{\mathrm{A}}\neq0$ if $j\in M^{*}$ and $\widehat{\beta}_{j}^{\mathrm{A}}=0$ if $j\in M^{*c}$. (b) indicates that Alasso's asymptotic distribution in the true active set is as if the oracle $M^{*}$ is known. In this nonstationary regression, Alasso's adaptiveness is maintained through the proper choice of $\widehat{\tau}_{j}=|\widehat{\beta}_{j}^{\mathrm{ols}}|^{-\gamma}$. When the true coefficient is nonzero, $\widehat{\tau}_{j}$ delivers a penalty of a negligible order $\lambda_{n}/n^{1-\delta_{j}\gamma}\to0$, recovering the OLS limit theory. When the true coefficient is zero, $\widehat{\tau}_{j}$ imposes a heavier penalty of the order $\lambda_{n}/n^{1-\gamma}\rightarrow\infty$, thereby achieving consistent variable selection. The intuition in zou2006adaptive under deterministic design is generalized in our proof of the setting with nonstationary regressors.

Both $\bar{\delta}$ and $\gamma$ appear in the rate condition ((ref)) for generality. In practice, because $\bar{\delta}$ is an unknown feature of the DGP and $\gamma$ can be controlled by the user, we recommend setting $\gamma=1$, which simplifies ((ref)) as $ \lambda_{n} / n^{1-\bar{\delta}}+ 1/ \lambda_{n} \rightarrow0.$ A conservative choice of $\lambda_{n}$ diverging more slowly than polynomial orders of $n$, say, $\lambda_{n}\asymp\log\log n$, satisfies the above condition for all $\bar{\delta}>0$. The practical choice of the tuning parameter in the numerical simulations and the empirical application will be discussed in the corresponding sections.

{Adaptive LASSO with Mixed Roots }

In practice, we often encounter a multitude of candidate predictors that exhibit various dynamic patterns. Some are stationary, whereas others can be highly persistent and/or cointegrated. In this section, we discuss the theoretical properties of Alasso under a mixed persistence environment. We extend the model in Section (ref) to accommodate I(0) and (local) unit root regressors, with possible cointegration in the latter. The theory in this section provides general guidance for multivariate predictive regressions.

Model

We introduce three types of predictors into the model. A $1\times p_{c}$ cointegrated system $x_{i\cdot}^{c}=(x_{i1}^{c},...,x_{ip_{c}}^{c})$ has cointegration rank $p_{1}$, so that $p_{2}=p_{c}-p_{1}$ is the number of local unit roots in this cointegration system. Let $x_{i\cdot}^{c}=(x_{1i\cdot}^{c}, x_{2i\cdot}^{c})$ admit a triangular representation

where $R_{2}=I_{p_{2}}+c_{2}/n$ with $c_{2}=\mathrm{diag}\left( \check{c}_{1},...,\check{c}_{p_{2}}\right)$, $L$ is the lag operator, and the vector $v_{1i\cdot}$ is the cointegrating residual, whose initialization is an $O_{p}(1)$ stationary variable. Each local-to-unity parameter $\check{c}_{l}$ is finite, thus including both stationary and nonstationary local unit root regions. This local-to-unity specification includes the unit root process as a special case when $\check{c}_{l}=0$. The triangular representation Phillips1991cointinference,elliott1998 is a convenient and general form of cointegrated systems. xu2018testing also adopts this structure in predictive regressions.

Assume that $y_{i}$ is generated from the linear model

Each time series $z_{l}=(z_{1l},...,z_{nl})'$ is a stationary regressor. Each $x_{l}=(x_{1l},...,x_{nl})'$ is a local-to-unity process (initialized from $O_{p}(1)$ stationary process) such that $x_{i\cdot}=(x_{i1},...,x_{ip_{x}})$ satisfies $(I_{p_{x}}-R_{x}L)x_{i\cdot}^{\prime}=e_{i\cdot}^{\prime},$ where $R_{x}=I_{p_{x}}+c_{x}/n$ with $c_{x}=\mathrm{diag}\left(\tilde{c}_{1},...,\tilde{c}_{p_{x}}\right)$ for finite $\tilde{c}_{l}$.

Eq.((ref)) is infeasible because the cointegrating residual $v_{1i\cdot}$ is unobservable without a priori knowledge about the cointegration relationship. What we observe is the vector $x_{i\cdot}^{c}$ that contains cointegration groups. Substituting ((ref)) into ((ref)), we obtain a feasible regression equation

where $\phi_{2}^{\ast}=-A_{1}^{\prime}\phi_{1}^{\ast}$ and $\phi^{\ast}=\left(\phi_{1}^{\ast\prime},\phi_{2}^{\ast\prime}\right)^{\prime}$. Stacking the sample of $n$ observations, the infeasible and feasible regressions can be written as

where the variable $V_{1}$ is the $n\times p_{1}$ matrix that stacks $\left(v_{1i\cdot}\right)_{i=1}^{n}$, and $X_{1}^{c}$, $X_{2}^{c}$, $X^{c}$, $Z$, and $X$ are defined similarly. We explicitly define $\alpha^{\ast}=\alpha^{0\ast}$ and $\phi^{\ast}=\phi^{0\ast}$ as two coefficients independent of the sample size, which are associated with $Z$ and $X^{c}$, respectively. As in Section (ref), the coefficients associated with $X$ are specified as local-to-zero sequences $\beta^{\ast}=\beta_{n}^{\ast}=(\beta_{l}^{0\ast}/n^{\delta_{j}})_{l=1}^{p_{x}}$, where $\beta_{l}^{0\ast}$ is invariant to the sample size.

We assume a linear process for the innovation and cointegrating residual vectors. In contrast to the simplistic iid assumption in Section (ref), the linear process assumption is fairly general, including as special cases many practical dependent processes such as the stationary autoregressive and moving average processes. Let $v_{i\cdot}=\left(v_{1i\cdot},v_{2i\cdot}\right)$ and $p=p_{z}+p_{c}+p_{x}$.

{OLS}

As Alasso attaches a penalty term to the OLS criterion function, we first study the asymptotic distribution of the OLS estimator \[ \widehat{\theta}^{\mathrm{ols}}=(W^{\prime}W)^{-1}W^{\prime}y \] under the mixed roots, where $W = \left( Z, X^c, X \right)$ is the observed predictor matrix. To state the result, we define the true coefficients $\theta_{n}^{*}=\left(\alpha^{0\ast\prime},\phi_{1}^{0\ast\prime},\phi_{2}^{0\ast\prime},\beta_{n}^{\ast\prime}\right)^{\prime}$, where $\phi_{2}^{0*}=-A_{1}^{\prime}\phi_{1}^{0*}$, a diagonal normalizing matrix $R_{n}=\mathrm{diag}\left(\left(\sqrt{n}\boldsymbol{1}_{p_{z}+p_{1}}^{\prime},n\boldsymbol{1}_{p_{2}+p_{x}}^{\prime}\right)\right)$, and a rotation matrix $Q=

$.

Adaptive LASSO

Similar to Section (ref), Alasso for model ((ref)) is estimated as

where $\widehat{\tau}_{j}=|\widehat{\theta}_{j}^{\mathrm{ols}}|^{-\gamma}$. The literature on Alasso has established the oracle property in many models. caner2013alternative and kock2016consistent study Alasso's rate adaptiveness in a pure autoregressive setting with iid error processes. In their cases, the potential nonstationary regressor is the first-order lagged dependent variable, and the other regressors are stationary. Therefore, the components of different convergence rates are known in advance. We complement this line of nonstationary LASSO literature by allowing a general regression framework with mixed degrees of persistence. We also generalize the error processes to the commonly used dependent processes, which is important in practice. For example, the long-horizon return regression in Section (ref) requires this type of dependence in its error structure because of the overlapping return construction.

Surprisingly, Theorem (ref) shows that in the mixed root model, Alasso's oracle property holds partially, but not for all regressors. To discuss variable selection in this context, we introduce the following notations. We partition the index set of all regressors $\mathcal{M}=\left\{ 1,\ldots,p\right\} $ into four components: $\mathcal{I}_{0}$ (I(0) variables associated with $Z$), $\mathcal{C}_{1}$ (associated with $X_{1}^{c}$), $\mathcal{C}_{2}$ (associated with $X_{2}^{c}$), and $\mathcal{I}_{1}$ (nonstationary variables associated with $X$). Let $\mathcal{C}=\mathcal{C}_{1}\cup\mathcal{C}_{2}$ and $\mathcal{I}=\mathcal{I}_{0}\cup\mathcal{I}_{1}$. Let $M^{\ast}=\left\{ j:\theta_{j}^{0\ast}\neq0\right\} $ be the true active set for the feasible representation, and let $\widehat{M}^{\mathrm{A}}=\{j:\widehat{\theta}_{j}^{\mathrm{A}}\neq0\}$ be Alasso's estimated active set. Next, let $\mathcal{M}_{Q}=\mathcal{I}\cup\mathcal{C}_{1}$ be the set of coordinates that are invariant to the rotation in $Q$. Similarly, let $M_{Q}^{\ast}=M^{\ast}\cap\mathcal{M}_{Q}$ be the active set in the infeasible regression equation ((ref)), and let $M_{Q}^{\ast c}=\mathcal{M}\backslash M_{Q}^{\ast}$ be the corresponding inactive set. The DGP ((ref)) obviously implies $\mathcal{C}_{2}\cap M_Q^{\ast}=\emptyset$ and $\mathcal{C}_{2}\subseteq M_{Q}^{\ast c}$. For a generic index set $M\subseteq\mathcal{M}$, let $\mathrm{CoRk}\left(M\right)$ be the cointegration rank of the variables in $M$.

Were we informed of the oracle about the true active variables in $M^{*}$\ and the cointegration matrix $A_{1}$, we would transform the cointegrated variables into cointegrating residuals, discard the inactive variables, and then run OLS. Ideally, we would conduct an estimation with variables in $M_{Q}^{*}$ only. Such an oracle OLS shares the same asymptotic distribution as its Alasso counterpart in ((ref)). Outside of the active set $M_{Q}^{*}$, ((ref)) shows that all of the other (transformed) variables in $M_{Q}^{*c}$ consistently converge to zero.

The variable selection results in Theorem (ref) are novel and interesting. ((ref)) shows that variable selection is consistent for the pure I(0) and local unit root variables, which is in line with the well-known oracle property of Alasso. However, instead of confirming the oracle property, ((ref)) indicates that in the cointegration set $\mathcal{C}$, the selected $\widehat{M}^{\mathrm{A}}$ asymptotically contains the true active variables in $M^{\ast}$ but Alasso may over-select inactive variables. The inconsistency stems from the mismatch between the convergence rate of the initial estimator and the marginal behavior of a cointegrated variable viewed in isolation. Consider a pair of inactive cointegrating variables, such as $\left(x_{c}^{c},x_{h}^{c}\right)$ in Table (ref). The unknown cointegration relationship precludes transforming this pair into the cointegrating residual $v_{c}$. Without the rotation, OLS associated with the pair can only achieve the $\sqrt{n}$ rate according to ((ref)). The resulting penalty weights are insufficient to remove these variables that individually appear as nonstationary. In consequence, Alasso fails to eliminate both $x_{c}^{c}$ and $x_{h}^{c}$ wpa1. To the best of our knowledge, this is the first case of Alasso's variable selection inconsistency in an important empirical model.

Under condition ((ref)) all variables in the active cointegration groups have nonzero coefficients, and Alasso selects them asymptotically according to ((ref)). Despite potential variable over-selection in $\mathcal{C}$, ((ref)) brings relief: in the limit, the cointegration rank in Alasso's selected set, $\mathrm{CoRk}(\widehat{M}^{\mathrm{A}})$, must equal the active cointegration rank $\mathrm{CoRk}(M^{*})$. Note that under our agnostic perspective, we do not need to know or use testing procedures to determine the value of $\mathrm{CoRk}(M^{*})$. Let $C^{*c}=M^{*c}\cap\mathcal{C}$ be the index set of the inactive variables in $\mathcal{C}$. ((ref)) implies that the variables in $C^{*c}\cap\widehat{M}^{\mathrm{A}}$---Alasso's mistakenly selected inactive variables---cannot form cointegration groups. In mathematical expression, $P\left(\mathrm{CoRk}\left(C^{*c}\cap\widehat{M}^{\mathrm{A}}\right)=0\right)\to1.$

The intuition gleaned in these examples can be generalized to cointegration relationships involving more than two variables and multiple cointegrated groups. The proof of Theorem (ref) formalizes this argument by inspecting a linear combination of the corresponding Karush-Kuhn-Tucker (KKT) condition for the selected variables.

In the literature, Alasso achieved the oracle property in a single implementation, so there was no need to run it twice. In our predictive regression with mixed roots, it is possible that a single run of Alasso over-selects inactive variables in the cointegration system. Given that the cointegrating ties are all shattered wpa1 in ((ref)), further action can help fulfill the oracle property.

Twin Adaptive LASSO

When the sample size is sufficiently large, with high probability, the mistakenly selected inactive variables have no cointegration relationship, so they behave as isolated local-to-unity processes in the post-selection regression equation of $y$ on $\left(w_{j}\right)_{j\in\widehat{M}^{\mathrm{A}}}$, where $w_{j}=(w_{1j},\ldots,w_{nj})^{\prime}$ is the $j$-th time series regressors. This observation suggests the need to run a post-selection Alasso. We first obtain the post-Alasso initial OLS estimator \[ \widehat{\theta}^{\mathrm{po}}=\left(W_{\widehat{M}^{\mathrm{A}}}^{\prime}W_{\widehat{M}^{\mathrm{A}}}\right)^{-1}W_{\widehat{M}^{\mathrm{A}}}^{\prime}y. \] The post-selection OLS estimator $\widehat{\theta}_{j}^{\mathrm{po}}=O_{p}\left(n^{-1}\right)$ for the over-selected inactive cointegrated variables $j\in C^{*c}\cap\widehat{M}^{\mathrm{A}}$, rather than $O_{p}\left(n^{-1/2}\right)$, as for the first-around initial $\widehat{\theta}_{j}^{\mathrm{ols}}$. The resulting post-selection penalty level $\widehat{\tau}_{j}^{\mathrm{po}}=|\widehat{\theta}_{j}^{\mathrm{po}}|^{-\gamma}$ is sufficiently heavy to wipe out these redundant variables in another round of Alasso. We call this procedure TAlasso, for which $\widehat{\theta}^{\mathrm{TA}}=\left(\widehat{\theta}_{\widehat{M}^{A}}^{\mathrm{TA\prime}},\widehat{\theta}_{\mathcal{M}\backslash\widehat{M}^{A}}^{\mathrm{TA\prime}}\right)^{\prime}$, where $\widehat{\theta}_{\mathcal{M}\backslash\widehat{M}^{A}}^{\mathrm{TA}}=\boldsymbol{0}$ and

TAlasso asymptotically reclaims variable selection consistency for all types of variables. Let $\widehat{M}^{\mathrm{TA}}$ be the active set for TAlasso.

Faced with a variety of potential predictors with unknown orders of integration, we may not be able to sort them into different persistence categories in predictive regressions without potential testing errors smeekes2020unit. Our research provides a valuable practical algorithm to address this problem. TAlasso is the first estimator that achieves the desirable oracle property without requiring prior knowledge on the persistence of multivariate regressors. Despite the intensive study of LASSO in recent years, Theorems (ref) and (ref) are novel and important. With the cointegration system in the predictors, the former shows that Alasso does not automatically adapt to the behavior of a system of regressors. Nevertheless, it at least breaks all redundant cointegration groups so that its defects can be easily addressed by another round of Alasso. This solution echoes the repeated implementation of a machine learning procedure as in Phillips2019boosting.

Conventional LASSO with Mixed Roots

Originally, LASSO was proposed as a plain $l_{1}$-penalized regression without a sophisticated weighting scheme, motivated by the optimization problem's variable screening effect: some coefficients are estimated as exactly zeros in finite samples under a wide range of the tuning parameter tibshirani1996regression. As mentioned in the introduction, in this paper, we call this original estimator

Plasso. However, Plasso is scale-variant in the sense that if we change the unit of $w_{j}$ by multiplying it by a nonzero constant $c$, such a change is not reflected in the penalty term in ((ref)), so the Plasso estimator does not change proportionally to $\widehat{\theta}_{j}^{\mathrm{P}}/c$. To keep the estimation scale-invariant to the choice of the unit of $x_{j}$, researchers often scale-standardize LASSO as

where $\widehat{\sigma}_{j}=\sqrt{n^{-1}\sum_{i=1}^{n}\left(w_{ij}-\bar{w}_{j}\right)^{2}}$ is the sample standard deviation of $w_{j}$ and $\bar{w}_{j}$ is the sample average. In this paper, we call ((ref)) Slasso. Such standardization is the default option for LASSO in many statistical packages, such as the R package glmnet. Note that Alasso is also scale-invariant with $\widehat{\theta}^{\mathrm{init}}=\widehat{\theta}^{\mathrm{ols}}$ and $\gamma=1$.

Both Plasso and Slasso are special cases of $\min_{\theta}\left\{ \Vert y-X\theta\Vert_{2}^{2}+\lambda_{n}\sum_{j=1}^{p}\widehat{\tau}_{j}|\theta_{j}|\right\} $, where the former uses $\widehat{\tau}_{j}=1$ and the latter specifies $\widehat{\tau}_{j}=\widehat{\sigma}_{j}.$ When Slasso's scale standardization is carried out with stationary and weakly dependent regressors, each $\widehat{\sigma}_{j}^{2}$ converges in probability to its finite population variance. In the literature, the difference between the asymptotic properties of Plasso and Slasso are not of interest because Slasso uses another set of constant configurations for the penalty levels.

In sharp contrast, for those $w_{j}$ associated with $\mathcal{C}\cup\mathcal{I}_{1}$ that are individually nonstationary, by ((ref)), we have

where $d_{j}$ is a non-degenerate random variable (whose expression can be found in the proof of Corollary (ref)), so that the sample standard deviation $\widehat{\sigma}_{j}=O_{p}\left(\sqrt{n}\right)$ diverges in the limit. As a result, it imposes a much heavier penalty on the associated coefficients than the stationary time series does.

Adopting a standard argument for LASSO as in knight2000asymptotics, we derive the following results characterizing Plasso's asymptotic behavior with various choices of $\lambda_{n}$. For exposition, we define a function $D:\mathbb{R}^{3}\mapsto\mathbb{R}$ as $D\left(s,v,b\right)=s\left[v\cdot\mathrm{sgn}(b)I(b\neq0)+|v|I(b=0)\right]$, where $\mathrm{sgn}(b)=I(b>0) -I(b<0) $ is the sign function.

In Corollary (ref)(a), the tuning parameter is small and the limit distribution of Plasso is equivalent to that of OLS, so there is no variable screening effect. Similar to Plasso in the cross-sectional setting, there is no guarantee of consistent variable selection. The screening effect occurs when the tuning parameter $\lambda_{n}$ becomes larger. In view of the OLS rate of convergence in ((ref)), we can call those $\theta_{j}$ associated with $\mathcal{I}_{0}\cup\mathcal{C}$ the slow coefficients (at rate $\sqrt{n}$) and those associated with $\mathcal{I}_{1}$ the fast coefficients (at rate $n$). When $\lambda_{n}$ is raised to the magnitude in (b), the term $D(1,v_{j},\theta_{j}^{0*})$ causes variable screening among the slow coefficients but not the fast coefficients. If we further increase $\lambda_{n}$ to the level of (c), then the convergence rate of the slow coefficients is dragged down by the large penalty, but there is still no variable screening effect for the fast coefficients. The local-to-zero coefficient $\delta_{l}$ for those in $\mathcal{I}_{1}$ plays no role in these asymptotic results because $\widehat{\beta}_{l}$ converges to $\beta_{ln}^{*}$ at the fast rate $n$ regardless of $\delta_{l}$.

Corollary (ref) reveals a major drawback of Plasso in the mixed root model. Because it has one uniform penalty level for all variables, it is not adaptive to these various types of predictors. To induce variable screening in $\widehat{\theta}_{\mathcal{I}_{1}}^{\mathrm{P}}$, the tuning parameter must be ballooned to $\lambda/n\to c_{\lambda}\in(0,\infty]$, but the consistency of the slow coefficients would collapse under such a disproportionately heavy $\lambda_{n}$.

We now turn to Slasso. Unlike the behavior of the slow coefficients in Plasso, the penalty scheme of Slasso makes the cointegration components in $\widehat{\theta}^{\mathrm{S}}$ susceptible to variable screening at much smaller scales of the tuning parameter.

To summarize this section, neither Plasso nor Slasso induces consistent parameter estimation and variable screening effect simultaneously for all components of the mixed regressors. The weight $\widehat{\tau}_{j}$ is a constant for Plasso, while it exploits only the marginal variation of $w_{j}$ for Slasso. Plasso uses a single tuning parameter and does not adapt to the different orders of magnitude of the slow and fast coefficients. Slasso suffers from overwhelming penalties for coefficients associated with the cointegration groups. Because the cointegrated regressors are individually local unit root processes, they can only form a linear combination to produce a stationary time series when classified into a system. In contrast, the penalty of Alasso/TAlasso reflects the cointegration system via the initial/post-selection OLS.

Simulations

In this section, we examine via simulations the performance of the LASSO methods in forecasting and variable screening. We consider different sample sizes to demonstrate the quality of the asymptotic approximation in finite samples. Comparison is based on the one-period-ahead out-of-sample forecast.

Simulation Design

Following the settings in Sections (ref), (ref), and (ref), we consider four DGPs. Each DGP starts after a burn-in of $1000$ periods for the innovations to mimic stationary $O_{p}\left(1\right)$ realizations.

DGP 1 (Pure unit roots). This DGP corresponds to the pure unit root model in Section 2. Consider a linear model with 9 unit root predictors, the same number as the persistent regressions in Section (ref). Each regressor $x_{ij}$ is drawn from the random walk $x_{ij}=x_{i-1,j}+e_{ij}$, and the dependent variable $y_{i}$ is generated from $y_{i}=\gamma^{\ast}+x_{i\cdot}\beta_{n}^{\ast}+u_{i}$, where the innovations $\left(e_{i\cdot},u_{i}\right)\sim iid\;N\left(0,\Sigma\right)$. The intercept $\gamma^{\ast}=0.25$, the slope coefficient $\beta_{n}^{\ast}=\left(1,1,1,0,0,0,0,0,0\right)^{\prime}/\sqrt{n}$, and the covariance matrix $\Sigma$ is estimated from the data of welch2008comprehensive, as detailed in Supplement S1.

DGP 2 (Mixed roots and cointegration). This DGP is designed for the mixed root model in Section 3, which emulates the kitchen-sink approach in Section (ref). The dependent variable

where $\gamma^{\ast}=0.3$, $\alpha^{\ast}=\left(0.4,0,0\right)$, $\phi^{\ast}=\left(0.5,-0.5,0,0\right)$, and $\beta_{n}^{\ast}=\left(1/\sqrt{n},1/\sqrt{n},0,0,0\right)$. The vector $\xi_{i}=( z_{i\cdot}, v_{i\cdot}, e_{i\cdot},u_{i} )$ follows a VAR(1) process $\xi_{i}=\Phi\xi_{i-1}+\varepsilon_{i}$, where $\varepsilon_{i}\sim iid\;N\left(0,\Sigma_{\varepsilon}\right)$ in which $\Phi$ and $\Sigma_{\varepsilon}$ are empirically estimated using the data of welch2008comprehensive, as detailed in Supplement S1. $x_{i\cdot}^{c}\in\mathbb{R}^{4}$ is a vector I(1) process with cointegration rank $2$ based on the VECM, $\Delta x_{i\cdot}^{c}=\Gamma^{\prime}\Lambda x_{i-1,\cdot}^{c}+v_{i\cdot},$ where $\Lambda=

$ and $\Gamma=

$ are the cointegrating matrix and the loading matrix, respectively. $\left(x_{il}\right)_{l=1}^{5}$ are non-cointegrating random walks generated by $x_{il}=x_{i-1,l}+e_{il}$.

DGP 3 (Mixed roots with local-to-unity predictors) The coefficients and the innovation processes are the same as those in DGP 2, but we replace the I(1) processes in DGP 2 with local-to-unity processes. $x_{i\cdot}^{c}$ is a vector of a local-to-unity process with cointegration rank $2$ based on the triangular representation

$\left(x_{il}\right)_{l=1}^{5}$ are non-cointegrating local-to-unity processes generated by $x_{il}=\left(1- \check{c}_{l}/n\right)x_{i-1,l}+e_{il}$, $l=1,2,\cdots,5$, and $( \check{c}_{1},\check{c}_{2},\cdots,\check{c}_{5})=\left(7.4934,\,5.9152,\,8.9277,\,0,\,0\right)$, where these local-unity coefficients are also estimated following xu2018testing. With some local-to-unity parameters being zero, DGP 3 incorporates a mixture of I(0), I(1), and local-to-unity processes.

DGP 4 (Stationary autoregression). The autoregressive distributed lag (ARDL) model is a classical specification for time series regressions. In addition to including lags of $y_{i}$, it is common to accommodate lags of predictors in predictive regressions, such as in medeiros2016L1. The stationary dependent variable in the following equation is generated from the ARDL model \[ y_{i}=\gamma^{\ast}+\rho^{\ast}y{}_{i-1}+\sum_{l=1}^{4}\left(\phi_{l1}^{\ast}x_{il}^{c}+\phi_{l2}^{\ast}x_{i-1,l}^{c}\right)+\sum_{l=1}^{5}\left(\beta_{l1n}^{\ast}x_{il}+\beta_{l2n}^{\ast}x_{i-1,l}\right)+\sum_{l=1}^{3}\left(\alpha_{l1}^{\ast}z_{il}+\alpha_{l2}^{\ast}z_{i-1,l}\right)+u_{i}, \] where $\gamma^{\ast}=0.3$, $\rho^{\ast}=0.4$, $\phi_{1}^{\ast}=\left(0.75,-0.75,0,0\right)$, $\phi_{2}^{\ast}=\left(0,0,0,0\right)$, $\beta_{1n}^{\ast}=\left(1.5/\sqrt{n},1,0,1/\sqrt{n},0\right)$, $\beta_{2n}^{\ast}=\left(0,-1,0,0,0\right)$, $\alpha_{1}^{\ast}=\left(0.6,0.4\right)$, $\alpha_{2}^{\ast}=\left(0.8,0\right)$, and $\alpha_{3}^{\ast}=\left(0,0\right)$. The predictors $x_{i}^{c}$, $x_{i}$, and $z_{i}$ are generated in the same way as in DGP 3.

As we develop our theory with regressors of fixed dimension, OLS is a natural benchmark. Another benchmark is the oracle OLS estimator under infeasible information. The sample sizes in our exercise are $n=80$, $120,$ $200$, $400$, and $800$. We run $10,000$ replications for each sample size and each DGP.

Each shrinkage estimator relies on its tuning parameter $\lambda_{n}$, which is the appropriate rate multiplied by a constant $c_{\lambda}$. We use fivefold cross validation (CV), where the sample is temporally ordered and then partitioned into five consecutive blocks, to guide the choice of $c_{\lambda}$. To ensure that the tuning parameter changes according to the rate specified in the asymptotic theory, we set $n=100$ and run an exploratory simulation 100 times for each method that requires a tuning parameter. In each replication, we use the fivefold CV to obtain $c_{\lambda}^{\left(1\right)},\ldots,c_{\lambda}^{\left(100\right)}$. We then fix $c_{\lambda}=\mathrm{median}(c_{\lambda}^{\left(1\right)},\ldots,c_{\lambda}^{\left(100\right)})$ in the full-scale 10,000 replications. To specify the tuning parameters for other sample sizes, we multiply the constant $c_{\lambda}$ that we have calibrated from $n=100$ by the rates suggested by our asymptotic theory. We multiply $c_{\lambda}$ by $\sqrt{n}$ for Plasso and Slasso, and by $\sqrt{n}/(\log \log n)$ for Alasso and TAlasso.

Performance Comparison

Table (ref)(a) reports the out-of-sample prediction accuracy in terms of the mean prediction squared error (MPSE) $E\big[(y_{n}-\widehat{y}_{n})^{2}\big]$. By the simulation design, unpredictable variation arises from the variance of the idiosyncratic error $u_{i}$, which is 1.

Plasso and Slasso achieve variable screening and consistent estimation, as the predictors are unit root processes in DGP 1. They are slightly better in forecasting than Alasso/TAlasso when the sample size is small; however, as the sample size increases to $n=800$, TAlasso surpasses the conventional LASSO estimators, which suggests that variable selection is conducive to forecasting in a large sample. In DGPs 2--4 of mixed roots and cointegrated regressors, the settings are more difficult for the conventional LASSO to navigate. TAlasso is the best performer, followed by Alasso, which also outperforms the conventional LASSO methods by a non-trivial margin.

Table (ref)(b) summarizes the variable screening performance. Recall that $M^{\ast}=\{j:\theta_{j}^{\ast}\neq0\}$ is the set of relevant regressors and $\widehat{M}=\{j:\widehat{\theta}_{j}\neq0\}$ is the estimated active set. We define two success rates for variable screening:

Here, $SR_{1}$ is the percentage of correct selections from the active set, and $SR_{2}$ is the percentage of correct removals from the inactive set. We also report the overall success rate of classification into zero or nonzero coefficients, defined as

These expectations in $SR$, $SR_{1}$, and $SR_{2}$ are computed by the average across the simulation replications.

In terms of the overall selection measure $SR$, TAlasso is the most effective and Alasso takes the second place. As the sample size increases, TAlasso's success rates approach 100% in all DGPs, which supports variable selection consistency. While the difference in $SR_{1}$ among these methods becomes negligible when the sample size is large, the gain of TAlasso and Alasso stems largely from $SR_{2}$. The asymptotic theory suggests that $\lambda_{n}\asymp\sqrt{n}$ is too small for Plasso to eliminate $0$ coefficients corresponding to the nonstationary regressors. Plasso and Slasso achieve high $SR_{1}$ at the cost of low $SR_{2}$. In view of the results in Table (ref), the advantage of TAlasso's variable screening capability helps with forecast accuracy in the context of predictive regressions where many included regressors actually exhibit no predictive power.\footnote{In Table (ref)(b) Slasso has lower $SR_{2}$ than Plasso in a few instances. However, due to the presence of $\widehat{\tau}_{j}=\widehat{\sigma}_{j}=O_{p}\left(\sqrt{n}\right)$ in the penalty term, in asymptotics it imposes a heavier penalty on coefficients of I(1) regressors than Plasso does. The reason is that in our simulations we fix $c_{\lambda}^{\mathrm{P}}$ (for Plasso) and $c_{\lambda}^{\mathrm{S}}$ (for Slasso) by CV separately. CV selects the tuning parameter $c_{\lambda}$ favoring lower MPSE and adjusts $c_{\lambda}$ in finite sample. For example, in DGP 2, $c_{\lambda}^{\mathrm{P}}=0.6576$ is much larger than $c_{\lambda}^{\mathrm{S}}=0.1676$. If we fixed $c_{\lambda}^{\mathrm{P}}$ by CV and set $c_{\lambda}^{\mathrm{S}}=c_{\lambda}^{\mathrm{P}}$, Slasso would have a much higher $SR_{2}$.}

Finally, as shown in Table (ref), we check the variable screening of the inactive cointegration group $\left(x_{i3}^{c},x_{i4}^{c}\right)$ in DGPs 2--4, where the true coefficients $\phi_{3}^{0*}=\phi_{4}^{0*}=0$. Alasso is effective in preventing both redundant variables from remaining in the regression according to the third column, while in the second column the chance of omitting one variable is around 40-50%, indicating that it breaks the cointegration pair but cannot screen out both inactive regressors simultaneously. In contrast, TAlasso successfully identifies both redundant variables with significantly higher probabilities, as in the first column. Plasso and Slasso break down in variable selection consistency, as our theory suggests.

Empirical Application

Data

We apply the LASSO methods to welch2008comprehensive's dataset to predict stock returns. We focus on the improvement in terms of prediction error and variable screening.

The dataset of welch2008comprehensive is one of the most widely used in predictive regressions. \citet*{koo2016high} update this monthly data from January 1945 to December 2012, and we use the same time span. The dependent variable is excess return (ExReturn), defined as the difference between the continuously compounded return on the S&P 500 index and the three-month Treasury bill rate. The estimated AR(1) coefficient of the excess return is $0.149$, indicating weak persistence. The 12 financial and macroeconomic predictors are introduced and depicted in Figure (ref). Three variables, namely ltr, infl, and svar, oscillate around the mean, and nine variables are highly persistent with estimated AR(1) coefficients greater than $0.95$. The pairs (tms, dfr) and (\texttt{dp}, \texttt{dy}) are visibly moving in a synchronized pattern that suggests potential cointegration, and it is difficult to determine whether cointegration holds among (\texttt{dp}, \texttt{dy}, \texttt{ep}). \texttt{ep} fluctuates with \texttt{dy} before 2000, but the link dissolves after 2000, and the two series diverge in opposite directions during the Great Recession. The presence of stationary and persistent predictors fits the mixed roots environment studied in this paper, and our agnostic approach avoids decision errors resulting from hypothesis testing.

As recognized in the literature, the signal of persistent predictors may become stronger in long-horizon return predictions cochrane2009asset. In addition to the one-month-ahead short-horizon forecast, we construct the long-horizon excess return as the sum of the continuous compounded monthly excess return on the S&P 500 index \[ \text{LongReturn}_{i}=\sum_{k=i}^{i+12h-1}\text{ExReturn}_{k}, \] where $h$ is the length of the forecasting horizon, and $h=1$ indicates one year. We use $h=1/12$ (one month), 1/4 (three months), 1/2 (half a year), 1, 2, and 3 in our empirical exercises.

Performance

We forecast short- and long-horizon stock returns recursively with either a $10$-year or $15$-year rolling window. All 12 variables are made available in the predictive regression, which welch2008comprehensive call the kitchen sink model. The tuning parameters for the shrinkage estimators are determined by a 10-fold CV on MPSE with consecutive partitions in each estimation window. The CV method is a preferred choice for prediction purposes. As a robustness check, we also use the Bayesian information criterion (BIC) to decide the tuning parameters for Alasso and TAlasso.\footnote{BIC is a popular choice geared toward variable selection, but in our context, it is incompatible with the Plasso/Slasso, which cannot cope with variable screening and consistent estimation simultaneously.}

The forecast returns of the LASSO methods are shown in Figure (ref) along with the true realized return in gray for $h=1/12=0.083$ (short horizon), $1$ (median horizon), and $3$ (long horizon). When $h=0.083$, the realized excess return resembles white noise that is extremely difficult to forecast. When the horizon is extended to $h=3$, the dynamics of the long-run aggregated return become clearer. For most of the time, Alasso and TAlasso track the realized return closely.

Table (ref) quantifies the forecast error in terms of the out-of-sample RMPSE (root MPSE) and mean predicted absolute error (MPAE) $E\left[\left|y_{n}-\widehat{y}_{n}\right|\right]$. In addition to OLS, which includes all of the variables without any screening, we include random walk with drift (RWwD), i.e., the historical average of the excess returns, $\widehat{y}_{n}=(n-1)^{-1}\sum_{i=1}^{n-1}y_{i}$, as another benchmark that utilizes no information from regressors. The results show that OLS loses in the short horizon, whereas RWwD suffers in the long horizon, indicating the ineffectiveness of the all-in and all-out approaches. Variable screening is essential for balanced performance in this empirical example. Alasso/TAlasso forecasts are more precise than Plasso/Slasso in general. When the horizon is $h=2$ or $h=3$, TAlasso can achieve the smallest RMPSE, and Alasso is stronger than Plasso and Slasso by a substantial margin. The results are robust when the tuning parameters are chosen by either CV or BIC.

There is an exceptional case of $h=1$ with a 15-year rolling window. In this particular instance, Alasso fails to foresee the recovery trend after the financial crisis in 2008, and another round of Alasso further worsens TAlasso. As shown in Figure (ref), when $h=1$, Plasso's forecast coincides with the movement of the realization during the recovery period after 2008, whereas Alasso/TAlasso moves in the opposite direction. Whereas large deviation exacerbates RMPSE, under MPAE, the gap between Plasso and Alasso/TAlasso in this case is narrowed or even reversed. Under the 10-year rolling window, all of the methods encounter difficulty around the financial crisis, and the difference between TAlasso and Slasso is negligible. Thus, we view the unsatisfactory RMPSE of Alasso/TAlasso here as an adverse case under the specific rolling window.

In terms of prediction performance, it is known that eliminating irrelevant predictors is more important than including relevant predictors; see, for example, ploberger2003empirical. The inclusion of the irrelevant predictors could be detrimental in forecasting contexts, and including irrelevant nonstationary predictors in predictive regression can be harmful because stock returns are stationary. In this sense, Alasso and TAlasso provide more conservative variable selection in predictive regressions. In Figure (ref), the instance with $h=1/2$ and the 15-year rolling window is used to illustrate the estimated coefficients under CV. The shrinkage methods select different variables over the estimation windows, indicating the evolution of the predictive models across time. Alasso and TAlasso throw out more variables than Plasso or Slasso and hence deliver more parsimonious models. For example, they completely eliminate the variables ltr and infl.

To highlight the potentially cointegrated variables, Figure (ref)(a) presents the Alasso and TAlasso estimates of (dp, dy) and (dfr, tms) by CV, and in Figure (ref)(b), BIC decides the tuning parameters and then produces the estimates. The table ahead of the subfigures lists the fraction of the active Alasso estimates annihilated by TAlasso over the rolling windows. In more than a quarter of the cases when Alasso selects dy or tms, TAlasso further penalizes the corresponding coefficient to be $0$ in the second round. In general, BIC tends to freeze out more variables than CV, especially for \texttt{dy} and \texttt{tms}. Although the numbers vary, similar patterns are found in other combinations of forecast horizons and rolling window lengths.

Conclusion

We explore LASSO procedures in the presence of stationary, nonstationary, and cointegrated predictors. Although it does not enjoy the well-known oracle property, Alasso breaks the links between inactive cointegration groups, and its repeated implementation TAlasso recovers the oracle property because of the differentiated penalty on the zero and nonzero coefficients. TAlasso is adaptive to a system of multiple predictors with various degrees of persistence, unlike Plasso's uniform penalty or Slasso's penalty based only on the marginal variation of each predictor. Moreover, TAlasso saves the effort of sorting the predictors according to their degrees of persistence, so we can be agnostic regarding the time series properties of the predictors. The automatic penalty adjustment of TAlasso guarantees consistent model selection and the optimal rate of convergence. Such desirable properties may improve the out-of-sample prediction in practice under complex predictive environments with a mixture of regressors.

To focus on the mixed root setting, we adopt the simplest asymptotic framework with a fixed $p$ and $n\to\infty$ to demonstrate the clear contrast between OLS, Alasso, TAlasso, Plasso, and Slasso. This asymptotic framework is in line with the state of the art of the predictive regression studies in financial econometrics \citep*{kostakis2014robust,phillips2016robust,xu2018testing}. However, the large number of potential regressors available in the era of big data indicate the need for a theoretical extension to allow for an infinite number of regressors in the limit. As the restricted eigenvalue condition \citep*{bickel2009simultaneous} is unsuitable in our context, in which the nonstationary part of the Gram matrix does not degenerate, a few papers onatski2018alternative,onatski2020spurious,zhang2019identifying,bykhovskaya2020cointegration point to the spectral analysis of a large random matrix bai2010spectral as the foundation of asymptotic analysis. In future research, we will apply the new technical apparatus to deal with the minimal eigenvalue of the Gram matrix to generalize the insight gleaned from low-dimensional asymptotics to high-dimensional asymptotics.

Another line of related literature concerns uniformly valid inference and forecasting after LASSO model selection; see \citet*{belloni2018valid,belloni2018uniformly} or \citet*{hirano2017forecasting}, for example. These papers allow for model selection error by LASSO and provide valid inference or prediction by introducing local limit theory with only small departures from the true models. Combining these recent developments with our current LASSO theory is anther exciting direction for future research.

Acknowledgements

We would like to thank the Editor, Serena Ng, the Associate Editor, and two anonymous referees for their thoughtful comments, which substantially improved this paper. We thank Mehmet Caner, Zongwu Cai, Yoosoon Chang, Changjin Kim, Bonsoo Koo, Zhipeng Liao, Tassos Magdalinos, Joon Park, Hashem Pesaran, Peter Phillips, Kevin Song, Jing Tao, Keli Xu, Jun Yu, and the seminar participants at Kansas, Indiana, Purdue, UBC, UW, Duke, KAEA, IPDC, AMES, and IAAE conferences for helpful comments. We also thank Bonsoo Koo for sharing the data for the empirical application. Shi acknowledges the financial support of the Hong Kong Research Grants Council No. 24614817 and No. 14500118. All remaining errors are ours.