EconBase
← Back to paper

Regularized Estimation of High-Dimensional Vector AutoRegressions with Weakly Dependent Innovations

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

53,848 characters · 10 sections · 56 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Regularized Estimation of High-Dimensional Vector Autoregressions with Weakly Dependent Innovations

\address{Center for Statistics and Machine Learning, Princeton University} \email{[email removed] }

\address{Department of Economics, Pontifical Catholic University of Rio de Janeiro} \email{[email removed]}

\address{School of Applied Mathematics (EMap), Getulio Vargas Foundation, Rio de Janeiro (FGV-RJ)} \email{[email removed]}

abstractThere has been considerable advance in understanding the properties of sparse regularization procedures in high-dimensional models. In time series context, it is mostly restricted to Gaussian autoregressions or mixing sequences. We study oracle properties of LASSO estimation of weakly sparse vector-autoregressive models with heavy tailed, weakly dependent innovations with virtually no assumption on the conditional heteroskedasticity. In contrast to current literature, our innovation process satisfy an $L^1$ mixingale type condition on the centered conditional covariance matrices. This condition covers $L^1$-NED sequences and strong ($\alpha$-) mixing sequences as particular examples. \newline \newline JEL: C32, C55, C58. \newline \newline Keywords: high-dimensional time series, LASSO, VAR, mixing.

\onehalfspace

Introduction

Modeling multivariate time series data is an important and vibrant area of research. Applications range from economics and finance, as in cS1980, Bauer2011, Chiriac2011, or vaR2016, to air pollution and ecological studies hG2013,kbElRdP2013,mSsBkE2017. Among alternatives, the Vector Autoregressive (VAR) model is certainly one of the most successful in modeling temporal evolution of vectors, networks, and matrices. See hL1991 or cW2015 for comprehensive textbook introductions.

The advances in data collection and storage have created data sets with large numbers of time series (Big Data), where the number of model parameters to be estimated may exceed the number of available data observations. A common approach to dealing with high-dimensional data is to impose additional structure in the form of (approximate) sparsity and estimate the parameters by some shrinkage method. Examples of estimation techniques range from Bayesian estimation with “spike-and-slab” priors to sparsity-inducing shrinkage, such as the least absolute and shrinkage estimator (LASSO) and its many extensions. See sMgR2019 for a nice survey on Bayesian VARs or abKmcMgfV2019 for a review on penalized regressions applied to time-series models.

Our Contributions

In this paper we study non-asymptotic properties of high-dimensional VAR models and their parameter estimates using equation-wise (row-wise or node-wise) LASSO. We show that, with high probability, estimated and population parameter vectors are close to each other in the Euclidean norm and discuss restrictions on the rate which the number of parameters can increase as the sample size diverges.

The importance of our results relies on the fact that our non-asymptotic guarantees serve as a fundamental ingredient for the derivation of asymptotic properties of penalized estimators in high-dimensional VAR models, as in rAsSiW2020. In particular, our results apply with minimal restrictions on the conditional variance model, allowing, for instance, large-dimensional multivariate linear processes in the variance. Moreover, auxiliary results proved in this paper are of independent interested and can, for instance, be used to derive finite bounds for other type of penalization such as group/structured lasso, elastic-net, SCAD or non-convex penalties.

The data are assumed to be generated from a covariance-stationary and weakly sparse VAR model, where the innovations are martingale difference with sub-Weibull tails and conditional covariance matrix satisfying a $L^1$ mixingale assumption. An important feature is that the resulting process $\{{\bf y}_t\}$ is not necessarily mixing. Mixing assumptions can be notoriously difficult to show and we avoid it in this paper. Nevertheless, it follows that our conditions cover strong mixing innovations as a particular case.

These conditions contemplate VAR models with conditional heteroskedasticity as in lBsLjR2006,fBfFrS2011 or stochastic volatility as in sCyOmA2009.

Literature review

Some consistency results on model estimation and selection of high-dimensional VAR processes were obtained by sSpB2011, though under much stronger assumptions, such as Gaussianity. pLmW2012, sBgM2015 {and aKlC2015} developed powerful concentration inequalities that enabled them to establish consistency under weaker conditions and prove that these conditions hold with high probability. In particular, sBgM2015 established consistency of $\ell_1$-penalized least squares and maximum likelihood estimators of the coefficients of high-dimensional Gaussian VAR processes and related the estimation and prediction error to the complex dependence structure of VAR processes. Other estimation approaches, including Bayesian approaches, are discussed by raDpZtZ2016. kMpPlS2019 proposed a factor-augmented large dimensional VAR and studied finite sample properties and provide estimation results. However, they assume independent and identically distributed errors. More recently, kWaTzL2017 derived finite-sample guarantees for the LASSO in a misspecified VAR model. Authors assume the series is either $\beta$-mixing process with sub-Weibull marginal distributions or $\alpha$-mixing Gaussian processes. {Finally, rAsSiW2020 develop theoretical results for point estimation and inference in near epoch dependent time series using desparsified lasso under high level conditions on the generating process.}

Organization of the Paper

The paper is organized as follows. In Section (ref) we define the model and the main assumptions in the paper. In Section (ref) we discuss examples of applications of our results. The theoretical results are presented in Section (ref), while in Section (ref) we provide a discussion of our findings and conclude the paper. All technical proofs are relegated to the Appendix.

Notation

Throughout the paper we use the following notation. For a vector ${\bf b} = (b_1, ..., b_k)'\in\mathbb{R}^k$ and $p\in[1,\infty]$, $ |{\bf b}|_p$ denotes its $\ell_p$ norm, i.e. $|{\bf b}|_p = (\sum_{i=1}^k|b_i|^p)^{1/p}$ for $p\in[1,\infty)$ and $|{\bf b}|_\infty = \max_{1\leq i\leq k} |b_i|$. We also define $|{\bf b}|_0 = \sum_{i=1}^kI(b_i\ne 0)$. For a random variable $X$, $\|X\|_p = (\mathbb{E}|X|^p)^{1/p}$ for $p\in[1,\infty)$ and $\|X\|_\infty = \inf\{a\in\mathbb{R}:\Pr(|X|\ge a)=0\}$. For a $m\times n$ matrix $\boldsymbol{A}$ with elements $a_{ij}$, we denote ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1 = \max_{1\le j\le n}\sum_{i=1}^m|a_{ij}|$, ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty = \max_{1\le i\le m}\sum_{j=1}^n|a_{ij}|$, the induced $\ell_\infty$ and $\ell_1$ norms respectively, and the maximum elementwise norm ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\max} = \max_{i,j} |a_{ij}|$. Also $\Lambda_{\min}(\boldsymbol A)$ and $\Lambda_{\max}(\boldsymbol A)$ denotes the minimum and maximum eigenvalues of the square matrix $\boldsymbol A$, respectively.

Model setup and Assumptions

Let $\{{\bf y}_t = (y_{t,1},...,y_{t,n})':t\in\mathbb{Z}\}$ be a vector stochastic process defined in some fixed probability space taking values on $\mathbb{R}^n$ given by

equation[equation omitted — 124 chars of source]

where ${\bf u}_t = (u_{t,1},...,u_{t,n})'$ is a zero-mean vector of innovations and $\boldsymbol{A}_1,\ldots,\boldsymbol{A}_p$, are $n\times n$ parameter matrices. The dimension $n:= n_T$ and order $p:= p_T$ of the process are allowed to increase with the number of observations $T$. Write the vector-autoregressive (VAR) process (ref) using its first-order representation:

equation[equation omitted — 118 chars of source]

where $\tilde{{\bf y}}_t = ({\bf y}_t',...,{\bf y}_{t-p+1}')'$, $\tilde{{\bf u}}_t = ({\bf u}_t',\boldsymbol{0}',...,\boldsymbol{0}')'$, and \[ \boldsymbol{F}_T = \left[

array[array omitted — 484 chars of source]

\right]. \]

Consider now the following assumptions.

assumption[A1] All roots of the reverse characteristic polynomial $\boldsymbol{\mathcal{A}}(z) = \boldsymbol{I}_n-\sum_{i=1}^p\boldsymbol{A}_jz^j$ lie outside the unit disk for each $p,n\in\mathbb{N}$, and there exist $\bar{c}_\Phi >0 $, $c_\phi>0 $ and $0<\gamma_1\le 1$ such that for all $m\in\mathbb{N}$ \begin{equation} \sum_{k=m}^\infty\left|\boldsymbol\phi_{k,i}\right|_1 \le \bar{c}_\Phi e^{-c_\phi m^{\gamma_1}}, \end{equation} uniformly in $1\leq i\leq n$, where $\boldsymbol\Phi_k := \boldsymbol J'\boldsymbol{\boldsymbol F}_T^k\boldsymbol J = (\boldsymbol\phi_{k,1},...,\boldsymbol\phi_{k,n})'$, $\boldsymbol{F}_T$ denote the companion matrix and $\boldsymbol J = (\boldsymbol{I}_n,\boldsymbol{0}_n,...,\boldsymbol{0}_n)'$.
assumption[A2] The sequence $\{({\bf u}_t,\mathcal{F}_t)\}_t$ is a covariance stationary (for each $T\in\mathbb{N}$) martingale difference process where the filtration $\{\mathcal{F}_t\}_t$ includes the natural filtration of $\{{\bf u}_t\}$ . The smallest and largest eigenvalues of $\boldsymbol\Sigma:=\mathbb{E}({\bf u}_1 {\bf u}_1')$ are bounded away from $0$ and $\infty$ respectively, uniformly in $T\in\mathbb{N}$. Furthermore, for all ${\bf b}_1,{\bf b}_2 \in \{\boldsymbol v\in\mathbb{R}^n:|\boldsymbol v|_1\le 1\}$ and $ m\in\mathbb{N}$, \[ \mathbb{E}\left|\mathbb{E}[{\bf b}_1'({\bf u}_t{\bf u}_t' - \boldsymbol\Sigma){\bf b}_2|\mathcal{F}_{t-m}]\right| \le a_1e^{-a_2m^{\gamma_2}}, \] for some $a_1,a_2>0$ and $0<\gamma_2\le 1$, uniformly in $1\leq t \leq T$ and $T\in\mathbb{N}$.
assumption[A3] For all ${\bf b}\in\{ \boldsymbol v \in \mathbb{R}^n:|\boldsymbol v|_1\le 1\}$ and all $0<x<\infty$, $\Pr(|{\bf b}'{\bf u}_t|>x)\le 2 e^{-|x/c_\alpha|^\alpha}$ for some $\alpha>0$, $0<c_\alpha<\infty$, uniformly in $1\leq t \leq T$ and $T\in\mathbb{N}$.

Assumption (A1) requires that the VAR process is stable and admits an infinite-order vector moving average, VMA($\infty$), representation for all $n$ and $p$ as

equation[equation omitted — 191 chars of source]

Furthermore, the coefficients of the MA($\infty$) representations of each $\{y_{i,t}\}$, $i=1,..,n$, are absolutely summable with exponentially decaying rate. This condition is satisfied in standard VAR($p$) models, where $n$ and $p$ are fixed. In models that $n$ is large, Lemma (ref) in Appendix (ref) shows that condition (ref) is satisfied if $\sum_{k=1}^p{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A}_k \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty<1$ and further regularity conditions on the size of the coefficients. Finally, notice that under (A1) it is also true that $\max_{k,i}|\boldsymbol\phi_{k,i}|_\infty \le \bar{c}_\Phi$, which means that the coefficients $\{\boldsymbol\Phi_k\}$ are uniformly upper bounded under the maximum entry-wise norm.

Assumption (A2) requires the error process to be a martingale difference process and satisfy a very weak dependence condition on its conditional variance. The former restricts the model to be correctly specified in the mean. Nevertheless, this assumption is standard in the literature and we are able to derive results covering a broad range of data generating processes and conditional dependence measures. The latter is the $L^1$ projective dependence measure appearing in weakdependence. Note that (1) strong mixing (or $\alpha$-mixing) sequences with exponential decay of the mixing coefficient satisfy this condition jD1994; and (2) uniform mixing sequences ($\phi$-mixing) and $\beta$-mixing sequences are also strong mixing, but the converse is not true rB2005. If we denote the centered outer product series $\boldsymbol{v}_t = \mathsf{vech}\,({\bf u}_t{\bf u}_t'-\boldsymbol\Sigma)$, this assumptions requires that $\{\boldsymbol{v}_t\}$ is $L^1$ mixingale. It means that stochastic process with $L^r$ bounded, $L^1$ near-epoch dependent, centered outer product series $\boldsymbol{v}_t$ are also contemplated in this setting dA1988. Finally, Assumptions (A1) and (A2) combined ensure that $\{{\bf y}_t\}$ is second order stationary for each $n$ and $p$ hL2006.

Condition (A3) imposes restrictions on the tail behavior of the innovation process $\{{\bf u}_t\}$ that are shared by $\{{\bf y}_t\}$. More precisely, we impose moment conditions on all linear combinations ${\bf b}'{\bf u}_t$. Lemma (ref), in the appendix, shows that each $\{y_{i,t}\}$ ($i=1,...,n$) also share the same tail properties of $\{{\bf u}_t\}$. This condition is essential for defining the rate in which $n$ and $p$ increase with $T$. We focus on the case the tail decays at rate $O(e^{-cx^\alpha})$ for some $\alpha>0$, that is, $\{{\bf b}'{\bf u}_t\}$ is sub-Weibull with parameter $\alpha$ studied in kWaTzL2017. Note that when $\alpha\ge 1$ and $\alpha\ge 2$ we have the sub-exponential and sub-Gaussian tails respectively. However, when $\alpha\in(0,1)$ the moment generating function does not exist at any point and and these variables are usually called heavy tailed.

It is convenient to write the model in stacked form. Let ${\bf x}_t = ({\bf y}_{t-1}',\ldots,{\bf y}_{t-p}')'$ be the $np\times 1$ vector of regressors and ${\bf X} = ({\bf x}_1,...,{\bf x}_T)'$ the $T\times np$ matrix of covariates. Let ${\bf Y}_i = (y_{i,1},...,y_{i,T})'$ be the $T\times 1$ vector of observations for the $i^{th}$ element of ${\bf y}_t$, and ${\bf U}_i = (u_{i,1},...,u_{i,T})'$ the corresponding vector of innovations. Denote $\boldsymbol\beta_i$ the $np\times 1$ vector of coefficients corresponding to equation $i$. Then, model (ref) is equivalent to

equation[equation omitted — 107 chars of source]

We now make additional assumptions concerning model (ref).

assumption[A4] The true parameter vectors $\boldsymbol\beta_i$, $i=1, \ldots, n$, satisfy $\sum_{j=1}^{np}|\beta_{i,j}|^q \le R_q$ for some $0\le q<1$ and $0<R_q<\infty$ where $R_q:=R_{q,T}$ is allowed to depend on the sample size $T$.
assumption[A5] For each $T\in\mathbb{N}$, the smallest eigenvalue of $\boldsymbol\Gamma := T^{-1}\mathbb{E}({\bf X}'{\bf X})$ is greater than a positive constant $\sigma_\Gamma^2$ that might depend on $T$.

Assumption (A4) imposes weak sparsity of the coefficients, in a sense that most of them are small. This condition is slightly stronger than we need in a sense that we may have distinct $q_i$ and $R_{q,i}$ for each equation. In the case $q=0$ we have sparsity in the standard sense, meaning that $R_0 = s$, the number of non-zero coefficients. In practice, we estimate a sparse model that truncates all coefficients close to zero. This assumption is standard for weak sparsity, see sNpRmWbY2012[section 4.3] and yHrT2019[Assumption 1] for an application in time series setting.

Assumption (A5) is often used in the sparse estimation literature aKlC2015,mcMeM2015,yHrT2019. sBgM2015 (Proposition 2.3) derived bounds for $\Lambda_{\min}(\boldsymbol\Gamma)$ and $\Lambda_{\max}(\boldsymbol\Gamma)$ using properties of the block Toeplitz matrix $\boldsymbol\Gamma$ and its generating function, the cross-spectral density of the generating VAR$(p)$ process:

equation[equation omitted — 376 chars of source]

where $\boldsymbol{\mathcal{A}}^*$ is the conjugate transpose of $\boldsymbol{\mathcal{A}}$, the reverse characteristic polynomial, defined in Assumption (A1). sBgM2015[Proposition 2.2] shows that under (A1), \[ {\max_{|z|=1}\Lambda_{\max}(\boldsymbol{\mathcal{A}}^*(z)\boldsymbol{\mathcal{A}}(z))} <\left[1+\frac{\sum_{k=1}^p({\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A}_k \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1+{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A}_k \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty)}{2}\right]^2. \] Hence, (A5) is satisfied if, for instance, $\Lambda_{\min}(\boldsymbol\Sigma)>0$, $\sum_{k=1}^p{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A}_k \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1<\infty$ and $\sum_{k=1}^p{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol{A}_k \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty<\infty$.

Examples

In this section we illustrate processes satisfying Assumptions (A2) and (A3). In the first two examples we discuss sufficient conditions involving mixing and near epoch dependent sequences, traditionally found in the literature. In the final two examples, we discuss variance process admitting an AR($\infty$) representation.

example[Strong mixing sequences] Let $\{{\bf u}_t\}$ denote a martingale difference, strong mixing sequence with coefficients $\alpha_m< b_1\exp(-b_2m^{\gamma_2})$ and common covariance matrix $\boldsymbol\Sigma$ with eigenvalues bounded away from zero and infinity, uniformly in $n$. It follows that $r_t = {\bf b}_1'{\bf u}_t{\bf u}_t'{\bf b}_2$ is also strong mixing of same size and, from jD1994, $\mathbb{E}[r_t-\mathbb{E}(r_t)|\mathcal{F}_{t-m}]\le a_1 \exp(-a_2m^{\gamma_2})$, for constants $a_1$ and $a_2$.
example[$L^1$ near-epoch dependent process] Let $\{{\bf u}_t\}$ denote a weakly stationary, martingale difference sequence. Suppose ${\bf b}'\boldsymbol{v}_t = {\bf b}'\mathsf{vech}\,({\bf u}_t{\bf u}_t'-\boldsymbol\Sigma)$ is a centered, $L^1$-NED sequence on $\mathcal{F}_t = \sigma\langle \epsilon_t, \epsilon_{t-1}, ...\rangle$, where $\{\epsilon_t\}$ is $\alpha$-mixing with coefficients $\alpha_m \le c_1 \exp(c_2 m^{\gamma_1})$, for all ${\bf b}\in\{{\bf b}\in\mathbb{R}^{n(n+1)/2}:|{\bf b}|_1\le 1\}$. It means that there are finite constants $\{d_t\}$ and $\{\psi_m\}$ such that \[ \mathbb{E}\left| {\bf b}'(\boldsymbol{v}_t-\mathbb{E}[\boldsymbol{v}_t|\mathcal{F}_{t-m:t}])\right| \le d_t\psi_m, \] where $\mathcal{F}_{t-m:t} = \sigma\langle \epsilon_t,...,\epsilon_{t-m}\rangle$ and $\psi_m \le \exp(c_3 m^{\gamma_2})$. Under Assumption (A3), it follows from kWaTzL2017 and Hölder inequality that for any $r<\infty$ \[ \|{\bf b}'\boldsymbol{v}_t\|_r\le |{\bf b}|_1^r\max_{1\le i\le j\le n}\|u_{it}u_{jt}\|_r\le\max_{1\le i\le n}\|u_{it}\|_{2r}\le c_4 r^{1/\alpha}. \] Finally, it follows from dA1988 that Assumption (A2) holds with $a_1 \ge (2\max_t d_t + c_4 r^{1/\alpha})(e^{c_3/2^{\gamma_2}}+6c_1e^{c_2(r-1)/r2^{\gamma_2}})$ and $a_2 \le (c_3\wedge c_2 (r-1)/r)/2^{\gamma_2}$.
example[Linear process in the variance] Let $\{\boldsymbol{v}_t,\boldsymbol{\epsilon}_t\}$ denote a sequence of centered independently and identically distributed, sub-Weibull random variables with parameter (at least) $2\alpha$, taking values in $\mathbb{R}^{2n}$ with identity covariance matrix. Let ${\bf u}_t = \boldsymbol{H}_t^{1/2}\boldsymbol{v}_t$ where $\boldsymbol{H}_t^{1/2}$ is the lower diagonal Cholesky decomposition of $\boldsymbol{H}_t$ and \[ \boldsymbol{h}_t = \mathsf{vech}\,(\boldsymbol{H}_t) = \boldsymbol c + \sum_{j=1}^\infty \boldsymbol \Psi_j\boldsymbol\eta_{t-j}. \] Here, $\mathsf{vech}\,(\boldsymbol{M})$ stacks the lower diagonal elements of matrix $\boldsymbol{M}$, $\boldsymbol c$ is a vector of constants and $\boldsymbol\eta_{t} = \mathsf{vech}\,(\boldsymbol{\epsilon}_t\boldsymbol{\epsilon}_t')$. For all $\tilde{{\bf b}}\in\{{\bf b}\in \mathbb{R}^{n(n+1)/2}:|{\bf b}|_1 \le 1\}$, $\{\boldsymbol \Psi_j\}$ satisfy $\sum_{j=m}^\infty|\tilde{{\bf b}}'\boldsymbol \Psi_j|_1\lesssim e^{-a_2 m^{\gamma_2}}$. We first show $\{{\bf u}_t\}$ is weakly stationary martingale difference with respect to $\mathcal{F}_{t-1} = \sigma\langle(\boldsymbol{v}_{t-j},\boldsymbol{\epsilon}_{t-j}): j=1,2,...\rangle$. The process $\{{\bf u}_t\}$ is $\mathcal{F}_{t}$ measurable and satisfy $\mathbb{E}[{\bf u}_t|\mathcal{F}_{t-1}] = \boldsymbol{H}_t^{1/2}\mathbb{E}[\boldsymbol{v}_t|\mathcal{F}_{t-1}] = \boldsymbol 0$. Its covariance matrix is \begin{align*} \mathbb{E}[{\bf u}_t{\bf u}_t'] = \mathbb{E}[\boldsymbol{H}_t^{1/2}\mathbb{E}(\boldsymbol{v}_t\boldsymbol{v}_t'|\mathcal{F}_{t-1}) (\boldsymbol{H}_t^{1/2})']=\mathbb{E}[\boldsymbol{H}_t]. \end{align*} Now, $\mathbb{E}[\boldsymbol{h}_t] = \boldsymbol c +\sum_{j=1}^\infty\boldsymbol\Psi_j\mathbb{E}\boldsymbol\eta_{t-j} = \boldsymbol c +\sum_{j=1}^\infty\boldsymbol\Psi_j\mathsf{vech}\,(\boldsymbol{I}_n)=\boldsymbol\Sigma$, where $\mathbb{E}(\boldsymbol\eta_t) = \mathsf{vech}\,(\mathbb{E}(\boldsymbol{\epsilon}_t\boldsymbol{\epsilon}_t')) = \mathsf{vech}\,(\boldsymbol{I}_n)$ for all $t$. For constant vectors ${\bf b}_1,{\bf b}_2\in\{{\bf b}\in \mathbb{R}_n:|{\bf b}|\le 1\}$, \begin{align*} \mathbb{E}[{\bf b}_1'({\bf u}_t{\bf u}_t'-\boldsymbol\Sigma){\bf b}_2|\mathcal{F}_{t-1}] &= {\bf b}_1'(\boldsymbol{H}_t-\mathbb{E} \boldsymbol{H}_t){\bf b}_2\\ &= \tilde{{\bf b}}'(\boldsymbol{h}_t-\mathbb{E} \boldsymbol{h}_t) \\ &= \sum_{j=1}^\infty\tilde{{\bf b}}'\boldsymbol\Psi_j(\boldsymbol\eta_{t-j}-\mathbb{E} \boldsymbol\eta_{t-j}), \end{align*} where $\tilde{{\bf b}}\in\{{\bf b}\in \mathbb{R}^{n(n+1)/2}:|{\bf b}|_1 \le 1\}$. It follows that \begin{align*} \mathbb{E}\left|\mathbb{E}[{\bf b}_1'({\bf u}_t{\bf u}_t'-\boldsymbol\Sigma){\bf b}_2|\mathcal{F}_{t-m}]\right| &= \mathbb{E}\left|\sum_{j=1}^\infty\tilde{{\bf b}}'\boldsymbol\Psi_j\mathbb{E}(\boldsymbol\eta_{t-j}-\mathbb{E} \boldsymbol\eta_{t-j}|\mathcal{F}_{t-m})\right|\\ &=\left\|\sum_{j=m}^\infty\tilde{{\bf b}}'\boldsymbol\Psi_j(\boldsymbol\eta_{t-j}-\mathbb{E} \boldsymbol\eta_{t-j})\right\|_1\\ &\le 2 \left(\sum_{j=m}^\infty|\tilde{{\bf b}}'\boldsymbol\Psi_j|_1\right)\max_{|{\bf b}|_1\le 1}\|{\bf b}'\boldsymbol{\epsilon}_t\|_2^2, \end{align*} where in the last line we use the same arguments of Lemma (ref) in the appendix, followed by the triangle inequality. Then, Assumption (A2) is satisfied under the condition that $\sum_{j=m}^\infty|\tilde{{\bf b}}'\boldsymbol\Psi_j|_1\lesssim e^{-a_2 m^{\gamma_2}}$ and $\|{\bf b}'\boldsymbol{\epsilon}_t\|_2\le c_2<\infty$. It follows from kWaTzL2017 that $\{\boldsymbol u_t\}$ is sub-Weibull with parameter $\alpha$ if $\operatorname*{\mathsf{sup}}_{d\ge_1}d^{-1/\alpha}\|{\bf b}'{\bf u}_t\|_d\leq c_\alpha<\infty$. For any $d\ge 1$, \begin{align*} d^{-1/\alpha}\|{\bf b}'{\bf u}_t\|_d &= d^{-1/\alpha}\|{\bf b}'{\bf u}_t{\bf u}_t'{\bf b}\|_{d/2}^{1/2}\\ &= d^{-1/\alpha}\|{\bf b}'\boldsymbol{H}_t^{1/2}\boldsymbol{v}_t\boldsymbol{v}_t'(\boldsymbol{H}_t^{1/2}){\bf b}\|_{d/2}^{1/2}\\ &= d^{-1/\alpha}\left\|{\bf b}'\boldsymbol{H}_t{\bf b}\times \frac{{\bf b}'\boldsymbol{H}_t^{1/2}\boldsymbol{v}_t\boldsymbol{v}_t'(\boldsymbol{H}_t^{1/2}){\bf b}}{{\bf b}'\boldsymbol{H}_t{\bf b}}\right\|_{d/2}^{1/2}\\ &\le d^{-1/\alpha}\left\{\mathbb{E}\left(\left|{\bf b}'\boldsymbol{H}_t{\bf b}\right|^{d/2}\mathbb{E}\left[\left|\frac{{\bf b}'\boldsymbol{H}_t^{1/2}\boldsymbol{v}_t\boldsymbol{v}_t'(\boldsymbol{H}_t^{1/2}){\bf b}}{{\bf b}'\boldsymbol{H}_t{\bf b}}\right|^{d/2}\Big\vert\mathcal{F}_{t-1}\right]\right)\right\}^{1/d}\\ &\le d^{-1/\alpha}\left\{\mathbb{E}\left(\left|{\bf b}'\boldsymbol{H}_t{\bf b}\right|^{d/2}\operatorname*{\mathsf{sup}}_{\boldsymbol\delta'\boldsymbol\delta = 1}\mathbb{E}[(\boldsymbol\delta'\boldsymbol{v}_t\boldsymbol{v}_t' \boldsymbol\delta)^{d/2}|\mathcal{F}_{t-1}]\right)\right\}^{1/d}\\ &= d^{-1/2\alpha}\left\|{\bf b}'\boldsymbol{H}_t{\bf b}\right\|_{d/2}^{1/2}\operatorname*{\mathsf{sup}}_{\boldsymbol\delta'\boldsymbol\delta = 1}d^{-1/\alpha}\|\boldsymbol\delta'\boldsymbol{v}_t\|_d,\\ &\le \left(d^{-1/\alpha}\left\|{\bf b}'\boldsymbol{H}_t{\bf b}\right\|_{d/2}\right)^{1/2} c_{2\alpha} \end{align*} where the two last lines follow because $\{\boldsymbol{v}_t\}$ is independent and sub-Weibull process with parameter $2\alpha$. There is a $\tilde{{\bf b}}\in\{{\bf b}\in\mathbb{R}^{n(n-1)/2}:|{\bf b}|_1\le 1\}$ such that \begin{align*} d^{-1/\alpha}\left\|{\bf b}'\boldsymbol{H}_t{\bf b}\right\|_{d/2} &= d^{-1/\alpha}\|\tilde{{\bf b}}'\boldsymbol{h}_t\|_{d/2}\\ &= d^{-1/\alpha}\left\|\tilde{{\bf b}}'\boldsymbol c+\sum_{j=1}^{\infty}\tilde{{\bf b}}'\boldsymbol\Psi_j\boldsymbol\eta_{t-j}\right\|_{d/2}\\ &\le d^{-1/\alpha}\left|\tilde{{\bf b}}'\boldsymbol c\right|_{d/2}+d^{-1/\alpha}\left\|\sum_{j=1}^{\infty}\tilde{{\bf b}}'\boldsymbol\Psi_j\boldsymbol\eta_{t-j}\right\|_{d/2}\\ &\le d^{-1/\alpha}\left|\tilde{{\bf b}}'\boldsymbol c\right\|t_{d/2}+ 2 \left(\sum_{j=1}^\infty|\tilde{{\bf b}}'\boldsymbol\Psi_j|_1\right)\max_{|{\bf b}|_1\le 1}\left(d^{-1/2\alpha}\|{\bf b}'\boldsymbol{\epsilon}_t\|_{d}\right)^2\\ &\le \left|\tilde{{\bf b}}'\boldsymbol c\right|_{d/2} + c_{2\alpha}^2. \end{align*} Combining these bounds, process $\{{\bf u}_t\}$ is sub-Weibull with parameter $\alpha$, satisfying Condition (A3).
example[{Stochastic covariance}] Let \[ {\bf y}_t = \sum_{i=1}^p\boldsymbol{A}_i{\bf y}_{t-i} + \boldsymbol{H}_t^{1/2}\boldsymbol{v}_t, \quad \boldsymbol{H}_{t+1} = \boldsymbol{C}_0 + \boldsymbol{\Psi} \boldsymbol{H}_t \boldsymbol{\Psi}' + \boldsymbol{\epsilon}_t\boldsymbol{\epsilon}_t', \] where $\boldsymbol{v}_t\overset{iid}{\sim}\mathsf{N}(\boldsymbol 0,\boldsymbol{I}_n)$, $\boldsymbol{\epsilon}_t \overset{iid}\sim (\boldsymbol 0,\boldsymbol{I}_n)$ is sub-Weibull with parameter $2\alpha$, and processes $\{\boldsymbol{v}_t\}$ and $\{\boldsymbol{\epsilon}_t\}$ are independent. Let $\mathcal{F}_t = \sigma\langle(\boldsymbol{v}_{t-j},\boldsymbol{\epsilon}_{t-j}):j=0,1,2,\ldots\rangle$. Process $\{\boldsymbol{H}_t\}$ is a matrix process characterizing the stochastic covariance of $\boldsymbol{u}_t = \boldsymbol{H}_t^{1/2}\boldsymbol{v}_t$ evolves according to a matrix autoregressive process. The intercept $\boldsymbol{C}_0$ is symmetric and positive definite matrix and the eigenvalues of $\boldsymbol{\Psi}$ are inside the unity circle. We also assume ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty < 1$, which implies that the largest singular value of $\boldsymbol\Psi$ is smaller than one. Under these conditions, the process $\{\boldsymbol{H}_t\}$ is ensured to be positive definite and stationary. To see the latter, vectorize the process to obtain $\boldsymbol{h}_t = \mathsf{vec}\,(\boldsymbol{H}_t)$, $\boldsymbol{c}_0 = \mathsf{vec}\,(\boldsymbol{C}_0)$, $\boldsymbol\eta_t = \mathsf{vec}\,(\boldsymbol\epsilon_t\boldsymbol\epsilon_t')$, $\bar{\boldsymbol\Psi} = \boldsymbol\Psi \otimes \boldsymbol\Psi'$ and \[ (\boldsymbol{I}-\bar{\boldsymbol{\Psi}} L)\boldsymbol{h}_{t+1} = \boldsymbol{c}_0 + \boldsymbol\eta_t \Leftrightarrow \boldsymbol{h}_{t+1} = \sum_{j=0}^\infty\bar{\boldsymbol\Psi}^j\boldsymbol{c}_0 + \sum_{j=0}^\infty\bar{\boldsymbol\Psi}^j\boldsymbol\eta_{t-j}. \] The process is stationary because eigenvalues of $\bar{\boldsymbol\Psi}$ are products of eigenvalues of $\boldsymbol\Psi$, which is also inside the unity circle. We are exactly in the setting of previous example. We have to show $\sum_{j=m+1}^\infty|{\bf b}'\bar{\boldsymbol\Psi}^j|_1 \lesssim e^{-a_2m^{\gamma_2}}$ for all ${\bf b}\in R^{n^2}:|{\bf b}|_1 = 1$. The term on the left hand side is bounded by $\sum_{j>m}{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \bar{\boldsymbol\Psi}^j \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1$ and each ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \bar{\boldsymbol\Psi}^j \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1\le{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi^j \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi^j \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty\le ({\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty)^{j}$. Under the assumption that ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Psi \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty = c_{\max} < 1$, $\sum_{j=m+1}^\infty|{\bf b}'\bar{\boldsymbol\Psi}^j|_1 \le \frac{c_{\max}}{1-c_{\max}}e^{-m\log(1/c_{\max})}$ meaning that condition is satisfied with $\gamma_2 = 1$ and $a_2 = \log(1/c_{\max})$. Finally, the unconditional covariance of ${\bf u}_t$ is \[ \boldsymbol{\Sigma} := \mathbb{E}[\boldsymbol H_t] = \sum_{j=0}^\infty \boldsymbol\Psi^j(\boldsymbol C_0+\boldsymbol I_n){\boldsymbol\Psi'}^j. \] The smallest eigenvalue $\min_{\boldsymbol\delta'\boldsymbol\delta = 1}\boldsymbol\delta'\boldsymbol\Sigma\boldsymbol\delta \ge \Lambda_{\min}(\boldsymbol C_0) + 1$ and largest eigenvalue of $\max_{\boldsymbol\delta'\boldsymbol\delta = 1}\boldsymbol\delta'\boldsymbol\Sigma\boldsymbol\delta \le (\rho(\boldsymbol C_0)+1)(1-\rho(\boldsymbol\Psi)^2)^{-1}$, where $\Lambda_{\min}(\boldsymbol A)$ is the smallest eigenvalue of $\boldsymbol A$ and $\rho(\boldsymbol A)$ is the spectral radius of $\boldsymbol A$.

LASSO estimation bounds

Let $\mathcal{L}_T(\boldsymbol\beta_i) = \frac{1}{T}|{\bf Y}_i-{\bf X}\boldsymbol\beta_i|_2^2$ denote the empirical squared risk, for each $i=1,...,n$. We estimate $\boldsymbol\beta_i$, $i=1,\ldots,n$, equation-wise using the LASSO procedure

equation[equation omitted — 220 chars of source]

where $\lambda_i$ are positive regularization parameters. For ease of exposition we assume $\lambda_1 = \cdots = \lambda_n = \lambda$. It is well known that $\boldsymbol\beta^*_i = \arg\min_{\boldsymbol\beta_i}\mathbb{E}\left\{\mathcal{L}_T(\boldsymbol\beta_i)\right\}$ are the population parameters in (ref), under stated conditions.

We follow the steps in sNpRmWbY2012 to derive error bounds for the equation-wise LASSO estimator. First define the pair of subspaces $\mathcal{M}(S) = \{\boldsymbol{u}\in\mathbb{R}^{np}|u_i=0, i\in S^c\}$ and its orthogonal complement $\mathcal{M}^\perp(S)= \{\boldsymbol{u}\in\mathbb{R}^{np}|u_i=0, i\in S\}$, where $S\subseteq\{1,\ldots,np\}$. Set $\boldsymbol{u}_\mathcal{M}$ and $\boldsymbol{u}_{\mathcal{M}^\perp}$ the projection of $\boldsymbol{u}$ on $\mathcal{M}(S)$ and $\mathcal{M}^\perp(S)$, respectively. Clearly, for any $\boldsymbol{u}\in\mathbb{R}^{np}$, $|u|_1 = |\boldsymbol{u}_\mathcal{M}|_1+|\boldsymbol{u}_{\mathcal{M}^\perp}|_1$. We say $|\cdot|_1$ is decomposable with respect to the pair $(\mathcal{M}(S),\mathcal{M}^\perp(S))$ for any set $S\subset\{1,\ldots,np\}$.

We have to show two conditions to obtain a finite sample estimation error bound for the parameter vectors. The first condition is known as restricted strong convexity (RSC) and restricts the geometry of the loss function around the optimum $\boldsymbol\beta^*$ and is related to the Restricted Eigenvalue vGpB2009. The second condition is known as deviation bound and restricts the size of the $\operatorname*{\mathsf{sup}}$-norm of the gradient $\nabla\mathcal{L}_T(\boldsymbol\beta^*)$. These conditions are shown to be satisfied in a set of large probability defined in Propositions (ref) and (ref).

definition[Deviation Bound (DB)] The deviation bound condition holds when the regularization parameter $\lambda$ satisfies $\{\lambda\ge 2|{\bf X}'\boldsymbol{U}_i/T|_\infty\}$ for all $i=1,...,n$.

Note that one may adopt individual $\lambda_i$s for each equation, in which above definition should be modified adequately.

definition[Restricted Strong Convexity (RSC)] Define $\mathbb{C}(\boldsymbol\beta^*,\mathcal{M},\mathcal{M}^\perp) = \{\boldsymbol\Delta\in\mathbb{R}^{np}||\boldsymbol\Delta_{\mathcal{M}^\perp}|_1\le 3|\boldsymbol\Delta_\mathcal{M}|_1+4|\boldsymbol\beta^*_{\mathcal{M}^\perp}|_1\}$. The restricted strong convexity holds for parameters $\kappa_\mathcal{L}$ and $\tau_\mathcal{L}$ if for any $\boldsymbol\Delta\in\mathbb{C}$, \[ \frac{\boldsymbol\Delta'{\bf X}'{\bf X}\boldsymbol\Delta}{T} \ge \kappa_\mathcal{L}|\boldsymbol\Delta|_2^2 - \tau_\mathcal{L}^2(\boldsymbol\beta^*). \]

sNpRmWbY2012[Section 4] show these conditions are satisfied by many loss functions and penalties. sBgM2015 show that both DB and RSC are satisfied by Gaussian VAR($p$) models in high dimensions.

If both DB and RSC hold with large probability, sNpRmWbY2012[Theorem 1] provides an $\ell_2$ estimation bound for $\widehat{\boldsymbol\beta_i}$. Our goal is to show that the error bounds are valid for each $\boldsymbol\Delta_i = \widehat{\boldsymbol\beta_i} - \boldsymbol\beta_i^*$, $i=1,\ldots,n$ at the same time.

Lemma (ref) characterizes the solutions of the optimization program in (ref). We require further notation. Define $\mathbb{C}_i := \mathbb{C}(\boldsymbol\beta_i^*,\mathcal{M}_{i,\eta},\mathcal{M}_{i,\eta}^\perp)$ for a pair of subsets $\mathcal{M}_{i,\eta} = \mathcal{M}(S_{i,\eta})$ and $\mathcal{M}^\perp_{i,\eta} = \mathcal{M}^\perp(S_{i,\eta})$, where $S_{i,\eta} = \{j\in\{1,...,pn\}||\beta_{i,j}|>\eta\}$ and $S_{i,\eta}^c = \{j\in\{1,...,pn\}||\beta_{i,j}|\le\eta\}$. These sets represent the active parameters under weak sparsity. In Theorem (ref) we set $\eta = \lambda/\sigma_\Gamma^2$ to derive our results.

lemmaSuppose $\{{\bf y}_t\}$ is generated from (ref) and Assumptions (A1), (A2) and (A3) are satisfied. Set \begin{equation} \lambda>\tau^*(\epsilon+\log(Tn^2p))^{2/\alpha}\sqrt{\frac{\epsilon+\log(n^2p)}{T}}, \end{equation} where $\epsilon>0$ and $\tau^*>0$ depends on $\tau$, $\alpha$ and $\bar{c}_\Phi$. Then, if $T>\epsilon+\log(n^2p)$, the event $\left\{\forall i=1,...,n:~\widehat{\boldsymbol\beta}_i - {\boldsymbol\beta}^* \in \mathbb{C}_i \right\}$ holds with probability at least $1-10e^{-\epsilon}$.

Lemma (ref) shows that under restrictions on $\lambda$ the solutions to the optimization program in (ref) lie inside the star-shaped sets $\mathbb{C}_i$ with high probability, as the sample size increases. It restricts the directions in which we should control the variation of our estimators. Next result shows the deviation bound holds with high probability for appropriate choice of $\lambda$. To formalize the idea, let

equation[equation omitted — 166 chars of source]

denote the event “DB holds for equation $i$ with regularization parameter $\lambda$.”

proposition[Deviation Bound] Suppose that $\{{\bf y}_t\}$ is generated from (ref), Assumptions (A1), (A2) and (A3) are satisfied and $T>\epsilon+\log(n^2p)$ for some $\epsilon>0$. Set the penalty parameters $\lambda$ as in (ref). Then, \[\Pr\left(\bigcup_{i=1}^n\mathcal{D}_i^c(\lambda)\right) \le \pi_1(\epsilon):=10e^{-\epsilon}. \]

Suppose $\epsilon = \log(np)$, $n^2p>T$. The regularization parameter $\lambda$ satisfies \[ \lambda \gtrsim [\log(np)]^{2/\alpha}\sqrt{\frac{\log(np)}{T}}, \] and $\pi_1(\lambda) \propto 1/n^2p$. This regularization parameter is $O([\log(np)]^{2/\alpha})$ larger, in rate, than one obtained in kWaTzL2017. Their results relied heavily in ${\bf y}_t$ being a $\beta$-mixing sequence in a sense that the concentration inequality derived in fMmPeR2011 depends on it. In our case, the dependence is characterized by the conditional variance of the innovation process and coefficients $\Phi_1,\Phi_2,...$, and we are not aware of "tight" concentration inequalities that hold under these assumptions. Nevertheless, for fixed $n$, it is possible to show that the concentration inequality for sub-Weibull martingales in Lemma (ref) is tight xFiGqL2012large.

Let $\boldsymbol\Gamma_T = {\bf X}'{\bf X}/T$ denote the scaled Gram matrix and $\boldsymbol{\Gamma}$ its expected value. We show that if each element in $\boldsymbol\Gamma_T$ is sufficiently close to its expectation, and Assumptions (A4) and (A5) hold, then RSC is satisfied with high probability.

lemma[Restricted Strong Convexity] Suppose Assumptions (A4) -- (A5) hold and that ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol\Gamma_T-\boldsymbol\Gamma \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\max}\leq \frac{\sigma^2_\Gamma\eta^q}{64R_q}$. Then, for any $\boldsymbol\Delta_i\in\mathbb{C}_i$ for $1\leq i\leq n$ \begin{equation} \boldsymbol\Delta_i'\boldsymbol\Gamma_T\boldsymbol\Delta_i\ge \frac{\sigma_\Gamma^2}{2}|\boldsymbol\Delta_i|_2^2 - \frac{\sigma_\Gamma^2}{2}R_q\eta^{2-q}. \end{equation}

To show RSC holds with high probability for all $i=1,..,n$ at the same time, we have to bound the event

equation[equation omitted — 239 chars of source]

where $a = \frac{\sigma^{2(1-q)}_\Gamma\lambda^q}{64R_q}$. If we assume distinct $R_{q,i}$ and $q_i$ for each equation, we should work with $\cap_i\mathcal{B}_i$ and $\mathcal{B}_i$ defined accordingly.

propositionSuppose Assumptions (A1), (A2) and (A3) hold. If \[ p< \frac{T^{\gamma_1\wedge\gamma_2}}{(\frac{2}{\gamma_1\wedge\gamma_2+1})(2 + \frac{1.4}{2\gamma_1 c_\phi\wedge a_2})}, \] and \[ a \ge \sqrt{\frac{2(1+\xi)^{1+2/\alpha}\tau^2[\log(npT)]^{1+2/\alpha}}{T}}, \] for some $\xi>0$, then $\Pr(\mathcal{B}^c(a)) \le \pi_2(a)$, where \[ \begin{split} \pi_2(a) & := \frac{2}{(np)^\xi T^{1+\xi}} + \frac{8}{(np)^\xi T^\xi}\\ & +\frac{n^2}{a}\left(b_1e^{-c_\phi\wedge a_2(T/2)^{\gamma_2\wedge\gamma_1}} + b_8e^{-2\gamma_1c_\phi(T/2)^{\gamma_1}}\right). \end{split} \]

This bound controls the proximity between the empirical and population covariance matrices. Similar concentration inequalities were derived by aKlC2015, pLmW2012 and mcMeM2016ER. Their results, however, cannot be applied in our setting. Explicit expressions for the constants $b_1$, $b_8$ and $\tau$ in Proposition (ref) are found in Lemma (ref). Also, one may replace $\epsilon$ by its lower bound to remove dependence.

This concentration guided the choice of dependence condition used in this work. Traditionally one uses either a Hanson-Wright inequality or a Bernstein or Hoeffding type inequality to bound the empirical covariance around its mean. We write the centered Gram matrix $\boldsymbol\Gamma_T-\boldsymbol\Gamma$ as a sum of martingales and a dependence term. The martingales are handled using a Bernstein type bound and the dependence term is handled using both assumptions (A1) and (A2). Combined, they imply a sub-Weibull type decay on expected value of dependence term.

Finally, we use the bounds $\pi_1(\cdot)$ and $\pi_2(\cdot)$ in Proposition (ref) and Proposition (ref) respectively, to derive an upper bound for the prediction error and for the difference between the lasso parameter estimates and the true parameters in the $\ell_2$ norm.

theoremSuppose assumptions (A1) -- (A3) hold. Under conditions of Proposition (ref), there exists $T_0>0$ such that for all $T\ge T_0$, \[ |\boldsymbol X(\widehat{\boldsymbol\beta}_i - {\boldsymbol\beta_i}^*)/T|_2^2 \le 12|{\boldsymbol\beta_i}^*|_1\lambda , \quad i=1,...,n, \] with probability at least $1-\pi_1(\epsilon)$ for $\epsilon>0$. Suppose further that assumptions (A4) and (A5) hold. Set $\eta = \lambda/\sigma^2_\Gamma$. Under conditions of Propositions (ref) and (ref), there exists $T_0>0$ such that for all $T\ge T_0$, \[ |\widehat{\boldsymbol\beta}_i - {\boldsymbol\beta_i}^*|_2^2 \le (44+2\lambda)R_q\left(\frac{\lambda}{\sigma_\Gamma^2}\right)^{2-q},\quad i=1,...,n, \] with probability at least $1-\pi_1(\epsilon) -\pi_2\left(\frac{\sigma^{2(1-q)}_\Gamma\lambda^q}{64R_q}\right)$.

Theorem (ref) states that, with high probability, estimated and population parameter vectors are close to each other in the Euclidean norm. It requires that Propositions (ref) and (ref) hold jointly, meaning that $\lambda$, $R_q$ and $\sigma_\Gamma^2$ must satisfy rate conditions. We show that if the size of 'small' coefficients and smallest eigenvalue $\sigma_\Gamma^2$ of $\Gamma$ are restricted, then the rate of $\lambda$ in Proposition (ref) is unaffected. For $\epsilon = \log(np)$ and $T<np^2$ Proposition (ref) requires, after simplification, \[\lambda\ge \tau^*\log(np)^{2/\alpha}\sqrt{\frac{4\log(np)}{T}},\] for some constant $\tau^*$. Replacing $a$ by $\frac{\sigma^{2(1-q)}_\Gamma\lambda^q}{64R_q}$ in Proposition (ref) we obtain \[ \lambda^q \gtrsim \left(\log(np)^{2/\alpha}\sqrt{\frac{\log(np)}{T}}\right)\times \left(\frac{R_q}{\sigma_\Gamma^{2(1-q)}\log(np)^{1/\alpha}}\right). \] However, it is not necessarily a constraint in the rate of $\lambda$. Propositions (ref) and (ref) will hold jointly for $T$ sufficiently large for $0\le q<1$ if \[ \frac{R_q}{\sigma_\Gamma^{2(1-q)}} = o\left(\log(np)^{(2q-1)/\alpha}\left({\frac{T}{\log(np)}}\right)^{(1-q)/2}\right). \] In other words, if the small parameters are not too large and smallest eigenvalue of $\boldsymbol\Sigma$ is not too small as a function of $T$.

Simulation

In order to evaluate the finite-sample performance of the LASSO estimator in large VARs with stochastic volatility, we consider the following model:

equation[equation omitted — 473 chars of source]

We consider 1,000 replications of model (ref) with $T=100,300$ observations and the number of series is a function of the sample size, i.e., $n=c\times T$, where $c=\{1,2,3\}$. $\boldsymbol{A}_1$ and $\boldsymbol{A}_4$ are two block-diagonal matrices with blocks of dimension $5\times 5$ and entries equal to 0.15 and -0.1 respectively. $\boldsymbol{C}_0$ and $\boldsymbol\Psi$ are two diagonal matrices with diagonal elements all equal to 1e-5 and 0.8, respectively.

This model satisfies Assumptions (A1) -- (A5). First, Assumption (A1) is satisfied because the model is block diagonal, with fixed block size. Assumptions (A2) and (A3) follows as in Example (ref) given diagonal specification of $\boldsymbol C_0$ and $\boldsymbol\Psi$. Assumption (A4) is satisfied with $q=0$ and $R=10$. Finally, assumption (A5) holds because lower bound in equation (ref), is strictly positive ($\sum_{i=1}^4({\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol A_i \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1+{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \boldsymbol A_i \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty)<\infty$ and $\Lambda_{\min}(\boldsymbol\Sigma)>0$, independently of $n$ and $p$).

Table (ref) reports the simulation results of the LASSO estimation of model (ref) for different combinations of sample size and number of variables. The penalty parameter of the LASSO is selected by the BIC. Panel (a) reports the mean squared error of the estimation of the VAR parameters. Panel (b) reports the ratio of the mean squared out-of-sample forecasting error of the LASSO with respect to the Oracle forecast.

table[table omitted — 1,120 chars of source]

Discussion

This work provides finite sample $\ell_2$ error bounds for the equation-wise LASSO parameters estimates of a weakly sparse, high-dimensional, VAR($p$) model, with dependent and heavy tailed innovation process. It covers a large collection of specifications as illustrated in section (ref).

A distinctive feature of this work is that the dependence structure of the innovations are characterized by a very weak projective dependence condition that is naturally verifiable in settings where one is interested in the conditional variance of the process. The series of innovations is not necessarily mixing nor the resulting time series $\{{\bf y}_t\}$.

Our bounds hold under a heavy tailed setting in a sense that we do not require the moment generating function to exist. Despite the tails in $\{{\bf y}_t\}$ being sub-Weibull as in kWaTzL2017, we are not able to recover the same rates and lower bound for the regularization parameter $\lambda$. The reason is that kWaTzL2017 bounds rely heavily on the concentration inequality for mixing sequences in fMmPeR2011. Given the weak projective dependence adopted, we chose to use a martingale concentration and overcome all together the issue of using the dependence metric for deriving the concentration bound. Nevertheless, we believe the loss in efficiency is minimal. Close inspection of proof of Lemma (ref) shows that the loss of efficiency is concentrated in bounding the tail. It amounts to an extra $\log(T)$ term, which does not change the rates under assumption that $T < n^2p$, eventually.

A limitation of this work is the restriction that the model is correctly specified in the mean, in a sense that innovations are martingale differences. Nevertheless, this assumption is standard in the literature and we are able to derive results covering a broad range of data generating processes and conditional dependence measures. The martingale difference condition cannot be relaxed at this moment as our deviation bound depends on it. Furthermore, we do not require strong sparsity in a sense that near zero coefficients are effectively treated as zero as long as they are concentrated in some slowly increasing $\ell_q$ ball ($0\leq q<1$) around the origin.

Results in this paper can be extended to polynomial tails. The strategy is to replace the martingale concentration in Lemma (ref), used to prove Propositions (ref) and (ref) by \[\Pr\left(\max_{1\le i\le n}|\sum_{t=1}^T\xi_{it}|>Ta\right)\le \frac{nK}{(a\sqrt{T})^ {d}},\] whenever $\|\xi_{it}\|_d<\infty$. If available under our dependence conditions, one could employ a Fuk-Nagaev type inequality. Nevertheless, it follows that under appropriate changes to concentration rates, equation-wise LASSO estimators also admit oracle bounds. A direct consequence is that moments conditions on mCxC2002 and cHaP2009factor,cHaP2009garch are directly applicable.

Despite working with a relatively simple structure and estimation model, the machinery can be applied to more complex settings. The key points are showing that the empirical covariance concentrates around its mean in terms of its maximum entry-wise norm and the concentration inequality for large dimensional, sub-Weibull martingales. Following development of sNpRmWbY2012, the results may be naturally extended to structured regularization with node-wise regression and replacing using the Frobenius norm for system estimation. Finally, the high-dimensional VAR specification encompasses large dimensional vector-panels among other models.

In rAsSiW2020, authors consider a near epoch dependent time series. This condition covers misspecification and non-Gaussian, conditionally heteroskedastic models, such as GARCH innovations. The main focus of the authors is on inference using the desparsified lasso, but estimation error bounds are also derived. Assumptions are in the same line of mcMeM2015, where concentration bounds are assumed to hold in probability. In contrast, we focus on finite sample error bounds with relatable dependence condition on the conditional covariance. The heavy lifting in our paper is to derive the required concentration inequalities. Effectively, one could use our results to provide theoretical justification for the concentration bounds under particular model specifications. On the other hand, rAsSiW2020 provides theoretical justification for desparsified inference in some large dimensional models discussed in our paper.