Data-Driven Tuning Parameter Selection for High-Dimensional Vector Autoregressions

Data-Driven Tuning Parameter Selection for High-Dimensional Vector Autoregressions

\noindentKeywords: High-dimensional time series, vector autoregressive model, $\ell_1$-penalized estimation, data-driven tuning parameter selection.

Introduction

Multivariate time series play a fundamental role in many areas of research. The quintessential approach to modelling these is based on the linear VAR model; see, e.g., lutkepohl_new_2005. In order to fully capture the dynamics of the time series at hand and to decrease the risk of omitting variables, one frequently includes many explanatory variables. These considerations result in a high-dimensional VAR model wherein the number of parameters can exceed the number of observations. Consequently, the parameters can no longer be estimated by least squares, and as an alternative there has been a surge of research on Lasso-type estimators [tibshirani_regression_1996]. For instance, assuming independent and identically distributed (i.i.d.) Gaussian innovations in the VAR, consistency and oracle inequalities for the Lasso have been studied by han_transition_2013, basu_regularized_2015, kock_oracle_2015 and davis_sparse_2016, among others. The assumption of Gaussian innovations has been relaxed in several papers, including song_large_2011, wong_lasso_2020 and miao_var_2022. In particular, wong_lasso_2020 derive consistency results for the Lasso in VAR models under suitable mixing and moment conditions, whereas miao_var_2022 derive rates of convergence and oracle properties for factor augmented VAR models, relying on results for weakly dependent processes [in the sense of wu_nonlinear_2005]. masini_VAR_2022 establish oracle inequalities in the case of martingale difference innovations that are not necessarily mixing. chernozhukov_lasso_2021 consider Lasso-driven inference for time series and spatial models, while adamek_lasso_2022 consider inference based on the desparsified Lasso when the data-generating process allows near-epoch dependence. Further related papers are han_direct_2015, guo_var_2016 and wu_performance_2016.

To implement Lasso-type estimators one must choose the penalization level $\lambda$, and the theoretical results listed above are derived under a suitable choice of this level. However, the “right” choice of this tuning parameter typically depends on unknown population quantities such as mixing coefficients or other coefficients quantifying the dependence structure of the data generating process {[}wong_lasso_2020, babii_timeseries_2022, masini_VAR_2022{]}, the population covariance matrix of the observed variables and innovations {[}kock_oracle_2015, MEDEIROS_highdim_2016{]}, the population coefficient matrix {[}basu_regularized_2015{]}, or other quantities depending on the data generating process [chernozhukov_lasso_2021, adamek_lasso_2022, miao_var_2022]. Consequently, despite the large amount of research on penalized estimation of VAR models, there is currently no firm guidance on how to choose $\lambda$ in practice. Instead, one often resorts to information criteria [e.g. kock_oracle_2015, MEDEIROS_highdim_2016 and masini_VAR_2022], cross-validation [e.g. wong_lasso_2020, babii_timeseries_2022 and miao_var_2022] or other methods without providing theoretical guarantees for these. Crucially, such choices of penalization do not necessarily satisfy the conditions imposed in the theoretical results. Therefore, strictly speaking, the theoretical guarantees provided are only valid for fortuitous choices of $\lambda$. In this paper we resolve this problem by proposing a data-driven way of choosing $\lambda$ along with prediction and estimation error guarantees for the resulting weighted Lasso estimator of the parameters in large VAR models. As upper bounds on the estimation error play a crucial role in establishing the validity of inference based on debiasing in high-dimensional models [see javanmard_confidence_2014, van_de_geer_asymptotically_2014 and zhang_confidence_2014], our results open the door for inference in VAR models with a data-driven choice of $\lambda$.

The penalization algorithm we study originates from belloni_sparse_2012, who consider regressions in high dimensions with independent data. There are several challenges in adapting this algorithm to time series data: First, the analysis of belloni_sparse_2012 relies on the independence as well as certain high-level conditions on the explanatory variables and their relation to the model errors. In VAR models one cannot impose such conditions as their validity is completely determined by the model (which also generates the explanatory variables). Instead, imposing only primitive conditions, we carefully take into account the inherent dependence in the VAR. Second, as heavy tails are omnipresent in many time series, we allow for the possibility of certain types of heavy-tailed innovation distributions, that is, so-called sub-Weibull innovations. To accommodate dependence and heavy tails, we establish a novel maximal inequality for centered sums of dependent sub-Weibull random variables, that may be of independent interest. The validity of the penalty loadings (i.e. Lasso weights) proposed in the algorithm is then established by proving that these are close to certain infeasible ideal loadings that are constructed by means of blocking-based self-normalization. In order to verify that the ideal loadings are well-behaved, we rely on recent moderate deviation theory for self-normalized block-sums of weakly dependent processes [in the wu_nonlinear_2005 functional dependence sense] derived by chen_self-normalized_2016 and gao_refined_2022.

To alleviate the shrinkage bias introduced by the Lasso, one often refits the parameters of variables selected by the Lasso using least squares. The performance of the resulting post-Lasso depends crucially on the Lasso variable selection. As this selection depends on the chosen level of penalization, the post-Lasso itself depends on this tuning parameter. We show that the post-Lasso following the weighted Lasso implied by our data-driven tuning parameter choice obeys the same performance guarantees as the weighted Lasso. Note, however, that these guarantees are established under an additional assumption (not needed for the weighted Lasso) on certain sparse eigenvalues of the population regressor covariance matrix.

Finally, the square-root Lasso (henceforth: sqrt-Lasso) of belloni2011square is a popular alternative to the Lasso for i.i.d. data in the absence of conditional heteroskedasticity. Apart from the work of sqrtJTSA, who study a linear model with fixed regressors but dependent error terms, the theoretical performance of the sqrt-Lasso has not been studied for time series and, in particular, it is not clear how to choose its tuning parameter. We resolve this problem by providing a fully data-driven implementation of the sqrt-Lasso with performance guarantees matching those of the other estimators studied in the absence of conditional heteroskedasticity.

Outline

In Section (ref) we present the weighted Lasso estimator and describe the data-driven tuning parameter selection. Section (ref) presents the assumptions and performance guarantees for the weighted Lasso. Sections (ref) and (ref) cover the post- and sqrt-Lasso, respectively. Sections (ref) and (ref) contain simulations and an empirical illustration, respectively. Proofs, technical lemmas, and additional simulation output are contained in the supplementary appendices.

Notation

For $k\in\mathbb N$, we write $[k]:=\{1,\dotsc,k\}$. For $\boldsymbol{\delta}\in\mathbb R^{k}$, we denote its $\ell_{r}$-norm, by $\norm{\boldsymbol{\delta}}_{\ell_{r}}:=(\sum_{j=1}^k|\delta_j|^r)^{1/r},r\in[1,\infty)$, and $\norm{\boldsymbol{\delta}}_{\ell_{\infty}}:=\max_{j\in[k]}|\delta_j|$. For $\mathcal{S}\subseteq [k]$ non-empty, $\boldsymbol{\delta}_\mathcal{S}\in\mathbb R^{|\mathcal{S}|}$ is the subvector of $\boldsymbol{\delta}$ picked out by by $\mathcal{S}$. When applied to a real matrix $\mathbf{A}$, the aforementioned norms are understood as the induced (operator) norms. We use $\mathbf{A}_{k:k+l,m:m+n}$ to denote the submatrix picked out by rows $k$ through $k+l$ and columns $m$ through $m+n$. For a square matrix $\mathbf{A}\in\mathbb{R}^{k\times k}$, we denote its spectrum (set of eigenvalues) by $\Lambda(\mathbf{A})$ and write $\rho(\mathbf{A}):=\max\{|\lambda|:\lambda\in\Lambda(\mathbf{A})\}$ for its spectral radius (with $\left|\cdot\right|$ the complex modulus). If $\mathbf{A}\in\mathbb R^{k\times k}$ is symmetric, we write $\Lambda_{\min}(\mathbf{A})$ and $\Lambda_{\max}(\mathbf{A})$ for the minimum and maximum eigenvalues, respectively.

For a random scalar $X$, we denote its $L_r$-norm by $\norm{X}_{r}:=(\mathrm{E}[|X|^{r}])^{1/r}$, $r\in[1,\infty)$, with $\mathrm{E}[\cdot]$ denoting the expectation operator. For $\alpha\in(0,\infty)$, we define the sub-Weibull($\alpha$) norm of $X$ as $\norm{X}_{\psi_{\alpha}}:=\sup_{r\in[1,\infty)}r^{-1/\alpha}\norm{X}_{r}$. A random scalar $X$ is said to be sub-Weibull$(\alpha)$ if $\norm{X}_{\psi_{\alpha}}<\infty$.\footnote{One may show that the space of sub-Weibull$(\alpha)$ random variables is complete with respect to $\left\Vert \cdot \right\Vert _{\psi_{\alpha}}$. Consequently, it holds that the norm $\left\Vert \cdot \right\Vert _{\psi_{\alpha}}$ is countably sub-additive, which we will use repeatedly.} For a $k$-dimensional random vector $\boldsymbol{X}$, we define its joint sub-Weibull($\alpha$) norm $\Vert \boldsymbol{X}\Vert _{\psi_{\alpha}}:=\sup\{\Vert \boldsymbol{u}^{\top}\boldsymbol{X}\Vert _{\psi_{\alpha}}:\left\Vert \boldsymbol{u}\right\Vert _{\ell_{2}}=1\}$ and call $\boldsymbol{X}$ jointly sub-Weibull($\alpha$) if $\left\Vert \boldsymbol{X}\right\Vert _{\psi_{\alpha}}<\infty$.

For non-random numbers $a_n$ and positive numbers $b_n,n\in\mathbb N,$ we write $a_n=o(1)$ if $a_n\to0$ as $n\to\infty$, and $a_n\lesssim b_n,$ if the sequence $a_n/b_n$ is bounded by a constant. For random variables $V_n$ and positive numbers $b_n,$ we write $V_n\lesssim_{\mathrm P} b_n,$ if the sequence $V_n/b_n$ is bounded in probability. We take $n\geqslant 3$ and $p\geqslant 2$ throughout the manuscript and reserve the word “constant” for non-random quantities that do not depend on $n$.

Model and Penalization Algorithm

We study the $q$\textsuperscript{th}-order $(q\in\mathbb N)$ VAR model given by

yielding a stochastic process $\{\boldsymbol{Y}_t\}_{t\in\mathbb Z}$ which is (strictly) stationary under the assumptions in Section (ref). Here $\boldsymbol{Y}_{t}:=(Y_{t,1},\dotsc,Y_{t,p})^{\top}$ is a random vector of length $p$, $\{\boldsymbol{\Theta}_{0j}\}_{j=1}^q$ are $p\times p$ (unknown) coefficient matrices, and with $\boldsymbol{\varepsilon}_{t}:=(\varepsilon_{t,1},\dotsc,\varepsilon_{t,p})^{\top}$, $\{\boldsymbol{\varepsilon}_t\}_{t\in\mathbb Z}$ is a sequence of innovations. Given observations $\{\boldsymbol{Y}_{t}\}_{t=-(q-1)}^{n}$ from ((ref)), the objective is to estimate $\{\boldsymbol{\Theta}_{0j}\}_{j=1}^q$, while allowing (but not requiring) that the number of elements $p^{2}q$ in $\{\boldsymbol{\Theta}_{0j}\}_{j=1}^q$ is larger than the (effective) sample size $n$. Throughout we take $p\geqslant2$ and $n\geqslant3$.\footnote{We consider the VAR model in ((ref)), allowing the number of output variables $p$, the (common) distribution of $\boldsymbol{\varepsilon}_{t}$, and, hence, that of $\boldsymbol{Y}_{t}$ to depend on the sample size $n$. That is, we consider an array $\{\{\boldsymbol{Y}_{t}^{(n)}\}_{t\in\mathbb{Z}}\}_{n\in\mathbb N}$ of stochastic processes, each process $\{\boldsymbol{Y}_{t}^{(n)}\}_{t\in\mathbb{Z}}$ presumed strictly stationary, in which each $\boldsymbol{Y}_{t}^{(n)}=(Y_{t,1}^{(n)},\dotsc,Y_{t,p_{n}}^{(n)})^\top,t\in\mathbb Z$, is a random element of $\mathbb R^{p_{n}}$. To ease notation, we henceforth suppress the $n$ superscript.}

The process in ((ref)) may be written in companion form,

for \[ \underbrace{\boldsymbol{Z}_{t}}_{pq\times1}:=

, \quad \underbrace{\widetilde{\boldsymbol{\Theta}}_{0}}_{pq\times pq}:=

, \quad\underbrace{\widetilde{\boldsymbol{\varepsilon}}_{t}}_{pq\times1}:=

, \] and $\mathbf{I}_p\in\mathbb R^{p\times p}$ being the identity matrix. In the special case of $q=1$, we interpret $\widetilde{\boldsymbol{\Theta}}_{0}$ as $\boldsymbol{\Theta}_{01}$ and $\widetilde{\boldsymbol{\varepsilon}}_{t}$ as $\boldsymbol{\varepsilon}_{t}$. For $i\in[p]$, the $i$\textsuperscript{th} row in ((ref)) is given by

where the vector $\boldsymbol{\beta}_{0i}^{\top}$ of length $pq$ is the $i$\textsuperscript{th} row of the companion matrix $\widetilde{\boldsymbol{\Theta}}_{0}$. In a low-dimensional setting where $\boldsymbol{\beta}_{0i}^{\top}$ is of fixed length (and short), estimation can be done by equationwise least squares with the same $n\times pq$ regressor matrix $\mathbf{X}:=[\boldsymbol{Z}_{0}:\cdots:\boldsymbol{Z}_{n-1}]^{\top}$ for all $i\in[p]$. However, when $pq>n$, the $pq\times pq$ Gram matrix $\mathbf{X}^{\top}\mathbf{X}$ has reduced rank, and Lasso-type estimators have been studied as an alternative under various sparsity assumptions. Nevertheless, as discussed in the introduction, the theoretical guarantees hitherto established for these are valid only for specific tuning parameter choices depending on unknown population quantities. As a result---despite the sensitivity of shrinkage estimators to the tuning parameter choice---one must currently resort to choosing this parameter by methods without theoretical guarantees when implementing shrinkage estimators in time series models. We provide a solution to this problem in the context of a weighted Lasso by proposing an algorithm for tuning parameter selection and explicitly incorporate the data-driven tuning parameter choice into our theoretical guarantees.

A weighted Lasso estimator satisfies

where $\widehat Q_i$ is the sample average squared error loss function

for $i\in [p]$, $\lambda\in(0,\infty)$ is a penalty level and each $\widehat{\boldsymbol{\Upsilon}}_{i}:=\mathrm{diag}(\widehat{\upsilon}_{i,1},\dotsc,\widehat{\upsilon}_{i,pq})$ is a diagonal matrix of data-dependent penalty loadings $\widehat\upsilon_{i,j}\in(0,\infty),j\in[pq]$. Both the penalty level and loadings will be specified in Algorithm (ref) below. Although the minimization problem in ((ref)) may have multiple solutions, we sometimes refer to such a $\widehat{\boldsymbol{\beta}}_{i}$ as the Lasso estimator. The results established below apply to any such (measurable) minimizer.\footnote{By the measurable selection theorem 6.7.22 in pfanzagl, $\widehat{\boldsymbol{\beta}}_{i}$ in (ref) can be chosen measurable (provided the $\widehat{\boldsymbol{\Upsilon}}_{i}$ are chosen measurable).}

Note that much of the literature on the Lasso and its variants employs the same penalty loadings for all regressors (e.g. $\widehat{\upsilon}_{i,j}\equiv1$). Equal weighting implicitly treats the regressors as either being on the same scale or having been brought onto the same scale by some preliminary transformation of the data, which is then typically abstracted from in the theoretical analysis. Taking serious the effect of such a preliminary transformation is not trivial as it can alter the already intricate dependence structure of the time series at hand. By incorporating data-dependent loadings, our algorithm is not subject to this caveat.

Adapting belloni_sparse_2012 for independent data to our setting, we consider the following data-driven penalization. Let

and let $K\in\mathbb{N}_{0}$ denote a fixed number of loading updates.\footnote{Although Theorem (ref) is valid for all choices of $c>1$, our simulations (Section (ref)) indicate that choices just above $1$ lead to the best performance in practice. Furthermore, as long as $c$ is close to $1$, the performance of Algorithm (ref) is not sensitive to the exact choice. This is documented in Figure (ref) in Section (ref) of the supplementary appendix, which also shows that the same remarks on the choice of $c$ are valid for the post- and sqrt-Lasso. $c=1.1$ was suggested in belloni_sparse_2012.} Following belloni_sparse_2012, we fix $K$ at 15 in our simulations and empirical illustration. To state Algorithm (ref) below, let $\Phi$ be the standard Gaussian cumulative distribution function.

Note that Algorithm (ref) does not require any knowledge of the degree of dependence in the process $\{\boldsymbol{Y}_{t}\}_{t\in\mathbb Z}$ (as quantified by mixing coefficients or other unknown population quantities). This feature of the penalization is in contrast to current recommendations in the literature.

To discuss error rates for the implied Lasso estimator, let $\left\Vert \cdot\right\Vert _{2,n}$ be the prediction norm \[ \Vert \boldsymbol{\delta}\Vert_{2,n}:=\sqrt{\frac{1}{n}\sum_{t=1}^{n}\left(\boldsymbol{Z}^{\top}_{t-1}\boldsymbol{\delta}\right)^{2}} \] of the vector $\boldsymbol{\delta}\in\mathbb R^{pq}$, and let

denote the sparsity number, defined as the largest cardinality of the support of $\boldsymbol{\beta}_{0i}$ across $i\in[p]$. (Without loss of generality, we take $s\geqslant 1$.) Then, under the primitive conditions given in Section (ref), any Lasso estimators $\widehat\boldsymbol{\beta}_i:=\widehat\boldsymbol{\beta}_{i}(\lambda_n^{\ast},\widehat\boldsymbol{\Upsilon}_i^{(K)}),i\in[p]$, in ((ref)) based on the penalty level $\lambda_n^{\ast}$ in ((ref)) and the (final) penalty loadings $\widehat\boldsymbol{\Upsilon}_i^{(K)}$ arising from Algorithm (ref) satisfy

The formal statement, including a set of sufficient conditions, is given in Theorem (ref) below. Importantly, the performance guarantees in (ref)--(ref) match those currently available based on infeasible penalty level choices.

Assumptions and Performance Guarantees for Algorithm (ref)

For the remainder of the paper we invoke (a subset of) the following assumptions.

Assumption (ref).(ref) implies that the innovations $\boldsymbol{\varepsilon}_t$ are causal and strictly stationary but allows them to be $\overline{q}$-dependent. Here $\overline{q}=0$ corresponds to independent errors, as imposed in, e.g., kock_oracle_2015, wong_lasso_2020 and miao_var_2022. Assumption (ref).(ref) states that $\boldsymbol{\varepsilon}_{0}$ belongs to the sub-Weibull family of distributions with tail parameter (at least) $\alpha$, as imposed in wong_lasso_2020 and masini_VAR_2022. Special cases are the sub-exponential $(\alpha=1)$ and sub-Gaussian $(\alpha=2)$ families. If $\boldsymbol{\varepsilon}_{0}$ is Gaussian---as imposed in, e.g., han_transition_2013, basu_regularized_2015, kock_oracle_2015 and davis_sparse_2016---then Assumption (ref).(ref) holds if, in addition, $\Lambda_{\max}(\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}})$ is bounded from above, with $\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}}$ being the covariance matrix of $\boldsymbol{\varepsilon}_{0}$. Assumption (ref).(ref) parallels part of masini_VAR_2022 and miao_var_2022.

Assumption (ref).(ref) is standard and combined with Assumption (ref) ensures that the process $\{\boldsymbol{Y}_{t}\}_{t\in\mathbb Z}$ in ((ref)) [as well as $\{\boldsymbol{Z}_{t}\}_{t\in\mathbb Z}$ in ((ref))] is strictly stationary. In particular, $\{\boldsymbol{Y}_{t}\}_{t\in\mathbb Z}$ has the moving average representation \[ \boldsymbol{Y}_{t}=\sum_{\ell=0}^{\infty}(\widetilde{\boldsymbol{\Theta}}_{0}^{\ell})_{1:p,1:p}\boldsymbol{\varepsilon}_{t-\ell}, \] where the series converges absolutely almost surely (and in sub-Weibull norm).

Assumption (ref).(ref) is milder than masini_VAR_2022, which imposes an exponential decay on the $\ell_1$-norms of the rows of $(\widetilde{\boldsymbol{\Theta}}_{0}^{h})_{1:p,1:p}$ as $h$ grows (as opposed to our decay placed on the $\ell_2$-norms). Likewise, the assumption is also milder than miao_var_2022, which imposes an exponential decay on the spectral norm of $(\widetilde{\boldsymbol{\Theta}}_{0}^{h})_{1:p,1:p}$ (exceeding the row-wise $\ell_2$-norms). \footnote{See \Citet[Lemma 4]{masini_VAR_2022} for further conditions sufficient for Assumption (ref).(ref).}

Assumption (ref).(ref) states that the lag-order $q$ is constant. This appears to be reasonable for most practical purposes as $q$ is often thought of as small relative to $n$.\footnote{This point was also made by kock_oracle_2015, although their analysis allows for growing $q$.} Indeed, in our application in Section (ref) we have $n=758$ and as a result of monthly sampling $q=12$ often suffices. We use that $q$ does not depend on $n$ to exploit recent moderate deviation results in gao_refined_2022 for self-normalized sums based on geometric moment contracting processes. Assumption (ref) therefore seems hard to relax with currently available methods.

Assumptions (ref) and (ref) suffice to show that the population covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{Z}}:=\mathrm{E}[\boldsymbol{Z}_{0}\boldsymbol{Z}_{0}^\top]$ of $\boldsymbol{Z}_0$ exists (in $\mathbb R^{pq\times pq}$).\footnote{See Lemma (ref) in the appendix for details.} We further assume:

Assumption (ref) ensures that $\boldsymbol{\Sigma}_{\boldsymbol{Z}}$ has full rank and is implied by, e.g., miao_var_2022. As demonstrated in basu_regularized_2015, Assumption (ref) follows from additional (boundedness) conditions on $\widetilde{\boldsymbol{\Theta}}_{0}$, although these are not necessary.\footnote{For the purpose of Theorem (ref), we note that Assumption (ref) is used in the proof of Lemma (ref). Inspection of its proof reveals that one can replace Assumption (ref) by a milder requirement that only certain restricted eigenvalues of $\boldsymbol{\Sigma}_{\boldsymbol{Z}}$ are bounded away from zero. Restricted eigenvalue conditions appear in, e.g., kock_oracle_2015 and MEDEIROS_highdim_2016.}

Assumption (ref) imposes sparsity on each row $\boldsymbol{\beta}_{0i}^\top$ of the horizontally concatenated parameter matrices $[\boldsymbol{\Theta}_{01}:\cdots:\boldsymbol{\Theta}_{0q}]$ as in, e.g., kock_oracle_2015, chernozhukov_lasso_2021, and masini_VAR_2022. To help digest the growth and sparsity conditions in Assumption (ref), note that in the special case of sub-Gaussian innovations $(\alpha=2)$ and exponential decay in the $\ell_2$-norms of the rows of $(\widetilde{\boldsymbol{\Theta}}_{0}^{h})_{1:p,1:p}$ as $h$ grows $(\tau=1)$, the exponents equal $C(\alpha,\tau)=6.66\dotsc$ and $1/\tau+\widetilde{C}(2, \alpha)=3$.

We can now state the properties of the weighted Lasso estimator with data-driven tuning parameter selection based on Algorithm (ref).

The estimation and prediction error rates guaranteed by Theorem (ref) match those currently available in the literature building on infeasible penalty parameter choices. Algorithm (ref) is thus the first data-driven tuning parameter selection method for which theoretical guarantees have been provided in the context of dependent data. Upper bounds on the estimation error are also crucial ingredients in establishing valid inference based on the debiasing/desparsification methodology originating in javanmard_confidence_2014, van_de_geer_asymptotically_2014 and zhang_confidence_2014. This methodology has been used in the context of time series in adamek_lasso_2022. Theorem (ref) therefore opens the door to inference in high-dimensional VAR models with data-driven penalty parameter choice.

In establishing Theorem (ref), we cannot call upon the large existing body of maximal and concentration inequalities for independent random variables. Instead, we establish a new maximal inequality for sums of dependent random variables (Lemma (ref) in the appendix). This result is then used to show that the penalty loadings of Algorithm (ref) are close to certain ideal, yet infeasible, blocking-based loadings (Lemma (ref)). Note that although these ideal loadings involve blocking, Algorithm (ref) does not, meaning that no choice of block size is needed in practice. The ideal loadings can, in turn, be used as self-normalizing factors in an application of the recent moderate deviation theorems of gao_refined_2022 for dependent random variables. Using these results, we show that $\lambda_n^*$ is a high-probability upper bound on the maximum of such self-normalized sums. Section (ref) in the appendix explains further steps and challenges involved in establishing Theorem (ref).

Performance Guarantees for the Post-Lasso

To alleviate shrinkage bias, the post-Lasso $\widetilde{\boldsymbol{\beta}}_{i}(\lambda,\widehat{\boldsymbol{\Upsilon}}_{i})$, say, uses least squares to refit the coefficients of the variables $\widehat{T}_i=\widehat{T}_i\del[1]{\widehat{\boldsymbol{\beta}}_i(\lambda,\widehat{\boldsymbol{\Upsilon}}_i)}:=\operatorname{supp}\del[1]{\widehat{\boldsymbol{\beta}}_i(\lambda,\widehat{\boldsymbol{\Upsilon}}_i)}$ selected by a Lasso with generic tuning parameter $\lambda$ (and generic penalty loadings $\widehat{\boldsymbol{\Upsilon}}_i$) in each equation. Here $\operatorname{supp}(\bm{b})=\cbr[0]{j\in[pq]:b_j\neq 0}$ denotes the support of $\bm{b}\in\mathbb R^{pq}$. Therefore, the not necessarily unique, post-Lasso satisfies\footnote{The post-Lasso $\widetilde{\boldsymbol{\beta}}_{i}$ need not be unique since:

Nevertheless, just like Theorem (ref), Theorem (ref) below applies to any post-Lasso estimator and guarantees uniqueness with probability tending to one.}

All coefficients classified as zero by the Lasso are also set to zero by the post-Lasso, i.e. $\widetilde{\boldsymbol{\beta}}_{i,\widehat{T}_i^c}=\bm{0}_{|\widehat{T}_i^c|}$. In particular, $\widetilde{\boldsymbol{\beta}}_{i}=\bm{0}_{pq}$ if $\widehat{T}_i=\emptyset$.

Since the set of variables selected by the “first stage” Lasso $\widehat{T}_i=\widehat{T}_i\del[1]{\widehat{\boldsymbol{\beta}}_i(\lambda,\widehat{\boldsymbol{\Upsilon}}_i)}$ depends on the choice of $\lambda$, the performance of the “second stage” post-Lasso also depends on this tuning parameter. Thus, to obtain a data-driven implementation of the post-Lasso with theoretical performance guarantees, one needs a data driven tuning parameter selection for the first stage Lasso with performance guarantees. For the penalty choice $\lambda_n^*$ and loadings $\widehat{\boldsymbol{\Upsilon}}_i^{(K)}$ stemming from Algorithm (ref) we have provided such guarantees for the first stage Lasso $\widehat{\boldsymbol{\beta}}_i(\lambda_n^*,\widehat{\boldsymbol{\Upsilon}}_i^{(K)})$ in Theorem (ref) for any fixed $K\in\mathbb N_0$.

Writing $\mathbf{X}_{I}\in\mathbb R^{n\times|I|}$ for the matrix consisting of all columns of $\mathbf{X}$ with column index in $I\subseteq [pq]$ and letting $\bm{y}_i=(Y_{1,i},\hdots,Y_{n,i})^\top$, the post-Lasso estimator of the coefficients of the variables retained by the Lasso can be uniquely written as

if $\widehat T_i\neq \emptyset$ and $(\mathbf{X}_{\widehat T_i})^\top \mathbf{X}_{\widehat T_i}$ has full rank $(|\widehat{T}_i|)$. To state Theorem (ref) below, let $\|\boldsymbol{\delta}\|_{\ell_0}:=|\operatorname{supp}(\boldsymbol{\delta})|$ be the $\ell_0$-“norm” of $\boldsymbol{\delta}$, denote

and define the $m$-sparse eigenvalue of a symmetric $pq\times pq$ matrix $\mathbf{A}$ as

Imposing that $\limsup_{n\to\infty}\phi_{\max}(s\ln(n), \boldsymbol{\Sigma}_{\boldsymbol{Z}})<\infty$ in addition to Assumptions (ref)--(ref), the post-Lasso following the Lasso implemented with the penalty parameter $\lambda_n^\ast$ and loadings $\widehat\boldsymbol{\Upsilon}_i^{(K)}$ from Algorithm (ref) obeys the following performance guarantees:

When $\lambda$ is chosen appropriately, Theorem (ref) provides the same asymptotic guarantees for the post-Lasso as Theorem (ref) does for the first stage Lasso. However, note that Theorem (ref) imposes the assumption $\limsup_{n\to\infty}\phi_{\max}(s\ln(n), \boldsymbol{\Sigma}_{\boldsymbol{Z}})<\infty$, which is not imposed for the Lasso in Theorem (ref). An assumption of this type was employed already for independent data in \citet*[Corollary 1]{belloni_sparse_2012} in order to establish rates of convergence in estimation and prediction error for the post-Lasso. This additional assumption is used to (i) upper bound the number of irrelevant variables retained by the Lasso and (ii) to control the omitted variable bias of the second stage least squares estimator (as the first stage Lasso may not retain all relevant variables). Thus, although the post-Lasso performs well in many of our simulations (see Section (ref)) with $\lambda_n^*$ as dictated by Algorithm (ref), it does impose an additional assumption compared to the Lasso. Moreover, the same simulations show that post-Lasso can perform worse than the latter if this additional assumption is not satisfied.

Performance Guarantees for the Sqrt-Lasso

In the context of i.i.d. data the sqrt-Lasso was proposed by belloni2011square to construct a shrinkage estimator for which a theoretically justifiable tuning parameter choice does not depend on the (unknown) scale of the error term(s).\footnote{The sqrt-Lasso coincides with the scaled Lasso of sun_scaled_2012 after concentrating the error standard deviation out of their objective function.} For this reason, it has become a popular alternative to the Lasso. We propose a data-driven implementation of the sqrt-Lasso for high-dimensional VAR models, which also accounts for the scales of the regressors $\bm{Z}_{t-1}$ through penalty loadings $\widehat{\boldsymbol{\Upsilon}}_{i}$ not present in the original sqrt-Lasso.\footnote{belloni2011square treat the regressors as fixed/conditioned on throughout (ibid., p. 793). This view allows one to normalize their scales to any desired values. In the context of time series such normalization is not innocuous, as regressor pre-processing can substantially alter the dependence structure. Hence, care needs to be taken in handling the regressor scales when implementing the sqrt-Lasso for time series.} In particular, we study the following sqrt-Lasso

and set the penalty level $\lambda=\lambda_n^*/2$ [cf. (ref)] and loadings $\widehat{\boldsymbol{\Upsilon}}_{i}=\dot{\boldsymbol{\Upsilon}}=\mathrm{diag}(\dot\upsilon_1,\hdots,\dot\upsilon_{pq})$ with $\dot\upsilon_j:=(n^{-1}\sum_{t=1}^nZ^2_{t-1,j})^{1/2}$, $j\in[pq]$.\footnote{In contrast to the loadings $\{\widehat{\boldsymbol{\Upsilon}}_i\}_{i\in[p]}$ from Algorithm (ref), the sqrt-Lasso loadings $\dot{\boldsymbol{\Upsilon}}$ do not depend on $i$.} Crucially, the implementation of the proposed version of the sqrt-Lasso does not depend on any population unknowns---it is fully data-driven.

The proposed sqrt-Lasso $\dot\boldsymbol{\beta}_i(\lambda_n^\ast/2,\dot\boldsymbol{\Upsilon})$ that adjusts to the scales of the regressors $\boldsymbol{Z}_{t-1}$ through the choice of loadings $\dot{\boldsymbol{\Upsilon}}$ obeys the same performance bounds as the (post-)Lasso. The condition $\overline q=0$ used in Theorem (ref) implies that the $\{\boldsymbol{\varepsilon}_t\}_{t\in\mathbb Z}$ are i.i.d. and, in particular, conditionally homoskedastic given the past outcomes. Conditional homoskedasticity is a crucial element in the analysis of the sqrt-Lasso for even i.i.d. observations, cf. belloni2011square. Just like for i.i.d. data, we shall see in the simulations in Section (ref) that the sqrt-Lasso can be inferior to the Lasso in the presence of conditional heteroskedasticity.

Simulations

We next explore the finite-sample behavior of the methods studied in a list of experiments.

Designs

Seven experimental designs (labelled A--G) are considered, all of which have zero intercepts. The number of Monte Carlo replications is $1{,}000$ and the effective sample sizes are $n\in\{100, 200, \dotsc, 1000\}$. We consider system sizes $p\in\{16,32,64,128\}$, the largest of which was chosen to mimic the 127 series used in our empirical illustration (Section (ref)).

Designs \hyperlink{Design A}{A}, \hyperlink{Design B}{B} and \hyperlink{Design C}{C} below are inspired by kock_oracle_2015. For these designs, the innovations $\{\boldsymbol{\varepsilon}_{t}\}_{t\in\mathbb Z}$ are i.i.d. mean zero Gaussian with a diagonal covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}}=\sigma_{\varepsilon}^2\mathbf{I}_{p}$ and $\sigma_{\varepsilon}=0.1$. The remaining details of these designs are:

As in Design \hyperlink{Design A}{A}, Designs \hyperlink{Design D}{D}, \hyperlink{Design E}{E} and \hyperlink{Design F}{F} all involve a VAR(1) with $\boldsymbol{\Theta}_{01}=(0.5)\mathbf{I}_p$ and $\sigma_\varepsilon=0.1$. They differ from Design \hyperlink{Design A}{A} in terms of $\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}}$, the heavy-tailedness of the $\boldsymbol{\varepsilon}_t$, or allowing for conditional heteroskedasticity. Design \hyperlink{Design G}{G} studies a local-to-unity variant of Design \hyperlink{Design A}{A}. The details of Designs \hyperlink{Design D}{D}--\hyperlink{Design G}{G} are as follows:

Implementation and Performance Measure

The Lasso, post-Lasso, and sqrt-Lasso are implemented as follows:

We specify $c$ and $\gamma_n$ as in (ref) and $K=15$.

For all procedures, the intercept is (correctly) enforced to be zero and the number of lags included is the smallest one ensuring that the model is correctly specified (i.e. four in Design \hyperlink{Design C}{C} and one for all other designs). The number of parameters to be estimated therefore ranges from $256$ in a VAR(1) of size $p=16$ to 65{,}536 in a VAR(4) of size $p=128.$

All simulations are carried out in R with user-written functions for each of the three above-mentioned estimators, taking as (basic) inputs the time series data (a matrix of dimensions $(n+q)\times p$) and desired lag length $(q)$.\footnote{Each Lasso problem (ref) is solved using glmnet. For the sqrt-Lasso, we downloaded the Matlab\textsuperscript{\textregistered} implementation of the coordinatewise method from Alexandre Belloni's website and translated it to R.}

Results

We study the maximum row-wise $\ell_{2}$-estimation error $\max_{i\in[p]}\Vert\widehat{\boldsymbol{\beta}}_i-\boldsymbol{\beta}_{0i}\Vert_{\ell_{2}}$ for which we report the average across the 1,000 Monte Carlo replications relative to that of the weighted Lasso (henceforth: Lasso). Thus, numbers less than one mean a procedure outperforms the Lasso. Figure (ref) plots the relative errors as a function of the sample size $n$, system size $p$ and design. The raw (non-relative) estimation errors can be found in Figure (ref) of Section (ref) in the appendix.

In Design \hyperlink{Design A}{A} each equation has only one relevant variable, which is uncorrelated with the irrelevant ones. Thus, the Lasso does well in terms of model selection and, as a result, the post-Lasso works very well here (for all $p$). Design \hyperlink{Design B}{B} is qualitatively similar, but note that the sqrt- and post-Lasso are less precise than the Lasso for $p=128$ and $n$ sufficiently small. In the block-diagonal Design \hyperlink{Design C}{C}, the Lasso and the sqrt-Lasso perform similarly, but for $p/n$ large the post-Lasso performs worse.

Design \hyperlink{Design D}{D} is challenging in the sense that the single relevant variable in each equation is highly correlated with the irrelevant ones. This makes model selection difficult and explains why the post-Lasso does relatively worse in this setting. The sqrt-Lasso also suffers for $p/n$ large. Design \hyperlink{Design E}{E} adds heavy tails to Design \hyperlink{Design D}{D}. This leaves the relative estimation errors unaffected.

Design \hyperlink{Design F}{F} confirms that the sqrt-Lasso may suffer substantially under heteroskedasticity (in particular for large $p$), cf. the discussion surrounding Theorem (ref). Finally, Design \hyperlink{Design G}{G} modifies Design \hyperlink{Design A}{A} to a near unit root design. This results in a reversal of the relative performance of the post- and sqrt-Lasso.

Empirical Illustration

We apply the methods from the above simulations to forecast a large set of macroeconomic variables using the Federal Reserve Economic Data (FRED) monthly data (MD) database. This database is maintained and regularly updated by the Federal Reserve Bank of St. Louis and described in detail on Michael W. McCracken's website.\footnote{\href{https://research.stlouisfed.org/econ/mccracken/fred-databases/}{https://research.stlouisfed.org/econ/mccracken/fred-databases/}. See also mccracken_FRED_2016.} The data is pre-processed in a standard manner using the Matlab\textsuperscript{\textregistered} code on McCracken's website leaving us with $758$ observations on $p=127$ macroeconomic variables covering March 1959 through April 2022.\footnote{The data pre-processing amounts to carrying out (deterministic) stationarity inducing transformations (prepare_missing.m), then removing outliers (remove_outliers.m), and finally replacing missing values with the unconditional average of the corresponding series (as in the initialization of factors_em.m).}

Forecasting

In each of the last $120$ months (i.e. ten years) of the sample we forecast the $p=127$ variables one month ahead. Specifically, we estimate VAR$(q)$ models of orders $q\in[12]$ using a rolling estimation window of size $n=758-120-12=626,$ and create one-month-ahead out-of-sample forecasts as $\widehat{\boldsymbol{Y}}_{t+1}:=\widehat{\boldsymbol{\mu}}+\sum_{j=1}^{q}\widehat{\boldsymbol{\Theta}}_{j}\boldsymbol{Y}_{t+1-j}$, where $\widehat{\boldsymbol{\mu}}$ and $\{\widehat{\boldsymbol{\Theta}}_{j}\}_{j=1}^q$ have not seen $\boldsymbol{Y}_{t+1}$. The three methods (weighted Lasso, post-Lasso and sqrt-Lasso) are implemented as in the simulations in Section (ref), but we now include an unpenalized intercept in each equation to account for non-zero means of the variables.\footnote{See Section (ref) for a modification of Algorithm (ref), which includes unpenalized intercepts.} For each method, each $q$, and each month $t\in\left\{ 639,640,\dotsc,758\right\} $ to be forecast, we calculate the forecast errors $\widehat{\boldsymbol{Y}}_{t}-\boldsymbol{Y}_{t}.$ Due to the different scaling of the variables, we then calculate the inverse-variance-weighted squared forecast error (IVWSFE) \[ \mathrm{IVWSFE}_{t}:=\sum_{i=1}^{p}(\widehat{Y}_{t,i}-Y_{t,i})^{2}/\widehat{\sigma}_{Y_{i}}^{2}, \] where $\widehat{\sigma}_{Y_{i}}^{2}$ denotes the sample variance of the $i$\textsuperscript{th} variable (over the entire pre-processed series).

Figure (ref) shows both the average and 95\textsuperscript{th} percentile of the 120 IVWSFEs of the weighted Lasso, post-Lasso and sqrt-Lasso, respectively. To facilitate comparison, these measures are put relative to that of the weighted Lasso with $q=1$. In terms of the average IVWSFE error, the sqrt-Lasso does slightly better than both the weighted Lasso and post-Lasso, no matter the choice of lag length $q$, hovering between 96 and 97 pct. of VAR(1) Lasso. However, when looking at the 95\textsuperscript{th} percentile IVWSFE, the picture is somewhat reversed, in that the post-Lasso here outperforms the sqrt-Lasso (and weighted Lasso). Thus, these methods cannot be ranked in terms of the quality of their forecasts, in general.\footnote{We also experimented with equation-by-equation ordinary least squares (OLS). With a lag order of one, OLS led to an average IVWSFE 40 percent higher than our weighted Lasso benchmark (for $q=1$). For higher lag orders, the implied design matrix repeatedly fell short of full rank, thus preventing a meaningful comparison.}