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.
135,141 characters · 23 sections · 94 citation commands
LASSO Inference for High Dimensional Predictive Regressions
\thispagestyle{empty}
Key words: Data-rich environment, Forecast, Hypothesis testing, LASSO, Local unit root, Shrinkage
JEL code: C22, C53, C55
The evaluation of economic and financial predictability has attracted widespread interest from theoretical and applied researchers for many decades. In today's era of big data, we have unprecedented access to a vast amount of digital information about the economy. Recent advancements in inference with high dimensional data have uncovered new empirical patterns in predictive practices using large-scale datasets with temporal features.
This paper aims at a plain quest: in a high dimensional linear predictive regression model, where the number of potential regressors is larger than the sample size, how can one conduct valid inference for a regressor of primary interest? No research has solved this question before. The challenges are twofold. First, predictive regressions were mainly studied in the low-dimensional context. A defining feature of predictive regression theory lies in persistent regressors stambaugh1999predictive, which upend standard inference based on the standard $t$-statistics. Second, one must estimate the coefficient via some regularization methods to cope with high dimensionality. For example, when the underlying true regression model is sparse, LASSO tibshirani1996regression is the off-the-shelf method. It is well known that the LASSO estimator is biased toward zero due to absolute-value shrinkage and exhibits a nonstandard asymptotic distribution that is distinct from the normal distribution. If we intend to provide an asymptotically normally distributed estimator to facilitate standard statistical inference, we must simultaneously combat two evils: the Stambaugh bias due to persistent regressors and the shrinkage bias caused by the LASSO penalty.
The above diagnosis hints at a plausible solution path. In low-dimensional predictive regressions, phillips2009econometric's IVX method leverages a self-generated instrument to alleviate the regressor persistence, thereby overcoming the Stambaugh bias. In high dimensional cross sectional regressions, zhang2014confidence's desparsified LASSO (Dlasso) constructs a score vector for the parameter of interest and removes the shrinkage bias via an auxiliary regression. Each method provides an asymptotic normal estimator in its respective environment.
Can we combine these two methods into a single procedure to address inference in high dimensional predictive regressions? We find that the answer is both no and yes. “No” is in the sense that a naive combination of the two does not lead to desirable results. “Yes”, on the other hand, is established upon a deep understanding of the mechanisms of both components and their adaptation to the context. This research culminates in a new IVX-Desparsified LASSO (XDlasso) estimator that is free from both biases and has an asymptotic normal distribution.
With a predictor of interest in mind, the construction of XDlasso, detailed in Algorithm (ref) in Section (ref), is summarized as follows. First, a workhorse estimator is needed to lay the groundwork for a high dimensional predictive regression. When both nonstationary and stationary predictors are present in the regression, mei2022lasso has recently established the consistency of the standardized LASSO (Slasso), making it a natural candidate for the workhorse estimator. Beyond the pure unit roots considered in MS24, we extend the characterization of nonstationary regressors by generalizing our framework to allow for local unit root (LUR) processes. The convergence rates of Slasso for LURs align with those for pure unit roots. However, the technical proofs for LURs are more involved than those for pure unit roots. We address the complexity arising from high dimensional predictors with both LURs and stationary regressors, and derive the convergence rates of the initial Slasso estimator.
Second, the common practice of generating the instrument in IVX is insufficient --- the IV must be scale-standardized to have the stochastic order aligned with all other predictors in Slasso. The standardized IV serves as the target variable for the auxiliary Slasso regression in Dlasso to estimate the shrinkage bias. The XDlasso estimator of the parameter of interest is defined as the initial Slasso estimator plus the bias-correction term, and the companion $t$-statistic is employed for statistical inference by comparing it to critical values from the standard normal distribution.
We further establish the asymptotic normality of our proposed XDlasso estimator and the convergence rate of its standard error. Specifically, the XDlasso estimator is $\sqrt{n}$-consistent for a stationary regressor while its convergence accelerates for an LUR regressor. Moreover, to conduct simultaneous inference for multiple parameters of interest, we develop a Wald statistic with an asymptotic $\chi^{2}$ distribution based on XDlasso. This Wald test is valid even when the parameters involve both stationary and nonstationary regressors.
To tackle persistent regressors, the self-generated IVX instrument in the second step is the key ingredient. The generated IV is less persistent than the nonstationary regressors modeled as LURs. This important feature enables us to decorrelate the IV from other covariates in the auxiliary LASSO regression, so that the resulting XDlasso estimator possesses these two properties: (i) It is free from the Stambaugh bias and thus enjoys asymptotic normality; (ii) It reduces the order of shrinkage bias to make it correctable. In contrast, the ordinary Dlasso encounters a spurious auxiliary regression, failing to correct the bias arising from persistent regressors (see Section (ref) for details). More importantly, XDlasso inference does not require a priori knowledge of the persistence of the regressor of interest and is thus immune to pretesting bias. To the best of our knowledge, this is the first methodology to handle the inferential problem in high dimensional predictive regressions ($p\gg n$) with nonstationary predictors. This is also the first paper that extends the IVX technique into the high dimensional framework.
Monte Carlo simulations show that XDlasso successfully removes the bias for inference on the coefficient of a nonstationary regressor, but the ordinary Dlasso fails to do so. Our procedure is applied to the high dimensional macroeconomic FRED-MD dataset mccracken2016fred with both stationary and persistent variables, to study two important macro-finance problems: financial market return predictability and the Phillips curve in macroeconomics.
Literature review. With the advent of big data, machine learning methods have spread to time series topics such as nonstationarity phillips2021boosting,smeekes2021automated,mei2024boosted, cointegration testing onatski2018alternative,zhang2019identifying,bykhovskaya2022cointegration, and structural breaks deshpande2023online,tu2023penetrating. This paper builds on several strands of literature. First, LASSO is one of the most studied methods in recent years, with well-developed theory in high dimension bickel2009simultaneous. It is well received and used for economic applications; see belloni2012sparse, shi2016estimation, caner2018asymptotically, and babii2022machine, to name a few. In recent years, the properties of LASSO are studied in various topics in high dimensional time series, including nonstationary time series models koo2020high,lee2022lasso and inference based on the heteroskedasticity and autocorrelation consistent (HAC) estimation babii2021highdimensionalgrangercausalitytests,babii2024high. None of these works has considered hypothesis testing problem for high dimensional predictive regressions with both LURs and stationary regressors.
Hypothesis testing after LASSO is challenging because of the shrinkage bias. To validate hypothesis testing in high dimensions, zhang2014confidence, van2014asymptotically, and javanmard2014confidence have developed the desparsified (debiased) LASSO estimators under the independently and identically distributed (i.i.d.) setting. adamek2022lasso generalize the Dlasso inference to high dimensional stationary time series. We follow this line of desparsified LASSO literature thanks to its convenience, which requires a baseline regression and an auxiliary regression only. On the other hand, chernozhukov2018's double machine learning (DML) is a more general theoretical framework of debiased inference, widely used in cross sectional data where sample-splitting is readily implementable. However, none of the aforementioned works has devised any inferential procedure for high dimensional nonstationary time series. hecq2023inference apply a post-double selection procedure to test the Granger causality in high dimensional nonstationary vector autoregressive models with cointegrated data. In contrast, our procedure relies on desparsified LASSO without variable selection.
The other strand is the vast literature on predictive regressions. As highlighted by campbell2006efficient and jansson2006optimal, non-standard distortion in the asymptotic distribution arises from persistent regressors. The peculiar asymptotic distributions invalidate the standard inferential procedures. There have been multiple proposals for valid inference, for example, the Bonferroni method campbell2006efficient, the conditional likelihood method jansson2006optimal, the linear projection method cai2014testing, the weighted empirical likelihood approach zhu2014predictive,liu2019unified,yang2021unified, and the implication-based inference xu2020testing. Some of these methods are designed for univariate predictive regressions; it would be difficult to extend them to the high dimensional case, where regularization is required to handle many parameters. On the other hand, phillips2009econometric's IVX estimator gained its popularity by recovering asymptotic normality, enabling valid inference for mean regressions kostakis2015robust,kostakis2018taking,phillips2013predictive,phillips2016robust,yang2020testing,demetrescu2023extensions and quantile regressions lee2016predictive,fan2019predictive,cai2023new,liu2023unified with low dimensional regressors. IVX recovers asymptotic normality by projecting the persistent regressor onto a self-generated IV.
Layout. The rest of the paper is organized as follows. Section (ref) introduces the high dimensional predictive regression model with a mixture of stationary and nonstationary regressors and proposes XDlasso. Section (ref) establishes the theoretical results, justifying the size and power of the XDlasso inference procedure. Section (ref) carries out simulation studies that corroborate the theory. Section (ref) applies XDlasso inference to two macro-finance empirical examples. Technical proofs are relegated to the Online Appendices.
Notations. We set up the notation before the formal discussion. We define $\boldsymbol{1}\left\{ \cdot\right\} $ as the indicator function, and $\Delta$ as the difference operator so that $\Delta x_{t}=x_{t}-x_{t-1}$. The set of natural numbers, integers, and real numbers are denoted as $\mathbb{N}$, $\mathbb{Z}$, and $\mathbb{R}$, respectively. For some $n\in\mathbb{N}$, the integer set $\{1,2,\cdots,n\}$ is denoted as $[n]$, and the space of $n$-dimensional vectors is denoted as $\mathbb{R}^{n}.$ For $x=(x_{t})_{t\in[n]}\in\mathbb{R}^{n}$, the $L_{1}$-norm is $\left\Vert x\right\Vert _{1}=\sum_{t=1}^{n}\left|x_{t}\right|$, and the sup-norm is $\|x\|_{\infty}=\sup_{t\in[n]}|x_{t}|$. Let $0_{n}$ be an $n\times1$ zero vector, $1_{n}$ be an $n\times1$ vector of ones, and $I_{n}$ be the $n\times n$ identity matrix. For a generic matrix $B,$ let $B_{ij}$ be the $(i,j)$-th element, and $B^{\top}$ be its transpose. Let $\|B\|_{\infty}=\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”. We use $\stackrel{\mathrm{p}}{\to}$ and $\stackrel{\mathrm{d}}{\to}$ to denote convergence in probability and in distribution, respectively. For any time series $\{a_{t}\}_{t=1}^{n}$, we use $\bar{a}$ to denote its sample mean $n^{-1}\sum_{t=1}^{n}a_{t}$. For any time series $\{a_{t}\}$ and $\{b_{t}\}$, we say they are asymptotically uncorrelated if their sample correlation coefficient $\frac{\sum_{t=1}^{n}(a_{t}-\bar{a})(b_{t}-\bar{b})}{\sqrt{\sum_{t=1}^{n}(a_{t}-\bar{a})^{2}\sum_{t=1}^{n}(b_{t}-\bar{b})^{2}}}\stackrel{\mathrm{p}}{\to}0$ as $n\to\infty.$
Suppose that a time series of the outcome $y_{t}$ is generated by the following linear predictive regression:\footnote{For simplicity of exposition, an intercept in (ref) is omitted, without loss of generality. As explained by MS24, the intercept in LASSO can be handled by the well-known Frisch-Waugh-Lovell theorem. In practical implementation --- throughout all simulations and empirical exercises in this paper --- we keep an unpenalized intercept in the model.}
where the error term $u_{t}$ is a stationary martingale difference sequence (m.d.s.) with mean zero and conditional variance $\sigma_{u}^{2}$. We consider two types of regressors with different stochastic properties. Firstly, the $p_{x}\times1$ vector $X_{t}=(x_{1,t},\cdots,x_{p_{x},t})^{\top}$ collects the LURs:
where $e_{t}=(e_{1,t},\cdots,e_{p_{x},t})^{\top}$ is a $p_{x}$-dimensional vector of stationary time series. The AR(1) coefficient $\rho_{j}^{*}$ in ((ref)) is close to 1 when the sample size $n$ is large, specified as
where $c_{j}^{*}\in\mathbb{R}$ is allowed to be negative, positive, or zero. Therefore, our framework accommodates nonstationary regressors that are locally integrated ($c_{j}^{*}<0$), unit roots ($c_{j}^{*}=0$), and locally explosive ($c_{j}^{*}>0$). Secondly, stationary regressors are stored in the $p_{z}\times1$ vector $Z_{t}=(z_{1,t},\cdots,z_{p_{z},t})^{\top}$. The two types of regressors are combined into a long vector $W_{t}=(X_{t}^{\top},Z_{t}^{\top})^{\top}=(w_{1,t},\cdots,w_{p,t})^{\top}$ of $p$ ($=p_{x}+p_{z}$) elements, and the associated coefficients are placed into $\theta^{*}=(\beta^{*\top},\gamma^{*\top})^{\top}\in\mathbb{R}^{p}$. Following the literature, we refer to $W_{t}$, which has multiple degrees of persistence, as mixed root regressors. For simplicity, let the initial value $\|W_{t=0}\|_{\infty}=O_{p}(1)$. Define the sparsity index $s=\sum_{j=1}^{p}\boldsymbol{1}\{\theta_{j}^{*}\neq0\}$ as the number of nonzero components in the coefficient vector $\theta^{*}$.
As in the default R program option glmnet::glmnet(x,y), it is a common practice in LASSO to scale-standardize each regressor $w_{j,t}$ by its sample standard deviation (s.d.) $\widehat{\sigma}_{j}=\sqrt{\frac{1}{n}\sum_{t=1}^{n}(w_{j,t-1}-\bar{w}_{j})^{2}},$ where $\bar{w}_{j}=n^{-1}\sum_{t=1}^{n}w_{j,t-1}$ is the sample mean. Let the diagonal matrix $D={\rm diag}\left(\widehat{\sigma}_{1},\widehat{\sigma}_{2},\cdots,\widehat{\sigma}_{p}\right)$ store the sample standard deviations. The standardized LASSO (Slasso) estimator is
The Slasso estimator is scale-invariant: if the regressor $w_{j,t-1}$ is multiplied by a nonzero constant $m$, then the $j$-th coefficient estimator changes proportionally to $\widehat{\theta}_{j}^{{\rm S}}/m.$ The standardization renders the magnitudes of LURs into the same order as those of stationary regressors, so that the same LASSO tuning parameter $\lambda$ in ((ref)) is valid for both stationary and persistent regressors. In contrast, the plain LASSO (Plasso) with the matrix $D$ in ((ref)) replaced with the identity matrix is scale-variant. What is worse, equipped with a single tuning parameter $\lambda$, Plasso favors the LURs of a larger order, and shrinks the coefficients of stationary regressors with a smaller order all the way to zero --- thus becoming inconsistent.
When data are i.i.d. or stationary, LASSO is subject to shrinkage bias as well as a nonstandard asymptotic distribution, which cannot be used for standard inference fu2000asymptotics. This motivates zhang2014confidence to bring forth the desparsified LASSO to correct the bias and recover asymptotic normality. With high dimensional LURs, $\widehat{\theta}^{{\rm S}}$ is subject to not only the shrinkage bias, but also the Stambaugh bias due to persistence, which further distorts the standard t-statistics inference. In this section, we examine both the shrinkage bias and the Stambaugh bias, and propose XDlasso for correcting both biases.
We are interested in inference on a null hypothesis $\mathbb{H}_{0}:\theta_{j}^{*}=\theta_{0,j}$ for a $j\in[p]$, a prevalent practice in empirical studies. In a low dimensional linear regression where $p$ is fixed, the Frisch-Waugh-Lovell theorem yields the following formulation of the ordinary least squares (OLS) estimator: \[ \widehat{\theta}_{j}^{{\rm OLS}}=\frac{\sum_{t=1}^{n}w_{j,t-1}^{\perp}y_{t}}{\sum_{t=1}^{n}w_{j,t-1}^{\perp}w_{j,t-1}}, \] where $w_{j}^{\perp}=(w_{j,0}^{\perp},\dots,w_{j,n-1}^{\perp})^{\top}$ is the OLS residual from regressing $w_{j,t}$ on all other regressors $W_{-j,t}=(w_{k,t})_{k\neq j}$. OLS induces a large variance as $p$ gets large, and becomes infeasible when $p>n$. Now, consider replacing the OLS residual $w_{j}^{\perp}$ by a generic score vector $r_{j}=(r_{j,0},\dots,r_{j,n-1})^{\top}$ to construct an estimator of $\theta_{j}^{*}$ that is linear in $y_{t}$ in the form of $\widehat{\theta}_{j}^{({\rm lin})}=\frac{\sum_{t=1}^{n}r_{j,t-1}y_{t}}{\sum_{t=1}^{n}r_{j,t-1}w_{j,t-1}}.$ Since \[ y_{t}=w_{j,t-1}\theta_{j}^{*}+W_{-j,t-1}^{\top}\theta_{-j}^{*}+u_{t} \] where $\theta_{-j}^{*}=(\theta_{k}^{*})_{k\neq j}$ is the vector of coefficients excluding the $j$-th entry, the generic estimator $\widehat{\theta}_{j}^{({\rm lin})}$ can be decomposed into \[ \widehat{\theta}_{j}^{({\rm lin})}=\theta_{j}^{*}+\frac{\sum_{t=1}^{n}r_{j,t-1}u_{t}}{\sum_{t=1}^{n}r_{j,t-1}w_{j,t-1}}+\frac{\sum_{t=1}^{n}r_{j,t-1}W_{-j,t-1}^{\top}\theta_{-j}^{*}}{\sum_{t=1}^{n}r_{j,t-1}w_{j,t-1}}=:\theta_{j}^{*}+{\rm N}_{j}+{\rm B}_{j}, \] where ${\rm N}_{j}$ is the noise component that determines the asymptotic distribution of $\widehat{\theta}_{j}^{({\rm lin})}$, and ${\rm B}_{j}$ is the potential bias due to the choice of $r_{j}$.
For OLS, the score vector $r_{j}=w_{j}^{\perp}$ is orthogonal to the column space of $W_{-j,\cdot}:=(W_{-j,0},\dots,W_{-j,n-1})^{\top}$, under which ${\rm B}_{j}=0$ and no bias is present. The bias term ${\rm B}_{j}$ pops up whenever $r_{j}$ is not orthogonal to $W_{-j,\cdot}$, which happens if we add a penalty to the OLS objective function. With the LASSO penalty at place, we call ${\rm B}_{j}$ shrinkage bias.
Following zhang2014confidence, we replace the unknown parameter $\theta_{-j}^{*}$ by the feasible workhorse Slasso estimator $\widehat{\theta}_{-j}^{{\rm S}}=(\widehat{\theta}_{k}^{{\rm S}})_{k\neq j}$ to obtain
to compensate for ${\rm B}_{j}$. Equivalently, $\widehat{\theta}_{j}$ can be written as
where $\widehat{u}_{t}=y_{t}-W_{t-1}^{\top}\widehat{\theta}^{{\rm S}}$ is the Slasso residual. Though LASSO may shrink $\widehat{\theta}_{j}^{{\rm S}}$ all the way to exactly zero, the second term in ((ref)) is continuously distributed and therefore $\widehat{\theta}_{j}$ is a desparsified version of $\widehat{\theta}_{j}^{{\rm S}}$. A straightforward calculation yields
where $\widehat{{\rm B}}_{j}-{\rm B}_{j}$ is the approximation error of the shrinkage bias.
Let $\omega_{j}$ denote the standard deviation of ${\rm N}_{j}$. To secure asymptotic normality for $\widehat{\theta}_{j}$, we need an appropriate score vector $r_{j}$ such that as $n\to\infty$:
These two results will furnish $\widehat{\theta}_{j}$ with asymptotic normality, validating inference based on the standard $t$-statistic.
The result (R1) is well established by zhang2014confidence under i.i.d. data when the score $r_{j,t}$ is taken as the LASSO residual of regressing $w_{j,t}$ on all other regressors. However, when $w_{j,t}$ is an LUR, this practice leads to a highly persistent $r_{j,t}$ and ruins the asymptotic normality of ${\rm N}_{j}$. This is the Stambaugh bias --- a highly persistent regressor produces an asymptotically non-normally distributed OLS estimator skewed away from zero. To safeguard (R1), we must seek a score that is less persistent than an LUR.
To retain (R2), we again examine the estimation error of the shrinkage bias
Note that $|{\rm \widehat{B}}_{j}-{\rm B}_{j}|\leq\|g_{j}\|_{\infty}\|\widehat{\theta}_{-j}^{{\rm S}}-\theta_{-j}^{*}\|_{1}.$ Slasso's $L_{1}$ estimation error $\|\widehat{\theta}_{-j}^{{\rm S}}-\theta_{-j}^{*}\|_{1}$ is invariant with the choice of score vector $r_{j,t}$, and diminishes if $\widehat{\theta}^{{\rm S}}$ is consistent. Thus, it suffices to control the order of $\|g_{j}\|_{\infty}$ to achieve (R2). The explicit expression of $g_{j}$ in ((ref)) suggests that a weak correlation in sup-norm between $W_{-j,t}$ and the score $r_{j,t}$ (relative to the correlation between $w_{j,t}$ and $r_{j,t}$) helps.
With these routes in mind, we devise XDlasso in the following section. It proceeds with two key steps: (i) conducting an IVX transformation to get a new variable less persistent than the LUR regressors to eliminate the Stambaugh bias; (ii) running an auxiliary LASSO regression to construct an $r_{j,t}$ to remove the shrinkage bias.
For inference of $\theta_{j}^{*}$ in low dimensions, phillips2009econometric's IVX method generates an instrument by quasi-differencing $w_{j,t}$. When $w_{j,t}$ is an LUR, this self-generated instrument is mildly integrated. The mitigation of persistence will remove the Stambaugh bias when the sample size passes to infinity, and thus recovers asymptotic normality of the IVX estimator. Due to the coexistence of Stambaugh bias and shrinkage bias, the wisdom of IVX cannot be directly transplanted into the high dimensional case. Instead, we integrate IVX with the idea of desparsified LASSO after comprehending the mechanism of the biases illustrated in Section (ref).
Specifically, IVX adopts the following instrumental variable:
where $\rho_{\zeta}\in(0,1)$ is a user-determined tuning parameter. Define the s.d. of the instrument as $\widehat{\varsigma}_{j}=\sqrt{n^{-1}\sum_{t=0}^{n-1}\left(\zeta_{j,t}-\bar{\zeta}_{j}\right)^{2}.}$ We unify the scale for LURs and stationary regressors by standardizing the instrument with its s.d.:
For low dimensional predictive regressions, the IVX literature has established asymptotic normality of ${\rm N}_{j}$ in ((ref)) taking $r_{j,t}=\widetilde{\zeta}_{j,t}$; see phillips2009econometric. When $w_{j,t}$ is an LUR, the IV $\zeta_{j,t}$ is mildly integrated and less persistent than an LUR. Furthermore, recall that the vector of regressors $W_{-j,t}$ includes either LURs or stationary regressors. Due to different degrees of persistence, the mildly integrated IV $\widetilde{\zeta}_{j,t}$ and the regressors $W_{-j,t}$ are asymptotically uncorrelated. Thus, when $w_{j,t}$ is an LUR, we can choose the score vector as $r_{j,t}=\widetilde{\zeta}_{j,t}$ to deliver a small order of $\|g_{j}\|_{\infty}$. If $w_{j,t}$ is stationary, however, this score vector fails. When $n$ is large, the instrumental variable $\zeta_{j,t}$ behaves similarly as the stationary regressor $w_{j,t}$, and thus its correlation to the high dimensional stationary regressors in $W_{-j,t}$ is not negligible.
The above analysis implies that we must decorrelate the score vector with the other regressors $W_{-j,t}$ to reduce the order of $g_{j}$ to control for the approximation error $\widehat{{\rm B}}_{j}-{\rm B}_{j}$. This decorrelation will produce a unified testing approach for both LURs and stationary regressors. To this end, we construct a residual score vector $\widehat{r}_{j}=(\widehat{r}_{j,0},\cdots,\widehat{r}_{j,n-1})^{\top}$ by the following auxiliary LASSO regression
with the LASSO tuning parameter $\mu_{j}$ and $D_{-j}={\rm diag}(\{\widehat{\sigma}_{k}\}_{k\neq j})$. In low dimensional multivariate predictive regressions, IVX transforms each regressor into a less persistent instrumental variable parallel to ((ref)), and constructs a two-stage least squares estimator using all the self-generated instrumental variables. In contrast, we only transform the variable of interest $w_{j,t}$ and estimate one auxiliary regression ((ref)).
The score vector $\widehat{r}_{j}$ in ((ref)) accommodates stationary and LUR regressors. Recall that when $w_{j,t}$ is stationary, the magnitude of the instrument $\zeta{}_{j,t}$ behaves similarly as $w_{j,t}$. Thus, the score $\widehat{r}_{j}$ is asymptotically equivalent to the standardized residual of the LASSO regression of $w_{j,t}$ on $W_{-j,t}$. The latter is proportional to the score in zhang2014confidence for cross-sectional data, which is also used in adamek2022lasso for stationary time series. When $w_{j,t}$ is an LUR, the instrument $\zeta_{j,t}$ is mildly integrated and has a different degree of persistence from $W_{-j,t}$. Therefore, $\zeta_{j,t}$ is asymptotically uncorrelated to each regressor in $W_{-j,t}$. For notational conciseness, define $\widetilde{W}_{-j,t}=(D_{-j})^{-1}W_{-j,t}$ and $\widetilde{\varphi}^{(j)}=D_{-j}\widehat{\varphi}^{(j)}$. The analysis above suggests $\widetilde{\varphi}^{(j)}\approx0$ so that $\widehat{r}_{j,t}=\widetilde{\zeta}_{j,t}-\widetilde{W}_{-j,t}^{\top}\widetilde{\varphi}^{(j)}\approx\widetilde{\zeta}_{j,t}$. Recall from the discussion right after ((ref)) that the standardized instrument $\widetilde{\zeta}_{j,t}$ is a valid score process for an LUR regressor, and thus the score $\widehat{r}_{j,t}$ can also remove the biases asymptotically. As a result, the residual of the auxiliary LASSO regression provides a unified construction of the score for either stationary or nonstationary $w_{j,t}$. It allows practitioners to conduct hypothesis testing on the coefficient in high dimensional predictive regression regardless of the order of integration of $w_{j,t}$. We maintain an agnostic attitude about the persistence of the regressors; in practical implementation, we need not distinguish these two types.
Following ((ref)), XDlasso is constructed as
with the standard error
where $\widehat{\sigma}_{u}^{2}=n^{-1}\sum_{t=1}^{n}\widehat{u}_{t}^{2}$. For simplicity, we focus on homoskedastic errors in the main text. In Section B.3 of the appendix, we provide a formula of the heteroskedasticity-robust standard error, and show that our method is robust to conditional heteroskedasticity in both Monte Carlo simulations (Section B.3) and empirical applications (Section C.2).
With the point estimator ((ref)) and the associated standard error in ((ref)), we perform the $t$-test for the null hypothesis $\mathbb{H}_{0}:\theta_{j}^{*}=\theta_{0,j}$. We summarize the testing procedure using the XDlasso estimator below.
\singlespacing \fbox{
}
\onehalfspacing
The testing procedure in Algorithm (ref) refines the conventional predictive regression inference using modern high dimensional inference techniques. In the following, we elaborate on the necessity of the IVX transformation by explaining the drawback of Dlasso under high dimensional LURs.
Given i.i.d. data, Dlasso proceeds with the following residual as the score vector:
where the dependent variable in ((ref)) is also scaled to be kept in line with ((ref)) and ((ref)). For simplicity, in this subsection we temporarily restrict all nonstationary regressors to be pure unit roots ($\rho_{j}^{*}=1$). We show that the standard Dlasso procedure fails to correct the bias in this special case of LURs, thereby invalidating the inference in general cases.
First, ((ref)) is a spurious regression as in granger1974spurious. In a low dimensional regression where all regressors have unit roots, the standardized least squares counterpart of the regression ((ref)) follows a functional central limit theorem (FCLT)
where $\mathcal{B}=(\mathcal{B}_{1},\dots,\mathcal{B}_{p})^{\top}$ is a $p$-dimensional Brownian motion, $\sigma_{j}^{*}=\sqrt{\int_{0}^{1}\mathcal{B}_{j}^{2}-(\int_{0}^{1}\mathcal{B}_{j})^{2}}$, and $\mathcal{B}_{-j}=(\mathcal{B}_{k})_{k\neq j}$. In other words, the OLS estimator converges in distribution to a nondegenerate random variable. This is the well-known spurious regression phenomenon in unit root regressions. While ((ref)) suggests that the limit target coefficients are random and continuously distributed, the asymptotics for the LASSO estimator ((ref)) in high dimensions depends on the sparsity of the target coefficients. The randomness of coefficients is contradictory to the sparsity required in LASSO.
Second, the score vector in ((ref)) cannot remove the Stambaugh bias. Even in low dimensions, the least squares residual $\widehat{r}_{j,t}^{{\rm OLS}}=\widetilde{w}_{j,t}-\widetilde{W}_{-j,t}^{\top}\widetilde{\psi}^{(j,{\rm OLS})}$remains highly persistent due to spurious regression. Using $\widehat{r}_{j,t}^{{\rm OLS}}$ as the score vector in low dimensions would therefore keep the Stambaugh bias in the noise component ${\rm N}_{j}$. This issue in low dimension is also present in the score vector $r_{j,t}$ construction in ((ref)) by high dimensional LASSO.
In contrast, with the help of the IVX transformation, the score vector in ((ref)) achieves (R1) and (R2). Figure (ref) provides an illustrative simulation to compare Dlasso with XDlasso. We set $\beta_{1}^{*}$ (the coefficient of the first unit root) as zero, $n=300$, $(p_{x},p_{z})=(150,300)$, and use i.i.d. innovations (See Eq. ((ref))). We generate the data following the DGP ((ref)), with a mixture of stationary and nonstationary regressors as in Section (ref). Figure (ref) displays the histograms of the $t$-statistics over 2000 replications. The density of XDlasso $t$-statistic is well approximated by $\mathcal{N}(0,1)$, whereas the Dlasso $t$-statistic suffers from a substantial bias.
We have devised the XDlasso inference for a scalar coefficient $\theta_{j}^{*}$. In low dimensional predictive regression, the IVX estimator is applicable to multiple coefficients, which jointly follows an asymptotic multivariate normal distribution. Therefore, a Wald statistic by IVX is available to jointly test the predictability of multiple regressors, provided there are a finite number of them. This test statistic is shown to be valid, even when the parameters of interest involve both stationary and nonstationary regressors.
The validity of Wald test extends to XDlasso in high dimensional predictive regression. Specifically, suppose that we are interested in a subset of regressors indexed by $\mathcal{J}\subset[p]$ with a fixed cardinality $|\mathcal{J}|$. In this case, the XDlasso estimators $\widehat{\theta}_{\mathcal{J}}^{{\rm XD}}=(\widehat{\theta}_{j}^{{\rm XD}})_{j\in\mathcal{J}}$ defined in ((ref)) asymptotically follow a multivariate normal distribution. Let the null hypothesis $\mathbb{H}_{0}:\theta_{\mathcal{J}}^{*}=\theta_{0,\mathcal{J}}$ involve $|\mathcal{J}|$ restrictions, where $\theta_{0,\mathcal{J}}$ is a $|\mathcal{J}|$-dimensional vector. We construct the following Wald statistic
where $\widehat{\Omega}_{\mathcal{J}}^{{\rm XD}} = (\widehat{\Omega}_{j,k}^{{\rm XD}})_{j,k\in\mathcal{J}}$ estimates the covariance matrix of $\widehat{\theta}_{\mathcal{J}}^{{\rm XD}}$, with \[ \widehat{\Omega}_{j,k}^{{\rm XD}}=\widehat{\sigma}_{u}^{2}\dfrac{\sum_{t=1}^{n}\widehat{r}_{j,t-1}\widehat{r}_{k,t-1}}{\sum_{t=1}^{n}\widehat{r}_{j,t-1}w_{j,t-1}\cdot\sum_{t=1}^{n}\widehat{r}_{k,t-1}w_{k,t-1}} \] measuring the covariance between $\widehat{\theta}_{j}^{{\rm XD}}$ and $\widehat{\theta}_{k}^{{\rm XD}}$. Under the null hypothesis, the Wald statistic in ((ref)) will follow an asymptotic $\chi^{2}$ distribution with the degree of freedom $|\mathcal{J}|$, enabling joint inference on $\theta_{\mathcal{J}}^{*}$.
This section develops the limit theory to shed light on the asymptotic behaviors of the XDlasso estimator. Unsurprisingly, this paper's assumptions share similarities with those in MS24. We state our theoretical assumptions and then highlight and explain the differences between the assumptions in these two papers. Regarding the asymptotic framework, we define 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 asymptotic statements, we will explicitly send $n\to\infty$, and it is understood that $p(n)\to\infty$ as $n\to\infty$, while $s\left(n\right)$ is allowed to be either fixed or divergent. Recall that the stationary vector $e_{t}\in\mathbb{R}^{p_{x}}$ is the innovation of the LUR regressors. We assume that the stationary high dimensional vector $v_{t}=(e_{t}^{\top},Z_{t}^{\top})^{\top}$ is generated by the innovations $\varepsilon_{t}=(\varepsilon_{k,t})_{k\in[p]}$ via a linear transformation
where $\Phi$ is a $p\times p$ deterministic matrix. Let $\mathcal{F}_{t}$ denote the $\sigma$-field generated by $\{u_{s},\varepsilon_{s}\}_{s\leq t}$.
The following Assumption (ref) imposes restrictions on the $\alpha$-mixing coefficients that characterize the time dependence of the innovations $u_{t}$ and $\varepsilon_{t}$. For any two $\sigma$-fields $\mathcal{A}$ and $\mathcal{B}$, define $\alpha(\mathcal{A},\mathcal{B})=\sup_{A\in\mathcal{A},B\in\mathcal{B}}|\Pr\left(A\cap B\right)-\Pr(A)\Pr(B)|$ and $\alpha(d)=\sup_{s\in\mathbb{Z}}\alpha(\sigma(\{u_{t},\varepsilon_{t}\}_{t\leq s}),\sigma(\{u_{t},\varepsilon_{t}\}_{t\geq s+d}))$.
The following Assumption (ref) depicts the contemporary correlation of $v_{t}$ defined as ((ref)), as well as the constants in the local-to-unity AR coefficients specified in ((ref)). Define $\Omega=\Phi\Phi^{\top}$ where $\Phi$ has appeared in ((ref)).
We specify the user-determined parameter in ((ref)) as
with absolute constants $C_{\zeta}>0$ and $\tau\in(0,1)$. The choice of $\tau$ determines the persistence of the IV $\zeta_{j,t}$, which will be elaborated in Remark (ref). The following Assumption (ref) characterizes the number of regressors $p$ and the sparsity index $s$ relative to the sample size $n$.
In the following, we expound the differences between the assumptions above and those in MS24. First, Assumption (ref) imposes the m.d.s. and conditional homoskedasticity conditions for the error term $u_{t}$. Three remarks are in order to justify these two conditions.
Second, MS24 assume the innovations follow linear processes, and impose the mixing condition and the lower bounded long run variances through the coefficients in the linear process. In contrast, our Assumption (ref) does not assume any specific form of the linear process for $u_{t}$ and $\varepsilon_{t}$, but directly imposes the exponentially decaying rate for the mixing coefficient and the lower bound of the long run variances.
Third, Assumption (ref)(a) follows MS24 by allowing $p$ to diverge at a polynomial rate of $n$; it can be extended to an exponential rate of $n$ at the cost of expositional complications. Assumption (ref)(b) imposes a more restrictive condition for the sparsity index $s$, compared to $s=o(n^{1/4})$ in MS24. This is understandable as asymptotic normality is more delicate and demanding than consistency. This condition ensures that the shrinkage bias is accurately estimated to achieve the result (R2) in Section (ref), so that XDlasso asymptotically follows a normal distribution centered at the true coefficient; also see the second term of ((ref)) below.
Assumptions (ref)--(ref) will be sufficient for the consistency of Slasso with high dimensional LURs and stationary regressors. As the bias correction of XDlasso is mounted on the workhorse estimator Slasso, the consistency of the latter is a prerequisite for the ensuing maneuver.
The leading case of persistent regressors is LUR in the low-dimensional predictive regressions campbell2006efficient. LUR includes the unit root as a special case, and is thus more general in modeling nonstationary behaviors. lee2022lasso study the variable selection properties by the adaptive LASSO under a finite number of LUR regressors, and MS24 cover the consistency of Slasso under high dimensional unit roots. In this paper, the first theoretical result extends Slasso's consistency in the latter paper to incorporate the LUR processes in the former one. This generalization calls for sophisticated arguments, as briefed in Remark (ref) about the technical challenges of LURs.
The consistency of Slasso is founded on two building blocks. The first one, which is essential and challenging, is the restricted eigenvalue (RE) condition of the Gram matrix of the standardized regressors. For any $L>1$, the RE of any $p\times p$ matrix $\Sigma$ is defined as
where $\mathcal{R}(L,s)=\{\delta\in\mathbb{R}^{p}\backslash\{0_{p}\}:\|\delta_{\mathcal{M}^{c}}\|_{1}\leq L\|\delta_{\mathcal{M}}\|_{1},\ \text{for all }|\mathcal{M}|\leq s\}.$ The generic matrix $H$ is a placeholder and varies in different contexts. Let $\widehat{\Sigma}=W^{\top}W/n$ be the sample Gram matrix of all regressors. In the context of Slasso, we consider $\Sigma=\widehat{\Sigma}$ and $H=D$ along with the scale standardization in ((ref)). The choice of the constant $L$ is related to the procedures of technical proofs and does not impact the rate of convergence. Following the common practice bickel2009simultaneous, we set $L=3$ as a convenient choice, and simplify the notation as $\widehat{\kappa}_{D}=\kappa_{D}(\widehat{\Sigma},3,s)$. The quantity $\widehat{\kappa}_{D}$ will appear, according to Lemma 1 in MS24, in the denominator of Slasso's convergence rates. Therefore, a lower bound for $\widehat{\kappa}_{D}$ is essential for the consistency of Slasso.
The second condition for Slasso's consistency is the deviation bound (DB) of the cross-product between the error term $u_{t}$ in ((ref)) and the standardized regressors. The theoretical order of the tuning parameter $\lambda$ must be no smaller than that of $\|n^{-1}\sum_{t=1}^{n}D^{-1}W_{t-1}u_{t}\|_{\infty}$ to avoid overfitting. On the other hand, an excessively large $\lambda$ causes over shrinkage and damages consistency. A tight upper bound of $\|n^{-1}\sum_{t=1}^{n}D^{-1}W_{t-1}u_{t}\|_{\infty}$ is therefore indispensable.
Next, we establish the RE and DB conditions for high dimensional mixed roots, and highlight their similarities to and differences from those in MS24. We then leverage them to derive the convergence rates of Slasso.
Since LURs share similar asymptotic behavior with unit roots, the orders of the RE (Eq. ((ref))) and DB (Eq. ((ref))) are the same as those in Proposition 3 of MS24. Nevertheless, the technical proofs for LURs are challenging.
With RE and DB, we formulate the consistency of the initial workhorse Slasso estimator in the following lemma.
Lemma (ref) provides the consistency of the standardized coefficients, which is necessary to establish the asymptotic normality of XDlasso. The rates of tuning parameters are technical devices for proofs. In practice, we recommend cross-validation to select the tuning parameter $\lambda$.
While Lemma (ref) generalizes the convergence of Slasso from unit roots to LUR regressors, brand new theory must be developed for the auxiliary regression ((ref)), again by the workhorse Slasso when $w_{j,t}$ is either LUR or stationary. Importantly, as we maintain an agnostic attitude about the persistence of the variable of interest, which is the dependent variable of the auxiliary regression, the theory must be applicable in a unified manner to accommodate two very different types of time series. Without loss of generality, let $\mathcal{M}_{x}=[p_{x}]$ and $\mathcal{M}_{z}=[p]\backslash[p_{x}]$ denote the integer sets indexing the locations of LURs and stationary regressors, respectively.
We first examine the case when $j\in\mathcal{M}_{x}.$ The standardized instrument $\widetilde{\zeta}_{j,t}$ has a different degree of persistence from $X_{-j,t}$ and $Z_{t}$, due to the different orders of integration. In the low dimensional framework, the FCLT in phillips2016robust yields the following rate of convergence in OLS:
In high dimensional models, the sample Gram matrix is rank-deficient and the FCLT no longer works. Thanks to the $L_{1}$ penalization, Slasso has a comparable local-to-zero order as ((ref)), which is shown in Proposition (ref) below.
We then turn to $j\in\mathcal{M}_{z}$, and slightly abuse the notation $Z_{-j,t}$ to denote the vector of stationary regressors excluding $w_{j,t}$. Define
as the linear projection of $w_{j,t}$ on $Z_{-j,t}$, and normalize it by the standard deviation of the IV as
The “pseudo-true” model for the LASSO regression ((ref)) is
where $\text{\ensuremath{\widetilde{\eta}}}_{j,t}=(\zeta_{j,t}-Z_{-j,t}^{\top}\varphi_{0z}^{(j)*})/\widehat{\varsigma}_{j}$. Note that when $w_{j,t}$ is stationary, the IV $\zeta_{j,t}$ is close to $w_{j,t}$ under a large sample size, and thus asymptotically uncorrelated with the LUR regressors $X_{t}$. As a result, the coefficients associated with $X_{t}$ in the pseudo-true model ((ref)) are zero.
In addition, the error term $\widetilde{\eta}_{j,t}$ is close to the stationary time series $(w_{j,t}-Z_{-j,t}^{\top}\varphi_{0z}^{(j)*})/\widehat{\varsigma}_{j}$ and thus asymptotically uncorrelated to the nonstationary regressors $X_{t}$ due to different persistence. Furthermore, the coefficient $\varphi_{0z}^{(j)*}$ satisfies $\mathbb{E}\left(Z_{-j,t}(w_{j,t}-Z_{-j,t}^{\top}\varphi_{0z}^{(j)*})\right)=0$. Therefore, $\widetilde{\eta}_{j,t}$ is also asymptotically uncorrelated to the stationary regressors $Z_{-j,t}$, thereby ensuring the consistency of the Slasso estimator $\widehat{\varphi}^{(j)}$ in ((ref)). We impose the following assumption on the coefficient $\varphi_{0z}^{(j)*}$.
To bound the LASSO estimation errors, we need sparsity of not only the main regression ((ref)), but also the auxiliary regression ((ref)). This is similar to the high dimensional sparse instrumental variable regression; see zhu2018sparse, gold2020inference, and fan2025heteroscedasticity. In Assumption (ref), we slightly abuse the sparsity index $s$ to bound the number of nonzero coefficients in the vector $\varphi_{0z}^{(j)*}$. For simplicity, we directly impose this high-level sparsity assumption on $\varphi_{0z}^{(j)*}$, which can be deduced under the commonly used sparsity restriction on the precision matrix $\left(\mathbb{E}\left(Z_{t}Z_{t}^{\top}\right)\right)^{-1}$; see the definition of $s_{j}$ and the conditions in Theorem 2.1 of zhang2017simultaneous. Finally, the upper bound of the $L_{1}$-norm of $\varphi_{0z}^{(j)*}$ controls the variance of the error term in the pseudo-true model ((ref)). Parallel conditions naturally hold for $j\in\mathcal{M}_{x}$ as the pseudo-true coefficients are zero in view of the local-to-zero OLS estimate displayed in ((ref)).
The following proposition formally lays out the convergence rate of the auxiliary estimator ((ref)).
The consistency in the main equation and the auxiliary regression has paved the way for statistical inference. We will analyze the asymptotic size and power of the concerning test statistics in the next section.
The desirable rate of the auxiliary regression built in Proposition (ref) guarantees the following decomposition of the $t$-statistic:
where ${\rm sgn}_{j}=\frac{\left|\sum_{t=1}^{n}\widehat{r}_{j,t-1}w_{j,t-1}\right|}{\sum_{t=1}^{n}\widehat{r}_{j,t-1}w_{j,t-1}}$ is either 1 or $-1$ with probability one as $\sum_{t=1}^{n}\widehat{r}_{j,t-1}w_{j,t-1}$ is continuously distributed. The first term of the above expression is a counterpart of $\mathrm{N}_{j}/\omega_{j}$ in the discussion of generic desparsifying argument in ((ref)). The second term is the convergence rate of the approximation error of the shrinkage bias analogous to $(\widehat{{\rm B}}_{j}-{\rm B}_{j})/\omega_{j}$. Under Assumption (ref), the second term on the right-hand side of ((ref)) is asymptotically negligible, thereby yielding the following asymptotic normality for any $j\in[p]$.
The asymptotic normality in Theorem (ref) is our main theoretical result that delivers the valid asymptotic size of the hypothesis testing for $\mathbb{H}_{0}:\theta_{j}^{*}=\theta_{0,j}$ using the $t$-statistic $t_{j}^{{\rm XD}}$ in ((ref)). The following Theorem (ref) provides the convergence rate of the estimated standard error, which characterizes the asymptotic power of the test.
The convergence rates displayed in Theorem (ref) for the two types of regressors are coherent with the results in low dimensional predictive regressions. When $j\in\mathcal{M}_{x}$, the standard error converges faster than the rate $1/n^{\delta_{j}}$ for any $\delta_{j}\in(0,(1+\tau)/2)$. Thus, for the null hypothesis $\mathbb{H}_{0}:\beta_{j}^{*}=0$, the hypothesis testing based on the $t$-statistic $t_{j}^{{\rm XD}}$ is consistent under the alternative with $\beta_{j}^{*}=c/n^{\delta_{j}}$ for some $c>0$ over this class of $\delta_{j}$. The range $\delta_{j}\in(0,(1+\tau)/2)$ includes the important $1/\sqrt{n}$ rate of coefficients for the LUR regressor $X_{j,t}$, under which ${\rm var}(X_{j,t}\beta_{j}^{*})=O(1).$ The $1/\sqrt{n}$ factor thus balances the larger order of LUR regressors $X_{t}$ and the standard $O_{p}(1)$ stochastic order of $y_{t}$ phillips2015halbert,lee2022lasso, and the test by XDlasso is consistent under this class of alternatives. When $j\in\mathcal{M}_{z}$, XDlasso achieves the standard $\sqrt{n}$-consistency.
Finally, when we are interested in a joint null hypothesis $\mathbb{H}_{0}:\theta_{\mathcal{J}}^{*}=\theta_{0,\mathcal{J}}$ for the coefficients indexed by $\mathcal{J}\subset[p]$ with a fixed cardinality, we can use the Wald statistic in ((ref)) for this test.
Theorem (ref) is a natural extension from the $t$-statistic for a one-dimensional univariate hypothesis to a simultaneous multivariate hypothesis. Notice that the Wald statistic is valid even if $\mathcal{J}$ includes a mixture of LUR regressors and stationary ones. It allows the researcher to maintain the agnostic attitude when conducting the joint hypothesis testing.
Let us summarize the insights gained from the theoretical development. The fundamental principle behind the Dlasso method, as discussed in zhang2014confidence, is based on the Frisch-Waugh-Lovell theorem. This theorem purges the influence of other control variables to overcome shrinkage bias. In the context of persistent predictors, the two-stage least squares approach developed in magdalinos2009limit is used to eliminate Stambaugh bias. Each piece of XDlasso is a machine learning version of a classical idea. To adapt these existing procedures to high dimensional predictive regressions, we must rely on the consistency of Slasso, with its extension to LURs. This consistency is crucial for both the main regression and the auxiliary regression. Moreover, as mentioned in Remark (ref), the construction of the IV and the choice of $\tau$ must be tailored to balance the two types of predictors. This approach is original in that such a restriction is unique to high dimensional models and has not been studied before even in conventional low dimensional settings.
In this section, we evaluate the performance of the proposed XDlasso inference procedure by comparing its test size and power with those of Dlasso. Although Dlasso and other existing LASSO bias-correction procedures are not designed to handle persistent regressors, we include this comparison to highlight the value added by the IVX transformation. Furthermore, to demonstrate the robustness of the XDlasso approach, we benchmark it using tests based on infeasible estimators, where an oracle reveals the locations of the nonzero coefficients.
We consider the linear predictive regression model in ((ref)). The vector of the stationary components, denoted as $v_{t}=\left(u_{t},e_{t}^{\top},Z_{t}^{\top}\right)^{\top}$, is generated by the two different processes:
where $R_{n}=\text{diag}(0,0.3,0.3,\ldots,0.3)$. Under this $R_{n}$ the error term $u_{t}$ in the main regression ((ref)) remains i.i.d., satisfying the m.d.s. condition in Assumption (ref), while the (local) unit root innovations $e_{t}$ and the stationary regressors $Z_{t}$ are AR(1) processes. The covariance matrix $\Sigma=(\Sigma_{ij})_{i,j\in[p+1]}$ is specified as \[ \Sigma_{ij}=
\] The persistent regressors $X_{t}$ are generated by
where $\rho^{*}=(1,1-1/n,1+1/n,1,1-1/n,1+1/n,\cdots)^{\top}\in\mathbb{R}^{p_{x}}$. Recall $0_{p}$ is a $p$-dimensional zero vector and $1_{p}$ is a $p$-dimensional vector of ones. The true coefficient vectors are:
The specification involves four active LUR regressors and four active stationary regressors.\footnote{In Section B.1, we consider a setting with more nonzero coefficients and find robust performance of XDlasso.} The $1/\sqrt{n}$ scaling balances the regression by normalizing the coefficients of the LUR regressors. We test the hypothesis $\mathbb{H}_{0}:\beta_{1}^{*}=0$, $\mathbb{H}_{0}:\gamma_{1}^{*}=0$, and the joint null hypothesis $\mathbb{H}_{0}:\beta_{1}^{*}=\gamma_{1}^{*}=0$, respectively, under the sample sizes $n\in\{200,300,400,500,600\}$ and the dimensionality pairs $(p_{x},p_{z})\in\{(50,100),(100,150),(150,300)\}.$ We conduct 2000 replications in each setting.
We compare the finite sample performance of XDlasso, as described in Algorithm (ref), and Dlasso with the score vector in (ref) and (ref). In addition, we consider two infeasible testing procedures as benchmarks. Using the known active set of regressors, we conduct IVX inference (IVX oracle) and the standard $t$-test based on the OLS estimator (OLS oracle), employing only the regressor of interest and the active regressors, which form a low-dimensional predictive regression model. We set $C_{\zeta}=5$ and $\tau=0.5$ for the parameter $\rho_{\zeta}$ specified as ((ref)). As discussed in Remark (ref), the choice $\tau=$ 0.5 admits the weakest sparsity condition, and thus effectively improves the finite sample performance.
Both XDlasso and Dlasso involve Slasso, where the selection of tuning parameters $\lambda$ and $\mu$ affects finite sample performance. In our experiments, we employ the block 10-fold cross-validation (CV), splitting the sample into 10 equally sized chronologically ordered consecutive blocks for validations. Though the unconditional variances of nonstationary regressors vary in different chronological blocks, the standardization in the Slasso estimators like ((ref)) and ((ref)) account for such variation. This explains the robustness of the block CV to nonstationary time series, as shown by our simulation results.
In Theorem (ref), the tuning parameters are specified as constants multiplied by the appropriate rates of convergence determined by the sample size $n$, dimensionality $p$ and the mixing condition constant $r$. As a benchmark, we also calibrate the tuning parameters following lee2022lasso to examine the validity of the theoretical orders of tuning parameters specified in Theorem (ref). Specifically, we perform $500$ pilot replications for each DGP, with $n_{0}=400$ , $\left(p_{x0},p_{z0}\right)=\left(100,150\right)$, and $p_{0}=p_{x0}+p_{z0}$. In each replication $q=1,2,\dots,500$, we use the 10-fold cross-validation to choose the tuning parameters $\lambda^{(q)}$ and $\mu^{(q)}$, and calibrate the constants as \[ C_{{\rm m}}^{\left(q\right)}=\lambda^{(q)}n_{0}^{1/2}/(\log p_{0})^{\frac{3}{2}+\frac{1}{2r}},\ \ C_{{\rm a}}^{\left(q\right)}=\mu^{(q)}n_{0}^{[(1-\tau)\wedge\tau]/2}/(\log p_{0})^{2+\frac{1}{2r}}, \] where $r=1$ and $\tau=0.5$ are chosen in the simulation. We then fix $\widehat{C}_{\star}=\mathrm{median}(C_{\star}^{\left(1\right)},\ldots,C_{\star}^{\left(500\right)})$ for $\star\in\left\{ {\rm m},{\rm a}\right\} $ in the full-scale experiments. The tuning parameters are then set as $\widehat{\lambda}=\widehat{C}_{{\rm m}}(\log p)^{\frac{3}{2}+\frac{1}{2r}}/\sqrt{n}$ and $\widehat{\mu}=\widehat{C}_{{\rm a}}(\log p)^{2+\frac{1}{2r}}/\sqrt{n^{(1-\tau)\wedge\tau}}$ as in Theorem (ref).
We first investigate the empirical size of different testing methods at a 5% nominal significance level when the true coefficients $\beta_{1}^{*}=0$ and $\gamma_{1}^{*}=0$. Tables (ref) and (ref) report the empirical sizes under the IID and AR(1) innovations, respectively. Foremost among the findings is that XDlasso effectively controls the empirical size for both $\beta_{1}^{*}$, associated with a unit root regressor, and $\gamma_{1}^{*}$, associated with a stationary regressor. This performance stands in sharp contrast to that of Dlasso and OLS oracle, which exhibit severe size distortions for $\beta_{1}^{*}$. Such distortions can be attributed to the failure to account for the Stambaugh bias arising from nonstationarity. Furthermore, the results yield noteworthy insights regarding the tuning parameter selection. The empirical size of XDlasso with both cross-validated and calibrated tuning parameters is close to the nominal level. This result not only validates the asymptotic rates of tuning parameters specified in Theorem (ref) but also supports CV as a feasible data-driven tuning parameter selection method in practice. In addition, we investigate the efficiency and robustness of XDlasso in comparison with alternatives. When compared to the unbiased but infeasible “IVX oracle” estimator, the confidence intervals produced by XDlasso are only slightly wider. In contrast to the estimators solely for low-dimensional data, XDlasso demonstrates robustness by accommodating high dimensional covariates without compromising much in efficiency. Lastly, the empirical sizes for testing the joint null hypothesis $\mathbb{H}_{0}:\beta_{1}^{*}=\gamma_{1}^{*}=0$, as reported in Table (ref), are also well controlled around the nominal level across setups. The results validate the theoretical result in Theorem (ref) and is consistent with the findings for testing $\mathbb{H}_{0}:\beta_{1}^{*}=0$ and $\mathbb{H}_{0}:\gamma_{1}^{*}=0$ individually.
We now turn to the empirical power of the XDlasso inference. Figure (ref) plots power curves for the null hypotheses $\mathbb{H}_{0}:\beta_{1}^{*}=0$ and $\mathbb{H}_{0}:\gamma_{1}^{*}=0$ under varying true coefficient values. In our analysis, we vary either $\beta_{1}^{*}$ or $\gamma_{1}^{*}$ from $0$ to $0.5$. All remaining coefficients are held fixed as specified in (ref). Across various configurations, XDlasso exhibits increasingly high power against the null hypothesis as the true coefficient moves away from 0 and the sample size $n$ increases. Furthermore, Figure (ref) compares the power for testing $\beta_{1}^{*}$ and $\gamma_{1}^{*}$, which reveals that the power associated with a unit root regressor surpasses that of a stationary regressor. This observation provides finite sample evidence to support the theoretical results in Theorem (ref). Specifically, it corroborates the faster convergence of standard error for unit root regressors compared with stationary regressors, thereby inducing higher power in the hypothesis testing.
In recent years, high dimensional macroeconomic data have been extensively used to forecast key variables of interest; see smeekes2018macroeconomic, medeiros2021forecasting, and giannone2021economic, for example. Researchers have primarily focused on the point estimation of forecast. There has been limited empirical exploration of statistical inference on the predictive power of specific predictors, due to a lack of suitable toolkits. This section showcases two empirical applications using our proposed XDlasso inference method. We utilize the monthly data of 112 U.S. macroeconomic variables spanning from January 1960 to April 2025, sourced from the FRED-MD dataset by mccracken2016fred.
In financial economics, testing the predictive power of valuation ratios, particularly the earnings-price ratio, on stock returns has been subject to widespread discussions. Much of the literature tests the predictability using univariate predictive regressions (e.g., welch2008comprehensive, zhu2014predictive, and goyal2024comprehensive), but inference can be sensitive to model misspecification arising from omitted variables. Controlling for high dimensional covariates is therefore necessary not only to enhance out-of-sample prediction but also to mitigate the omitted variables for credible and accurate inference. The literature on predictive regression for stock returns has focused on identifying which variables possess significant predictive power for future returns. For better out-of-sample prediction, recent literature has documented gains in forecasting performance from incorporating high-dimensional covariates into predictive regressions; see, for instance, smeekes2018macroeconomic, gu2020empirical and medeiros2021forecasting. The increasing popularity of statistical learning methods with high dimensional predictors calls for suitable inference methods beyond univariate predictive regressions. The goal is to provide investors with signals to adjust their portfolios based on changes in these predictive variables. In this section, we investigate the predictability of stock returns using the log earnings-price ratio in a data-rich environment with high dimensional mixed-root control variables.
Our analysis focuses on predicting the monthly return of the S&P 500 index, calculated as $\text{Return}_{t}=\log(P_{t})-\log(P_{t-1}),$ where $P_{t}$ refers to S&P 500 (S&P's Common Stock Price Index: Composite). The primary predictor of interest is the log earnings-price ratio. We obtain it from the variable S&P PE ratio (S&P's Composite Common Stock: Price-to-Earnings Ratio), denoted as $\text{PE}_{t}$. The log earnings-price ratio is calculated by inverting the original price-earnings ratio as $\text{logEP}_{t}=\log(1/\text{PE}_{t}).$
Figure (ref) displays the time series plot of the monthly return of the S&P 500 index and the log earnings-price ratio from January 1960 to April 2025. The log earnings-price ratio exhibits persistent patterns, while the S&P 500 return appears stationary. We further report the AR(1) coefficient estimates and Augmented Dickey--Fuller (ADF) test $p$-values of both series under different sample periods in Table (ref). The S&P 500 return is evidenced to be stationary under the full sample with an ADF test $p$-value below 1%, rejecting the null hypothesis of nonstationarity. Conversely, the log earnings-price ratio shows high persistence with an AR(1) coefficient estimate equal to $0.993$. Given the $p$-value of $0.074$, nonstationarity is not rejected at the 5% significance level. Note that the nonstationary log earnings-price ratio can predict the stationary monthly return of S&P 500 in our model, since we allow for a local-to-zero coefficient to balance the different scales between a stationary outcome and a nonstationary regressor phillips2015halbert. Theorem (ref) and the paragraph that follows illuminate that inference by XDlasso has the power to detect a wide range of local-to-zero alternatives.
Return predictability can be time-varying and sporadic tu2023penetrating. campbell2006efficient discover that the return predictability test results can vary depending on the inclusion of the samples after 1994, and the break is revisited by zhu2014predictive. Following this empirical finding, we divide the sample into two periods: pre-1994 (January 1960 - December 1993) and post-1994 (January 1994 - April 2025). The S&P 500 return demonstrates stationarity across both subperiods with ADF $p$-values below 1%. In contrast, the persistence of the log earnings-price ratio depends on the sample period considered. Nonstationarity is evident in the Pre-1994 period with a large ADF test $p$-value $0.776$, while the null hypothesis of nonstationarity is rejected at the 1% significance level in the Post-1994 period. The ambiguity in stationarity of the log earnings-price ratio motivates the use of the XDlasso procedure, a unified approach for both stationary and nonstationary regressors without prior knowledge of their persistence.
Besides the monthly return of the S&P 500 index and the log earnings-price ratio, our analysis incorporates high dimensional covariates including all other 110 macroeconomic variables from the FRED-MD dataset. These variables comprise a mixture of stationary and nonstationary time series. In practice, it is common for empirical analysts to transform potentially nonstationary time series into stationary ones to avoid challenges arising from nonstationarity. To facilitate such stationarization, FRED assigns a transformation code (TCODE) to each variable denoting the recommended transformation.
In our analysis, we follow the common practice of using the TCODE to transform the 110 time series, and subsequently use the transformed variables as covariates. Nevertheless, we highlight that these elementary transformations are not a silver bullet in taming nonstationarity. We perform the ADF test for each of the transformed variables for different sample periods. A nontrivial proportion (9.1%) of the transformed variables still demonstrate nonstationarity based on the ADF test at the 5% significance level for both Pre-1994 and Post-1994 subperiods. The high persistence of the log earnings-price ratio and other covariates suggests the necessity of XDlasso.
As highlighted by smeekes2020unit, the predictive performance in regressions using FRED-MD data is sensitive to the transformations. To assess the robustness of our results, we also conduct our analysis using the original (untransformed) FRED-MD time series as covariates. To reduce the impact of highly nonstationary series, we exclude I(2) variables that require second differencing for stationarity according to their TCODE classification.\footnote{Appendix C.1 further investigates robustness by considering: (i) excluding nonstationary variables based on integrated orders determined by the bootstrap sequential testing procedure of smeekes2015bootstrap, as reported in smeekes2020unit, and (ii) applying only logarithmic transformations, as indicated by TCODE, without differencing.}
We study the one-month-ahead regression $\text{Return}_{t}=\alpha^{*}+\theta_{1}^{*}\times\text{logEP}_{t-1}+W_{-1,t-1}^{\top}\theta_{-1}^{^{*}}+u_{t},$ where $W_{-1,t-1}$ denotes a high dimensional vector that collects all control variables. We carry out hypothesis testing on the key parameter of interest $\theta_{1}^{*}$ that measures the predictive power of the log earnings-price ratio in forecasting the S&P 500 monthly return. Given the stationary pattern of the dependent variable, we additionally consider specifications where $W_{-1,t-1}$ includes the lagged dependent variable $\text{Return}_{t-1}$. For each model, we apply the wild bootstrapped automatic variance ratio test kim2009automatic to the Slasso residuals as a diagnostic test for the martingale difference sequence (m.d.s.) condition in Assumption (ref).\footnote{This practice serves as a heuristic diagnostic. Demonstrated by our simulation results in Appendix B.4, the variance ratio test on the Slasso residuals $\widehat{u}_{t}$ tends to severely over-reject the m.d.s. condition for $u_{t}$. To the best of our knowledge, at present the literature has no valid testing procedure yet for m.d.s. in high dimensional predictive regressions. It is an open question for future research.}
In addition to XDlasso and Dlasso, we conduct the IVX inference in a simple regression setup, omitting all control variables, to demonstrate the necessity of incorporating high dimensional control variables in practice. We compare the results of XDlasso to those of Dlasso for the high dimensional regression and IVX for a simple regression $\text{Return}_{t}=\alpha^{*}+\theta_{1}^{*}\times\text{logEP}_{t-1}+u_{t}.$
Table (ref) presents the point estimates and testing results for $\theta_{1}^{*}$. IVX delivers negative point estimates across all sample periods and detects significant effects in the full sample period and Pre-1994 subperiod at the 5% significance level. However, the resulting negative relation between stock return and earning performance contradicts the economic mechanism between the two variables. A large number of potential confounding variables are present in the dataset, which may be the culprit in producing the counterintuitive result from the simple regression.
Dlasso, which accounts for high dimensionality but not nonstationarity, generally reverses the sign of the estimates compared to the simple regression. For the full sample, Dlasso consistently yields positive but statistically insignificant coefficients. Zooming into the Pre-1994 period, Dlasso detects strong evidence of predictability with transformed data, with coefficients significant at the 5% level. This result persists with the untransformed data, though the significance weakens to the 10% level when the lagged dependent variable is included. Accounting for the potential nonstationarity in the predictors, XDlasso provides notably different estimates and standard errors from Dlasso. XDlasso reveals no significant predictive power of log earnings-price ratio for stock return in the full sample and the Post-1994 period. Moreover, for the Pre-1994 period, the evidence of predictability is weakened to the 10% significance level with transformed data, and becomes entirely insignificant with untransformed data. Our empirical findings with XDlasso are consistent with the general recognition that there is little predictability in the financial market, thereby supporting the efficient market hypothesis. The divergence between Dlasso and XDlasso underscores our theoretical prediction: XDlasso mitigates the Stambaugh bias arising from highly persistent regressors while Dlasso does not. Despite the short confidence intervals, Dlasso results can be misleading when nonstationary time series is present, as clearly shown by the illustrative simulation in Section (ref).
The variance ratio tests on the Slasso residuals provide a diagnostic check on the m.d.s. condition required by our asymptotic theory. The results exhibit stark differences across data specifications: with untransformed data and without the lagged dependent variable, the tests strongly reject the m.d.s. condition ($p$-values ranging from 0.001 to 0.046), suggesting potential model misspecification. Including the lagged return or using transformed data substantially improves the diagnostic results, with $p$-values generally exceeding 0.05. This pattern suggests that both data transformation and including dynamics help satisfy the underlying assumptions, thereby boosting credibility for the XDlasso results under these specifications.
It is essential for monetary policymakers to understand the relationship between unemployment and inflation. As engemann2020phillips pointed out, “The Federal Reserve has a dual mandate to promote maximum sustainable employment and price stability.” First alluded to by fisher1926statistical,fisher1973discovered though, phillips1958relation popularized the Phillips curve --- a negative relationship between the level of unemployment and the change rate of money wage rates. There has been a prolonged debate about whether the unemployment is a credible barometer for inflation among not only modern economic studies, but also policymakers.\footnote{Mary Daly, San Francisco Fed President, delivered at daly2019new a negative view on the Phillips curve that \textquotedblleft As for the Phillips curve\ldots most arguments today center around whether it\textquoteright s dead or just gravely ill. Either way, the relationship between unemployment and inflation has become very difficult to spot.\textquotedblright John Williams, New York Fed President, expressed a different opinion that \textquotedblleft The Phillips curve is the connective tissue between the Federal Reserve\textquoteright s dual mandate goals of maximum employment and price stability. Despite regular declarations of its demise, the Phillips curve has endured. It is useful, both as an empirical basis for forecasting and for monetary policy analysis.\textquotedblright See engemann2020phillips for more details.} Empirical findings suggest that inflation rate can be either positively or negatively correlated with unemployment, depending on the shocks to the economies, the policies, and the lag orders niskanen2002death,gordon2011history,gordon2013phillips. Given the ongoing debate, we revisit the Phillips curve in a predictive regression framework utilizing the FRED-MD dataset.
The inflation rate, as the outcome variable in the predictive regression, is calculated by $\text{Inflation}_{t}=\left(\log({\rm CPI}_{t})-\log({\rm CPI}_{t-1})\right)\times100,$ where ${\rm CPI}_{t}$ denotes the Consumer Price Index for All Urban Consumers: All Items (CPI). The unemployment rate, as the predictor of interest, denoted as $\text{Unrate}_{t}$, is retrieved in its original form under the name UNRATE.
Similar to our first empirical application in Section (ref), we follow benati2015long to delineate three subperiods in addition to the full sample: Pre-Volcker (January 1960 - July 1979), Volcker and Greenspan (August 1979 - January 2006), and Bernanke, Yellen, and Powell (February 2006 - April 2025). These subperiods correspond to different eras in U.S. monetary policy, each named after the Federal Reserve chairperson who presided during that time.\footnote{Pre-Volcker (before August 1979): This period was characterized by high and volatile inflation, with the Federal Reserve lacking a clear nominal anchor. Volcker and Greenspan (August 1979--January 2006): Volcker, who served as the Federal Reserve Chairman from August 1979 to August 1987, implemented aggressive anti-inflation measures, notably raising interest rates to historically high levels. Greenspan succeeded Volcker and continued to focus on maintaining price stability during his tenure, which contributed to a period of low and stable inflation known as the “Great Moderation”. Bernanke, Yellen, and Powell (February 2006--December 2019): Bernanke's tenure as Fed Chairman was marked by the Great Recession and the implementation of unconventional monetary policies, such as quantitative easing, aimed at stimulating the economy and preventing deflation. Yellen and Powell carried on these policies.} This periodization allows us to examine how the relationship between inflation and unemployment may have evolved across different policy regimes and economic conditions.
Figure (ref) plots the inflation and unemployment rates over our sample period. Visual inspection suggests that the unemployment rate is more persistent than the inflation rate. Table (ref) further reports the AR(1) coefficient estimates and ADF test $p$-values of the inflation and the unemployment rate under different sample periods. The inflation rate appears stationary in most periods, with AR(1) coefficient estimates ranging from 0.54 to 0.64. However, the inflation rate during the pre-Volcker period shows a slight upward trend, and nonstationarity is indicated by the ADF test with a $p$-value of 0.22. In contrast, the inflation rate is found to be stationary in the other two subperiods. The unemployment rate, on the other hand, appears highly persistent with AR(1) coefficient estimates close to 1. The ADF test rejects the null hypothesis of nonstationarity at a 10% significance level for both the full sample and the Volcker-Greenspan period. However, during the Pre-Volcker and Bernanke-Yellen-Powell periods, there is strong evidence of nonstationarity in the unemployment rate. The stationarity of both inflation and the unemployment rate in different periods is again ambiguous, which prompts the use of XDlasso from an agnostic perspective.
In addition to the unemployment rate, we incorporate all other 110 macroeconomic variables from the FRED-MD dataset as controls after the TCODE transformation. Still, a significant proportion, from 11% to 17% across subperiods, of the transformed variables exhibit nonstationarity according to the ADF test at the 5% significance level. Consistent with Section (ref), we also perform the analysis using the original (untransformed) FRED-MD time series as covariates. To further address potential concerns regarding the high correlation between control variables and the unemployment rate, we additionally present results excluding variables from the labor market group, as classified by FRED-MD.
We study, in our predictive regression framework, a one-month-ahead regression $\text{Inflation}_{t}=\alpha^{*}+\theta_{1}^{*}\times\text{Unrate}_{t-1}+W_{-1,t-1}^{\top}\theta_{-1}^{^{*}}+u_{t},$ where $W_{-1,t-1}$ denotes high-dimensional covariates. Given the ambiguous persistence of the inflation rate discussed in the previous section, including the lagged dependent variable in our model may introduce further technical complications. Therefore, we do not recommend its inclusion in this analysis, and will focus on the traditional predictive regression setting. We carry out hypothesis testing on the key parameter of interest $\theta_{1}^{*}$ that measures the predictive power of the unemployment rate in forecasting inflation. The predictive form is of particular interest among policymakers in leveraging the relationship as a practical tool.
Table (ref) reports the point estimates and standard errors for $\theta_{1}^{*}$ using IVX, Dlasso, and XDlasso across sample periods and specifications. In the benchmark setup with the transformed data, the diagnostic check rejects the m.d.s. condition except for the Pre-Volcker period. Nevertheless, Dlasso, which ignores nonstationarity, delivers significantly positive coefficients for both the full sample and the Pre-Volcker period, in striking contrast to XDlasso and IVX. Using untransformed data significantly alleviates concerns about model misspecification, with $p$-values of the diagnostic test greater than 5% in all cases. Both XDlasso and IVX consistently find no significant predictive power of unemployment for inflation across all specifications and time periods. The empirical findings add new insight to the recent debates on the Phillips curve and echo Mankiw2024's latest pessimistic remark: “The large confidence intervals for the natural rate, together with the apparent futility of this Holy Grail search, lead me to think that we should not expect much from the Phillips curve as a guide for forecasting inflation or for judging the stance of policy.” On the other hand, with untransformed data, Dlasso yields significantly negative coefficients at the 1% level for the full sample and at the 10% level in the Volcker-Greenspan period, without controlling for other labor market variables. The unstable results of Dlasso across setups highlight the necessity of accounting for nonstationarity in the inference.
This paper proposes XDlasso to overcome the difficulties in hypothesis testing for high dimensional predictive regressions with stationary and nonstationary regressors. XDlasso fuses the IVX technique from time series econometrics and the debiasing technique from the high dimensional statistics, thereby reducing the order of biases to make them readily correctable. We establish the asymptotic normality and convergence rate of XDlasso. The validity of our methods is further evidenced by simulation studies and empirical applications.
We thank the editor, the associate editor, and two referees for their constructive comments that led to substantial improvement of this paper. We thank Andrii Babii, Yongmiao Hong, Degui Li, Yiu Lim Lui, Kenwin Maung, Jun Yu, Wenyang Zhang, and Yichong Zhang for helpful comments and suggestions. Gao acknowledges the partial support from Southern Methodist University (SMU) University Research Council (URC) Grants. Mei acknowledges the partial financial support from the Start-up Research Grant (SRG) from University of Macau (Project No. SRG2025-00062-FBA). Shi acknowledges the partial financial support from the National Natural Science Foundation of China (Project No. 72425007, No. 72133002).
\onehalfspacing
\setcounter{footnote}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{equation}{0} \setcounter{thm}{0} \setcounter{lem}{0} \setcounter{rem}{0} \setcounter{prop}{0}