# Rate-Optimal Robust Estimation of High-Dimensional Vector Autoregressive
Models

Di Wang and Ruey S. Tsay
Booth School of Business, University of Chicago

###### Abstract

High-dimensional time series data appear in many scientific areas in the
current data-rich environment. Analysis of such data poses new challenges to
data analysts because of not only the complicated dynamic dependence between
the series, but also the existence of aberrant observations, such as missing
values, contaminated observations, and heavy-tailed distributions. For high-
dimensional vector autoregressive (VAR) models, we introduce a unified
estimation procedure that is robust to model misspecification, heavy-tailed
noise contamination, and conditional heteroscedasticity. The proposed
methodology enjoys both statistical optimality and computational efficiency,
and can handle many popular high-dimensional models, such as sparse, reduced-
rank, banded, and network-structured VAR models. With proper regularization
and data truncation, the estimation convergence rates are shown to be almost
optimal in the minimax sense under a bounded $(2+2\epsilon)$-th moment
condition. When $\epsilon\geq 1$, the rates of convergence match those
obtained under the sub-Gaussian assumption. Consistency of the proposed
estimators is also established for some $\epsilon\in(0,1)$, with minimax
optimal convergence rates associated with $\epsilon$. The efficacy of the
proposed estimation methods is demonstrated by simulation and a U.S.
macroeconomic example.

Keywords: Autocovariance, high-dimensional time series, minimax optimal,
robust statistics, truncation

##  1 Introduction

###  1.1 High-dimensional vector autoregression

Vector autoregressive (VAR) models are arguably the most commonly used
multivariate time series models in practice; see, e.g., Lütkepohl, (2005),
Tsay, (2013), and the references therein. Applications of the model can be
found in a wide range of fields, such as economics and finance (Wu and Xia,,
2016), time-course functional genomics (Michailidis and d’Alché Buc,, 2013),
and neuroimaging (Gorrostieta et al.,, 2012). Consider a $p$-dimensional zero-
mean VAR model of order $d$, i.e., VAR($d$) model,

$\mathbf{y}_{t}=\mathbf{A}_{1}\mathbf{y}_{t-1}+\mathbf{A}_{2}\mathbf{y}_{t-2}+\cdots+\mathbf{A}_{d}\mathbf{y}_{t-d}+\boldsymbol{\varepsilon}_{t},~{}~{}~{}~{}t=1,2,\dots,T,$
(1)

where $\mathbf{y}_{t}=(y_{1t},\dots,y_{pt})^{\top}\in\mathbb{R}^{p}$ is the
observed time series, $\mathbf{A}_{j}\in\mathbb{R}^{p\times p}$ is the lag-$j$
coefficient matrix, and
$\boldsymbol{\varepsilon}_{t}=(\varepsilon_{1t},\dots,\varepsilon_{pt})^{\top}\in\mathbb{R}^{p}$
is a serially uncorrelated white noise innovation. We assume that all
solutions of the determinant equation $|\mathbf{A}(z)|=0$ are outside the unit
circle, where
$\mathbf{A}(z)=\mathbf{I}_{p}-\mathbf{A}_{1}z-\cdots-\mathbf{A}_{d}z^{d}$ is
referred to as the AR matrix polynomial in $z$. In modern applications, the
dimension $p$ is often large. However, since the number of coefficient
parameters is $p^{2}d$ and those coefficients are often highly correlated, an
unrestricted VAR($d$) model is likely to encounter the difficulty of over-
parameterization and, hence, cannot provide reliable estimates nor accurate
forecasts without further restrictions.

Estimation consistency of high-dimensional VAR models is achievable under some
structural assumptions, provided that certain regularity conditions are
satisfied. For example, if the coefficient matrices have an unobserved low-
dimensional structure, such as sparsity or low-rankness, the structure-
inducing regularization methods, including Lasso (Basu and Michailidis,,
2015), Dantzig selector (Han et al.,, 2015), and nuclear norm penalty
(Negahban and Wainwright,, 2011), give consistent estimates under the Gaussian
assumption of the time series. Recently, Zheng and Raskutti, (2019), Zheng and
Cheng, (2021), and Wang et al., 2021b  developed novel technical tools to
relax the distributional assumption from Gaussian to sub-Gaussian.

However, in real applications, time series data often contain aberrant
observations, which can occur in many ways, such as missing values,
measurement error contamination, and heavy-tailed distribution. Those aberrant
observations, if overlooked, can lead to biased estimates, erroneous
inference, and sub-optimal forecasts. The situation can easily be further
exacerbated when the dimension $p$ is large. It is, therefore, important to
study robust estimation of high-dimensional VAR models. In addition, the true
data generating process is unlikely to follow a VAR model. Model uncertainty,
including using a high-dimensional VAR($d$) model as an approximation to the
true model, also deserves a careful investigation.

Recently there have been emerging interests in studying high-dimensional VAR
models with non-i.i.d. and/or non-sub-Gaussian innovations. For example, Wu
and Wu, (2016) studied theoretical properties of Lasso and constrained
$\ell_{1}$ minimization estimators for a VAR model with weakly correlated and
heavy-tailed $\boldsymbol{\varepsilon}_{t}$. Wong et al., (2020) investigated
the estimation and prediction performance of Lasso for the sub-Weibull time
series data under a $\beta$-mixing condition. Both theoretical and numerical
results in the literature show that the performance of standard $\ell_{1}$
regularized estimators deteriorates substantially when the data have heavy
tails. For robust estimation of high-dimensional heavy-tailed time series
data, Qiu et al., (2015) developed a quantile-based Dantzig selector for the
class of elliptical VAR processes. Han et al., (2020) proposed a robust
estimation method for high-dimensional sparse generalized linear models with
temporal dependent covariates. However, the existing literature on robust
estimation for time series data focuses on sparse models. To the best of our
knowledge, there is no unified solution to address the robust estimation
problem for a large class of high-dimensional VAR models.

Our proposed robust estimation procedure is built on two key ingredients: the
constrained Yule–Walker estimator and the robust autocovariance matrix
estimator. The first ingredient, the constrained Yule–Walker estimator,
provides a general and flexible estimation framework for two classes of high-
dimensional models, namely the approximately low-dimensional VAR models and
linear-restricted VAR models. The second ingredient, the robust autocovariance
matrix estimator, is easy to implement by truncating the time series data.
However, based on the specific model structure, we need to adapt the data
truncation methods. With a large class of distributions having bounded second
or higher-order moments, the proposed estimators are shown to be consistent
under high-dimensional scaling. To be specific, we summarize the main
contributions of this paper as follows:

  * (i)

For various high-dimensional VAR models, the paper provides a simple and
general estimation procedure robust to model misspecification, heavy-tailed
noise contamination, and conditional heteroskedasticity. Our proposal can
handle many popular high-dimensional VAR models, such as sparse, reduced-rank,
banded, and network VAR models. An efficient and scalable alternating
direction method of multipliers (ADMM) algorithm is developed.

  * (ii)

The proposed methodology enjoys statistical optimality in the minimax sense.
Our theoretical framework deals with heavy-tailed distributions with bounded
$(2+2\epsilon)$-th moments for any $\epsilon>0$. It results in a phase
transition on the rates of convergence: for $\epsilon\geq 1$, the estimator
achieves the same convergence rates as those obtained under the Gaussian or
sub-Gaussian distribution, while consistency is also established with a slower
rate for $\epsilon\in(0,1)$. By establishing the matching minimax lower
bounds, we show that the proposed estimators for high-dimensional VAR models
and autocovariance matrices are rate-optimal.

###  1.2 Related literature

This work is related to a huge body of literature on the robust estimation of
high-dimensional regression and covariance matrices. The early developments in
robust statistics were pioneered by Huber, (1964) and Hampel, (1971, 1974);
see also the overview by Hampel, (2001) and the references therein. For time
series data, robust estimation methods were developed for univariate and
multivariate ARMA models (Martin,, 1981; Muler et al.,, 2009; Muler,, 2013).
For high-dimensional i.i.d. data, inspired by Catoni, (2012), the non-
asymptotic deviation analysis for heavy-tailed variables and robust
$M$-estimators for high-dimensional regression were proposed by Fan et al.,
(2017), Loh, (2017), Sun et al., (2020), Wang et al., (2020), and many others.
These recent robust estimation methods and their theoretical results were
further extended to high-dimensional covariance and precision matrix
estimation problems, such as Avella-Medina et al., (2018), Minsker, (2018),
and Zhang, (2021). Another research line of robust regression is the least
absolute deviation loss, and more generally, quantile loss, which have been
studied extensively; see, e.g., Wang et al., (2007), Belloni and Chernozhukov,
(2011), and Wang et al., (2012). Recently, Fan et al., (2021) and Ke et al.,
(2019) independently proposed truncated variants of the sample covariance,
which are easy to implement and motivate us to develop the robust estimators
of high-dimensional autocovariance matrices.

The upper and lower bound development for high-dimensional estimation problems
under heavy-tailed distributions is another emerging and important research
topic. Under heavy-tailed distributions, the phase transition phenomenon in
the rate of convergence was previously discovered by Bubeck et al., (2013),
Avella-Medina et al., (2018), Sun et al., (2020), Tan et al., (2022), and
others. For the lower bound development, the phase transition phenomenon was
first established by Devroye et al., (2016) for univariate mean estimation
problem, and was extended to the fixed and high-dimensional linear regression
problems in Sun et al., (2020). Compared with the existing literature, our
minimax lower bound results are established under a much more complicated
setting. First, all existing lower bound results are developed for i.i.d.
data, but we allow weak serial dependency and establish the minimax lower
bounds under strong mixing conditions. Second, for the lower bounds of the
linear regression problem in Sun et al., (2020), a finite moment condition is
imposed on the random error terms while the covariates are assumed to be
light-tailed. Without any assumption on data generating process, we consider
the finite moment condition on the observed time series data that serve as
both predictors and responses in the autoregressive models, which is
fundamentally different from the conventional linear regression problem.

###  1.3 Notation and outline

We start with some notations used in the paper. Let $C$ denote a generic
positive constant, which is independent of the dimension and sample size. For
any two real-valued sequences $x_{k}$ and $y_{k}$, $x_{k}\gtrsim y_{k}$ if
there exists a $C>0$ such that $x_{k}\geq Cy_{k}$ for all $k$. In addition, we
write $x_{k}\asymp y_{k}$ if $x_{k}\gtrsim y_{k}$ and $y_{k}\gtrsim x_{k}$.
For any two real numbers $x$ and $y$, let $x\wedge y$ denote their minimum.
Throughout the paper, we use bold lowercase letters to denote vectors. For any
vector $\mathbf{v}=(v_{1},v_{2},\dots,v_{d})^{\top}$, denote its $\ell_{q}$
norm as $\|\mathbf{v}\|_{q}=(\sum_{i=1}^{d}v_{i}^{q})^{1/q}$ for $1\leq
q<\infty$, its maximum norm as $\|\mathbf{v}\|_{\infty}=\max_{1\leq i\leq
d}|v_{i}|$, and its $\ell_{0}$ norm as
$\|\mathbf{v}\|_{0}=\sum_{i=1}^{d}1\\{v_{i}\neq 0\\}$. We use bold uppercase
letters to denote matrices. For any matrix
$\mathbf{M}=[\mathbf{m}_{1},\dots,\mathbf{m}_{d_{2}}]\in\mathbb{R}^{d_{1}\times
d_{2}}$, where $\mathbf{m}_{j}\in\mathbb{R}^{d_{1}}$ is the $j$-th column of
$\mathbf{M}$, denote its $\ell_{p,q}$ norm as
$\|\mathbf{M}\|_{p,q}=(\sum_{j=1}^{d_{2}}\|\mathbf{m}_{j}\|_{p}^{q})^{1/q}$,
for any $1\leq p,q\leq\infty$, and its $i$th largest singular value as
$\sigma_{i}(\mathbf{M})$, for all $i=1,2,\dots,d_{1}\wedge d_{2}$. For a given
matrix $\mathbf{M}\in\mathbb{R}^{d_{1}\times d_{2}}$, we let
$\|\mathbf{M}\|_{\textup{F}}$, $\|\mathbf{M}\|_{\textup{op}}$, and
$\|\mathbf{M}\|_{\textup{nuc}}$ denote its Frobenius norm, operator norm, and
nuclear norm, respectively, where
$\|\mathbf{M}\|_{\textup{F}}=\|\mathbf{M}\|_{2,2}$,
$\|\mathbf{M}\|_{\textup{op}}=\sigma_{1}(\mathbf{M})$, and
$\|\mathbf{M}\|_{\textup{nuc}}=\sum_{i=1}^{d_{1}\wedge
d_{2}}\sigma_{i}(\mathbf{M})$. For any two matrices $\mathbf{M}_{1}$ and
$\mathbf{M}_{2}$, their Kronecker product is
$\mathbf{M}_{1}\otimes\mathbf{M}_{2}$. For any subspace
$\mathcal{M}\subset\mathbb{R}^{d_{1}\times d_{2}}$ and any matrix
$\mathbf{M}\in\mathbb{R}^{d_{1}\times d_{2}}$, denote by $\mathcal{M}^{\perp}$
and $\mathbf{M}_{\mathcal{M}}$ the orthogonal complement of $\mathcal{M}$ and
the projection of $\mathbf{M}$ onto $\mathcal{M}$, respectively.

The rest of the paper is organized as follows. Section 2 studies constrained
Yule–Walker estimators for two classes of VAR models: VAR with approximately
low-dimensional structure and VAR with linear restrictions. Section 3 develops
robust autocovariance estimators tailored for various low-dimensional
structures. Minimax lower bounds for both high-dimensional VAR estimation and
autocovariance matrix estimation are investigated in Section 4. Computational
algorithms and implementation details are discussed in Section 5. Section 6
presents some simulation results, and Section 7 shows an empirical application
of U.S. macroeconomic data. All technical proofs and detailed algorithms are
relegated to the Appendices. The codes and data can be found at
https://github.com/diwangstat/RobustAR.

##  2 Constrained Yule–Walker Estimation

###  2.1 Yule–Walker equation

Consider a general mean-zero and covariance stationary process
$\\{\mathbf{y}_{t}\\}_{t=1}^{T}$, where $\mathbf{y}_{t}\in\mathbb{R}^{p}$. In
many applications, it is common to model the time series data and predict
their future values using a linear VAR model. The VAR($d$) model in (1) can be
rewritten as

$\mathbf{y}_{t}=\mathbf{A}\mathbf{x}_{t}+\boldsymbol{\varepsilon}_{t},$ (2)

where
$\mathbf{x}_{t}=(\mathbf{y}_{t-1}^{\top},\dots,\mathbf{y}_{t-d}^{\top})^{\top}\in\mathbb{R}^{pd}$
is the predictor vector and
$\mathbf{A}=[\mathbf{A}_{1},\mathbf{A}_{2},\cdots,\mathbf{A}_{d}]\in\mathbb{R}^{p\times
pd}$ is the combined parameter matrix. Throughout this paper, we consider that
the lag order $d$ is fixed. Given the stationarity of $\mathbf{y}_{t}$, the
parameter matrix of interest $\mathbf{A}^{*}\in\mathbb{R}^{p\times pd}$ is
defined as the minimizer of the risk function

$\mathbf{A}^{*}:=\underset{\mathbf{A}\in\mathbb{R}^{p\times
pd}}{\arg\min}~{}\mathbb{E}\left[\|\mathbf{y}_{t}-\mathbf{A}\mathbf{x}_{t}\|_{2}^{2}\right].$
(3)

For any integer $\ell$, denote the lag-$\ell$ autocovariance matrix of
$\mathbf{y}_{t}$ by
$\mathbf{\Gamma}_{\ell}=\mathbb{E}[\mathbf{y}_{t}\mathbf{y}_{t-\ell}^{\top}]$.
By simple algebra,

$\begin{split}\|\mathbf{y}_{t}-\mathbf{A}\mathbf{x}_{t}\|_{2}^{2}=&\|\mathbf{y}_{t}-\textup{vec}(\mathbf{x}_{t}^{\top}\mathbf{A}^{\top})\|_{2}^{2}=\|\mathbf{y}_{t}-(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\textup{vec}(\mathbf{A}^{\top})\|_{2}^{2}\\\
=&~{}\mathbf{y}_{t}^{\top}\mathbf{y}_{t}-2\mathbf{y}_{t}^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\textup{vec}(\mathbf{A}^{\top})+\textup{vec}(\mathbf{A}^{\top})^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t}\mathbf{x}_{t}^{\top})\textup{vec}(\mathbf{A}^{\top})\\\
=&~{}\mathbf{y}_{t}^{\top}\mathbf{y}_{t}-2\textup{vec}(\mathbf{x}_{t}\mathbf{y}_{t}^{\top})^{\top}\textup{vec}(\mathbf{A}^{\top})+\textup{vec}(\mathbf{A}^{\top})^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t}\mathbf{x}_{t}^{\top})\textup{vec}(\mathbf{A}^{\top}),\end{split}$
(4)

which implies that $\mathbf{\Sigma}_{1}=\mathbf{A}^{*}\mathbf{\Sigma}_{0}$ or
$\mathbf{A}^{*}=\mathbf{\Sigma}_{1}\mathbf{\Sigma}_{0}^{-1}$ if
$\mathbf{\Sigma}_{0}$ is invertible, where
$\mathbf{\Sigma}_{1}:=\mathbb{E}[\mathbf{y}_{t}\mathbf{x}_{t}^{\top}]=[\mathbf{\Gamma}_{1},\mathbf{\Gamma}_{2},\dots,\mathbf{\Gamma}_{d}]$
and

$\mathbf{\Sigma}_{0}:=\mathbb{E}[\mathbf{x}_{t}\mathbf{x}_{t}^{\top}]=\begin{bmatrix}\mathbf{\Gamma}_{0}&\mathbf{\Gamma}_{1}&\cdots&\mathbf{\Gamma}_{d-1}\\\
\mathbf{\Gamma}_{1}^{\top}&\mathbf{\Gamma}_{0}&\cdots&\mathbf{\Gamma}_{d-2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{\Gamma}_{d-1}^{\top}&\mathbf{\Gamma}_{d-2}^{\top}&\cdots&\mathbf{\Gamma}_{0}\end{bmatrix}.$
(5)

This relationship between the parameter matrix $\mathbf{A}^{*}$ and
autocovariance matrices is well known as the multivariate Yule–Walker equation
for VAR($d$) models (Lütkepohl,, 2005; Tsay,, 2013). When the data generating
mechanism of $\mathbf{y}_{t}$ is not the VAR($d$) process, the Yule–Walker
equation still holds if VAR is used as a running model.

The existence of aberrant observations, such as missing values or data
contaminated by measurement errors, is ubiquitous in high-dimensional
macroeconomic, environmental, and genetic data, among many other applications.
In a VAR process, when some of the variables are removed as they cannot be
observed or measured, the rest of the available variables generally would not
follow any finite-order VAR model (Lütkepohl,, 2005). In addition, if the true
signal of the time series follows a VAR process and the observed time series
$\mathbf{y}_{t}$ is a contaminated version of the true signal plus a
measurement error, $\mathbf{y}_{t}$ no longer follows a VAR process. Hence,
the VAR model assumption commonly used in real applications is questionable,
but the Yule–Walker equation which holds without any data generating mechanism
assumption can be utilized to construct estimation methodology for the general
covariance stationary time series.

In addition to the sub-Gaussian innovation assumption, another key limitation
in the existing literature is the i.i.d. assumption for the innovation series
$\boldsymbol{\varepsilon}_{t}$. The conditional heteroskedasticity is often
observed in financial time series, and violation of the homogeneous error
assumption can affect the performance of standard estimation methods. Another
major advantage of our moment-based methodology is that the i.i.d. assumption
for $\boldsymbol{\varepsilon}_{t}$ can be relaxed to a serially uncorrelated,
but weakly stationary condition. Specifically, our theoretical analysis allows
some serial dependence in $\boldsymbol{\varepsilon}_{t}$, dependence between
$\boldsymbol{\varepsilon}_{t}$ and $\mathbf{x}_{t}$, and conditional
heteroskedasticity in $\boldsymbol{\varepsilon}_{t}$. Note that this setting
of innovations includes many important conditional heteroskedasticity models,
such as GARCH models (Engle,, 1982; Bollerslev,, 1986) and double AR models
(Ling,, 2004; Zhu et al.,, 2018).

Based on the Yule–Walker equation, we propose a constrained minimization
estimation approach, named the constrained Yule–Walker estimation, for two
classes of high-dimensional VAR($d$) models. In Section 2.2, we consider the
VAR($d$) model, where $\mathbf{A}^{*}$ in (3) can be approximated by a matrix
with some latent low-dimensional structure. In Section 2.3, we consider the
VAR($d$) model whose coefficients are subject to some linear restrictions.
These two classes of models are shown to encompass many commonly-used high-
dimensional VAR models in practice.

###  2.2 Approximately low-dimensional VAR

We first consider the situations where the dimension $p$ is relatively large
compared to the sample size $T$ and the parameter matrix of interest
$\mathbf{A}^{*}$ in (3) can be well approximated by a matrix with certain
types of low-dimensional structure. The approximately low-dimensional
structure, such as weak sparsity and approximate low-rankness, is general and
natural in high-dimensional time series modeling.

Based on the Yule–Walker equation, we consider a general class of constrained
estimation:

$\mathbf{\widehat{A}}=\underset{\mathbf{A}\in\mathbb{R}^{p\times
pd}}{\arg\min}~{}\mathcal{R}(\mathbf{A})~{}~{}~{}\text{such
that}~{}~{}\mathcal{R}^{*}(\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{A}\mathbf{\widetilde{\Sigma}}_{0})\leq\lambda,$
(6)

where $\mathcal{R}(\cdot)$ is a matrix norm as the regularization function,
$\mathcal{R}^{*}(\cdot)$ is its dual norm as the constraint function,
$\lambda$ is the constraint parameter, and $\mathbf{\widetilde{\Sigma}}_{1}$
and $\mathbf{\widetilde{\Sigma}}_{0}$ are robust autocovariance estimators to
be specified later in Section 3. When $\lambda$ is sufficiently small, the
constraint function $\mathcal{R}^{*}(\cdot)$ guarantees that the sample
version of the Yule–Walker equation holds roughly, and the regularizer
$\mathcal{R}(\cdot)$ induces the low-dimensional structure and improves the
estimation efficiency.

Following the framework for high-dimensional analysis in Negahban et al.,
(2012), we consider a decomposable regularizer $\mathcal{R}(\cdot)$. For a
generic low-dimensional structure, define
$\mathcal{M}\subseteq\overline{\mathcal{M}}\subset\mathbb{R}^{p\times pd}$,
where $\mathcal{M}$ is referred to as the model subspace to represent the
specific model constraints; for instance, it can be the subspace of low-rank
matrices (see Example 1). The orthogonal complement of
$\overline{\mathcal{M}}$, denoted by $\overline{\mathcal{M}}^{\perp}$, is the
associated perturbation subspace to capture the deviation from the model
subspace and is adopted to measure the error of approximation to the low-
dimensional structure. We assume that the true value $\mathbf{A}^{*}$ can be
decomposed into its projections onto $\mathcal{M}$ and
$\overline{\mathcal{M}}^{\perp}$, i.e.,
$\mathbf{A}^{*}=\mathbf{A}^{*}_{\mathcal{M}}+\mathbf{A}^{*}_{\overline{\mathcal{M}}^{\perp}}$.

We call a regularizer $\mathcal{R}(\cdot)$ decomposable with respect to a pair
of subspaces $(\mathcal{M},\overline{\mathcal{M}}^{\perp})$, if for any
$\mathbf{W}_{1}\in\mathcal{M}$ and
$\mathbf{W}_{2}\in\overline{\mathcal{M}}^{\perp}$,

$\mathcal{R}(\mathbf{W}_{1}+\mathbf{W}_{2})=\mathcal{R}(\mathbf{W}_{1})+\mathcal{R}(\mathbf{W}_{2}).$
(7)

Many of the commonly-used convex regularizers, such as the nuclear norm for
low-rankness, are shown to be decomposable; see Negahban et al., (2012) for
more details of decomposable regularizers. To measure the magnitude of the
low-dimensional structure, define by $\phi(\overline{\mathcal{M}})$ the
associated constant such that
$\mathcal{R}(\mathbf{W})\leq\phi(\overline{\mathcal{M}})\mathcal{R}^{*}(\mathbf{W})$
for any $\mathbf{W}\in\overline{\mathcal{M}}$. In our theoretical analysis, we
consider another auxiliary matrix norm $\mathcal{C}(\cdot)$ such that
$(\mathcal{R}(\cdot),\mathcal{R}^{*}(\cdot),\mathcal{C}(\cdot))$ satisfy that
$\mathcal{R}^{*}(\mathbf{W}_{1}\mathbf{W}_{2})\leq\mathcal{R}^{*}(\mathbf{W}_{1})\mathcal{C}(\mathbf{W}_{2})$
and
$\mathcal{R}^{*}(\mathbf{W}_{1}\mathbf{W}_{2})\leq\mathcal{R}(\mathbf{W}_{1})\mathcal{R}^{*}(\mathbf{W}_{2})$
for any symmetric matrix $\mathbf{W}_{2}$ and any compatible matrix
$\mathbf{W}_{1}$. The following example illustrates the suitable choice of
$\mathcal{R}(\cdot)$, $\mathcal{R}^{*}(\cdot)$, and
$(\mathcal{M},\overline{\mathcal{M}}^{\perp})$ for a reduced-rank VAR model.

######  Example 1 (Reduced-rank VAR).

The reduced-rank VAR model is an important approach to modeling high-
dimensional time series by imposing a low-rank structure on $\mathbf{A}$; see
also Velu and Reinsel, (2013), Basu et al., (2019), and Wang et al., 2021a .
To induce low-rankness, we use the nuclear norm as the regularizer and the
operator norm as the constraint function; that is,
$\mathcal{R}(\cdot)=\|\cdot\|_{\textup{nuc}}$ and
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\textup{op}}$. By the submultiplicative
property of the operator norm, $\mathcal{C}(\cdot)=\|\cdot\|_{\textup{op}}$.
For the approximately low-rank matrix $\mathbf{A}^{*}$, denote by
$\mathcal{U}$ and $\mathcal{V}$ the subspace spanned by its leading left and
right singular vectors, respectively. Define the model subspace

$\mathcal{M}(\mathcal{U},\mathcal{V})=\\{\mathbf{W}\in\mathbb{R}^{p\times
pd}:\text{col}(\mathbf{W})\subseteq\mathcal{U},~{}\text{col}(\mathbf{W}^{\top})\subseteq\mathcal{V}\\},$
(8)

and the perturbation space

$\overline{\mathcal{M}}^{\perp}(\mathcal{U},\mathcal{V})=\\{\mathbf{W}\in\mathbb{R}^{p\times
pd}:\text{col}(\mathbf{W})\perp\mathcal{U},~{}\text{col}(\mathbf{W}^{\top})\perp\mathcal{V}\\}.$
(9)

For the general mean-zero and covariance stationary process
$\\{\mathbf{y}_{t}\\}_{t=1}^{T}$, if we model the data by a VAR($d$) model in
(1) and apply the constrained Yule–Walker estimator in (6) with robust
autocovariance estimators $\mathbf{\widetilde{\Sigma}}_{0}$ and
$\mathbf{\widetilde{\Sigma}}_{1}$, we have the following upper bounds for the
estimation error between the estimated matrix $\mathbf{\widehat{A}}$ and the
true value $\mathbf{A}^{*}$ defined in (3).

######  Proposition 1.

For the constrained Yule–Walker estimator in (6), suppose that
$\mathcal{R}^{*}(\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\Sigma}_{0})\leq\zeta_{0}$,
$\mathcal{R}^{*}(\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{\Sigma}_{1})\leq\zeta_{1}$,
and $\lambda\geq\zeta_{0}\mathcal{R}(\mathbf{A}^{*})+\zeta_{1}$. Then, for
$(\mathcal{M},\overline{\mathcal{M}}^{\perp})$,

$\begin{split}\mathcal{R}^{*}(\mathbf{\widehat{A}}-\mathbf{A}^{*})&\leq
2\hskip 0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})\lambda,\\\
\mathcal{R}(\mathbf{\widehat{A}}-\mathbf{A}^{*})&\leq 4\hskip
0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})\phi(\overline{\mathcal{M}})\lambda+2\mathcal{R}(\mathbf{A}^{*}_{\overline{\mathcal{M}}^{\perp}}),\\\
\text{and}~{}~{}\|\widehat{\mathbf{A}}-\mathbf{A}^{*}\|_{\textup{F}}^{2}&\leq
8\hskip
0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})^{2}\phi(\overline{\mathcal{M}})\lambda^{2}+4\hskip
0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})\mathcal{R}(\mathbf{A}^{*}_{\overline{\mathcal{M}}^{\perp}})\lambda.\end{split}$
(10)

This proposition presents the deterministic estimation error bounds in
$\mathcal{R}^{*}(\cdot)$, $\mathcal{R}(\cdot)$, and
$\|\cdot\|_{\textup{F}}^{2}$, given that the estimators
$\mathbf{\widetilde{\Sigma}}_{0}$ and $\mathbf{\widetilde{\Sigma}}_{1}$
satisfy some regularity conditions. The estimation error bounds are related to
$\mathcal{R}(\mathbf{A}^{*})$ and $\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})$, and
these terms could be bounded or diverge slowly as the dimension $p$ increases;
see more discussions for each specific model in Section 3. When the true
parameter matrix $\mathbf{A}^{*}$ admits a low-dimensional approximation
rather than an exact structure,
$\mathcal{R}(\mathbf{A}^{*}_{\overline{\mathcal{M}}^{\perp}})$ is strictly
positive, and the second terms in the second and third upper bounds are
related to the approximation error. In this case, the upper bounds in this
proposition hold uniformly for a class of low-dimensional structures with
respect to $(\mathcal{M},\overline{\mathcal{M}}^{\perp})$, and the class of
upper bounds can be minimized by balancing the estimation error on
$\overline{\mathcal{M}}$ and the approximation error on
$\overline{\mathcal{M}}^{\perp}$.

Moreover, in some cases, it is of interest to consider the approximate low
dimensionality in a more structured manner. In particular, based on the multi-
response nature, the VAR model in (2) can be split into $p$ sub-models:

$y_{it}=\mathbf{a}_{i}^{\top}\mathbf{x}_{t}+\varepsilon_{it}$ (11)

where $\mathbf{a}_{i}^{\top}$ is the $i$-th row of $\mathbf{A}$, for
$i=1,\dots,p$. For these sub-models, as each $\mathbf{a}_{i}$ is a
$pd$-dimensional vector, we may consider some low-dimensional structures, such
as sparsity, on $\mathbf{a}_{i}$. Then, we consider the constrained estimation
framework for each $\mathbf{a}_{i}$:

$\widehat{\mathbf{a}}_{i}=\operatorname*{arg\,min}_{\mathbf{a}\in\mathbb{R}^{pd}}\mathcal{R}(\mathbf{a})~{}~{}\text{such
that}~{}~{}\mathcal{R}^{*}(\widetilde{\boldsymbol{\sigma}}_{1i}-\widetilde{\mathbf{\Sigma}}_{0}\mathbf{a})\leq\lambda,$
(12)

where $\mathcal{R}(\cdot)$ is a vector norm, $\mathcal{R}^{*}(\cdot)$ is its
dual norm, and $\widetilde{\boldsymbol{\sigma}}_{1i}^{\top}$ is the $i$-th row
of $\widetilde{\mathbf{\Sigma}}_{1}$. For simplicity, we consider the same
constraint parameter $\lambda$ for all sub-problems. In this framework, these
$p$ sub-problems can be handled in parallel and thus can be solved
efficiently.

In an analogous fashion, we consider a pair of subspaces of $\mathbb{R}^{pd}$,
$(\mathcal{M}_{i},\overline{\mathcal{M}}_{i}^{\perp})$, for each
$\mathbf{a}_{i}$, where
$\mathbf{a}^{*}_{i}=(\mathbf{a}^{*}_{i})_{\mathcal{M}_{i}}+(\mathbf{a}^{*}_{i})_{\overline{\mathcal{M}}_{i}^{\perp}}$
and
$\mathcal{R}(\mathbf{v}_{1}+\mathbf{v}_{2})=\mathcal{R}(\mathbf{v}_{1})+\mathcal{R}(\mathbf{v}_{2})$
for any $\mathbf{v}_{1}\in\mathcal{M}_{i}$ and
$\mathbf{v}_{2}\in\overline{\mathcal{M}}_{i}^{\perp}$, for all $i=1,\dots,p$.
In addition, define by $\phi(\overline{\mathcal{M}}_{i})$ the constant such
that
$\mathcal{R}(\mathbf{v})\leq\phi(\overline{\mathcal{M}}_{i})\mathcal{R}^{*}(\mathbf{v})$
for any $\mathbf{v}\in\overline{\mathcal{M}}_{i}$. Moreover, for theoretical
analysis, we consider the auxiliary matrix norms
$\overline{\mathcal{R}}^{*}(\cdot)$ and $\mathcal{C}(\cdot)$ satisfying that
$\mathcal{R}^{*}(\mathbf{W}\mathbf{v})\leq\overline{\mathcal{R}}^{*}(\mathbf{W})\mathcal{R}(\mathbf{v})$
and
$\mathcal{R}^{*}(\mathbf{W}^{\top}\mathbf{v})\leq\mathcal{C}(\mathbf{W})\mathcal{R}^{*}(\mathbf{v})$
for any compatible matrix $\mathbf{W}$ and vector $\mathbf{v}$. The following
example shows that the $\ell_{1}$ regularization can be applied to the sparse
VAR model.

######  Example 2 (Sparse VAR).

The sparse VAR model assumes that each $\mathbf{a}_{i}$ is strictly or weakly
sparse, so that each response variable $y_{it}$ is only related to a subset of
lagged values $y_{j,t-\ell}$. Though the overall sparsity in $\mathbf{A}$ has
been considered in many existing literature, such as Basu and Michailidis,
(2015) and Kock and Callot, (2015), the sparsity structure in rows seems more
natural in the time series context. For example, in the diagonal VAR(1) model,
the total degree of sparsity grows with $p$, while the sparsity level in each
row remains fixed. To induce sparsity in each row, as in Han et al., (2015),
we consider $\mathcal{R}(\cdot)=\|\cdot\|_{1}$,
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\infty}$,
$\overline{\mathcal{R}}^{*}(\cdot)=\|\cdot\|_{\infty}$, and
$\mathcal{C}(\cdot)=\|\cdot\|_{1,\infty}$. For the weakly sparse
$\mathbf{a}_{i}$, $i=1,\dots,p$, denote by $S_{i}\subset\\{j:1\leq j\leq
pd\\}$ any index set corresponding to those coefficients significantly distant
from zero in $\mathbf{a}_{i}$. For each $i$, define the model subspace
associated with the chosen index set $S_{i}$ as

$\mathcal{M}_{i}(S_{i})=\\{\mathbf{v}\in\mathbb{R}^{pd}:\mathbf{v}_{j}=0~{}\text{for
all}~{}j\notin S_{i}\\},$ (13)

and the perturbation subspace

$\overline{\mathcal{M}}_{i}^{\perp}(S_{i})=\\{\mathbf{v}\in\mathbb{R}^{pd}:\mathbf{v}_{j}=0~{}\text{for
all}~{}j\in S_{i}\\}.$ (14)

For the general mean-zero and covariance stationary process
$\\{\mathbf{y}_{t}\\}_{t=1}^{T}$, if we separate the VAR($d$) model to $p$
sub-problems in (11) and apply the separated constrained Yule–Walker estimator
in (12) with robust autocovariance estimators
$\mathbf{\widetilde{\Sigma}}_{0}$ and $\widetilde{\boldsymbol{\sigma}}_{1i}$,
we have the following upper bounds for the estimation error between
$\mathbf{\widehat{a}}_{i}$ and the true value $\mathbf{a}_{i}^{*}$.

######  Proposition 2.

For the separated Yule–Walker estimators in (12), suppose that
$\overline{\mathcal{R}}^{*}(\widetilde{\mathbf{\Sigma}}_{0}-\mathbf{\Sigma}_{0})\leq\zeta_{0}$,
$\mathcal{R}^{*}(\widetilde{\boldsymbol{\sigma}}_{1i}-\boldsymbol{\sigma}_{1i})\leq\zeta_{1i}$,
and $\lambda\geq\max_{1\leq i\leq
p}[\zeta_{0}\mathcal{R}(\mathbf{a}_{i}^{*})+\zeta_{1i}]$. Then, for
$(\mathcal{M}_{i},\overline{\mathcal{M}}_{i}^{\perp})$,

$\begin{split}\mathcal{R}^{*}(\widehat{\mathbf{a}}_{i}-\mathbf{a}_{i}^{*})&\leq
2\hskip 0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})\lambda,\\\
\mathcal{R}(\widehat{\mathbf{a}}_{i}-\mathbf{a}_{i}^{*})&\leq 4\hskip
0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})\phi(\overline{\mathcal{M}}_{i})\lambda+2\mathcal{R}(\mathbf{A}^{*}_{\overline{\mathcal{M}}_{i}^{\perp}}),\\\
\text{and}~{}~{}\|\widehat{\mathbf{a}}_{i}-\mathbf{a}_{i}^{*}\|_{2}^{2}&\leq
8\hskip
0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})^{2}\phi(\overline{\mathcal{M}}_{i})\lambda^{2}+4\hskip
0.72229pt\mathcal{C}(\mathbf{\Sigma}_{0}^{-1})\mathcal{R}(\mathbf{A}^{*}_{\overline{\mathcal{M}}_{i}^{\perp}})\lambda\end{split}$
(15)

for $i=1,\dots,p$.

For the approximately low-rank and sparse VAR models discussed in Examples 1
and 2, we can obtain a sharp upper bound by appropriately selecting the matrix
rank or the degree of sparsity. More concrete results are given in Section 3.

###  2.3 Linear-restricted VAR

In modern high-dimensional time series modeling, domain knowledge and prior
information are commonly available in addition to the time series data. In
many such cases, prior information can be formulated as linear restrictions on
the parameter matrix, and in this subsection, we apply the constrained
Yule–Walker estimation approach to the high-dimensional VAR models with linear
restrictions.

As discussed in Tsay, (2013), a general linear parameter restriction can be
expressed as
$\textup{vec}(\mathbf{A}^{\top})=\mathbf{C}\boldsymbol{\theta}+\boldsymbol{\gamma}$,
where $\mathbf{C}\in\mathbb{R}^{p^{2}d\times r}$ is a prespecified constraint
matrix, $\boldsymbol{\gamma}\in\mathbb{R}^{p^{2}d}$ is a known constant
vector, and $\boldsymbol{\theta}\in\mathbb{R}^{r}$ is the unknown parameter
vector. Note that the linear constraint representation is not unique; that is,
for any nonsingular matrix $\mathbf{O}\in\mathbb{R}^{r\times r}$,
$\mathbf{C}\boldsymbol{\theta}=(\mathbf{C}\mathbf{O})(\mathbf{O}^{-1}\boldsymbol{\theta})$.
Thus, we require that $\mathbf{C}$ is a tall orthonormal matrix, i.e.,
$\mathbf{C}^{\top}\mathbf{C}=\mathbf{I}_{r}$. Otherwise, we can apply singular
value decomposition or QR decomposition to $\mathbf{C}$.

For simplicity, we consider the case where $\boldsymbol{\gamma}=\mathbf{0}$.
In this case, we can rewrite the VAR($d$) model in (2) to

$\mathbf{y}_{t}=(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{C}\boldsymbol{\theta}+\boldsymbol{\varepsilon}_{t}.$
(16)

Similarly to the discussions in Section 2.1, for the general stationary time
series $\mathbf{y}_{t}$ which may not strictly follow a linear-restricted VAR
process, the parameter of interest $\boldsymbol{\theta}^{*}$ is defined as

$\boldsymbol{\theta}^{*}=\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\mathbb{R}^{r}}\mathbb{E}\left[\|\mathbf{y}_{t}-(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{C}\boldsymbol{\theta}\|_{2}^{2}\right].$
(17)

It implies the linear-restricted version of the Yule–Walker equation

$\boldsymbol{\theta}^{*}=[\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{\Sigma}_{0})\mathbf{C}]^{-1}\mathbf{C}^{\top}\textup{vec}(\mathbf{\Sigma}_{1}^{\top})~{}~{}\text{or}~{}~{}\mathbf{C}^{\top}\textup{vec}(\mathbf{\Sigma}_{1}^{\top})=\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{\Sigma}_{0})\mathbf{C}\boldsymbol{\theta}^{*},$
(18)

provided that
$\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{\Sigma}_{0})\mathbf{C}$ is
invertible.

Let
$\mathbf{\Omega}=\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{\Sigma}_{0})\mathbf{C}\in\mathbb{R}^{r\times
r}$ and
$\boldsymbol{\omega}=\mathbf{C}^{\top}\textup{vec}(\mathbf{\Sigma}_{1}^{\top})\in\mathbb{R}^{r}$.
For $\mathbf{\Omega}$ and $\boldsymbol{\omega}$, suppose that we have robust
estimators $\widetilde{\mathbf{\Omega}}$ and
$\widetilde{\boldsymbol{\omega}}$. Then, we consider the constrained
Yule–Walker estimator

$\boldsymbol{\widehat{\theta}}=\arg\min\|\boldsymbol{\theta}\|_{\infty}~{}~{}~{}\textup{subject
to}~{}~{}\|\mathbf{\widetilde{\Omega}}\boldsymbol{\theta}-\boldsymbol{\widetilde{\omega}}\|_{\infty}\leq\lambda,$
(19)

where $\lambda$ is the constraint parameter and the $\ell_{\infty}$ norm
provides both element-wise regularization on $\boldsymbol{\theta}$ and
element-wise constraint on the sample version of Yule–Walker equation.
Decompose the restriction matrix
$\mathbf{C}=[\mathbf{B}_{1}^{\top},\dots,\mathbf{B}_{p}^{\top}]^{\top}$, where
each $\mathbf{B}_{i}$ is a $pd\times r$ matrix. Let
$\mathbf{B}=[\mathbf{B}_{1},\dots,\mathbf{B}_{p}]\in\mathbb{R}^{pd\times pr}$
and then the estimator of $\mathbf{A}$ is
$\mathbf{\widehat{A}}=(\mathbf{I}_{p}\otimes\boldsymbol{\widehat{\theta}})^{\top}\mathbf{B}^{\top}$.

The linear-restricted model encompasses several important recent developments
in high-dimensional vector autoregression.

######  Example 3 (Banded VAR).

The banded VAR(1) proposed by Guo et al., (2016) has the following banded
coefficient structure:

$\mathbf{A}_{ij}=0,~{}~{}\text{for all }|i-j|>k_{0},$ (20)

where $k_{0}$ is called the bandwidth parameter. The banded structure is a
special sparse structure, so we can formulate it to a linear constraint
$\textup{vec}(\mathbf{A}^{\top})=\mathbf{C}\boldsymbol{\theta}$, where each
column of $\mathbf{C}$ is a coordinate vector corresponding to one nonzero
index in $\textup{vec}(\mathbf{A}^{\top})$.

######  Example 4 (Network VAR).

Zhu et al., (2017) proposed the network VAR model for the network time series
data $\mathbf{y}_{t}$, which has an observable network structure with the
adjacency matrix $\mathbf{W}$. The network VAR model assumes the form

$\mathbf{y}_{t}=\beta_{1}\mathbf{y}_{t-1}+\beta_{2}\mathbf{W}\mathbf{y}_{t-1}+\boldsymbol{\varepsilon}_{t}.$
(21)

In other words, it is a VAR(1) model with
$\mathbf{A}=\beta_{1}\mathbf{I}_{p}+\beta_{2}\mathbf{W}$, and the associated
coefficient vector can be written as
$\textup{vec}(\mathbf{A}^{\top})=\mathbf{C}\boldsymbol{\theta}$, where
$\mathbf{C}=(p^{-1/2}\textup{vec}(\mathbf{I}_{p}),\|\mathbf{W}\|_{\textup{F}}^{-1}\textup{vec}(\mathbf{W}^{\top}))$,
and
$\boldsymbol{\theta}=(p^{1/2}\beta_{1},\|\mathbf{W}\|_{\textup{F}}\beta_{2})^{\top}$.
In this parameterization, $\mathbf{C}$ is an orthonormal matrix.

For a general zero-mean and covariance stationary time series
$\\{\mathbf{y}_{t}\\}_{t=1}^{T}$, if we model the data by the linear-
restricted VAR($d$) model in (16) and apply the linear-restricted Yule–Walker
estimator in (19) with the robust estimators $\mathbf{\widetilde{\Omega}}$ and
$\boldsymbol{\widetilde{\omega}}$, we have the upper bounds for the estimation
error between the estimated parameter vector $\boldsymbol{\widehat{\theta}}$
and the true vector $\boldsymbol{\theta}^{*}$ defined in (17).

######  Proposition 3.

For the constrained Yule–Walker estimator in (19), suppose that
$\|\mathbf{\widetilde{\Omega}}-\mathbf{\Omega}\|_{1,\infty}\leq\zeta_{1}$,
$\|\boldsymbol{\widetilde{\omega}}-\boldsymbol{\omega}\|_{\infty}\leq\zeta_{2}$,
and $\lambda\geq\zeta_{1}\|\boldsymbol{\theta}^{*}\|_{\infty}+\zeta_{2}$.
Then,

$\|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta}^{*}\|_{\infty}\leq
2\lambda\|\mathbf{\Omega}^{-1}\|_{1,\infty}.$ (22)

This proposition presents the $\ell_{\infty}$ norm estimation error bounds of
the low-dimensional parameter $\boldsymbol{\widehat{\theta}}$, which is a new
result for the linear-restricted VAR models. Similarly to Proposition 1, the
upper bound is proportional to $\|\mathbf{\Omega}^{-1}\|_{1,\infty}$ that
could be related to the dimension $p$, depending on the specific linear
restrictions. The $\ell_{\infty}$ upper bound directly implies that
$\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\|_{2}\leq
2\sqrt{r}\lambda\|\mathbf{\Omega}^{-1}\|_{1,\infty}$, and some sharper bounds
may be obtained based on the the special structure of $\mathbf{C}$. More
concrete results and discussions are provided in Section 3.3.

######  Remark 1.

For the upper bound result in Proposition 3, the terms
$\|\boldsymbol{\theta}^{*}\|_{\infty}$ and
$\|\mathbf{\Omega}^{-1}\|_{1,\infty}$ may vary or even diverge with the
dimension $p$. Indeed, the term $\|\boldsymbol{\theta}^{*}\|_{\infty}$ for the
network VAR model in Example 3 naturally diverges at a rate $\sqrt{p}$. The
explicit rate of $\lambda\|\mathbf{\Omega}^{-1}\|_{1,\infty}$ depends on the
specific linear restriction structure and further assumptions. More
discussions about these terms are given in Section 3.3 for some specific
models.

The deterministic upper bounds in Propositions 1-3 imply that the performance
of parameter estimation of the constrained Yule–Walker estimators hinges on
the accuracy of autocovariance matrix estimation. It suffices to find reliable
robust autocovariance estimators, which will be introduced in the next
section.

##  3 Robust Autocovariance Estimation

###  3.1 Element truncation estimator

Fan et al., (2021) and Ke et al., (2019) proposed a simple robust estimation
approach via appropriate truncation on data and showed that the truncated
covariance estimator can achieve the optimal rate as that under the sub-
Gaussian distribution for independent samples. We first apply the truncation
method to estimate the elements of autocovariance matrices.

We consider the element-wise truncated data
$\mathbf{y}_{t}^{\text{E}}(\tau)=(y_{1t}(\tau),y_{2t}(\tau),\dots,y_{pt}(\tau))^{\top}$,
where $y_{it}(\tau)=\text{sign}(y_{it})(\tau\wedge|y_{it}|)$, for $1\leq i\leq
p$, and $\tau>0$ is the truncation parameter. Based on this truncation scheme,
$\mathbf{\Gamma}_{\ell}$ can be estimated by
$\mathbf{\widetilde{\Gamma}}_{\ell}^{\text{E}}(\tau)=T^{-1}\sum_{t=1}^{T}\mathbf{y}_{t}^{\text{E}}(\tau)\mathbf{y}_{t-\ell}^{\text{E}}(\tau)^{\top}$,
for any integer $\ell\geq 0$. The corresponding autocovariance estimators are

$\widetilde{\mathbf{\Sigma}}_{1}^{\text{E}}(\tau)=[\widetilde{\mathbf{\Gamma}}^{\text{E}}_{1}(\tau),\widetilde{\mathbf{\Gamma}}^{\text{E}}_{2}(\tau),\dots,\widetilde{\mathbf{\Gamma}}^{\text{E}}_{d}(\tau)],$
(23)

and

$\widetilde{\mathbf{\Sigma}}_{0}^{\text{E}}(\tau)=\begin{bmatrix}\widetilde{\mathbf{\Gamma}}^{\text{E}}_{0}(\tau)&\widetilde{\mathbf{\Gamma}}^{\text{E}}_{1}(\tau)&\cdots&\widetilde{\mathbf{\Gamma}}^{\text{E}}_{d-1}(\tau)\\\
\widetilde{\mathbf{\Gamma}}^{\text{E}}_{-1}(\tau)&\widetilde{\mathbf{\Gamma}}^{\text{E}}_{0}(\tau)&\cdots&\widetilde{\mathbf{\Gamma}}^{\text{E}}_{d-2}(\tau)\\\
\vdots&\vdots&\ddots&\vdots\\\
\widetilde{\mathbf{\Gamma}}^{\text{E}}_{-d+1}(\tau)&\widetilde{\mathbf{\Gamma}}^{\text{E}}_{-d}(\tau)&\cdots&\widetilde{\mathbf{\Gamma}}^{\text{E}}_{0}(\tau)\end{bmatrix},$
(24)

where
$\widetilde{\mathbf{\Gamma}}^{\text{E}}_{-j}(\tau)=\widetilde{\mathbf{\Gamma}}^{\text{E}}_{j}(\tau)^{\top}$,
for $j\geq 1$. The element-wise truncation can control the deviation of
$y_{it}y_{jt}$ or $y_{it}y_{j,t-\ell}$ and the truncation parameter $\tau$
allows us to balance the tradeoff between the truncation bias and robustness.

Note that no data generating mechanism assumption is imposed, and we adopt the
$\alpha$-mixing condition to quantify the serial dependency. For any
stochastic process $\\{\mathbf{y}_{t}\\}_{t=-\infty}^{\infty}$, the lag-$\ell$
$\alpha$-mixing dependence coefficient is defined as
$\alpha(\ell):=\sup_{s\in\mathbb{Z}}\alpha(\\{\mathbf{y}_{t}\\}_{t=-\infty}^{s},\\{\mathbf{y}_{t}\\}_{t=s+\ell}^{\infty})$,
where

$\alpha(\\{\mathbf{y}_{t}\\}_{t=-\infty}^{s},\\{\mathbf{y}_{t}\\}_{t=r}^{\infty})=\sup|\mathbb{P}(A\cap
B)-\mathbb{P}(A)\mathbb{P}(B)|,$ (25)

the supremum being taken over all events
$A\in\sigma(\\{\mathbf{y}_{t}\\}_{t=-\infty}^{s})$ and
$B\in\sigma(\\{\mathbf{y}_{t}\\}_{t=r}^{\infty})$, and $\sigma(\cdot)$ is the
sigma field generated by the process. To derive the theoretical guarantees of
the autocovariance estimation, we have the following assumptions.

######  Assumption 1.

The process $\\{\mathbf{y}_{t}\\}$ is weakly stationary and $\alpha$-mixing
with the mixing coefficients $\alpha(\ell)=O(r^{\ell})$, where $r=r(p)$ is a
sequence possibly depending on $p$ such that $0\leq r\leq\bar{r}$ for some
constant $\bar{r}<1$.

######  Assumption 2.

For $1\leq i\leq p$, $\mathbb{E}[|y_{it}|^{2+2\epsilon}]\leq M_{2+2\epsilon}$,
for some $\epsilon\in(0,1]$.

Assumption 1 states the weak stationarity and geometrically decayed
$\alpha$-mixing, rather than assuming that the true data generating process
(DGP) of $\mathbf{y}_{t}$ is a VAR($d$) process. If $\mathbf{y}_{t}$ truly
follows a VAR($d$) model in (1), Assumption 1 holds by the stationarity and
geometric ergodicity of the VAR process; see Proposition 2 in Liebscher,
(2005). Assumption 2 relaxes the commonly-used sub-Gaussian condition in the
existing literature to the bounded $(2+2\epsilon)$-th moment condition, where
$M_{2+2\epsilon}$ may diverge to infinity with the dimension $p$. For the
strong mixing time series, we denote $n_{\text{eff}}=T/\log(T)^{2}$ as the
effective sample size representing the number of effectively independent
samples from $T$ observations. Note that $n_{\text{eff}}$ is of almost the
same rate as $T$ since $T^{1-\delta}=o(n_{\text{eff}})$ for any $\delta>0$.

######  Proposition 4.

Under Assumptions 1 and 2, if $T\gtrsim\log(p^{2}d)$ and

$\tau\asymp\left[\frac{M_{2+2\epsilon}n_{\textup{eff}}}{\log(p^{2}d)}\right]^{\frac{1}{2+2\epsilon}},$
(26)

then, with probability at least $1-C\exp[-C\log(T)\log(p^{2}d)]$,

$\begin{split}&\|\widetilde{\mathbf{\Sigma}}_{0}^{\textup{E}}(\tau)-\mathbf{\Sigma}_{0}\|_{\infty}\lesssim\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p^{2}d)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}},\\\
\text{and~{}~{}}&\|\widetilde{\mathbf{\Sigma}}_{1}^{\textup{E}}(\tau)-\mathbf{\Sigma}_{1}\|_{\infty}\lesssim\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p^{2}d)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}}.\end{split}$
(27)

This proposition presents the $\ell_{\infty}$ bounds for the element-wise
truncation autocovariance estimators
$\mathbf{\widetilde{\Sigma}}^{\text{E}}_{0}(\tau)$ and
$\mathbf{\widetilde{\Sigma}}^{\text{E}}_{1}(\tau)$. If the truncation
parameter $\tau$ is chosen appropriately and the moment bound
$M_{2+2\epsilon}$ is fixed, the convergence rate of autocovariance matrix
estimation in the $\ell_{\infty}$ norm is
$O([\log(p^{2}d)/n_{\text{eff}}]^{\epsilon/(1+\epsilon)})$. When $\epsilon=1$,
the data has a bounded fourth moment and the convergence rates of the
autocovariance matrices scale as $\sqrt{\log(p^{2}d)/n_{\text{eff}}}$. When
$\epsilon\in(0,1)$, the convergence rates have a smooth phase transition
phenomenon, decreasing from $[\log(p^{2}d)/n_{\text{eff}}]^{1/2}$ to
$[\log(p^{2}d)/n_{\text{eff}}]^{\epsilon/(1+\epsilon)}$.

######  Remark 2.

For robust covariance estimators based on data truncation or shrinkage, the
existing theoretical analysis in Fan et al., (2021) and Ke et al., (2019)
focused on the data with a bounded fourth moment. Our results in Proposition 4
relax the fourth moment condition to the $(2+2\epsilon)$-th moment condition
and can effectively handle a much larger class of distributions. In addition,
Avella-Medina et al., (2018) proposed rank-based and adaptive Huber regression
methods for robust covariance estimation and obtained a similar phase
transition in the upper bound with respect to $\epsilon$.

For the weakly sparse VAR model in Example 2,
$\mathbf{\widetilde{\Sigma}}_{0}^{\text{E}}(\tau)$ and
$\mathbf{\widetilde{\Sigma}}_{1}^{\text{E}}(\tau)$ can be used as robust
autocovariance estimators in the constrained Yule–Walker estimation in (12)
with $\mathcal{R}(\cdot)=\|\cdot\|_{1}$ and
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\infty}$. The resulting estimator for the
$i$-th row of $\mathbf{A}$ is denoted as
$\widehat{\mathbf{a}}_{i}(\lambda,\tau)$ and the sparse estimator is
$\widehat{\mathbf{A}}_{\text{S}}(\lambda,\tau)=[\widehat{\mathbf{a}}_{1}(\lambda,\tau),\dots,\widehat{\mathbf{a}}_{p}(\lambda,\tau)]^{\top}$.

Define an $\ell_{q}$-“ball” with radius $s_{q}$ as
$\mathbb{B}_{q}(s_{q})=\\{\mathbf{A}\in\mathbb{R}^{p\times pd}:\max_{1\leq
i\leq p}\sum_{j=1}^{pd}|\mathbf{A}_{ij}|^{q}\leq s_{q}\\}$. Note that when
$q=0$, $\mathbb{B}_{0}(s_{0})$ is the set of all $p$-by-$pd$ matrices whose
rows are at most $s_{0}$-sparse. For $q>0$, $\mathbb{B}_{q}(s_{q})$ requires
that the absolute values of entries in $\mathbf{A}$ decay sufficiently fast,
which is more general than the exact sparsity assumption. Given the weak
sparsity in the rows of $\mathbf{A}^{*}$, it is natural to assume that
$\|\mathbf{A}^{*\top}\|_{1,\infty}$ is bounded.

We are ready to present the theoretical properties of
$\widehat{\mathbf{A}}_{\text{S}}(\lambda,\tau)$.

######  Theorem 1 (Sparse VAR upper bounds).

Suppose that $\mathbf{A}^{*}\in\mathbb{B}_{q}(s_{q})$ for some $q\in[0,1)$,
$\|\mathbf{A}^{*\top}\|_{1,\infty}\leq C<\infty$, and Assumptions 1 and 2
hold. If $T\gtrsim\log(p^{2}d)$,

$\tau\asymp\left[\frac{M_{2+2\epsilon}n_{\textup{eff}}}{\log(p^{2}d)}\right]^{\frac{1}{2+2\epsilon}}\quad\text{and}\quad\lambda\asymp\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p^{2}d)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}},$
(28)

then, with probability at least $1-C\exp[-C\log(T)\log(p^{2}d)]$,

$\begin{split}\|\mathbf{\widehat{A}}_{\textup{S}}(\lambda,\tau)-\mathbf{A}^{*}\|_{\infty}&\lesssim\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p^{2}d)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}},\\\
\|\mathbf{\widehat{A}}^{\top}_{\textup{S}}(\lambda,\tau)-\mathbf{A}^{*\top}\|_{1,\infty}&\lesssim
s_{q}\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}^{1-q}\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p^{2}d)}{n_{\textup{eff}}}\right]^{\frac{(1-q)\epsilon}{1+\epsilon}},\\\
\text{and}~{}~{}\|\mathbf{\widehat{A}}^{\top}_{\textup{S}}(\lambda,\tau)-\mathbf{A}^{*\top}\|_{2,\infty}&\lesssim\sqrt{s_{q}}\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}^{1-\frac{q}{2}}\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p^{2}d)}{n_{\textup{eff}}}\right]^{\left(1-\frac{q}{2}\right)\frac{\epsilon}{1+\epsilon}}.\end{split}$
(29)

This theorem presents non-asymptotic estimation upper bounds in various matrix
norms under the weakly sparse structure. If both $M_{2+2\epsilon}$ and
$\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}$ are fixed, the $\|\cdot\|_{\infty}$,
$\|\cdot\|_{1,\infty}$, and $\|\cdot\|_{2,\infty}$ upper bounds scale as
$[\log(p^{2}d)/n_{\text{eff}}]^{\epsilon/(1+\epsilon)}$,
$s_{q}[\log(p^{2}d)/n_{\text{eff}}]^{(1-q)\epsilon/(1+\epsilon)}$, and
$\sqrt{s_{q}}[\log(p^{2}d)/n_{\text{eff}}]^{(1-q/2)\epsilon/(1+\epsilon)}$,
respectively. If the data have a bounded fourth moment with $\epsilon=1$, the
first two rates of convergence match those of Gaussian sparse VAR in Han et
al., (2015). When the true model has a strict sparsity structure with $q=0$
and sparsity level $s_{0}$, based on the $\|\cdot\|_{2,\infty}$ bound, the
sample size requirement is $n_{\text{eff}}\gtrsim s_{0}\log(p^{2}d)$. That is,
if the sparsity level in each row of $\mathbf{A}^{*}$ is finite, the dimension
$p$ is allowed to be exponentially large compared to $n_{\text{eff}}$. If the
data do not have a bounded fourth moment but only a bounded $(2+2\epsilon)$-th
moment for some $\epsilon\in(0,1)$, the proposed estimator is still consistent
but the rate of convergence in the $\ell_{2,\infty}$ norm decreases from
$\sqrt{s_{q}}[\log(p^{2}d)/n_{\text{eff}}]^{1/2}$ to
$\sqrt{s_{q}}[\log(p^{2}d)/n_{\text{eff}}]^{\epsilon/(1+\epsilon)}$.

######  Remark 3.

The boundedness of $\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}$ is widely assumed
in the existing literature of the Dantzig selector for sparse VAR models; see,
for example, Han et al., (2015) and Wu and Wu, (2016). Indeed, if
$\mathbf{y}_{t}$ follows a VAR(1) process in (1), the covariance of
$\mathbf{y}_{t}$ has an explicit form
$\mathbf{\Sigma}_{0}=\sum_{i=0}^{\infty}(\mathbf{A}^{*})^{i}\mathbf{\Sigma}_{\boldsymbol{\varepsilon}}(\mathbf{A}^{*\top})^{i}$,
where $\mathbf{\Sigma}_{\boldsymbol{\varepsilon}}$ is the covariance matrix of
$\boldsymbol{\varepsilon}_{t}$, and $\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}$
can be verified to be independent of $p$ given that $\mathbf{A}^{*}$ and
$\mathbf{\Sigma}_{\boldsymbol{\varepsilon}}$ satisfy certain structures. For
example, if both $\mathbf{A}^{*}$ and
$\mathbf{\Sigma}_{\boldsymbol{\varepsilon}}$ are block diagonal with a fixed
block size, $\mathbf{\Sigma}_{0}^{-1}$ is also block diagonal and its
$\ell_{1,\infty}$ norm is independent of $p$. If $\mathbf{A}^{*}$ has a banded
structure as in Example 3, $(\mathbf{\Sigma}_{0}^{-1})_{ij}$ shrinks to zero
exponentially as $|i-j|$ increases to infinity, and hence the boundedness of
$\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}$ can be expected.

######  Remark 4.

The robust estimation of sparse linear regression has been investigated in Fan
et al., (2017), Sun et al., (2020), and Fan et al., (2021), among many others.
We highlight some main differences between conventional linear regression and
vector autoregression in terms of robust estimation. First, for linear
regression models, the heavy-tailed distribution condition is considered on
the random errors; that is, the response is conditionally heavy-tailed given
the predictors. In our analysis, as no data generating mechanism assumption is
imposed, the heavy-tailed distributional assumption, i.e. the moment condition
in Assumption 2, is directly considered on the observed data $\mathbf{y}_{t}$.
Second, when the design matrix of the linear regression model satisfies some
regulatory conditions, the $(1+\epsilon)$-th moment condition is considered
for the random errors in Sun et al., (2020). However, for the VAR models, as
both predictor and response are simultaneously heavy-tailed, a more stringent
$(2+2\epsilon)$-th moment condition has to be imposed in our theoretical
analysis. Though the distributional assumptions are different in these two
problems, the similar phase transition phenomenon can be obtained for two
regimes $\epsilon\in(0,1)$ and $\epsilon\geq 1$.

###  3.2 Vector truncation estimator

In this subsection, we propose and study an autcovariance estimator that is
robust in the operator norm. Intuitively, in order to achieve the operator
norm robustness, we need to control the spectrum of
$\mathbf{x}_{t}\mathbf{x}_{t}^{\top}$ and
$\mathbf{y}_{t}\mathbf{x}_{t}^{\top}$. Note that
$\|\mathbf{x}_{t}\mathbf{x}_{t}^{\top}\|_{\textup{op}}=\|\mathbf{x}_{t}\|_{2}^{2}$
and
$\|\mathbf{y}_{t}\mathbf{x}_{t}^{\top}\|_{\textup{op}}=\|\mathbf{y}_{t}\|_{2}\|\mathbf{x}_{t}\|_{2}$.
Therefore, we propose the truncation method to the whole vector
$\mathbf{y}_{t}$ in the $\ell_{2}$ norm. To construct robust estimators for
$\mathbf{\Sigma}_{0}$ and $\mathbf{\Sigma}_{1}$ of a VAR($d$) model, we
consider the vector-truncated responses and predictors
$\mathbf{y}^{\text{V}}_{t}(\tau_{1})=(\tau_{1}\wedge\|\mathbf{y}_{t}\|_{2})\mathbf{y}_{t}/\|\mathbf{y}_{t}\|_{2}$
and
$\mathbf{x}^{\text{V}}_{t}(\tau_{2})=(\tau_{2}\wedge\|\mathbf{x}_{t}\|_{2})\mathbf{x}_{t}/\|\mathbf{x}_{t}\|_{2}$,
for $1\leq t\leq T$, where two truncation parameters $\tau_{1}$ and $\tau_{2}$
are adopted as $\mathbf{y}_{t}$ and $\mathbf{x}_{t}$ are of different
dimension when $d>1$. For the case with $d=1$, we can use a single parameter
by setting $\tau_{1}=\tau_{2}$. Based on the vector truncation of the data,
the corresponding truncation autocovariance estimators are defined as

$\widetilde{\mathbf{\Sigma}}^{\text{V}}_{0}(\tau_{2})=\frac{1}{T}\sum_{t=1}^{T}\mathbf{x}_{t}^{\text{V}}(\tau_{2})\mathbf{x}_{t}^{\text{V}}(\tau_{2})^{\top},~{}~{}\text{and}~{}~{}\widetilde{\mathbf{\Sigma}}^{\text{V}}_{1}(\tau_{1},\tau_{2})=\frac{1}{T}\sum_{t=1}^{T}\mathbf{y}_{t}^{\text{V}}(\tau_{1})\mathbf{x}_{t}^{\text{V}}(\tau_{2})^{\top}.$
(30)

######  Remark 5.

Fan et al., (2021) proposed a robust covariance estimator based on the vector-
wise truncation in the $\ell_{4}$ norm and investigated its theoretical
properties under the finite fourth moment condition. However, as we are
interested in the relaxed setting with the finite $(2+2\epsilon)$-th moment,
the truncation in the $\ell_{4}$ norm is not appropriate in our setting. Ke et
al., (2019) proposed the spectrum-wise truncation for robust estimation of
covariance matrix, and our method shares the same idea.

For the vector-wise truncation estimator, we adopt the $\beta$-mixing
condition to quantify the serial dependency. Specifically, for any stochastic
process $\\{\mathbf{y}_{t}\\}_{t=-\infty}^{\infty}$, the lag-$\ell$
$\beta$-mixing dependence coefficient is defined as
$\beta(\ell):=\sup_{s\in\mathbb{Z}}\beta(\\{\mathbf{y}_{t}\\}_{t=-\infty}^{s},\\{\mathbf{y}_{t}\\}_{t=s+\ell}^{\infty})$,
where

$\beta(\\{\mathbf{y}_{t}\\}_{t=-\infty}^{s},\\{\mathbf{y}_{t}\\}_{t=r}^{\infty})=\sup\frac{1}{2}\sum_{i=1}^{I}\sum_{j=1}^{J}\left|\mathbb{P}(A_{i}\cap
B_{j})-\mathbb{P}(A_{i})\mathbb{P}(B_{j})\right|,$ (31)

the supremum being taken over all pairs of partitions
$\\{A_{1},\dots,A_{I}\\}$ and $\\{B_{1},\dots,B_{J}\\}$ such that
$A_{i}\in\sigma(\\{\mathbf{y}_{t}\\}_{t=-\infty}^{s})$ and
$B_{j}\in\sigma(\\{\mathbf{y}_{t}\\}_{t=r}^{\infty})$ . To study the
theoretical properties of the vector-wise truncation estimators, we need the
following assumptions.

######  Assumption 3.

The process $\\{\mathbf{y}_{t}\\}$ is weakly stationary and $\beta$-mixing
with the mixing coefficients $\beta(\ell)=O(r^{\ell})$, where $r=r(p)$ is a
sequence possibly depending on $p$ such that $0\leq r\leq\bar{r}$ for some
constant $\bar{r}<1$.

######  Assumption 4.

For any $\mathbf{v}\in\mathbb{R}^{pd}$ such that $\|\mathbf{v}\|_{2}=1$ and
some $\delta>0$,
$\mathbb{E}[(\mathbf{x}_{t}^{\top}\mathbf{v})^{2+2\epsilon+\delta}]\leq
M_{2+2\epsilon}$.

Note that the $\beta$-mixing condition in Assumption 3 is slightly stronger
than the $\alpha$-mixing condition in Assumption 1 for the element-wise
truncation estimators. Here we adopt the $\beta$-mixing condition in order to
apply the Bernstein-type inequality for $\beta$-mixing random matrices
developed by Banna et al., (2016). For the VAR process with independent and
identically distributed $\boldsymbol{\varepsilon}_{t}$, the absolute
regularity with geometrically decayed $\beta$-mixing coefficients is
equivalent to the geometric ergodicity (Liebscher,, 2005). The bounded
$(2+2\epsilon+\delta)$-th moment condition in Assumption 4 is also slightly
stronger than the bounded $(2+2\epsilon)$-th moment condition in Assumption 2.
Here $\delta$ is only for the technical purpose and can be arbitrarily small.

######  Proposition 5.

Under Assumptions 3 and 4, if $T\gtrsim pd$,

$\tau_{1}\asymp(p^{\epsilon}M_{2+2\epsilon}n_{\textup{eff}})^{\frac{1}{2+2\epsilon}}~{}~{}\text{and}~{}~{}\tau_{2}\asymp(p^{\epsilon}d^{\epsilon}M_{2+2\epsilon}n_{\textup{eff}})^{\frac{1}{2+2\epsilon}},$
(32)

then, with probability at least $1-C\exp[-C\log(T)]$,

$\begin{split}\|\mathbf{\widetilde{\Sigma}}_{0}^{\textup{V}}(\tau_{2})-\mathbf{\Sigma}_{0}\|_{\textup{op}}&\lesssim\left(\frac{pM_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}},\\\
\text{and}~{}~{}\|\mathbf{\widetilde{\Sigma}}_{1}^{\textup{V}}(\tau_{1},\tau_{2})-\mathbf{\Sigma}_{1}\|_{\textup{op}}&\lesssim\left(\frac{pM_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}}.\end{split}$
(33)

This proposition presents the operator norm bounds for the vector truncation
autocovariance estimators
$\mathbf{\widetilde{\Sigma}}_{0}^{\text{V}}(\tau_{2})$ and
$\mathbf{\widetilde{\Sigma}}_{1}^{\text{V}}(\tau_{1},\tau_{2})$. Similarly to
the element truncation estimator in Section 3.1, we define
$n_{\text{eff}}=T/\log(T)^{2}$. If the moment bound $M_{2+2\epsilon}$ is
fixed, the operator norm convergence rates of both autocovariance matrix
estimators scale as $(pd/n_{\text{eff}})^{\epsilon/(1+\epsilon)}$. When the
time series data has a bounded fourth moment, that is $\epsilon=1$, our
results are almost the same as those of covariance matrix estimators for
i.i.d. data in Ke et al., (2019) and Fan et al., (2021). Smooth transition on
the operator norm convergence rate is observed when the fourth moment
condition is relaxed to the $(2+2\epsilon)$-th moment condition.

Based on the rates of the vector truncation autocovariance matrices, we derive
the estimation rate of the nuclear norm constrained Yule–Walker estimator for
the weakly low-rank VAR model. Define an $\ell_{q}$-“ball” for the singular
values with radius $r_{q}$ as
$\widetilde{\mathbb{B}}_{q}(r_{q})=\\{\mathbf{A}\in\mathbb{R}^{p\times
p}:\sum_{i=1}^{p}\sigma_{i}^{q}(\mathbf{A})\leq r_{q}\\}$. When the singular
values of $\mathbf{A}^{*}$ are weakly sparse, it is natural to assume that
$\|\mathbf{A}^{*}\|_{\textup{nuc}}$, the sum of singular values of
$\mathbf{A}^{*}$, is bounded. Denote by
$\mathbf{\widehat{A}}_{\text{RR}}(\lambda,\tau_{1},\tau_{2})$ the reduced-rank
constrained Yule–Walker estimator with the regularizer
$\mathcal{R}(\cdot)=\|\cdot\|_{\textup{nuc}}$, the constraint function
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\textup{op}}$, the constraint parameter
$\lambda$, and the vector truncation autocovariance estimators in (30) with
truncation parameters $\tau_{1}$ and $\tau_{2}$.

######  Theorem 2 (Reduced-rank VAR upper bounds).

Suppose that $\mathbf{A}^{*}\in\widetilde{\mathbb{B}}_{q}(r_{q})$ for some
$q\in[0,1)$, $\|\mathbf{A}^{*}\|_{\textup{nuc}}\leq C<\infty$, and Assumptions
3 and 4 hold. If $T\gtrsim p$,

$\tau_{1}\asymp(p^{\epsilon}M_{2+2\epsilon}n_{\textup{eff}})^{\frac{1}{2+2\epsilon}},~{}~{}\tau_{2}\asymp(p^{\epsilon}d^{\epsilon}M_{2+2\epsilon}n_{\textup{eff}})^{\frac{1}{2+2\epsilon}},~{}~{}\text{and}~{}~{}\lambda\asymp\left(\frac{pM_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}},$
(34)

then, with probability at least $1-C\exp[-C\log(T)]$,

$\begin{split}\|\mathbf{\widehat{A}}_{\textup{RR}}(\lambda,\tau_{1},\tau_{2})-\mathbf{A}^{*}\|_{\textup{op}}&\lesssim\left(\frac{pM_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}},\\\
\|\mathbf{\widehat{A}}_{\textup{RR}}(\lambda,\tau_{1},\tau_{2})-\mathbf{A}^{*}\|_{\textup{nuc}}&\lesssim
r_{q}\|\mathbf{\Sigma}_{0}^{-1}\|^{1-q}_{\textup{op}}\left(\frac{pM_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{(1-q)\epsilon}{1+\epsilon}},\\\
\text{and}~{}~{}\|\mathbf{\widehat{A}}_{\textup{RR}}(\lambda,\tau_{1},\tau_{2})-\mathbf{A}^{*}\|_{\textup{F}}&\lesssim\sqrt{r_{q}}\|\mathbf{\Sigma}_{0}^{-1}\|^{1-\frac{q}{2}}_{\textup{op}}\left(\frac{pM_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\left(1-\frac{q}{2}\right)\frac{\epsilon}{1+\epsilon}}.\end{split}$
(35)

This theorem presents the non-asymptotic estimation upper bounds in the
operator norm, nuclear norm, and Frobenius norm, respectively. If
$M_{2+2\epsilon}$ is fixed, the convergence rates of
$\widehat{\mathbf{A}}_{\text{RR}}$ in the operator norm and Frobenius norm
scale as $(p/n_{\text{eff}})^{\epsilon/(1+\epsilon)}$ and
$\sqrt{r_{q}}(p/n_{\text{eff}})^{(1-q/2)\epsilon/(1+\epsilon)}$. Specifically,
when $\epsilon=1$, that is the time series have a bounded fourth moment, the
Frobenius norm convergence rate of the robust estimator is nearly the same as
those obtained by the standard nuclear norm penalized estimators for Gaussian
VAR model in Negahban and Wainwright, (2011) and Basu et al., (2019). When the
time series only has a bounded $(2+2\epsilon)$-th moment for some
$\epsilon\in(0,1)$, the estimation error rates decrease from
$\sqrt{r_{q}}(p/n_{\text{eff}})^{1/2}$ to
$(p/n_{\text{eff}})^{\epsilon/(1+\epsilon)}$. Note that though the estimation
convergence rates decrease, if $r_{q}$ is fixed, the sample size requirement
$n_{\text{eff}}\gtrsim pd$ for estimation consistency remains unchanged when
the moment condition is relaxed.

######  Remark 6.

The boundedness of $\|\mathbf{\Sigma}_{0}^{-1}\|_{\textup{op}}$ is equivalent
to that the smallest eigenvalue of $\mathbf{\Sigma}_{0}$ is bounded away from
zero. If $\mathbf{y}_{t}$ follows a stationary VAR model in (1), this
condition can be guaranteed if the smallest eigenvalue of
$\mathbf{\Sigma}_{\boldsymbol{\varepsilon}}$ is bounded away from zero, where
$\mathbf{\Sigma}_{\boldsymbol{\varepsilon}}$ is the covariance matrix of
$\boldsymbol{\varepsilon}_{t}$.

###  3.3 Linear-restricted truncation estimator

For the linear-restricted VAR model in (16), the linear transformations of
autocovariance matrices are defined as

$\mathbf{\Omega}=\mathbb{E}[\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t})(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{C}]~{}~{}\text{and}~{}~{}\boldsymbol{\omega}=\mathbb{E}[\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t})\mathbf{y}_{t}].$
(36)

Motivated by the element and vector-wise truncation in the previous
subsections, in order to robustly estimate each element of $\mathbf{\Omega}$,
we apply the truncation to
$(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{c}_{i}$, where each
$\mathbf{c}_{i}\in\mathbb{R}^{p^{2}d}$ is the $i$-th column of $\mathbf{C}$.
When each $\mathbf{c}_{i}$ is highly sparse, e.g., $\mathbf{c}_{i}$ is an
coordinate vector for the banded VAR in Example 3, the vector
$(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{c}_{i}$ is also highly
sparse. In this case,let $S_{i}$ be the non-zero index set of
$(\mathbf{I}_{p}\otimes\mathbf{x}_{t})\mathbf{c}_{i}$, and
$\mathbf{w}_{it}=[(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{c}_{i}]_{S_{i}}$
and $\mathbf{z}_{it}=(\mathbf{y}_{t})_{S_{i}}$ be the sub-vectors. It is
obvious to check that
$\mathbf{c}_{i}^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t})(\mathbf{I}_{p}\otimes\mathbf{x}_{t}^{\top})\mathbf{c}_{j}=\mathbf{w}_{it}^{\top}\mathbf{w}_{jt}$
and
$\mathbf{c}_{i}^{\top}(\mathbf{I}_{p}\otimes\mathbf{x}_{t})\mathbf{y}_{t}=\mathbf{w}_{it}^{\top}\mathbf{z}_{it}$,
for $1\leq i,j\leq r$.

Therefore, we consider the vector-wise truncation
$\widetilde{\mathbf{w}}_{it}(\tau_{1})=(\tau_{1}\wedge\|\mathbf{w}_{it}\|_{2})\mathbf{w}_{it}/\|\mathbf{w}_{it}\|_{2}$
and
$\widetilde{\mathbf{z}}_{it}(\tau_{2})=(\tau_{2}\wedge\|\mathbf{z}_{it}\|_{2})\mathbf{z}_{it}/\|\mathbf{z}_{it}\|_{2}$,
and the elements of $\mathbf{\Omega}$ and $\boldsymbol{\omega}$ can be
estimated by the linear-restricted truncated data

$\widetilde{\mathbf{\Omega}}_{ij}(\tau_{1})=\frac{1}{T}\sum_{t=1}^{T}\widetilde{\mathbf{w}}_{it}(\tau_{1})^{\top}\widetilde{\mathbf{w}}_{jt}(\tau_{1})~{}~{}\text{and}~{}~{}\widetilde{\boldsymbol{\omega}}_{i}(\tau_{1},\tau_{2})=\frac{1}{T}\sum_{t=1}^{T}\widetilde{\mathbf{w}}_{it}(\tau_{1})^{\top}\widetilde{\mathbf{z}}_{it}(\tau_{2}),$
(37)

for $1\leq i,j\leq r$.

######  Remark 7.

For the special matrix $\mathbf{C}$ whose columns are coordinate vectors
corresponding to the non-zero entries in $\textup{vec}(\mathbf{A}^{\top})$,
the linear-restricted estimators with $\tau_{1}=\tau_{2}$ are equivalent to
the linear transformations of element truncation estimators, namely
$\widetilde{\mathbf{\Omega}}(\tau_{1})=\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\widetilde{\mathbf{\Sigma}}_{0}^{\text{E}}(\tau_{1}))\mathbf{C}$
and
$\widetilde{\boldsymbol{\omega}}(\tau_{1},\tau_{1})=\mathbf{C}^{\top}\textup{vec}(\widetilde{\mathbf{\Sigma}}_{1}^{\text{E}}(\tau_{1})^{\top})$.

For the linear-restricted models, we consider the following moment conditions.

######  Assumption 5.

For $1\leq i\leq r$, $\mathbb{E}[\|\mathbf{w}_{it}\|_{2}^{2+2\epsilon}]\leq
M_{1,2+2\epsilon}$ and $\mathbb{E}[\|\mathbf{z}_{it}\|_{2}^{2+2\epsilon}]\leq
M_{2,2+2\epsilon}$.

The moment conditions in Assumption 5 can be viewed as an extension of the
element-wise moment condition in Assumption 2. In addition, due to the
normalization of the matrix $\mathbf{C}$, the moment bounds
$M_{1,2+2\epsilon}$ and $M_{2,2+2\epsilon}$ are allowed to vary with the
dimension $p$. More discussions on the moment bounds are presented below for
each specific linear-restricted model. For the linear-restricted model in (16)
with $r>1$, the general results of the autocovariance estimators are given as
follows.

######  Proposition 6.

Under Assumptions 1 and 5, if $T\gtrsim\log(r)$,

$\tau_{1}\asymp\left[\frac{M_{1,2+2\epsilon}n_{\text{eff}}}{\log(r)}\right]^{\frac{1}{2+2\epsilon}},~{}~{}\text{and}~{}~{}\tau_{2}\asymp\left[\frac{M_{2,2+2\epsilon}n_{\text{eff}}}{\log(r)}\right]^{\frac{1}{2+2\epsilon}},$
(38)

then, with probability at least $1-C\exp[-C\log(r)\log(T)]$,

$\begin{split}\|\widetilde{\mathbf{\Omega}}(\tau_{1})-\mathbf{\Omega}\|_{\infty}&\lesssim\left[\frac{M_{1,2+2\epsilon}^{1/\epsilon}\log(r)}{n_{\text{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}},\\\
\text{and}~{}~{}\|\widetilde{\boldsymbol{\omega}}(\tau_{1},\tau_{2})-\boldsymbol{\omega}\|_{\infty}&\lesssim\left[\frac{M_{1,2+2\epsilon}^{1/(2\epsilon)}M_{2,2+2\epsilon}^{1/(2\epsilon)}\log(r)}{n_{\text{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}}.\end{split}$
(39)

This proposition presents the estimation upper bounds for
$\widetilde{\mathbf{\Omega}}(\tau_{1})$ and
$\widetilde{\boldsymbol{\omega}}(\tau_{2})$ in the $\ell_{\infty}$ norm.
Following this proposition, we first consider the banded VAR model with the
bandwidth $k_{0}$ in Example 3, denote by
$\widehat{\boldsymbol{\theta}}_{\text{B}}(\lambda,\tau)$ and
$\widehat{\mathbf{A}}_{\text{B}}(\lambda,\tau)$ the linear-restricted
constrained Yule–Walker estimators with robust autocovariance estimators
$\widetilde{\mathbf{\Omega}}(\tau)$ and
$\widetilde{\boldsymbol{\omega}}(\tau,\tau)$. Note that as each
$\mathbf{c}_{i}$ is a coordinate vector, both $\mathbf{w}_{it}$ and
$\mathbf{z}_{it}$ are 1-dimensional, and the Assumption 5 reduces to the
element-wise moment bound in Assumption 2. The following estimation upper
bounds can be derived.

######  Theorem 3 (Banded VAR upper bounds).

For the banded VAR model with a bounded bandwidth $k_{0}$ in Example 3,
suppose that $\|\boldsymbol{\theta}^{*}\|_{\infty}\leq C<\infty$, and
Assumptions 1 and 2 hold. If $T\gtrsim\log(p)$,

$\tau\asymp\left[\frac{M_{2+2\epsilon}n_{\textup{eff}}}{\log(p)}\right]^{\frac{1}{1+\epsilon}},~{}~{}\text{and}~{}~{}\lambda\asymp\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p)}{n_{\textup{eff}}}\right]^{\frac{1}{1+\epsilon}},$
(40)

then, with probability at least $1-C\exp[-C\log(p)\log(T)]$,

$\begin{split}\|\boldsymbol{\widehat{\theta}}_{\textup{B}}(\lambda,\tau)-\boldsymbol{\theta}^{*}\|_{\infty}&\lesssim\|\mathbf{\Omega}^{-1}\|_{1,\infty}\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}},\\\
\|\mathbf{\widehat{A}}_{\textup{B}}(\lambda,\tau)-\mathbf{A}^{*}\|_{\textup{op}}&\lesssim\|\mathbf{\Omega}^{-1}\|_{1,\infty}\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}},\\\
\text{and}~{}~{}\|\boldsymbol{\widehat{\theta}}_{\textup{B}}(\lambda,\tau)-\boldsymbol{\theta}^{*}\|_{2}=\|\mathbf{\widehat{A}}_{\textup{B}}(\lambda,\tau)-\mathbf{A}^{*}\|_{\textup{F}}&\lesssim\sqrt{p}\|\mathbf{\Omega}^{-1}\|_{1,\infty}\left[\frac{M_{2+2\epsilon}^{1/\epsilon}\log(p)}{n_{\textup{eff}}}\right]^{\frac{\epsilon}{1+\epsilon}}.\end{split}$
(41)

When $M_{2+2\epsilon}$ and $\|\mathbf{\Omega}^{-1}\|_{1,\infty}$ are fixed,
the upper bounds in this theorem imply that asymptotically the rates of
convergence in the operator norm and Frobenius norm scale as
$[\log(p)/n_{\text{eff}}]^{\epsilon/(1+\epsilon)}$ and
$[p\log(p)/n_{\text{eff}}]^{\epsilon/(1+\epsilon)}$, respectively. In other
words, the sample size requirement in the operator norm is
$n_{\text{eff}}\gtrsim\log(p)$. When $\epsilon=1$, our results are comparable
and consistent with those in Guo et al., (2016) up to a logarithm factor.
However, in Guo et al., (2016), they imposed the i.i.d. condition and sub-
exponential tail condition for $\\{\boldsymbol{\varepsilon}_{t}\\}$ when
$p\gtrsim T$. Hence, under high-dimensional scaling, both the operator norm
and Frobenius norm convergence rates of the robust constrained Yule–Walker
estimator for heavy-tailed data with a bounded fourth moment are almost the
same as those of the standard ordinary least squares under the sub-exponential
tail condition. In sum, the convergence rates obtained under the sub-
exponential distribution can be achieved by the robust estimator under a much
relaxed fourth moment condition, and we also establish the estimation
consistency and the rates of convergence under the $(2+2\epsilon)$-th moment
condition.

######  Remark 8.

The stationarity of the VAR(1) model requires that the largest eigenvalue of
$\mathbf{A}^{*}$ in terms of absolute value is strictly smaller than one, and
hence the absolute value of $\boldsymbol{\theta}^{*}_{i}$ is expected to be
bounded. As
$\mathbf{\Omega}=\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{\Sigma}_{0})\mathbf{C}$
and the columns of $\mathbf{C}$ are all coordinate vectors,
$\|\mathbf{\Omega}^{-1}\|_{1,\infty}$ is bounded by
$\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}$, and the boundedness of
$\|\mathbf{\Sigma}_{0}^{-1}\|_{1,\infty}$ has been discussed in Remark 3.

Next, we consider another case of the linear-restricted VAR models, the
network VAR model in Example 4. In this model, $\mathbf{C}$ is a $p^{2}$-by-2
matrix with $\mathbf{c}_{1}=p^{-1/2}\textup{vec}(\mathbf{I}_{p})$ and
$\mathbf{c}_{2}=\|\mathbf{W}\|_{\textup{F}}^{-1}\textup{vec}(\mathbf{W}^{\top})$
being the restrictions associated with the nodal effect and network effect,
respectively. Accordingly,
$\boldsymbol{\theta}=(p^{1/2}\beta_{1},\|\mathbf{W}\|_{\textup{F}}\beta_{2})^{\top}$.
As $\mathbf{C}$ needs to be normalized fro identification purpose, we assume
that the unscaled parameters, $\beta_{1}$ and $\beta_{2}$ are constants
independent of $p$. In addition, we assume that each node in the network is
only connected to a fixed number of nodes, so the network adjacency matrix
satisfies $\|\mathbf{W}\|_{\textup{F}}\asymp\sqrt{p}$.

According to the linear restriction matrix $\mathbf{C}$ of the network VAR
model, it can be verified that $\mathbf{w}_{1t}=p^{-1/2}\mathbf{y}_{t-1}$,
$\mathbf{w}_{2t}=\|\mathbf{W}\|_{\textup{F}}^{-1}\mathbf{W}\mathbf{y}_{t-1}$,
and $\mathbf{z}_{1t}=\mathbf{z}_{2t}=\mathbf{y}_{t}$. Hence, if
$\max_{i=1,2}\mathbb{E}\|\mathbf{w}_{it}\|_{2}^{2+2\epsilon}=M_{2+2\epsilon}$
which is independent of $p$, then,
$\mathbb{E}\|\mathbf{z}_{1t}\|_{2}^{2+2\epsilon}=\mathbb{E}\|\mathbf{z}_{2t}\|_{2}^{2+2\epsilon}\leq
p^{1+\epsilon}M_{2+2\epsilon}$. For the network VAR model, denote by
$\widehat{\boldsymbol{\theta}}_{\text{N}}(\lambda,\tau_{1},\tau_{2})$ the
constrained Yule–Walker estimator with the constraint parameter $\lambda$ and
the autocovariance estimators $\widetilde{\mathbf{\Omega}}(\tau_{1})$ and
$\widetilde{\boldsymbol{\omega}}(\tau_{1},\tau_{2})$. The estimation upper
bounds of the network VAR model can be derived as follows.

######  Theorem 4 (Network VAR upper bounds).

For the network VAR model in Example 4 with
$\|\mathbf{W}\|_{\textup{F}}\asymp\sqrt{p}$ and
$\|\boldsymbol{\theta}^{*}\|_{\infty}\asymp\sqrt{p}$, under Assumptions 1 and
5 with $M_{1,2+2\epsilon}=M_{2+2\epsilon}$ and
$M_{2,2+2\epsilon}=p^{1+\epsilon}M_{2+2\epsilon}$, if

$\tau_{1}\asymp(M_{2+2\epsilon}n_{\textup{eff}})^{\frac{1}{2+2\epsilon}},~{}~{}\tau_{2}\asymp\sqrt{p}(M_{2+2\epsilon}n_{\textup{eff}})^{\frac{1}{2+2\epsilon}},~{}~{}\text{and}~{}~{}\lambda\asymp\sqrt{p}\left(\frac{M_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}},$
(42)

then, with probability at least $1-C\exp[-C\log(T)]$,

$\begin{split}\|\boldsymbol{\widehat{\theta}}_{\textup{N}}(\lambda,\tau_{1},\tau_{2})-\boldsymbol{\theta}^{*}\|_{\infty}&\lesssim\|\mathbf{\Omega}^{-1}\|_{1,\infty}\sqrt{p}\left(\frac{M_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}},\\\
\text{and}~{}~{}\|\mathbf{\widehat{A}}_{\textup{N}}(\lambda,\tau_{1},\tau_{2})-\mathbf{A}^{*}\|_{\textup{F}}=\|\boldsymbol{\widehat{\theta}}_{\textup{N}}(\lambda,\tau_{1},\tau_{2})-\boldsymbol{\theta}^{*}\|_{2}&\lesssim\|\mathbf{\Omega}^{-1}\|_{1,\infty}\sqrt{p}\left(\frac{M_{2+2\epsilon}^{1/\epsilon}}{n_{\textup{eff}}}\right)^{\frac{\epsilon}{1+\epsilon}}.\end{split}$
(43)

The upper bounds in this theorem are new for the network VAR model under the
non-asymptotic scheme. Due to the normalization for the linear restriction
matrix $\mathbf{C}$, both the true value $\|\boldsymbol{\theta}^{*}\|_{2}$ and
estimated value $\|\boldsymbol{\widehat{\theta}}\|_{2}$ will diverge at a rate
of $\sqrt{p}$ as the dimension $p$ increases to infinity. Hence, the
$\ell_{2}$ rate of convergence for the normalized coefficients
$p^{-1/2}\boldsymbol{\theta}$ scales as
$n_{\text{eff}}^{-\epsilon/(1+\epsilon)}$, independent of the dimension $p$,
and our result with $\epsilon=1$ is comparable to that obtained under the sub-
Gaussian condition in Zheng and Cheng, (2021). The convergence rates under a
bounded $(2+2\epsilon)$-th moment condition are established under a phase
transition pattern.

##  4 Minimax Lower Bounds

In this section, we investigate theoretical properties on the lower bounds of
estimation tasks for heavy-tailed high-dimensional time series data, including
estimation of the high-dimensional VAR models and the autocovariance matrices.
The lower bound analysis shows that the rates of convergence obtained in
Section 3 are optimal in the minimax sense: there exists a distributional
setting for the time series process for which the upper bounds obtained cannot
be improved without further distributional assumptions.

###  4.1 Lower bounds of VAR estimation

Similarly to the upper bound analysis, no assumption on the data generating
mechanism is imposed on the time series data in the lower bound analysis. We
consider the VAR($d$) model and denote the distribution of
$\\{\mathbf{y}_{t}\\}_{t=1}^{T}$ as $\mathbb{P}$. By the Yule–Walker equation,
the true value of the VAR model coefficient matrix is defined as
$\mathbf{A}^{*}(\mathbb{P})=\mathbf{\Sigma}_{1}(\mathbb{P})\mathbf{\Sigma}_{0}^{-1}(\mathbb{P})$,
where
$\mathbf{\Sigma}_{0}(\mathbb{P})=\mathbb{E}_{\mathbb{P}}[\mathbf{x}_{t}\mathbf{x}_{t}^{\top}]$
and
$\mathbf{\Sigma}_{1}(\mathbb{P})=\mathbb{E}_{\mathbb{P}}[\mathbf{y}_{t}\mathbf{x}_{t}^{\top}]$.

The first case we consider is the sparse VAR model, where all row vectors of
$\mathbf{A}^{*}$ are strictly sparse. For any $M>0$, $\epsilon\in(0,1]$ and
$r\in(0,1)$, let $\mathcal{P}_{\text{E}}(M,\epsilon,r)$ denote the class of
all joint distributions for the $\alpha$-mixing stochastic process
$\mathbf{v}_{t}=(v_{1t},\dots,v_{pt})^{\top}$ with the element-wise moment
condition, such that they satisfy $\max_{1\leq i\leq
p}\mathbb{E}|v_{it}|^{2+2\epsilon}=M$ and
$\alpha(\\{\mathbf{v}_{t}\\}_{t=-\infty}^{s},\\{\mathbf{v}_{t}\\}_{t=s+\ell}^{\infty})=O(r^{\ell})$
for any integer $s$ and lag order $\ell>0$. The moment parameter $M$ is
imposed on the element of $\mathbf{v}_{t}$, which is consistent with
Assumption 2.

For $\mathbf{A}^{*}$ belonging to the sparse set $\mathbb{B}_{0}(s_{0})$
defined in Section 3.1 and the process $\mathbf{y}_{t}$ following the
distribution $\mathbb{P}\in\mathcal{P}_{\text{E}}(M,\epsilon,r)$, we have the
minimax lower bound for the VAR estimation.

######  Theorem 5 (Sparse VAR lower bound).

For any $M>0$, $\epsilon\in(0,1]$ and $r\in(0,1)$, suppose that
$T/\log_{r^{-2}}T\geq 2s_{0}$ and the joint distribution $\mathbb{P}$ of
$\mathbf{y}_{t}$ belongs to $\mathcal{P}_{\textup{E}}(M,\epsilon,r)$. Then,
for any estimator
$\widehat{\mathbf{A}}:=\widehat{\mathbf{A}}(\\{\mathbf{y}_{t}\\}_{t=1}^{T})$
which depends on the observations from $\mathbf{y}_{1}$ to $\mathbf{y}_{T}$,

$\inf_{\widehat{\mathbf{A}}}\sup_{\begin{subarray}{c}\mathbb{P}\in\mathcal{P}_{\textup{E}}(M,\epsilon,r),\\\
\mathbf{A}^{*}(\mathbb{P})\in\mathbb{B}_{0}(s_{0})\end{subarray}}\mathbb{E}\left[\|\widehat{\mathbf{A}}^{\top}-\mathbf{A}^{*}(\mathbb{P})^{\top}\|_{2,\infty}\right]\gtrsim\sqrt{s_{0}}\|\mathbf{\Sigma}_{0}^{-1}(\mathbb{P})\|_{1,\infty}\left[\frac{M^{1/\epsilon}\log(T)}{T}\right]^{\frac{\epsilon}{1+\epsilon}}.$
(44)

The minimax lower bound in this theorem matches the $\ell_{2,\infty}$ upper
bound in Theorem 1 with $q=0$ up to a logarithm factor that is negligible
compared with $T$. Hence, for the strictly sparse VAR model, the constrained
Yule–Walker estimator $\widehat{\mathbf{A}}_{\text{S}}(\lambda,\tau)$
introduced in Section 3.1 is nearly rate-optimal when the regularization
parameter $\lambda$ and truncation parameter $\tau$ are properly chosen. In
addition, the effective sample size in the lower bound is $T/\log(T)$, where
the $\log(T)$ factor can be referred to as the price paid for the serial
dependency in the data.

Moreover, it is clear that the banded VAR model with the bandwidth $k$ is a
special case of the strictly sparse VAR model with the sparsity level
$s_{0}=2k+1$. When $k$ is fixed, Theorem 5 directly implies that the minimax
lower bound of the banded VAR model is of rate
$\|\mathbf{\Omega}^{-1}\|_{1,\infty}[M^{1/\epsilon}\log(T)/T]^{\epsilon/(1+\epsilon)}$
in the operator norm, matching the upper bound of the banded VAR in Theorem 3.

The second case considered is the reduced-rank VAR model where
$\mathbf{A}^{*}$ is of low rank. For any $M>0$, $\epsilon\in(0,1]$ and
$r\in(0,1)$, let $\mathcal{P}_{\text{V}}(M,\epsilon,r)$ denote the class of
all joint distributions for the $\alpha$-mixing stochastic process
$\mathbf{v}_{t}=(v_{1t},\dots,v_{pt})^{\top}$ with the vector-wise moment
condition, such that
$\sup_{\mathbf{u}\in\mathbb{R}^{p},\|\mathbf{u}\|_{2}=1}\mathbb{E}|\mathbf{v}_{t}^{\top}\mathbf{u}|^{2+2\epsilon}=M$
and
$\alpha(\\{\mathbf{v}_{t}\\}_{-\infty}^{s},\\{\mathbf{v}_{t}\\}_{s+\ell}^{\infty})=O(r^{\ell})$
for any integer $s$ and lag order $\ell>0$. The moment condition defined on
the whole vector is similar to Assumption 3 for the reduced-rank VAR model,
but the technical term $\delta$ is omitted in the lower bound analysis.

For $\mathbf{A}^{*}$ belonging to the low-rank set
$\widetilde{\mathbb{B}}_{0}(r_{0})$ defined in Section 3.2 and the process
$\mathbf{y}_{t}$ following the distribution
$\mathbb{P}\in\mathcal{P}_{\text{V}}(M,\epsilon,r)$, the lower bound in terms
of Frobenius norm is derived.

######  Theorem 6 (Reduced-rank VAR lower bound).

For any $M>0$, $\epsilon\in(0,1]$ and $r\in(0,1)$, suppose that $p\geq 20$,
$T>p/128$, and the joint distribution $\mathbb{P}$ of $\mathbf{y}_{t}$ belongs
to $\mathcal{P}_{\textup{V}}(M,\epsilon,r)$. Then, for any estimator
$\widehat{\mathbf{A}}:=\widehat{\mathbf{A}}(\\{\mathbf{y}_{t}\\}_{t=1}^{T})$
which depends on the observations from $\mathbf{y}_{1}$ to $\mathbf{y}_{T}$,

$\inf_{\widehat{\mathbf{A}}}\sup_{\begin{subarray}{c}\mathbb{P}\in\mathcal{P}_{\textup{V}}(M,\epsilon,r),\\\
\mathbf{A}^{*}(\mathbb{P})\in\widetilde{\mathbb{B}}_{0}(r_{0})\end{subarray}}\mathbb{E}\left[\|\widehat{\mathbf{A}}^{\top}-\mathbf{A}^{*}(\mathbb{P})^{\top}\|_{\textup{F}}\right]\gtrsim\sqrt{r_{0}}\|\mathbf{\Sigma}_{0}^{-1}(\mathbb{P})\|_{\textup{op}}\left(\frac{pM^{1/\epsilon}}{T}\right)^{\frac{\epsilon}{1+\epsilon}}.$
(45)

This theorem presents the minimax lower bound of the reduced-rank VAR model
with the exact low-rankness in terms of Frobenius norm, which matches the
upper bound in Theorem 2 up to a logarithm factor.

######  Remark 9.

The minimax lower bounds in Theorems 5 and 6 are developed for the VAR models
with an exact low-dimensional structure, such as the strict row-wise sparsity
and exact low-rankness. The lower bounds for the high-dimensional regression
and covariance estimation with the weakly sparse coefficients have been
investigated by Raskutti et al., (2011), Cai and Zhou, (2012), and others
under the $\ell_{q}$ ball constraints. However, the theoretical techniques in
their proofs rely heavily on the Gaussian distributional assumption, and hence
cannot be applied to the heavy-tailed setting. The minimax lower bounds for
the VAR estimation under the weak sparsity and approximate low-rankness are
left for future research.

###  4.2 Lower bounds of autocovariance estimation

For any process $\\{\mathbf{y}_{t}\\}$ following distribution
$\mathbb{P}\in\mathcal{P}_{\text{E}}(M,\epsilon,r)$, let
$\mathbf{x}_{t}=(\mathbf{y}_{t-1}^{\top},\dots,\mathbf{y}_{t-d}^{\top})^{\top}$,
and a by-product of our lower bound analysis is the minimax lower bounds for
the autocovariance estimation. First, we develop the minimax lower bounds in
the $\ell_{\infty}$ norm for estimating the autocovariance matrices
$\mathbf{\Sigma}_{0}(\mathbb{P})=\mathbb{E}_{\mathbb{P}}[\mathbf{x}_{t}\mathbf{x}_{t}^{\top}]$
and
$\mathbf{\Sigma}_{1}(\mathbb{P})=\mathbb{E}_{\mathbb{P}}[\mathbf{y}_{t}\mathbf{x}_{t}^{\top}]$.

######  Proposition 7.

For any $M>0$, $\epsilon\in(0,1]$ and $r\in[0,1)$, suppose that the joint
distribution $\mathbb{P}$ of $\mathbf{y}_{t}$ belongs to
$\mathcal{P}_{\textup{E}}(M,\epsilon,r)$. Then, for any covariance estimator
$\widehat{\mathbf{\Sigma}}_{k}:=\widehat{\mathbf{\Sigma}}_{k}(\\{\mathbf{y}_{t}\\}_{t=1}^{T})$
which depends on the observations from $\mathbf{y}_{1}$ to $\mathbf{y}_{T}$,

$\inf_{\widehat{\mathbf{\Sigma}}_{k}}\sup_{\mathbb{P}\in\mathcal{P}_{\textup{E}}(M,\epsilon,r)}\mathbb{E}\left[\|\widehat{\mathbf{\Sigma}}_{k}-\mathbf{\Sigma}_{k}(\mathbb{P})\|_{\infty}\right]\gtrsim\left[\frac{M^{1/\epsilon}\log(T)}{T}\right]^{\frac{\epsilon}{1+\epsilon}},\quad
k=0,1.$ (46)

The minimax lower bound matches the upper bounds in Proposition 4 up to a
logarithmic factor $[\log(p^{2}d)\log(T)]^{\epsilon/(1+\epsilon)}$, indicating
the nearly minimax optimality of the proposed shrinkage estimator
$\widetilde{\mathbf{\Sigma}}_{k}^{\text{E}}(\tau)$ with $\tau$ chosen
properly. Compared with the upper and lower bounds of the covariance matrix
estimation for the i.i.d. data in Avella-Medina et al., (2018) and Devroye et
al., (2016), our minimax lower bound involves a factor
$\log(T)^{\epsilon/(1+\epsilon)}$, which is due to the weak serial dependency
under the geometrically decayed $\alpha$-mixing condition. It is also
noteworthy that the $\log(T)$ factor can also be found in the upper bound
analysis in Zhang, (2021) under the functional dependence measures (Wu,,
2005).

In addition, this minimax lower bound can easily be extended to the linear
transformations of autocovariance matrices
$\mathbf{\Omega}=\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\mathbf{\Sigma}_{0})\mathbf{C}$
and $\boldsymbol{\omega}=\mathbf{C}^{\top}\textup{vec}(\mathbf{\Sigma}_{1})$,
where $\mathbf{C}$ consists of coordinate vector columns. In other words,
Proposition 7 implies that the upper bounds in Proposition 6 are rate-optimal
up to a logarithm factor.

Moreover, the minimax lower bound in the operator norm is also established.

######  Proposition 8.

For any $M>0$, $\epsilon\in(0,1]$ and $r\in[0,1)$, suppose that $p\geq 20$,
$T\geq p/128$, and the joint distribution $\mathbb{P}$ of $\mathbf{y}_{t}$
belongs to $\mathcal{P}_{\textup{V}}(M,\epsilon,r)$. Then, for any covariance
estimator
$\widehat{\mathbf{\Sigma}}_{k}:=\widehat{\mathbf{\Sigma}}_{k}(\\{\mathbf{y}_{t}\\}_{t=1}^{T})$
which depends on the observations from $\mathbf{y}_{1}$ to $\mathbf{y}_{T}$,

$\inf_{\widehat{\mathbf{\Sigma}}_{k}}\sup_{\mathbb{P}\in\mathcal{P}_{\textup{V}}(M,\epsilon,r)}\mathbb{E}\left[\|\widehat{\mathbf{\Sigma}}_{k}-\mathbf{\Sigma}_{k}(\mathbb{P})\|_{\textup{op}}\right]\gtrsim\left(\frac{pM^{1/\epsilon}}{T}\right)^{\frac{\epsilon}{1+\epsilon}},\quad
k=0,1.$ (47)

This lower bound matches the upper bound in Proposition 5 up to a logarithm
factor $\log(T)^{\epsilon/(1+\epsilon)}$, indicating that the vector-wise
truncated estimators $\widetilde{\mathbf{\Sigma}}^{\text{V}}_{k}(\tau)$
introduced in Section 3.2 are nearly rate-optimal when $\tau$ is chosen
properly.

######  Remark 10.

To the best of our knowledge, for covariance or autocovaraince estimation
problems, the minimax lower bound in Proposition 8 is the first operator norm
lower bound result with the phase transition phenomenon. This lower bound is
sharper than that in Fan et al., (2021) by constructing a special multivariate
discrete distribution with the $(2+2\epsilon)$-th moment condition; see
Appendix B for details.

##  5 Algorithm and Implementation

###  5.1 Linearized ADMM algorithm

The constrained Yule–Walker estimators in (6), (12) and (19) are convex
optimization problems as each of them consists of a convex objective function
and a convex constraint function. For the constrained problem in (6), define
the constraint set
$\mathbb{C}_{\mathcal{R}^{*}}(\lambda)=\\{\mathbf{M}\in\mathbb{R}^{p\times
pd}:\mathcal{R}^{*}(\mathbf{M})\leq\lambda\\}$ and the optimization problem
can be rewritten as

$\min_{\mathbf{A},\mathbf{D}}\mathcal{R}(\mathbf{A}),\quad\text{subject to
}\mathbf{D}=\mathbf{A}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\widetilde{\Sigma}}_{1}~{}\text{and}~{}\mathbf{D}\in\mathbb{C}_{\mathcal{R}^{*}}(\lambda).$
(48)

The augmented Lagrangian form is given as

$\mathcal{L}_{\rho}(\mathbf{A},\mathbf{D};\mathbf{W})=\mathcal{R}(\mathbf{A})-\langle\mathbf{W},\mathbf{A}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{D}\rangle+\frac{\rho}{2}\|\mathbf{A}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{D}\|_{\textup{F}}^{2},$
(49)

subject to $\mathbf{D}\in\mathbb{C}(\lambda)$, where $\mathbf{W}$ is the
Lagrangian multiplier and $\rho$ is the regularization parameter. The
augmented Lagrangian form can be solved by the alternating direction method of
multipliers (Boyd et al.,, 2011, ADMM) with the iterative updates of
$(\mathbf{A},\mathbf{D},\mathbf{W})$:

$\begin{split}\mathbf{A}^{(j+1)}&=\underset{\mathbf{A}}{\operatorname*{arg\,min}}~{}\mathcal{R}(\mathbf{A})+\frac{\rho}{2}\|\mathbf{A}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{D}^{(j)}-\mathbf{W}^{(j)}/\rho\|_{\textup{F}}^{2},\\\
\mathbf{D}^{(j+1)}&=\underset{\mathbf{D}\in\mathbb{C}_{\mathcal{R}^{*}}(\lambda)}{\operatorname*{arg\,min}}\|\mathbf{D}+\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{A}^{(j+1)}\mathbf{\widetilde{\Sigma}}_{0}+\mathbf{W}^{(j)}/\rho\|_{\textup{F}}^{2},\\\
\mathbf{W}^{(j+1)}&=\mathbf{W}^{(j)}+\rho(\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{A}^{(j+1)}\mathbf{\widetilde{\Sigma}}_{0}+\mathbf{D}^{(j+1)}).\end{split}$
(50)

The efficiency of the ADMM algorithm depends largely on the complexity of the
resulting subproblems. Note that $\mathbf{A}$-update in the ADMM algorithm is
a regularized least squares problem of $\mathbf{A}$, which might not have an
explicit solution and can be computationally expensive.

To alleviate the computational burden in $\mathbf{A}$-update, we use an
inexact proximal method to minimize

$\begin{split}&\widetilde{\mathcal{L}}_{\rho}^{(j)}(\mathbf{A},\mathbf{D}^{(j)};\mathbf{W}^{(j)})\\\
=&\mathcal{R}(\mathbf{A})+\frac{\rho\mu}{2}\left\|\mathbf{A}-\left[\mathbf{A}^{(j)}-\frac{2}{\mu}(\mathbf{A}^{(j)}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\widetilde{\Sigma}}_{1}-\mathbf{D}^{(j)}-\mathbf{W}^{(j)}/\rho)\mathbf{\widetilde{\Sigma}}_{0}\right]\right\|_{\textup{F}}^{2}\end{split}$
(51)

where $\mu$ is another prespecified parameter. As in Wang and Yuan, (2012),
$\widetilde{\mathcal{L}}^{(j)}_{\rho}(\mathbf{A},\mathbf{D}^{(j)};\mathbf{W}^{(j)})$
can be viewed as a linearized approximation of
$\mathcal{L}_{\rho}(\mathbf{A},\mathbf{D}^{(j)};\mathbf{W}^{(j)})$.

For the separated constrained Yule–Walker estimation in (12), the linearized
ADMM algorithm can be developed in similar fashion. For the $\ell_{1}$ norm
and nuclear norm regularized estimators, in $\mathbf{D}$-update and the
proximal variant of $\mathbf{A}$-update, the soft-thresholding operator and
truncation operator can be applied to the elements and singular values,
respectively. The proposed ADMM algorithm can be illustrated as a particular
application of the inexact Bregman ADMM algorithm (Wang and Banerjee,, 2014).
To guarantee the global convergence, $\mu/2$ needs to be greater than the
largest eigenvalue of $\mathbf{\widetilde{\Sigma}}_{0}^{2}$.

In summary, the ADMM algorithm has explicit updates and can be solved
efficiently. The details of the linearized ADMM algorithm for sparse and
reduced-rank VAR models are presented in Appendix C.

###  5.2 Linear and semidefinite programmings

Some specific cases of the constrained Yule–Walker estimator in (6) can be
formulated as linear programmings (LP) or semidefinite programs (SDP), and can
thus be efficiently solved by any of the standard LP or SDP solvers.

As discussed in Example 2, for sparse VAR models, a separated constrained
minimization framework in (12) is proposed with
$\mathcal{R}(\cdot)=\|\cdot\|_{1}$ and
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\infty}$. Following the lead of Dantzig
selector (Candès and Tao,, 2007), each of the sub-problems can be recast to an
LP:

$\begin{split}\text{minimize}&~{}\mathbf{1}_{pd}^{\top}\mathbf{u}\\\
\text{subject
to}&~{}-\mathbf{u}\leq\mathbf{a}_{i}\leq\mathbf{u}~{}\text{and}~{}-\lambda\cdot\mathbf{1}_{pd}\leq\widetilde{\boldsymbol{\sigma}}_{1i}-\widetilde{\mathbf{\Sigma}}_{0}\mathbf{a}_{i}\leq\lambda\cdot\mathbf{1}_{pd},\end{split}$
(52)

where the optimization variables are $\mathbf{u}\in\mathbb{R}^{pd}$ and
$\mathbf{a}_{i}\in\mathbb{R}^{pd}$. As discussed in Candès and Tao, (2007),
there is a large class of efficient algorithms for solving such problems.

Similarly, the constrained Yule–Walker estimator in (6), with
$\mathcal{R}(\cdot)=\|\cdot\|_{\textup{nuc}}$ and
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\textup{op}}$, can be used to estimate the
reduced-rank VAR models. As introduced by Candès and Plan, (2011), the nuclear
norm minimization problem can be recast to a SDP with a linear matrix
inequality constraint:

$\begin{split}\text{minimize}&~{}[\text{tr}(\mathbf{W}_{1})+\text{tr}(\mathbf{W}_{2})]/2\\\
\text{subject to}&~{}\begin{bmatrix}\mathbf{W}_{1}&\mathbf{A}\\\
\mathbf{A}^{\top}&\mathbf{W}_{2}\\\ \end{bmatrix}\succeq
0,~{}\text{and}~{}\begin{bmatrix}\lambda\mathbf{I}_{p}&\widetilde{\mathbf{\Sigma}}_{1}-\mathbf{A}\widetilde{\mathbf{\Sigma}}_{0}\\\
\widetilde{\mathbf{\Sigma}}_{1}^{\top}-\widetilde{\mathbf{\Sigma}}_{0}\mathbf{A}^{\top}&\lambda\mathbf{I}_{pd}\end{bmatrix}\succeq
0,\end{split}$ (53)

with optimization variables $\mathbf{W}_{1}\in\mathbb{R}^{p\times p}$,
$\mathbf{W}_{2}\in\mathbb{R}^{pd\times pd}$, and
$\mathbf{A}\in\mathbb{R}^{p\times pd}$. Interior point methods or first-order
methods, such as splitting cone solver, can be applied to efficiently solve
large-scale SDP.

For linear-restricted VAR models, the $\ell_{\infty}$ optimization problem of
the constrained Yule–Walker estimator in (19) can also be formulated as a LP:

$\begin{split}\text{minimize}&~{}t\\\ \text{subject
to}&~{}-t\cdot\mathbf{1}_{r}\leq\boldsymbol{\theta}\leq
t\cdot\mathbf{1}_{r},~{}\text{and}~{}-\lambda\cdot\mathbf{1}_{r}\leq\mathbf{\widetilde{\Omega}}\boldsymbol{\theta}-\boldsymbol{\widetilde{\omega}}\leq\lambda\cdot\mathbf{1}_{r},\end{split}$
(54)

with optimization variables $t\in\mathbb{R}$ and
$\boldsymbol{\theta}\in\mathbb{R}^{r}$, and can be solved by standard LP
solvers.

As the data truncation and robust autocovariance estimator calculation are
computationally cheap, the operation time of ADMM, LP and SDP algorithms
mainly depends on the dimension $p$ and AR order $d$. As shown in the
simulation experiments in Section 6, the proposed ADMM algorithm is much more
efficient than the LP and SDP solvers, especially when $p$ is large.

###  5.3 Tuning parameter selection

For the constrained Yule–Walker estimators, both the robustification
parameter(s) $\tau$ (or $\tau_{1}$ and $\tau_{2}$) and constraint parameter
$\lambda$ need to adapt properly to the dimension $p$, sample size $T$ and
moment condition bound $M_{2+2\epsilon}$ to achieve optimal trade-off between
the truncation bias and tail robustness. However, based on the conditions and
convergence rates in Section 3, the optimal value of both parameters rely on
some unknown parameters $\epsilon$ and $M_{2+2\epsilon}$. In the literature on
robust estimation of i.i.d. data, cross-validation is an intuitive data-driven
method to select robustification and regularization parameters. To accommodate
the intrinsically ordered nature of time series data, we use a rolling
forecasting validation, one of the standard approaches to tuning parameter
selection for time series data.

To simultaneously select robustification and constraint parameters, a two-
dimensional or multidimensional grid of tuning parameters is constructed. The
search of the constraint parameter $\lambda$ starts from a predetermined
$\lambda_{\max}$ and decreases in log-linear increments. For the
robustification parameter $\tau$, grid search on the data-driven interval
$[\tau_{\min},\tau_{\max}]$ is suggested, where $\tau_{\min}$ and
$\tau_{\max}$ are empirical quantiles of the data. For instance, to tune
$\tau$ for the element truncation estimators, we may set $\tau_{\min}$ and
$\tau_{\max}$ to be the median and maximum of $|y_{it}|$, for all $1\leq i\leq
p$ and $1\leq t\leq T$, respectively. If at least one of $p$ and $T$ is very
large, we may calculate the quantiles of some randomly sampled $|y_{it}|$.

For the dataset with time series observations from $t=1$ to $T$, we determine
a validation period $[T_{\text{val}},T]$ whose length $T-T_{\text{val}}$ is
supposed to be slightly smaller than $T$. Given the constraint and
robustification parameters $(\lambda,\tau)$, at any time point
$t\in[T_{\text{val}},T-1]$, we use all the historical data to calculate the
constrained Yule–Walker estimates $\mathbf{\widehat{A}}(\lambda,\tau)$ and the
corresponding one-step-ahead forecast
$\mathbf{\widehat{A}}(\lambda,\tau)\mathbf{x}_{t}$. The performance evaluation
is based on the mean squared forecast error (MSFE)

$\text{MSFE}(\lambda,\tau;T_{\text{val}},T)=\frac{1}{T-T_{\text{val}}}\sum_{t=T_{\text{val}}+1}^{T}\|\mathbf{y}_{t}-\mathbf{\widehat{A}}(\lambda,\tau)\mathbf{x}_{t}\|_{2}^{2},$
(55)

and the tuning parameters are selected by minimizing the
$\text{MSFE}(\lambda,\tau;T_{\text{val}},T)$ over the grid points.

##  6 Simulation Study

###  6.1 VAR estimation

In this subsection, we compare the performance of the robust procedure and two
standard estimation procedures for sparse, reduced-rank, and banded VAR models
in finite samples. For each model, we consider three innovation settings:
standardized $t_{2.1}$-distributed innovations, standardized log-normal
innovations, and standard Gaussian innovations. They represent heavy-tailed
symmetric distributions with finite second moments, heavy-tailed asymmetric
distributions with finite fourth moments, and light-tailed symmetric
distribution, respectively. The data are drawn from the VAR(1) process
$\mathbf{y}_{t}=\mathbf{A}^{*}\mathbf{y}_{t-1}+\boldsymbol{\varepsilon}_{t}$
with the following structures of $\mathbf{A}^{*}$.

DGP-1: Sparse VAR with $\mathbf{A}^{*}$ specified as $\mathbf{A}^{*}_{ij}=0.5$
if $i=j$, $\mathbf{A}^{*}_{ij}=0.4$ if $i-j=1$, $\mathbf{A}^{*}_{ij}=-0.4$ if
$i-j=-1$, and $\mathbf{A}^{*}_{ij}=0$ otherwise. We consider $p\in\\{30,50\\}$
and $T\in\\{100,150,200,250,300\\}$.

DGP-2: Reduced-rank VAR with $\mathbf{A}^{*}$ having a low-rank singular value
decomposition $\mathbf{A}^{*}=\mathbf{U}\mathbf{D}\mathbf{V}^{\top}$.
Specifically, we consider a fixed $\mathbf{D}=\text{diag}(1.5,0.8)$, and
generate random orthonormal $\mathbf{U}\in\mathbb{R}^{p\times 2}$ and
$\mathbf{V}\in\mathbb{R}^{p\times 2}$ from the first two leading singular
vectors of a random $p$-by-$p$ Gaussian matrix. We consider $p\in\\{30,50\\}$
and $T\in\\{300,350,400,450,500\\}$.

DGP-3: Banded VAR with a banded $\mathbf{A}^{*}$ having bandwidth $k=2$.
Specifically, $\mathbf{A}^{*}_{ij}=0.5$ if $i=j$; $\mathbf{A}^{*}_{ij}=0.3$ if
$i-j=1$; $\mathbf{A}^{*}_{ij}=-0.3$ if $i-j=-1$; $\mathbf{A}^{*}_{ij}=0.2$ if
$|i-j|=2$; and $\mathbf{A}^{*}_{ij}=0$ for all $|i-j|>2$. We consider
$p\in\\{50,100\\}$ and $T\in\\{100,150,200,250,300\\}$.

We first specify the robust estimation procedure, denoted as ROB, for these
three DGPs. For DGP-1, as discussed in Example 2, we apply the split
estimation method in (12) with $\mathcal{R}(\cdot)=\|\cdot\|_{1}$,
$\mathcal{R}(\cdot)=\|\cdot\|_{\infty}$, and the element-wise truncation
autocovariance estimators $\widetilde{\mathbf{\Sigma}}_{0}^{\text{E}}(\tau)$
and $\widetilde{\mathbf{\Sigma}}_{1}^{\text{E}}(\tau)$ specified in Section
3.1. For DGP-2, as discussed in Example 1, we apply the estimation method in
(6) with $\mathcal{R}(\cdot)=\|\cdot\|_{\textup{nuc}}$,
$\mathcal{R}^{*}(\cdot)=\|\cdot\|_{\textup{op}}$, and the vector-wise
truncation autocovariance estimators
$\widetilde{\mathbf{\Sigma}}_{0}^{\text{V}}(\tau)$ and
$\widetilde{\mathbf{\Sigma}}_{1}^{\text{V}}(\tau)$ specified in Section 3.2.
For DGP-3, as discussed in Example 3, we apply the estimation method in (19)
with the element-wise truncation autocovariance estimators
$\widetilde{\mathbf{\Omega}}^{\text{E}}(\tau)$ and
$\widetilde{\boldsymbol{\omega}}^{\text{E}}(\tau)$ in Section 3.3. The tuning
parameters $\lambda$ and $\tau$ are selected simultaneously by the MSFE in
Section 5.3.

Alternatively, if the possibility of heavy-tailed distribution is ignored, we
consider the sample autocovariance estimators

$\widehat{\mathbf{\Sigma}}_{0}=\frac{1}{T}\sum_{t=1}^{T}\mathbf{y}_{t}\mathbf{y}_{t}^{\top},~{}~{}\widehat{\mathbf{\Sigma}}_{1}=\frac{1}{T-1}\sum_{t=2}^{T}\mathbf{y}_{t}\mathbf{y}_{t-1}^{\top},~{}~{}\widehat{\mathbf{\Omega}}=\mathbf{C}^{\top}(\mathbf{I}_{p}\otimes\widetilde{\mathbf{\Sigma}}_{0})\mathbf{C},~{}~{}\text{and}~{}~{}\mathbf{C}^{\top}\textup{vec}(\widehat{\mathbf{\Sigma}}_{1}^{\top}),$
(56)

and plug them in the Yule–Walker estimators. We denote this estimation
procedure as the regular Yule–Walker (RYW) method. In the high-dimensional VAR
literature, another popular estimation procedure is the regularized least
squares (RLS) method. For DGP-1 and 2, we apply the regularized method

$\widehat{\mathbf{A}}_{\text{RLS}}=\operatorname*{arg\,min}_{\mathbf{A}}\frac{1}{T}\sum_{t=1}^{T}\|\mathbf{y}_{t}-\mathbf{A}\mathbf{y}_{t-1}\|_{2}^{2}+\lambda\mathcal{R}(\mathbf{A})$
(57)

with $\mathcal{R}(\cdot)=\|\cdot\|_{1}$ and $\|\cdot\|_{\textup{nuc}}$,
respectively. For DGP-3, we consider the least squares (LS) estimator

$\widehat{\boldsymbol{\theta}}_{\text{LS}}=\operatorname*{arg\,min}_{\boldsymbol{\theta}}\frac{1}{T}\|\mathbf{y}_{t}-(\mathbf{I}_{p}\otimes\mathbf{y}_{t-1}^{\top})\mathbf{C}\boldsymbol{\theta}\|_{2}^{2}.$
(58)

To evaluate the estimation performance and to confirm the theoretical results
in Theorems 1-3, for $\widehat{\mathbf{A}}-\mathbf{A}^{*}$, we calculate the
$\ell_{2,\infty}$ norm for sparse VAR models, the Frobenius norm for reduced-
rank VAR models, and the operator norm and $\ell_{\infty}$ norm for banded VAR
models. For the three DGPs, the estimation errors are calculated by averaging
500 replications for each setting and presented in Figures 1, 2 and 3,
respectively.

Figure 1: Estimation errors of robust (ROB), regularized Yule–Walker (RYW),
and regularized least squares (RLS) methods in the $\ell_{2,\infty}$ norm v.s.
sample size $T$ for various dimension and innovation distribution settings in
DGP-1. Figure 2: Estimation errors of robust (ROB), regularized Yule–Walker
(RYW), and regularized least squares (RLS) methods in the Frobenius norm v.s.
sample size $T$ for various dimension and innovation distribution settings in
DGP-2. Figure 3: Estimation errors of robust (ROB), regularized Yule–Walker
(RYW), and least squares (LS) methods in the operator norm v.s. sample size
$T$ for various dimension and innovation distribution settings in DGP-3.

For the three types of VAR models considered, we can observe from Figures 1-3
that the proposed robust estimation method yields much small statistical
errors than both standard estimation methods under the heavy-tailed
innovations, i.e. the standardized $t_{2.1}$ distributed innovation and log-
normal innovation. In particular, under the $t_{2.1}$ innovation setting, the
estimation performance of both standard estimation methods is not stable.
Under the setting of Gaussian innovation, the proposed robust estimator
produces almost the same or even slightly smaller estimation errors than the
standard methods, indicating that it does not hurt to use the robust procedure
under the light-tailed innovation setting. In sum, the simulation results
generally confirm the theoretical convergence rates and demonstrate the
robustness of the proposed estimators against heavy-tailed distributions.

We also compare the performance of two computational algorithms discussed in
Section 5, using the simulated data with a fixed $T=300$ and varying dimension
$p\in\\{25,50,75,100,125\\}$ from the Gaussian settings of DGP-1 and 2,
respectively. We obtain almost identical estimates from the two algorithms
under the same choices of $\lambda$ and $\tau$. Regarding computational time,
we record the averaged time consumed over 100 replications for each algorithm.
Figure 4 demonstrates that the ADMM is much more efficient than the LP and SDP
solvers, especially when the dimension $p$ is large.

Figure 4: Computational time (in seconds) of ADMM and LP/SDP algorithms v.s.
dimension $p$ in DGP-1 and DGP-2.

###  6.2 Autocovariance matrix estimation

We also conduct a simulation experiment to compare the truncation-based
autocovariance matrix estimator with the standard sample autocovariance
estimator. Two data generating processes, DGP-1 and DGP-2 in Subsection 6.1,
are employed. Similarly to the previous experiment, we also consider three
innovation settings for each DGP. The element and vector truncation estimators
are considered for these two processes. Based on 1000 replications, the
average estimation errors of autocovariance matrices in terms of
$\ell_{\infty}$ and operator norms are presented in Figures 5 and 6 for the
two processes, respectively.

Figure 5: Estimation errors of
$\widetilde{\mathbf{\Sigma}}_{0}^{\textup{E}}(\tau)$ and
$\widetilde{\mathbf{\Sigma}}_{1}^{\textup{E}}(\tau)$ in the $\ell_{\infty}$
norm v.s. sample size $T$ for different dimension and innovation distribution
settings in DGP-1. Figure 6: Estimation errors of
$\widetilde{\mathbf{\Sigma}}_{0}^{\textup{E}}(\tau)$ and
$\widetilde{\mathbf{\Sigma}}_{1}^{\textup{E}}(\tau)$ in the operator norm v.s.
sample size $T$ for different dimension and innovation distribution settings
in DGP-2.

For these two DGPs, the truncation-based autocovariance matrix estimators
yield much smaller estimation errors than the sample autocovariance estimator
under the heavy-tailed innovation settings. Similarly to the previous
experiment, under the $t_{2.1}$ innovation setting, the performance of the
sample autocovariance estimator is not stable, but the smoothly decreasing
pattern of the proposed robust estimator can be observed. Under the Gaussian
innovation setting, the performance of both estimators are nearly the same.
The numerical results in this experiment confirm the rates of convergence in
Propositions 4 and 5 and verify the robustness of the proposed autocovariance
estimators.

##  7 Real Data Example

In this section, we apply the proposed robust estimation procedure to a real
data set consisting of 40 quarterly macroeconomic variables of the United
States from the third quarter of 1959 to fourth quarter of 2007, with 194
observations for each variable. Following Stock and Watson, (2009), all
variables are seasonally adjusted except for financial variables, transformed
by differencing or log differencing to be stationary, and standardized to zero
mean and unit standard deviation. These macroeconomic variables capture many
aspects of the U.S. economy, including GDP, National Association of Purchasing
Manager indices, industrial production, pricing, employment, credit and
interest rate. Both factor model and VAR model have been applied to these
series in empirical econometric analyses for structural analysis and
forecasting; see Stock and Watson, (2009), Koop, (2013) and Wang et al., 2021a
. The kurtosis of the standardized variables are plotted in Figure 7, with the
blue and red dashed lines indicating the kurtosis of normal distribution and
$t_{6}$ distribution, respectively. Among 40 variables, 36 and 11 variables
have larger kurtosis than standard normal distribution and $t_{6}$
distribution, providing some evidence of the existence of heavy-tailed
distributions in this macroeconomic dataset.

Figure 7: Histogram of kurtosis of 40 macroeconomic variables. Blue and red
dashed line mark the theoretical kurtosis of standard normal distribution and
$t_{6}$ distribution.

Following Koop, (2013), we apply a VAR(4) model to these macroeconomic time
series. If no structural assumption is imposed on the parameter matrix, we can
directly use the ordinary least squares estimator by minimizing the least
squares loss function, namely

$\mathbf{\widehat{A}}_{\text{OLS}}=\operatorname*{arg\,min}_{\mathbf{A}\in\mathbb{R}^{40\times
160}}\frac{1}{T}\sum_{t=1}^{T}\|\mathbf{y}_{t}-\mathbf{A}\mathbf{x}_{t}\|_{2}^{2},$
(59)

where
$\mathbf{x}_{t}=(\mathbf{y}_{t-1}^{\top},\cdots,\mathbf{y}_{t-4}^{\top})^{\top}$.
If the low-rank structure on $\mathbf{A}$ is considered and the rank $r$ is
pre-specified, the rank-constrained least squares estimator is formulated as

$\mathbf{\widehat{A}}_{\text{RRR}}=\operatorname*{arg\,min}_{\mathbf{A}\in\mathbb{R}^{40\times
160},\text{rank}(\mathbf{A})\leq
r}\frac{1}{T}\sum_{t=1}^{T}\|\mathbf{y}_{t}-\mathbf{A}\mathbf{x}_{t}\|_{2}^{2}.$
(60)

As shown in Section 6, if the low-rank structure is not directly imposed, we
can use the nuclear norm regularized least squares (NN-RLS), nuclear norm
regularized Yule–Walker (NN-RYW), and the proposed nuclear norm robust (NN-
ROB) methods. For the sparse VAR model, the $\ell_{1}$ regularized least
squares (L1-RLS), $\ell_{1}$ regularized Yule–Walker (L1-RYW), and $\ell_{1}$
robust (L1-ROB) methods can be applied as discussed in Section 6.

The performance of these eight methods are compared via out-of-sample
forecasting errors. From the first quarter of 1993 ($t=135$) to the fourth
quarter 2007 ($t=194$), we fit the VAR model by different methods utilizing
all the historical data available until time $t-1$ and obtain the one-step-
head forecast $\mathbf{\widehat{A}}\mathbf{x}_{t}$. Then, we calculate the
rolling forecasting errors
$\widehat{\boldsymbol{\delta}}_{t}=\mathbf{\widehat{A}}\mathbf{x}_{t}-\mathbf{y}_{t}$
in the $\ell_{2}$ and $\ell_{\infty}$ norms, respectively, and summarized the
their means and medians in Table 1.

Table 1: Mean and median of out-of-sample rolling forecasting errors of various methods for the full, sparse, and reduced-rank VAR(4) models. Model | Full | Sparse |  | Reduced-rank
---|---|---|---|---
Method | OLS | L1-RLS | L1-RYW | L1-ROB |  | RRR | NN-RLS | NN-RYW | NN-ROB
mean of $\|\widehat{\boldsymbol{\delta}}_{t}\|_{2}$ | 22.63 | 5.64 | 5.82 | 4.17 |  | 19.67 | 8.82 | 8.91 | 6.52
median of $\|\widehat{\boldsymbol{\delta}}_{t}\|_{2}$ | 14.57 | 5.21 | 5.20 | 3.81 |  | 9.62 | 5.39 | 5.34 | 3.97
mean of $\|\widehat{\boldsymbol{\delta}}_{t}\|_{\infty}$ | 9.29 | 2.16 | 2.42 | 1.66 |  | 6.42 | 3.44 | 3.27 | 2.49
median of $\|\widehat{\boldsymbol{\delta}}_{t}\|_{\infty}$ | 6.42 | 1.96 | 2.01 | 1.43 |  | 3.06 | 2.16 | 2.09 | 1.80

As shown in Table 1, the OLS estimator produces the largest out-of-sample
prediction errors as the full VAR(4) model is highly over-parameterized, while
low-dimensional and sparse structure can alleviate the problem of overfitting
and significantly improve the forecasting performance. Both L1-ROB and NN-ROB
estimators have smaller out-of-sample forecasting errors than their standard
counterparts, indicating the necessity of robust estimation for large-scale
time series data and the satisfactory performance of the proposed methodology.
The L1-ROB estimator performs best among all estimators as it produces a
parsimonious sparse parameter matrix and prevents overfitting effectively.

## References

  * Avella-Medina et al., (2018) Avella-Medina, M., Battey, H. S., Fan, J., and Li, Q. (2018).  Robust estimation of high-dimensional covariance and precision matrices.  Biometrika, 105:271–284.
  * Banna et al., (2016) Banna, M., Merlevède, F., and Youssef, P. (2016).  Bernstein-type inequality for a class of dependent random matrices.  Random Matrices: Theory and Applications, 5:1650006.
  * Basu et al., (2019) Basu, S., Li, X., and Michailidis, G. (2019).  Low rank and structured modeling of high-dimensional vector autoregressions.  IEEE Transactions on Signal Processing, 67:1207–1222.
  * Basu and Michailidis, (2015) Basu, S. and Michailidis, G. (2015).  Regularized estimation in sparse high-dimensional time series modevls.  Annals of Statistics, 43:1535–1567.
  * Belloni and Chernozhukov, (2011) Belloni, A. and Chernozhukov, V. (2011).  $\ell_{1}$-penalized quantile regression in high-dimensional sparse models.  Annals of Statistics, 39:82–130.
  * Bollerslev, (1986) Bollerslev, T. (1986).  Generalized autoregressive conditional heteroskedasticity.  Journal of Econometrics, 31:307–327.
  * Boyd et al., (2011) Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011).  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine learning, 3:1–122.
  * Bubeck et al., (2013) Bubeck, S., Cesa-Bianchi, N., and Lugosi, G. (2013).  Bandits with heavy tail.  IEEE Transactions on Information Theory, 59:7711–7717.
  * Cai et al., (2011) Cai, T., Liu, W., and Luo, X. (2011).  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association, 106:594–607.
  * Cai and Zhou, (2012) Cai, T. T. and Zhou, H. H. (2012).  Optimal rates of convergence for sparse covariance matrix estimation.  Annals of Statistics, 40:2389–2420.
  * Candès and Tao, (2007) Candès, E. and Tao, T. (2007).  The dantzig selector: Statistical estimation when p is much larger than n.  Annals of Statistics, 35:2313–2351.
  * Candès and Plan, (2011) Candès, E. J. and Plan, Y. (2011).  Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements.  IEEE Transactions on Information Theory, 57:2342–2359.
  * Catoni, (2012) Catoni, O. (2012).  Challenging the empirical mean and empirical variance: a deviation study.  Annales de l’IHP Probabilités et Statistiques, 48:1148–1185.
  * Devroye et al., (2016) Devroye, L., Lerasle, M., Lugosi, G., and Oliveira, R. I. (2016).  Sub-gaussian mean estimators.  Annals of Statistics, 44:2695–2725.
  * Doukhan, (1994) Doukhan, P. (1994).  Mixing.  In Mixing, pages 15–23. Springer.
  * Engle, (1982) Engle, R. F. (1982).  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  Econometrica, 50:987–1007.
  * Fan et al., (2017) Fan, J., Li, Q., and Wang, Y. (2017).  Estimation of high dimensional mean regression in the absence of symmetry and light tail assumptions.  Journal of the Royal Statistical Society: Series B, 79:247–265.
  * Fan et al., (2021) Fan, J., Wang, W., and Zhu, Z. (2021).  A shrinkage principle for heavy-tailed data: High-dimensional robust low-rank matrix recovery.  Annals of Statistics, 49:1239–1266.
  * Fan and Yao, (2008) Fan, J. and Yao, Q. (2008).  Nonlinear time series: nonparametric and parametric methods.  Springer Science & Business Media.
  * Gorrostieta et al., (2012) Gorrostieta, C., Ombao, H., Bédard, P., and Sanes, J. N. (2012).  Investigating brain connectivity using mixed effects vector autoregressive models.  NeuroImage, 59:3347–3355.
  * Guo et al., (2016) Guo, S., Wang, Y., and Yao, Q. (2016).  High-dimensional and banded vector autoregressions.  Biometrika, 103:889–903.
  * Hampel, (1971) Hampel, F. R. (1971).  A general qualitative definition of robustness.  The Annals of Mathematical Statistics, 42:1887–1896.
  * Hampel, (1974) Hampel, F. R. (1974).  The influence curve and its role in robust estimation.  Journal of the American Statistical Association, 69:383–393.
  * Hampel, (2001) Hampel, F. R. (2001).  Robust statistics: A brief introduction and overview.  In Research report/Seminar für Statistik, Eidgenössische Technische Hochschule (ETH), volume 94. Seminar für Statistik, Eidgenössische Technische Hochschule.
  * Han et al., (2015) Han, F., Lu, H., and Liu, H. (2015).  A direct estimation of high dimensional stationary vector autoregressions.  Journal of Machine Learning Research, 16:3115–3150.
  * Han et al., (2020) Han, Y., Tsay, R. S., and Wu, W. B. (2020).  Robust estimation of high dimensional generalized linear models for temporal dependent data.  Technical report.
  * Huber, (1964) Huber, P. J. (1964).  Robust estimation of a location parameter.  The Annals of Mathematical Statistics, 35:73–101.
  * Ke et al., (2019) Ke, Y., Minsker, S., Ren, Z., Sun, Q., and Zhou, W.-X. (2019).  User-friendly covariance estimation for heavy-tailed distributions.  Statistical Science, 34:454–471.
  * Kock and Callot, (2015) Kock, A. B. and Callot, L. (2015).  Oracle inequalities for high dimensional vector autoregressions.  Journal of Econometrics, 186:325–344.
  * Koop, (2013) Koop, G. M. (2013).  Forecasting with medium and large bayesian vars.  Journal of Applied Econometrics, 28:177–203.
  * Liebscher, (2005) Liebscher, E. (2005).  Towards a unified approach for proving geometric ergodicity and mixing properties of nonlinear autoregressive processes.  Journal of Time Series Analysis, 26:669–689.
  * Ling, (2004) Ling, S. (2004).  Estimation and testing stationarity for double-autoregressive models.  Journal of the Royal Statistical Society: Series B, 66:63–78.
  * Loh, (2017) Loh, P.-L. (2017).  Statistical consistency and asymptotic normality for high-dimensional robust M-estimators.  Annals of Statistics, 45:866–896.
  * Lütkepohl, (2005) Lütkepohl, H. (2005).  New introduction to multiple time series analysis.  Springer Science & Business Media.
  * Martin, (1981) Martin, R. D. (1981).  Robust methods for time series.  In Applied time series analysis II, pages 683–759. Elsevier.
  * Massart, (2007) Massart, P. (2007).  Concentration inequalities and model selection: Ecole d’Eté de Probabilités de Saint-Flour XXXIII-2003.  Springer.
  * Michailidis and d’Alché Buc, (2013) Michailidis, G. and d’Alché Buc, F. (2013).  Autoregressive models for gene regulatory network inference: Sparsity, stability and causality issues.  Mathematical Biosciences, 246:326–334.
  * Minsker, (2018) Minsker, S. (2018).  Sub-gaussian estimators of the mean of a random matrix with heavy-tailed entries.  Annals of Statistics, 46:2871–2903.
  * Muler, (2013) Muler, N. (2013).  Robust estimation for vector autoregressive models.  Computational Statistics & Data Analysis, 65:68–79.
  * Muler et al., (2009) Muler, N., Pena, D., and Yohai, V. J. (2009).  Robust estimation for arma models.  Annals of Statistics, 37:816–840.
  * Negahban and Wainwright, (2011) Negahban, S. and Wainwright, M. J. (2011).  Estimation of (near) low-rank matrices with noise and high-dimensional scaling.  Annals of Statistics, 39:1069–1097.
  * Negahban et al., (2012) Negahban, S. N., Ravikumar, P., Wainwright, M. J., and Yu, B. (2012).  A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers.  Statistical Science, 27:538–557.
  * Qiu et al., (2015) Qiu, H., Xu, S., Han, F., Liu, H., and Caffo, B. (2015).  Robust estimation of transition matrices in high dimensional heavy-tailed vector autoregressive processes.  In International Conference on Machine Learning, pages 1843–1851. PMLR.
  * Raskutti et al., (2011) Raskutti, G., Wainwright, M. J., and Yu, B. (2011).  Minimax rates of estimation for high-dimensional linear regression over $\ell_{q}$-balls.  IEEE Transactions on Information Theory, 57:6976–6994.
  * Stock and Watson, (2009) Stock, J. H. and Watson, M. (2009).  Forecasting in dynamic factor models subject to structural instability.  The Methodology and Practice of Econometrics. A Festschrift in Honour of David F. Hendry, 173:205.
  * Sun et al., (2020) Sun, Q., Zhou, W.-X., and Fan, J. (2020).  Adaptive huber regression.  Journal of the American Statistical Association, 115:254–265.
  * Tan et al., (2022) Tan, K. M., Sun, Q., and Witten, D. (2022).  Sparse reduced rank huber regression in high dimensions.  Journal of the American Statistical Association.  To appear.
  * Tsay, (2013) Tsay, R. S. (2013).  Multivariate time series analysis: with R and financial applications.  John Wiley & Sons.
  * Velu and Reinsel, (2013) Velu, R. and Reinsel, G. C. (2013).  Multivariate reduced-rank regression: theory and applications, volume 136.  Springer Science & Business Media.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-dimensional statistics: A non-asymptotic viewpoint, volume 48.  Cambridge University Press.
  * (51) Wang, D., Lian, H., Zheng, Y., and Li, G. (2021a).  High-dimensional vector autoregressive time series modeling via tensor decomposition.  Journal of the American Statistical Association.  To appear.
  * (52) Wang, D., Zheng, Y., and Li, G. (2021b).  High-dimensional low-rank tensor autoregressive time series modeling.  arXiv preprint arXiv:2101.04276.
  * Wang and Banerjee, (2014) Wang, H. and Banerjee, A. (2014).  Bregman alternating direction method of multipliers.  Advances in Neural Information Processing Systems, 4:2816–2824.
  * Wang et al., (2007) Wang, H., Li, G., and Jiang, G. (2007).  Robust regression shrinkage and consistent variable selection through the LAD-Lasso.  Journal of Business & Economic Statistics, 25:347–355.
  * Wang et al., (2020) Wang, L., Peng, B., Bradic, J., Li, R., and Wu, Y. (2020).  A tuning-free robust and efficient approach to high-dimensional regression.  Journal of the American Statistical Association, 115:1700–1714.
  * Wang et al., (2012) Wang, L., Wu, Y., and Li, R. (2012).  Quantile regression for analyzing heterogeneity in ultra-high dimension.  Journal of the American Statistical Association, 107:214–222.
  * Wang and Yuan, (2012) Wang, X. and Yuan, X. (2012).  The linearized alternating direction method of multipliers for dantzig selector.  SIAM Journal on Scientific Computing, 34:A2792–A2811.
  * Wong et al., (2020) Wong, K. C., Li, Z., and Tewari, A. (2020).  Lasso guarantees for $\beta$-mixing heavy-tailed time series.  Annals of Statistics, 48:1124–1142.
  * Wu and Xia, (2016) Wu, J. C. and Xia, F. D. (2016).  Measuring the macroeconomic impact of monetary policy at the zero lower bound.  Journal of Money, Credit and Banking, 48:253–291.
  * Wu, (2005) Wu, W. B. (2005).  Nonlinear system theory: Another look at dependence.  Proceedings of the National Academy of Sciences, 102:14150–14154.
  * Wu and Wu, (2016) Wu, W.-B. and Wu, Y. N. (2016).  Performance bounds for parameter estimates of high-dimensional linear models with correlated errors.  Electronic Journal of Statistics, 10:352–379.
  * Zhang, (2021) Zhang, D. (2021).  Robust estimation of the mean and covariance matrix for high dimensional time series.  Statistica Sinica, 31:797–820.
  * Zheng and Raskutti, (2019) Zheng, L. and Raskutti, G. (2019).  Testing for high-dimensional network parameters in auto-regressive models.  Electronic Journal of Statistics, 13:4977–5043.
  * Zheng and Cheng, (2021) Zheng, Y. and Cheng, G. (2021).  Finite-time analysis of vector autoregressive models under linear restrictions.  Biometrika, 108:469–489.
  * Zhu et al., (2018) Zhu, Q., Zheng, Y., and Li, G. (2018).  Linear double autoregression.  Journal of Econometrics, 207:162–174.
  * Zhu et al., (2017) Zhu, X., Pan, R., Li, G., Liu, Y., and Wang, H. (2017).  Network vector autoregression.  Annals of Statistics, 45:1096–1123.

Supplementary material for
“Rate-Optimal Robust Estimation of High-Dimensional Vector Autoregressive
Models”

Di Wang and Ruey S. Tsay

Booth School of Business, University of Chicago

This supplementary material provides all technical proofs of the theoretical
results in the main paper, some auxiliary lemmas, and details of the
linearized ADMM algorithm. Specifically, Appendices A and B present the proofs
of upper and lower bound results, respectively, and some related auxiliary
lemmas. Appendix C shows the ADMM algorithms for the $\ell_{1}$ regularized
sparse VAR model and the nuclear norm regularized reduced-rank VAR model.

##  Appendix A Proofs of Upper Bound Results

We present the proofs of the deterministic upper bounds (Propositions 1–3) in
Section A.1, the error bounds of the robust autocovariance estimators
(Propositions 4–6) in Section A.2, and the error bounds of the robust VAR
estimators (Theorems 1–4) in Section A.3, respectively. The auxiliary lemmas
are relegated to Section A.4.

We start with some notation used in the Appendix. Denote by
$\mathbb{S}^{p}=\\{\mathbf{v}\in\mathbb{R}^{p}:\|\mathbf{v}\|_{2}=1\\}$ the
$p$-dimensional sphere with unit radius in the Euclidean norm. For any
symmetric matrix $\mathbf{M}$, let $\lambda_{\max}(\mathbf{M})$ and
$\lambda_{\min}(\mathbf{M})$ be its largest and smallest eigenvalues.

###  A.1 Proof of Propositions 1–3

###### Proof of Proposition 1.

The general idea of the proof follows from that of Dantzig selector (Candès
and Tao,, 2007) and constrained minimization estimation (Cai et al.,, 2011;
Han et al.,, 2015).

By the conditions in Proposition 1,
$\mathcal{R}^{*}(\widetilde{\mathbf{\Sigma}}_{0}-\mathbf{\Sigma}_{0})\leq\zeta_{0}$
and
$\mathcal{R}^{*}(\widetilde{\mathbf{\Sigma}}_{1}-\mathbf{\Sigma}_{1})\leq\zeta_{1}$.
Let $\mathbf{\Delta}=\widehat{\mathbf{A}}-\mathbf{A}^{*}$ be the estimation
error. We first show that the true value $\mathbf{A}^{*}$ is a feasible
solution to the optimization problem in (6). Note that as
$\mathcal{R}^{*}(\cdot)$ is convex,

$\begin{split}&\mathcal{R}^{*}(\mathbf{A}^{*}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\widetilde{\Sigma}}_{1})=\mathcal{R}^{*}(\mathbf{A}^{*}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\Sigma}_{1}+\mathbf{\Sigma}_{1}-\mathbf{\widetilde{\Sigma}}_{1})\\\
=&\mathcal{R}^{*}(\mathbf{A}^{*}\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{A}^{*}\mathbf{\Sigma}_{0}+\mathbf{\Sigma}_{1}-\mathbf{\widetilde{\Sigma}}_{1})\\\
\leq&\mathcal{R}^{*}(\mathbf{A}^{*}[\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\Sigma}_{0}])+\mathcal{R}^{*}(\mathbf{\Sigma}_{1}-\mathbf{\widetilde{\Sigma}}_{1})\leq\mathcal{R}(\mathbf{A}^{*})\mathcal{R}^{*}(\mathbf{\widetilde{\Sigma}}_{0}-\mathbf{\Sigma}_{0})+\zeta_{1}\\\
\leq&\zeta_{0}\mathcal{R}(\mathbf{A}^{*})+\zeta_{1}\leq\lambda,\end{split}$
(61)

where the first inequality follows from the triangle inequality, and the
second inequality follows from the properties of $\mathcal{R}^{*}(\cdot)$ and
$\mathcal{R}(\cdot)$.

Therefore, $\mathbf{A}^{*}$ is feasible in the optimization equation, and
hence $\mathcal{R}(\widehat{\mathbf{A}})\leq\mathcal{R}(\mathbf{A}^{*})$.
Then, by the triangle inequality, we have

$\begin{split}&\mathcal{R}^{*}(\mathbf{\Delta})=\mathcal{R}^{*}(\widehat{\mathbf{A}}-\mathbf{\Sigma}_{1}\mathbf{\Sigma}_{0}^{-1})\\\
=&\mathcal{R}^{*}([\widehat{\mathbf{A}}\mathbf{\Sigma}_{0}-\mathbf{\Sigma}_{1}]\mathbf{\Sigma}_{0}^{-1})\\\
=&\mathcal{R}^{*}([\widehat{\mathbf{A}}\mathbf{\Sigma}_{0}-\widehat{\mathbf{A}}\widetilde{\mathbf{\Sigma}}_{0}+\widehat{\mathbf{A}}\widetilde{\mathbf{\Sigma}}_{0}-\widetilde{\mathbf{\Sigma}}_{1}+\widetilde{\mathbf{\Sigma}}_{1}-\mathbf{\Sigma}_{1}]\mathbf{\Sigma}_{0}^{-1})\\\