Sparse High-Dimensional Vector Autoregressive Bootstrap

Sparse High-Dimensional Vector Autoregressive Bootstrap

\onehalfspacing

\doublespacing

Introduction

We introduce theory for bootstrapping the distribution of high-dimensional means of sparse, finite order, stable vector autoregressive (VAR) processes. For an $N$-dimensional vector of time series $\boldsymbol{x}_t=(x_{1,t},\dots,x_{N,t})^\prime$, we provide an approximation for the distribution of $\max\limits_{1\leq j\leq N}\left\lvert\frac{1}{\sqrt{T}}\sum\limits_{t=1}^T x_{j,t}\right\rvert$, where the number of variables $N$ is potentially much larger than the sample size $T$, and can asymptotically grow faster than $T$. This prototypical statistic is commonly considered in high-dimensional settings, see e.g. the closely related work of CCK13, CCK17, zhang2017gaussian, chernozhukov2020nearly, giessing2020bootstrapping, or the review by chernozhukov2022high, who investigate the properties of this estimator for independent data. In this paper, we extend these results to high-dimensional linear processes, including stable VARs. Related work in time series settings include zhang2018gaussian, who provide Gaussian approximations in the general framework of functional dependence of Wu05.

The VAR sieve bootstrap is well-known in the low-dimensional time series bootstrapping literature, see e.g. paparoditis1996, Park2002, ChangPark2003, MeyerKreiss15, and Section 12.2 of kilian2017structural. It fits a VAR to the time series data, resamples the residuals of the estimated VAR, and re-applies the VAR recursively to place the dependence back into the bootstrap sample. Under appropriate conditions, the VAR sieve bootstrap allows for valid inference. We extend this approach to high dimensions where the VAR is estimated by the lasso (Tibshirani96) or another sparse estimation method, and use a multiplier (or wild) bootstrap to resample the residuals. Our work is related to that of Trapani2013, bi2021ar and Krampe19. The two former papers assume a dense structure on the data, and apply the VAR sieve bootstrap to a low-dimensional set of factors. The latter consider a sparse setting, providing bootstrap inference for desparsified estimators of VAR coefficients. We assume a data-generating process (DGP) similar to the one considered in Krampe19.

All theoretical results in this paper are established under two different sets of assumptions on the errors. First, we assume the errors have sub-gaussian moments, which generally allows $N$ to grow at an exponential rate of $T$. Second, we assume that the errors have some finite number of absolute moments, which effectively restricts the growth of $N$ to some polynomial rate of $T$. In (ref), we introduce the multiplier bootstrap for sparsely estimated high-dimensional VARs. In (ref), we start by providing a high-dimensional central limit theorem (HDCLT) for linear processes in (ref), which may be of independent interest. In (ref), we introduce the stable VAR model, and show that under consistent estimation, the long run covariance structure is recovered with high probability. (ref) provides a consistency result for the covariance matrix. In (ref), we show that the bootstrap's behaviour is asymptotically similar to that of the original sample. In particular, (ref) provides a HDCLT for the bootstrap process which mirrors that of (ref), and (ref) shows consistency of the bootstrap. (ref) then shows how these results can be used to establish validity of inference in VARs estimated by the lasso.

Notation. For a random variable $x$, $\left\lVertx\right\rVert_{L_p}=\left(\mathbb{E}\left\lvertx\right\rvert^p\right)^{1/p}$, $\left\lVertx\right\rVert_{\psi_2}=\inf\left\lbrace c>0:\mathbb{E}\exp(\left\lvertx\right\rvert^2/c^2)\leq 2\right\rbrace$ denote the $L_p$ and Orlicz norms. For any $N$ dimensional vector $\boldsymbol{x}$, $\left\Vert \boldsymbol{x}\right\Vert_p=\left(\sum\limits_{j=1}^{N}\left\vert x_j\right\vert^p\right)^{1/p}$ denotes the $p$-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_p = \max_{\left\lVert\boldsymbol{x}\right\rVert_p = 1} \left\lVert\boldsymbol{A} \boldsymbol{x}\right\rVert_p$ for any $p \in [0, \infty]$ and $\left\lVert\boldsymbol{A}\right\rVert_{\max}=\max\limits_{i,j}\left\vert a_{i,j}\right\vert$. $\Lambda_{\min}(\boldsymbol{A})$ and $\Lambda_{\max}(\boldsymbol{A})$ denote the smallest and largest eigenvalues of $\boldsymbol{A}$, and $\rho(\boldsymbol{A})$ the spectral radius of $\boldsymbol{A}$, i.e. the largest absolute eigenvalue of $\boldsymbol{A}$, or equivalently $\rho(\boldsymbol{A})=\lim\limits_{k\to\infty}\left\lVert\boldsymbol{A}^k\right\rVert^{1/k}$ for any induced norm $\left\lVert\cdot\right\rVert$. For $\boldsymbol{A}$ a square matrix, we let its zero-th power $\boldsymbol{A}^0=\boldsymbol{I}$. 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_{i,t}$). When they are used, it should be understood that there exists some constant $C$ or sequence $\eta_T\to0$ such that the given statement holds.

Vector Autoregressive Bootstrap

We introduce our proposed bootstrap procedure for sparsely estimated high-dimensional VARs and subsequently discuss how it can be used to perform inference on high-dimensional time series.

Bootstrap for High-Dimensional VARs

Let $\boldsymbol{x}_t$ be an $N$-dimensional time series process. We assume the data is generated by a stable, finite order, high-dimensional VAR($K$) model

with autoregressive parameter matrices $\boldsymbol{A}_k \ (k=1,\ldots,K)$, independent errors $\boldsymbol{\epsilon}_t $ with $\mathbb{E}\boldsymbol{\epsilon}_t=\boldsymbol{0}$ and covariance matrix $\boldsymbol{\Sigma}_{\epsilon} :=\frac{1}{T}\sum\limits_{t=1}^T\mathbb{E}\boldsymbol{\epsilon}_t\boldsymbol{\epsilon}_t^\prime$, and $\boldsymbol{x}_t=\boldsymbol{\epsilon}_t=\boldsymbol{0}$ for $t<1$. We can re-write (ref) as a collection of linear equations

where $\boldsymbol{a}_{j,k}$ is the $j$th row of $\boldsymbol{A}_k$, $\boldsymbol{\beta}_{j}=(\boldsymbol{a}_{j,1},\dots,\boldsymbol{a}_{j,K})^\prime$, and $\mathcal{X}_t=(\boldsymbol{x}_{t-1}^\prime,\dots,\boldsymbol{x}_{t-K}^\prime)^\prime$. We denote data stacked into a matrix as $\underset{T\times N}{\boldsymbol{X}}=(\boldsymbol{x}_{1}^\prime, \dots, \boldsymbol{x}_{T}^\prime)^\prime$. The lasso estimator of equation $j$ is defined as

where $\lambda_j$ is a tuning parameter that determines the degree of penalization in equation $j$, and can be selected independently in each equation. For tuning parameter selection, one could use e.g. the theoretically founded method of kock2024data, the iterative plug-in procedure described in Section 5.1 of adamek2021lasso, or information criteria.

Once all equations $j=1,\ldots, N$ are estimated by the lasso, we collect the VAR coefficient estimates as follows

Our object of interest is the scaled high-dimensional mean $$ Q= \max\limits_{1\leq j\leq N}\left\lvert\frac{1}{\sqrt{T}}\sum\limits_{t=1}^{T} x_{j,t}\right\rvert$$ of the sparse VAR. To approximate its distribution, we apply the VAR multiplier bootstrap summarized in Algorithm (ref). When $B$ is sufficiently large, the CDF of $Q$ can be approximated by the quantiles of the ordered statistics $Q^{*(1)}, \dots, Q^{*(B)}$. Note that while we derive results for the maximum absolute mean, this bootstrap procedure is equally valid for statistics such as $\max\limits_{1\leq j\leq N}\frac{1}{\sqrt{T}}\sum\limits_{t=1}^Tx_{j,t}$ or $\min\limits_{1\leq j\leq N}\frac{1}{\sqrt{T}}\sum\limits_{t=1}^Tx_{j,t}$, which would allow for one-sided tests, or tests with an asymmetric rejection region.

Bootstrap Inference on (Approximate) Means

Statistics such as the scaled mean $Q$ are useful in high-dimensional settings, since they allow us to simultaneously test a high-dimensional set of hypotheses. For example, let $\mu_{j}=\mathbb{E} x_{j,t}$ be the means of a high-dimensional stable autoregressive process, and assume we are interested in testing the hypothesis

Under the null hypothesis, this process follows (ref), which allows us to directly test the null using the quantiles of $Q^{*(1)}, \dots, Q^{*(B)}$. Specifically, one would reject the null at significance level $\alpha$ if $Q>Q^{*(B[1-\alpha])}$. To know for which means the null can be rejected, one can use the stepdown procedure of romano2005exact, as detailed in the description in Section 5 of CCK13. Importantly, this procedure is asymptotically exact -- non-conservative -- as it takes into account the possible correlations between statistics, instead of using the conservative worst case of independence.

More generally, this bootstrap procedure can be used to test any high-dimensional set of hypotheses, provided its test statistic can be expressed as an approximate mean, that is, $\frac{1}{\sqrt{T}}\sum\limits_{t=1}^Tx_{j,t}+o_p(1)$. While we do not formally consider this extension here, we can adapt the arguments in Section 5 of CCK13 (which do not rely on independent data) to establish this result in our context as well. This opens up the way for applications to statistics that are much more general than just sample means, as many statistics of practical interest, such as (high-dimensional) regression estimates, can be written in this form. Our results therefore form a first step towards a more general bootstrap theory for high-dimensional inference using VAR models on statistics that can be well-approximated by the mean of a linear process.

HDCLT for Linear Processes

In this section, we establish a high-dimensional CLT for linear processes, which is a useful result in its own right, but also a vital building block to establish theoretical results for the bootstrap. We therefore give it a self-contained treatment in this section, before applying it to the VAR process in (ref) and covering the theory for the bootstrap in the following sections.

Under appropriate invertibility conditions, it is well-known that the VAR process in (ref) can be written in the following infinite order vector moving average (VMA) form

where $\mathcal{B}(z)=\sum\limits_{k=0}^{\infty}\boldsymbol{B}_k z^k=\left(\boldsymbol{I}-\sum\limits_{k=1}^K\boldsymbol{A}_k z^k\right)^{-1}$, and $L$ is the lag operator. We derive a Gaussian approximation for linear processes of the form in (ref), which builds on and extends similar approximations for independent and identically distributed (i.i.d.) processes by chernozhukov2020nearly and others (see (ref)).

Specifically, we show that the distribution of $\max\limits_{1\leq j\leq N}\left\lvert\frac{1}{\sqrt{T}}\sum\limits_{t=1}^T x_{j,t}\right\rvert=\left\lVert\frac{1}{\sqrt{T}}\sum\limits_{t=1}^T\boldsymbol{x}_t\right\rVert_{\infty}$ can be asymptotically approximated by $\left\lVert\boldsymbol{z}\right\rVert_{\infty}$, with $\boldsymbol{z}\sim N(\boldsymbol{0},\boldsymbol{\Sigma})$ and $\boldsymbol{\Sigma}$ an appropriate covariance matrix. This result parallels well-known results in low-dimensional settings, where scaled means of linear processes converge in distribution to a Gaussian random variable as $T\to\infty$. However, in our high-dimensional setting, we consider the case where $N$ and $T$ diverge simultaneously, and $\left\lVert\frac{1}{\sqrt{T}}\sum\limits_{t=1}^T\boldsymbol{x}_t\right\rVert_{\infty}$ does not converge to a well defined limit; the maximum over a growing number of elements generally also grows. As such, we instead show that their distributions grow closer together asymptotically, in the sense that the Kolmogorov distance between between $\left\lVert\frac{1}{\sqrt{T}}\sum\limits_{t=1}^T\boldsymbol{x}_t\right\rVert_{\infty}$ and $\left\lVert\boldsymbol{z}\right\rVert_\infty$ converges to 0. Even though to our knowledge, there does not exist a closed-form expression for the CDF of $\left\lVert\boldsymbol{z}\right\rVert_{\infty}$, it can be approximated for any $N$ by Monte Carlo simulation, making it a useful asymptotic approximation in practice.

The broad sketch of our proof is as follows. We use the Beveridge-Nelson decomposition to write

where $\tilde{\mathcal{B}}(z)=\sum\limits_{j=0}^\infty\sum\limits_{k=j+1}^\infty\boldsymbol{B}_k z^j$. The first term is a scaled sum of independent errors with covariance matrix $\boldsymbol{\Sigma}:=\mathcal{B}(1)\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}\mathcal{B}(1)^\prime$, $\sigma^2_{j}:=\boldsymbol{\Sigma}_{(j,j)}$, and can therefore be approximated by a Gaussian maximum thanks to chernozhukov2020nearly when $\boldsymbol{\Sigma}$ is non-degenerate and the $\boldsymbol{\epsilon}_t$'s satisfy certain moment conditions (see (ref)). The second term is an asymptotically negligible leftover under certain summability conditions on the VMA coefficient matrices $\boldsymbol{B}_k$ (see (ref)). Formally, we make the following assumptions:

We derive our results under two different moment assumptions. In (ref).(ref) we require that the errors are uniformly sub-gaussian over $j$ and $t$; or in (ref).(ref) that the moments possess some number ($m$) of finite absolute moments. By equation (2.15) in vershynin2019high, (ref).(ref) follows automatically for all $m$ from (ref).(ref), making the latter a considerably less stringent assumption. Under these assumptions, (ref) provides an upper bound on the Kolmogorov distance between our statistic of interest and a Gaussian maximum:

Under (ref).(ref), convergence of this upper bound to 0 depends on the size of the terms $\tilde{S}$ and $S_2$, and the relative growth rates of $N$ and $T$. As $N$ only enters in logs compared to $\sqrt{T}$ in the denominator, it is possible to have $N$ grow at some exponential rate of $T$. Under (ref).(ref), $N$ enters the numerator at a polynomial rate through the sequence $d_{N}$; this effectively restricts the growth rate of $N$ to some polynomial of $T$, though it can still grow faster than $T$ when $m$ is sufficiently large. Our results under these two sets of assumptions therefore mainly differ (apart from the different proof strategies required for each case), in this regard: if exponential growth of $N$ is desirable, we need finite exponential moments of $\boldsymbol{\epsilon}_t$; whereas if polynomial growth of $N$ is sufficient, we only need finite polynomial moments of $\boldsymbol{\epsilon}_t$.

Application to VAR Models

(ref) is a key building block in our derivations for the bootstrap, as it can be applied to our VAR in (ref) under appropriate conditions. In this section, we explain our assumptions on the VAR process, and on the consistency properties of lasso estimation. While the lasso is our running example, the following theoretical results do not rely on the lasso specifically, and are equally valid for any other estimation method which satisfies our consistency conditions. We return to the lasso in (ref), where we show examples of it satisfying these conditions.

For the following exposition, it is useful to define the companion matrix

of the VAR in (ref). This matrix allows us to re-write the VAR($K$) as a VAR(1) with

and allows for a simple expression for the corresponding VMA coefficients in (ref): $\boldsymbol{B}_k=\boldsymbol{J}\mathds{A}^k\boldsymbol{J}^\prime$, where $\underset{N\times KN}{\boldsymbol{J}}=\left(\boldsymbol{I},\boldsymbol{0},\dots,\boldsymbol{0}\right)$.\footnote{See page 279 of paparoditis1996.} This inversion is only possible if the VAR is invertible.

(ref) is based on Assumption 1(ii) of Krampe19, and its purpose is twofold. First, it allows us to derive summability properties for the quantities $\tilde{S}$ and $S_q$ in (ref), since $\left\lVert\boldsymbol{B}_j\right\rVert_{\infty}\leq \left\lVert\mathds{A}^j\right\rVert_{\infty}\leq \psi_N\theta^j$. Second, it implies that the VAR process in (ref) is stable, since $\rho(\mathds{A})=\lim\limits_{k\to\infty}\left\lVert\mathds{A}^k\right\rVert_{\infty}^{1/k}\leq\lim\limits_{k\to\infty}\left(\psi_N\theta^k\right)^{1/k}=\theta$, and it can therefore be inverted into a VMA. Based on this inequality, it is also clear that when $k$ is large, $\left\lVert\mathds{A}^k\right\rVert_{\infty}\approx \rho(\mathds{A})^k\leq \theta^k$, i.e., the powers of $\mathds{A}$ will eventually converge at an approximately exponential rate. The magnitude of $\psi_N$ controls the magnitude of $\left\lVert\mathds{A}\right\rVert_\infty$, which may be substantially larger than 1 even in VAR models with low persistence. The growth rate of $\psi_N$ controls how quickly $\left\lVert\mathds{A}^k\right\rVert_\infty$ approaches $\theta^k$, as the dimension of $\mathds{A}$ increases. Sequences of VAR models which require $\psi_N$ to grow were (to our knowledge) first highlighted in liu2021bernstein, who relate the growth of $\psi_N$ to spatial dependence, as opposed to temporal dependence tied to $\theta$.

While our results allow for DGPs with $\psi_N$ growing, it should be noted that such DGPs suffer in terms of convergence rates required for bootstrap validity, and many are already implicitly excluded by (ref). To illustrate this, consider a VAR(1) where $\left\lVert\mathds{A}^j\right\rVert_{\infty}$ grows with $N$. In many cases this leads to $\mathcal{B}(1)=\sum\limits_{j=0}^{\infty} \mathds{A}^j$ growing with $N$ as well, resulting in $\sigma_j^2$ growing. However, this is not always the case, and (ref) in Appendix (ref) shows a DPG which satisfies (ref) while requiring $\psi_N$ to grow exponentially with $N$.

Next, we make the following assumptions about consistency of the estimators $\hat\mathds{A}$, and the residuals $\hat\boldsymbol{\epsilon}_t$:

While we leave the sequences $\xi_{N,T}$ and $\phi_{N,T}$ unspecified and derive later results in terms of these sequences, the reader may think of them as $\xi_{N,T}$ converging at a rate close to $\frac{1}{\sqrt{T}}$ and $\phi_{N,T}$ close to $\frac{1}{T}$ for reasonable estimators. Regarding the assumption that $\psi_N\xi_{N,T}\leq \bar{C}(1-\theta)^2$, a sufficient condition to satisfy this is that $\psi_N\xi_{N,T}\to0$ and $N,T$ are sufficiently large. However, this formulation highlights that our requirements on $\xi_{N,T}$ -- and therefore on the estimation error $\left\lVert\hat{\mathds{A}}-\mathds{A}\right\rVert_{\infty}$ -- are stricter for VARs with large temporal and/or spatial dependence. We elaborate more on these rates when using the lasso in (ref).

In our proof strategy, we make use of the probabilistic sets denoted by calligraphic letters $\mathcal{P}$ to $\mathcal{U}$. They describe events involving functions of the random variables $\boldsymbol{x}_t$ and $\boldsymbol{\epsilon}_t$, and can therefore only hold with a certain probability. For the sets $\mathcal{P}$ and $\mathcal{Q}$, we assume that they hold with probability converging to 1 as $N,T\to \infty$. For the other sets, they are chosen in such a way that we can show they hold with probability converging to 1 under our assumptions. For example, relevant to this section are the sets

and

The different subscripts of these sets indicate for which version of (ref) they are intended. We show they hold with high probability in (ref). Note that many of our intermediate results are phrased as non-random bounds on random quantities, which hold on these sets, i.e., these bounds hold with probability 1 conditionally on these random events occurring. For the main result in (ref), we then show that the probability of all these random events occurring jointly converges to 1, such that these non-random bounds hold asymptotically.

The main result of this section concerns the consistency of our estimate of $\boldsymbol{\Sigma}$, namely $\hat\boldsymbol{\Sigma}:=\hat\mathcal{B}(1)\hat\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}\hat\mathcal{B}(1)^\prime$, with $\hat\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}:=\frac{1}{T}\sum\limits_{t=1}^T\hat\boldsymbol{\epsilon}_t\hat\boldsymbol{\epsilon}_{t}^\prime$, $\hat\mathcal{B}(z)=\boldsymbol{I}+\sum\limits_{k=1}^{\infty}\hat\boldsymbol{B}_k z^{k}$, $\hat\mathcal{B}(z)=\boldsymbol{I}+\sum\limits_{k=1}^{\infty}\hat\boldsymbol{B}_k z^{k}$. Unsurprisingly, the form of $\hat\boldsymbol{\Sigma}$ mirrors that of $\boldsymbol{\Sigma}$, since we apply the same Beveridge-Nelson decomposition in (ref) to the bootstrap process. To do so, the estimated VAR is required to be invertible, i.e. $\rho(\hat\mathds{A})<1$; we show that this is the case with probability converging to 1 in (ref).4. This justifies our suggested invertibility correction in (ref), since it is asymptotically negligible. In finite samples one can perform this correction by, for example, checking if $\rho(\hat\mathds{A})>0.999$, and if so, multiplying each element of $\hat\mathds{A}$ by $0.999/\rho(\hat\mathds{A})$. In (ref) we establish a covariance closeness result which plays a crucial role in showing consistency of our proposed bootstrap method in the next section.

Bootstrap Consistency

In this section, we introduce some of the bootstrap-related notation, and flesh out the exact properties of the processes $\boldsymbol{x}_t^*$ and $\boldsymbol{\epsilon}_t^*$. In (ref), we then give a Gaussian approximation for the bootstrap process, mirroring (ref). Finally, (ref) provides the main result of bootstrap consistency.

As is customary in the bootstrap literature, we define the following bootstrap conditional notation: Let $\mathbb{P}^*\left(\cdot\right)$ denote the bootstrap probability conditional on the sample $\boldsymbol{X}$, and $\mathbb{E}^*\left(\cdot\right)$ the expectation with respect to $\mathbb{P}^*$, and similarly ley $\left\lVertx\right\rVert_{\psi_2}^*:=\inf\left\lbrace c>0:\mathbb{E}^*\exp(\left\lvertx\right\rvert^2/c^2)\leq 2\right\rbrace$ and $\left\lVertx\right\rVert_{L_p}^*:=\left(\mathbb{E}^*\left\lvertx\right\rvert^p\right)^{1/p}$ denote the corresponding conditional norms. We let

and $\boldsymbol{x}_t^*$ built from $\boldsymbol{\epsilon}_t^*$

where $\boldsymbol{A}_{k}^*:=\hat\boldsymbol{A}_k$. By construction, the bootstrap processes $\boldsymbol{x}_t^*$ and $\boldsymbol{\epsilon}_t^*$ then follow a VAR process mirroring (ref), and can be inverted under appropriate conditions to a VMA process mirroring (ref): $\hat\boldsymbol{B}_k=\boldsymbol{J}\hat\mathds{A}^k\boldsymbol{J}^\prime$, where $\underset{N\times KN}{\boldsymbol{J}}=\left(\boldsymbol{I},\boldsymbol{0},\dots,\boldsymbol{0}\right)$. This then also leads to the bootstrap versions of $\tilde{S}$ and $S_q$, and the following bootstrap equivalent of (ref).

Since $\boldsymbol{z}$ in (ref) is the same as in (ref), we can combine both theorems and a telescopic sum argument to bound the distance between distributions of $\left\lVert\frac{1}{\sqrt{T}}\sum_{t=1}^T\boldsymbol{x}_t\right\rVert_{\infty}$ and $\left\lVert\frac{1}{\sqrt{T}}\sum_{t=1}^T\boldsymbol{x}_t^*\right\rVert_{\infty}$, giving us bootstrap consistency in the following theorem.

Bootstrap Consistency for VAR Estimation by the Lasso

The application of our proposed bootstrap method requires that the lasso satisfies (ref) with sequences $\psi_{N}$, $\xi_{N,T}$, and $\phi_{N,T}$ such that the bound in (ref) converges to 0. In this section, we show that this is the case under both options of (ref), and under both weak and exact row-wise sparsity of the underlying VAR.

As described in (ref) we propose to estimate the VAR equation-by-equation, using the lasso estimators in (ref). Our goal is therefore to find bounds on $\max\limits_{j}\left\lVert\hat\boldsymbol{\beta}_j-\boldsymbol{\beta}_j\right\rVert_1$ and $\max\limits_{j}\frac{1}{T}\left\lVert\hat\boldsymbol{\epsilon}_{j}-\boldsymbol{\epsilon}_j\right\rVert_2^2=\max\limits_{j}\frac{1}{T}\sum\limits_{t=1}^T\left[(\hat\boldsymbol{\beta}_j-\boldsymbol{\beta}_j)^\prime\mathcal{X}_t\right]^2$. For this purpose, we will be using error bounds in Corollary 1 of our previous work in adamek2021lasso, though similar error bounds have been derived in different contexts by other authors; see e.g. bickel2009simultaneous, KockCallot2015, MedeirosMendes16, and masini2019regularized. Next, we will elaborate on the assumptions under which these error bounds hold.

For Assumption 1 of adamek2021lasso, we have $\mathbb{E}\boldsymbol{x}_t=0\implies \mathbb{E}\mathcal{X}_t=0$ by the structure of (ref), and $\mathbb{E}\boldsymbol{x}_t\epsilon_{j,t}=0,~\forall j$, by independence of the errors. We then need to assume that $\max\limits_{j,t}\mathbb{E}\left\lvertx_{j,t}\right\rvert^m\leq C$ in addition to (ref).(ref) in this paper to ensure the first part of the assumption is satisfied. This high-level assumption on moments of $x_{j,t}$ can also be shown to hold under more primitive conditions, such as a moment condition on linear combinations of the errors, $\max\limits_{\left\lVert\boldsymbol{u}\right\rVert_2\leq 1, t}\mathbb{E}\left\lvert\boldsymbol{u}^\prime\boldsymbol{\epsilon}_{t}\right\rvert^m\leq C$, and a new summability condition on the rows of $\boldsymbol{B}_k$, $\max\limits_{j}\sum\limits_{k=0}^\infty\left\lVert\boldsymbol{b}_{j,k}\right\rVert_2^m\leq C$:

Note that $m$ in this paper corresponds to $2\bar{m}$ in adamek2021lasso. Under an additional assumption that $\psi_N\leq C$,\footnote{This additional assumption is in line with e.g. kock2024data who require this in their Assumption 2.(2) to obtain error bounds on the lasso.} (ref) ensures that the NED assumption is satisfied uniformly across equations and as $N$ grows. The VMA coefficients decay at an exponential rate, therefore satisfying any polynomial decay rate on the NED sequence, and the assumption is satisfied for any arbitrarily large $d$. Assumption 2 of adamek2021lasso requires that the rows of $\mathds{A}$ are weakly sparse, in the sense that $\left\lVert\boldsymbol{\beta}_j\right\rVert_r^r=\left\lVert[\mathds{A}]_{j,\cdot}\right\rVert_r^r\leq s_{r,j}$ for some $0\leq r<1$. Assumption 3 of adamek2021lasso requires that the covariance matrix of the regressors satisfies a form of compatibility condition; for simplicity, we can assume that $\Lambda_{\min}\left(\frac{1}{T}\sum\limits_{t=1}^T\mathbb{E}\mathcal{X}_{t}\mathcal{X}_{t}^\prime\right)$ is bounded away from zero, which is sufficient to satisfy the condition simultaneously for all equations. For an example of conditions when this is satisfied, see Equation 6 of masini2019regularized. Under these conditions, we have by Corollary 1 of adamek2021lasso that

with probability converging to 1 under appropriate restrictions on the $\lambda_j$, detailed in Theorem 1 of adamek2021lasso. Note that these restrictions are a function of the dependence (NED size $d$) and sparsity ($s_{r,j}$) within each equation, so in order to satisfy (ref), these properties should hold uniformly across equations.

To further simplify this result, we can use the asymptotic setup of Example C.1 of adamek2021lasso where $N$, $\lambda_j$, and $s_{r,j}$ grow at a polynomial rate of T. While that example provides the full details on the tradeoff between $r$, the number of moments, and the growth rates of $s_{r,j}$ and $N$ relative to $T$, here we fix $r=1/2$ and $s_{r,j}\sim T^{1/8},~\forall j$ for illustrative purposes.

While (ref) shows an example of conditions for bootstrap consistency using the finite absolute moments in (ref).(ref), the stronger assumption of sub-gaussian moments in (ref).(ref) allows for faster growth of $N$ relative to $T$. In this scenario, we can consider the error bounds in Theorem 2 of KockCallot2015,

with $\lambda_j=C\ell_{T}^{5/2}\ell_{N}^{2}\ell_{K}\ell_{N^2K}^{1/2}\sigma_T^2/\sqrt{T}$. Note that $\sigma_T^2$ denotes the largest variance among all $\epsilon_{j,t}$ and $x_{j,t}$, so we once again make the high level assumption that $\max\limits_{j,t}\mathbb{E} x_{j,t}^2\leq C$. To obtain these bounds, we need the additional assumption that the errors are Gaussian, so $\boldsymbol{\epsilon}_t\overset{\text{iid}}{\sim} N(\boldsymbol{0},\boldsymbol{\Sigma}_{\epsilon})$, which implies (ref).(ref). Additionally, they consider the case of exact sparsity, with $\sum\limits_{k=1}^{KN}\mathds{1}_{\{\left\lvert\left[\mathds{A}\right]_{j,k}\right\rvert>0\}}\leq s_{0,j}$. Finally, $\kappa_{j}$ play a similar role to the compatibility constant in Assumption 2 of adamek2021lasso, and are bounded away from 0 when $\Lambda_{\min}\left(\frac{1}{T}\sum\limits_{t=1}^T\mathbb{E}\mathcal{X}_{t}\mathcal{X}_{t}^\prime\right)\geq 1/C$, see the discussion on page 7 of KockCallot2015 for details. Regarding the growth rates of $N$ and $s_{0,j}$, we take a similar example to Theorem 3 of KockCallot2015, with $N\sim e^{(T^a)}$ and $s_{0,j}\leq C T^b$.

Simulations

To evaluate the finite sample performance of our proposed method, our simulation study covers a variety of DGPs on which we compare size and power with other bootstrap methods typically used in a high-dimensional time series setting.

Setup

We implement our proposed VAR multiplier bootstrap with two different ways of selecting the lasso penalty. First, we estimate the VAR with the penalty chosen by the Bayesian information criterion jointly over all equations (VAR-BIC). Second, we use the theoretically founded data-driven method of kock2024data (VAR-TF). For both methods the number of lags $K$ is chosen as the informative upper bound in Section 5 of HecqMargaritellaSmeekes2023, as mentioned in (ref). For details, see Algorithm (ref) in Appendix C. Additionally, we leave the diagonal elements of the VAR coefficient matrices unpenalized in the lasso estimation. We believe this is good common practice with lasso VAR estimation, because a series' own lags are often more important than those of other series for explaining the dynamic properties. This approach is similar to the “Own-Other” hierarchical penalties in Wilms2020Hlag or the Minnesota prior in Bayesian VAR estimation. To guarantee stability of the estimated VAR, we apply the finite sample correction mentioned in (ref): If $\rho(\hat\mathds{A})>0.999$, we multiply each element of $\hat\mathds{A}$ by $0.999/\rho(\hat\mathds{A})$.

As a benchmark, we also show results for the `oracle' method, which does no VAR estimation, and generates bootstrap samples using the true VAR coefficients (VAR-oracle).

In addition to the VAR-based bootstrap, we consider two block-based bootstrap methods: the block wild/multiplier bootstrap (BWB) based on e.g. Shao2011 or ZhangCheng2014bootstrapping, and the moving block bootstrap based on e.g. PALM2011 or Smeekes2015 (MBB). For both block-based bootstraps, we use a block length using the automatic bandwidth estimator for the Bartlett kernel in andrews1991heteroskedasticity.

We study four DGPs used by other work in this field. Specifically, we take inspiration from KockCallot2015, Krampe19, Barigozzi2024FNETS. In all DGPs, we consider every combination of $T\in\left\lbrace 50,100,200,500\right\rbrace$, and $N\in\left\lbrace 20,40,100,200 \right\rbrace$. To estimate size, we generate the data with population mean $0$ for each variable. The nominal level is $\alpha=0.05$, and for better readability, all size plots are truncated at a rejection rate of 0.5. For power, we add a nonzero constant $\mu$ to a proportion $p$ of variables, such that the first $Np$ variables have mean $\mu$ and the remaining $N(1-p)$ variables have mean $0$. We consider $p=0.5$ for all DGPs, and choose $\mu$ separately for each DGP according to an initial calibration exercise, such that the power is relatively low (around 25%) for $N=20$, $T=50$. In DGP1, we also investigate the effects on power of increasing $p$ to $0.9$, and doubling $\mu$.

DGP1: Diagonal VAR(1)

This DGP is based on Experiment A of KockCallot2015:

where $\boldsymbol{A}=\text{diag}(0.5,\dots,0.5)$ and $\boldsymbol{\Sigma}_\epsilon=\text{diag}(0.01,\dots,0.01)$. This DGP satisfies (ref) with $\Lambda_{\min}\left(\boldsymbol{\Sigma}\right)=\max\limits_{1\leq j\leq N}\sigma_j^2=0.04$ for all $N$, (ref).(ref) with Gaussian errors and (ref) with $\theta=0.5$, $\psi_N=1$. This DGP is the “best-case” setup for our proposed method because the lasso generally performs well in sparse models, and all the true non-zero parameters in this DGP are left unpenalized.

Regarding the size in the top row of (ref), we generally see the VAR-based methods achieve correct, slightly conservative size. With the exception of $N=100,~T=100$, VAR-BIC and VAR-TF perform very similarly, being slightly more conservative than the oracle method. They are generally more conservative at larger $N$, but improve and reach close to nominal size as $T$ increases. At $N=100,~T=100$, BIC tends to select a very low value of the tuning parameter, often reaching the lower edge of the grid. This results in models with almost no regularization, excessive variance, and poor performance of VAR-BIC. This phenomenon is also observed in later DGPs, so this seems to be a somewhat pervasive issue with BIC. Both block-based bootstrap methods have comparable performance, reaching size between 5 and 15%. This large size is most pronounced at low $N$, though we see improvement with growing $T$. At $N=200$, both methods exceed 5% only slightly, with the BWB outperforming the MBB.

Power is given in the bottom three rows of (ref). We see similar patterns across all three settings: For all methods, power grows considerably with $T$, and slightly with $N$, and reaches close to 100% at $N=200,~T=500$. The VAR-based methods have slightly lower power than the oracle method, and the block-based methods beat the oracle. This is not necessarily an indictment against the VAR-based methods, as the block-based methods do not achieve size control. The abnormal behavior of the BIC is also reflected in the power, reaching 100% at $N=100,~T=100$. Comparing between the three settings, we see that increasing the nonzero proportion from $p=0.5$ to $p=0.9$ increases the power only slightly, by around 5-15 percentage points. Doubling the mean from $\mu=0.0175$ to $\mu=0.035$ had a much larger impact, more than doubling the power in most cases. This is not a surprising pattern, given that the test statistic is based on the maximum of means.

DGP2: Block-diagonal VAR(1)

DGP2 is based on Example 1 of Krampe19. It follows (ref) with $\boldsymbol{A}$ and $\boldsymbol{\Sigma}_\epsilon$ having a block-diagonal structure. The blocks are $20\times20$ in both cases; their precise definition\footnote{We use the $\xi=0.6$ version of this DGP.} can be found in Appendix D of Krampe19, and we provide a visual overview of the pattern within blocks in (ref) in Appendix C.\footnote{Note that we do not shuffle the indices of variables like Krampe19.} This DGP satisfies (ref) with $\Lambda_{\min}\left(\boldsymbol{\Sigma}\right) \approx 0.0782$ and $\max\limits_{1\leq j\leq N}\sigma_j^2 \approx 38.322$ for all $N$ (in multiples of 20), (ref).(ref) with Gaussian errors and (ref) with $\theta=0.8$, $\psi_N\approx 3.121$. We expect our proposed method to perform well in this DGP: most of the structure within the blocks is on the unpenalized diagonal, and is quite sparse even in the last 6 rows.

In terms of size (top row of (ref)), all methods other than the oracle are generally oversized. Between the VAR-based methods, the VAR-TF has better size than VAR-BIC except at $T=500$ and $N=100,200$, where VAR-TF performs the worst with around $15\%$ size. Except the $T=500$ case, VAR-TF has the best size, with performance close to the oracle at $N=100,200$. For low $T$, VAR-BIC's performance changes significantly over different $N$, with size around $15\%$ at $N=20$, but well below $5\%$ at $N=200$.

Given that the oracle method has the correct size, the relatively poor performance of the VAR-based methods is largely due to estimation. Estimation is challenging in this DGP because of high persistency with $\rho(\mathds{A})=0.8$, as VAR estimates can be heavily biased in such cases, even when using least squares estimation. A classic solution to this issue in low-dimensional settings is the double bootstrap of Killian1998doublebootstrap; it is an interesting avenue of future research to investigate whether the results would improve using a similar approach in our setting. The block-based methods both have similar performance, with size around $15-20\%$ at $T=50$, and reaching $5-10\%$ at $T=500$. The high persistence of this DGP also hampers the block-based methods, since they need long blocks to accurately capture the dependence.

Regarding the power results displayed in the bottom row of (ref), we see large improvements with growing $T$, and changes over $N$ are broadly in line with the changes in rejection rates seen in the size plots.

DGP3: Weakly sparse VAR(1)

This DGP is based on Experiment D of KockCallot2015. It follows (ref) with $\boldsymbol{A}$ having a Toeplitz structure and exponentially decaying off-diagonals, $a_{ij}=(-1)^{\left\lverti-j\right\rvert} \rho^{\left\lverti-j\right\rvert+1}$, $\rho=0.3$. $\boldsymbol{\Sigma}_\epsilon$ is the same as in DGP1. For (ref), $\boldsymbol{\Sigma}$ changes as $N$ grows, but its properties stabilize at $\Lambda_{\min}\left(\boldsymbol{\Sigma}\right) \approx 0.0234$ and $\max\limits_{1\leq j\leq N}\sigma_j^2 \approx 0.0142$. (ref) is satisfied with $\theta=0.6$ and $\psi_N=1$. While this DGP is not sparse in the exact sense, it is weakly sparse with elements far from the diagonal taking values very close to zero. The lasso will inevitably set most parameters equal to zero, but we do not expect this to have a large impact on performance, since the effect of these near-zero parameters on the dynamic properties of the process is negligible.

We see a similar pattern in the size (top row of (ref)) as for DGP1: the VAR-based methods perform similarly, being slightly conservative, except a few cases where VAR-BIC fails. The block-based methods are oversized again, with size around $10\%$.

For power (bottom row of (ref)) the pattern is also similar to DGP1: power generally increases greatly with $T$, and slightly with $N$, and the relative power of different methods is in line with the differences in size.

DGP4: Factor model with sparse idiosyncratic component

This DGP is based on the simulation setup (E1)+(C2) in Appendix E.1 of Barigozzi2024FNETS:

where the entries of $\boldsymbol{\lambda}_{i,\ell}\in\mathds{R}^2$ are generated as i.i.d. standard Gaussian, $\boldsymbol{D}=\boldsymbol{D}_0\cdot 0.7 /\Lambda_{\max}(\boldsymbol{D}_0)$, where $\boldsymbol{D}\in\mathds{R}^{2\times 2}$ has off-diagonal elements generated i.i.d. from $U[0,0.3]$ and diagonal elements generated from $U[0.5,0.8]$. The $w_i$ are such that the sample estimate of $\text{Var}(\chi_{i,t})/\text{Var}(\xi_{i,t})=1,~\forall i$. To generate $\boldsymbol{A}$, first $\boldsymbol{A}_0$ is generated, with its entries drawn i.i.d. from $Bernoulli(1/N)\cdot 0.275$. Then, if $\Lambda_{\max}(\boldsymbol{A}_0)\leq 0.9$, $\boldsymbol{A}=\boldsymbol{A}_0$; otherwise $\boldsymbol{A}=\boldsymbol{A}_0\cdot 0.9/\Lambda_{\max}(\boldsymbol{A}_0)$. This DGP does not fit the VAR structure in (ref), and Assumptions (ref)-(ref) do not hold. The process is stationary, but if a VAR representation exists, it is likely not sparse due to the factor structure. We expect our proposed method to perform more poorly relative to the block-based bootstrap methods, since it is an adverse setting for the lasso. Note that since the DGP is not a VAR model, the oracle method is not implemented for this DGP.

Contrary to our expectations, size results in the top row of (ref), demonstrate good performance of the VAR-based methods, especially compared to the block-based methods. They are slightly oversized at around $10\%$ for $T=50$, but are close to nominal for larger $T$. On the other hand, the block-based methods are oversized across the board, with size at $20\%$ at $T=50$ and only decreasing to $10\%$ at $T=500$. Power in the bottom row of (ref) shows improvements with increasing $T$ as for the other DGPs, and not much change with $N$. The relative powers of the different methods is broadly in line with the size differences.

Conclusion

In this paper, we introduce a VAR multiplier bootstrap procedure which approximates the distribution of scaled high-dimensional means, using the lasso to estimate the VAR. We motivate the usefulness of this procedure as a tool for inference in high-dimensional time series, allowing for non-conservative simultaneous testing of a large set of hypotheses. We show that the bootstrap is consistent under two different moment assumptions on the errors: sub-gaussian moments, and a finite number of absolute moments. Under the former, $N$ can grow at an exponential rate of $T$. Under the latter, $N$ can only grow at a polynomial rate of $T$, with the growth rate of $N$ limited by the number of absolute moments available.

We provide guidance for estimating the VAR bootstrap model by the lasso as a running example. We show that the lasso satisfies appropriate error bounds for consistency of the bootstrap distribution, under the assumption that the underlying VAR process is (row-wise) sparse. In our examples, we derive explicit limits on the growth rate of $N$ relative to $T$ thereby allowing for exact and weak sparsity of the VAR.

To establish the consistency of the VAR multiplier bootstrap, we derive a Gaussian approximation for the maximum mean of a linear process, which may be of independent interest. Our results can be applied to more complex statistics than simple means, and we believe that extending this method to inference for linear model coefficients is an interesting avenue for future research. Our simulation results show generally good performance of the lasso-VAR-based bootstrap, with the exception of highly persistent DGPs. We believe that another interesting extension would be a bias-corrected version of the bootstrap to improve performance in highly persistent DGPs.

\numberwithin{lemma}{section} \numberwithin{equation}{section}