High Dimensional Generalised Penalised Least Squares

High Dimensional Generalised Penalised Least Squares

} \affil[1]{Essex Business School , University of Essex} } } \affil[2]{ King's Business School , King's College London}

Introduction

Research in high-dimensional statistics and econometrics has witnessed a surge, because the dimensionality of available datasets, models and associated parameter spaces has grown massively, in relation to the sample size. After the seminal work of tibshirani1996regression, the Lasso has become a focus of this continuously growing literature, since it conducts \textcolor{black}{simultaneously} model \textcolor{black}{estimation and} selection. More recent work establishes the asymptotic behaviour of the Lasso estimator as well as its model selection consistency, known as the oracle property, see, for example, meinshausen2009lasso, raskutti2010restricted, van2009conditions, van2008high and zhang2008sparsity.

Although the majority of the theoretical work is centred around the asymptotic behaviour of the Lasso and relevant estimators in generalised linear models, the main interest still lies in simple designs, such as i.i.d. or fixed covariates, while limited work has been done towards the direction of regularised linear models with time series components. Early work of wang2007regression suggests a linear model with autoregressive error terms and their resulting Lasso estimator satisfies a knight2000asymptotics-type asymptotic property, {where} the number of covariates cannot be larger than the sample size. hsu2008subset consider the Lasso estimator under a Vector Autoregressive (VAR) process and derive their asymptotic results following the same setting as the former. Various other papers derive asymptotic results ensuring effective model selection for regularised estimators: nardi2011autoregressive consider autoregressive structure on the covariates, basu2015regularized consider stochastic regressions and transition matrix estimation in VAR models, and kock2015oracle consider models with stationary VAR covariates. These papers are indicative of the work in the high-dimensional literature under more general conditions, but use restrictive assumptions in the error term in order to get the oracle property.

There is significant research in high-dimensional econometrics which allows for more relaxed assumptions: kock_2016 shows that the adaptive Lasso is oracle efficient in stationary and non-stationary autoregressions, masini2019regularized consider linear time-series models with non-Gaussian errors, and 10.1214/19-AOS1840 consider sparse Gaussian VAR models exploring the efficiency of $\beta$-mixing assumptions, used to bound the prediction and estimation error of sub-Weibull covariates. Another class of papers focuses entirely on high-dimensional financial econometrics: babii2020machine and babii2021machine consider prediction and now-casting with panel data and high-dimensional time series respectively, sampled at different frequencies using the sparse-group Lasso.

Although, the Lasso provides an efficient avenue to estimation and variable selection in high-dimensional datasets, one is unable to perform inference on the parameters estimated. leeb2005model have proven that ignoring the model selection step in the Lasso, leads to invalid \textcolor{black}{uniform} inference. Recent developments hinge on post-selection inference, see, for example, berk2013valid and taylor2015statistical, among others. The incentive is that post-selection methods provide valid inference on the non-zero coefficients obtained after model selection, which is typically carried out at a first step, using the Lasso. Although post-\textit{Lasso} inference guarantees valid confidence intervals, it is subject to the model selection made prior to that step, which can be {misleading} when the conditions for a "well-behaved" model are not met, e.g. \emph{i.i.d.} errors and/or covariates. The latter facilitates the necessity for \textit{post-Lasso} inference to allow for more relaxed assumptions, as well as uniformity in the limit theory for general penalised models. A way to relax the \textcolor{black}{\textit{i.i.d.}/Gaussianity} assumption has been proposed by tian2017asymptotics and tibshirani2018uniform, who consider a bootstrap approach for asymptotically valid testing. However, \emph{i.i.d.} conditions on the covariates and errors need to be {assumed} in order examine the large sample properties of the method.

Alternative approaches that allow inference on the true parameters without the limitation of a "well-behaved" model, have been developed. These are based on debiased or desparsified versions of the Lasso, see, for example, javanmard2014confidence, van2014asymptotically, zhang2014confidence and on certain assumptions on the covariance matrix, such as fixed design, {i.i.d}-ness or sub-Gaussianity. Further, extensions to a time series framework have been introduced in the literature, see, for example, wang2020lasso consider the covariates and error terms of a high-dimensional model to be temporally and cross-sectionally dependent and apply bootstrap methods towards estimation and inference of the true parameters. babii2020inference use debiased sparse-group Lasso for inference in a lower dimensional group of parameters. kock2016oracle studies high-dimensional correlated random effects panel data models, where they allow for correlation between time invariant covariates and unobserved heterogeneity, as under fixed effects. \textcolor{black}{Finally, adamek2020lasso extend the \emph{desparsified Lasso} to a time series setting, using near-epoch dependence assumptions, allowing for non-Gaussian, serially correlated and heteroscedastic processes. Notice though that throughout the literature, no link has been made on the use of Generalised Least Squares (\emph{GLS}) type of inference to account for non-spherical errors, which is the main focus of this paper. }

Inference on Ridge regression instead of the Lasso has been considered by zhang2020ridge, where they propose a wild bootstrap algorithm to construct confidence regions and perform hypothesis testing for a linear combination of parameters. A similar approach has been considered by zhang2021debiased, allowing for inference under non-stationary and heteroscedastic errors.

In this paper, we contribute to the ongoing literature of high-dimensional inference under general conditions. We propose a penalised GLS-type estimator which utilizes estimated autocovariances of the residual in a linear regression, where the residual is allowed to follow a general autoregressive process. Lasso works as a preliminary estimator and is shown to be asymptotically consistent under mild assumptions on the covariates and error processes. We perform uniform inference via the debiased GLS Lasso, imposing mild restrictions on the autocorrelation of the error term. In addition, we relax assumptions on the error and most importantly the covariates processes, commonly used in the existing literature, for example fixed design, $i.i.d.$, sub-Gaussianity, by allowing them to be stationary $\alpha$-mixing processes.

The remainder of this paper is organised as follows. Sections (ref) -- (ref) present the model, proposed methodology and theoretical results. Section (ref) describes the regularisation parameter tuning, while Section (ref) contains the Simulation Study. Section (ref) concludes. Proofs and additional simulation results are relegated to the Supplementary Material.

Setup and notation

For any vector $\boldsymbol{x}\in \mathbb{R}^{n}, $ we denote the $\ell_{p}$-, $\ell_{\infty}$- and $\ell_0$- norms, as $\left\lVert \boldsymbol{x} \right\rVert_{p}=\left(\sum_{i=1}^{n} | x_{i} |^p\right)^{1/p}$, $\lVert \boldsymbol{x} \rVert_{\infty} = \max_{i=1,\ldots,n}|x_{i}|$, $\lVert \boldsymbol{x}\rVert_{0}= \sum_{i=1}^{n}\boldsymbol{\mathrm{1}}_{(x_{i}\neq 0)} = \rm{supp}(\boldsymbol{\beta})$, respectively, and $\rm{sign}$$(\boldsymbol{x})$ denotes the sign function applied element-wise on $\boldsymbol{x}$. We use "$\to_{P}$" to denote convergence in probability. For two deterministic sequences $a_n$ and $b_n$ we define asymptotic proportionality, "$\asymp$", by writing $a_n \asymp b_n$ if there exist constants $0 < a_1 \leq a_2$ such that $a_1b_n \leq a_n \leq a_2b_n$ for all $n \geq 1$. For any set $A$, $ |A|$ denotes its cardinality, while $A^c$ denotes its complement. For a real number $a$, $\lfloor a \rfloor$ denotes the largest integer no greater than $a$.

Theoretical considerations

We consider the following high-dimensional linear regression model,

where $\{\boldsymbol{x}_t\} = \left \{ \left(x_{t,1},\ldots,x_{t,p} \right)' \right \}$ is a $p$-dimensional vector-valued stationary time series, $ \boldsymbol{\beta}$ is a $p\times 1$ vector of unknown parameters and $u_{t}$ is a general second order stationary process. We make the following assumption on the processes $\left \{u_{t} \right \}$ and $\left \{ \boldsymbol{x}_t \right \}$:

Under Assumption (ref), $\{u_t\}$ admits an $AR\left(q\right)$ representation, which is a more general setting compared to the standard one, e.g. $u_{t}\sim i.i.d. $ We propose to estimate the unknown parameters associated with (ref) feasibly within the scope of a high-dimensional model. More specifically $ \left \{u_{t} \right \}$ is a second-order stationary process satisfying

Notice that a stationary $AR(p)$ process with i.i.d. increments such as (ref), in which the marginal distribution of $\{ \varepsilon_{t} \}$ has unbounded support satisfies the $\alpha$-mixing assumption. To further specify the conditions under which $\{ u_{t} \}$ is strong mixing, one may refer to theorem 12.4 of davidson1994stochastic, describing random sequences, which by definition are strong mixing and further on Chapter 14 of davidson1994stochastic.

To describe the method in detail, let \[\widetilde{u }_t = y_t - \boldsymbol{x}_t'\boldsymbol{\widetilde{\beta}} , \quad t=1,\ldots, T,\] where $\boldsymbol{ \widetilde{\beta}}$ is a preliminary Lasso estimate of $\boldsymbol{\beta}$. Further, let $\boldsymbol{\widehat{{\phi}}} = (\widehat{\phi}_1,\ldots, \widehat{\phi}_q) '$ be the $q^{\text{th}}$ order OLS estimator of the autoregressive parameters for $\{\widetilde{u}_{t} \}$ obtained as the the solution to the minimisation of \[ T^{-1}\sum_{t=q+1}^T \left(\widetilde{u}_t - \phi_{1} \widetilde{u}_{t-1} - \cdots - \phi_{q} \widetilde{u}_{t-q}\right)^2, \vspace{-1mm} \] over $\boldsymbol{\phi} = (\phi_1, \ldots, \phi_q)'$. Then, an asymptotically valid feasible and penalised GLS estimate of $\boldsymbol{\beta}$ can be obtained as the solution to the following problem

where $ \widetilde{\boldsymbol{X}}= \boldsymbol{\widehat{L} X} $, $ \widetilde{\boldsymbol{y}}= \boldsymbol{\widehat{L} y} $, $\boldsymbol{y}=\left( y_{1}, \ldots, y_{T} \right)'$ and $\widehat{\boldsymbol{L}}$ is a $\left(T-q \right)\times T$ matrix defined as:

One can use the scalar representation of $\widetilde{\boldsymbol{y}}$ and vector representation of $\widetilde{\boldsymbol{X}}$, to obtain

The loss function in (ref) corresponds to the $\ell_{1}$-penalised loss function using the estimates of $\boldsymbol{\widehat{\phi}}$, where the additional penalty is added to the least squares objective. Note that the asymptotic consistency of $\boldsymbol{\widehat{\beta}}$ can be established even prior to the estimation of $\boldsymbol{{\phi}}$. We prove this result in Lemma (ref). \textcolor{black}{Further, following the same argument that we used to define $\widehat{\boldsymbol{ L }}$ in (ref), we define ${\boldsymbol{L}}$, a $\left(T-q \right)\times T$ matrix, with $\widehat{\phi}_{j}$ replaced by $\phi_{j}$, denoting the matrix of $\boldsymbol{\phi}$ used to derive the infeasible GLS-Lasso estimates, when the degree of autocorrelation in $u_{t}$ is known. }

We highlight now the Assumptions on the covariates and errors, necessary to derive our theoretical results.

Our objective is to incorporate a GLS estimator to address the issue of serial-autocorrelation in the error term, $u_{t}$, that will enable sharper inference following the paradigm of the debiased Lasso. Such methods have been examined in the past in low dimensional cases, where $T>>p$, see, for example, amemiya1973generalized, and kapetanios2016semiparametric. We extend this framework to allow for $p>>T$, ensuring that $\widetilde{u}_t$ is asymptotically consistent using a penalised estimate of $\boldsymbol{\beta}$, while we limit our research to finite autoregressions. \textcolor{black}{ Note that theoretical results can be established even when $\{u_{t}\}$ admits an $AR(\infty)$ representation, however certain additional assumptions must be made about $\{u_t\}$, and the theoretical proofs would come at an increased mathematical cost. }

The use of penalised models highlights the necessity of imposing sparsity conditions on the parameter set, $\boldsymbol{\beta}$, which in turn allow for a degree of misspecification in the model. The latter can be approximated by a sparse linear model following certain boundedness conditions on the smallest eigenvalues of the sample variance-covariance matrix, $\left(\boldsymbol{x}_t\boldsymbol{x}_t'\right)/T$.

For some $\lambda\geq 0 $, we define the following index set

with cardinality $s_{0}=\left| S_{0} \right|$. Under an appropriate choice of the regularisation parameter, $\lambda$, $S_{0}$ contains all $ \beta_i$ that are "sufficiently large", while for $\lambda=0$, $s_0$ is $\lVert \boldsymbol{\beta}\rVert_{0}$. We make use of this set in the following assumption, which forms the compatibility condition, as seen in bickel2009simultaneous, raskutti2010restricted and Chapter 6 of buhlmann2011statistics. We make the following assumption under $|S_{0}|$.

Further, in Lemma (ref), we show that the oracle inequalities hold prior to estimating $\boldsymbol{\widehat{\phi}}$. Notice that in (ref) we use the population covariance matrix instead of the sample, $T^{-1}{\boldsymbol{x}_t\boldsymbol{x}_t'}$, for which an argument similar to (ref) holds and is formalised in the Supplementary Material, see Lemma 10.7. Notice that (ref) is the restricted eigenvalue condition on the population covariance, $E(\boldsymbol{x}_{t}\boldsymbol{x}_{t}')$, while $\phi_{0}^{-2}$ is the compatibility constant. Lemma (ref) is of paramount importance and is required in order to establish the following corollary and Theorem (ref).

A feasible penalised GLS

We start with the main theorem that provides non- asymptotic guarantees for the estimation and prediction errors of the Lasso under a modified compatibility condition, i.e. Assumption (ref), which restricts the smallest eigenvalue of the variance-covariance matrix. \textcolor{black}{The results of the next theorem are based on the deviation inequality, which illustrates that as long as the co-ordinates of $\boldsymbol{X'{ \boldsymbol{L}}'}\boldsymbol{Lu}/T$ are uniformly concentrated around zero and the quantity $\lVert \boldsymbol{X'{ \boldsymbol{L}}'}\boldsymbol{Lu}/T \rVert_{\infty}$ is sharply bounded, using $\boldsymbol{\phi}$ instead of $\boldsymbol{\widehat{\phi}}$, (ref) and (ref) hold w.p.a. one.}

The following Corollary serves as a consequence of Theorem (ref), providing asymptotic rates for the bounded processes in (ref) and (ref).

{

Point-wise valid inference based on the GLS Lasso

A natural avenue to inference, having obtained the estimates of $\boldsymbol{\beta}$ from (ref) is to re-estimate the parameters of the estimated active set, $\widehat{S}_{0}= \{ i : \; \widehat\beta_{i}\neq 0 \}$, via OLS. Formally, let $\boldsymbol{\ddot{\beta}}_{S_0}$ be the vector whose $i^{th}$ element equals the least square re-estimate for all $i \in \widehat{S}_{0}$ and zero otherwise, while $\boldsymbol{\widehat{\beta}}_{S_{0}}$ denotes the oracle assisted least squares estimates only including the relevant variables, those indexed by ${S}_{0}$. The following theorem shows that this indeed leads to point-wise valid confidence bands for the non-zero entries of $\boldsymbol{\beta}$, i.e. those indexed by $S_{0} $, since $\widehat{S} _{0}=S_{0}$ asymptotically.

Uniformly valid inference based on the GLS Lasso

We define the following notation: Let $$\boldsymbol{\Sigma}_{xu} = E[(\boldsymbol{{X}}'\boldsymbol{L}'\boldsymbol{Lu})(\boldsymbol{{X}}'\boldsymbol{L}'\boldsymbol{Lu})'], \quad \boldsymbol{\widehat{\Sigma}}_{xu} =(\boldsymbol{\widetilde{X}}'\boldsymbol{\widehat{L}}\boldsymbol{u})(\boldsymbol{\widetilde{X}}'\boldsymbol{\widehat{L}}\boldsymbol{u})'/T.$$

We extend the debiased Lasso estimator, introduced in javanmard2014confidence, and van2014asymptotically \textcolor{black}{to accommodate non-spherical errors while relaxing the Assumptions made for $\boldsymbol{x}_t$, to perform valid inference in high-dimensional linear models via a penalised GLS-type estimator, introduced in (ref).} The key idea of the debiased Lasso is to invert the Karush-Kuhn-Tucker (KKT) conditions of the Lasso solution defined in (ref). To proceed under correlated errors we invert the KKT conditions of (ref). Then,

where $\widehat{\boldsymbol{\kappa}}$ is obtained from the subgradient of the $\ell_{1}$-norm at $\boldsymbol{\widehat{\beta}}$, the last term of (ref), such that

The solution of the normal equations, can be further written as:

If $p>T$, then $\widehat{\boldsymbol{\Sigma}}$ is not invertible, so we seek a different avenue on approximating its inverse. Let $\widehat{\boldsymbol{\Theta}}$ be an approximate inverse of $\widehat{\boldsymbol{\Sigma}}$. Then, (ref) can be rewritten as

where $\boldsymbol{\delta}$ is the error resulting from the approximation of $\widehat{\boldsymbol{{\Theta}}}$, which is shown to be asymptotically negligible, in Theorem (ref). We use the fact that $ {T}^{-1}\widetilde{\boldsymbol{X}}'(\widetilde{\boldsymbol{y}} - \widetilde{\boldsymbol{X}}\boldsymbol{\widehat{\beta}}) = \lambda\widehat{\boldsymbol{\kappa}}$ and define the debiased GLS Lasso estimator

In essence, (ref) is the same as in equation (5) of van2014asymptotically, but using $(\widetilde{ \boldsymbol{y}},\widetilde{ \boldsymbol{X}})$ instead of $(\boldsymbol{y}, \boldsymbol{X})$. \textcolor{black}{Notice that, with (ref), we obtain the asymptotic pivotal quantity, $ ({\boldsymbol{\widehat{b }}} - \boldsymbol{\beta})$, since $\boldsymbol{\delta}=o_{P}(1)$, arriving at the following expression} \textcolor{black}{

The error $(\widehat{\boldsymbol{\Theta}}\widetilde{\boldsymbol{X}}{'}{\boldsymbol{\widehat{L}}\boldsymbol{u}})/\sqrt{T}$ }is asymptotically Gaussian and $\boldsymbol{\delta}$ is asymptotically negligible, as shown in Theorem (ref).

Construction of $\widehat{\boldsymbol{\Theta}}$

The most common practice of constructing the approximate inverse of $\widehat{\boldsymbol{\Sigma}}$, is given by node-wise regressions on the design matrix $\widetilde{\boldsymbol{X}}$. The main idea is similar to van2014asymptotically, but the asymptotic properties of the method change in the presence of serially correlated errors. Let $\widetilde{\boldsymbol{X}}$ be a $T\times p$ design matrix, and $\widetilde{\boldsymbol{X}}_{-i}$ a design sub-matrix, missing the $i^{th}$ column, for $i=1,\ldots,p$.

We first consider the population nodewise regressions, obtained by the linear projections

where $\boldsymbol{x}_{t}^{\star} = \boldsymbol{x}_t-\sum_{j=1}^q \phi_j\boldsymbol{x}_{t-j}$, and $( \tau^{\star}_{i})^{2}= E({{\upsilon_{i}^{\star}}})^{2} \label{ti2}$. Under Assumption (ref) and by consequence of Theorem 14.1 of davidson1994stochastic, we make the following Assumption for $ j=1,\ldots, p-1,\; i=1,\ldots,p$:

A feasible way to estimate $\boldsymbol{\gamma}_{i} $ would be to use $\boldsymbol{\widetilde{ \boldsymbol{x}}}_t$ instead of $\boldsymbol{x}^{\star}_{t}$ in (ref). Let $\widehat{\boldsymbol{\gamma}}_{i}$ be the estimates of the node-wise regressions, obtained as the solution of the following problem:

with components $\widehat{\boldsymbol{\gamma}}_{i}=\left\{ \widehat{{\gamma}}_{i,\kappa},\; \kappa\neq i,\; \kappa=1,\ldots,p \right \}$. Using the latter we define the following $p\times p$ matrix:

which is the solution to the series of the $p$ node-wise regressions, each of them tuned with a different regularisation parameter $\lambda_{i}\asymp \sqrt{\log{p}/T}$. We define $\widehat{\boldsymbol{\Theta}}$ as:

where

It remains to show that $\widehat{\boldsymbol{\Theta}}$ is a suitable approximation of the inverse of $\boldsymbol{\widehat{\Sigma}}$, which we explore in Proposition 10.1 of the Supplementary Material. Similarly to van2014asymptotically and kock2016oracle, we define the $i^{th}$ row of $\boldsymbol{\widehat{\Theta}}$ as $\boldsymbol{\widehat{\Theta}}_{i}$, a $1\times p$ vector and analogously $\widehat{\boldsymbol{C}}_i$. Hence, $\boldsymbol{\widehat{\Theta}}_i= \widehat{\boldsymbol{C}}_i\;\widehat{\tau}_{i}^{-2}$. The first order (KKT) conditions for the node-wise Lasso regressions in (ref), imply that

where $\boldsymbol{e}_{i}$ is a $p\times 1$ unit vector. Notice that the bound in (ref) is controlled by $\widehat{\tau}_i^{2}$, which depends explicitly on the solutions of the node-wise lasso regressions, $\widehat{\boldsymbol{\gamma}}_i$. In Lemma 10.2 of the Supplementary Material, we show that $\widehat{\boldsymbol{\gamma}}_i$ is consistently estimated, and $\lVert\boldsymbol{\gamma}_i-\widehat{\boldsymbol{\gamma}}_i\rVert_{\ell}, \; \ell=\{1,2\}$, attains a non-asymptotic bound. We further show that (ref) holds a.s., in Proposition 10.1 of the Supplementary Material.

Asymptotic properties of $\widehat{\boldsymbol{\Theta}}$

In order to show that $\sqrt{T}{\boldsymbol{\Theta}}\widetilde{\boldsymbol{X}}{'}{\boldsymbol{\widehat{L}u}} $ is asymptotically Gaussian, where $\widetilde{\boldsymbol{X}}$ is defined prior to (ref), one needs to explore the limiting behaviour of the approximate inverse $\widehat{ {\boldsymbol{\Theta}}}$. Consider ${\boldsymbol{\Sigma}} = E( \boldsymbol{\widetilde{x}}_t \boldsymbol{\widetilde{x}}_t')$ and $\boldsymbol{\Theta} $ its approximate inverse. Following a similar partition of $\boldsymbol{\Sigma}^{-1}$ as yuan2010high, we show that $\widehat{\boldsymbol{\Theta}}$ is asymptotically equivalent to ${\boldsymbol{\Theta}}$ under certain sparsity assumptions (formalised in Lemma 10.4 of the Supplementary Material). We write

where $\Theta_{i,i} $ is the $i^{th}$ element of the diagonal of $\boldsymbol{\Theta}$, $\Sigma_{i,i}$ is the $i^{th}$ element of the diagonal of $\boldsymbol{\Sigma}$, $\boldsymbol{\Sigma}_{i,-i}$ is the $1\times (p-1)$ vector of the design, obtained by removing the $i^{th }$ column of the $i^{th}$ row, $\boldsymbol{\Sigma}_{-i,-i}$ is the sub-matrix of $\boldsymbol{\Sigma}$ resulting from the removal of the $i^{th}$ row and column, $\boldsymbol{\Sigma}_{-i,i} $ is the $(p-1)\times 1$ vector of the design, obtained by removing the $i^{th}$ row of the same column and $\boldsymbol{\Theta}_{i,-i}$ is the $1\times (p-1)$ vector of $\boldsymbol{\Theta}$, obtained by removing the $i^{th }$ column of the $i^{th}$ row. \\ Let (ref) be the $p$ node-wise regressions. By linear projections, let its first order conditions be

Replacing $\widetilde{\boldsymbol{X}}_{-i,-i}^{-1}\;\widetilde{\boldsymbol{x}}_{-i,i}$ of (ref) into (ref), we have that

which illustrates that the $i^{{th}}$ row of $\boldsymbol{\Theta}$ is sparse, if and only if, \textcolor{black}{${\boldsymbol{\widetilde{\gamma}}_{i}}$} is sparse. Furthermore, let the error, $\boldsymbol{\psi}_{i}= \widetilde{\boldsymbol{x}}_{i}-\widetilde{\boldsymbol{X}}_{-i}\textcolor{black}{\boldsymbol{\widetilde{\gamma}}_{i}}$, then

Then, using Assumption (ref), and the definition of $\textcolor{black}{\boldsymbol{\widetilde{\gamma}}_{i}}$ as a penalised quadratic objective function we have that

for some small positive constant $c$. In light of Theorem (ref), it is sensible to note that $\widehat{\boldsymbol{\gamma}}_{i}$, defined in (ref) is close to its population counterpart, $\textcolor{black}{\boldsymbol{\widetilde{\gamma}}_{i}}$. We then can write

Notice that $\Theta_{i,-i}=-\textcolor{black}{\boldsymbol{\widetilde{\gamma}}_{i}}\tau_{i}^{-2}$. Therefore, replacing the terms in (ref), with their population counterparts, we can write $\boldsymbol{{M}}^{2} = \text{diag}(\tau_{1}^{2},\ldots, \tau_{p}^{2})$, where matrices $\boldsymbol{T}, \;\boldsymbol{C}$ are defined similarly to $\widehat{\boldsymbol{T}},\; \widehat{\boldsymbol{C}}$ replacing $\widehat{\boldsymbol{\gamma}}_{i}$ with $\textcolor{black}{\boldsymbol{\widetilde{\gamma}}_{i}}$, which in turn is asymptotically close to the true/population value of $\boldsymbol\gamma_i$ as it is shown in Lemma 10.4 of the Supplementary Material. Further, we illustrate in Lemma 10.4 of the Supplementary Material that $(\widehat{\tau}^{2}_i- \tau_{i}^2)\to_{P} 0$, such that, $\widehat{\boldsymbol{\Theta}}\to_{P} \boldsymbol{\Theta}$ a.s. for $p\to \infty$.

Asymptotic normality

We explore the asymptotic normality of the debiased GLS Lasso, and establish uniformly valid confidence intervals based on the theoretical results considered throughout Section (ref).

Following the results of Theorem (ref), we introduce asymptotic point-wise confidence intervals for the true parameter, $\beta_{i}$, $1\leq i\leq p$, given by

$\boldsymbol{\widehat{\Sigma}}_{xu} =T^{-1} (\boldsymbol{\widetilde{X}}'\boldsymbol{\widehat{L}}\boldsymbol{u})(\boldsymbol{\widetilde{X}}'\boldsymbol{\widehat{L}}\boldsymbol{u})'$, $z_{\alpha/2}:=\Phi^{-1}\left(1-\alpha/2\right)$, $\alpha $ is the confidence level, and $\Phi(\cdot)$ is the standard normal CDF, }such that

and $\widehat{\sigma}_{u}^2$ is a consistent estimate of the variance of ${u}_{t}$. Thus, one can perform inference on $\beta_{i}$, for some $z\in \mathbb{R}$, using the following form

for any $z\in \mathbb{R}$. Note that (ref) is implied in Theorem (ref) for a single-dimensional component, ${\beta}_i$. It can be generalised to hold trivially for a group $G$, of components, $S_G\subseteq\{1,\ldots,p\}$ which can be allowed to be large, see, for example, van2014asymptotically.

Selection of the optimal regularisation parameter

The selection of the regularisation parameter used to obtain $\boldsymbol{ \widehat{\beta}}$, $\boldsymbol{\widetilde{\beta}}$, and $\boldsymbol{\widehat{\gamma}}_i$, defined in (ref), (ref) and (ref) respectively, is carried out by a $k$-fold-cross validation scheme, where each fold can be considered as a block. Due to the autocorrelation present in the errors, we do not randomly shuffle the data within the blocks, retaining the structure of dependence intact. Each of the blocks $B_{1}, \ldots, B_{k}$ are non-overlapping and consequent to each other.

However, different types of cross validation schemes can be considered. racine2000consistent proposed a cross-validation approach carried out in blocks and highlighting near-independence of the training and validation sub-samples, as an extension to the $h$-block method of burman1994cross, which proved consistent for general stationary processes. The former encompasses our method of cross-validation, the leave-one-out method and the $h$-block method of burman1994cross, and can be an alternative in selecting the optimal regularisation parameter in our framework as well.

Simulation study

We use simulations to verify the theoretical properties of the proposed methodology in finite samples. We assess the performance of our method against the debiased Lasso. We generate an array of samples from the following model

where $\boldsymbol{x}_{t}, \; \varepsilon_{t}\sim\rm{i.i.d.} \, \mathcal{N}(0,1)$ and the sample sizes $T,p=\{100, 200, 500\}$. The active set has cardinality $s_{0} = |S_{0}|= 3$, where $S_{0}=\{ i:\beta_{i}\neq0, \; i=1,\ldots,p \}$, while $S_{0}^{c}= \{i:\beta_{i}=0, \; i=s_0+1,\ldots, p \}$. We choose the sparsity level $s_0$, according to the results of van2014asymptotically, which are indicative of good performance for the debiased Lasso. Simulation results for $s_0=7$ are relegated to the Supplementary Material. Each of the sets, $S_{0}, \;S_{0}^{c}$ assume the form: $ S_{0} = \{1,2,\ldots,s_{0}\}\equiv\{v_{1},\ldots,v_{s_{0}}\}, $ where $v_{1},\ldots,v_{s_{0}}$ is a realisation of random draws of $S_{0}$ without replacement from $\{1,\ldots,p\}$. The parameters, $\beta_{i}\in S_0$ are simulated from the $U[0,1]$ distribution at each replication, while the autoregressive parameter $\phi$ used for the simulation of $u_{t}$, takes values $\phi=[0,0.5,0.8,0.9]$, where $\phi=0$ indicates that $u_{t}\sim\rm{i.i.d.} \, \mathcal{N}(0,1)$.

We base our findings on 1000 Monte Carlo replications, and evaluate the following performance measures:

where $z=s_{0}, \; S= S_{0}$ when $\beta_{i}\in {S_{0}}$, $z=p-s_{0},\; S=S_{0}^{c}$, when $\beta_{i}\in {S_{0}^{c}}$, $\forall\; i=1,\ldots,p$ and CI$_{i}$ is a two-sided confidence interval for either $\beta_{i}\in {S_{0}}$ or $\beta_{i}\in {S_{0}^{c}}$, denoted in (ref). Further, we test the null hypothesis $H_{0} : \beta_i = 0$, using the following t-statistic

where, for all $i=1,\ldots, p$, $\widehat{S}_i\sim t^{(T-1)}$, and $ t^{(T-1)}$ the student's-$t$ distribution with $T-1$ degrees of freedom, $\beta_{i|H_0}=0$ under the null, $\widehat{b}_{i}$ is defined in (ref), $\widehat{\boldsymbol{\Sigma}}={T}^{-1}\widetilde{ \boldsymbol{x}}_t\widetilde{ \boldsymbol{x}}_t'$ is the sample covariance matrix, and $\boldsymbol{\widehat{\Theta}}$ is defined in Section (ref). Note that in the case where $i\in S_{0}$, we report the empirical power, while for $i\in S_{0}^{c}$, we report the empirical size

All tests are carried out at $\alpha=5\%$ significance level and confidence intervals are at the $1-\alpha=95\%$ level. Prior to the estimation of $\widehat{b}_{i}$, the estimate of ${\widehat{\phi}}$, and a value for the regularisation parameter are required. In Section (ref), we describe the method used for the estimation of ${\widehat{\phi}}$, along with the testing scheme used for the selection of the optimal lag order, $q$, of the $AR(q)$ model. Further, we use the $k$-fold cross-validation scheme described in Section (ref), where $k=10$, to select the optimal regularisation parameter, $\lambda$ to obtain $\boldsymbol{ \widehat{\beta}}$, $\boldsymbol{\widetilde{\beta}}$ and $\boldsymbol{\widehat{\gamma}}_i$ in (ref), (ref) and (ref), respectively. The results are displayed in Table (ref).

We also examine the performance of the Lasso, GLS Lasso, debiased Lasso and debiased GLS Lasso estimators in terms of average root mean squared error (RMSE) throughout 1000 replications of model (ref). The results are displayed in Table (ref).

The first panel of Table (ref) reports the ratio of the average RMSE of the Lasso estimator over the RMSE of the GLS Lasso. The second panel of Table (ref), reports the ratio of the average RMSE of the debiased Lasso estimator over the RMSE of debiased GLS Lasso. Entries larger than 1 indicate superiority of the competing model (GLS \emph{Lasso}). Highlighted are the entries corresponding to the RMSE of the {\emph{GLS} Lasso} in Panel I and the \emph{debiased \emph{GLS} Lasso} in Panel II.

The evidence is compelling, since both the GLS Lasso and debiased GLS Lasso have, in general, the smallest RMSE compared with their counterpart respectively. The most pronounced cases of improvement are when the autocorrelation is strong, i.e. $\phi=0.9$ for all possible sample sizes, where the RMSE of the proposed method is almost (and /or more) two times smaller than the one recorded by the benchmark, Lasso. Notice that in the case of $\phi=0$ GLS Lasso and Lasso perform equally well.

In the first panel of Table (ref) we present the average coverage rates (as $\text{AvgCov}$ $S_0$ or $S_{0}^{c}$) and average length of the confidence intervals, of the debiased GLS Lasso compared to the debiased Lasso. Notice that $1-\text{AvgCov} S_0^{c} $ is the size of the test in (ref). In the second panel of Table (ref) we facilitate the comparison of the two methods by reporting the size adjusted power of the test. By size-adjusted power, we mean that if the test is found to reject at a rate $\widehat{\alpha}>0.05$ under the null, then the power of the test is adjusted by finding the $0.95$ empirical quantile of the statistic in (ref). We then report the rate that exceeds this quantile. Alternatively, if $\widehat\alpha\leq 0.05$, then no adjustment is needed.

It is noticeable that when autocorrelation is present in the error term, debiased Lasso underperforms severely in terms of all measures reported. More specifically, as $T$ increases the benefits of debiased GLS Lasso are more clear, the coverage rates reported for our method are approaching the desirable rate of $0.95$, while for the debiased Lasso there are noticeable deviations from the nominal rate. The latter is more severe for higher serial autocorrelation reported, e.g. $\phi=[0.8,0.9]$.

Further, debiased GLS Lasso reports narrower confidence intervals, indicating that the variance of the proposed estimator is significantly smaller compared to the competing method. In addition, our method appears to be more powerful and correctly sized throughout the different cases reported. The best performing case reported is in Table (ref), for $(p,T)=(200,500)$ and $\phi=0.9$, where the size is $0.05$ and for $(p,T)= (500,500)$, using the same $\phi$, the reported size is $0.049$.

In Figures (ref)--(ref) we graph the estimation and prediction error rates of GLS Lasso, defined in (ref) and (ref) respectively, of Lasso, defined in (ref) and (ref), and of the sub-optimal GLS Lasso considering a $\lambda$, in the preliminary estimate $\boldsymbol{\widetilde{\beta}}$ which maximizes rather than minimising the following test-loss function w.r.t. the regularisation parameter, $\lambda$, for sample sizes $T=[100,500], \; p=[100,200,500]$ and $u_{t} = 0.9 u_{t-1}+\varepsilon_{t},$

where $B_k$ is one block (fold) for $k=1,\ldots,10$, $\boldsymbol{\widetilde{f}}_{\lambda}^{-i}\left(\boldsymbol{{X}}_{j}\right) =\boldsymbol{X}_{j} \boldsymbol{\widetilde{\beta}}$ excluding block $i$.

Notice that the prediction and estimation error for the GLS Lasso are bounded from below as Theorem (ref) indicates, while a lower bound for the Lasso is provided in Lemma (ref). In empirical problems, cross-validation does not always yield an optimal selection of $\lambda$, because one has to select between true model recovery (interpretable model) and a parsimonious (highly regularised) model. The latter can potentially be an issue because the selection of $\lambda$ controls the lower bound of both the prediction and estimation error, see, for example, (ref)--(ref). To that end, GLS Lasso provides an asymptotic guarantee that even when the selection of $\lambda$ is sub-optimal, it is preferable to using Lasso. Indicative cases of such behaviour are $T=100,p=200, 500$ and $T=100, p=100, 200, 500$. Furthermore, in the remaining sets of cases, of a sub-optimal selection of $\lambda$, GLS Lasso still reports smaller errors than Lasso but by a smaller margin.

As a summary of Figure (ref), GLS Lasso can be more attractive under relaxed assumptions due to the fact that neither the estimation, nor the prediction error deteriorate severely with different selections of $\lambda$. Further, even when the selection of $\lambda$ for the preliminary estimate is sub-optimal, GLS Lasso appears to be, in the majority of the cases, the best method to use as a preliminary vehicle to inference.

In Figure (ref) we graph the density estimates of $\sqrt{T}(\boldsymbol{\beta} - \boldsymbol{\widehat{{b}}} )$ corresponding to debiased GLS Lasso starting from the left panel, the second panel corresponds to debiased Lasso and the third panel illustrates the quantile-quantile (QQ) plots of the aforementioned. Figure (ref) corresponds to the case where the sample size is the largest, i.e. $T=500$ and the number of regressors is $p=200$ for all $\phi= [0, 0.5, 0.8, 0.9]$. A full set of graphs considering all the cases that are reported in the Tables (ref)--(ref) are available upon request.

In the first row of graphs in Figure (ref) where $\phi=0.9$, it is observed that extreme dependence between $u_t$ and $u_{t-1}$ is related to overdispersion of the distribution of $\sqrt{T}(\boldsymbol{\beta} - \boldsymbol{\widehat{{b}}} )$, corresponding to the debiased Lasso estimator, while in the same case, departures from the standard normal are observed, with the empirical data showing evidence of heavy tails. On the contrary, the empirical distribution of the debiased GLS Lasso estimator is closer to the Standard Normal, as it is indicative from the QQ-plot of the same row. However, this effect fades as $\phi$ takes values closer to $0$, indicative of the case where $u_{t}\sim \; i.i.d.\mathcal{N}(0,1). $ In the second row of Figure (ref), debiased Lasso appears slightly skewed left with heavy tails to persist indicative of departures from the standard normal distribution, while debiased GLS Lasso appears to approach the standard normal distribution, similar behaviour for both methods is reported in the third row. Finally, in the case where $\phi= 0 $ the two methods coincide supporting the evidence in Tables (ref)--(ref).

Discussion

This paper provides a complete inferential methodology for high dimensional linear regressions with serially correlated error terms. We propose GLS-type of estimators in the high-dimensional framework, derive a debiased GLS Lasso estimator, and establish its asymptotic normality. Theoretical results lead to valid confidence intervals regardless of the dimensionality of the parameters and the degree of autocorrelation, under stationarity of the error and covariate processes. We derive non-asymptotic bounds using Bernstein-type of inequalities derived in DGK, allowing for a degree of flexibility in our model compared to what classical assumptions typically allow. Finally, our method explores the underlying structure of the process in the error term, introducing a flexible way of dealing with autocorrelation in high-dimensional models, without a prior knowledge of neither its existence nor its structure. We support this statement both theoretically and with our simulation results, where we show that even when autocorrelation is absent the debiased Lasso and debiased GLS Lasso are asymptotically equivalent.

There are a number of interesting avenues for future work. For example, a natural extension would be to generalise our method using not only $\ell_1$-regularised models or $\ell_2$ (e.g. zhang2021debiased), but a combination of the two, e.g. the elastic net, covering a wider class of penalised models. Further, since our assumptions allow for a wide class of dependent stochastic processes, one could consider heteroskedastic error processes, see, for example, the framework in cck2022, where

where $h_t$ is a persistent scaling factor, smoothly varying in time, while $\epsilon_t\sim\; i.i.d.$

{\smallAcknowledgements: }{For discussions and/or comments we thank Jianqing Fan, Dacheng Xiu, Zacharias Psaradakis, Tassos Magdalinos and Yunyi Zhang.}

\setcounter{page}{1} \counterwithin*{equation}{section} \counterwithin{table}{section} \emptythanks \emptythanks

This Supplement provides proofs of the theoretical results given in the text of the main paper. It is organised as follows: Section (ref) provides proofs of Lemma \textcolor{red}{1}, Corollary \textcolor{red}{1}, Theorem \textcolor{red}{1} and Corollary \textcolor{red}{2} of the main paper. Section (ref) provides the proofs of Theorem \textcolor{red}{2} and Theorem \textcolor{red}{3} of the main paper. Section (ref) provides proofs of auxiliary technical lemmas. Section (ref) contains supplementary simulation material.

Formula numbering in this supplement includes the section number, e.g. $(8.1)$, and references to lemmas are signified as “Lemma 8.\#", “Lemma 9.\#", e.g. Lemma 8.1. Formula numbering in the Notation section is signified with the letter N.\#, e.g. (N.1). Theorem references to the main paper are signified, e.g. as Theorem $ \textcolor{red}{1}$, while equation references are signified as, e.g. \textcolor{red}{(1)}, \textcolor{red}{(2)} (denoted in red).

In the proofs, $C, d, e$ stand for generic positive constants which may assume different values in different contexts.

{ Notation: For any vector $\boldsymbol{x}\in \mathbb{R}^{n}, $ we denote the $\ell_{p}$-, $\ell_{\infty}$- and $\ell_0$- norms, as $\left\lVert \boldsymbol{x} \right\rVert_{p}=\left(\sum_{i=1}^{n} | x_{i} |^p\right)^{1/p}$, $\lVert \boldsymbol{x} \rVert_{\infty} = \max_{i=1,\ldots,n}|x_{i}|$ and $\lVert \boldsymbol{x}\rVert_{0}= \sum_{i=1}^{n}\boldsymbol{\mathrm{1}}_{(x_{i}\neq 0)} = \rm{supp}(\boldsymbol{\beta})$, respectively. Throughout this supplement, $\mathbb{R}$ denotes the set of real numbers. We denote the cardinality of a set $S_0$ by $s_{0}=|S_0|$, while $S_{0}^c$ denotes its complement. For any $\boldsymbol{x}\in \mathbb{R}^{n}$, $\rm{sign}(x)$ denotes the sign function applied to each component of $\boldsymbol{x}$, where $\mathrm{sign}(x)=0$ if $x=0$. We use "$\to_{P}$" to denote convergence in probability. For two deterministic sequences $a_n$ and $b_n$ we define asymptotic proportionality, "$\asymp$", by writing $a_n \asymp b_n$, if there exist constants $0 < a_1 \leq a_2$ such that $a_1b_n \leq a_n \leq a_2b_n$ for all $n \geq 1$. For a real number $a$, $\lfloor a \rfloor$ denotes the largest integer no greater than $a$. { For a symmetric matrix $\boldsymbol{B}$, we denote its minimum and maximum eigenvalues by $\lambda_{\min}\left( \boldsymbol{B} \right)$ and $\lambda_{\max}\left( \boldsymbol{B} \right)$ respectively. For a general (not necessarily square) $n\times m$ matrix $\boldsymbol{B}=(b_{ij})$, $\Vert \boldsymbol{B}\Vert_1$, $\Vert \boldsymbol{B}\Vert_F$, $\Vert \boldsymbol{B}\Vert_{\infty}$ are its $\ell_{1}$-, $\ell_{F}$-, $\ell_{\infty}$-norm respectively. In particular, $\Vert \boldsymbol{B}\Vert_1= \max_j\sum_{i}|b_{ij}|$, $\Vert \boldsymbol{B}\Vert_F = [\rm{tr}(\boldsymbol{B}\boldsymbol{B}')]^{1/2}$, $\Vert \boldsymbol{B}\Vert_{\infty} = \max_i\sum_{j} |b_{ij}|$ and $\rm{tr}(\boldsymbol{B}) $$= \sum_{i=1}^{n}b_{ii}$.} We further write that, $\widetilde{{\boldsymbol{\beta}}}$ is the Lasso estimator obtained from the solution of \textcolor{red}{(13)} in the main paper, $\boldsymbol{\widehat{\beta}}$ is the GLS Lasso estimator obtained from the solution of \textcolor{red}{(3)}, given some estimates of the autoregressive parameters ${ \boldsymbol{\widehat{\phi}}} = ( \widehat{\phi}_{1}, \ldots, \widehat{\phi}_{q} )'$ and $\boldsymbol{\beta}$ is a $p\times 1$ parameter vector containing the true parameters.

We define $\widetilde{\boldsymbol{x}}_t= \boldsymbol{x}_t- \sum_{ j=1}^{q}\widehat\phi_{j}\boldsymbol{x}_{t-j}$, and its matrix counterpart, $ \widetilde{\boldsymbol{X}}= \boldsymbol{\widehat{L} X} $, where $\widehat{\boldsymbol{L}}$ is a $\left(T-q \right)\times T$ matrix defined as:

Further, following the same argument that we used to define $\widehat{\boldsymbol{ L }}$, we define ${\boldsymbol{L}}$, a $\left(T-q \right)\times T$ matrix, with $\widehat{\phi}_{j}$ replaced by $\phi_{j}$, denoting the autoregressive parameters, $\forall\; j=1,\ldots,q$. As a consequence, we define ${\boldsymbol{X} }^{\star} = {\boldsymbol{L}}\boldsymbol{ X }$.

We shall frequently refer to the $\alpha$-mixing Assumption \textcolor{red}{2} and its property $\alpha_{k}\leq c_{*}\phi^{k}, \; k\geq 1$ of $\{\boldsymbol{x}_{t}\}$ and $\{u_{t } \}$ of Section \textcolor{red}{2} of the main paper. The r.v.'s $\boldsymbol{x}_{t}, \; {u}_{t}$ have thin-tailed distributions with properties denoted in \textcolor{red}{(5)}--\textcolor{red}{(6)} of the main paper.

Proofs of Lemma 1, Corollary 1, Theorem 1 and Corollary 2 of the main paper

This section contains the proofs of the results of Section \textcolor{red}{2} of the main paper.

}

Proofs of Theorems 2 and 3 of the main paper

This Section provides proofs of Theorems \textcolor{red}{2}, and \textcolor{red}{3} of the main paper.

Proofs of auxiliary results

This section contains auxiliary technical lemmas used in the proofs of Section (ref).