Lasso Inference for High-Dimensional Time Series

Lasso Inference for High-Dimensional Time Series

\onehalfspacing

Introduction

In this paper we propose methods for performing uniformly valid inference on high-dimensional time series regression models. Specifically, we establish the uniform asymptotic normality of the desparsified lasso method vandeGeer14 under very general conditions, thereby allowing for inference in high-dimensional time series settings that encompass many econometric applications. That is, we establish validity for potentially misspecified time series models, where the regressors and errors may exhibit serial dependence, heteroskedasticity and fat tails. In addition, as part of our analysis we derive new error bounds for the lasso Tibshirani96, on which the desparsified lasso is based.

Although traditionally approaches to high-dimensionality in econometric time series have been dominated by factor models BaiNg08survey,StockWatson11, shrinkage methods have rapidly been gaining ground. Unlike factor models where dimensionality is reduced by assuming common structures underlying regressors, shrinkage methods assume a certain structure on the parameter vector. Typically, sparsity is assumed, where only a small, unknown subset of the variables is thought to have “significantly non-zero” coefficients, and all the other variables have negligible -- or even exactly zero -- coefficients. The most prominent among shrinkage methods exploiting sparsity is the lasso proposed by Tibshirani96, which adds a penalty on the absolute value of the parameters to the least squares objective function. This penalty ensures that many of the coefficients will be set to zero and thus variable selection is performed, an attractive feature that helps to make the results of a high-dimensional analysis interpretable. Due to this feature, the lasso and its many extensions are now standard tools for high-dimensional analysis Hesterberg08,Vidaurre13,Hastie15.

Much effort has been devoted to establish error bounds for lasso-based methods to guarantee consistency for prediction (e.g., greenshtein2004persistence, buehlmann2006boosting) and estimation of a high-dimensional parameter (e.g., bunea2007sparsity, zhang2008sparsity, bickel2009simultaneous, meinshausen2009lasso, Huang08). While most of these advances have been made in frameworks with independent and identically distributed (IID) data, early extensions of lasso-based methods to the time series case can be found in Wang07, Hsu08. These authors, however, only consider the case where the number of variables is smaller than the sample size. Various papers (e.g., NardiRinaldo11; KockCallot2015 and BasuMichailidis15) let the number of variables increase with the sample size, but often require restrictive assumptions (for instance Gaussianity) on the error process when investigating theoretical properties of lasso-based estimators in time series models.

Exceptions are MedeirosMendes16, WuWu2016, masini2019regularized, and wong2020lasso. MedeirosMendes16 consider the adaptive lasso for sparse, high-dimensional time series models and show that it is model selection consistent and has the oracle property, even when the errors are non-Gaussian and conditionally heteroskedastic. WuWu2016 consider high-dimensional linear models with dependent non-Gaussian errors and/or regressors and provide asymptotic theory for the lasso with deterministic design. To this end, they adopt the functional dependence framework of Wu05. masini2019regularized focus on weakly sparse high-dimensional vector autoregressions for a class of potentially heteroskedastic and serially dependent errors, which encompass many multivariate volatility models. The authors derive finite sample estimation error bounds for the parameter vector and establish consistency properties of lasso estimation. wong2020lasso derive nonasymptotic inequalities for estimation error and prediction error of the lasso without assuming any specific parametric form of the DGP. The authors assume the series to be either $\alpha$-mixing Gaussian processes or $\beta$-mixing processes with sub-Weibull marginal distributions thereby accommodating settings with heavy-tailed non-Gaussian errors.

While one of the attractive feature of lasso-type methods is their ability to perform variable selection, this also causes serious issues when performing inference on the estimated parameters. In particular, performing inference on a (data-driven) selected model, while ignoring the selection, causes the inference to be invalid. This has been discussed by, among others, LeebPoetscher05 in the general context of model selection and LeebPoetscher08 for shrinkage estimators. As a consequence, recent statistical literature has seen a surge in the development of so-called post-selection inference methods that circumvent the problem induced by model selection; see for example the literature on selective inference FST15,LSST16 and simultaneous inference PoSI13,bachoc2020.

In the context of lasso-type estimation, methods have been developed based on the idea of orthogonalizing the estimation of the parameter of interest to the estimation (and potential incorrect selection) of the other parameters. BCH14,CHS15ARE propose a post-double-selection approach that uses a Frisch-Waugh partialling out strategy to achieve this orthogonalization by selecting important covariates in initial selection steps on both the dependent variable and the variable of interest, and show this approach yields uniformly valid and standard normal inference for independent data. In a related approach, JavanmardMontanari14; vandeGeer14 and ZhangZhang14 introduce debiased or desparsified versions of the lasso that achieve uniform validity based on similar principles for IID Gaussian data. Extensions to the time series case include chernozhukov2021timeandspace who provide desparsified simultaneous inference on the parameters in a high-dimensional regression model allowing for temporal and cross-sectional dependency in covariates and error processes, Krampe18 who introduce bootstrap-based inference for autoregressive time series models based on the desparsification idea, HMS19 who use the post-double-selection procedure of BCH14 for constructing uniformly valid Granger causality test in high-dimensional VAR models, and ghysels20 who use a debiased sparse group lasso for inference on a low dimensional group of parameters.

In this paper, we contribute to the literature on shrinkage methods for high-dimensional time series models by providing novel theoretical results for both point estimation and inference via the desparsified lasso. We consider a very general time series-framework where the regressors and errors terms are allowed to be non-Gaussian, serially correlated and heteroskedastic, and the number of variables can grow faster than the time dimension. Moreover, our assumptions allow for both correctly specified and misspecified models, thus providing results relevant for structural interpretations if the overall model is specified correctly, but not limited to this.

We derive error bounds for the lasso in high-dimensional, linear time series models under mixingale assumptions and a weak sparsity assumption on the parameter vector. Our setting generalizes the one from MedeirosMendes16, who require a martingale difference sequence (m.d.s.) assumption -- and hence correct specification -- on the error process. Moreover, we relax the traditional sparsity assumption to allow for weak sparsity, thereby recognizing that the true parameters are likely not exactly zero. The error bounds are used to establish estimation and prediction consistency even when the number of parameters grows faster than the sample size.

We extend the error bounds to the nodewise regressions performed in the desparsified lasso, where each regressor (on which inference is performed) is regressed on all other regressors. Note that, contrary to the setting with independence over time, these nodewise regressions are inherently misspecified in dynamic models with temporal dependence. As such our error bounds are specifically derived under potential misspecification. We then establish the asymptotic normality of the desparsified lasso under general conditions. As such, we ensure uniformly valid inference over the class of weakly sparse models. This result is accompanied by a consistent estimator for the long run variance, thereby providing a complete set of tools for performing inference in high-dimensional, linear time series models. As such, our theoretical results accommodate various financial and macro-economic applications encountered by applied researchers.

The remainder of this paper is structured as follows. (ref) introduces the time series setting and assumptions thereof. In (ref), we derive an error bound for the lasso ((ref)) that forms the basis for the nodewise regressions performed for the desparsfied lasso. In (ref), we establish the theory that allows for uniform inference with the desparsified lasso. (ref) contains a simulation study examining the small sample performance of the desparsified lasso, and (ref) concludes. The main proofs and preliminary lemmas needed for (ref) are contained in Appendix (ref), while Appendix (ref) contains the results and proofs on (ref). Appendix C contains supplementary material.

A word on notation. For any $N$ dimensional vector $\boldsymbol{x}$, $\left\Vert \boldsymbol{x}\right\Vert_r=\left(\sum\limits_{i=1}^{N}\left\vert x_i\right\vert^r\right)^{1/r}$ denotes the $L_r$-norm, with the familiar convention that $\left\lVert\boldsymbol{x}\right\rVert_0 = \sum_{i} 1 (\left\lvertx_i\right\rvert>0)$ and $\left\Vert \boldsymbol{x}\right\Vert_{\infty}=\max\limits_{i}\left\vert x_i\right\vert$. For a matrix $\boldsymbol{A}$, we let $\left\lVert\boldsymbol{A}\right\rVert_r = \max_{\left\lVert\boldsymbol{x}\right\rVert_r = 1} \left\lVert\boldsymbol{A} \boldsymbol{x}\right\rVert_r$ for any $r \in [0, \infty]$ and $\left\lVert\boldsymbol{A}\right\rVert_{\max}=\max\limits_{i,j}\left\vert a_{i,j}\right\vert$. We use $\overset{p}{\to}$ and $\overset{d}{\to}$ to denote convergence in probability and distribution respectively. Depending on the context, $\sim$ denotes equivalence in order of magnitude of sequences, or equivalence in distribution. We frequently make use of arbitrary positive finite constants $C$ (or its sub-indexed version $C_i$) whose values may change from line to line throughout the paper, but they are always independent of the time and cross-sectional dimension. Similarly, generic sequences converging to zero as $T\to\infty$ are denoted by $\eta_T$ (or its sub-indexed version $\eta_{T,i}$). We say a sequence $\eta_T$ is of size $-x$ if $\eta_T=O\left(T^{-x-\varepsilon}\right)$ for some $\varepsilon>0$.

The High-Dimensional Linear Model

Consider the linear model

where $\boldsymbol{x}_t=\left(x_{1,t},\dots, x_{N,t}\right)'$ is a $N\times 1$ vector of explanatory variables, $\boldsymbol\beta^0$ is a $N\times 1$ parameter vector and $u_t$ is an error term. Throughout the paper, we examine the high-dimensional time series model where $N$ can be larger than $T$.

We impose the following assumptions on the processes $\{\boldsymbol{x}_t\}$ and $\{u_t\}$.

(ref)(ref) ensures that the error terms are contemporaneously uncorrelated with each of the regressors, and that the process has finite and constant unconditional moments. One can think of $\boldsymbol{s}_{T,t}$ in (ref)(ref) as an underlying shock process driving the regressors and errors in $\boldsymbol{z}_t$, where we assume $\boldsymbol{z}_t$ to depend almost entirely on the “near epoch” of $s_{T,t}$.\footnote{Since $\boldsymbol{z}_t$ grows asymptotically in dimension, it is natural to let the dimension of $\boldsymbol{s}_{T,t}$ grow with $T$, though this is not theoretically required. Although, like $\boldsymbol{s}_{T,t}$, technically our stochastic process $\boldsymbol{z}_t$ is a triangular array due to dimension $N$ increasing with $T$, in the remainder of the paper we suppress the dependence on $T$ for notational convenience.}

Near epoch dependence of $\boldsymbol{z}_t$ can be interpreted as $\boldsymbol{z}_t$ being “approximately” mixing, in the sense that it can be well-approximated by a mixing process. The NED framework in (ref) therefore allows for very general forms of dependence that are often encountered in econometrics applications including, but not limited to, strong mixing processes McLeish75, linear processes including ARMA models, various types of stochastic volatility and GARCH models hansen1991garch, and nonlinear processes davidson2002establishing. Moreover, NED holds in cases where mixing has well-known failures for common processes, such as the AR(1) process discussed in andrews1984non. These properties have made NED a very popular tool for modelling dependence in econometrics Davidson02.\footnote{To make the paper self-contained, we include formal definitions on NED and mixingales in Appendix (ref).}

To our knowledge, our paper is the first to utilize the NED framework for establishing uniformly valid high-dimensional inference. wong2020lasso consider time series models with $\beta$-mixing errors, which has the advantage of allowing for general forms of dynamic misspecification resulting in serially correlated error terms, but, as discussed above, rules out several relevant data generating processes, and is in addition typically difficult to verify. Alternative approaches that avoid mixing assumptions are found in ghysels20, who consider $\tau-$dependence, as well as WuWu2016 and chernozhukov2021timeandspace, who use functional dependence for modeling the dependence allowed in regressors and innovations. Finally, masini2019regularized use an m.d.s. assumption on the innovations in combination with sub-Weibull tails and a mixingale assumption on the conditional covariance matrix. The m.d.s. assumption of MedeirosMendes16 and masini2019regularized however does not allow for dynamic misspecification of the full model. Importantly, the NED assumption on $u_t$ does allow for misspecified models as well, in which case we view $\boldsymbol{\beta}_0$ as the coefficients of the pseudo-true model when restricting the class of models to those linear in $\boldsymbol{x}_t$. In particular, it allows one to view (ref) as simply the linear projection of $y_t$ on all the variables in $\boldsymbol{x}_t$, with $\boldsymbol{\beta}^0$ in that case representing the corresponding best linear projection coefficients. In such a case $\mathbb{E}\left[u_t\right]=0$ and $\mathbb{E}\left[u_t x_{j,t}\right]=0$ hold by construction, and the additional conditions of (ref) can be shown to hold under weak further assumptions. On the other hand, $u_t$ is not likely to be an m.d.s. in that case. As will be explained later, allowing for misspecified dynamics is crucial for developing the theory for the nodewise regressions underlying the desparsified lasso.

It is important to note that we do not consider $\boldsymbol{\beta}^0$ as the projection coefficients of the (lasso) selected model, but only of the full, pseudo-true, model. Our approach simply allows for the possibility of the full model being misspecified, for instance if the econometrician has missed relevant confounders in the initial dataset. This does not imply a “failure” of our lasso inference method, but rather a failure of the econometrician in setting up the initial model.\footnote{Of course, the misspecification may be intentional, as even in dynamically misspecified models, the parameter of interest can still have a structural meaning. One example is the local projections of jorda2005estimation, where $h$-step ahead predictive regressions with generally serially correlated error terms are performed.} Allowing for such misspecification is crucial for the nodewise regressions we consider in Section (ref) which are simply projections of one explanatory variable on all the others, and therefore inherently misspecified.

We further elaborate on misspecification in (ref), after we present two examples of correctly specified common econometric time series DGPs.

(ref) demonstrate that (ref) is sufficiently general to include common time series models in econometrics. While these examples are equally well covered by other commonly used assumptions such as the martingale difference sequence (m.d.s) framework chosen in MedeirosMendes16 or masini2019regularized, we opt for the more general NED framework, as it additionally covers many relevant cases -- in particular for our nodewise regressions -- where properties such as m.d.s.\ fail. The following examples provide simple illustrations of these cases.

The key condition to apply the lasso successfully is that the parameter vector $\boldsymbol{\beta}_0$ is (at least approximately) sparse. We formulate this in (ref) below.

(ref) implies that ${\boldsymbol{\beta}}^0$ is sparse with the degree of sparsity governed by both $r$ and $s_r$. Without further assumptions on $r$ and $s_r$, (ref) is not binding, but as will be seen later, the allowed rates will interact with other DGP parameters creating binding conditions. (ref) generalizes the common assumption of exact sparsity taking $r=0$ (see e.g., MedeirosMendes16; vandeGeer14; chernozhukov2021timeandspace; ghysels20), which assumes that there are only a few (at most $s_0$) non-zero components in $\boldsymbol\beta^0$, to weak sparsity (see e.g., vandeGeer19). This allows us to have many non-zero elements in the parameter vector, as long as they are sufficiently small. It follows directly from the formulation in (ref) that, given the compactness of the parameter space, exact sparsity of order $s_0$ implies weak sparsity with $r >0$ of the same order (up to a fixed constant). In general, the smaller $r$ is, the more restrictive the assumption. The relaxation to weak sparsity is straightforward and follows from elementary inequalities (see e.g., Section 2.10 of vandeGeer2016book and the proof of (ref)).

For $\lambda\geq0$, define the weak sparsity index set

and complement set $S^c_\lambda=\left\lbrace1,\dots,N\right\rbrace\setminus S_\lambda$. With an appropriate choice of $\lambda$, this set contains all `sufficiently large' coefficients; for $\lambda=0$ it contains all non-zero parameters. We need this set in the following condition, which formulates the standard compatibility conditions needed for lasso consistency SHD11.

The compatibility constant in (ref) is an upper bound on the minimum eigenvalue of ${\boldsymbol{\Sigma}}$, so this condition is considerably weaker than assuming ${\boldsymbol{\Sigma}}$ to be positive definite. We formulate the compatibility condition in (ref) on the population covariance matrix rather than directly on the sample covariance matrix $\hat{\boldsymbol{\Sigma}}:=\boldsymbol{X}^\prime\boldsymbol{X}/T$, see e.g., the restricted eigenvalue condition in MedeirosMendes16 or Assumption (A2) in chernozhukov2021timeandspace. Verifying this assumption on the population covariance matrix is generally more straightforward than directly on the sample covariance matrix.\footnote{Though note that BasuMichailidis15 show in their Proposition 3.1 that the restricted eigenvalue condition holds with high probability under general time series conditions when $\boldsymbol{x}_t$ is a stable process with full-rank spectral density and $T$ is sufficiently large. Their Proposition 4.2 includes a stable VAR process as an example.}

Finally, note that the compatibility assumption for the weak sparsity index set $S_\lambda$ is weaker than (and implied by) its equivalent for $S_0$, see Lemma 6.19 in SHD11, and that the strictness of this assumption depends on the choice of the tuning parameter $\lambda$.

Error Bound and Consistency for the Lasso

In this section, we derive a new error bound for the lasso in a high-dimensional time series model. The lasso estimator Tibshirani96 of the parameter vector $\boldsymbol\beta^0$ in Model (ref) is given by

where $\boldsymbol{y}=(y_1, \ldots, y_T)^\prime$ is the $T\times 1$ response vector, $\boldsymbol{X}=\left(\boldsymbol{x}_1,\dots,\boldsymbol{x}_T\right)'$ the $T\times N$ design matrix and $\lambda>0$ a tuning parameter. Optimization problem (ref) adds a penalty term to the least squares objective to penalize parameters that are different from zero.

When deriving this error bound, one typically requires that $\lambda$ is chosen sufficiently large to exceed the empirical process $\max\limits_{j}\left\lvert\frac{1}{T}\sum_{t=1}^Tx_{j,t}u_t\right\rvert$ with high probability. To this end, we define the set $\mathcal{E}_{T}(z):=\left\lbrace\max\limits_{j\leq N,l\leq T}\left\lvert\sum\limits_{t=1}^{l}u_t x_{j,t}\right\rvert\leq z\right\rbrace$, and establish the conditions under which $\mathbb{P}\left(\mathcal{E}_{T}(T\lambda/4)\right)\to 1$. In addition, since we formulate the compatibility condition in (ref) on the population covariance matrix, we need to show that $\boldsymbol{\Sigma}$ and $\hat\boldsymbol{\Sigma}$ are sufficiently close under the DGP assumptions. To this end, we define the set $\mathcal{CC}_T(S):=\left\lbrace\left\lVert\hat\boldsymbol{\Sigma}-\boldsymbol{\Sigma}\right\rVert_{\max}\leq C/\left\lvertS\right\rvert\right\rbrace$, and show that $\mathbb{P}\left(\mathcal{CC}_T(S_{\lambda})\right)\to1$. (ref) then presents both results.

(ref) thus establishes that the sets $\mathcal{E}_{T}(T\lambda/4)$ and $\mathcal{CC}_T(S_\lambda)$ hold with high probability. Each set has a condition under which its probability converges to 1, which follow from (ref) respectively. For the set $\mathcal{E}_{T}(T\lambda/4)$, the condition $\lambda\geq{ C\ln(\ln( T))}^{1/m}\frac{N^{1/m}}{\sqrt{T}}$ is required. The $\ln(\ln( T))$ appearing throughout the theorem is chosen arbitrarily as a sequence which grows slowly as $T\to\infty$; we only need some sequence tending to infinity sufficiently slowly. The details can be found in the proof of (ref). For the set $\mathcal{CC}_T(S_\lambda)$, we need to distinguish the cases $0<r<1$ and $r=0$ due to the way the size of the sparsity index set in (ref) is bounded. For $0<r<1$, a lower bound on $\lambda$ is imposed which is stricter than the one for the empirical process, hence only that bounds appears in (ref). For $r=0$, the conditions do not depend on $\lambda$ hence both bounds appear in (ref).

(ref) directly yields an error bound for the lasso in high-dimensional time series models by standard arguments in the literature, see e.g., Chapter 2 of vandeGeer2016book. The proofs of Lemmas A.6 and A.7 in the Supplementary Appendix C.1 provide details.

Under the additional assumption that $\lambda^{1-r}s_{r}\to0$, these error bounds directly establish prediction and estimation consistency. The bounds in (ref) thereby put implicit limits on the divergence rate of $N$, and $s_{r}$ relative to $T$. In particular, the term offsetting the divergence in $N$, and $s_{r}$ is of polynomial order in $T$. The order of the polynomial, and therefore the restriction on the growth of $N$ and $s_{r}$, is determined by the moments $m$ and dependence parameter $d$; the higher the number of moments $m$ and the larger the dependence parameter $d$, the fewer restrictions one has on the allowed polynomial growth of $N$ and $s_r$. In the limit, if $m$ and $d$ tend to infinity (all moments exist and the data are mixing), the order of the polynomial restriction on $N$ tends to infinity, thereby approaching exponential growth. A similar trade off between the allowed growth of $N$ and the existence of moments was found in MedeirosMendes16. In Example C.1 we study in greater detail how the different rates interact, thereby providing an overview of the restrictions under different scenarios.

While (ref) is a useful result in its own right, it is vital to derive the theoretical results for the desparsified lasso, which we turn to next.

Uniformly Valid Inference via the Desparsified Lasso

We use the desparsified lasso to perform uniformly valid inference in general high-dimensional time series settings. After briefly reviewing the desparsified lasso, we formulate the assumptions needed in (ref). The asymptotic theory is then derived in (ref) for inference on low-dimensional parameters of interest, and (ref) for inference on a high-dimensional parameters.

The desparsified lasso vandeGeer14 is defined as

where $\hat{\boldsymbol{\beta}}$ is the lasso estimator from (ref) and $\hat{{\boldsymbol{\Theta}}}:=\hat{\boldsymbol{\Upsilon}}^{-2}\hat{\boldsymbol{\Gamma}}$ is a reasonable approximation for the inverse of $\hat\boldsymbol{\Sigma}$. By de-sparsifying the initial lasso, the bias in the lasso estimator is removed and uniformly valid inference can be obtained. The matrix $\hat{\boldsymbol{\Gamma}}$ is constructed using nodewise regressions; regressing each column of $\boldsymbol{X}$ on all other explanatory variables using the lasso. Let the lasso estimates of the $j=1,\dots,N$ nodewise regressions be

where the $T\times(N-1)$ matrix $\boldsymbol{X}_{-j}$ is $\boldsymbol{X}$ with its $j$th column removed. Their components are given by $\hat{\boldsymbol{\gamma}}_j=\left\lbrace\hat\gamma_{j,k}:k=\{1,\dots,N\}\setminus j\right\rbrace$. Stacking these estimated parameter vectors row-wise with ones on the diagonal gives the matrix

We then take $\hat{\boldsymbol{\Upsilon}}^{-2}:=\text{diag}\left(1/\hat{\tau}_1^2,\dots,1/\hat{\tau}_N^{2}\right)$, where $\hat{\tau}_j^2:=\frac{1}{T} \left\lVert \boldsymbol{x}_j-\boldsymbol{X}_{-j} \hat{\boldsymbol{\gamma}}_j\right\rVert_2^2 + 2\lambda_j \left\lVert\hat{\boldsymbol{\gamma}}_j\right\rVert_1$.

We use the index set $H\subseteq \left\lbrace 1,\dots,N\right\rbrace$ with cardinality $h=\left\lvertH\right\rvert$ to denote the set of variables whose coefficients we wish to perform inference on. In this case computational gains can be obtained with respect to the nodewise regressions, as we only need to obtain the sub-vector of the desparsified lasso corresponding to $\hat{\boldsymbol{b}}_H:=\hat{\boldsymbol{\beta}}_H + \hat{\boldsymbol{\Theta}}_H \boldsymbol{X}({\boldsymbol{y}}-\boldsymbol{X}\hat{{\boldsymbol{\beta}}})$, with the subscript $H$ indicating that we only take the respective rows of $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{\Theta}}$. To compute $ \hat{\boldsymbol{\Theta}}_H$, one only needs to compute $h$ nodewise regressions instead of $N$, which can be a considerable reduction for small $h$ relative to large $N$.

Assumptions

Consider the population nodewise regressions defined by the linear projections

with $\tau_j^2:=\frac{1}{T}\sum\limits_{t=1}^T\mathbb{E}\left[v_{j,t}^2\right]$. Note that by construction, it holds that $\mathbb{E}\left[v_{j,t}\right]=0,\ \forall t, j$ and $\mathbb{E}\left[v_{j,t}x_{k,t}\right]=0,\ \forall t, k\neq j$. We first present (ref), which allow us to extend (ref) to the nodewise lasso regressions.

(ref) requires the errors $v_{j,t}$ from the nodewise linear projections to have bounded moments of an order greater than fourth. By the properties of NED processes, we use (ref) to establish mixingale properties of the products $v_{j,t}u_t=:w_{j,t}$ and $w_{j,t}w_{k,t-l}$ in (ref), which are used extensively in the derivation of the desparsified lasso's asymptotic distribution.

(ref)(ref), similar to (ref), requires weak sparsity of the nodewise regressions, not exact sparsity. The latter could be problematic, as it would imply many of the regressors to be uncorrelated. In contrast, weak sparsity is a plausible alternative, see e.g.,\ (ref). Importantly, the weak sparsity of the nodewise regressions is fully determined by the model and hence should be verified. Below, we provide concrete examples where the weak sparsity assumption holds.

(ref)(ref) requires the population covariance matrix to be positive definite, with its smallest eigenvalue bounded away from zero, and to have finite variances. (ref)(ref) implies the compatibility condition and thus replaces (ref) in (ref), with $\Lambda_{min}$ fulfilling the role of $\phi_{{\boldsymbol{\Sigma}}}^2$. It also implies that the explanatory variables, including the irrelevant ones, cannot be linear combinations of each other even as we let the number of variables tends to infinity. Although this is a considerable strengthening of (ref), it is important to realize this assumption is still made on the population matrix instead of the sample version, and may therefore still hold in fairly general, high-dimensional models. For example, BasuMichailidis15 provide a lower bound for $\Lambda_{\min}$ in VAR models on their Proposition 2.3, which can be shown to be bounded away from zero under realistic conditions, see also masini2019regularized (masini2019regularized, p.\ 6). Similarly, this assumption can be shown to hold in factor models under minimal assumptions on the idiosyncratic errors (see (ref) below).

(ref) allow us to apply (ref) to the nodewise regressions. Specifically, if the conditions on $\lambda$ formulated in (ref) hold for both $\underset{\bar{}}{\lambda}:=\min\limits_{j\in H}\lambda_j$ and $\bar{\lambda}:=\max\limits_{j\in H}\lambda_j$, the error bounds -- with $\bar{s}_r:=\max\limits_{j\in H} s^{(j)}_{r}$ substituted for $s_r$ -- apply to the nodewise regressions as well. As we generally need the error bounds to hold uniformly over all relevant nodewise regressions as well as the initial regression, we combine these bounds and state our results on the quantities

which simplifies many of the final expressions. While some conditions could be weakened if we keep them in terms of $\bar{\lambda}$ or $\bar{s}_r$ explicitly, this would be at the expense of more conditions and readability, and therefore we opt against it.

Inference on low-dimensional parameters

In this section we establish the uniform asymptotic normality of the desparsified lasso focusing on low-dimensional parameters of interest. We consider testing $P$ joint hypotheses of the form $\boldsymbol{R}_{N}\boldsymbol{\beta}^0=\boldsymbol{q}$ via a Wald statistic, where $\boldsymbol{R}_{N}$ is an appropriate $P\times N$ matrix whose non-zero columns are indexed by the set $H:=\left\lbrace j:\sum_{p=1}^P\vert r_{N,p,j}\vert>0\right\rbrace$ of cardinality $h:=\vert H\vert$. As can be seen from the lemmas in Appendix (ref), all our results up to application of the central limit theorem allow for $h$ to increase in $N$ (and therefore $T$). In (ref) we first focus on inference on a finite set of parameters, such that we can apply a standard central limit theorem under the assumptions listed above. An alternative, high-dimensional approach under more stringent conditions is considered in (ref).

Given our time series setting, the long-run covariance matrix

where $\boldsymbol{w}_t=(v_{1,t}u_t,\dots,v_{N,t}u_t)'$, enters the asymptotic distribution in (ref). $\boldsymbol{\Omega}_{N,T}$ can equivalently be written as $\boldsymbol{\Omega}_{N,T}=\boldsymbol{\Xi}(0)+\sum\limits_{l=1}^{T-1}(\boldsymbol{\Xi}(l)+\boldsymbol{\Xi}^\prime(l))$, where $\boldsymbol{\Xi}(l) = \frac{1}{T}\sum\limits_{t=l+1}^T\mathbb{E}\left[ \boldsymbol{w}_t \boldsymbol{w}_{t-l}^\prime\right]$.

In order to estimate the asymptotic variance $\boldsymbol{\Psi}$, we suggest to estimate $\boldsymbol{\Omega}_{N,T}$ with the long-run variance kernel estimator

where $\hat{\boldsymbol{\Xi}}(l)=\frac{1}{T-l}\sum\limits_{t=l+1}^{T}\hat{\boldsymbol{w}}_t \hat{\boldsymbol{w}}_{t-l}^\prime$ with $\hat{w}_{j,t}=\hat{v}_{j,t}\hat{u}_t$, the kernel $K(\cdot)$ can be taken as the Bartlett kernel $K(l/Q_T) = \left(1-\frac{l}{Q_T}\right)$ NeweyWest87 and the bandwidth $Q_T$ should increase with the sample size at an appropriate rate. A similar heteroskedasticity and autocorrelation consistent (HAC) estimator was considered by ghysels20, though under a different framework of dependence. In (ref), we show that $\hat{\boldsymbol{\Psi}} = \boldsymbol{R}_{N}({\hat{\boldsymbol{\Upsilon}}}^{-2}\hat{\boldsymbol{\Omega}}{\hat{\boldsymbol{\Upsilon}}}^{-2}){\boldsymbol{R}'_{N}}$ is a consistent estimator of $\boldsymbol{\Psi}$ in our NED framework.

Note that here we restrict $\boldsymbol{R}_N$ such that the number of hypotheses $P$ may not grow faster than the number of parameters of interest $h$, but $h$ may grow with $T$ at a controlled rate. (ref) therefore allows for variance estimation of an increasing number of estimators. We believe the restrictions on $P$ are reasonable, as they apply to the most commonly performed hypothesis tests in practice, such as joint significance tests (where $\boldsymbol{R}_N$ is the identity matrix), or tests for the equality of parameter pairs.

As a natural implication of (ref), (ref) gives an asymptotic distribution result for a quantity composed exclusively of estimated components.

(ref) allows one to perform a variety of hypothesis tests. For a significance test on a single variable $j$, for instance, take $\boldsymbol{R}_{N}$ as the $j$th basis vector. Then, inference on $\beta^0_j$ of the form $\mathbb{P}\left(\frac{\sqrt{T}(\hat b_j-{{\beta}}_j^0)}{\sqrt{\hat\omega_{j,j}/\hat\tau^4_j}}\leq z\right)-\boldsymbol{\Phi}(z)=o_p(1),\quad \forall z\in\mathds{R}$, can be obtained where $\boldsymbol{\Phi}(\cdot)$ is the standard normal CDF. One can then obtain standard confidence intervals $CI(\alpha):=\left[\hat{b}_j-z_{\alpha/2}\sqrt{\frac{\hat{\omega}_{j,j}/\hat{\tau}_j^4}{T}}, \ \hat{b}_j+z_{\alpha/2}\sqrt{\frac{\hat{\omega}_{j,j}/\hat{\tau}_j^4}{T}}\right]$, where $z_{\alpha/2}:=\boldsymbol{\Phi}^{-1}(1-\alpha/2)$, with the property that $\sup\limits_{{\boldsymbol{\beta}}^0\in\boldsymbol{B}(s_r)}\left\vert\mathbb{P}\left(\beta_j^0\in CI(\alpha)\right)-(1-\alpha)\right\vert=o_p(1)$. For a joint test with $P$ restrictions on $h$ variables of interest of the form $\boldsymbol{R}_{N} \boldsymbol{\beta}^0 = \boldsymbol{q}$, one can construct a Wald type test statistic based on (ref), and compare it to the critical value $F_P^{-1}(1-\alpha)$. Note that these results can also be used to test for nonlinear restrictions of parameters via the Delta method CB2002.

As the bounds and convergence rates as displayed in full generality in (ref) may be hard to interpret, we investigate in (ref) how the conditions of (ref) can be satisfied in a simplified asymptotic setup, thereby illustrating how the different growth rates interact. As for (ref), the conditions on $\lambda$ effectively require that $Q_T$, $N$, and $s_{r,\max}$ grow at a polynomial rate of $T$, which we exploit in (ref) to simplify the conditions.

Inference on high-dimensional parameters

The reason for considering $h\leq C$ in (ref) lies entirely in the application of the central limit theorem. However, while inference on a finite set of parameters covers many cases of interest in practice, it does not allow for simultaneous inference on all parameters. We therefore next consider inference on a growing number of parameters (or hypotheses). We follow the approach pioneered by chernozhukov2013gaussian to consider tests which can be formulated as a maximum over individual tests, and apply a high-dimensional CLT for the maximum of a random vector of increasing length. zhang2017gaussian and zhang2018gaussian provide such a CLT for high-dimensional time series, with serial dependence characterized through the functional dependence framework of Wu05, while chernozhukov2019inference derive a similar result under general $\beta$-mixing conditions. In more recent work, chang2021central derive a high-dimensional CLT for $\alpha$-mixing processes, that we base our result on. Recalling that a process which is NED on an $\alpha$-mixing process can be well-approximated by a mixing process, this mixing condition remains conceptually close to, if more stringent than, our NED framework.\footnote{Ideally one would directly have a high-dimensional CLT available for NED processes, such that it would directly fit to our assumptions. However, such a result is, to our knowledge, currently not available in the literature. While such a result would clearly be very interesting to obtain, this is left for future research given the intricacies needed to derive it.} We therefore build on their results to provide distributional results for high-dimensional inference in (ref). While the core of the proof directly follows by applying the CLT of chang2021central, one still needs to integrate this with the results from (ref) on the consistency of the covariance matrix, as well as adapting the CLT to our estimators. We therefore believe it is worthwhile to state this as a formal result in (ref). Correspondingly, we now strengthen our assumptions as follows.

(ref)(ref) implies (ref)(ref). (ref)(ref) states that the NED process $\boldsymbol{z}_t$ can be well-approximated by an $\alpha$-mixing process; clearly this holds when it is itself $\alpha$-mixing. More specifically, the sequence is NED on itself, such that (ref)(ref) is satisfied for any positive $d$. Furthermore, the exponential decay of the $\alpha$-mixing coefficients is stricter than our restrictions on $\boldsymbol{s}_{T,t}$. Similarly, the sub-gaussian moments in (ref)(ref) imply that all finite moments in (ref)(ref) and (ref) exist, so $m$ may be arbitrarily large.

Unlike (ref), (ref) allows one to simultaneously test a growing number of hypotheses, while controlling for family-wise error rate, for example by the stepdown method described in Section 5 of chernozhukov2013gaussian. One such test is an overall test of significance, with the null hypothesis $\boldsymbol{\beta}^0=\boldsymbol{0}$; in this case $P=h=N$ and $\boldsymbol{R}_N=\boldsymbol{I}$. Note that although $\mathbb{P}\left(\max\limits_{1\leq p\leq P}\hat g_p\leq z\right)$ cannot be calculated analytically, it can easily be approximated with arbitrary accuracy by simulation.

Due to the stronger assumptions in (ref), we can relax the conditions on the growth rates of $N$ and $s_{r,\max}$ compared to (ref) and (ref). In particular, the size of $a$ and $\mathcal{H}$ are not restricted, meaning that $N$ and $h$ can grow at an arbitrarily large polynomial rate of $T$. The conditions on $s_{r,\max}$ can also be relaxed so it can grow up to a rate of $\sqrt{T}$, depending on $r$. This corresponds to our analysis in (ref) when we let $m$ and $d$ tend to infinity.

Analysis of Finite-Sample Performance

We analyze the finite sample performance of the desparsified lasso by means of simulations. We start by discussing tuning parameter selection in (ref). We then discuss three simulation settings: a high-dimensional autoregressive model with exogenous variables (in (ref)), a factor model (in (ref)), and a weakly sparse VAR model (in (ref)). In (ref) and (ref), we compute coverage rates of confidence intervals for single hypothesis tests. In (ref), we perform a multiple hypothesis test for Granger causality.

Tuning parameter selection

While the previous sections give some theoretical restrictions on the tuning parameter choice, these results cannot be used in practice since its value depends on properties of the underlying model that are unobservable. In this section, we provide a feasible recommendation to select the tuning parameters (in both the original regression and nodewise regressions) in a data-driven way.

In particular, we adapt the iterative plug-in procedure (PI) used in, for instance, BCCH12,BCH14,belloni2017program to a time series setting. We build on the theoretical relation between the tuning parameter and the empirical process in (ref), namely the restriction that $\frac{1}{T}\left\lVert\boldsymbol{X}^\prime \boldsymbol{u}\right\rVert_{\infty}\leq C\lambda$ needs to hold with high probability, to guide the choice of $\lambda$. For large $N$ and $T$, $\frac{1}{T}\left\lVert\boldsymbol{X}^\prime \boldsymbol{u}\right\rVert_{\infty}$ can be approximated by the maximum over an $N$-dimensional multivariate Gaussian distribution with covariance matrix $\Omega_{N,T}^{(\mathcal{E})} = \mathbb{E} \left[ \frac{1}{T} \boldsymbol{X}^\prime \boldsymbol{u} \boldsymbol{u}^\prime \boldsymbol{X} \right]$.\footnote{Under minimal extra assumptions (sub-Gaussian moments for $\boldsymbol{x}_t$, and minimum eigenvalue of the long-run covariance matrix bounded away from 0), (ref) substantiates the validity of this approximation.} One may therefore approximate its quantiles by simulating from a multivariate Gaussian with covariance matrix a consistent estimate $\hat{\Omega}^{(\mathcal{E})}$ of $\Omega_{N,T}^{(\mathcal{E})}$.

Our time series setting requires the usage of a consistent long-run variance estimator, which is provided by (ref). We therefore take $\hat{\Omega}^{(\mathcal{E})}$ as in (ref) with $\hat\boldsymbol{\Xi}^{(\mathcal{E})}(l)=\frac{1}{T-l}\sum\limits_{t=l+1}^T \boldsymbol{x}_t\hat u_t\hat u_{t-l}\boldsymbol{x}_{t-l}^\prime$. We set the number of lags in the long-run covariance estimator as the automatic bandwidth estimator in andrews1991heteroskedasticity, specifically $Q_T=\left\lceil1.1447(\hat\alpha(1)T)^{1/3}\right\rceil$, with $\hat\alpha(1)$ computed based on an AR(1) model, as detailed in eq. (6.4) therein. As the estimates $\hat{u}_t$ require a choice of $\lambda$, we iterate the algorithm until the chosen $\lambda$ converges. Full details are provided in Supplementary Appendix C.5. Throughout all simulations, the lasso estimates are obtained through the coordinate descent algorithm GLMnet applied to standardized data.

Autoregressive model with exogenous variables

Inspired by the simulation studies in KockCallot2015 (Experiment B) and MedeirosMendes16, we take the following DGP

where $\boldsymbol{x}_t$ is a $(N-1)\times 1$ vector of exogenous variables. In this simulation design (and the following ones), we consider different values of the time series length $T=\left\lbrace100,200,500,1000\right\rbrace$ and number of regressors $N=\left\lbrace101,201,501,1001\right\rbrace$. For this data generating process, we take $\rho=0.6$, $\beta_j=\frac{1}{\sqrt{s}}(-1)^j$ for $j=1,\dots,s$, and zero otherwise. For $N=101,201$ we set $s=5$ and $s=10$ for $N=501,1001$. The autoregressive parameter matrices $\boldsymbol{A}_1$ and $\boldsymbol{A}_4$ are block-diagonal with each block of dimension $5\times5$. Within each matrix, all blocks are identical with typical elements of 0.15 and -0.1 for $\boldsymbol{A}_1$ and $\boldsymbol{A}_4$ respectively. Due to the misspecification of nodewise regressions, there is induced autocorrelation in the nodewise errors $v_{j,t}$. However, the block diagonal structure of $\boldsymbol{A}_1$ and $\boldsymbol{A}_4$ keeps the sparsity of nodewise regressions constant asymptotically.

We consider different processes for the error terms $u_t$ and $\boldsymbol{\nu}_t$:

For all designs, we evaluate whether the 95% confidence intervals corresponding to $\rho$ and $\beta_1$ cover their true values at the correct rates. The intervals are constructed as $\left[\hat{\rho}\pm z_{0.025}\sqrt{\frac{\hat{\omega}_{1,1}/\hat{\tau}_1^4}{T}}\right]$ and $\left[\hat{\beta}_1\pm z_{0.025}\sqrt{\frac{\hat{\omega}_{2,2}/\hat{\tau}_2^4}{T}}\right]$. These results are obtained based on 2,000 replications. The rates at which the intervals contain the true values are reported in (ref).

We start by discussing the results for the model with Gaussian errors (Model A). Coverage for $\rho$ is close to the nominal level of 95% for all combinations of $N$ and $T$, with some combinations producing slightly conservative results. The coverage rates for $\beta_1$ are worse than for $\rho$. This is likely due to the fact that the exogenous variables $\boldsymbol{x}_t$ within the same block are strongly correlated to each other which negatively impacts the performance of the lasso.

Turning to the results for the model with GARCH errors (Model B), similar finite sample coverage rates are obtained. We do see a small increase in the mean interval width, which is to be expected given the heteroskedastic error structure. With correlated errors (Model C), we again observe consistent coverage rates near the nominal level for $\rho$. Interestingly, the coverage rates for $\beta_1$ appear considerably better than in Models A and B, though in most cases still remaining below the nominal rate at around 90%. We also observe higher mean interval widths than Model A, which is due to larger variance of $\boldsymbol{x}_t$ induced by the cross-sectional covariance of the errors.

In Supplementary Appendix C.6 we provide details on an examination of various selection methods for tuning parameters through heat maps for the coverage levels, which also shed some further light on the relatively poor performance for $\beta_1$ compared to $\rho$ visible for models A and B. In addition to selection by our PI method, we indicate selection by the BIC, the AIC, and the EBIC as in chen2012extended, with $\gamma=1$.\footnote{For additional stability in the high-dimensional settings, we restrict the BIC, AIC, and EBIC to only select models with at most $T/2$ nonzero parameters, though this restriction appears to be binding for the AIC only.} We summarize the main findings below. First, notice that there are regions with coverage close to the nominal level in nearly all scenarios and combinations of $N$ and $T$, suggesting that good coverage could be achieved by selecting the tuning parameters well. Second, across all scenarios, PI generally tends to result in coverage rates closest to the nominal coverage of 95%. As expected, the AIC produces, overall, the least sparse solutions, the EBIC the sparsest and BIC lies in between. PI lies mostly between the BIC and EBIC. Third, there is a region of relatively low coverage for large values of the tuning parameter in the initial and nodewise regressions (see the top right corner of the heat maps). This occurs more pronouncedly for $\beta_1$ than for $\rho$ and especially for $T=1000$. Since PI tends to select near this region, it partly explains why its coverage is worse for $\beta_1$. The relatively better coverage of $\beta_1$ in Model C is matched by this region being much less prominent. Given that the regions of good coverage are in different places for $\rho$ and $\beta_1$, using the BIC or EBIC for generally smaller or larger $\lambda$ would not lead to consistently better coverage across scenarios.\footnote{To confirm this analysis, we also performed the simulations results for all three setups using selection of $\lambda$ by BIC (the best performing information criterion); in line with the heat maps, the coverage rates for BIC are generally somewhat worse than for PI. Results are available upon request.}

Factor model

We take the following factor model

where $\boldsymbol{x}_t$ is a $N\times 1$ vector generated by the AR(1) factor $f_t$. We take $\boldsymbol{\beta}$ as in (ref) with $s$ increased by one to match the number of non-zero parameters. The $N \times 1$ vector of factor loadings $\boldsymbol{\Lambda}$ is chosen with the first $s$ entries (corresponding to the variables with non-zero entries in $\boldsymbol{\beta}$) set to 0.5, and the remaining entries $\Lambda_i=(i-s+1)^{-1}$. This choice of weakly sparse factor loadings ensures that the nodewise regressions are weakly sparse too, as shown in (ref). By letting the large loadings coincide with the non-zero entries in $\boldsymbol{\beta}$, we ensure that there is a large potential for incurring (omitted variable) bias in the estimates, and thus that this DGP provides a serious test for the desparsified lasso.

We investigate whether the confidence interval for $\beta_1$, $\left[\hat{\beta}_1\pm z_{0.025}\sqrt{\frac{\hat{\omega}_{1,1}/\hat{\tau}_2^4}{T}}\right]$, covers the true value at the correct rate. Results are reported in (ref). Coverage rates improve with growing values of $N$ and $T$, with empirical coverages of approximately 85% for small $N$ and $T$, and increasing towards the nominal level when either $N$ or $T$ increases. This result is therefore in line with our theoretical framework, and provides a relevant practical setting in which the desparsified lasso is appropriate to use even if exact sparsity is not present.

Weakly sparse VAR(1)

Inspired by KockCallot2015 (Experiment D), we consider the VAR(1) model

with $\boldsymbol{z}_{t}$ a $(N/2)\times 1$ vector. We focus on testing whether $x_t$ Granger causes $y_t$ by fitting a a VAR(2) model, such that we have a total of $N$ explanatory variables per equation. The $(j,k)$-th element of the autoregressive matrix $A_1^{(j,k)}=(-1)^{\vert j-k\vert}\rho^{\vert j-k\vert+1}$, with $\rho=0.4$. To measure the size of the test, we set $A_1^{(1,2)}=0$; to measure the power of the test, we keep its regular value of $-\rho^2$. Weak sparsity holds\footnote{The weak sparsity measure is $\sum\limits_{j=1}^{N}\vert\rho^{j}\vert^{r}$ with asymptotic limit $\frac{\rho^{r}}{1-\rho^{r}}<\infty$, trivially satisfying $B=0$.} under our choice of the autoregressive parameters, but exact sparsity is violated by having half of the parameters non-zero. Note that the desparsified lasso is convenient for estimating the full VAR equation-by-equation, since all equations share the same regressors, and $\hat{\boldsymbol\Theta}$ needs to be computed only once. For our Granger causality test, however, only a single equation needs to be estimated.

We test whether $x_t$ Granger causes $y_t$ by regressing $y_t$ on the first and second lag of $\boldsymbol{z}_t$. To this end, we test the null hypothesis $A^{(1,2)}_{1}=A^{(1,2)}_2=0$ by using the Wald test statistic in (ref), with $\hat{\boldsymbol{b}}_H=\left(0,\hat{A}^{(1,2)}_{1},0\dots0,\hat{A}^{(1,2)}_{2},0\dots0\right)'$, $H=\left\lbrace2,N/2+1\right\rbrace$, and $\hat{A}^{(1,2)}_{1}$, $\hat{A}^{(1,2)}_{2}$ obtained by regressing $y_t$ on $\left(\boldsymbol{z}_{t-1}',\boldsymbol{z}_{t-2}'\right)'$. We reject the null hypothesis when the statistic exceeds $\chi^2_{2,0.05}\approx5.99$.

We start by discussing the size of the test in (ref). Overall, the empirical sizes exceed the nominal size of 5%, with performance generally not improving for larger sample sizes. In particular, rejection rates slightly deteriorate for larger $N$. However, the observed changes in performance across $N$ and $T$ are rather small and may be due to simulation randomness. The power of the test increases with both $N$ and $T$, reaching 1 at $T=1000$ regardless of the value for $N$.

To improve the finite-sample performance of the method, a natural extension would be to consider the bootstrap for constructing confidence intervals as opposed to asymptotic theory. Bootstrap-based inference for desparsified lasso methods in high dimensions has already been explored by several authors, for example dezeure2017high in the IID setting, and in time series by Krampe18, chernozhukov2019inference and chernozhukov2021timeandspace. In particular, block or block multiplier bootstrap methods, which would allow one to capture serial dependence nonparametrically, would fit our setup well. The block bootstrap has the additional advantage of correcting the finite-sample performance of statistics based on long-run variance estimators, which might be a factor for our tests as well GoncalvesVogelsang11. However, due to the lack of theory about such bootstrap methods, and the associated selection of tuning parameters like the block length, for high-dimensional NED processes, we do not consider such methods here. The development of such theory would be a highly relevant and interesting topic for future research.

Conclusion

We provide a complete set of tools for uniformly valid inference in high-dimensional stationary time series settings, where the number of regressors $N$ can possibly grow at a faster rate than the time dimension $T$. Our main results include (i) an error bound for the lasso under a weak sparsity assumption on the parameter vector, thereby establishing parameter and prediction consistency; (ii) the asymptotic normality of the desparsified lasso under a general set of conditions, leading to uniformly valid inference for finite subsets of parameters; (iii) asymptotic normality of a maximum-type statistic of a growing, high-dimensional, number of tests, valid under more stringent conditions, thereby also permitting simultaneous inference over a potentially large number of parameters, and (iv) a consistent Bartlett kernel Newey-West long-run covariance estimator to conduct inference in practice.

These results are established under very general conditions, thereby allowing for typical settings encountered in many econometric applications where the errors may be non-Gaussian, autocorrelated, heteroskedastic and weakly dependent. Crucially, this allows for certain types of misspecified time series models, such as omitted lags in an AR model.

Through a small simulation study, we examine the finite sample performance of the desparsified lasso in popular types of time series models. We perform both single and joint hypothesis tests and examine the desparsified lasso's robustness to, amongst others, regressors and error terms exhibiting serial dependence and conditional heteroskedasticity, and a violation of the sparsity assumption in the nodewise regressions. Overall our results show that good coverage rates are obtained even when $N$ and $T$ increase jointly. The factor model design shows that the desparsified lasso remains applicable when the exact sparsity assumption of the nodewise regressions is violated. Finally, Granger causality tests in the VAR are slightly oversized, but empirical sizes generally remain close to the nominal sizes, and the test's power increases with both $N$ and $T$.

There are several extensions to our approach that are interesting to consider. The development of a high-dimensional central limit theorem for NED processes would allow to weaken the dependence conditions needed for establishing simultaneous, high-dimensional inference. Similarly, using sample splitting would likely allow for weakening sparsity assumptions. Finally, improvements in finite sample performance may be achieved by bootstrap procedures. All of these extensions would require the development of novel theory, and thus provide challenging but worthwhile avenues for future research.

Acknowledgements

We thank the editor, associate editor and three referees for their thorough review and highly appreciate their constructive comments which substantially improved the quality of the manuscript.

The first and second author were financially supported by the Netherlands Organization for Scientific Research (NWO) under grant number 452-17-010. The third author was supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sk\lodowska-Curie grant agreement No 832671. Previous versions of this paper were presented at CFE-CM Statistics 2019, NESG 2020, Bernoulli-IMS One World Symposium 2020, (EC)$^2$ 2020, and the 2021 Maastricht Workshop on Dimensionality Reduction and Inference in High-Dimensional Time Series. We gratefully acknowledge the comments by participants at these conferences. In addition, we thank Etienne Wijler for helpful discussions. All remaining errors are our own.

\numberwithin{lemma}{section} \numberwithin{equation}{section} \numberwithin{example}{section} \numberwithin{definition}{section}