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

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

64,005 characters

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




\title{Data-Driven Tuning Parameter Selection for High-Dimensional Vector Autoregressions
\thanks{First version: October 2022. We thank Robert Adamek, Peter Boswijk, Abhishek Chakrabortty, Lan Gao, Frank Kleibergen, Bent Nielsen, Anders Rahbek and Mikkel S{\o}lvsten for insightful comments and discussion. We also thank participants at the Aarhus
Workshops in Econometrics I and IV, the  2024 Bristol econometrics study group conference, the 35\textsuperscript{th} (EC)\textsuperscript{2} conference, and seminar participants at University of Amsterdam for suggestions and discussion.
Moreover, we thank Zifan Li for providing details about the simulation
experiment in \citet{wong_lasso_2020}. Pedersen gratefully acknowledges
financial support from Independent Research Fund Denmark (FSE 0133-00162B).}
}

\author{Anders Bredahl Kock\footnote{University of Oxford  and Aarhus Center for Econometrics (ACE); e-mail:
\tt{[email removed]}.} \and Rasmus S{\o}ndergaard Pedersen\footnote{University of Copenhagen and Danish Finance Institute; e-mail:
\tt{[email removed]}.} \and Jesper
Riis-Vestergaard S{\o}rensen\footnote{University of Copenhagen and Aarhus Center for Econometrics (ACE); e-mail:
\tt{[email removed]}.}}

\maketitle

\begin{abstract}
Lasso-type estimators are routinely used to estimate high-dimensional time series models. The theoretical guarantees established for these estimators typically require the penalty level to be chosen in a suitable fashion often depending on \emph{unknown} population quantities. Furthermore, the resulting estimates and the number of variables retained in the model depend crucially on the chosen penalty level. However, there is currently no theoretically founded guidance for this choice in the context of high-dimensional time series. Instead, one resorts to selecting the penalty level in an ad hoc manner using, e.g., information criteria or cross-validation. We resolve this problem by considering estimation of the perhaps most commonly employed multivariate time series model, the linear vector autoregressive (VAR) model, and propose versions of the Lasso, post-Lasso, and square-root Lasso estimators with penalization chosen in a \emph{fully data-driven} way. The theoretical guarantees that we establish for the resulting estimation and prediction errors match those currently available for methods based on infeasible choices of penalization. We thus provide a first solution for choosing the penalization in high-dimensional time series models.

\end{abstract}

\medskip\noindent\textbf{Keywords:} High-dimensional time series, vector autoregressive model, $\ell_1$-penalized estimation, data-driven tuning parameter selection.



\section{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., \citet{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 [\citet{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 \citet{han_transition_2013}, \citet{basu_regularized_2015}, \citet{kock_oracle_2015} and \citet{davis_sparse_2016}, among
others. The assumption of Gaussian innovations has been
relaxed in several papers, including \citet{song_large_2011}, \citet{wong_lasso_2020} and \citet{miao_var_2022}. In
particular, \citet{wong_lasso_2020} derive consistency results for
the Lasso in VAR models under suitable mixing and moment conditions, whereas \citet{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
\citet{wu_nonlinear_2005}]. \citet{masini_VAR_2022} establish oracle inequalities
in the case of martingale difference innovations that are not necessarily
mixing. \citet{chernozhukov_lasso_2021} consider Lasso-driven inference for time series and spatial models, while \citet{adamek_lasso_2022} consider inference based on the desparsified Lasso when the data-generating process allows near-epoch dependence. Further related papers are \citet{han_direct_2015}, \citet{guo_var_2016} and \citet{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 \emph{unknown} population quantities such as mixing coefficients or other coefficients quantifying
the dependence structure of the data generating process {[}\citet{wong_lasso_2020},
\citet{babii_timeseries_2022}, \citet{masini_VAR_2022}{]}, the population
covariance matrix of the observed variables and innovations {[}\citet{kock_oracle_2015},
\citet{MEDEIROS_highdim_2016}{]}, the population coefficient matrix
{[}\citet{basu_regularized_2015}{]}, or other quantities depending
on the data generating process [\citet{chernozhukov_lasso_2021}, \citet{adamek_lasso_2022}, \citet{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.~\citet{kock_oracle_2015}, \citet{MEDEIROS_highdim_2016} and
\citet{masini_VAR_2022}], cross-validation [e.g. \citet{wong_lasso_2020}, \citet{babii_timeseries_2022} and \citet{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 \citet{javanmard_confidence_2014}, \citet{van_de_geer_asymptotically_2014} and \citet{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 \citet{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 \citet{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 \citet{wu_nonlinear_2005} functional dependence sense]
derived by \citet{chen_self-normalized_2016} and \citet{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~\cite{belloni2011square} is a popular alternative to the Lasso for i.i.d.~data in the absence of conditional heteroskedasticity. Apart from the work of~\cite{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.



\subsection*{Outline}

In Section \ref{sec:Model-and-Estimation} we present the weighted Lasso estimator and describe the data-driven tuning parameter selection. Section \ref{sec:Ass+Res} presents the assumptions and performance guarantees for the weighted Lasso. Sections~\ref{sec:postlasso} and \ref{sec:SqrtLasso} cover the post- and sqrt-Lasso, respectively. Sections \ref{sec:Simulations} and \ref{sec:Empirical-Illustration} contain simulations and an empirical illustration, respectively. Proofs, technical lemmas, and additional simulation output are contained in the supplementary appendices.

\subsection*{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$.

\section{Model and Penalization Algorithm\label{sec:Model-and-Estimation}}
We study the $q$\textsuperscript{th}-order $(q\in\mathbb N)$ VAR model given by
\begin{equation}
\boldsymbol{Y}_{t}=\sum_{j=1}^{q}\boldsymbol{\Theta}_{0j}\boldsymbol{Y}_{t-j}+\boldsymbol{\varepsilon}_{t},\quad t\in\mathbb{Z},\label{eq:VARvector}
\end{equation}
yielding a stochastic process $\{\boldsymbol{Y}_t\}_{t\in\mathbb Z}$ which is (strictly) stationary under the assumptions in Section~\ref{sec:Ass+Res}. 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{eq:VARvector}),
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{eq:VARvector}), 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{eq:VARvector}) may be written in companion form,
\begin{equation}
\boldsymbol{Z}_{t}=\widetilde{\boldsymbol{\Theta}}_{0}\boldsymbol{Z}_{t-1}+\widetilde{\boldsymbol{\varepsilon}}_{t},\quad t\in\mathbb Z,\label{eq:VARcompanion}
\end{equation}
for
\[
\underbrace{\boldsymbol{Z}_{t}}_{pq\times1}:=\begin{pmatrix}
\boldsymbol{Y}_{t}\\
\boldsymbol{Y}_{t-1}\\
\vdots\\
\boldsymbol{Y}_{t-(q-1)}
\end{pmatrix},
\quad
\underbrace{\widetilde{\boldsymbol{\Theta}}_{0}}_{pq\times pq}:=\begin{pmatrix}\boldsymbol{\Theta}_{01} & \boldsymbol{\Theta}_{02} & \cdots &  & \boldsymbol{\Theta}_{0q}\\
\mathbf{I}_{p} & \boldsymbol{0}_{p\times p} & \cdots &  & \boldsymbol{0}_{p\times p}\\
\boldsymbol{0}_{p\times p} & \mathbf{I}_{p} & \boldsymbol{0}_{p\times p} &  & \vdots\\
\vdots & \ddots & \ddots & \ddots \\
\boldsymbol{0}_{p\times p} & \cdots & \boldsymbol{0}_{p\times p} & \mathbf{I}_{p} & \boldsymbol{0}_{p\times p}
\end{pmatrix},
\quad\underbrace{\widetilde{\boldsymbol{\varepsilon}}_{t}}_{pq\times1}:=\begin{pmatrix}\boldsymbol{\varepsilon}_{t}\\
\boldsymbol{0}_{p\times 1}\\
\vdots\\
\boldsymbol{0}_{p\times 1}
\end{pmatrix},
\]
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{eq:VARvector})
is given by
\begin{equation}
Y_{t,i}=\boldsymbol{Z}_{t-1}^{\top}\boldsymbol{\beta}_{0i}+\varepsilon_{t,i},\label{eq:VARequation}
\end{equation}
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
\begin{equation}
\widehat{\boldsymbol{\beta}}_{i}:=\widehat{\boldsymbol{\beta}}_{i}(\lambda,\widehat{\boldsymbol{\Upsilon}}_{i})\in\operatornamewithlimits{argmin}\limits_{\boldsymbol{\beta}\in\mathbb R^{pq}}\left\{ \widehat Q_i(\boldsymbol{\beta}) +\frac{\lambda}{n}\Vert\widehat{\boldsymbol{\Upsilon}}_{i}\boldsymbol{\beta}\Vert_{\ell_{1}}\right\},\label{eq:LASSOVector}
\end{equation}
where~$\widehat Q_i$ is the sample average squared error loss function
\begin{equation}\label{eq:Qihat}
\widehat Q_i(\boldsymbol{\beta}):=\frac{1}{n}\sum_{t=1}^{n}\left(Y_{t,i}-\boldsymbol{Z}_{t-1}^{\top}\boldsymbol{\beta}\right)^{2},\quad\boldsymbol{\beta}\in\mathbb R^{pq},
\end{equation}
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{alg:Data-Driven-Penalization} below. Although the minimization problem in (\ref{eq:LASSOVector}) may have multiple solutions, we sometimes refer to such a $\widehat{\boldsymbol{\beta}}_{i}$ as \textit{the} Lasso estimator. The results established below apply to any such (measurable) minimizer.\footnote{By the measurable selection theorem 6.7.22~in \citet{pfanzagl},~$\widehat{\boldsymbol{\beta}}_{i}$ in~\eqref{eq:LASSOVector} 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 \citet[Algorithm A.1]{belloni_sparse_2012} for independent data to our setting,
we consider the following data-driven penalization. Let
\begin{equation}
c:=1.1\quad\text{and}\quad\gamma_{n}=0.1/\ln(\max\{n,pq\})\label{eq:tuning_c_and_gamma}
\end{equation}
and let $K\in\mathbb{N}_{0}$ denote a fixed number of loading updates.\footnote{Although Theorem \ref{thm:Rates-for-LASSO-data-driven-loadings} is valid for all choices of~$c>1$, our simulations (Section \ref{sec:Simulations}) 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{alg:Data-Driven-Penalization} is not sensitive to the exact choice. This is documented in Figure~\ref{fig:cchoice} in Section \ref{sec:appsim} 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~\cite{belloni_sparse_2012}.} Following \citet{belloni_sparse_2012}, we fix $K$ at 15 in our simulations and
empirical illustration. To state Algorithm \ref{alg:Data-Driven-Penalization} below, let $\Phi$ be the standard Gaussian cumulative distribution function.
\begin{mdframed}
\begin{lyxalgorithm}
[\textbf{Data-Driven Penalization}]\label{alg:Data-Driven-Penalization}\ \\
\textbf{Initialize}: Specify the penalty level~$\lambda$ in~\eqref{eq:LASSOVector} as
\begin{equation}
\lambda^{\ast}_n := 2c\sqrt{n}\Phi^{-1}\big(1-\gamma_n/(2p^{2}q)\big),\label{eq:PenaltyLevelPractice}
\end{equation}
and specify the initial penalty loadings as
\begin{equation}
\widehat{\upsilon}_{i,j}^{\left(0\right)} := \sqrt{\frac{1}{n}\sum_{t=1}^{n}Y_{t,i}^{2}Z_{t-1,j}^{2}},\quad i\in[p],\quad j\in[pq].\label{eq:PenaltyLoadingsInitial}
\end{equation}
Use $\lambda^{\ast}_n$ and $\widehat{\boldsymbol{\Upsilon}}_{i}^{\left(0\right)}:=\mathrm{diag}(\widehat{\upsilon}_{i,1}^{\left(0\right)},\dotsc,\widehat{\upsilon}_{i,pq}^{\left(0\right)})$
to compute a Lasso estimate $\widehat{\boldsymbol{\beta}}_{i}^{\left(0\right)}:=\widehat{\boldsymbol{\beta}}_{i}(\lambda^{\ast}_n,\widehat{\boldsymbol{\Upsilon}}_{i}^{\left(0\right)})$
via (\ref{eq:LASSOVector}) for each $i\in[p]$. Store the residuals $\widehat{\varepsilon}_{t,i}^{\left(0\right)}:=Y_{t,i}-\boldsymbol{Z}_{t-1}^{\top}\widehat{\boldsymbol{\beta}}_{i}^{\left(0\right)},t\in[n],i\in[p]$, and set $k=1$.

\noindent\textbf{Update:}
While $k\leqslant K$, specify the penalty loadings as
\begin{equation}
\widehat{\upsilon}_{i,j}^{\left(k\right)}:=\sqrt{\frac{1}{n}\sum_{t=1}^{n}(\widehat{\varepsilon}_{t,i}^{\left(k-1\right)})^2Z_{t-1,j}^{2}},\quad i\in[p],\quad j\in[pq].\label{eq:PenaltyLoadingsRefined}
\end{equation}
Use $\lambda_n^{\ast}$ and $\widehat{\boldsymbol{\Upsilon}}_{i}^{\left(k\right)}:=\mathrm{diag}(\widehat{\upsilon}_{i,1}^{\left(k\right)},\dotsc,\widehat{\upsilon}_{i,pq}^{\left(k\right)})$
to compute a Lasso estimate $\widehat{\boldsymbol{\beta}}_{i}^{\left(k\right)}:=\widehat{\boldsymbol{\beta}}_{i}(\lambda_n^{\ast},\widehat{\boldsymbol{\Upsilon}}_{i}^{\left(k\right)})$
via (\ref{eq:LASSOVector}) for each $i\in[p]$. Store the residuals $\widehat{\varepsilon}_{t,i}^{\left(k\right)}:=Y_{t,i}-\boldsymbol{Z}_{t-1}^{\top}\widehat{\boldsymbol{\beta}}_{i}^{\left(k\right)},t\in[n],i\in[p]$,
and increment $k\leftarrow k+1.$
\end{lyxalgorithm}
\end{mdframed}
\noindent Note that Algorithm \ref{alg:Data-Driven-Penalization} 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
\begin{equation}
s:=\max_{i\in[p]}\sum_{j=1}^{pq}\textbf{1}(\beta_{0i,j}\neq0),\label{eq:def_s_sparsity}
\end{equation}
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{sec:Ass+Res}, any
Lasso estimators $\widehat\boldsymbol{\beta}_i:=\widehat\boldsymbol{\beta}_{i}(\lambda_n^{\ast},\widehat\boldsymbol{\Upsilon}_i^{(K)}),i\in[p]$, in (\ref{eq:LASSOVector}) based on the penalty level $\lambda_n^{\ast}$ in (\ref{eq:PenaltyLevelPractice})
and the (final) penalty loadings $\widehat\boldsymbol{\Upsilon}_i^{(K)}$ arising from Algorithm \ref{alg:Data-Driven-Penalization} satisfy
\begin{align}
\max_{i\in[p]}\Vert\widehat{\boldsymbol{\beta}}_{i}-\boldsymbol{\beta}_{0i}\Vert_{2,n} & \lesssim_{\mathrm{P}}\sqrt{\frac{s\ln\left(pn\right)}{n}},\label{eq:lasso_rates_1}\\
\max_{i\in[p]}\Vert\widehat{\boldsymbol{\beta}}_{i}-\boldsymbol{\beta}_{0i}\Vert_{\ell_{1}} & \lesssim_{\mathrm{P}}\sqrt{\frac{s^{2}\ln\left(pn\right)}{n}}\quad\text{and}\label{eq:lasso_rates_2}\\
\max_{i\in[p]}\Vert\widehat{\boldsymbol{\beta}}_{i}-\boldsymbol{\beta}_{0i}\Vert_{\ell_{2}} & \lesssim_{\mathrm{P}}\sqrt{\frac{s\ln\left(pn\right)}{n}}.\label{eq:lasso_rates_3}
\end{align}
The formal statement, including a set of sufficient conditions, is given
in Theorem \ref{thm:Rates-for-LASSO-data-driven-loadings} below. Importantly, the performance guarantees in~\eqref{eq:lasso_rates_1}--\eqref{eq:lasso_rates_3} match those currently available based on \emph{infeasible} penalty level choices.

\begin{rem}[Including Intercepts]
In applications one typically includes model intercepts in an unpenalized manner. In Section \ref{subsec:IncludingIntercepts}, we modify Algorithm \ref{alg:Data-Driven-Penalization} for this purpose.\hfill$\diamondsuit$
\end{rem}

\section{Assumptions and Performance Guarantees for Algorithm \ref{alg:Data-Driven-Penalization}}\label{sec:Ass+Res}

For the remainder of the paper we invoke (a subset of) the following assumptions.
\begin{assumption}[\textbf{Innovations}] \label{assu:Innovations} There is a $\overline p\in\mathbb N$, a constant $\overline q\in\mathbb N_0$, a measurable mapping $\boldsymbol{F}:\mathbb R^{\overline p\cdot\overline q}\to\mathbb R^p$, an i.i.d.~process $\{\boldsymbol{\eta}_t\}_{t\in\mathbb Z}$, each $\boldsymbol{\eta}_t$ taking values in $\mathbb R^{\overline p}$, and constants $a_1,a_2,\alpha\in(0,\infty)$ such that: \begin{inparaenum}[(1)]\item\label{enu:EpsIID} $\boldsymbol{\varepsilon}_t=\boldsymbol{F}(\boldsymbol{\eta}_{t-\overline q},\dotsc,\boldsymbol{\eta}_t)$ for each $t\in\mathbb Z$.
\item\label{enu:EpsSubWeibullBeta} $\Vert \boldsymbol{\varepsilon}_{0}\Vert _{\psi_{\alpha}}\leqslant a_{1}$. \item\label{enu:EigValsOfSigmaEpsBddAwayFromZero} With $\mathcal{F}^{\boldsymbol{\eta}}_t$ denoting the $\sigma$-field generated by $\{\boldsymbol{\eta}_u\}_{u=-\infty}^{t}$, for each $t\in\mathbb Z$, $\mathrm{E}[\boldsymbol{\varepsilon}_t \mid \mathcal{F}^{\boldsymbol{\eta}}_{t-1}]=\mathbf{0}_p$ almost surely and $\min_{i\in[p]}\mathrm{E}[\varepsilon_{t,i}^2 \mid \mathcal{F}^{\boldsymbol{\eta}}_{t-1}]\geqslant a_{2}$ almost surely.
\end{inparaenum}
\end{assumption}
Assumption \ref{assu:Innovations}.\ref{enu:EpsIID} 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., \citet{kock_oracle_2015}, \citet[Example 1]{wong_lasso_2020} and \citet{miao_var_2022}. Assumption \ref{assu:Innovations}.\ref{enu:EpsSubWeibullBeta} states that $\boldsymbol{\varepsilon}_{0}$ belongs to the sub-Weibull family of distributions with tail parameter (at least) $\alpha$, as imposed in \citet{wong_lasso_2020} and \citet{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., \citet{han_transition_2013}, \citet{basu_regularized_2015}, \citet{kock_oracle_2015} and \citet{davis_sparse_2016}---then Assumption \ref{assu:Innovations}.\ref{enu:EpsSubWeibullBeta}
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{assu:Innovations}.\ref{enu:EigValsOfSigmaEpsBddAwayFromZero} parallels part of \citet[Assumption A2]{masini_VAR_2022} and  \citet[Assumption A.1(i)]{miao_var_2022}.


\begin{assumption}
[\textbf{Companion Matrix}]\label{assu:Companion}\textbf{ }\begin{inparaenum}[(1)]\item\label{enu:SpectralRadius}~The spectral
radius $\rho(\widetilde{\boldsymbol{\Theta}}_{0})<1$.
\item\label{enu:RowNormDecay} There are constants $b_1, b_2\in(0,\infty)$ and $\tau\in(0,1]$ such that
$\max_{i\in[p]}\Vert(\widetilde{\boldsymbol{\Theta}}_{0}^{h})_{i,1:p}\Vert_{\ell_{2}}\leqslant b_{1}\mathrm{e}^{-b_{2}h^\tau}$
for all $h\in\mathbb{N}$.
\item\label{enu:lag-order} The lag-order $q$ is constant.
\end{inparaenum}
\end{assumption}

Assumption \ref{assu:Companion}.\ref{enu:SpectralRadius} is standard and combined with Assumption \ref{assu:Innovations} ensures that the process $\{\boldsymbol{Y}_{t}\}_{t\in\mathbb Z}$ in (\ref{eq:VARvector}) [as well as
$\{\boldsymbol{Z}_{t}\}_{t\in\mathbb Z}$ in (\ref{eq:VARcompanion})] 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{assu:Companion}.\ref{enu:RowNormDecay} is milder than \citet[Assumption A1]{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 \citet[Assumption A.1(vi)]{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{assu:Companion}.\ref{enu:RowNormDecay}.}

Assumption \ref{assu:Companion}.\ref{enu:lag-order} 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 \citet{kock_oracle_2015}, although their analysis allows for growing $q$.} Indeed, in our application in Section \ref{sec:Empirical-Illustration} 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 \citet{gao_refined_2022} for self-normalized sums based on geometric moment contracting processes. Assumption \ref{assu:Companion} therefore seems hard to relax with currently available methods.

Assumptions \ref{assu:Innovations} and \ref{assu:Companion} 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{lem:ZandEpsZNormBnds} in the appendix for details.} We further assume:
\begin{assumption}
[\textbf{Covariance}]\label{assu:CovarianceZ}
There is a constant $d\in(0,\infty)$ such that
$\Lambda_{\min}(\boldsymbol{\Sigma}_{\boldsymbol{Z}})\geqslant d$.
\end{assumption}
Assumption \ref{assu:CovarianceZ} ensures that $\boldsymbol{\Sigma}_{\boldsymbol{Z}}$ has full rank and is implied by, e.g., \citet[Assumption A.1]{miao_var_2022}. As demonstrated in \citet[Section 3]{basu_regularized_2015}, Assumption \ref{assu:CovarianceZ} follows from additional (boundedness) conditions on~$\widetilde{\boldsymbol{\Theta}}_{0}$, although these are not necessary.\footnote{For the purpose of Theorem \ref{thm:Rates-for-LASSO-data-driven-loadings}, we note that Assumption \ref{assu:CovarianceZ} is used in the proof of Lemma \ref{lem:Restricted-Eigenvalue-Bound}. Inspection of its proof reveals that one can replace Assumption \ref{assu:CovarianceZ} by a milder requirement that only certain \textit{restricted} eigenvalues of $\boldsymbol{\Sigma}_{\boldsymbol{Z}}$ are bounded away from zero. Restricted eigenvalue conditions appear in, e.g., \citet{kock_oracle_2015} and \citet[DESIGN(3)a]{MEDEIROS_highdim_2016}.}
\begin{assumption}
[\textbf{Growth and Sparsity}]\label{assu:RowSparsity} For  constants $\alpha\in(0,\infty)$  and $\tau\in(\textstyle{\frac{1}{2}},1]$ satisfying respectively Assumptions \ref{assu:Innovations}.\ref{enu:EpsSubWeibullBeta} and \ref{assu:Companion}.\ref{enu:RowNormDecay}, we have $(\ln p)^{C(\alpha,\tau)}=o(n)$ and $s^2(\ln(pn))^{1/\tau+\widetilde{C}(2,\alpha)}=o(n)$, where $s$ is given in (\ref{eq:def_s_sparsity}), $\widetilde{C}(x,\alpha):=\max\{(2x/\alpha),(x+\alpha)/\alpha\}$ and
\[
C(\alpha,\tau):=\max\left\{\frac{\widetilde{C}(4,\alpha)(1+4\tau)}{4\tau-1},\frac{1+4\tau}{4\tau-2},\frac{1+4\tau}{\tau}\right\}.
\]
\end{assumption}
\begin{comment}
{\color{red}
\begin{itemize}
    \item J: Perhaps add reference to \citet[Assumption 8]{wong_lasso_2020}, which similarly restricts the dependence (there measured in terms of geometric decay in beta-mixing coefficients).
\end{itemize}
}

\end{comment}
Assumption \ref{assu:RowSparsity}
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., \citet{kock_oracle_2015},
\citet{chernozhukov_lasso_2021}, and \citet{masini_VAR_2022}.~To help digest the growth and sparsity conditions in Assumption \ref{assu:RowSparsity}, 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{alg:Data-Driven-Penalization}.

\begin{thm}[\textbf{Convergence Rates for Lasso with Data-Driven Penalization}]\label{thm:Rates-for-LASSO-data-driven-loadings}
Let Assumptions \ref{assu:Innovations}, \ref{assu:Companion}, \ref{assu:CovarianceZ} and \ref{assu:RowSparsity} hold and fix~$K\in\mathbb N_0$. Then Lasso estimators $\widehat\boldsymbol{\beta}_i(\lambda_n^\ast,\widehat\boldsymbol{\Upsilon}_i^{(K)}),i\in[p]$, arising from the data-driven penalization in Algorithm \ref{alg:Data-Driven-Penalization} satisfy the rates in~(\ref{eq:lasso_rates_1})--(\ref{eq:lasso_rates_3}).
\end{thm}

The estimation and prediction error rates guaranteed by Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings} match those currently available in the literature building on infeasible penalty parameter choices. Algorithm~\ref{alg:Data-Driven-Penalization} 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~\citet{javanmard_confidence_2014}, \citet{van_de_geer_asymptotically_2014} and \citet{zhang_confidence_2014}. This methodology has been used in the context of time series in \citet{adamek_lasso_2022}. Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings} therefore opens the door to inference in high-dimensional VAR models with data-driven penalty parameter choice.

In establishing Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings}, 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{lem:MaximalInequalityFunctionalDependence} in the appendix). This result is then used to show that the penalty loadings of Algorithm~\ref{alg:Data-Driven-Penalization} are close to certain ideal, yet infeasible, blocking-based loadings (Lemma~\ref{lem:AsymptoticValidityDataDrivenPenaltyLoadings}). Note that although these ideal loadings involve blocking, Algorithm~\ref{alg:Data-Driven-Penalization} 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 \cite{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{sec:proofsteps} in the appendix explains further steps and challenges involved in establishing Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings}.


\begin{rem}[Penalty Loading Updates]
The error rates (\ref{eq:lasso_rates_1})--(\ref{eq:lasso_rates_3}) guaranteed by Theorem \ref{thm:Rates-for-LASSO-data-driven-loadings} are valid for any fixed $K\in\mathbb N_0$ in Algorithm \ref{alg:Data-Driven-Penalization}. Preliminary and unreported simulations indicate that there are typically finite-sample benefits from updating the loadings (as opposed to no updating, $K=0$). The same exercises indicate that, it suffices to use~$K=5$ updates, as neither the loadings nor the parameter estimates change much for larger~$K$.\hfill$\diamondsuit$
\end{rem}

\section{Performance Guarantees for the Post-Lasso}\label{sec:postlasso}
To alleviate shrinkage bias, the \textit{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:
\begin{enumerate}[(i)]
    \item The first step Lasso estimator need not be unique and the elements of~$\mathrm{argmin}_{\boldsymbol{\beta}\in\mathbb R^{pq}}\{\widehat Q_i(\boldsymbol{\beta}) +(\lambda/n)\Vert\widehat{\boldsymbol{\Upsilon}}_{i}\boldsymbol{\beta}\Vert_{\ell_{1}}\}$
    need not select the same variables, cf.~\cite{tibshirani2013lasso}. As a results, the post-Lasso depends on which Lasso estimator one picks.
    \item Even for a given~$\widehat T_i$ it need not be the case that $\mathrm{argmin}_{\boldsymbol{\beta}\in\mathbb R^{pq}:\operatorname{supp}(\boldsymbol{\beta})\subseteq \widehat{T}_i}\widehat{Q}_{i}(\boldsymbol{\beta})$
    is unique as~$(\mathbf{X}_{\widehat T_i})^\top \mathbf{X}_{\widehat T_i}$ may not have full rank.
\end{enumerate}
Nevertheless, just like Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings}, Theorem~\ref{thm:postLasso} below applies to any post-Lasso estimator and guarantees uniqueness with probability tending to one.}
\begin{equation}\label{eq:PostLassom}
\widetilde{\boldsymbol{\beta}}_{i}\in\operatornamewithlimits{argmin}\limits_{\mathclap{\substack{\boldsymbol{\beta}\in\mathbb R^{pq}:\\ \operatorname{supp}(\boldsymbol{\beta})\subseteq \widehat{T}_i}}} \widehat{Q}_{i}\del[0]{\boldsymbol{\beta}},\qquad i\in[p].
\end{equation}
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{alg:Data-Driven-Penalization} we have provided such guarantees for the first stage Lasso~$\widehat{\boldsymbol{\beta}}_i(\lambda_n^*,\widehat{\boldsymbol{\Upsilon}}_i^{(K)})$ in Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings} 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
\begin{align}\label{eq:PLdef}
	\widetilde{\boldsymbol{\beta}}_{i,\widehat T_i}
	=
	\del[1]{(\mathbf{X}_{\widehat T_i})^\top \mathbf{X}_{\widehat T_i}}^{-1}(\mathbf{X}_{\widehat T_i})^\top \bm{y}_{i},
\end{align}
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{thm:postLasso} below, let $\|\boldsymbol{\delta}\|_{\ell_0}:=|\operatorname{supp}(\boldsymbol{\delta})|$ be the $\ell_0$-``norm'' of $\boldsymbol{\delta}$, denote
\begin{align}\label{eq:Dm}
	\mathfrak{D}(m):=\cbr[1]{\boldsymbol{\delta}\in\mathbb R^{pq}:\Vert\boldsymbol{\delta}\Vert_{\ell_0}\leqslant m, \Vert\boldsymbol{\delta}\Vert_{\ell_2}= 1}\qquad\text{for } m\in[1,\infty),
\end{align}
and define the~\emph{$m$-sparse eigenvalue} of a symmetric~$pq\times pq$ matrix~$\mathbf{A}$ as
\begin{equation}
	\phi_{\max}(m,\mathbf{A}):=\max_{\boldsymbol{\delta}\in\mathfrak{D}(m)}\boldsymbol{\delta}^\top\mathbf{A}\boldsymbol{\delta}.\label{eq:largestmsparseeigenvalue}
\end{equation}
Imposing that $\limsup_{n\to\infty}\phi_{\max}(s\ln(n), \boldsymbol{\Sigma}_{\boldsymbol{Z}})<\infty$ in addition to Assumptions~\ref{assu:Innovations}--\ref{assu:RowSparsity}, the post-Lasso following the Lasso implemented with the penalty parameter~$\lambda_n^\ast$ and loadings~$\widehat\boldsymbol{\Upsilon}_i^{(K)}$ from Algorithm~\ref{alg:Data-Driven-Penalization} obeys the following performance guarantees:
\begin{thm}[\textbf{Convergence Rates for Post-Lasso with Data-Driven Penalization}]\label{thm:postLasso}
	Let Assumptions \ref{assu:Innovations}, \ref{assu:Companion}, \ref{assu:CovarianceZ} and \ref{assu:RowSparsity} hold and fix~$K\in\mathbb N_0$. Assume that~$\limsup_{n\to\infty}\phi_{\max}(s\ln(n), \boldsymbol{\Sigma}_{\boldsymbol{Z}})<\infty$. Then post-Lasso estimators~$\widetilde{\boldsymbol{\beta}}_{i}(\lambda_n^\ast,\widehat\boldsymbol{\Upsilon}_i^{(K)}),i\in[p]$, from~\eqref{eq:PostLassom} following Lasso estimators $\widehat\boldsymbol{\beta}_i(\lambda_n^\ast,\widehat\boldsymbol{\Upsilon}_i^{(K)}),i\in[p]$, arising from the data-driven penalization in Algorithm \ref{alg:Data-Driven-Penalization}, satisfy the rates in~\eqref{eq:lasso_rates_1}--\eqref{eq:lasso_rates_3}.
	\end{thm}
	\emph{When~$\lambda$ is chosen appropriately}, Theorem~\ref{thm:postLasso} provides the same asymptotic guarantees for the post-Lasso as Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings} does for the first stage Lasso. However, note that Theorem~\ref{thm:postLasso} imposes the assumption~$\limsup_{n\to\infty}\phi_{\max}(s\ln(n), \boldsymbol{\Sigma}_{\boldsymbol{Z}})<\infty$, which is \emph{not} imposed for the Lasso in Theorem~\ref{thm:Rates-for-LASSO-data-driven-loadings}. 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{sec:Simulations}) with~$\lambda_n^*$ as dictated by Algorithm \ref{alg:Data-Driven-Penalization}, 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.

\section{Performance Guarantees for the Sqrt-Lasso}\label{sec:SqrtLasso}
In the context of i.i.d.~data the sqrt-Lasso  was proposed by~\cite{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 \textit{scaled Lasso} of \citet{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{\citet{belloni2011square} treat the regressors as fixed/conditioned on throughout (\textit{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
\begin{align}\label{eq:srtLasso}
	\dot{\boldsymbol{\beta}}_{i}:=\dot{\boldsymbol{\beta}}_{i}(\lambda,\widehat{\boldsymbol{\Upsilon}}_{i})\in\operatornamewithlimits{argmin}\limits_{\boldsymbol{\beta}\in\mathbb R^{pq}}\left\{ \del[1]{\widehat Q_i(\boldsymbol{\beta})}^{1/2} +\frac{\lambda}{n}\Vert\widehat{\boldsymbol{\Upsilon}}_{i}\boldsymbol{\beta}\Vert_{\ell_{1}}\right\}
\end{align}
and set the penalty level~$\lambda=\lambda_n^*/2$ [cf.~\eqref{eq:PenaltyLevelPractice}] 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{alg:Data-Driven-Penalization}, the sqrt-Lasso loadings~$\dot{\boldsymbol{\Upsilon}}$ do not depend on~$i$.}
Crucially, the implementation of the proposed version of the sqrt-Lasso does \emph{not} depend on any population unknowns---it is fully data-driven.
\begin{thm}[\textbf{Convergence Rates for Sqrt-Lasso with Data-Driven Penalization}]\label{thm:Rates-for-sqrtLASSO-data-driven-loadings}
Let Assumptions \ref{assu:Innovations}, \ref{assu:Companion}, \ref{assu:CovarianceZ} and \ref{assu:RowSparsity} hold with $\overline q=0$. Then sqrt-Lasso estimators $\dot\boldsymbol{\beta}_i(\lambda_n^\ast/2,\dot\boldsymbol{\Upsilon}),i\in[p]$, arising from~\eqref{eq:srtLasso} satisfy the rates in~(\ref{eq:lasso_rates_1})--(\ref{eq:lasso_rates_3}).
\end{thm}
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{thm:Rates-for-sqrtLASSO-data-driven-loadings} 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.~\citet{belloni2011square}. Just like for i.i.d.~data, we shall see in the simulations in Section \ref{sec:Simulations} that the sqrt-Lasso can be inferior to the Lasso in the presence of conditional heteroskedasticity.

\section{Simulations\label{sec:Simulations}}
We next explore the finite-sample behavior of the methods studied in a list of experiments.

\subsection{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{sec:Empirical-Illustration}).

Designs \hyperlink{Design A}{A}, \hyperlink{Design B}{B} and \hyperlink{Design C}{C} below are inspired by \citet[Section 5]{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:
\begin{itemize}
\item\textbf{Design \hypertarget{Design A}{A}:} The data-generating process (DGP) is a VAR(1) with diagonal matrix $\boldsymbol{\Theta}_{01}=(0.5)\mathbf{I}_{p}$
implying a spectral radius $\rho(\boldsymbol{\Theta}_{01})=0.5$. This setting
is very row sparse as $s=1.$
\item\textbf{Design \hypertarget{Design B}{B}:} The DGP is a VAR(1) with coefficient matrix having the Toeplitz structure $\boldsymbol{\Theta}_{01,i,j}=(-1)^{|i-j|}(0.4)^{1+|i-j|},(i,j)\in[p]^2$. Since~$s=p$ this design leads to a violation of exact/strong sparsity, but the entries decay exponentially fast in magnitude as one moves away from the
diagonal. It holds that $\rho(\boldsymbol{\Theta}_{01})=0.9$.
\item\textbf{Design \hypertarget{Design C}{C}:} The DGP is a VAR(4) with both~$\boldsymbol{\Theta}_{01}$
and~$\boldsymbol{\Theta}_{04}$  block-diagonal matrices with diagonal blocks of size $4\times4$ with all entries equal
to $0.15$ and $-0.1,$ respectively. $\boldsymbol{\Theta}_{02}=\boldsymbol{\Theta}_{03}=\mathbf{0}_{p\times p}$. Thus, $s=8$ and the spectral radius of the companion matrix is~$\rho(\widetilde{\boldsymbol{\Theta}}_{0})=0.9.$\footnote{We deviate from the $5\times5$ blocks in \citet[Experiment B, p.~333]{kock_oracle_2015}
to ensure that the block sizes are divisors of~$p$.}
\end{itemize}
 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:
\begin{itemize}
\item\textbf{Design \hypertarget{Design D}{D}:} The $\{\boldsymbol{\varepsilon}_{t}\}_{t\in\mathbb Z}$ are i.i.d.~centered Gaussian with covariance matrix  $\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon},i,j}=\sigma_{\varepsilon}^2\cdot(0.9)^{\mathbf{1}(i\neq j)},(i,j)\in[p]^2$. This covariance structures implies strong correlation among the regressors. Specifically, the largest $m$-sparse eigenvalue $\phi_{\max}(m,\boldsymbol{\Sigma}_{\boldsymbol{Y}})$ [see \eqref{eq:largestmsparseeigenvalue}] grows linearly with $m$. Thus, $\lim_{n\to\infty}\phi_{\max}(s\ln(n),\boldsymbol{\Sigma}_{\boldsymbol{Y}})=\infty$, i.e.~an assumption used to establish performance guarantees for the post-Lasso in Theorem \ref{thm:postLasso} is violated.
\item\textbf{Design \hypertarget{Design E}{E}:} The
$\{\boldsymbol{\varepsilon}_{t}\}_{t\in\mathbb Z}$ are i.i.d.~with $\boldsymbol{\varepsilon}_t\sim(\sigma_{\varepsilon}/\sqrt{5/3})\cdot t_5(\mathbf{0}_p,[(0.9)^{\mathbf{1}(i\neq j)}])$, where $t_5(\mathbf{0}_p,[(0.9)^{\mathbf{1}(i\neq j)}])$ denotes the multivariate Student distribution with five degrees of freedom, all-zero locations, and scales $(0.9)^{\mathbf{1}(i\neq j)},(i,j)\in[p]^2$. This design conflicts with Assumption \ref{assu:Innovations} in that the innovation tails are heavier than sub-Weibull.
\item\textbf{Design \hypertarget{Design F}{F}:} This design investigates the effect of conditional heteroskedasticity. Specifically, let $\{\boldsymbol{\eta}_t\}_{t\in\mathbb Z}$ be i.i.d.~with $\boldsymbol{\eta}_t\sim \mathrm{N}(\mathbf{0}_p,\mathbf{I}_p)$. Let $\boldsymbol{\varepsilon}_t = \boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}}^{1/2}(\boldsymbol{\eta}_{t-1}) \boldsymbol{\eta}_t$, where $\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}}(\boldsymbol{\eta}_{t-1}):=\mathrm{diag}(\{\sigma_{\varepsilon,i}^2(\boldsymbol{\eta}_{t-1})\}_{i\in[p]})$ with conditional standard deviations taking the form
\[
\sigma_{\varepsilon,i}(\boldsymbol{\eta}_{t-1}):=\sigma_{\varepsilon}\cdot\begin{cases}
\mathrm{e}^{-1.5|\eta_{t-1,i}|+1.5|\eta_{t-1,i+1}|}, & \text{if}\;i\in[p-1],\\
\mathrm{e}^{-1.5|\eta_{t-1,p}|+1.5|\eta_{t-1,1}|}, & \text{if}\;i=p.
\end{cases}
\]
This design violates Assumption \ref{assu:Innovations} since the $\{\boldsymbol{\varepsilon}_{t}\}_{t\in\mathbb Z}$ have conditional variances that can become arbitrarily small/large. Assumption \ref{assu:Innovations} can be enforced by censoring the~$\sigma_{\varepsilon,i}(\boldsymbol{\eta}_{t-1})$ from above and below by an arbitrarily large (small) number.
\item\textbf{Design \hypertarget{Design G}{G}:} The $\{\boldsymbol{\varepsilon}_{t}\}_{t\in\mathbb Z}$ are as in Design \hyperlink{Design D}{D} but $\boldsymbol{\Theta}_{01}=(1-5/n)\mathbf{I}_p$. This is a local to unit root design which violates Assumption \ref{assu:Companion} since $\sup_{n\in\mathbb N}\max_{i\in[p]}\|(\boldsymbol{\Theta}_{01}^{h})_{i,1:p}\|_{\ell_{2}}=\sup_{n\in\mathbb N}(1-5/n)^{h}=1$ for all $h\in\mathbb N$.
\end{itemize}


\subsection{Implementation and Performance Measure}
The Lasso, post-Lasso, and sqrt-Lasso are implemented as follows:
\begin{itemize}
    \item The weighted Lasso  $\widehat\boldsymbol{\beta}_i(\lambda_n^\ast,\widehat\boldsymbol{\Upsilon}_i^{(K)})$ in \eqref{eq:LASSOVector} tuned via Algorithm~\ref{alg:Data-Driven-Penalization}.
    \item The post-Lasso $\check\boldsymbol{\beta}_i(\lambda_n^\ast,\check\boldsymbol{\Upsilon}_i^{(K)})$ in \eqref{eq:postestim} tuned via Algorithm~\ref{alg:Data-Driven-Penalization-with-Refitting}, refitting in every step.\footnote{Theorem~\ref{alg:Data-Driven-Penalization-with-Refitting} provides performance guarantees for this algorithm that are identical to those in Theorem~\ref{thm:postLasso} for the post-Lasso which refits $\widehat\boldsymbol{\beta}_i(\lambda_n^\ast,\widehat\boldsymbol{\Upsilon}_i^{(K)})$ from the final step of Algorithm~\ref{alg:Data-Driven-Penalization} alone.}
    \item The sqrt-Lasso $\dot{\boldsymbol{\beta}}_{i}(\lambda_n^\ast/2,\dot{\boldsymbol{\Upsilon}})$ in \eqref{eq:pensqrtLasso}.
\end{itemize}
We specify $c$ and $\gamma_n$ as in \eqref{eq:tuning_c_and_gamma} 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 \texttt{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 \eqref{eq:LASSOVector} is solved using \texttt{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 \texttt{R}.}

\subsection{Results}\label{subsec: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{fig:SimulationsRelativeErrors} 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{fig:SimulationsAbsoluteErrors} of Section \ref{sec:appsim} in the appendix.

\begin{figure}\label{fig:SimulationsRelativeErrors}\caption{Average Estimation Errors relative to Weighted Lasso}
\centering{}\includegraphics[width=.90\textwidth]{img/Figure_rel_mean_max_ell2_errors_p_by_design_1000_MC_100_to_1000_n_16_to_128_p}
\end{figure}

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{thm:Rates-for-sqrtLASSO-data-driven-loadings}. 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.


\section{Empirical Illustration\label{sec: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 \citet{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 (\texttt{prepare\_missing.m}), then removing
outliers (\texttt{remove\_outliers.m}), and finally replacing missing
values with the unconditional average of the corresponding series
(as in the initialization of \texttt{factors\_em.m}).}

\subsection{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{sec:Simulations}, but we now include an unpenalized intercept in each equation to account for non-zero means of the variables.\footnote{See Section \ref{subsec:IncludingIntercepts} for a modification of Algorithm \ref{alg:Data-Driven-Penalization}, 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{fig:FRED-MD:-Forecasting-Performance} 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 \emph{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 \emph{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.}

\begin{figure}[htb]
\caption{FRED-MD: Forecasting Performance\label{fig:FRED-MD:-Forecasting-Performance}}

\centering{}\includegraphics[width=0.495\textwidth]{img/Figure_FRED_rel_mean_ivwsfe_N_120_qmax_12.png}
\includegraphics[width=0.495\textwidth]{img/Figure_FRED_rel_q95_ivwsfe_N_120_qmax_12.png}
\end{figure}


\bibliographystyle{apalike}
\bibliography{My_Library}