arXiv:1708.01505v1 [math.ST] 3 Aug 2017

Lasso Guarantees for β-Mixing Heavy Tailed Time
Series
Wong, Kam Chung
kamwong@umich.com

Tewari, Ambuj
tewaria@umich.com

Abstract
Most existing theoretical results for the lasso, however, require the
samples to be iid. Recent work has provided guarantees for the lasso
assuming that the time series is generated by a sparse Vector AutoRegressive (VAR) model with Gaussian innovations. Proofs of these
results rely critically on the fact that the true data generating mechanism (DGM) is a finite-order Gaussian VAR. This assumption is quite
brittle: linear transformations, including selecting a subset of variables, can lead to the violation of this assumption. In order to break
free from such assumptions, we derive non-asymptotic inequalities for
estimation error and prediction error of the lasso estimate of the best
linear predictor without assuming any special parametric form of the
DGM. Instead, we rely only on (strict) stationarity and geometrically
decaying β-mixing coefficients to establish error bounds for the lasso for
subweibull random vectors. The class of subweibull random variables
that we introduce includes subgaussian and subexponential random
variables but also includes random variables with tails heavier than an
exponential. We also show that, for Gaussian processes, the β-mixing
condition can be relaxed to summability of the α-mixing coefficients.
Our work provides an alternative proof of the consistency of the lasso
for sparse Gaussian VAR models. But the applicability of our results
extends to non-Gaussian and non-linear times series models as the examples we provide demonstrate.

1

Introduction

High dimensional statistics is a vibrant area of research in modern statistics
and machine learning [7, 20]. The interplay between computational and
statistical aspects of estimation in high dimensions has led to a variety of
1

eﬃcient algorithms with statistical guarantees including methods based on
convex relaxation (see, e.g., [8, 35]) and methods using iterative optimization
techniques (see, e.g., [4, 1, 15]). However, the bulk of existing theoretical
work focuses on iid samples. The extension of theory and algorithms in high
dimensional statistics to time series data, where dependence is the norm
rather than the exception, is just beginning to occur. We brieﬂy summarize
some recent work in Section 1.1 below.
Our focus in this paper is to give guarantees for ℓ1 -regularized least squares
estimation, or lasso [20], that hold even when there is temporal dependence
in data. We build upon the recent work of [3] who took a major step forward
in providing guarantees for lasso in the time series setting. They considered
Gaussian VAR models with ﬁnite lag (see Example 1) and deﬁned a measure
of stability using the spectral density, which is the Fourier transform of the
autocovariance function of the time series. Then they showed that one can
derive error bounds for lasso in terms of their measure of stability. Their
bounds are an improvement over previous work [33, 26, 19] that assumed
operator norm bounds on the transition matrix. These operator norm conditions are restrictive even for VAR models with a lag of 1 and never hold if
the lag is strictly larger than 1! Therefore, the results of [3] hold in greater
generality than previous work. But they do have limitations.
A key limitation is that [3] assume that the VAR model is the true data generating mechanism (DGM). Their proof techniques rely heavily on having
the VAR representation of the stationary process available. The VAR model
assumption, while popular in many areas, can be restrictive since the VAR
family is not closed under linear transformations: if Zt is a VAR process
then CZt may not expressible as a ﬁnite lag VAR [27]. We later provides
examples (Examples 4 and 2) of VAR processes where leaving out a single
variable breaks down the VAR assumption. What if we do not assume that
Zt is a ﬁnite lag VAR process but simply that it is stationary? Under stationarity (and ﬁnite 2nd moment conditions), the best linear predictor of Zt
in terms of Zt−d , . . . , Zt−1 is well deﬁned even if Zt is not a lag d VAR. If
we assume that this best linear predictor involves sparse coeﬃcient matrices,
can we still guarantee consistent parameter estimation? Our paper provides
an aﬃrmative answer to this important question.
We provide ﬁnite sample parameter estimation and prediction error bounds
for lasso in two cases: (a) for stationary processes with subweibull marginals
and geometrically decaying β-mixing coeﬃcients (Section 4), and (b) for
stationary Gaussian processes with suitably decaying α-mixing coeﬃcients

2

(Section 3). It is well known that guarantees for lasso follow if one can
establish restricted eigenvalue (RE) conditions and provide deviation bounds
(DB) for the correlation of noise with the regressors (see the Master Theorem
in Section 2.3 below for a precise statement). Therefore, the bulk of the
technical work in this paper boils down to establishing, with high probability,
that RE and DB conditions hold under the the Gaussian α-mixing and the
subweibull β-mixing assumptions respectively (Propositions 7, 8, 3, 2). Note
that RE conditions were previously shown to hold under the iid assumption
by [41] for Gaussian random vectors and by [44] for sub-Gaussian random
vectors.

1.1

Summary of Recent Work on High Dimensional Time
Series.

While we discussed the work of [3] – since ours is closely related to theirs –
we wish to emphasize that several other researchers have recently published
work on statistical analysis of high dimensional time series. [46], [55] and
[2] give theoretical guarantees assuming that RE conditions hold. As [3]
pointed out, it takes a fair bit of work to actually establish RE conditions
in the presence of dependence. [10, 11, 12] use high dimensional time series
for global macroeconomic modeling. Alternatives to lasso that have been
explored include the Dantzig selector [19], quantile based methods for heavytailed data [40], quasi-likelihood approaches [49], and two-stage estimation
techniques [13]. [16] considered the case of multiple sequences of univariate α-mixing heavy-tailed dependent data. Under a stringent condition on
the auto-covariance structure (please refer to Appendix D for details), the
paper established ﬁnite sample ℓ2 consistency in the real support for penalized least squares estimators. In addition, under mutual incoherence type
assumption , it provided sign and ℓ∞ consistency. An AR(1) example was
given as an illustration. Both [49] as well as [23] establish oracle inequalities for the lasso applied to time series prediction. [49] provided results
not just for lasso but also for estimators using penalties such as the SCAD
penalty. Also, instead of assuming Gaussian errors, it is only assumed that
fourth moments of the errors exist. [23] provided non-asymptotic lasso error
and prediction error bounds for stable Gaussian VARs. Both [45] and [30]
considered subexponential designs. [45] studied lasso on iid subexponential
designs and provide ﬁnite sample bounds. [30] studied adaptive lasso for
linear time series models and provide sign consistency results. [52] provided
theoretical guarantees for lasso in linear regression models with autoregres3

sive errors. Other structured penalties beyond the ℓ1 penalty have also been
considered [38, 37, 18, 36]. [57], [29], [53] and [9] consider estimation of the
covariance (or precision) matrix of high dimensional time series. [29] and [32]
both highlight that autoregressive (AR) estimation, even in univariate time
series, leads to high dimensional parameter estimation problems if the lag is
allowed to be unbounded.

1.2

Organization of the Paper

Section 2 introduces our notation, presents the key assumptions and also
records some useful facts needed later. Then we present two sets of high
probability guarantees for lower restricted eigenvalue and deviation bound
conditions: Section 3 covers α-mixing Gaussian time series. Note that αmixing is a weaker notion than β-mixing and all the parameter dependences
are explicit. It is followed by Section 4 which covers β-mixing time series
with subweibull observations and we make the dependence on the subweibull
norm explicit.
We present ﬁve examples, two involving α-mixing Gaussian processes, three
β-mixing sub-Gaussian vectors, and one β-mixing subweibull vectors. They
are presented along with the corresponding theoretical results to illustrate
potential applications of the theory. Examples 1 and 2 concern applications
of the results in Section 3. We consider VAR models with Gaussian innovations when the model is correctly or incorrectly speciﬁed. In Examples 3, 4,
and 5, we focus on the case of subweibull random vectors. We consider
VAR models with subweibull innovations when the model is correctly or incorrectly speciﬁed (Examples 3 and 4). In addition, we go beyond linear
models and introduce non-linearity in the DGM in Example 5.
These examples serve to illustrate that our theoretical results for lasso on
high dimensional dependent data estimation extend beyond the classical linear Gaussian setting and provides guarantees potentially in the presence of
one or more of the following scenarios: model mis-speciﬁcation, heavy tailed
non-Gaussian innovations and nonlinearity in the DGM.

2

Preliminaries

p
Consider a stochastic process of pairs (Xt , Yt )∞
t=1 where Xt ∈ R , Yt ∈
q
R , ∀t. We will be interested in time series prediction where, given a time

4

series (Zt )∞
t=1 , we might be interested in predicting Yt = Zt using Xt =
(Zt−d , . . . , Zt−1 ). As such, we cannot, and will not, assume that the pairs
(Xt , Yt ) are iid. We will assume that the sequence (Xt , Yt )Tt=1 is strictly stationary. We will be interested in estimating the best linear predictor of Yt
in terms of Xt . That is, our parameter matrix of interest Θ⋆ ∈ Rp×q , is
Θ⋆ = arg min E[ Yt − Θ′ Xt
Θ∈Rp×q

2
].
2

(2.1)

Note that Θ⋆ does not depend on t because of stationarity. Given T pairs
(Xt , Yt )Tt=1 of observations, deﬁne the matrices that collect the observations
together:
Y = (Y1 , Y2 , . . . , YT )′ ∈ RT ×q

X = (X1 , X2 , . . . , XT )′ ∈ RT ×p .

(2.2)

A matrix referred to in our analysis, but not available to an estimator, is the
matrix of residuals:
W := Y − XΘ⋆ .

(2.3)

b ∈ Rp×q be deﬁned
Let the ℓ1 -penalized least squares, or lasso, estimator Θ
as
b = arg min 1 k vec(Y − XΘ)k22 + λT kvec(Θ)k .
Θ
(2.4)
1
Θ∈Rp×q T

2.1

Notation

For scalars a and b, deﬁne shorthands a ∧ b := min{a, b} and a ∨ b :=
max{a, b}. For a symmetric matrix M, let λmax (M) and λmin (M) denote
its maximum and minimum eigenvalues respectively. For any matrix let M,
r(M), |||M|||, |||M|||
p ∞ , and′ |||M|||F denote its spectral radius maxi {|λi (M)|},
λmax (M M), entry-wise ℓ∞ norm maxi,j |Mi,j |, and Frobeoperator norm
p
nius norm tr(M′ M) respectively. For any vector v ∈ Rp , kvkq denotes
P
its ℓq norm ( pi=1 |vi |q )1/q . Unless otherwise speciﬁed, we shall use k·k to
denote
the ℓ2 norm. For any vector v ∈ Rp , we use kvk0 and kvk∞ to denote
Pp
i=1 1{vi 6= 0} and maxi {|vi |} respectively. Similarly, for any matrix M,
|||M|||0 = kvec(M)k0 where vec(M) is the vector obtained from M by concatenating the rows of M . We say that matrix M (resp. vector v) is s-sparse
if |||M|||0 = s (resp. kvk0 = s). We use v ′ and M′ to denote the transposes
of v and M respectively. When we index a matrix, we adopt the following
conventions. For any matrix M ∈ Rp×q , for 1 ≤ i ≤ p, 1 ≤ j ≤ q, we deﬁne
5

M[i, j] ≡ Mij := e′i Mej , M[i, :] ≡ Mi: := e′i M and M[:, j] ≡ M:j := Mej
where ei is the vector with all 0s except for a 1 in the ith coordinate. The
set of integers is denoted by Z.
For a lag l ∈ Z, we deﬁne the auto-covariance matrix w.r.t. (Xt , Yt )t as
Σ(l) = Σ(X;Y ) (l) := E[(Xt ; Yt )(Xt+l ; Yt+l )′ ]. Note that Σ(−l) = Σ(l)′ . Simi′ ],
larly, the auto-covariance matrix of lag l w.r.t. (Xt )t is ΣX (l) := E[Xt Xt+l
′
and w.r.t. (Yt )t is ΣY (l) := E[Yt Yt+l ]. The cross-covariance matrix at lag l
′ ]. Note the diﬀerence between Σ
is ΣX,Y (l) := E[Xt Yt+l
(X;Y ) (l) and ΣX,Y (l):
the former is a (p + q) × (p + q) matrix, the latter is a p × q matrix. Thus,
Σ(X;Y ) (l) is a matrix consisting of four sub-matrices. Using Matlab-like notation, Σ(X;Y ) (l) = [ΣX , ΣX,Y ; ΣY,X , ΣY ]. As per our convention, at lag 0,
we omit the lag argument l. For example, ΣX,Y denotes ΣX,Y (0) = E[Xt Yt′ ].
′
Finally, let Γ̂ := XTX be the empirical covariance matrix.

2.2

Sparsity, Stationarity and Zero Mean Assumptions

The following assumptions are maintained throughout; we will make additional assumptions speciﬁc to each of the subweibull and Gaussian sceb − Θ⋆ .
narios. Our goal is to provide ﬁnite sample bounds on the error Θ
We shall present theoretical guarantees on the ℓ2 parameter estimation erb − Θ⋆ )k2 and also the associated (in-sample) prediction error
ror k vec(Θ
b − Θ⋆ )′ Γ̂(Θ
b − Θ⋆ ) .
(Θ
F

Assumption 1. The matrix Θ⋆ is s-sparse, i.e., kvec(Θ⋆ )k0 ≤ s.

Assumption 2. The process (Xt , Yt ) is strictly stationary: i.e., ∀τ, n ≥ 1,
d

((X1 , Y1 ), · · · , (Xn , Yn )) = ((X1+τ , Y1+τ ), · · · , (Xn+τ , Yτ +n )).
d

where “ =” denotes equality in distribution.
Assumption 3. The process (Xt , Yt ) is centered; i.e., ∀t, E(Xt ) = 0p×1 ,
and E(Yt ) = 0q×1 .

2.3

A Master Theorem

We shall start with what we call a “master theorem" that provides nonasymptotic guarantees for lasso estimation and prediction errors under two
well-known conditions, viz. the restricted eigenvalue (RE) and the deviation
6

bound (DB) conditions. Note that in the classical linear model setting (see,
e.g., [21, Ch 2.3]) where sample size is larger than the dimensions (n > p), the
conditions for consistency of the ordinary least squares(OLS) estimator are
P

as follows: (a) the empirical covariance matrix X′ X/T → Q and Q invertible,
i.e., λmin (Q) > 0, and (b) the regressors and the noise are asymptotically
uncorrelated, i.e., X′ W/T → 0.

In high-dimensional regimes, [5], [26] and [34] have established similar consistency conditions for lasso. The ﬁrst one is the restricted eigenvalue (RE)
condition on X′ X/T (which is a special case, when the loss function is the
squared loss, of the restricted strong convexity (RSC) condition). The second
is the deviation bound (DB) condition on X′ W. The following lower RE and
DB deﬁnitions are modiﬁed from those given by [26].
Definition 1 (Lower Restricted Eigenvalue). A symmetric matrix Γ ∈ Rp×p
satisﬁes a lower restricted eigenvalue condition with curvature α > 0 and
tolerance τ (T, p) > 0 if,
∀v ∈ Rp , v ′ Γv ≥ α kvk22 − τ (T, p) kvk21 .
Definition 2 (Deviation Bound). Consider the random matrices X ∈ RT ×p
and W ∈ RT ×q deﬁned in (2.2) and (2.3) above. They are said to satisfy the
deviation bound condition if there exist a deterministic multiplier function
Q(X, W, Θ⋆ ) and a rate of decay function R(p, q, T ) such that,
1
T

X′ W

∞

≤ Q(X, W, Θ⋆ )R(p, q, T ).

We now present a master theorem that provides guarantees for the ℓ2 parameter estimation error and for the (in-sample) prediction error. The proof,
given in Appendix A builds on existing result of the same kind [5, 26, 34]
and we make no claims of originality for either the result or for the proof.
Theorem 1 (Estimation and Prediction Errors). Consider the lasso estimab defined in (2.4). Suppose Assumption 1 holds. Further, suppose that
tor Θ
Γ̂ := X′ X/T satisfies the lower RE(α, τ ) condition with α ≥ 32sτ and X′ W
satisfies the deviation bound. Then, for any λT ≥ 4Q(X, W, Θ⋆ )R(p, q, T ),
we have the following guarantees:
√
b − Θ⋆ ) ≤ 4 sλT /α,
vec(Θ

b − Θ⋆ )′ Γ̂(Θ
b − Θ⋆ )
(Θ
7

2

F

≤

32λ2T s
.
α

(2.5)
(2.6)

With this master theorem at our disposal, we just need to establish the validity of the restricted eigenvalue (RE) condition and deviation bound (DB)
conditions for stationary time series by making appropriate assumptions.
We shall do that without assuming any parametric form of the data generating mechanism. Instead, we will impose appropriate tail conditions on
the random vectors Xt , Yt and also assume that they satisfy some type of
mixing condition. Speciﬁcally, in Section 4, we will consider β-mixing subweibull random vectors (we deﬁne these below in Section 4.1). Next, in
Section 3, we will consider α-mixing Gaussian random vectors. Classically,
mixing conditions were introduced to generalize classic limit theorems in
probability beyond the case of iid random variables [42]. Recent work on
high dimensional statistics has established the validity of RE conditions in
the iid Gaussian [41] and iid sub-Gaussian cases [44]. One of the main contributions of our work is to extend these results in high dimensional statistics
from the iid to the mixing case.

2.4

A Brief Overview of Mixing Conditions

Mixing conditions [6] are well established in the stochastic processes literature as a way to allow for dependence in extending results from the iid
case. The general idea is to ﬁrst deﬁne a measure of dependence between
two random variables X, Y (that can vector-valued or even take values in a
Banach space) with associated sigma algebras σ(X), σ(Y ). For example,
α(X, Y ) = sup{|P (A ∩ B) − P (A)P (B)| : A ∈ σ(X), B ∈ σ(Y )}.

Then for a stationary stochastic process (Xt )∞
t=−∞ , one deﬁnes the mixing
coeﬃcients, for l ≥ 1,
α(l) = α(X−∞:t , Xt+l:∞ ).
We say that that the process is mixing, in the sense just deﬁned, when
α(l) → 0 as l → ∞. The particular notion we get using the α measure of
dependence above is called “α-mixing". It was ﬁrst used by [42] to extend
the central limit theorem to dependent random variables. There are other,
stronger notions of mixing, such as ρ-mixing and β-mixing that are deﬁned
using the dependence measures:
ρ(X, Y ) = sup{Cov(f (X), g(Y )) : Ef = Eg = 0, Ef 2 = Eg 2 = 1}
I

J

1 XX
β(X, Y ) = sup
|P (Si ∩ Tj ) − P (Si )P (Tj )|
2
i=1 j=1

8

where the last supremum is over all pairs of partitions {A1 , . . . , AI } and
{B1 , . . . , BI } of the sample space Ω such that Ai ∈ σ(X), Bj ∈ σ(Y ) for all
i, j. The ρ-mixing and β-mixing conditions do not imply each other but each,
by itself, implies α-mixing [6]. For stationary gaussian processes, ρ-mixing
is equivalent to α-mixing (see Fact 3 below).
The β-mixing condition has been of interest in statistical learning theory
for obtaining ﬁnite sample generalization error bounds for empirical risk
minimization [50, Sec. 3.4] and boosting [24] for dependent samples. There is
also work on estimating β-mixing coeﬃcients from data [28]. The usefulness
of β-mixing lies in the fact that by using a simple blocking technique, that
goes back to the work of [56], one can often reduce the situation to the
iid setting. At the same time, many interesting processes such as Markov
and hidden Markov processes satisfy a β-mixing condition [50, Sec. 3.5].
To the best of our knowledge, however, there are no results showing that
RE and DB conditions holds under mixing conditions. Next we ﬁll this
gap in the literature. Before we continue, we note an elementary but useful
fact about mixing conditions, viz. they persist under arbitrary measurable
transformations of the original stochastic process.
Fact 1. The range of values that the α-, β- and ρ-mixing coeﬃcients can
take on are bounded(see e.g. [6]): Consider the probability space (Ω, F, P ),
for any two sigma ﬁelds A, B ∈ F, we have
0 ≤ α(A, B) ≤ 1/4,

0 ≤ β(A, B) ≤ 1,

0 ≤ ρ(A, B) ≤ 1

Fact 2. Suppose a stationary process {Ut }Tt=1 is α, ρ, or β-mixing. Then
the stationary sequence {f (Ut )}Tt=1 , for any measurable function f (·), also
is mixing in the same sense with its mixing coeﬃcients bounded by those of
the original sequence.

3

Gaussian Processes under α-Mixing

Here we will study Gaussian processes and the α-mixing condition which
is weaker than that of the β-mixing. We make the following additional
assumption.
Assumption 4 (Gaussianity). The process (Xt , Yt ) is a Gaussian process.
Assume (Xt , Yt )Tt=1 satisﬁes Assumptions 2, 3, and 4. Note that Xt ∼
N (0, ΣX ) and Yt ∼ N (0, ΣY ). To control dependence over time, we will
assume α-mixing, the weakest notion among α, ρ and β mixing.
9

Assumption P
5 (α-Mixing). The process (Xt , Yt ) is an α-mixing process.
Let Sα (T ) := Tl=0 α(l). If α(l) is summable, we let α̃ := lim Sα (T ) < ∞.
T →∞

We will use the following useful fact [22, p. 111] in our analysis.
Fact 3. For any stationary Gaussian process, the α and ρ mixing coeﬃcients
are related as follows:
∀l ≥ 1, α(l) ≤ ρ(l) ≤ 2πα(l).

Proposition 2 (Deviation Bound, Gaussian Case). Suppose Assumptions
2–3 and 4–5 hold. Then, there exists a deterministic
positive constant c̃, and
q

a free parameter b > 0, such that, for T ≥ b+1
c̃ log(pq), we have


X′ W
≤ Q(X, W, Θ⋆ )R(p, q, T ) ≥ 1 − 8 exp(−b log(pq))
P
T
∞
where

r





(b + 1)
⋆ 2
Q(X, W, Θ ) = 8π
|||ΣX ||| 1 + max kΘ:i k2 + |||ΣY |||
1≤i≤p
c̃
r
log(pq)
.
R(p, q, T ) = Sα (T )
T
Remark 1. Note that the free parameter b serves to trade-oﬀ between the
success probability on the one hand and the sample size threshold and multiplier function Q on the other. A large b increases the success probability
but worsen the sample size threshold and the multiplier function.
⋆

Proposition 3 (RE, Gaussian Case). Suppose Assumptions 2–3 and 4–5
log(p)
hold. Then, for some constant c > 0, when the sample size T ≥ c42e
min{1,η2 } ,

then we have, with probability at least 1 − 2 exp − 2c T min{1, η 2 } , that for
every vector v ∈ Rp ,
|v ′ Γ̂v| > α1 kvk22 − τ1 (T, p)kvk21 ,

(3.1)

where
1
λmin (ΣX ),
2
λmin (ΣX )
.
η=
108πSα (T )λmax (ΣX )
α1 =

τ1 (T, p) = α/⌈c

T
min{1, η 2 }⌉,
4 log(p)

and

Remark 2. Note that, in Theorem 1, it is advantageous to have a large α
and a smaller τ so that the convergence rate is fast and the initial sample
threshold for the result to hold is small. The result above, therefore, clearly
shows that is advantageous to have a well-conditioned ΣX .
10

3.1

Estimation and Prediction Errors

Substituting the RE and DB constants from Propositions 2-3 into Theorem 1
immediately yields the following guarantees
Corollary 4 (Lasso Guarantee for Gaussian Vectors under α-Mixing). Suppose Assumptions 2–5 hold. Let c, c̃ be fixed constants from Propositions 3
and 2 and let b be free parameter. Then, for sample size
)
(
r
b+1
log(p)
max {42e, 128s} , log(pq)
T ≥ max
c min{1, η 2 }
c̃
where η =

λmin (ΣX )
108πSα (T )λmax (ΣX )


we have, with probability at least 1−2 exp − 2c T min{1, η 2 } −8 exp(−b log(pq)),
the lasso error bounds (2.5) and (2.6) hold with
1
α = λmin (ΣX )
2
λT = 4Q(X, W, Θ⋆ )R(p, q, T )
where

r





(b + 1)
⋆ 2
Q(X, W, Θ ) = 8π
|||ΣX ||| 1 + max kΘ:i k2 + |||ΣY ||| ,
1≤i≤p
c̃
r
log(pq)
.
R(p, q, T ) = Sα (T )
T
Remark 3. If the α-mixing coeﬃcients are summable,qi.e., Sα (T ) ≤ α̃ <
⋆

). Also, the
∞, ∀T , then we get the usual convergence rate of O( log(pq)
T
threshold sample size is O(s log(pq)). This is in agreement with what is
happens in the iid Gaussian case. When α(l) is not summable then both
the initial sample threshold required for the guarantee to be
q valid as well

as the rate of error decay deteriorate. The latter becomes O( Sα (T )Tlog(pq) ).
√ 
We see that as long as Sα (T ) ∈ o
T , we still have consistency. In the
ﬁnite order stable VAR case considered by [3], the α-mixing coeﬃcients are
geometrically decaying and hence summable (see Example 1 for details).

3.2

Examples

We illustrate applicability of our theory in Section 3 using the examples
below.
11

Example 1 (Gaussian VAR). Transition matrix estimation in sparse stable
VAR models has been considered by several authors in recent years [14, 19,
46]. The lasso estimator is a natural choice for the problem.
Formally a ﬁnite order Gaussian VAR(d) process is deﬁned as follows. Con+d
sider a sequence of serially ordered random vectors (Zt )Tt=1
, Zt ∈ Rp that
admits the following auto-regressive representation:
Zt = A1 Zt−1 + · · · + Ad Zt−d + Et

(3.2)

where each Ak , k = 1, . . . , d is a sparse non-stochastic coeﬃcient matrix in
Rp×p and innovations Et are p-dimensional random vectors from N (0, Σǫ ).
Assume λmin (Σǫ ) > 0 and λmax (Σǫ ) < ∞.


Pd
k
Assume that the VAR(d) process is stable; i.e. det Ip×p − k=1 Ak z 6=

′
)′ and Yt := Zt+d for
0, ∀ |z| ≤ 1. Now, we identify Xt := (Zt′ , · · · , Zt−d+1
t = 1, . . . , T .

We can verify (see Appendix E.1 for details) that Assumptions 1–5 hold.
Note that Θ⋆ = (A1 , . . . , Ad )′ ∈ Rdp×p . As a result, Propositions 2 and
3, and thus Corollary 4 follow and hence we have all the high probabilistic
guarantees for lasso on Example 1. This shows that our theory covers the
stable Gaussian VAR models for which [3] provided lasso errors bounds.
We state the following convenient fact because it allows us to study any ﬁnite
order VAR model by considering its equivalent VAR(1) representation. See
Appendix E.1 for details.
Fact 4. Every VAR(d) process can be written in VAR(1) form (see e.g. [27,
Ch 2.1]).
Therefore, without loss of generality, we can consider VAR(1) model in the
ensuing Examples.
Example 2 (Gaussian VAR with Omitted Variable). We study OLS estimator of a VAR(1) process when there are endogenous variables omitted. This
arises naturally when the underlying DGM is high-dimensional but not all
variables are available/observable/measurable to the researcher to do estimation/prediction. This also happens when the researcher mis-speciﬁes the
scope of the model.
Notice that the system of the retained set of variables is no longer a ﬁnite
order VAR(and thus non-Markovian). There is model mis-speciﬁcation and
12

this example serves to illustrate that our theory is applicable to models
beyond the ﬁnite order VAR setting.
+1
Consider a VAR(1) process (Zt , Ξt )Tt=1
such that each vector in the sequence
is generated by the recursion below:

(Zt ; Ξt ) = A(Zt−1 ; Ξt−1 ) + (EZ,t−1 ; EΞ,t−1 )
where Zt ∈ Rp , Ξt ∈ R, EZ,t ∈ Rp , and EΞ,t ∈ R are partitions of the random
vectors (Zt , Ξt ) and Et into p and 1 variables. Also,


AZZ AZΞ
A :=
AΞZ AΞΞ
is the coeﬃcient matrix of the VAR(1) process with AZΞ 1-sparse, AZZ psparse and r(A) < 1. Et := (EX,t−1 ; EZ,t−1 ) for t = 1, . . . , T + 1 are iid draws
from a Gaussian white noise process.
We are interested in the OLS 1-lag estimator of the system restricted to the
set of variables in Zt . Recall that


2
Θ⋆ := arg min E Zt − B′ Zt−1 2
B∈Rp×p

Now, set Xt := Zt and Yt := Zt+1 for t = 1, . . . , T . It can be shown that
(Θ⋆ )′ = AZZ + AZΞ ΣΞZ (0)(ΣZ )−1 . We can verify that Assumptions 1–
5 hold. See Appendix E.2 for details. As a result, Propositions 2 and 3,
and thus Corollary 4 follow and hence we have all the high probabilistic
guarantees for lasso on this non-Markovian example.

4

Subweibull Random Vectors under β-Mixing

Existing analyses of lasso mostly assume subgaussian or subexponential tails.
These assumptions ensure that the moment generating function exists, at
least for some values of the free parameter. Non-existence of the moment
generating function is often taken as the deﬁnition of having a heavy tail [17].
We now introduce a family of random variables that includes subgaussian
and subexponential random variables. In addition, it includes some heavy
tailed distributions.

13

4.1

Subweibull Random Variables and Vectors

The subgaussian and subexponential random variables are characterized by
the behavior of their tails. Among the several equivalent deﬁnitions of these
random variables, we recall the ones that are based on the growth behavior of
moments. Recall that a subgaussian (resp. subexponential) random variable
√
X can be deﬁned as one for which E(|X|p )1/p ≤ K p, ∀p ≥ 1 for some
constant K (resp. E(|X|p )1/p ≤ Kp, ∀p ≥ 1). A natural generalization of
these deﬁnitions that allows for heavier tails is as follows. Fix some γ > 0,
and require
kXkp := (E|X| p )1/p ≤ Kp1/γ , ∀p ≥ 1 ∧ γ
There are a few diﬀerent equivalent ways to imposing the condition above
including a tail condition that says that the tail is no heavier than that of a
Weibull random variable with parameter γ. That is the reason why we call
this family “subweibull(γ)”.
Lemma 5. (Subweibull properties) Let X be a random variable. Then
the following statements are equivalent for every γ > 0. The constants
K1 , K2 , K3 differ from each other at most a constant depending only on γ.
1. The tails of X satisfies
P (|X| > t) ≤ 2 exp {−(t/K1 )γ } , ∀t ≥ 0.
2. The moments of X satisfy,
kXkp := (E|X| p )1/p ≤ K2 p1/γ , ∀p ≥ 1.
3. The moment generating function of |X|γ is finite at some point; namely
E [exp (|X| /K3 )γ ] ≤ 2.
Remark 4. A similar tail condition is called “Condition C0” by [47]. However, to the best of our knowledge, this family has not been systematically
introduced. The equivalence above is related to the theory of Orlicz spaces
(see, for example, Lemma 3.1 in the lecture notes of [39]).
Definition 3. (Subweibull(γ) Random Variable and Norm). A random variable X that satisﬁes any property in Lemma 5 is called a subweibull(γ) random variable. The subweibull(γ) norm associated with X, denoted kXkψγ ,
14

is deﬁned to be the smallest constant such that the moment condition in
deﬁnition Lemma 5 holds. In other words,
kXkψγ := sup(E|X| p )1/p p−1/γ .
p≥1

It is easy to see that k·kψγ , being a pointwise supremum of norms, is indeed
a norm on the space of subweibull(γ) random variables.
Remark 5. It is common in the literature (see, for example [17]) to call a random variable heavy-tailed if its tail decays slower than that of an exponential
random variable. This way of distinguishing between light and heavy tails is
natural because the moment generating function for a heavy-tailed random
variable thus deﬁned fails to exist at any point. Note that, under such a
deﬁnition, subweibull(γ) random variables with γ < 1 include heavy-tailed
random variables.
In our theoretical analysis, we will often be dealing with squares of random
variables. The next lemma tells us what happens to the subweibull parameter
γ and the associated constant, under squaring.
Lemma 6. For any γ ∈ (0, ∞), if a random variable X is subweibull(2γ)
then X 2 is subweibull(γ). Moreover,
kX 2 kψγ ≤ 21/γ kXk2ψ2γ .
We now deﬁne the subweibull norm of a random vector to capture dependence among its coordinates. It is deﬁned using one dimensional projections
of the random vector in the same way as we deﬁne subgaussian and subexponential norms of random vectors.
Definition 4. Let γ ∈ (0, ∞). A random vector X ∈ Rp is said to be
a subweibull(γ) random vector if all of its one dimensional projections are
subweibull(γ) random variables. We deﬁne the subweibull(γ) norm of a
random vector as,
kXkψγ := sup kv ′ Xkψγ
v∈S p−1

where S p−1 is the unit sphere in Rp .
To prove our results, we need measures that control the amount of dependence in the observations across time as well as within a given time period.
15

Assumption 6. The process (Xt , Yt ) is geometrically β-mixing; i.e., there
exist constants c > 0 and γ1 > 0 such that
β(n) ≤ 2 exp(−c · nγ1 ), ∀n ∈ N.
Assumption 7. Each random vector in the sequences (Xt ) and (Yt ) follows
a subweibull(γ2 ) distribution with kXt kψγ2 ≤ KX , kYt kψγ2 < KY for t =
1, · · · , T .
Finally, we make an joint assumption on the allowed pairs γ1 , γ2 .
Assumption 8. Assume γ < 1 where


1
2 −1
γ :=
.
+
γ1 γ2
Remark 6. Note that the parameters γ1 and γ2 deﬁnes a diﬃculty landscape.
The "easy case" where γ1 ≥ 1 and γ2 ≥ 2 are already addressed in the
literature (e.g. in [54]). This paper serves to provide theoretical guarantees
for the diﬃcult scenario when the tail probability decays slowly (γ2 < 2)
and/or data exhibit strong temporal dependence (γ1 < 1) and hence extends
the literature to all spectrum of positive γ1 and γ2 .
Now, we are ready to provide high probability guarantees for the deviation
bound and restricted eigenvalue conditions.
Proposition 7 (Deviation Bound, β-Mixing Subweibull Case). Suppose Assumptions 1-3 and 6-8 hold. Let c′ > 0 be a constant and let K be defined
as
K := 22/γ2 (KY + KX (1 + |||Θ⋆ |||))2 .
2

−1

Then with sample size T ≥ C1 (log(pq)) γ , we have
!
r
1
log(pq)
′
P
≤ 2 exp(−c′ log(pq))
X W ∞ > C2 K
T
T
where the constants C1 , C2 depend only on c′ and the parameters γ1 , γ2 , c
appearing in Assumptions 6 and 7.
Proposition 8 (RE, β-mixing Subweibull Case). Suppose Assumptions 1-3
and 6-8 hold. Let
2
K := 22/γ2 KX
.

16

Then for sample size
(

 2−γ   1 )
1−γ
54K (2C1 log(p))1/γ
C2 1−γ
54K
T ≥ max
,
λmin (ΣX (0))
λmin (ΣX (0))
C1
we have with probability at least
(λmin (ΣX (0)))γ
1 − 2T exp {−c̃T } , where c̃ =
,
(54K)γ 2C1
γ

that for all v ∈ Rp ,
1
kXvk22 ≥ α kvk22 − τ kvk21 .
T


α
. Note that the constants
where α = 12 λmin (ΣX (0)) and τ = 2c̃
· log(p)
γ
T
C1 , C2 depend only on the parameters γ1 , γ2 , c appearing in Assumptions 6
and 7.

4.2

Estimation and Prediction Errors

Substituting the RE and DB constants from Propositions 7-8 into Theorem 1
immediately yields the following guarantee
Corollary 9 (Lasso Guarantee for Subweibull Vectors under β-Mixing).
Suppose Assumptions 1-3 and 6-8 hold. Let c′ , C1 , C2 , c̃ be constants as defined in Propositions 7-8, and let K := 22/γ2 (KY + KX (1 + |||Θ⋆ |||))2 ,
then for sample size
(
T ≥ max C1

2
−1
(log(pq)) γ ,

54K [2 max{8s/c̃, C1 } log(p)]1/γ
,
λmin (ΣX (0))

we have with probability at least
1 − 2T exp {−c̃T γ } − 2 exp(−c′ log(pq))
that the lasso error bounds (2.5) and (2.6) hold with
1
λmin (ΣX )
2
λT = 4Q(X, W, Θ⋆ )R(p, q, T )
α=

17



54K
λmin (ΣX (0))

 2−γ 
1−γ

C2
C1



1
1−γ

)

where
Q(X, W, Θ⋆ ) = C2 K,
r
log(pq)
R(p, q, T ) =
.
T

4.3

Examples

We explore applicability of our theory in Section 4 beyond just linear Gaussian processes using the examples below. Together, we will show that the
high probabilistic guarantees for lasso potentially extend to data generated
from DGM involving subgaussianity, model mis-speciﬁcation, and/or nonlinearity.
Example 3 (Subweibull VAR). We study a generalization of VAR, one
that has subweibull(γ2 ) realizations. Consider a VAR(1) model deﬁned as in
Example 1 except that we replace the Gaussian white noise innovations with
iid random vectors from some subweibull(γ2 ) distribution with a non-singular
covariance matrix Σǫ . Now, consider a sequence (Zt )t generated according
to the model. Then, each Zt will be a mean zero subweibull random vector.
′
)′ and Yt := Zt+d for t = 1, . . . , T .
Now, we identify Xt := (Zt′ , · · · , Zt−d+1
Assuming that Ai ’s are sparse, |||A||| < 1 and (Zt )t is stable, we can verify
(see Appendix E.1 for details) that Assumptions 1-3 and 6-8 hold. Note that
Θ⋆ = (A1 , . . . , Ad )′ ∈ Rdp×p . As a result, Propositions 7 and 8 follow and
hence we have all the high probability guarantees for lasso on Example 3.
This shows that our theory covers DGMs beyond just the stable Gaussian
processes.

Example 4 (VAR with Subweibull Innovations and Omitted Variable). Using the same setup as in Example 2 except that we replace the Gaussian
white noise innovations with iid random vectors from some subweibull(γ2 )
distribution with a non-singular covariance matrix Σǫ . Now, consider a sequence (Zt )t generated according to the model. Then, each Zt will be a mean
zero subweibull random vector.
Now, set Xt := Zt and Yt := Zt+1 for t = 1, . . . , T . Assume |||A||| < 1. It
can be shown that (Θ⋆ )′ = AZZ + AZΞ ΣΞZ (0)(ΣZ )−1 . We can verify that
Assumptions 1-3 and 6-7 hold. See Appendix E.2 for details. Therefore,
Propositions 7 and 8 and thus Corollary 9 follow and and hence we have
all the high probabilistic guarantees for subweibull random vectors from a
non-Markovian model.
18

Example 5 (Multivariate ARCH). We explore the generality of our theory
by considering a multivariate nonlinear time series model with subweibull
innovations. A popular nonlinear multivariate time series model in econometrics and ﬁnance is the vector autoregressive conditionally heteroscedastic
(ARCH) model. We chose the following speciﬁc ARCH model just for convenient validation of the geometric β-mixing property of the process; it may
potentially be applicable to a larger class of multivariate ARCH models.
+1
Let (Zt )Tt=1
be random vectors deﬁned by the following recursion, for any
constants c > 0, m ∈ (0, 1), a > 0, and A sparse with |||A||| < 1:

Zt = AZt−1 + Σ(Zt−1 )Et
Σ(z) := c · clipa,b (kzkm ) Ip×p

(4.1)

where Et are iid random vectors from some subweibull(γ2 ) distribution with
a non-singular covariance matrix Σǫ , and clipa,b (x) clips the argument x to
stay in the interval [a, b]. Consequently, each Zt will be a mean zero subweibull random vector. Note that Θ∗ = A′ , the transpose of the coeﬃcient
matrix A here.
Now, set Xt := Zt and Yt = Zt+1 for t = 1, . . . , T . We can verify (see
Appendix E.3 for details) that Assumptions 1-3 and 6-7 hold. Therefore,
Propositions 7 and 8, and thus Corollary 9 follow and and hence we have
all the high probabilistic guarantees for lasso on nonlinear models with subweibull innovations.

Acknowledgments
We thank Sumanta Basu and George Michailidis for helpful discussions, and
Roman Vershynin for pointers to the literature. We acknowledge the support
of NSF via a regular (DMS-1612549) and a CAREER grant (IIS-1452099).

References
[1] Alekh Agarwal, Sahand Negahban, and Martin J Wainwright. Fast
global convergence of gradient methods for high-dimensional statistical
recovery. The Annals of Statistics, 40(5):2452–2482, 2012.
[2] Pierre Alquier, Paul Doukhan, et al. Sparsity considerations for dependent variables. Electronic journal of statistics, 5:750–774, 2011.
19

[3] Sumanta Basu and George Michailidis. Regularized estimation in
sparse high-dimensional time series models. The Annals of Statistics,
43(4):1535–1567, 2015.
[4] Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding
algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183–202, 2009.
[5] Peter J Bickel, Ya’acov Ritov, and Alexandre B Tsybakov. Simultaneous
analysis of lasso and dantzig selector. The Annals of Statistics, pages
1705–1732, 2009.
[6] Richard C Bradley. Basic properties of strong mixing conditions. a survey and some open questions. Probability surveys, 2(2):107–144, 2005.
[7] Peter Bühlmann and Sara Van De Geer. Statistics for high-dimensional
data: methods, theory and applications. Springer Science & Business
Media, 2011.
[8] Venkat Chandrasekaran, Benjamin Recht, Pablo A Parrilo, and Alan S
Willsky. The convex geometry of linear inverse problems. Foundations
of Computational mathematics, 12(6):805–849, 2012.
[9] Xiaohui Chen, Mengyu Xu, and Wei Biao Wu. Covariance and precision matrix estimation for high-dimensional time series. The Annals of
Statistics, 41(6):2994–3021, 2013.
[10] Alexander Chudik and M Hashem Pesaran. Inﬁnite-dimensional VARs
and factor models. Journal of Econometrics, 163(1):4–22, 2011.
[11] Alexander Chudik and M Hashem Pesaran. Econometric analysis of high
dimensional VARs featuring a dominant unit. Econometric Reviews,
32(5-6):592–649, 2013.
[12] Alexander Chudik and M Hashem Pesaran. Theory and practice of
GVAR modelling. Journal of Economic Surveys, 2014.
[13] Richard A Davis, Pengfei Zang, and Tian Zheng. Sparse vector autoregressive modeling. arXiv preprint arXiv:1207.0520, 2012.
[14] Richard A Davis, Pengfei Zang, and Tian Zheng. Sparse vector autoregressive modeling. Journal of Computational and Graphical Statistics,
(just-accepted):1–53, 2015.
20

[15] David L Donoho, Arian Maleki, and Andrea Montanari. Messagepassing algorithms for compressed sensing. Proceedings of the National
Academy of Sciences, 106(45):18914–18919, 2009.
[16] JianQing Fan, Lei Qi, and Xin Tong. Penalized least squares estimation
with weakly dependent data. Science China Mathematics, 59(12):2335–
2354, 2016.
[17] Sergey Foss, Dmitry Korshunov, Stan Zachary, et al. An introduction to
heavy-tailed and subexponential distributions, volume 6. Springer, 2011.
[18] Shaojun Guo, Yazhen Wang, and Qiwei Yao. High dimensional and
banded vector autoregressions. arXiv preprint arXiv:1502.07831, 2015.
[19] Fang Han and Han Liu. Transition matrix estimation in high dimensional time series. In Proceedings of the 30th International Conference
on Machine Learning (ICML-13), pages 172–180, 2013.
[20] Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical
learning with sparsity: the lasso and generalizations. CRC Press, 2015.
[21] Fumio Hayashi. Econometrics princeton university press, 2000.
[22] IlÊźdar Abdulovič Ibragimov and Yurii Antolevich Rozanov. Gaussian
random processes. Springer, 1978.
[23] Anders Bredahl Kock and Laurent Callot. Oracle inequalities for high dimensional vector autoregressions. Journal of Econometrics, 186(2):325–
344, 2015.
[24] Sanjeev Kulkarni, Aurelie C Lozano, and Robert E Schapire. Convergence and consistency of regularized boosting algorithms with stationary β-mixing observations. In Advances in neural information processing
systems, pages 819–826, 2005.
[25] Eckhard Liebscher. Towards a uniﬁed approach for proving geometric
ergodicity and mixing properties of nonlinear autoregressive processes.
Journal of Time Series Analysis, 26(5):669–689, 2005.
[26] Po-Ling Loh and Martin J Wainwright. High-dimensional regression
with noisy and missing data: Provable guarantees with nonconvexity.
The Annals of Statistics, 40(3):1637–1664, 2012.
[27] Helmut Lütkepohl. New introduction to multiple time series analysis.
Springer Science & Business Media, 2005.
21

[28] Daniel J Mcdonald, Cosma R Shalizi, and Mark J Schervish. Estimating beta-mixing coeﬃcients. In International Conference on Artificial
Intelligence and Statistics, pages 516–524, 2011.
[29] Timothy L McMurry and Dimitris N Politis. High-dimensional autocovariance matrices and optimal linear prediction. Electronic Journal of
Statistics, 9:753–788, 2015.
[30] Marcelo C Medeiros and Eduardo F Mendes. ℓ1 -regularization of highdimensional time-series models with non-gaussian and heteroskedastic
errors. Journal of Econometrics, 191(1):255–271, 2016.
[31] Florence Merlevède, Magda Peligrad, and Emmanuel Rio. A bernstein type inequality and moderate deviations for weakly dependent sequences. Probability Theory and Related Fields, 151(3-4):435–474, 2011.
[32] Yuval Nardi and Alessandro Rinaldo. Autoregressive process modeling
via the lasso procedure. Journal of Multivariate Analysis, 102(3):528–
549, 2011.
[33] Sahand Negahban and Martin J Wainwright. Estimation of (near) lowrank matrices with noise and high-dimensional scaling. The Annals of
Statistics, pages 1069–1097, 2011.
[34] Sahand Negahban and Martin J Wainwright. Restricted strong convexity and weighted matrix completion: Optimal bounds with noise. The
Journal of Machine Learning Research, 13(1):1665–1697, 2012.
[35] Sahand N Negahban, Pradeep Ravikumar, Martin J Wainwright,
and Bin Yu. A uniﬁed framework for high-dimensional analysis
of M -estimators with decomposable regularizers. Statistical Science,
27(4):538–557, 2012.
[36] Rodrigue Ngueyep and Nicoleta Serban. Large vector auto regression
for multi-layer spatially correlated time series. Technometrics, 2014.
[37] William Nicholson, David Matteson, and Jacob Bien. VARX-L: Structured regularization for large vector autoregressions with exogenous
variables. arXiv preprint arXiv:1508.07497, 2015.
[38] William B Nicholson, Jacob Bien, and David S Matteson. Hierarchical
vector autoregression. arXiv preprint arXiv:1412.5250, 2014.

22

[39] Gilles Pisier. Subgaussian sequences in probability and fourier analysis,
2016. arXiv preprint arXiv:1607.01053v3.
[40] Huitong Qiu, Sheng Xu, Fang Han, Han Liu, and Brian Caﬀo. Robust
estimation of transition matrices in high dimensional heavy-tailed vector autoregressive processes. In Proceedings of the 32nd International
Conference on Machine Learning (ICML-15), pages 1843–1851, 2015.
[41] Garvesh Raskutti, Martin J Wainwright, and Bin Yu. Restricted eigenvalue properties for correlated gaussian designs. The Journal of Machine
Learning Research, 11:2241–2259, 2010.
[42] Murray Rosenblatt. A central limit theorem and a strong mixing condition. Proceedings of the National Academy of Sciences of the United
States of America, 42(1):43, 1956.
[43] Mark Rudelson and Roman Vershynin. Hanson-wright inequality and
sub-gaussian concentration. Electronic Communications in Probability,
18(82):1–9, 2013.
[44] Mark Rudelson and Shuheng Zhou. Reconstruction from anisotropic
random measurements. Information Theory, IEEE Transactions on,
59(6):3434–3447, 2013.
[45] Vidyashankar Sivakumar, Arindam Banerjee, and Pradeep K Ravikumar. Beyond sub-gaussian measurements: High-dimensional structured
estimation with sub-exponential designs. In Advances in Neural Information Processing Systems, pages 2206–2214, 2015.
[46] Song Song and Peter J Bickel. Large vector auto regressions. arXiv
preprint arXiv:1106.3915, 2011.
[47] Terence Tao and Van Vu. Random matrices: Sharp concentration of eigenvalues. Random Matrices: Theory and Applications,
2(03):1350007, 2013.
[48] Dag Tjøstheim. Non-linear time series and markov chains. Advances in
Applied Probability, pages 587–611, 1990.
[49] Yoshimasa Uematsu.
Penalized likelihood estimation in highdimensional time series models and uts application. arXiv preprint
arXiv:1504.06706, 2015.

23

[50] Mathukumalli Vidyasagar. Learning and generalisation: with applications to neural networks. Springer Science & Business Media, second
edition, 2003.
[51] Gabrielle Viennet. Inequalities for absolutely regular sequences: application to density estimation. Probability theory and related fields,
107(4):467–492, 1997.
[52] Hansheng Wang, Guodong Li, and Chih-Ling Tsai. Regression coefﬁcient and autoregressive order shrinkage and selection via the lasso.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(1):63–78, 2007.
[53] Zhaoran Wang, Fang Han, and Han Liu. Sparse principal component analysis for high dimensional vector autoregressive models. arXiv
preprint arXiv:1307.0164, 2013.
[54] Kam Chung Wong, Zifan Li, and Ambuj Tewari. Lasso guarantees for
time series estimation under subgaussian tails and β-mixing. arXiv
preprint arXiv:1602.04265 [stat.ML], 2017.
[55] W. B. Wu and Y. N. Wu. High-dimensional linear models with dependent observations. Electronic Journal of Statistics, 10(1):352–379,
2016.
[56] Bin Yu. Rates of convergence for empirical processes of stationary mixing sequences. The Annals of Probability, pages 94–116, 1994.
[57] Danna Zhang and Wei Biao Wu. Gaussian approximation for high dimensional time series. arXiv preprint arXiv:1508.07036, 2015.

A

Proof of Master Theorem

Proof of Theorem 1. We will break down the proof in steps.
b is optimal for 2.4 and Θ⋆ is feasible,
1. Since Θ
1
T

b
Y − XΘ

2

F

b
+ λT vec(Θ)

1

24

≤

1
|||Y − XΘ⋆ |||2F + λT kvec(Θ⋆ )k1
T

ˆ := Θ
b − Θ⋆ ∈ Rp×q
2. Let ∆
1
T

2

ˆ
X∆

F

Note that



2
⋆
′ ′
b
ˆ
≤ tr(∆ X W) + λT kvec(Θ )k1 − vec(Θ)
T
1

ˆ
vec(Θ⋆ + ∆)

1

ˆ S) }
− kvec(Θ⋆ )k1 ≥{kvec(Θ⋆S )k1 − vec(∆
1

ˆ Sc )
+ vec(∆
ˆ Sc )
= vec(∆

⋆

1

− kvec(Θ )k1

1

ˆ S)
− vec(∆

1

where S denote the support of Θ⋆ .
3. With RE constant α and
qtolerance τ , deviation bound constant Q(ΣX , ΣW )
and λT ≥ 2Q(ΣX , ΣW )
ˆ
α ∆

2

F

log(q)
T ,

we have

ˆ 2
−τ k vec(∆)k
1
1
|||X∆|||2F
T
2
ˆ ′ X′ W) + λT { vec(∆
ˆ S ) − vec(∆
ˆ Sc ) }
≤ tr(∆
T
1
1
q
X
2
ˆ Sc ) }
ˆ S ) − vec(∆
ˆ :k k1 k(X′ W):k k∞ + λT { vec(∆
k∆
≤
T
1
1
RE

≤

k=1

2
ˆ 1 X′ W
≤ k vec(∆)k
T

∞

ˆ S)
+ λT { vec(∆

1

ˆ Sc ) }
− vec(∆
1

DB

ˆ 1 Q(ΣX , ΣW )R(p, q, T ) + λT { vec(∆
ˆ S)
≤ 2k vec(∆)k

ˆ S)
ˆ 1 λN /2 + λT { vec(∆
≤ k vec(∆)k
3λT
ˆ S)
≤
vec(∆
2
ˆ
≤ 2λT vec(∆)

λT
ˆ Sc )
−
vec(∆
2
1

1

ˆ Sc ) }
− vec(∆
1

1

1

ˆ S ) ≥ vec(∆
ˆ Sc )
4. In particular, this says that 3 vec(∆
1
ˆ
ˆ
ˆ S ) ≤ 4√s vec(∆)
So vec(∆)
≤ 4 vec(∆
1

1

25

1

1

ˆ Sc ) }
− vec(∆
1

5. Finally, with α ≥ 32sτ ,
α
ˆ
vec(∆)
2

2

ˆ
≤ (α − 16sτ ) vec(∆)

F

ˆ
≤ α vec(∆)

2
F

2
F

ˆ 21
− τ k vec(∆)k

ˆ 1
≤ 2λT k vec(∆)k
√
ˆ F
≤ 2 sλT k∆k
6.
ˆ
vec(∆)

F

√
4λT s
≤
α

7. From step 4, we have
1
T

ˆ
X∆

2
F

√
ˆ
≤ 8λT s vec(∆)

Then, from step 6
1
ˆ
X∆
T

B

2
F

√
ˆ ≤ 32λ2 s/α
≤ 8λT s vec(∆)
T

Proofs for Gaussian Processes under α-Mixing

We will also need the following result to control spectral/operator norms of
matrices.
Fact 5 (Schur Test). For any matrix M, we have
|||M|||2 ≤ max kMi: k1 · max kM:j k1 .
i

j

Therefore, for any symmetric matrix M ∈ Rn×n , |||M||| ≤ max1≤i≤n kMi: k1 .

Claim 1. For any random vectors X ∈ Rn and Y ∈ Rn , we have


E XY ′

=



E Y X′

26

≤

|||ΣX ||| + |||ΣY |||
2

Proof. We have,




E XY ′
= E Y X′
:=

sup

kuk≤1, kvk≤1

=

sup
kuk≤1, kvk≤1

≤

sup
kuk≤1, kvk≤1

= sup
kuk≤1

=
≤

p

p



E u′ Y X ′ v



E (Y ′ u)(X ′ v)
p

E [(Y ′ u)2 ]

E [(Y ′ u)2 ] sup
kvk≤1

p

p

E [(X ′ v)2 ] by CauchyâĂŞSchwarz ineq.

E [(X ′ v)2 ]

p
|||E [XX ′ ]||| |||E [Y Y ′ ]|||

|||E [XX ′ ]||| + |||E [Y Y ′ ]|||
.
2

The proof of Proposition 3 relies on the following result.
Lemma 10. For a second order stationary ρ-mixing sequence of random
vectors {Xt }t , their l-th auto-covariance matrix can be bounded as follows:
|||ΣX (l)||| ≤ ρ(l)|||ΣX (0)|||, ∀ l ∈ Z.
Proof. Recall the deﬁnition of ρ-mixing, for random vectors X and Y on the
probability space (Ω, F, P), let A := σ(X) and B := σ(Y ).
ρ(X, Y ) := sup{cor(f (X), g(Y )) | f ∈ L2 (A), g ∈ L2 (B)}
≥ cor(u′ X, v ′ Y )

=p

|E[u′ Xv ′ Y

]|

∀ ﬁxed u, v
u, v non-zero

E(u′ X)2 E(v ′ Y )2

Hence, ∀u, v ﬁxed,

p
p
|u′ E[XY ′ ]v| ≤ ρ(X, Y ) E(u′ X)2 E(v ′ Y )2
r
r
sup |u′ E[XY ′ ]v| ≤ ρ(X, Y ) sup E(u′ X)2 sup E(v ′ Y )2
u

u,v

27

v

But,
′

E[XY ]

r
r
′
2
≡ sup |u E[XY ]v| ≤ ρ(X, Y ) sup E(u X) sup E(v ′ Y )2
′

′

u,v

u

v

For a stationary time series {Xt }, recall that ∀t, l
′
].
ΣX (l) := E[Xt Xt+l

By stationarity, ∀t, l

ρ(Xt , Xt+l ) = ρ(l).

Hence,
|||ΣX (l)||| ≤ ρ(l)|||ΣX (0)|||.

Proof of Proposition 2. Note that by Fact 2 α and ρ-mixing are equivalent
for stationary Gaussian processes. The proof will operate via arguments in
ρ-mixing coeﬃcients.
Recall

X′ W

∞

= max1≤i≤p,1≤j≤q |[X′ W]i,j | = max1≤i≤p,1≤j≤q X′:i W:j .

By Assumption (3), we have
EX:i = 0, ∀i and
EY:j = 0, ∀j

By ﬁrst order optimality of the optimization problem in (2.1), we have
EX′:i (Y − XΘ⋆ ) = 0, ∀i ⇒ EX:i ′ W:j = 0, ∀i, j
We know ∀i, j
X′:i W:j = X′:i W:j − E[X′:i W:j ]



1
=
kX:i + W:j k2 − E[kX:i + W:j k2 ] − kX:i k2 − E[kX:i k2 ] − kW:j k2 − E[kW:j k2 ]
2
1
1
1
kX:i + W:j k2 − E[kX:i + W:j k2 ] + kX:i k2 − E[kX:i k2 ] + kW:j k2 − E[kW:j k2 ]
≤
2
2
2
(B.1)

28

Therefore,


1
′
P
X:i W:j > 3t
T




1
1
2
2
2
2
kX:i + W:j k − E[kX:i + W:j k ] > t + P
kX:i k − E[kX:i k ] > t
≤P
2T
2T


1
2
2
+P
kW:j k − E[kW:j k ] > t
2T
This suggests proof strategy via controlling tail probability on each of the
terms kX:i k2 − E[kX:i k2 ] , kW:j k2 − E[kW:j k2 ] and kX:i + W:j k2 − E[kX:i + W:j k2 ] .
Assuming the conditions in Proposition 2, we can apply the Hanson-Wright
inequality (Lemma 11) on each of them because we know that
∀i, j

X:i ∼ N (0, ΣX:i ), i = 1, · · · , p, and
W:j := Y:j − [XΘ⋆ ]:j ∼ N (0, ΣW:j ), j = 1, · · · , q

since both {Xt }Tt=1 and {Yt }Tt=1 are centered Gaussian vectors.
So,
X:i + W:j ∼ N (0, ΣX:i + ΣW:j + ΣW:j , X:i + ΣX:i ,W:j )

We are ready to apply the tail bound on each term on the RHS of (B.1). By
Lemma (11), ∃ constant c > 0 such that ∀t ≥ 0
kX:i k2 − EkX:i k2
≤ t w.p. at least 1 − 2 exp(−cT min{t, t2 })
T |||ΣX:i |||
kW:j k2 − EkW:j k2
≤ t w.p. at least 1 − 2 exp(−cT min{t, t2 })
T ΣW:j
kX:i + W:j k2 − EkX:i + W:j k2
≤ t w.p. at least 1 − 2 exp(−cT min{t, t2 })
T ΣX:i +W:j
With Claim 1, the third inequality implies, for some c̃ > 0, that w.p. at least
1 − 8 exp(−c̃T min{t, t2 })
the following holds
kX:i + W:j k2 − EkX:i + W:j k2
kX:i + W:j k2 − EkX:i + W:j k2
≤t
≤
2T (|||ΣX:i ||| + ΣW:j )
T ΣX:i +W:j
29

Therefore,
X′:i W:j
≤ 3t w.p. at least 1 − 8 exp(−c̃T min{t, t2 })
3T (|||ΣX:i ||| + ΣW:j )
Appealing to the union bound over all i ∈ [1 · · · p] and j ∈ [1 · · · q], for any
∆
P[ max
X′:i W:j ≥ ∆] ≤ pqP[X′:i W:j ≥ ∆]
1≤i≤p,1≤j≤q

We can conclude that
X′:i W:j
≤ 3t] ≥ 1 − 8pq exp(−c̃T min{t, t2 })
1≤i≤p,1≤j≤q 3T (|||ΣX:i ||| + ΣW:j )

P[

max

q
Now, for a free parameter b > 0, chooset = (b+1)c̃Tlog(pq) , for T ≥ (b+1) c̃log(pq)
we have
"
#
r
X′ W
(b + 1) log(pq)
P
≤
max (|||ΣX:i ||| + ΣW:j )
1≤i≤p,1≤j≤q
T
c̃T
∞
≥ 1 − 8 exp[−b log(pq)]
Let’s ﬁnd out what

ΣW:j

and |||ΣX:i ||| are. Recall W = Y − XΘ⋆ . So,

ΣW:i = ΣY:i + ΣXΘ⋆:i + ΣY:i , XΘ⋆:i + ΣXΘ⋆:i , Y:i

(B.2)

By Claim 1, we have
|||ΣW:i ||| ≤ 2|||ΣY:i ||| + 2 ΣXΘ⋆:i

(B.3)

Let’s ﬁgure out each of the summands on the RHS of equation (B.3) above.


ΣXΘ⋆:i [l, k] := E (XΘ⋆:i )(XΘ⋆:i )′ [l, k]


= E e′l (XΘ⋆:i )(XΘ⋆:i )′ ek


= E (X′l,: Θ⋆:i )((Θ⋆:i )′ Xk,: )


= E (Θ⋆:i )′ Xl,: X′k,: Θ⋆:i


= (Θ⋆:i )′ EXl,: X′k,: Θ⋆:i
30

With the equality above,
≤ max

ΣXΘ⋆:i

1≤k≤T

≤2

T
X
l=0

(ΣXΘ⋆:i )[k, :]

by Fact 5

1

ρ(l) kΘ⋆:i k22 |||ΣX (0)|||

by Lemma 10

Therefore,
max

1≤i≤p

≤ max 2

ΣXΘ⋆:i

1≤i≤p

T
X
l=0

ρ(l) kΘ⋆:i k22 |||ΣX (0)||| = 2|||ΣX (0)|||

T
X
l=0

ρ(l) max kΘ⋆:i k22
1≤i≤p

Similarly,
max |||ΣY:i ||| ≤ 2|||ΣY (0)|||

T
X

ρ(l)

max |||ΣX:i ||| ≤ 2|||ΣX (0)|||

T
X

ρ(l)

1≤i≤p

and
1≤i≤p

l=0

l=0

So,. by inequality (B.3)
max |||ΣW:i ||| ≤ 4

1≤i≤q

T
X
l=0



ρ(l) |||ΣX ||| max

1≤i≤p

kΘ⋆:i k22

+ |||ΣY |||



Therefore,
max

1≤i≤p,1≤j≤q

(|||ΣX:i ||| +

ΣW:j ) ≤ 4

T
X
l=0

Finally, to state the ﬁnal result:
For a free parameter b > 0, chooset =
have with probability at least





ρ(l) |||ΣX ||| 1 + max

q

(b+1) log(pq)
,
c̃T

1≤i≤p

for T ≥

kΘ⋆:i k22



(b+1) log(pq)
c̃


+ |||ΣY |||

we

1 − 8 exp[−b log(pq)]
that
X′ W
T

∞

≤

r





T
(b + 1) log(pq) X
⋆ 2
4
ρ(l) |||ΣX ||| 1 + max kΘ:i k2 + |||ΣY |||
1≤i≤p
c̃T
l=0

31

Also, because of Fact 3, we have
r




X′ W
(b + 1) log(pq)
⋆ 2
≤
8πSα (T ) |||ΣX ||| 1 + max kΘ:i k2 + |||ΣY |||
1≤i≤p
T
c̃T
∞

Proof of Proposition 3. Note that, by Fact 3, α and ρ-mixing are equivalent
for stationary Gaussian processes. The proof will operate via arguments
involving ρ-mixing coeﬃcients.
For a ﬁxed unit test vector v ∈ Rp , kvk2 = 1, consider the Gaussian vector
Xv ∈ RT . To apply the Hanson-Wright inequality (Lemma (11)), we have
to upper bound the operator norm of the covariance matrix Q of Xv.
Q takes the form
 ′
v EX1 X1′ v · · ·

..
..

.
.

′ EX X ′ v
Q=
v
t 1


..

.
v ′ EXT X1′ v · · ·

v ′ EX1 Xj′ v
v ′ EXt Xt′ v

···
..

.
v ′ EXT X1′ v · · ·


v ′ EX1 XT′ v

..

.

′
′
v EXt XT v 


..

.
′
′
v EXT XT v

We can thus use Fact 5 and Lemma 10 to upper bound |||Q||| by
T
X

ρ(l)|||ΣX (0)|||

t=0

Now, we can apply Lemma 11 on any ﬁxed unit test vector v ∈ Rp , kvk2 = 1.
′
Recall Γ̂ := XTX ∈ Rp×p . Using Lemma 11, we have, ∀η > 0
P[|v ′ (Γ̂ − ΣX (0))v > η|||Q||||] ≤ 2 exp{−cT min(η, η 2 )} ⇒
P[v ′ (Γ̂ − ΣX (0))v > η

T
X
t=0

ρ(l)|||ΣX (0)|||] ≤ 2 exp{−cT min(η, η 2 )}

32

Using Lemma F.2 in [3], for any integer k > 0, we extend it to all vectors in
J(2k) := {v ∈ Rp : kvk ≤ 1, kvk0 ≤ 2k}:
#
"
T
X
ρ(l)|||ΣX (0)|||
P sup |v ′ (Γ̂ − ΣX (0))v| > η
v∈J(2k)

t=0

≤ 2 exp{−cT min{η, η 2 } + 2k min{log(p), log(

21ep
))}
2k

By Lemma 12 in [26], we further extend the bound to all ∀v ∈ Rp ,
′

|v (Γ̂ − ΣX (0))v| > 27η

T
X

ρ(l)|||ΣX (0)|||

t=0



kvk22

1
+ kvk21
k



21ep
))}
2k
m


T
X
1
ρ(l)|||ΣX (0)||| kvk22 + kvk21
|v ′ (Γ̂ − ΣX (0))v| ≤ 27η
k
t=0

w.p. ≤ 2 exp{−cT min(η, η 2 ) + 2k min(log(p), log(

w.p. > 1 − 2 exp{−cT min(η, η 2 ) + 2k min(log(p), log(

′

|v (Γ̂)v| > −27η

T
X

ρ(l)|||ΣX (0)|||

t=0



kvk22

1
+ kvk21
k



21ep
))}
2k
⇓

+ λmin (ΣX (0))kvk22

w.p. > 1 − 2 exp{−cT min(η, η 2 ) + 2k min(log(p), log(

21ep
))}
2k

Intuitively, we know the quadratic form of a Hermitian matrix should have
its magnitude bounded from below by its minimum eigenvalue. To achieve
. So, we have
that, pick η = PT λmin (ΣX (0))
54

t=0 ρ(l)λmax (ΣX (0))

|v ′ Γ̂v| >

λmin (ΣX (0))
1
λmin (ΣX (0))kvk22 −
kvk21
2
2k

w.p.
≥ 1 − 2 exp{−cT min(1, η 2 ) + 2k min(log(p), log(
because min(1, η 2 ) ≤ min(η, η 2 ).
33

21ep
))}
2k

Now, we choose k to make sure the ﬁrst component in the exponential domT
2
inates. For now, assume p ≥ 21ep
2k . Let k = ⌈c 4 log(p) min{1, η }⌉. Now,
choose T such that k ≥

21e
2

. Let T ≥

42e log(p)
,
c min{1,η2 }

where s is the sparsity.

log(p)
Finally, we have, for T ≥ s c42e
, with probability at least
min{1,η2 }

1 − 2 exp{−T

c
min{1, η 2 }}
2

the following holds
|v ′ Γ̂v| >
Also, let η̃ :=

C

1
λmin (ΣX (0))
λmin (ΣX (0))kvk22 −
kvk21
2
2k

λmin (ΣX (0))
108πSα (T )λmax (ΣX (0))

we can bound η with η̃ by Fact 3.

Hanson-Wright Inequality

The general statement of the Hanson-Wright inequality can be found in the
paper by [43, Theorem 1.1]. We use a form of the inequality which is derived
in the proof of Proposition 2.4 of [3] as an easy consequence of the general
result. We state the modiﬁed form of the inequality and the proof below for
completeness.
Lemma 11 (Variant of Hanson-Wright Inequality). If Y ∼ N (0n×1 , Qn×n ),
then there exists universal constant c > 0 such that for any η > 0,





1
kY k22 − E kY k22 > η|||Q||| ≤ 2 exp −cn min η, η 2 .
P
(C.1)
n
Proof. The lemma easily follows from Theorem 1.1 in [43]. Write Y =
Q1/2 X, where X ∼ N (0, I) and (Q1/2 )′ (Q1/2 ) = Q. Note that each component Xi of X is independent N (0, 1), so that kXi kψ2 ≤ 1. Then, by the
above theorem,




1
1
2
′
′
X QX − EX QX > η|||Q|||
kY k2 − Tr(Q) > η|||Q||| = P
P
n
n
"
(
)#
n2 η 2 |||Q||| nη|||Q|||
≤ 2 exp −c min
,
|||Q|||
|||Q|||2F

 2 
≤ 2 exp −c min η, η
since |||Q|||2F ≤ n|||Q|||2
Lastly, note that Tr(Q) = Tr(EY Y ′ ) = E Tr(Y Y ′ ) = E Tr(Y ′ Y ) = E Tr kY k2 =
E kY k2 .
34

D

D.1

Proofs for Subweibull Random Vectors under βMixing
Proof Related to Subweibull Properties

Proof. (of Lemma 5) Property 1 ⇒ Property 2: Since we can scale X by K1 ,
without loss of generality, we can assume K1 = 1. Then we have, for p ≥ γ,
Z ∞
p
P (|X|p ≥ u) du
E |X| =
0
Z ∞
=
P (|X| ≥ t) ptp−1 dt
using change of variable u = tp
Z0 ∞
γ
by Property 1
2e−t ptp−1 dt
≤
0
Z ∞
p−1 1−γ
2p
=
using change of variable v = tγ
e−v · v γ v γ dv
γ 0
Z
2p ∞ −v p/γ−1
e ·v
dv
=
γ 0
 
2p
p
=
·Γ
γ
γ
 
2p p p/γ
since Γ(x) ≤ xx , ∀x ≥ 1.
≤
γ γ
Therefore, for p ≥ γ,
(E |X|p )1/p ≤ 21/p (1/γ)1/p p1/p (p/γ)1/γ ≤ Cγ · p1/γ
where Cγ = 4(1/γ ∨ 1)(1/γ)1/γ . If γ ≤ 1, this covers all p ≥ 1. If γ > 1, we
have, for p = 1, . . . , ⌈γ⌉ − 1,
(E |X|p )1/p ≤ 21/p (1/γ)1/p p1/p

max

i=1,...,⌈γ⌉−1

Γ(i/γ)1/i ≤ Cγ′ ,

where Cγ′ = 4(1/γ ∨ 1) maxi=1,...,⌈γ⌉−1 Γ(i/γ)1/i . Therefore, for all p,
(E |X|p )1/p ≤ (Cγ ∨ Cγ′ ) · p1/γ .

35

Property 2 ⇒ Property 3: Without loss of generality, we can assume K2 = 1.
Using Taylor series expansion of exp(·), for some positive λ,


∞
γ p
X
E
[((λ
|X|)
)
]

E exp [(λ |X|)γ2 ] = E 1 +
p!
p=1
≤1+
=

∞
X

∞
X
(λγ γp)p
p=1

by Property 2 and Stirling’s approx.

(p/e)p

(eγλγ )p =

p=0

1
≤ 2,
1 − eγλγ

where the last inequality holds for any λ satisfying eγλγ ≤ 1/2, i.e., λ ≤
(2eγ)−1/γ . Therefore Property 2 holds with K3 = (2eγ)1/γ .
Property 3 ⇒ Property 1: Without loss of generality, we can assume K3 = 1.
For all t > 0,
P (|X| > t) = P (exp (|X|γ2 ) ≥ exp (tγ2 ))
≤ exp (− (tγ2 )) E exp (|X|γ2 )
γ2

≤ 2 exp (− (t ))

by Markov’s inequality
by Property 3.

Proof. (of Lemma 6) By deﬁnition,
X2

ψγ

= sup p−1/γ E X 2
p≥1



= sup p
p≥1

−1/(2γ)

p 1/p


1/2p 2
2p
E |X|

Now we make a change of variables p̃ := 2p. Then, we have,
X

1/γ

2
ψγ

=2

sup p̃
p̃≥2

1/γ

≤2

1/γ

=2





sup p̃
p̃≥1

−1/(2γ)

−1/(2γ)

sup p̃

−1/(2γ)

p̃≥1
1/γ

=2

kXk2ψ2γ .
36

1/p̃ 2

p̃
E |X|


1/p̃ 2
p̃
E |X|
1/p̃

E |X|p̃

!2

D.2

Subweibull Norm Under Linear Transformations

We will need the following result about changes to the subweibull norm under
linear transformations.
Lemma 12. Let X be a random vector and A be a fixed matrix. We have,
kAXkψγ ≤ |||A||| · kXkψγ
Proof. We have,
v ′ AX

kAXkψγ = sup

kvk2 ≤1

≤

ψγ

u′ X

sup
kuk2 ≤|||A|||

= |||A||| sup

D.3

(A′ v)′ X

kvk2 ≤1

ψγ

ψγ

u′ X

kuk2 ≤1

= sup

ψγ

= |||A||| kXkψγ .

Concentration Inequality for Sums of β-Mixing Subweibull Random Variables

We will state and prove a modiﬁed form of Theorem 1 of [31]. This concentration result will be used to prove the high probability guarantees on
the deviation bound (Lemma 7) and lower restricted eigenvalue (Lemma 8)
conditions.
Lemma 13. Let (Xj )Tj=1 be a strictly stationary sequence of zero mean
random variables that are subweibull(γ2 ) with subweibull constant K. Denote their sum by ST . Suppose their β-mixing coefficients satisfy β(n) ≤
2 exp(−cnγ1 ). Let γ be a parameter given by
1
1
1
=
+ .
γ
γ1 γ2
Further assume γ < 1. Then for T > 4, and any t > 1/T ,






(tT )γ
t2 T
ST
> t ≤ T exp − γ
+ exp − 2
P
T
K C1
K C2
where the constants C1 , C2 depend only on γ1 , γ2 and c.

37

(D.1)

Proof. Note that, in this proof, constants C, C1 , C2 , . . . can depend on c, γ1
and γ2 and C1 , C2 in the proof are not the same as the eventual constants
C1 , C2 that appear in the lemma statement.
Further, we will assume that K = 1. The general form then follows by scaling the random variables by 1/K and applying the lemma with t replaced
by t/K. The proof consists of two parts. First, we will state a concentration inequality of [31] and bound a certain parameter V appearing in their
inequality using the β-mixing assumption. Second, we will simplify the expression that we get directly from their concentration inequality to get a
more convenient form.
Step 1: Controlling the V parameter using β-mixing coefficients
First, recall that Theorem 1 of [31], under the condition of our lemma, gives




u2
uγ
+ exp −
P {|ST | > u} ≤ T exp −
C1
C2 (1 + T V )
(
!γ ))
(
u2
1
u(1−γ)
(D.2)
+ exp −
exp
C3 T
C4 log(u)

First of all, we need to control the quantity V that appears in the denominator of the second term of (D.2). V is a worst case measure of the partial
sum of the auto-covariances on the clipped dependent sequence (Xt )Tt=1 . It
is increasing in time horizon T and related to dimension p and sparsity s
and hence not an absolute constant. To the best of our knowledge, V is
not controllable under the weaker α-mixing condition. As [31] mention in
their Section 2.1.1, using results of [51], we have, for any β-mixing strictly
stationary sequence (Yt ) with geometrically decaying β-mixing coeﬃcients;
i.e.,
β(k) ≤ 2 exp {−ck γ1 } for any positive k
the associated quantity V can be upper bounded as
X
V ≤ EX12 + 4
E(Bk X12 )
k≥0

for some sequence (Bk ) with values in [0, 1] satisfying E(Bk ) ≤ β(k). In our
case, (Xt ) is stationary and we know that its ﬁnite moments exist because
38

of Assumption 7. Then,
X
V ≤ EX12 + 4
E(Bk X12 )
k≥0

Xq
≤ EX12 + 4
E(Bk2 )E(X14 )

Cauchy-Schwarz ineqeuality

k≥0

q
Xq
E(Bk2 )
= EX12 + 4 E(X14 )

all ﬁnite moments of X1 exist

k≥0

q
Xp
≤ EX12 + 4 E(X14 )
E(Bk )
k≥0

≤ C,

where the second to last inequality follows because p
Bk ∈ [0, 1] p
⇒ Bk2 ≤
B
inequality
β(k) ≤
p comes from the fact that E(Bk ) ≤
√k . The
 last
1 γ1
2 exp − 2 ck
⇒ ( E(Bk )) summable. Moreover since X1 is subweibull(γ2 )
with constant 1, both EX12 and E(X14 ) are bounded with constants depending
only on γ2 . Note that C depends on c, γ1 and γ2 .
Step 2: Deriving a Convenient form Eventually we will apply the
concentration inequality above with u = tT , and we will choose t such that
u = tT > 1. Under the condition that u > 1, we will now show that the
term appearing in the exponent in the third term in (D.2),
(u)(1−γ)
,
(log(u))

(D.3)

is larger than a γ-dependent constant. Along with the fact that V is a
constant in the second term, the second and third terms in (D.2) can then
be combined into one.
Let u > 1. Note that the expression (D.3) remains positive and blows up to
inﬁnity as u approaches 1 from above. Taking derivative with respect to u,
we obtain


d u(1−γ)
u−γ
1
=
(1 − γ) −
du (log(u))
log(u)
log(u)
1

Observe that the derivative is negative when u < u∗ = e 1−γ ; for u > u∗ , it
becomes positive again. Hence, the expression (D.3) reaches its minimum at
39

u∗ , where its value is,
 1 1−γ
e 1−γ
1
1−γ

= e(1 − γ),

which is positive since γ < 1.

D.4

Proofs of Deviation and RE Bounds

Proof. (of Proposition 7) Note that constants C1 , C2 , . . . can change from
line to line and depend only on γ1 , γ2 , c appearing in Assumption 6 and
Assumption 7, and on the constant c′ appearing in the high probability
guarantee.
Recall that W := Y − XΘ⋆ , and
max1≤i≤p,1≤j≤q |(X:i )′ W:j |.

X′ W

∞

= max1≤i≤p,1≤j≤q |[X′ W]i,j | =

By Assumption 3, we have
EX:i = 0, ∀i = 1, · · · , p and

EY:j = 0, ∀j = 1, · · · , q

By ﬁrst order optimality of the optimization problem in (2.1), we have
E(X:i )′ (Y − XΘ⋆ ) = 0, ∀i ⇒ E(X:i )′ W:j = 0, ∀i, j

We know ∀i, j
(X:i )′ W:j

= (X:i )′ W:j − E[(X:i )′ W:j ]



1
=
kX:i + W:j k2 − E[kX:i + W:j k2 ] − kX:i k2 − E[kX:i k2 ] − kW:j k2 − E[kW:j k2 ]
2
1
1
1
≤ kX:i + W:j k2 − E[kX:i + W:j k2 ] + kX:i k2 − E[kX:i k2 ] + kW:j k2 − E[kW:j k2 ]
2
2
2
Therefore,


1
′
P
(X:i ) W:j > 3t
T




1
1
2
2
2
2
kX:i + W:j k − E[kX:i + W:j k ] > t + P
kX:i k − E[kX:i k ] > t
≤P
2T
2T


1
+P
kW:j k2 − E[kW:j k2 ] > t
2T
40

We will control the tail probabilities for each of the terms kX:i k2 − E[kX:i k2 ] ,
kW:j k2 − E[kW:j k2 ] and kX:i + W:j k2 − E[kX:i + W:j k2 ] . Before we apply Lemma 13, we have to ﬁgure out their subweibull norms. We will ﬁrst
calculate the subweibull(γ2 ) norm of Xti , Wtj and Xti + Wtj . This will
immediate yield control of the subweibull(γ2 /2) norms of their squares via
Lemma 6.
Recall that
Wt: = Yt: − (XΘ⋆ )t:
= Yt: − Xt: Θ⋆

Therefore, we have,
by Deﬁnition 4

kWtj kγ2 ≤ kWt: kγ2

= kYt: − Xt: Θ⋆ kγ2

≤ kYt: kγ2 + kXt: Θ⋆ kγ2

≤ kYt: kγ2 + kXt: kγ2 |||Θ⋆ |||
⋆

≤ KY + |||Θ |||KX

k·kγ2 is a norm

by Lemma 12

by Assumption 7.

We also have,
kXti + Wtj kγ2 ≤ kXti kγ2 + kWtj kγ2

⋆

≤ KY + KX (1 + |||Θ |||) .

k·kγ2 is a norm

Using Lemma 6, we know that the subweibull(γ2 /2) constants of the squares
of Xti , Wtj and Xti + Wtj are all bounded by
K = 22/γ2 (KY + KX (1 + |||Θ⋆ |||))2 .
We now apply Lemma 13 three times with γ2 replaced by γ2 /2, to get, for
any t > 1/2T ,






1
(2tT )γ
4t2 T
′
(X:i ) W:j > 3t ≤ 3T exp − γ
+ 3 exp − 2
,
P
T
K C1
K C2
where γ = (1/γ1 + 2/γ2 )−1 is less than 1 by Assumption 8.
Now, taking a union bound over the pq possible values of i, j, gives us






1
(2tT )γ
4t2 T
′
P
X W ∞ > 3t ≤ 3T pq exp − γ
+ 3pq exp − 2
.
T
K C1
K C2
41

If we set,

(

t = K max C2

r

log(3pq) C1
,
(log(3T pq))1/γ
T
T

)

then the probability of the large deviation event above is at most
2 exp(−c′ log(3pq)).
Note that the constant c′ can be made arbitrarily large but aﬀects the constants C1 , C2 above.
In the expression
√ for t above, we want to ensure that two conditions are met.
First, the 1/ T term should dominate. That is, we want,
r
log(3pq)
C1
≥
(log(3T pq))1/γ ,
T
T

which, in turn, is implied by
r
C2
log(3pq)
≥
(log(3T ))1/γ
T
T

r

and
2

Both of these are met if T ≥ C3 (log(pq)) γ

−1

C2
log(3pq)
≥
(log(pq))1/γ
T
T

.

Finally, the condition t > 1/2T should be met.That is,
r
1
log(3pq)
C2
>
T
2T
which happens as soon as T ≥ C22 /4.
Proof. (of Lemma 8) Recall that X1 , · · · , Xt ∈ Rp are subweibull random
variables forming a β-mixing and stationary sequence.
Step I: Concentration for a fixed vector Now, ﬁx a unit vector v ∈
Rp , kvk2 = 1. Deﬁne real valued random variables Zt = v ′ Xt , t = 1, · · · , T .
Note that the β-mixing rate of (Zt ) is bounded by the same of (Xt ) by
≤ 22/γ2 kZt k2ψγ . Moreover,
Fact 2. From Lemma 6, we know that Zt2 ψ
γ2 /2

2

kZt kψ ≤ kXt kψγ . Therefore, we can invoke Lemma 13 for the sum ST (v) =
2
PT γ2 2
2 with γ replaced by γ /2, γ = (1/γ + 2/γ )−1 and K =
−
EZ
Z
2
2
1
2
t
t
t=1
2 to get the following bound, for T > 4 and t > 1/T ,
22/γ2 KX






ST (v)
(tT )γ
t2 T
P
> t ≤ T exp − γ
+ exp − 2
T
K C1
K C2
42

Step II: Uniform concentration over all vectors Let J(2k) denote the
set of 2k-sparse vector with Euclidean norm at most 1. Then, using union
bound arguments similar to those in Lemma F.2 of [3], we have
(
)
ST (v)
P
sup
> 3t
T
v∈J(2k)




t2 T
(tT )γ
+ k log(p) + exp − 2
+ k log(p) .
≤ exp log(T ) − γ
K C1
K C2
From the 2k-sparse set, we will extend our bound to all v ∈ Rp . To do so,
we will apply Lemma 12 in [26]. For k ≥ 1, with probability at least




(tT )γ
t2 T
1 − exp log(T ) − γ
+ k log(p) − exp − 2
+ k log(p)
(D.4)
K C1
K C2
the following holds uniformly for all v ∈ Rp


1
1
2
2
|ST (v)| ≥ 27t kvk2 + kvk1 .
T
k
Let Σ̂T (v) := T1 kXvk22 and note that EΣ̂T (v) = v ′ ΣX (0)v. Therefore,
1
T ST = Σ̂T (v) − EΣ̂T (v). Using these notations, the above inequality implies
that


1
2
2
′
Σ̂T (v) ≥ v (ΣX (0)) v − 27 · t kvk2 + kvk1
k


1
2
2
2
≥ λmin (ΣX (0)) kvk2 − 27 · t kvk2 + kvk1
k
27t
kvk21
= kvk22 (λmin (ΣX (0)) − 27t) −
k
1
λmin (ΣX (0))
= kvk22 λmin (ΣX (0)) −
kvk21 ,
2
2k
where the last line follows by picking t =

1
54 λmin (ΣX (0)).

Step III: Selecting parameters The only thing left is to set the parameter k appropriately. We want to set it so that


(tT )γ t2 T
2k log p = min
,
K γ C1 K 2 C2
43

so that the failure probability in (D.4) is at most 1 − 2T exp(−k log p). We
want the minimum above to be attained at the ﬁrst term which means we
want
  2−γ   1
K 1−γ C2 1−γ
T ≥
t
C1
Under this condition, we have
k=

(tT )γ
.
2K γ C1 log p

To ensure that k ≥ 1, we need
T ≥

54K (2C1 log(p))1/γ
λmin (ΣX (0))

To conclude, we have the following RE guarantee. For sample size
(

 2−γ   1 )
1−γ
C2 1−γ
54K (2C1 log(p))1/γ
54K
T ≥ max
,
λmin (ΣX (0))
λmin (ΣX (0))
C1
we have with probability at least

(λmin (ΣX (0)))γ
1 − 2T exp −c′ T γ , where c′ =
(54K)γ 2C1

we have, for all v ∈ Rp ,

Σ̂T (v) ≥ α kvk22 − τ kvk21

where
1
α = λmin (ΣX (0)),
2

E
E.1

τ=

α
·
2c′



log(p)
Tγ



.

Verification of Assumptions for the Examples
VAR

Note that every VAR(d) process has an equivalent VAR(1) representation
(see e.g. [27, Ch 2.1]) as
Z̃t = ÃZ̃t−1 + Ẽt
44

(E.1)

where




Zt
 Zt−1 


Z̃t :=  . 
 .. 
Zt−d+1 (pd×1)



 
Et
0
 
Ẽt :=  . 
 .. 

A1 A2 · · ·
 Ip
0
0

 0
Ip
and Ã := 
 ..
..
 .
.
0 (pd×1)
0
0 ···


Ad−1 Ad
0
0 

0
0 

..
.. 
.
. 
Ip
0 (dp×dp)
(E.2)

Because of this equivalence, justiﬁcation of Assumptions 5(Gaussian case)
and 6 (subweibull case) will operate through this corresponding augmented
VAR(1) representation.
For both Gaussian and subweibull VARs, Assumption 3 is true since the
sequences (Zt ) is centered. Second, Θ⋆ = (A1 , · · · , Ad ). So Assumption 1
follows from construction.
For the remaining Assumptions, we will consider the Gaussian and subweibull cases separately.
Gaussian VAR

(Zt ) satisﬁes Assumption 4 by model assumption.

To show that (Zt ) is α-mixing with summable coeﬃcients, we use the following facts together with the equivalence between (Zt ) and (Z̃t ) and Fact
2.
Since (Z̃t ) is stable, the spectral radius of Ã, r(Ã) < 1, hence Assumption 2
holds. Also the innovations Ẽ has ﬁnite ﬁrst absolute moment and positive
support everywhere. Then, according to (author?) [48, Theorem 4.4], (Z̃t )
is geometrically ergodic. Note here that Gaussianity is not required here.
Hence, it also applies to innovations from mixture of Gaussians.
Next, we present a standard result (see e.g. [25, Proposition 2]).
Fact 6. A stationary Markov chain {Zt } is geometrically ergodic implies
{Zt } is absolutely regular (a.k.a. β-mixing) with
β(n) = O(γ n ), γ ∈ (0, 1)
By the fact that β-mixing implies α-mixing (see Section 2.4) for a random
process, we know that α-mixing coeﬃcients decay geometrically and hence
is summable. So, Assumption 5 holds.

45

Subweibull VAR To show that (Zt ) satisﬁes Assumptions 2 and 6, we
establish that (Zt ) is geometrically ergodic. To show the latter, we use
Propositions 1 and 2 in [25] together with the equivalence between (Zt ) and
(Z̃t ) and Fact 2.
To apply Proposition 1 in [25], we check the three conditions one by one.
Condition (i) is immediate with m = 1, E = Rp , and µ is the Lebesgue
measure. For condition (ii), we set E = Rp , µ to be the Lebesgue measure,
and m̄ = ⌈inf u∈C,v∈A ku − vk2 ⌉ the minimum “distance" between the sets C
and A. Because C is bounded and A Borel, m̄ is ﬁnite. Lastly, for condition
(iii), we again let E = Rp , µ to be the Lebesgue measure, and now the
2EkẼt k
function Q(·) = k·k and the set K = {x ∈ Rp : kxk ≤
c ǫ } where
Ã . Then,

c=1−

• Recall from model assumption that
E

h

Z̃t+1

Ã

< 1; hence,


c
Ã kzk + E( E˜t+1 ) ≤ 1 −
kzk − ǫ,
2

i
Z̃t = z <

for all z ∈ E\K
• For all z ∈ K,
E

h

Z̃t+1

i
Z̃t = z <

Ã kzk + E( Ẽt+1 ) ≤

Ã

2E E˜t
cǫ

• For all z ∈ K,
0 ≤ kzk ≤

2E Ẽt
cǫ

Now, by Proposition 1 in [25], (Z̃t ) is geometrically ergodic; hence (Z̃t ) will
be stationary. Once it reaches stationarity, by Proposition 2 in the same
paper, the sequence will be β-mixing with geometrically decaying mixing
coeﬃcients. Therefore, Assumptions 2 and 6 hold.
We are left with checking Assumption 7. Let γ be the subweibull parameter
associated with (Et ).
Assume that the spectral radius of A is smaller than 1; i.e. r(A) < 1. This
is an equivalent notion of stability of VAR process.
46

By the deﬁnition of the spectral radius,
lim |||Am |||1/m = r(A) < 1

m→∞

In other words, there exists a positive integer k < ∞ such that

Ak

< 1.

By the recursive nature of the time series,
kZt kψγ ≤

Ak kZt−1 kψγ +

i=1

Ak−i kEt−k+i kψγ

Ak−i . Using stationarity, we have the

To simplify notation, let Ci :=
following
kZt kψγ ≤

k
X

kǫt kψγ

1 − |||Ak |||

k
X

ci

i=1

!

<∞

The last inequality follows because Ci < ∞, ∀i = 1, · · · , k.
Thus, the sequence (Zt ) satisﬁes Assumption 7.

E.2

VAR with Misspecification

Assumption 3 is immediate from model deﬁnitions. By the same arguments
as in Appendix E.1, (Zt , Ξt ) are stationary and so is the sub-process (Zt );
Assumption 2 holds. Again, (Zt , Ξt ) satisfy Assumption 5 (for Example 2)
and Assumption 6 (for Example 4) according to Appendix E.1. By Fact 2,
we have the same Assumptions hold for the respective sub-processes (Zt ) in
the respective cases.
To show that (Θ⋆ )′ = AZZ + AZΞ ΣΞZ (0)(ΣZ (0))−1 , consider the following
arguments. By Assumption 2, we have the auto-covariance matrix of the
whole system (Zt , Ξt ) as
Σ(Z, Ξ) =



ΣX (0) ΣXΞ (0)
ΣΞX (0) ΣΞ (0)



Recall our Θ⋆ deﬁnition from Eq. (2.1)

Θ⋆ := arg min E Zt − B ′ Zt−1
B∈Rp×p

47

2
2



Taking derivatives and setting to zero, we obtain
(Θ⋆ )′ = ΣZ (−1)(ΣZ )−1

(E.3)

Note that
ΣZ (−1) = Σ(Z, Ξ) (−1)[1 : p1 , 1 : p1 ]
′
= E (AZZ Zt−1 + AZΞ Ξt−1 + EZ,t−1 ) Zt−1

′
′
′
= E AZZ Zt−1 Zt−1
+ AZΞ Ξt−1 Zt−1
+ EZ,t−1 Zt−1

= AZZ ΣZ (0) + AZΞ ΣΞZ (0)



by Assumptions 2 and the fact that the innovations are iid.
Naturally,
(Θ⋆ )′ = AZZ ΣZ (0)(ΣZ (0))−1 +AZΞ ΣΞZ (0)(ΣZ (0))−1 = AZZ +AZΞ ΣΞZ (0)(ΣZ (0))−1
Remark 7. Notice that AZΞ is a column vector and suppose it is 1-sparse,
and AZZ is p-sparse, then Θ⋆ is at most 2p-sparse. So Assumption 1 can be
built in by model construction.
Remark 8. We gave an explicit model here where the left out variable Ξ
was univariate. That was only for convenience. In fact, whenever the set
of left-out variables Ξ aﬀect only a small set of variables Ξ in the retained
system Z, the matrix Θ⋆ is guaranteed to be sparse. To see that, suppose
Ξ ∈ Rq and AZΞ has at most s0 non-zero rows (and let AZZ to be s-sparse
as always), then Θ⋆ is at most (s0 p + s)-sparse.
Lastly, for Example 2, the sub-process (Zt ) is Gaussian because is obtained
from a linear transformation of (Zt , Ξt ) which is Gaussian; we have Assumption 4. For Example 4, note that Zt = M(Zt , Ξt ) where M = [Ip , 0; 0′ , 0] is
a sub-setting matrix that selects the ﬁrst p entries of a (p + 1)-dimensional
vector. Hence, the fact that Zt is subweibull follows from the same arguments in Appendix E.1 pertaining to establishing the subweibull property in
conjunction with applying Lemma 12 on Zt = M(Zt , Ξt ); so, Assumption 7
holds.
Remark 9. Any VAR(d) process has an equivalent VAR(1) representation
(Lutkepohl 2005). Our results extend to any VAR(d) processes.

E.3

ARCH

Verifying the Assumptions.
To show that Assumption 6 hold for a
process deﬁned by Eq. (4.1) we leverage on Theorem 2 from [25]. Note
48

that the original ARCH model in [25] assumes the innovations to have positive support everywhere. However, this is just a convenient assumption to
establish the ﬁrst two conditions in Proposition 1 (on which proof of Theorem 2 relies) from the same paper. ARCH model with innovations from
more general distributions (e.g. uniform) also satisﬁes the ﬁrst two conditions of Proposition 1 by the same arguments in the Subweibull paragraph
of Appendix E.1.
Theorem 2 tells us that for our ARCH model, if it satisﬁes the following conditions, it is guaranteed to be absolutely regular with geometrically decaying
β-coeﬃcients.
• Et has positive density everywhere on Rp and has identity covariance
by construction.
• Σ(z) = o(kzk) because m ∈ (0, 1).
•

Σ(z)−1

≤ 1/(ac), |det (Σ(z)) | ≤ bc

• r(A) ≤ |||A||| < 1
So, Assumption 6 is valid here. We check other assumptions next.
Mean 0 is immediate, so we have Assumption 3. When the Markov chain
did not start from a stationary distribution, geometric ergodicity implies
that the sequence is approaching the stationary distribution exponentially
fast. So, after a burning period, we will have Assumption 2 approximately
valid here.
The subweibull constant of Σ(Zt−1 )Et given Zt−1 = z is bounded as follows:
for every z,

kΣ(z)Et kψγ ≤ |||Σ(z)||| kEt kψγ

by Lemma 12

≤ Ke cb =: KE

where Ke := sup kEt kψγ
t

By the recursive relationship of (Zt )t , we have
kZt kψγ ≤ |||A||| kZt−1 kψγ + KE .
which yields the bound kZt kψγ ≤ KE /(1 − |||A|||) < ∞. Hence Assumption 7
holds.
49

We will show below that Θ⋆ = A′ . Hence, sparsity (Assumption 1) can be
built in when we construct our model 4.1.
Recall Eq. E.3 from Appendix E.2 that
Θ⋆ = ΣZ (−1)(ΣZ )−1
Now,
′
ΣZ (−1) = EZt Zt−1

= E (AZt−1 +
=

by stationarity
′
Σ(Zt−1 )Et ) Zt−1

= AΣZ +
= AΣZ +
= AΣZ +
= AΣZ

′
+ EΣ(Zt−1 )Et Zt−1
′
E[c clipa,b (kZt−1 km ) Et Zt−1
]
m
′
E[cEt Zt−1 clipa,b (kZt−1 k )]
 ′

cE [Et ] E Zt−1
clipa,b (kZt−1 km )

Eq. (4.1)

′
AEZt−1 Zt−1

i.i.d. innovations
Et mean 0,

where clipa,b (x) := min{max{x, a}, b} for b > a.
Since ΣZ is invertible, we have (Θ⋆ )′ = ΣZ (−1)(ΣZ )−1 = A.

50

