Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
112,110 characters · 15 sections · 78 citation commands
On LASSO for High Dimensional Predictive Regression
\thispagestyle{empty}
Key words: Cointegration, Forecast, Macroeconomics, Time series, Unit root
JEL code: C22, C53, C55
Ziwei Mei: [email removed]. Corresponding author: Zhentao Shi: [email removed]. Tel: (852) 3943 1432. Fax: (852) 2603 5805. Address: 928 Esther Lee Building, the Chinese University of Hong Kong, Shatin, New Territories, Hong Kong SAR, China. We thank Anna Bykhovskaya, Jinyuan Chang, Jianqing Fan, Yingying Li, Alexey Onatskiy, Whitney Newey, Liangjun Su, and Etienne Wijler for helpful comments.
\onehalfspacing
Machine learning is a rapidly evolving field that has significantly reshaped numerous academic disciplines. While statisticians often concentrate on scenarios where the sample comprises independently and identically distributed (i.i.d.) observations, econometricians pay special attention to settings where variables are gathered over time. Temporal dependence plays a crucial role in these data generating processes (DGP). In the context of parameter estimation, weakly dependent data, under certain technical conditions, bear resemblance to i.i.d. data as the time span increases. However, many established theoretical results under i.i.d. data are inapplicable to instances where time series exhibit high persistence.
Prediction is an important theme of empirical macroeconomics and finance. Although forecasting the stock market is notoriously challenging, recent advancements in machine learning offer some silver lining gu2020empirical; forecasting macroeconomic variables, on the other hand, is more feasible Stock2012,medeiros2021forecasting. Macroeconomic time series encompass a diverse range of dynamic patterns. GDP, industrial production index, exchange rates, and money supply all provide perspectives on the economy's current state and may offer insights into its future trajectory.
Linear predictive regression is a straightforward model. However, the presence of persistent regressors can pose a multitude of theoretical and practical challenges in its estimation and inference. In recent years, the advent of macroeconomic big data has sparked research interest in digesting numerous potential variables in macroeconomics ng2013variable. The least absolute shrinkage and selection operator (LASSO) tibshirani1996regression, an off-the-shelf machine learning method for linear regressions, is one of such tools. While LASSO and its variants have been extensively explored in statistics for i.i.d. data, only a handful of recent econometric papers have examined predictive regression in the context of nonstationary regressors with growing dimensions.
This paper serves as a stepping stone toward understanding LASSO in high dimensional predictive regressions with persistent variables. In particular, it considers the setting when a large number of unit root regressors are present, and in scenarios where the number of regressors ($p$) exceeds the sample size ($n$). We allow the innovation processes to be time dependent and non-Gaussian. Under these conditions, we introduce novel asymptotic rates for the deviation bound (DB) and the restricted eigenvalue (RE), which will be discussed in Sections (ref) and (ref). DB and RE are two pivotal conditions that govern LASSO's behavior. We establish convergence rates under the nonstationary time series setting, which differ from those of i.i.d. data buhlmann2011statistics and weakly dependent data kock2015oracle,medeiros2016,Mogliani2021.
To enhance the practical relevance of our theory, we expand the stylized regression model with all unit root regressors in two ways. Firstly, we examine not only tibshirani1996regression's original LASSO, which imposes the same penalty level to all coefficients, but also a LASSO variant that standardizes each regressor with its sample standard deviation (s.d.). We refer to the former as Plain LASSO (Plasso) and the latter as Standardized LASSO (Slasso), following lee2022lasso. While Plasso is more straightforward for theoretical analysis, Slasso is commonly the default in applications. The asymptotic theory for Slasso under i.i.d. data can be easily extended from Plasso, as each sample s.d. is expected to converge in probability to a positive constant. However, the sample s.d. of a unit root process, when divided by $\sqrt{n}$, converges in distribution to a non-degenerate stochastic integral, introducing additional randomness and altering the convergence rate. Secondly, we consider a mix of unit root, stationary, and cointegrated regressors. The researcher maintains an agnostic stance and includes all these variables in the regression without pre-testing to categorize them. We find that while Slasso maintains asymptotic guarantees, Plasso encounters multiple challenges. These enhancements enrich the theory and broaden the applicability.
We utilize the FRED-MD database to employ LASSO in forecasting the unemployment rate in the United States. Initially, we include all 121 variables from the database as predictors. Adhering to a standard practice in empirical studies, we transform each nonstationary time series into a stationary one, and then compare the forecast results with those based on the raw data without any transformation. Our findings reveal that Plasso underperforms in comparison to Slasso when given the same set of predictors. Moreover, the raw data are stronger than the stationarized data as they better match the persistence of the dependent variable. Subsequently, we experiment with a more comprehensive setting, incorporating four lags of each predictor along with lagged dependent variables and extracted factors. With a total of 504 regressors, Slasso further reduces the prediction error, suggesting that macroeconomic domain knowledge is beneficial in guiding initial specifications.
This paper adds to a burgeoning literature concerning many nonstationary time series. lee2022lasso explore variable selection of (adaptive) LASSO in the low dimensional setting where $p$ is fixed and highlights that some well-known LASSO properties for i.i.d. data collapse when faced with nonstationary data. Several papers consider the minimum eigenvalue or the RE of the Gram matrix of root unit processes as $p/\sqrt{n}\to0$, which we call the case of moderate dimension. koo2020high study a predictive regression with the unit root regression forming cointegration systems, and they leave the symbol of RE in the rate of convergence; they do not provide a lower bound for the RE. fan2023predictive work with quantile regressions in a similar setting with an assumed RE. An explicit rate that bounds the minimum eigenvalue is deduced in zhang2019identifying as a by-product of their exploration of cointegration systems, and smeekes2021automated use it to bound the RE in moderate dimensional predictive regressions.
When we were preparing this manuscript, wijler2022restricted independently derived the RE of high dimensional ($p>n$) unit root processes under the assumption that the innovations are i.i.d. Gaussian, as in kock2015oracle. Our paper differs from wijler2022restricted in the following aspects. First, wijler2022restricted uses non-asymptotic tail bounds based on sub-Gaussian distributions Wainwright2019high to obtain the rate of convergence of Plasso. Our approach, based on the non-asymptotic deviation inequalities for the maximum and minimum eigenvalues of Wishart random matrices Wainwright2019high offers sharper rates. Second, we go beyond i.i.d. Gaussian and accommodate sub-exponential and temporally dependent innovations by leveraging the technique of Koml\'{o}s-Major-Tusn\'{a}dy coupling Komlos1975 to achieve Gaussian approximation. Third, as in lee2022lasso our paper provides a comprehensive discussion of Plasso and Slasso in the setting with a mix of unit root, stationary time series, and cointegrated variables.
While high dimensional estimation counts on a well-behaved RE, testing problems often involve the maximum eigenvalue or a few large eigenvalues. Techniques of eigen-analysis for large random matrices are carried over into nonstationary time series by zhang2018clt for unit root tests, by onatski2018alternative and bykhovskaya2022asymptotics,bykhovskaya2022cointegration for cointegration tests, and by onatski2021spurious for principal component analysis and spurious regressions.
Besides a handful of papers mentioned above, machine learning grows fast in econometrics, covering i.i.d. data chernozhukov2017double,caner2018asymptotically, panel data su2016identifying,Shiforthcoming,shi2023forward, weakly dependent time series yousuf2021boosting,babii2022machine, and nonstationary time series phillips2021boosting,mei2022boosted,masini2022counterfactual, to name a few.
The rest of the paper is organized as follows. Section (ref) introduces LASSO and two variants in implementation, namely Plasso and Slasso. We put them into a unified framework, via a lemma that highlights the two key building blocks. In Section (ref), we first focus on the low-level assumptions for DB and RE. We then apply them to obtain the rates of convergence of Plasso and Slasso given pure unit root regressors, respectively. To better match practical circumstances, we further study LASSO given mixed regressors. Section (ref) carries out Monte Carlo simulations and the results corroborate the theoretical analysis. Section (ref) applies LASSO to predict the unemployment rate. Section (ref) concludes the paper. All technical proofs and additional simulation results are relegated to the Online Appendices.
Prior to formal presentation, we set up the notations. The set of natural numbers, integers, real numbers, and complex numbers are denoted as $\mathbb{N}$, $\mathbb{Z}$, $\mathbb{R}$, and $\mathbb{C}$, respectively. The integer set $\{1,2,\cdots,n\}$ is denoted as $[n]$ for some $n\in\mathbb{N}$. The integer floor function and ceiling function are denoted as $\left\lfloor \cdot\right\rfloor $ and $\left\lceil \cdot\right\rceil $, respectively. For an $n$-dimensional vector $x=(x_{t})_{t\in[n]}$, the $L_{2}$-norm is $\left\Vert x\right\Vert _{2}=\sqrt{\sum_{t=1}^{n}x_{t}^{2}}$, the $L_{1}$-norm is $\left\Vert x\right\Vert _{1}=\sum_{t=1}^{n}\left|x_{t}\right|$, and its sup-norm is $\|x\|_{\infty}=\sup_{t\in[n]}|x_{t}|$; we use “double dots” to denote the demeaned version $\ddot{x}=x-\bar{x}\cdot1_{n}$, where $\bar{x}=n^{-1}\sum_{t=1}^{n}x_{t}$ and $1_{n}$ is a vector of $n$ ones. Let $0_{n}$ be an $n\times1$ zero vector, and $I_{n}$ be the $n\times n$ identity matrix. For a generic index set $\mathcal{M}\subset[p]$ for some $p\in\mathbb{N}$, we use $\mathcal{M}^{c}=[p]\backslash\mathcal{M}$ to denote its complement, and $x_{\mathcal{M}}=\{x_{j}\}_{j\in\mathcal{M}}$ to denote the subvector of $x$ with coordinates located in $\mathcal{M}.$ For a generic matrix $B,$ let $B_{ij}$ be the $(i,j)$th element, and $B^{\top}$ be its transpose. Let $\|B\|_{\max}=\max_{i,j}|B_{ij}|$, and $\lambda_{\min}(B)$ and $\lambda_{\max}(B)$ be the minimum and maximum eigenvalues, respectively. Define $a\wedge b:=\min\left\{ a,b\right\} $, and $a\vee b:=\max\left\{ a,b\right\} $. An absolute constant is a positive, finite constant that is invariant with the sample size. The abbreviation “w.p.a.1” is short for “with probability approaching one”. “$a_{n}\stackrel{\mathrm{p}}{\preccurlyeq}b_{n}$” means that there is an absolute constant, say $c$, such that the event $\left\{ a_{n}\leq cb_{n}\right\} $ holds w.p.a.1. Symmetrically, “$a_{n}\stackrel{\mathrm{p}}{\succcurlyeq}b_{n}$” means “$b_{n}\stackrel{\mathrm{p}}{\preccurlyeq}a_{n}$”.
Let $W_{t}=(W_{jt})_{j\in[p]}$ be a $p$-vector of regressors. At time $n$, an econometrician is interested in using a linear combination $\alpha+W_{n}^{\top}\theta$ to predict a future outcome $y_{n+1}$. To learn the coefficients $\alpha$ and $\theta$, she collects historical data $Y=\left(y_{t}\right)_{t\in[n]}$ ($n\times1$ vector) and $W=\left(W_{0},W_{1},\ldots,W_{n-1}\right)^{\top}=(W_{t-1}^{\top})_{t\in[n]}$ ($n\times p$ matrix). When $p$ is close to $n$, or larger than $n$, LASSO tibshirani1996regression is one of the off-the-shelf estimation methods. It minimizes the sum of squared residuals plus an $L_{1}$ penalty
where the intercept $\alpha$ is not penalized. The superscript “$\mathrm{P}"$ of the estimator signifies Plasso. Prediction is made as $\widehat{y}_{n+1}^{{\rm P}}=\widehat{\alpha}^{{\rm P}}+W_{n}^{\top}\widehat{\theta}^{{\rm P}}$.
Plasso is not scale-invariant, meaning that if we multiply a non-zero constant $c_{j}$ to a regressor $W_{j,t-1}$, the corresponding LASSO estimate will not change proportionally to $\widehat{\theta}_{j}^{\mathrm{P}}/c_{j}$. Given that scale-invariance is a desirable property, a common practice --- like the default option of LASSO via glmnet::glmnet(x,y) in the R software --- scale-standardizes each regressor by its sample s.d. $\widehat{\sigma}_{j}=(n^{-1}\sum_{t=1}^{n}(W_{j,t-1}-\bar{W}_{j})^{2})^{1/2}$. Let $D={\rm diag}(\widehat{\sigma}_{1},\widehat{\sigma}_{2},\cdots,\widehat{\sigma}_{p})$ be the diagonal matrix that stores the sample s.d., and the Slasso estimator is
for which the prediction is made as $\widehat{y}_{n+1}^{{\rm S}}=\widehat{\alpha}^{{\rm S}}+W_{n}^{\top}\widehat{\theta}^{{\rm S}}$.
To analyze Plasso and Slasso under the same framework, we write
where $H$ is a positive definite placeholder: $H=I_{p}$ in Plasso, or $H=D$ in Slasso. In this paper, we will focus on the high dimensional component $\widehat{\theta}$, and the intercept is obviously $\widehat{\alpha}=\bar{Y}-\bar{W}\widehat{\theta}$ as it is unpenalized. Substitute $\widehat{\alpha}$ back to the criterion function in ((ref)), the $\theta$ component is numerical equivalent to
where $\ddot{Y}=Y-\bar{Y}1_{n}$ and similar demeaning applies to each column of $W$ to produce $\ddot{W}$.
The above is the numerical programming independent of the DGP. Now, suppose that the dependent variable is generated by
where $\left(\alpha^{*},\theta^{*}\right)$ are the true parameters. Sparsity means that most elements in $\theta^{*}$ are exactly zero. Let $\mathcal{S}=\{j\in[p]:\theta_{j}^{*}\neq0\}$ be the true active set, i.e., the location of the non-zero components, with its cardinality $s=|\mathcal{S}|.$
Well-known since bickel2009simultaneous, the two essential building blocks for the convergence of high dimensional LASSO are the DB (See the condition in Lemma (ref) below) and the RE. Let $\widehat{\Sigma}=\ddot{W}^{\top}\ddot{W}/n$ be the sample covariance matrix of all regressors. For some $L>0$, the restricted eigenvalue is defined, in our context, as
where $\mathcal{R}(L,s)=\{\delta\in\mathbb{R}^{p}\backslash\{0\}:\|\delta_{\mathcal{M}^{c}}\|_{1}\leq L\|\delta_{\mathcal{M}}\|_{1},\ \text{for all }|\mathcal{M}|\leq s\}.$ As our paper focuses on the rate of convergence, without loss of generality we follow buhlmann2011statistics and bickel2009simultaneous by setting $L=3$ as a convenient choice of the constant and use $\widehat{\kappa}_{H}=\kappa_{H}(\widehat{\Sigma},3,s)$ to simplify the notation. The following finite sample bounds hold for the generic LASSO estimator ((ref)).
The condition requires that the tuning parameter $\lambda$ in the LASSO estimation should be chosen above the deviation $\|n^{-1}\sum_{t=1}^{n}H^{-1}\ddot{W}_{t-1}u_{t}\|_{\infty}$, which is governed by the DGP of $W_{t}$ and $u_{t}$. This is the DB condition. The convergence rates of the LASSO estimator, signified by the right-hand side expressions in Lemma (ref), are determined by the sparsity index $s$, the RE $\widehat{\kappa}_{H}$, and the tuning parameter $\lambda$. We study in the next section the conditions under which we can establish desirable rates for the DB and RE, and then apply these two quantities to Plasso and Slasso for their rates of convergence.
This paper highlights unit root regressors. While $W_{j}$ is for a generic regressor with coefficient $\theta_{j}$, we denote the $j$th unit root regressor as $X_{j}=\left(X_{j0},\ldots,X_{j,t-1}\right)^{\top}$ and use $\beta_{j}$ as its coefficient. In this section we consider a DGP
where $X_{t}=\left(X_{1t},\ldots,X_{pt}\right)^{\top}$ is a vector of $p$ unit root processes $X_{t}=X_{t-1}+e_{t}$, and for simplicity let the initial value $\left\Vert X_{t=0}\right\Vert _{\infty}=O_{p}(1)$. We concatenate it with the error term $u_{t}$ into a $(p+1)$-vector $v_{t}=(e_{t}^{\top},u_{t})^{\top}$, and assume it is generated from
where $\ensuremath{\varepsilon_{t}=(\varepsilon_{jt})_{j\in[p+1]}}$ is a $(p+1)\times1$ random vector and $\Phi$ is a $(p+1)\times(p+1)$ deterministic matrix. For each $j\in[p+1]$, the shock
follows a linear process phillips1992asymptotics, which yields temporal dependence.
We will use low-level assumptions to build up the two high-level asymptotic properties DB and RE. Throughout this paper, we take the number of regressors $p=p\left(n\right)$ and the sparsity index $s=s\left(n\right)$ as deterministic functions of the sample size $n$. In formal asymptotic statements, we explicitly send $n\to\infty$ only, while it is understood that $p(n)\to\infty$ as $n\to\infty$ whereas $s\left(n\right)$ is either fixed or divergent.
We begin with the DB, which involves $n$ and $p$ only. As we allow high dimensionality in that $p>n$, the model ((ref)) must be regularized by assumptions. We first impose Assumption (ref) concerning the marginal distribution of the underlying shocks $\eta_{jt}$.
In Assumption (ref) we assume i.i.d. $\eta_{jt}$ with density $f$ over both $j$ and $t$, following zhang2019identifying and smeekes2021automated. It allows us to invoke concentration inequalities in the high-dimensional setting. Condition ((ref)) is known as the sub-exponential tail condition, which includes the familiar sub-Gaussian tail as a special case.
Assumption (ref) is concerning the coefficient $\psi_{jd}$ in the linear process ((ref)), which governs the temporal dependence of $\varepsilon_{jt}$. For any $z\in\mathbb{C}$, we denote the polynomial $\psi_{j}(z)=\psi_{j0}+\sum_{d=1}^{\infty}\text{\ensuremath{\psi_{jd}\cdot z^{d}}}=1+\sum_{d=1}^{\infty}\text{\ensuremath{\psi_{jd}\cdot z^{d}}}$, where without loss of generality we normalize $\psi_{j0}=1$.
Assumption (ref) is a sufficient condition for linear processes to satisfy the geometrically strong mixing ($\alpha$-mixing) condition (See Lemma (ref) in the Appendix). Finite-order strictly stationary ARMA processes are special cases of ((ref)) as they admit MA($\infty$) representations with exponentially decaying coefficients. The sub-exponential tail in Assumption (ref) and strong mixing in Assumption (ref) are common conditions in high dimensional time series regressions fan2011high,fan2013large,ding2021high.
The cross-sectional dependence across the regressors, encoded in $\Phi$, must be regularized as well. Assumption (ref) is concerning $\Phi$, which maps $\varepsilon_{t}$ into $v_{t}$ via ((ref)). Let $\Omega=\Phi\Phi^{\top}$.
Assumption (ref) (a) controls the magnitude of cross-sectional correlation. It rules out the unfavorable cases where innovations are very strongly correlated. Part (b), together with the condition ((ref)), guarantees the sub-exponential tail of $v_{jt}=\sum_{\ell=1}^{p+1}\Phi_{j\ell}\varepsilon_{jt}$ for all $j\in[p+1]$.
The above assumptions have been prepared for DB.
It is known that $n^{-1}\sum_{t=1}^{n}\ddot{X}_{j,t-1}u_{t}=O_{p}\left(1\right)$ as $n\to\infty$ for an individual unit root process $X_{j,t-1}$ and stationary error $u_{t}$ phillips1986understanding. Here to accommodate all $p$ unit root time series in a uniform matter, the DB grows at a mild speed $(\log p)^{1+\frac{1}{2r}}$, where the $r$ from Assumption (ref) governs the rate of diminishing temporal dependence. When $r$ is arbitrarily large, $\varepsilon_{t}$ will approach to temporal independence and the rate on the right-hand side of ((ref)) is reduced to $(\log p)^{1+\frac{1}{\infty}}=\log p$. The constant “4” on the left-hand side replicates the same constant required for $\lambda$ in Lemma (ref).
When $H$ in ((ref)) is an identity matrix, we study the RE $\widehat{\kappa}_{I}=\kappa_{I}(\widehat{\Sigma},3,s)$ for Plasso, associated with the sample Gram matrix $\widehat{\Sigma}=\ddot{X}^{\top}\ddot{X}/n$. In the i.i.d. case, it is easy to establish RE as the $(j,k)$th entry $\widehat{\Sigma}_{jk}$ converges in probability to the population covariance ${\rm cov}(X_{jt},X_{kt})$ for any fixed $j,k$, and then in high dimension we can apply concentration inequalities to construct a uniform bound for $\|\widehat{\Sigma}-\Sigma_{X}\|_{\max}=o_{p}(1)$, where $\Sigma_{X}:=\mathbb{E}(\widehat{\Sigma})$ has minimum eigenvalue bounded away from 0; See bickel2009simultaneous. This strategy does not carry over into non-stationary data. Recall that $X_{t}=X_{t-1}+e_{t}$ and define $\Sigma_{e}:=\mathbb{E}\left(e_{t}e_{t}^{\top}\right).$ After scaling by $1/n$, for each fixed pair $(j,k)$ the random variable \[ n^{-1}\widehat{\Sigma}_{jk}\stackrel{d}{\to}\mathcal{D}_{jk}=\int_{0}^{1}\mathcal{B}_{j}(r)\mathcal{B}_{k}(r)dr-\int_{0}^{1}\mathcal{B}_{j}(r)dr\int_{0}^{1}\mathcal{B}_{k}(r)dr \] where \textquotedblleft $\stackrel{d}{\to}$\textquotedblright denotes convergence in distribution, $\mathcal{B}_{j}$ and $\mathcal{B}_{k}$ are two Brownian motions, and the limiting distribution $\mathcal{D}_{jk}$ is a non-degenerate stable law. Since the diagonal elements $\mathcal{D}_{jj}$ has non-trivial probability in any small neighbor of zero, when the dimension $p$ accumulates $\widehat{\kappa}_{I}/n$ will shrink to 0. (\textcolor{black}{See Section (ref) for elaboration.}) This is in sharp contrast with the case of i.i.d. regressors, where the RE is bounded away from 0.
As the RE appears in the denominator of the error bounds in Lemma (ref), the convergence of LASSO requests that the RE shrinks to zero slowly enough. Lemma (ref) prepares an RE condition when the underlying innovations are i.i.d. normal, and Proposition (ref) allows non-Gaussian and time dependent innovations. This is one of the main theoretical contributions of this paper.
With a fixed $p$, phillips1990statistical show $\widehat{\Sigma}/n=\ddot{X}^{\top}\ddot{X}/n^{2}$ is positive-definite w.p.a.1\@. as $n\to\infty$. For unit root regressors the denominator under $\ddot{X}^{\top}\ddot{X}$ is $n^{2}$, instead of $n$ as in the i.i.d. case.\footnote{In high dimensional regressions, the rates of convergence under weakly dependent data are largely similar to those under the i.i.d. data. In the rest of the paper we mostly compare our results with what happens under the i.i.d. case for simplicity.} To align with this convention, we put $\widehat{\kappa}/n$ on the left-hand side of ((ref)). When $\varepsilon_{t}$ are i.i.d. normal, the right-hand side of ((ref)) gives a lower bound of RE proportional to $1/(s\log p)$. This result echoes smeekes2021automated where they establish $\lambda_{\min}(\widehat{\Sigma}/n)\stackrel{\mathrm{p}}{\succcurlyeq}1/p$ in the moderate dimensional case when $p/\sqrt{n}\to0$. Our ((ref)) replaces $\lambda_{\min}(\widehat{\Sigma}/n)$ by the restricted version $\widehat{\kappa}_{I}/n$ on the left-hand side, and replaces $p$ by $s\log p$ on the right-hand side.
The i.i.d. normality in Lemma (ref) is a strong assumption. To obtain an RE that accommodates more general innovations, we must control the relative magnitude among $n$, $p$ and $s$.
Assumption (ref) (a) allows $p$ to be of high dimension. The polynomial rate $n^{\nu}$ is for simplicity of presentation.\footnote{The proofs can still go through if we relax $p$ to grow at some exponential rate of $n$. But the speed of such rates will be peculiar to each of our convergence statement, thereby complicate the notations. For example, the right-hand side of ((ref)) would involve many specific terms.} In the meantime, by Assumption (ref) (b) the sparsity index cannot grow faster than $n^{1/4}$ when the innovations are non-normal.
Proposition (ref) substantially relaxes the distributional and dependence conditions by substituting the normality in Lemma (ref) with the sub-exponential tails in Assumptions (ref), and replacing i.i.d. with the mixing condition in Assumption (ref).
Given the two building blocks, DB and RE, we are ready to apply them to study the LASSO estimators.
Consider Plasso \[ \widehat{\beta}^{{\rm P}}=\text{\ensuremath{\arg\min_{\beta}\left\{ \dfrac{1}{n}\|\ddot{Y}-\ddot{X}\beta\|_{2}^{2}+\lambda\|\beta\|_{1}\right\} }} \] in the form of a special case of ((ref)). Parallel results to Lemma (ref) immediately follows.
It is well-known that with high dimensional i.i.d. data, Plasso's $L_{1}$ and $L_{2}$ estimation error bounds are $s\sqrt{\left(\log p\right)/n}$ and $\sqrt{s\left(\log p\right)/n}$, respectively, under standard conditions. Instead of the usual $\sqrt{n}$ for i.i.d. data, in the rates of convergence the denominators are $n$, yielding the familiar super-consistency when unit root regressors are present. The numerators, on the other hand, are multiplied by an extra factor $s(\log p)^{\frac{3}{2}+\frac{1}{2r}}$. This additional factor reflects the effect of the nonstationary time series, where $s\log p$ comes from the denominator of the lower bound of RE in ((ref)), and another $(\log p)^{\frac{1}{2}+\frac{1}{2r}}$ term stems from DB. When $r$ is arbitrarily large, the temporal dependence in $\varepsilon_{t}$ vanishes and the extra factor is reduced to $s\left(\log p\right)^{3/2}$. For example, in the special case of i.i.d. normal $\varepsilon_{t}$, in ((ref)) our convergence rate under the $L_{1}$-norm is $\dfrac{s^{2}}{n}(\log p)^{2}$ as $r=\infty$. This rate is faster than wijler2022restricted's Corollary 1, which shows \[ \|\ensuremath{\widehat{\beta}^{{\rm P}}-\beta}^{*}\|_{1}=O_{p}\left(\frac{s^{3}}{n^{1-\zeta_{1}}}(\log p)^{2}\right)=O_{p}\left(\frac{s^{2}}{n}(\log p)^{2}\times sn^{\zeta_{1}}\right) \] for any $\zeta_{1}>0$.
The tuning parameter $\lambda$ in Theorem (ref) involves an absolute constant $C_{\mathrm{DB}}$, which in turn depends on the absolute constants in the assumptions that are unknown in practice. Nevertheless, for all the left-hand side quantities in ((ref))--((ref)) to converge to zero in probability, it suffices if
which specifies a wide range of admissible rates for $\lambda$.
Plasso is the prototype of the $L_{1}$-penalized regression. In practice, Slasso is more often implemented in statistical software as scale-invariance is a desirable property. Again, we focus on the high dimensional coefficient $\beta$ in \[ \widehat{\beta}^{{\rm S}}:=\text{\ensuremath{\arg}\ensuremath{\min_{\beta}\left\{ \dfrac{1}{n}\|\ddot{Y}-\ddot{X}\beta\|_{2}^{2}+\lambda\|D\beta\|_{1}\right\} }}. \]
Proposition (ref) establishes the bounds for $\widehat{\sigma}_{\min}$ and $\widehat{\sigma}_{\max}$, which refresh the DB for Slasso as well as the RE $\widehat{\kappa}_{D}:=\kappa_{D}(\widehat{\text{\ensuremath{\Sigma}}},3,s)$, which is the restricted eigenvalue of the sample correlation coefficient matrix $D^{-1}\widehat{\Sigma}D^{-1}$ of the original data.
In terms of fitting, the leading term $s^{2}/n$ in ((ref)) is the same as that in ((ref)) for Plasso, up to an extra logarithm term. Super-consistency is preserved in $\widehat{\beta}^{{\rm S}}$ for the true original coefficient $\beta^{*}$. To counter the unknown absolute constant $\widetilde{C}_{{\rm DB}}$, it is sufficient to specify $\lambda$ as
for the consistency of the quantities in Theorem (ref).
In reality when we predict a target variable $y_{t}$ with many potential regressors, the regressors are most likely to have various dynamic patterns and we would not restrict ourselves by using nonstationary regressors exclusively. In the low dimensional case it is possible by pre-testing to classify variables into stationary and nonstationary ones, but the power of unit root tests are known to be weak in finite sample. What is worse, in high dimensional cases the individual test errors will accumulate in multiple testing procedures. For these reasons, lee2022lasso study LASSO with mixed-root regressors without pre-testing. That is, if we have a pool of mixed stationary and nonstationary regressors, we keep an agnostic view and throw them all into LASSO. It is in line with the attitude that a good machine learning method should adapt to the complex nature of the regression.
In this section we allow additional stationary regressors $Z_{t}$ in a DGP
The generic regressor $W_{t}$ and parameter $\theta^{*}$ in ((ref)) represent $(X_{t}^{\top},Z_{t}^{\top})^{\top}$ and $(\beta^{*\top},\gamma^{*\top})^{\top}$, respectively. Let $p_{x}$ be the length of $X_{t-1}$, and $p_{z}$ be the length of $Z_{j,t-1}$. We assume that the stationary components $Z_{t}$, $e_{t}$ and $u_{t}$ are potentially correlated in the following form
where we slightly abuse the notations to keep using $v_{t}$, $\varepsilon_{t}$ and $\Phi$, understanding that they are adapted to the mixed root case with the total number of regressors $p=p_{x}+p_{z}$.\footnote{The fact $E[\varepsilon_{t}]=0_{p+1}$ and ((ref)) imply that the stationary regressor $E\left[Z_{jt}\right]=0$. This restriction is merely for the conciseness of notation and there is no loss of generality. In the model ((ref)) the intercept $\alpha^{*}$ can absorb the non-zero means of the stationary regressors.} In addition to Assumptions (ref)--(ref) under the redefined symbols in this section, we impose one more condition.
Assumption (ref) is a necessary condition for identifying the coefficient $\gamma^{*}$ in ((ref)); otherwise $Z_{t-1}$ becomes endogenous and we must resort to external instrumental variables for identification and consistent estimation.
We study Slasso first in the mixed root case. We refresh $\widehat{\Sigma}=n^{-1}\ddot{W}^{\top}\ddot{W}$ and $\widehat{\kappa}_{D}=\kappa_{D}(\widehat{\Sigma},3,s)$ here. We have established in Proposition (ref) (a) that $\sqrt{n(\log p)^{-1}}\stackrel{\mathrm{p}}{\preccurlyeq}\widehat{\sigma}_{j}\stackrel{\mathrm{p}}{\preccurlyeq}\sqrt{n\log p}$ w.p.a.1. uniformly for all $j$ associated with the unit root variables. Under Assumptions (ref)--(ref) $\max_{j}\left|\widehat{\sigma}_{j}-\mathrm{s.d.}(Z_{j,t-1})\right|\stackrel{p}{\to}0$ for all $j$ associated the stationary variables, where $\mathrm{s.d.(\cdot)}$ denotes the population s.d. of a stationary regressor. The scale-standardization makes the DB and RE for the nonstationary variables comparable to those of the stationary variables up to some $\log p$ terms. As a result, convergence rates similar to Theorem (ref) follow.
For consistency, the admission rate for $\lambda$ in ((ref)) remains valid. Slasso provides provable rates of convergence for the mixed root case, thanks to the scale-standardization that aligns the unit root variables with the stationary ones. Up to some logarithmic term, the estimator $\widehat{\beta}^{{\rm S}}$ is super-consistent and $\widehat{\gamma}^{{\rm S}}$ maintains the standard $\sqrt{n}$ rate. Overall, $\widehat{\theta}^{\mathrm{S}}$ is consistent for $\theta^{*}$ under both the $L_{1}$ and $L_{2}$ norms.
Up to now we have addressed the asymptotic theory for a mix of $I(0)$ and $I(1)$ regressors. What happens if some nonstationary regressors are actually cointegrated? We discuss it in the next section.
To introduce the cointegrated variables into the predictive regression, we consider a cointegration system of $p_{c}$ observable variables with cointegration rank $p_{c1}$. We write the cointegration system into the triangular representation Phillips1991cointinference:
where $A^{(1)}$ ($p_{c1}\times p_{c2}$ matrix, where $p_{c2}:=p_{c}-p_{c1}$) stores the $p_{c1}$ cointegration vectors, and the cointegration error $v_{t}^{(1)}$ and the innovation $e_{t}^{(2)}$ are strictly stationary. Following lee2022lasso, we consider a model where the cointegration error $v_{t}^{(1)}$ enters the regression linearly via a coefficient $\phi_{1}^{*}$:
and $u_{t}$ is uncorrelated with $Z_{t-1}$ and $v_{t-1}^{(1)}$ to ensure the identification of their respective parameters $\gamma^{*}$ and $\phi_{1}^{*}$.
In practice the econometrician has no knowledge about the nature of the regressions a priori. Without the identities of $X_{t}^{\mathrm{co}(1)}$ and $X_{t}^{\mathrm{co}(2)}$ she cannot identify or estimate $A$, and thus $v_{t}^{(1)}$ is a latent variable, making ((ref)) an infeasible regression. She can, nevertheless, throw the $p$ ($=p_{x}+p_{z}+p_{c}$) observable regressors into the feasible regression
where $\phi_{2}^{*}=-A^{\top}\phi_{1}^{*}$ by substituting ((ref)) into ((ref)), $X_{t-1}^{\mathrm{co}}:=(X_{t-1}^{\mathrm{co}(1)\top},X_{t-1}^{\mathrm{co}(2)\top})^{\top}$ with associated parameter $\phi^{*}:=(\phi_{1}^{*\top},\phi_{2}^{*\top})^{\top}$, and $W_{t-1}$ collects all regressors with the corresponding parameter $\theta^{*}$ to fit into our framework ((ref)).
Penalized estimation methods face a generic challenge with cointegration systems in the regressors, as demonstrated in the following example.
We continue the example with LASSO's specific issue.
\setcounter{example}{0}
The above Example (ref) implies that in ((ref)) LASSO cannot achieve consistent estimation for the whole parameter $\theta^{*}$, which includes $\phi^{*}$ as a component. Instead, we should benchmark it with the a tailored regression of $y_{t}$ on $X_{t-1}$ and $Z_{t-1}$ only, where $\phi^{*}$ is suppressed to zero. In the DGP ((ref)) the latent variable $v_{t-1}^{(1)}$, if correlated with $Z_{t-1}$, will induce the well-known omitted variable bias in the population model
where $\omega^{*}:=\left[\mathbb{E}(Z_{t}Z_{t}^{\top})\right]^{-1}\mathbb{E}(Z_{t}v_{t}^{(1)\top})\phi_{1}^{*}$ is the projection coefficient of $v_{t-1}^{(1)\top}\phi_{1}^{*}$ onto the linear space spanned by $Z_{t-1}$, and the projection leads to the new regression coefficient $\gamma^{(1)*}=\gamma^{*}+\omega^{*}$ for $Z_{t-1}$ to ensure that the new residual $u_{t}^{(1)}=u_{t}+v_{t-1}^{(1)\top}\phi_{1}^{*}-Z_{t-1}^{\top}\omega^{*}$ is orthogonal to $Z_{t-1}$. From the perspective of prediction, this re-calibration of the population coefficient from $\gamma^{*}$ to $\gamma^{(1)*}$ is desirable in that the predictive power of $v_{t}^{(1)}$ can be partially absorbed by the observable $Z_{t-1}$ to reduce the variance of the error term as $\mathrm{var}\big(u_{t}^{(1)}\big)\leq\mathrm{var}\big(v_{t-1}^{(1)\top}\phi_{1}^{*}+u_{t}\big)$ by construction.
Now we present the formal asymptotic analysis of Slasso when cointegrated variables are present. We consider that the stationary components $v_{t}^{(1)},$ $e_{t}^{(2)}$, $e_{t}$, $Z_{t}$, and $u_{t}$ are potentially correlated in the form
understanding that $\Phi$ and $\varepsilon_{t}$ are redefined to adapt to the DGP ((ref)) and ((ref)). Define $\|A\|_{r1}:=\max_{j}\sum_{k}|A_{jk}|$ as the maximum row-wise $L_{1}$ norm. We further regularize the new coefficients $\gamma^{(1)*}$ and the cointegration matrix $A$. Let $\left\Vert \cdot\right\Vert _{0}$ be the cardinality of non-zero elements in a vector.
The restriction on $\|\beta^{*}\|_{0}+\|\gamma^{*(1)}\|_{0}$ controls the sparsity of the coefficients in ((ref)). The finite $L_{1}$-norm of $\|\omega^{*}\|_{1}+\|\phi_{1}^{*}\|_{1}$ governs the deviation bound for $u_{t}^{(1)}$. The restrictions on the cointegration matrix $A$ regularize the high dimensional cointegration system.
Given the discussion in Remark (ref) that Slasso is favored over Plasso when $X_{t-1}$ and $Z_{t-1}$ are present, we apply Slasso to ((ref)) and obtain the following results.
The in-sample fitting performance ((ref)) shows that Slasso for the feasible regression ((ref)) effectively learns the information in ((ref)). The parameter estimation performance in ((ref)) illustrates in terms of the $L_{1}$-norm that $\widehat{\beta}^{{\rm S}}$ is consistent for the pure unit root predictors $X_{t}$, whereas for the $I(0)$ regressors $\widehat{\gamma}^{{\rm S}}$ consistently estimates $\gamma^{*(1)}$ to absorb the information in $v_{t}^{(1)}$. As explained in Example (ref), Slasso shrinks $\widehat{\phi}^{{\rm S}}$ all the way to 0 due to the excessive penalty after scale-normalization. Such over-penalization violates the conditions in Lemma (ref) and we must devise a new technique to cope with the variable rotation in the proof (see Appendix (ref)). In summary, the Slasso estimator $\widehat{\theta}^{{\rm S}}$ converges in probability to the parameters in ((ref)), and the component associated with cointegrated variables is suppressed to 0 asymptotically.
In this section we carry out Monte Carlo simulations with mixed roots and pure unit roots.\footnote{We perform additional simulations based on DGP ((ref)) and ((ref)) with the presence of cointegrated variables. To save space, we defer the designs and the results to Section (ref) in the Appendix.} We first consider the DGP ((ref)) and generate the innovation $v_{t}=(e_{t}^{\top},Z_{t}^{\top},u_{t})^{\top}$ by a (vector) autoregressive (AR) process
where ${\bf 1}(\cdot)$ is the indicator function. The AR(1) coefficient is chosen to set the unconditional variance $0.84\Omega/(1-0.4^{2})=\Omega$, and the indicator function ensures Assumption (ref) with uncorrelated $Z_{j,t-1}$ and $u_{t}$. We consider $n\in\{120,240,360\}$ and $p=2n$. We try $p_{x}=\{0.5n,0.8n,1.2n,1.5n\}$ for each $n$, and thus $p_{z}=2n-p_{x}$ for each pair of $(n,p_{x})$$.$ The sparsity indices are $s_{x}=s_{z}=2\lceil\log n\rceil$, so that $s=4\lceil\log n\rceil$. We set the true coefficients of the stationary component $\gamma^{*}=(0.3\times[s_{z}]^{\top},0_{p_{z}-s_{z}}^{\top})^{\top}$. We specify two cases for the unit root regressors, which vary only in the coefficients $\beta_{(1)}^{*}=(n^{-1/2}1_{s_{x}}^{\top},0_{p_{x}-s_{x}}^{\top})^{\top}$ where the factor $n^{-1/2}$ bound the dependent variable $y_{t}$ to be of non-explosive, and $\beta_{(2)}^{*}=(1,n^{-1/2}1_{s_{x}-1}^{\top},0_{p_{x}-s_{x}}^{\top})^{\top}$where the coefficient of the first regressor is invariant with the sample size $n$. We label the data generated by the following coefficients as
For each DGP, we report the one-period-ahead out-of-sample root mean squared prediction error (RMSPE), defined as $\left\{ \mathbb{E}\left[(y_{n+1}-\widehat{y}_{n+1})^{2}\right]\right\} ^{1/2}$, and the parameter estimation root mean squared error, defined as $\left[\mathbb{E}\left(\Vert\widehat{\theta}-\theta^{*}\Vert{}_{2}^{2}\right)\right]^{1/2}$. The estimation method is either Plasso or Slasso. The expectations are approximated by the empirical average over 5000 replications. As a benchmark, we compare LASSO with an oracle estimator---the OLS with known active variables.
A key ingredient in implementing LASSO is the choice of the tuning parameter $\lambda$. One common data-driven approach is cross validation (CV). In our time series context, we cut $t\in[n]$ into 10 chronically ordered blocks and choose the $\lambda$ that minimizes the CV means squared error as each block serves as a validation dataset in turn whereas the other 9 blocks work as the training data. We refer to the $\lambda$ chosen by this time series 10-fold CV as “CV $\lambda$”. CV $\lambda$ is completely data-driven.
Alternatively, to evaluate our theoretical statement where $\lambda$ is specified as a constant multiplied by an expansion rate determined by $n$, $p$ and $r$, we follow lee2022lasso to use a small-scale experiment to calibrate an initial choice. We try 100 replications with $(n_{0},p_{x0})$, the smallest $n$ and $p_{x}$ considered in the simulations, save $\lambda$ in each replication according to the 10-fold CV described in the previous paragraph, and let $\widehat{\lambda}_{0}$ be the median of these $\lambda$'s. We then scale up $\widehat{\lambda}_{0}$ based on the theoretical expansion rate. We refer to this scheme as “calibrated $\lambda$”. For Slasso with an initial $\widehat{\lambda}_{0}^{(\mathrm{s})}$ obtained from $(n_{0},p_{x0},p_{z0})=(120,60,180)$, we use
to adhere to the rate $n^{-1/2}(\log p)^{\frac{3}{2}+\frac{1}{2r}}$ in Theorem (ref) when $p$ is proportional to $n$, where $r=1$ for the AR(1) innovation. Section (ref) has elaborated that Slasso enjoys theoretical guarantees whereas the convergence of Plasso with mixed regressors is unknown. As a numerical exercise we naively borrow the rate $(\log p)^{1+\frac{1}{2r}}$ in Theorems (ref) after obtaining the initial calibrated $\widehat{\lambda}_{0}^{(\mathrm{p})}$ and then calculate
Consistent with our theory, in Table (ref) we find that Slasso with the calibrated $\lambda$ performs well. In both DGP1 and DGP2 we observe that the prediction error and parameter estimation error decrease as $n$ increases. Similar error reduction is observed under CV $\lambda$. On the contrary, the simulation evidence suggests possible inconsistency of Plasso under either the calibrated $\lambda$ or the CV $\lambda$. Similar patterns are found in terms of mean absolute prediction error (MAPE) and mean absolute parameter estimation error in Table (ref) in the Appendix.
To better understand the unsatisfactory performance of Plasso, Table (ref) in the Appendix shows the percentage of variables selected from the active and inactive coefficients. For example, under CV $\lambda$ Plasso selects fewer active $\beta^{*}$ than Slasso, and the gap is particularly big when $n=360$. In the meantime, it makes more mistakes in selecting the inactive $\beta^{*}$. Plasso faces an inherent dilemma concerning the suitable tuning parameter levels of the stationary and nonstationary components, as discussed in Remark (ref). Nonstationary variables have larger variations and are more influential in prediction. In order to achieve the variable selection effect amongst the nonstationary variables, Plasso requires a large $\lambda$ as in ((ref)). Such emphasis in the nonstationary component imposes a heavy cost in the stationary component, where a non-trivial proportion of the active $\gamma^{*}$ is eliminated, although it also rules out almost all the inactive $\gamma^{*}$. This observation echoes the discussion in Remark (ref) about variable selection in relation to the scales.
The scale normalization in Slasso balances the two types of time series, which allows it to choose active variables in both components. Slasso improves upon Plasso in $\beta^{*}$ for both the active and inactive ones. Moreover, in Table (ref) Slasso produces nearly perfect variable selection in the active $\gamma^{*}$, and in the meantime it controls the estimation error in the inactive one, as shown in Table (ref) in the Appendix about the parameter estimation RMSE of each subset of the coefficients. In this table, the most prominent estimation error comes from the active $\gamma^{*}$ by Plasso, where its large $\lambda$ that accommodates the nonstationary component results in substantial shrinkage bias in the estimation.
For completeness, we check LASSO's performance under the prototype pure unit root case. We consider the same set of $n$ and $p_{x}$ following ((ref)):
where we simply remove all the stationary regressors from DGP1 and 2, respectively. The innovation of $v_{t}=(e_{t}^{\top},u_{t}^{\top})^{\top}$ is again generated according to ((ref)), where $\Omega_{ij}$ is update to $0.8^{|j-j'|}$ for all $i,j\in[p_{x}]$ to allow correlation between $e_{t}$ and $u_{t}$. Table (ref) reports the RMSPE and Table (ref) in the Appendix displays MAPE. Plasso is slightly stronger than Slasso, reflecting the tighter rates of convergence in Theorem (ref) than those in Theorem (ref) as Slasso involves extra randomness in $\widehat{\sigma}_{j}$.
Faced with multiple regressors, some applied econometricians may be inclined to avoid nonstationary regressors in view of the resulting nonstandard asymptotic inference; they may prefer transforming them into stationary ones. Whether we use the stationarized variable or the nonstationary original variable count on the true DGP. The advantage of nonstationary data arises from the super-consistency as the large variation of the stochastic trend can accelerate the rate of convergence, making the parameter estimation more accurate and thereby improving prediction. There is little compelling justification for excluding nonstationary variables a priori in predictive regressions.
We use the FRED-MD macroeconomic database McCracken2016 to check the predictability of the unemployment rate of the United States. Given that the data cover 6 decades from 1960:Jan to 2019:Dec, we adopt a rolling window of length 10 years, 20 years or 30 years, and to make the results comparable we set the entire testing sample as 1990:Jan to 2019:Dec.
The dependent variable, labeled as UNRATE in the database, is plotted in Figure (ref). It ranges from 3.4% to 10.8%, and peaks in the early 1980s recession and the 2008 Global Financial Crisis. It is a persistent time series. If we run a simple AR(1) regression in the entire sample, the AR coefficient 0.995 is close to unity.
We include as potential predictors all the other 121 variables in the database which have no missing values during the sample period. Each variable in FRED-MD is accompanied with a transformation code (TCODE), which suggests a way to transform the raw sequence into a stationary time series. There are 7 categories in total. For a generic scalar time series $(w_{t})_{t=1}^{n}$, the labels 1--7 correspond to the following transformations: (1) null (10 variables out of our 121 predictors); (2) $\Delta w_{t}$ (17 variables); (3) $\Delta^{2}w_{t}$ (none); (4) $\log(w_{t})$ (10 variables); (5) $\Delta\log(w_{t})$ (50 variables); (6) $\Delta^{2}\log(w_{t})$ (33 variables); (7) $\Delta(w_{t}/w_{t-1}-1)$ (1 variable). UNRATE is classified into (2). Figure (ref) draws a representative time series in each category. Obviously the dynamic patterns vary substantively. For example, the one labeled (5) exhibits a clear upward trend, the one labeled (6) shows exponential acceleration, and the one labeled (7) has a dramatic structural break after 2008. One option to avoid nonstationary time series is to stationarize all the raw sequences according to the TCODE. We call this practice stationarization transformation (ST).
Figure (ref) compares the scale of the variables with no transformation (NT) and those with ST. Each dot on the left panel represents the sample s.d. $\widehat{\sigma}_{j}$ (in logarithm base 10 along the y-axis) of each variable, ordered from lower to high for every TCODE marked along the x-axis. The right panel shows the histogram of all variables (the axis again in logarithm base 10). We observe enormous diversity in the upper sub-figures where the data are at their original scales. The smallest sample s.d. is about $10^{-2}$ and the biggest goes over $10^{6}$. Large-scale variables are particularly common in TCODE (4) and (5). Under ST, the variables are much more concentrated. In particular, ST pulls down considerably the scale of all variables that need to be stationarized.
We conduct 1, 2, or 3-month ahead out-of-sample prediction, denoted by $h=1,2$, or $3$. We set two simple benchmark models: (i) Random walk with drift (RWwD), where $\widehat{y}_{n+h}=y_{n}+\frac{h}{n}(y_{n}-y_{0})$; and (ii) AR model $\widehat{y}_{n+h}=\widehat{\pi}_{0,h}+\widehat{\pi}_{1,h}y_{n}+\cdots+\widehat{\pi}_{q,h}y_{n-q+1}$ where the AR coefficients are estimated by OLS and the number of lags $q$ is determined by the Bayesian information criterion. All these models use information up to time $n$. Table (ref) shows RMSPE averaged over the entire testing sample 1990:Jan--2019:Dec, and three testing sub-samples for each decade. Across the lengths of the rolling windows, a 30-year rolling window does not necessarily improve RMSPE, indicating potential model uncertainty over a long training sample. Across the testing sub-samples, RMSPE is the largest during 2000--2010, which includes the Great Recession. Across the forecast horizons, the estimation error increases as the horizon gets farther in the future.
When implementing LASSO, we use the data-driven 10-fold CV as introduced in Section (ref). Regarding the potential regressors, we first throw all the 121 predictors into the linear regression. With the rich mix of time series of various temporal patterns, we find in Table (ref) that Plasso is much worse than Slasso, and Slasso under NT outperforms the best benchmark model in most cases. If we transform all variables to ST according to TCODE, the outcomes deteriorate. Relative performance is similar when the error is measured by MAPE, reported in Table (ref) in the Appendix. These empirical results echo smeekes2020unit, who find that the best forecast strategy should be devised based on the nature of the target time series as well as the predictors; the information contained in NT often has an edge over ST.\footnote{smeekes2020unit carry out empirical exercises targeting several variables in FRED-MD and they also use Google Trend to nowcast Dutch unemployment. These empirical applications are elaborated in smeekes2018macroeconomic and smeekes2021automated.}
We plot the prediction errors $(\widehat{y}_{n+1}-y_{n+1})$ in Figure (ref) under NT and ST with $h=1$ and 20-year rolling window. The graphs are similar under other $h$ and rolling window lengths. The prediction errors based on Slasso with NT (upper right panel) fluctuate around 0 in a narrow range between $\pm0.5$, as the 121 regressors form a linear combination that predicts well the one-month-ahead unemployment rate. The errors produced by Plasso (upper left panel) remain persistent, swinging wildly between 2008 and 2016. Under ST (lower panels) the outliers are not as pronounced, but the prediction errors go beyond the range of $\pm1.2$ and appear persistent. Furthermore, under the same $h$ and $n$ we check LASSO's selected variables under NT. The FRED database classifies all time series into 8 categories based on economic implications, and UNRATE belongs to the Labor Market (LM) group. Table (ref) reports the top 10 most frequently selected variables over the rolling window estimation. Among these 10 variables that Slasso picks out, 9 are from LM group, which showcase the economic relevance of the variables chosen by Slasso. In contrast, Plasso is inclined to select the variables with large s.d.
That the linear combination of 121 predictors under Slasso can outperform the benchmarks is encouraging. It illustrates the value of a high dimensional model estimated by an off-the-shelf machine learning method. Next, we experiment with an augmented model. stock2002forecasting propose computing diffusion indices---the principal components from many potential predictors, and bai2008forecasting further add lagged dependent variables into predictive regressions. Following medeiros2021forecasting, we incorporate the lagged dependent variable and four diffusion indices, making 126 unique regressors, and to allow potential delayed effects we include four time lags of each predictor, totaling $126\times4=504$ regressors.
Columns under “504 predictors” in Table (ref) show the corresponding RMPSE. While Plasso remains worse than the simple benchmarks, we observe improvement in Slasso. First, under ST the additional lagged dependent variables and diffusion indices mitigate the imbalance in the predictive regression and therefore strengthen the performance of Slasso. Moreover, these additional regressors improve Slasso under NT, which is the overall best performer. It reduces the RMPSE in 29 out of the 36 instances relative to the counterpart with “121 predictors”, and beats the last column in most instances except $h=3$, where the errors mainly occur during 2000--2009. These results indicate that macroeconomic domain knowledge is instrumental in guiding the initial specification to determine the pool of regressors, and then Slasso takes care of the estimation of many coefficients. This fusion of field expertise and machine learning is more effective than simply LASSOing with all variables in the database.
This paper studies asymptotic properties of LASSO in predictive regressions where many nonstationary time series are present. We establish new bounds for the RE, which allows us to derive convergence rates for Plasso and Slasso. The consistency of Slasso is extended to the model of mixed stationary and nonstationary regressors, and it can further digest information from cointegrated variables. The simulations and the empirical application provide numerical evidence that supports the merits of Slasso, which we recommend for practice.
As a first step of exploration, this paper uses the unit root process as a representative of nonstationary time series. There are other popular models that characterize persistence, for example, local-to-unity and fractional integration. Future investigation of these nonstationary time series will generalize the theory and further guide practical implementation. Moreover, the theoretical results of this paper rely on a tuning parameter expanding at some rate based on the same size. It will be important to explore the behaviors of LASSO involving nonstationary regressors when the tuning parameter is selected by a data-driven method, such as the CV. Last but not least, a formal testing procedure for coefficients will be feasible in high dimensions if we debias the LASSO estimator.
\setcounter{footnote}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{equation}{0} \setcounter{thm}{0} \setcounter{lem}{0} \setcounter{rem}{0} \setcounter{cor}{0} \setcounter{prop}{0}