Electronic Journal of Statistics
Vol. 0 (2006) 1–8
ISSN: 1935-7524

arXiv:1604.04002v5 [] 30 Sep 2017

Sparse transition matrix estimation for
high-dimensional and locally stationary
vector autoregressive models∗
Xin Ding‡ , Ziyi Qiu†? , Xiaohui Chen†
‡ University

of British Columbia
of Illinois at Urbana-Champaign
? University of Chicago
e-mail: xin.ding@stat.ubc.ca; ziyiqiu@illinois.edu; xhchen@illinois.edu
† University

Abstract: We consider the estimation of the transition matrix in the highdimensional time-varying vector autoregression (TV-VAR) models. Our
model builds on a general class of locally stationary VAR processes that
evolve smoothly in time. We propose a hybridized kernel smoothing and `1 regularized method to directly estimate the sequence of time-varying transition matrices. Under the sparsity assumption on the transition matrix, we
establish the rate of convergence of the proposed estimator and show that
the convergence rate depends on the smoothness of the locally stationary
VAR processes only through the smoothness of the transition matrix function. In addition, for our estimator followed by thresholding, we prove that
the false positive rate (type I error) and false negative rate (type II error)
in the pattern recovery can asymptotically vanish in the presence of weak
signals without assuming the minimum nonzero signal strength condition.
Favorable finite sample performances over the `2 -penalized least-squares
estimator and the unstructured maximum likelihood estimator are shown
on simulated data. We also provide two real examples on estimating the dependence structures on financial stock prices and economic exchange rates
datasets.
MSC 2010 subject classifications: Primary 62M10, 62H12; secondary
91B84.
Keywords and phrases: vector autoregression, time-varying parameters,
locally stationary processes, kernel smoothing, high-dimension, sparsity.

Contents
1
2
3

4

Introduction . . . . . . . . . . . . . .
Method and algorithm . . . . . . . .
Asymptotic properties . . . . . . . .
3.1 Locally stationary VAR process
3.2 Rate of convergence . . . . . .
3.3 Pattern recovery . . . . . . . .
Simulation studies . . . . . . . . . .
4.1 Data generation . . . . . . . . .
4.2 Tuning parameter selection . .

. . . .
. . . .
. . . .
in L2
. . . .
. . . .
. . . .
. . . .
. . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

2
4
6
6
7
8
10
11
11

∗ Research partially supported by NSF grant DMS-1404891 and UIUC Research Board
Award RB15004.

1
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

4.3

Comparison results . . . . . . . . . . . . . . . . . . .
4.3.1 Estimation errors . . . . . . . . . . . . . . . .
4.3.2 Pattern recovery . . . . . . . . . . . . . . . .
5 Real data analysis . . . . . . . . . . . . . . . . . . . . . .
5.1 Finance data: stock prices . . . . . . . . . . . . . . .
5.2 Economic data: exchange rates . . . . . . . . . . . .
6 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . .
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A Detailed setup of the simulation studies . . . . . . . . . .
B Examples for four graph structures . . . . . . . . . . . . .
C ROC curves based on the true thresholding parameters u]
D Selected stocks and currencies in real data analysis . . . .

2

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

12
12
12
12
12
19
22
27
28
30
31
32
33

1. Introduction
Vector autoregression (VAR) is a basic tool in multivariate time series analysis
and it has been extensively used to model the cross-sectional and serial dependence in various applications from economics and finance [31, 3, 4, 10, 32, 14,
34, 13, 1]. There are two fundamental limitations of the vector-autoregressive
models. First, conventional methods to estimate the transition matrix of a VAR
model are based on the least squares (LS) estimator and the maximum likelihood estimator (MLE), in which the parameter estimation is consistent when
the sample size increases and the model size is fixed [8, 25]. Since the number
of parameters grows quadratically in the number of time series variables, the
VAR model typically includes no more than ten variables in many real applications [18]. However, due to the recent explosive data enrichment, analysis of
panel data with a few hundred variables is often encountered, in which the LS
and MLE are not suitable even for a moderate problem size. Second, stationarity plays a major role in the VAR model. Therefore, the stationary VAR does
not capture the time-varying underlying data generation structures, which have
been observed in a broad regime of applications in economics and finance [2, 27].
Motivated by the limitations of the VAR, this paper studies the estimation
problems of the time-varying VAR (TV-VAR) model for high-dimensional time
series data. Let Xd×n = (x1 , · · · , xn ) be a sequence of d-dimensional observations generated by a mean-zero TV-VAR of order 1 (TV-VAR(1))
xi = A(i/n)xi−1 + ei ,

(1.1)

where A(t), t ∈ [0, 1], is a d × d matrix-valued deterministic function consisting
of the transition matrices Ai := A(i/n) at evenly spaced time points and ei
are independent and identically distributed (iid) mean-zero random errors, i.e.
innovations. In this paper, our main focus is to estimate the transition matrices
Ai for the TV-VAR(1) and the extension to higher-order VAR is straightforward.
Indeed, for a general TV-VAR of order k ≥ 1
xi = Ai,1 xi−1 + Ai,2 xi−2 + · · · + Ai,k xi−k + ei ,
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

3

>
>
we can rewrite zi = (x>
at time i = k, k + 1, · · · , n, as a TVi , · · · , xi−k+1 )
VAR(1) in the augmented space

zi = Ai zi−1 + ẽi ,
where




Ai = 



Ai,1
Id×d
0d×d
..
.

Ai,2
0d×d
Id×d
..
.

···
···
···
..
.

Ai,k−1
0d×d
0d×d
..
.

Ai,k
0d×d
0d×d
..
.

0d×d

0d×d

···

Id×d

0d×d





,





ei

 0d×1

ẽi =  .
 ..
0d×1




,


Id×d is the d × d identity matrix and 0d×d0 is the d × d0 zero matrix. Then, we
need only to estimate the first d rows in Ai .
To make the estimation problem feasible for high-dimensional time series
when d is large, it is crucial to carefully regularize the coefficient matrices Ai .
The main idea is to use certain low-dimensional structures in Ai and the degree
of freedom of (1.1) can be dramatically reduced. In our problem, we assume
two key structures in Ai . First, since our goal is to estimate a sequence of
transition matrices, which can be viewed as the discrete version of a matrixvalued function, it is natural to impose the smoothness condition to A(t). In
this case, (1.1) is closely related to the locally stationary processes, a general
class of non-stationary processes proposed by [15]. Examples of other linear and
nonlinear locally stationary processes can be found in [11]. In particular, let
i0 = 1, · · · , n, be any time point and t0 = i0 /n. Then, under suitable regularity
conditions [16], there exists a stationary process x̃i (t0 ) such that for all j =
1, · · · , d,
|Xij − X̃ij (t0 )| = OP (|i/n − t0 | + 1/n),
where Xij and X̃ij (t0 ) are the j-th element of xi and x̃i (t0 ) respectively, and
x̃i (t0 ) = A(t0 )x̃i−1 (t0 ) + ei .
Therefore, xi is an approximately stationary VAR process x̃i (t0 ) in a small
neighborhood of t0 .
Second, at each time point i = 1, · · · , n, we need to estimate a d × d matrix.
There have been different structural options in literature, such as the sparse
[21, 35, 33, 30], banded [19], and low-rank [22] transition matrix, all of which
only considered the stationary VAR model. In this paper, we consider a sequence
of sparse transition matrices and we allow the sparsity patterns to change over
time. Note that our problem is also different from the emerging literature for
estimating the high-dimensional covariance matrix and its related functionals
of time series data [11, 12, 37, 5], since our goal is to directly estimate the data
generation mechanism specified by A(·).
To simultaneously address the two issues, we propose a hybridized method
of the nonparametric smoothing technique and `1 -regularization to estimate the
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

4

sparse transition matrix of the locally stationary VAR. The proposed method
is equivalent to solving a sequence of large numbers of linear programs and
therefore the estimates can be efficiently obtained by using the high-performance
(i.e. parallel) computing technology; see Section 2 for details. In Section 3, we
establish the rate of convergence under suitable assumptions on the smoothness
of the covariance function and the sparsity of the transition matrix function A(·).
Specifically, the dimension d is permitted to increase sub-exponentially fast in
the sample size n to obtain consistent estimation, i.e. d = o(exp(n)). In addition,
we also prove that when our estimator is followed by thresholding, type I and
type II errors in the pattern recovery asymptotically vanish in the presence
of weak signals. In contrast with the existing literature on consistent model
selection, we do not require the minimum nonzero signal strength to be bounded
away from zero. Simulation studies in Section 4 and two real data examples in
Section 5 demonstrate favorable performances and more interpretable results of
the proposed TV-VAR model. Technical proofs are deferred to Section 6.
We fix the notations that will be used in the rest of the paper. Let M be a
generic matrix, I ⊂ N is an index subset, [M ]∗I and [M ]I∗ be the submatrix
of M with columns and rows indexed by I, resp. Denote J = {1, · · · , d}. Let
ρ(M ) = sup|v|=1 |M v| be the spectral norm of M . |M |1 , |M |∞ and |M |F are
Pd
the entry `1 , `∞ and Frobenius norm of M , resp. |M |`1 = maxj≤d k=1 |Mjk |
Pd
and |M |`∞ = maxk≤d j=1 |Mjk | are the matrix `1 and `∞ norm of M , resp.
supp(M ) = {(`, j) : M`j 6= 0} is the support of M and |S| is the number
of nonzeros in the support set S. For q > 0 and a random variable (r.v.) X,
kXkq = (E|X|q )1/q and we say that X ∈ Lq iff kXkq < ∞. Write kXk = kXk2 .
Denote a ∧ b = min(a, b) and a ∨ b = max(a, b). If a and b are vectors, then the
maximum and minimum are interpreted element-wise.
2. Method and algorithm
For m ∈ Z, let Σi,m = Cov(xi , xi+m ) = E(xi x>
i+m ). Then
Σi−1,1 = Σi−1,0 A>
i

and

Σi,−1 = Ai Σi−1,0 .

(2.1)

Therefore, for any estimator, say Âi , that is reasonably close to the true coefficient matrices Ai , we must have (Σi−1,1 − Σi−1,0 Â>
i ) ≈ 0d×d and (Σi,−1 −
Âi Σi−1,0 ) ≈ 0d×d . A naive estimator for Ai would be constructed to invert the
sample versions of (2.1). Estimators of this kind do not have good statistical
properties in high-dimensions because of the ill-conditioning and the dependence
information in Ai is not directly used in estimation. If Ai is known to be sparse
as a priori, we may consider the following constrained minimization program
minimizeΛ∈Rd×d
subject to

|Λ|1

(2.2)

|Σi−1,1 − Σi−1,0 Λ|∞ ≤ τ
|Σi,−1 − Λ> Σi−1,0 |∞ ≤ τ.

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

5

Because Σi,1 and Σi,0 are unknown, the solution of (2.2) is an oracle estimator
and therefore it is not implementable in practice.
Let A(·) be the continuous version of Ai , i = 1, · · · , n, in (1.1), and fix a
t ∈ (0, 1). To estimate A(t), we first estimate Σi,1 and Σi,0 in (2.1) by their
empirical versions. Let
Σ̂t,` =

n
X

w(t, m)xm x>
m+` ,

` = −1, 0, 1,

(2.3)

m=1

be the kernel smoothed estimator of Σi,0 , Σi,1 and Σi,−1 , where w(t, m) are
the nonnegative weights. Here, we consider the Nadaraya-Watson smoothed
estimator
K( t−m/n
bn )
,
(2.4)
w(t, m) = Pn
t−`/n
`=1 K( bn )
where K(·) is a nonnegative bounded kernel function with support in [−1, 1] such
R1
that −1 K(v)dv = 1 and bn is the bandwidth parameter. We assume throughout
the paper that the bandwidth satisfies the natural conditions: bn = o(1) and
n−1 = o(bn ). Then, our estimator Â(t) is defined as the transpose of the solution
of
minimizeΛ∈Rd×d
subject to

|Λ|1

(2.5)

|Σ̂i−1,1 − Σ̂i−1,0 Λ|∞ ≤ τ
|Σ̂i,−1 − Λ> Σ̂i−1,0 |∞ ≤ τ.

Observe that (2.5) is equivalent to the following d optimization sub-problems
β̂ j (t)

=

subject to

arg min |u|1
u∈Rd

(2.6)

|[Σ̂i−1,1 ]∗j − Σ̂i−1,0 u|∞ ≤ τ
|[Σ̂i,−1 ]j∗ − u> Σ̂i−1,0 |∞ ≤ τ

in that [Â(t)]j∗ = β̂ j (t)> . Since the d sub-problems (2.6) can be independently
solved, we can efficiently compute the solution of (2.5) by parallelizing the optimizations of (2.6). In addition, each sub-problem in (2.6) can be recasted as a
linear program (LP)
minimizev,u∈Rd+
subject to

|v|1 + |w|1
−Σ̂i−1,0 v + Σ̂i−1,0 w ≤ τ − max{[Σ̂i−1,1 ]∗j , [Σ̂i−1,−1 ]>
∗j }
Σ̂i,0 v − Σ̂i−1,0 w ≤ τ + min{[Σ̂i,1 ]∗j , [Σ̂i,−1 ]>
∗j }

so that β̂ j (t) = v̂ − ŵ, where (v̂, ŵ) is the nonnegative solution of the LP. Note
that (2.6) can be efficiently solved by the simplex algorithm [26]. Compared with
the estimation of sparse transition matrix for the stationary VAR model [21],
our estimator requires to solve two sample versions of Yule-Walker equations in
(2.1) and (2.5) due to the non-stationarity and therefore two-sided constraints
of the estimator Â(t) are needed in (2.5).
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

6

3. Asymptotic properties
3.1. Locally stationary VAR process in L2
To establish an asymptotic theory for estimating the continuous matrix-valued
transition function A(·) of the non-stationary VAR, it is more convenient to
model the time series as realizations from a continuous d-dimensional process.
Following the framework of [17], x1 , · · · , xn are viewed as realizations of the
continuous d-dimensional process x(t) = (X1 (t), · · · , Xd (t))> , t ∈ [0, 1], on the
discrete grid ti = i/n, i = 1, · · · , n, i.e. x(i/n) = A(i/n)x((i − 1)/n) + e(i/n)
and the temporal dependence of xi is carried out on the rescaled time indices
ti .
Definition 3.1 (Locally stationary VAR process in L2 ). A mean-zero nonstationary VAR process {xi ; i = 1, · · · , n} in Rd of the form (1.1) is said to be
locally stationary in L2 , or weakly locally stationary, if for every i = 1, · · · , n and
t = i/n, there exists a mean-zero stationary VAR process {x̃m (t); m = 1, · · · , n}
given by
x̃m (t) = A(t)x̃m−1 (t) + em
(3.1)
such that for all m = 1, · · · , n,

max kX̃mj (t) − Xmj k ≤ C

1≤j≤d

1
m
−t +
n
n


(3.2)

and the constant C > 0 does not depend on d, m, and n.
Note that the approximating stationary VAR process in Definition 3.1 generally depends on t. As suggested by the Yule-Walker equations (2.1), given a
fixed time of interest, estimation of A(t) relies on an estimate of Σt,` . Consider
` = 0 and let Σ(t) = Σt,0 be the matrix-valued covariance function at lag-zero.
With a finite number of observations x1 , · · · , xn , it is unclear how to define the
covariance function Σ(·) off the n points ti = i/n. Nevertheless, the weakly locally stationary VAR processes provide a natural framework for extending Σ(ti )
to Σ(t) for all t ∈ (0, 1).
Since A(·) is continuous on [0, 1], the stationary VAR process in (3.1) is
defined for all t ∈ (0, 1). Let Σ̃(t) = E[x̃m (t)x̃m (t)> ]. By (3.2) and the CauchySchwarz inequality, we have for each j, k = 1, · · · , d,
|σjk (ti ) − σ̃jk (ti )| =

|E(Xij Xik ) − E(X̃ij (ti )X̃ik (ti ))|

≤

kXij k2 · kXik − X̃ik (ti )k + kX̃ik (ti )k · kXij − X̃ij (ti )k

≤

(kXij k2 + kX̃ik (ti )k) max kXij − X̃ij (ti )k

≤

1≤j≤d

−1

Cn

,

where the constant C here is uniform in j, k = 1, · · · , d and i = 1, · · · , n.
Therefore, we get
Σ(ti ) = Σ̃(ti ) + O(n−1 ).
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

7

Letting n → ∞ and by the continuity of A(·), we can extend Σ(t) = Σ̃(t) for all
t ∈ (0, 1). Similar extension can be done for Σt,` for ` = ±1. In Section 3.2, it
will be shown that the asymptotic theory of estimating A(t), t ∈ (0, 1) depends
on the smoothness of the weakly locally stationary VAR processes only through
the smoothness of Σ(t) and therefore A(t).
3.2. Rate of convergence
In this section, we characterize the rate of convergence of our estimator (2.5)
under various matrix norms. We assume d ≥ 2. To study the asymptotic properties of the proposed estimator, we make the following assumptions on model
(1.1).
1. The coefficient matrices are sparse: let 0 ≤ α < 1, for each i = 1, · · · , n,
Ai ∈ Gα (s, Md ) = Gα| (s, Md ) ∩ Gα− (s, Md ),

(3.3)

where
Gα| (s, Md )

d
n
o
X
d×d
|Θ`j |α ≤ s, |Θ|`1 ≤ Md ,
=
Θ∈R
: max

Gα− (s, Md )

d
n
o
X
=
Θ ∈ Rd×d : max
|Θ`j |α ≤ s, |Θ|`∞ ≤ Md .

j≤d

`≤d

`=1

j=1

2. The marginal and lag-one covariance matrix processes {Σ0 (t)}t∈[0,1] and
{Σ1 (t)}t∈[0,1] are smooth functions in Rd×d : Σ`,jk (t) ∈ C 2 ([0, 1]), ` = 0, 1,
where C 2 ([0, 1]) is the class of functions defined on [0, 1] that are twice
differentiable with bounded derivatives uniformly in j, k = 1, · · · , d.
3. The random innovations ei = (ei1 , · · · , eid )> have iid components and
2
sub-Gaussian tail: P(|eij | ≥ x) ≤ Ce−cx for all x > 0.
Before proceeding, we discuss the above assumptions. (3.3) requires that the
transition matrices are sparse in both columns and rows at all time points.
A similar matrix class defined by (3.3) was first proposed in [7] for symmetric
matrices and it has been widely used for estimating high-dimensional covariance
and precision matrix; see e.g. [9, 11]. If α = 0, then the maximum number of
nonzeros in columns and rows of Ai is at most s. Assumption 2 requires that
the marginal and lag-one covariance matrices evolve smoothly in time. The
smoothness is not defined directly on A(·) for the ease of theorem statements.
In view of (3.2), Assumption 2 is implied by the smoothness on Ai under extra
regularity conditions. For a generic matrix M (t) parameterized by t ∈ (0, 1), we
let Ṁ (t) and M̈ (t) be the first two element-wise derivatives of M (t) w.r.t. t.
Lemma 3.1. Assume that supt∈[0,1] |A(t)|`1 < 1. Then we have
|Σ̇(t)|∞ ≤

2|Ȧ(t)|∞ · |A(t)|`1 · |Σ(t)|`1
1 − |A(t)|2`1

(3.4)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

8

and
|Σ̈(t)|∞

≤

8|Ȧ(t)|∞ · |A(t)|2`1 · |Ȧ(t)|`1 · |Σ(t)|`1
(1 − |A(t)|2`1 )2
+

(3.5)

4 max{|Ȧ(t)|∞ , |Ä(t)|∞ } · max{|A(t)|`1 , |Ȧ(t)|`1 } · |Σ(t)|`1
.
1 − |A(t)|2`1

By Lemma 3.1, if supt∈[0,1] |A(t)|`1 < 1 and Ajk (·) ∈ C 2 ([0, 1]) such that
supt∈[0,1] |Ȧ(t)|`1 ∨ |Σ(t)|`1 ≤ K for some constant K, then Σjk (·) is also a
C 2 ([0, 1]) function and therefore Assumption 2 is fulfilled.
Assumption 3 specifies the tail probability of the innovations ei . In [21], ei ’s
follow iid N (0, Ψ) for some error covariance matrix Ψ. A simple transformation by Ψ−1/2 will reduce to the case that ei has iid components with the
standard normal distribution, a special case of Assumption 3, which covers the
sub-Gaussian innovations.
Theorem 3.2. Let ρ̄ = supi≥0 ρ(Ai ). Fix an i ∈ [nbn + 1, n(1 − bn ) + 1].
Suppose that Assumption 1,2,3 are satisfied, ρ̄ < 1 and Ai ∈ Gα (s, Md ). Then,
−1
with probability at least
p1−2d , we have the estimator Âi with tuning parameter
2
τ ≥ C(1 + Md )(bn + (log d)/(nbn )) obeys
|Âi − Ai |∞

≤

ρ(Âi − Ai ) ≤
−1

d

|Âi −

Ai |2F

≤

2τ |Σ−1
i−1,0 |`1 ,

(3.6)

1−α
C(α)s(τ |Σ−1
,
i−1,0 |`1 )

(3.7)

2−α
C(α)s(τ |Σ−1
.
i−1,0 |`1 )

(3.8)

From Theorem 3.2, the bandwidth bn  ((log d)/n)−1/5 gives the optimal rate
of convergence and the resulting tuning parameter τ ≥ C(1+Md )((log d)/n)−2/5 .
When d is fixed, it is known that the optimal bandwidth in the nonparametric
kernel density estimation is n−1/5 for twice continuously differentiable functions
under the mean integrated square error. So in the high-dimensional context, the
dimension only has a logarithm impact on the choice of optimal bandwidth.
3.3. Pattern recovery
We also study the recovery of time-varying patterns using the estimator (2.5).
Let Si = supp(Ai ) be the nonzero positions of Ai . If the nonzero entries of Ai are
small enough, then it is impossible to accurately distinguish the small nonzeros
and the zeros. Therefore, the best hope we can expect is that the nonzero entries
of Ai with large magnitudes can be well separated from the zeros in Ai . Let
u] = 2τ |Σ−1
i−1,0 |`1 ,

(3.9)

where τ is determined in Theorem 3.2. We use the thresholded version of (2.5)
as an estimator of Si
Ŝi = Tu] (Âi ) = {1(|Âi,jk | > u] )}dj,k=1 .

(3.10)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

9

By Theorem 3.2, the maximal fluctuation |Âi − Ai |∞ is controlled by u] with
probability at least 1 − 2d−1 . So if we apply a thresholding procedure for Âi at
the level u] , then we expect that the zeros and non-zeros with magnitudes larger
than 2u] in Ai can be identified by the thresholded support of Âi . Precisely, we
have the instantaneous partial recovery consistency.
Theorem 3.3. Assume that Assumption 1,2,3 are satisfied. Then, for any fixed
time points i ∈ [nbn + 1, n(1 − bn ) + 1], P(Ŝi ⊂ Si ) ≥ 1 − 2d−1 and P({(m, k) :
|Ai,mk | > 2u] } ⊂ Ŝi ) ≥ 1 − 2d−1 .
Theorem 3.3 states that, with high probability, the zeros in Ai can be identified and the nonzero entries in Ai with strong signal strength above 2u] can
be recovered by Ŝi . Therefore, the false positives (type I error) of the estimator (3.10) are asymptotically controlled; see Theorem 3.5 for precise statement.
However, Theorem 3.3 does not provide too much information regarding the false
negatives since there is no characterization of the signal strength in (0, 2u] ).
Let β > 0 and u0 ∈ (0, 1). We introduce the following d × d matrix class
Gα,β (s, Md , Ld )


|{(m, k) : |Θmk | ∈ (0, u)}|
β
≤ Ld u , ∀0 < u < u0 . (3.11)
=Gα (s, Md ) ∩ Θ :
|supp(Θ)|
For A ∈ Gα,β (s, Md , Ld ), the parameters β and Ld control the proportion of
small entries in the support of A. If β is large and Ld grows slowly, then the
fraction of weak signals in A is small and therefore the false negatives (type II
error) can also be well controlled. Below, we shall give such an example.
Example 3.1 (A spatial design). Let 0 < r < 1. Consider a d × d symmetric
matrix A = (Amk )d×d that is generated by the covariance function of a spatial process Z1 , Z2 , ..., Zd , which is a random vector that is observed at sites
h◦1 , ..., h◦d ∈ D in some spatial domain D. Assume that the covariance between
Zm and Zk , m, k = 1, · · · , d satisfies
(
1−r
f (π(h◦m , h◦k )) π(h◦m , h◦k ) < d 2
Amk := cov(Zm , Zk ) =
,
(3.12)
1−r
0
π(h◦m , h◦k ) ≥ d 2
where π(h◦m , h◦k ) is the distance between the sites h◦m and h◦k
π(h◦m , h◦k ) =

|m − k|
dr

(3.13)

and f is a real-value covariance function. Here, we consider the rational quadratic
covariance function [11, 36]
f (π(h◦m , h◦k )) =

1
,
(1 + π(h◦m , h◦k )2 )γ

γ > 1.

(3.14)

In this example, Amk = 0 if |m − k| ≥ d/2, so A is a banded matrix.
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

10

Lemma 3.4. For A in Example 3.1, then there exists a large enough constant
C > 0 depending only on r, α, γ such that A ∈ Gα,β (s, Md , Ld ), where

if 2γα = 1
 Cdr log d
Cdr
if 2γα > 1
s=
,
(3.15)

Cdr d(1−r)(1−2γα) if 0 < 2γα < 1

if γ = 1/2
 Cdr log d
Cdr
if γ > 1/2
Md =
,
(3.16)

Cdr d(1−r)(1−2γ) if 0 < γ < 1/2
and Ld = Cd2(1−r)γβ .
For any fixed distance parameter r ∈ (0, 1), the weak signal parameter Ld has
a natural dependence on γ. If γ is smaller, then the covariance function f decays
to zero slower and there is a less fraction of weak signals in A. This can allow Ld
to grow slowly in d. Note that class Gα,β (s, Md , Ld ) is much less stringent than
the widely used condition for support recovery and model selection in literature,
which requires that the minimal nonzero signal strength is uniformly bounded
away from zero [28]. To quantify the error in the pattern recovery, we use the
following two error rate measures.
Definition 3.2. The false positive rate (FPR) and false negative rate (FNR)
of Ŝi are defined as
FPRi =

|Ŝi ∩ Sic |
,
|Sic |

FNRi =

|Ŝic ∩ Si |
.
|Si |

(3.17)

By convention, if Sic = ∅, then FPRi = 0; if Si = ∅, then FNRi = 0.
If FPRi = FNRi = 0 with probability tending to one, which is a very strong
requirement, then we have the pattern recovery consistency P(Ŝi = Si ) → 1.
Below, we show that both FPR and FNR are asymptotically controlled in the
presence of weak signals.
Theorem 3.5. Assume that Assumption 1,2,3 are satisfied. Fix an i ∈ [nbn +
1, n(1 − bn ) + 1]. If Ai ∈ Gα (s, Md ), then we have P(FPRi = 0) ≥ 1 − 2d−1 .
If in addtiion Ai,m ∈ Gα,β (s, Md , Ld ), then we also have P(FNRi ≤ 2β Ld uβ] ) ≥
1 − 2d−1 .
Since u] = o(1), the FNR vanishes with probability tending to one if Ld =
o(u−β
] ).
4. Simulation studies
In this section, we present some numerical results on simulated datasets. We
compare the following five methods:
(i) Our TV-VAR estimator (2.5).
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

(ii)
(iii)
(iv)
(v)

The
The
The
The

11

stationary vector autoregressive model (Stat. VAR) [21].
time-varying lasso method [23].
time-varying ridge method [20].
time-varying maximum likelihood estimator (MLE).

Methods (iii), (iv) and (v) are extensions of the lasso [23], ridge regression [20]
and MLE to the time-varying setting by kernel smoothing. We include (ii), i.e.
the stationary VAR model, as the baseline method that ignores the dynamic
features in the transition matrices. (iii) is a competitor of (i). We solve (iii) with
the FISTA [6] and we solve (i) and (ii) with the simplex algorithm [26].
4.1. Data generation
We consider different setups of n = 100 and d = 20, 30, 40 and 50. For each setup
(n, d), the data are generated by the following procedure. First, the baseline coefficient matrices A01 and A02 are generated by using the sugm.generator() in
flare R package [24]. We consider four graph structures defined in flare: hub,
cluster, band and random for A01 and A02 . Examples of these four structures
are shown in Appendix B. Then we normalize ρ(A01 ) = 0.2, ρ(A02 ) = 1 and
smoothly interpolate on the intermediate values
Ai = (1 − ti )4 A01 + t2i A02 ,

ti = i/n for i = 1, · · · , n.

Following [21], we specify the innovation covariance matrix Ψ = Σ − A01 ΣA>
01
where Σ = Id . [21] showed that the choice of Σ does not affect the numeric
performance significantly. In our simulation studies, we use the Epanechnikov
kernel K(v) = 0.75(1 − v 2 )I(|v| ≤ 1) and fix the model order k = 1. The
bandwidth bn is set to be bn = 0.8n−1/5 for (i), (iii) and (iv).
4.2. Tuning parameter selection
For the tuning parameter selection in (i)–(iv), we propose a data-driven procedure that minimizes the one-step-ahead prediction errors as follows.
1. Choose a grid for the tuning parameter (say τ ) and the number n1 of
training data points.
2. For each τ , perform the one-step-ahead prediction on the testing set by
estimating At with Ât (τ ) and then predict Xt by X̂t = Ât (τ )Xt−1 , where
t = n1 +1, ..., n. Then calculate the prediction error at time t as Errt (τ ) =
||Xt − Ât (τ )Xt−1 ||2 .
3. Calculate the average errors over the last n − n1 time points
Err(τ ) =

1
n − n1

n
X

Errt (τ )

t=n1 +1

4. Select τ̂ that minimizes Err(τ ).
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

12

4.3. Comparison results
4.3.1. Estimation errors
For each setup, we run 100 simulations and for each simulation indexed by nn we
estimate the transition matrix at each time point indexed by t, where t = 1, .., n.
Then, we calculate the error errt,nn = |Ât,nn − At |ρ for t = 1, ..., n and nn =
1, ..., 100. We also report errors under `F , `1 and `∞ norms. Next, P
the averaged
b
1
errors over different time points are aggregated as errnn = b−a
t=a errt,nn
where a = [nbn ] + 1 and b = [n(1 − bn )] − 1. Finally, the averaged errors are
P100
1
reported over 100 simulations as err = 100
nn=1 errnn . The results are shown
in Table 1, 2, 3 and 4.
From Table 1 to Table 4, larger d often results in larger errors. In general, the
unstructured MLE performs the worst almost under all matrix norms. Although
the ridge method shrinks the coefficients in the transition matrix to zeros, the
values of those coefficients are not exactly zeros. Thus, the estimated transition
matrices by the ridge method are not sparse which lead to higher estimation
errors compared with the TV-VAR, time-varying lasso and stationary VAR. Stationary VAR always perform worse than TV-VAR and time-varying lasso, since
it cannot capture the dynamic structure. The time-varying lasso performs better
than any other methods except TV-VAR. The proposed TV-VAR performs the
best almost under all matrix norms.
4.3.2. Pattern recovery
For the TV-VAR, we also report the FPRt,l,nn and FNRt,l,nn in (3.17), where
l = 1, .., 30 which are the indexes of 30 possible tuning
P100parameters
Pb from 0.001 to
1
1
0.45. Then, we calculate the averaged FPRl = 100
{
nn=1 b−a
t=a FPRt,l,nn }
P100
Pb
1
1
and FNRl = 100 nn=1 { b−a t=a FNRt,l,nn }. Following [9], we set the u] to
10−3 , which is considered to be numerical nonzero. The ROC curves for all
possible values of the sparsity control parameter are plotted in Fig 6a, Fig 6b,
Fig 6c and Fig 6d. Based on the ROC curves, the TV-VAR method has better
discrimination power for band or random patterns than hub or cluster patterns.
We also calculate u] by using the true value of Σ in (3.9) and the resulting
ROC curves are shown in Appendix C. Those ROC curves are similar to those
based on u] = 10−3 .
5. Real data analysis
5.1. Finance data: stock prices
In this section, we compare the aforementioned estimators on a real financial
dataset. The dataset is from Yahoo! Finance (finance.yahoo.com). The data
matrix contains daily closing prices of 452 stocks that are consistently in the

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

13

Table 1
Comparison of estimation errors under different setups. The standard deviations are shown
in parentheses. Here `ρ , `F , `1 and `∞ are the spectral, Frobenius, `1 and `∞ matrix norm
resp. The pattern of transition matrix is ‘hub’.
d=20

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.407
(0.083)
0.395
(0.088)
0.329
(0.038)
0.239
(0.011)

1.102
(0.248)
1.747
(0.516)
0.778
(0.147)
0.292
(0.032)

0.453
(0.068)
0.446
(0.076)
0.348
(0.032)
0.228
(0.011)

1.530
(0.075)
1.532
(0.068)
0.604
(0.022)
0.342
(0.008)

2.770
(0.205)
2.762
(0.255)
1.157
(0.08)
0.568
(0.02)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.643
(0.125)
0.694
(0.135)
0.430
(0.046)
0.245
(0.013)

1.154
(0.295)
2.197
(0.652)
0.846
(0.141)
0.288
(0.026)

0.695
(0.111)
0.732
(0.109)
0.443
(0.04)
0.244
(0.014)

2.116
(0.087)
2.102
(0.073)
0.710
(0.024)
0.393
(0.006)

4.350
(0.302)
4.302
(0.302)
1.607
(0.112)
0.729
(0.02)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.711
(0.103)
0.812
(0.144)
0.460
(0.041)
0.253
(0.011)

1.221
(0.256)
2.868
(0.775)
0.958
(0.124)
0.296
(0.02)

0.800
(0.095)
0.872
(0.126)
0.486
(0.035)
0.259
(0.011)

2.594
(0.079)
2.586
(0.082)
0.793
(0.024)
0.429
(0.006)

6.183
(0.426)
6.198
(0.406)
2.168
(0.142)
0.907
(0.021)

d=30

`∞
`1
`ρ
`F

d=40

`∞
`1
`ρ
`F

d=50

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.827
(0.124)
1.019
(0.19)
0.504
(0.045)
0.266
(0.014)

1.158
(0.219)
3.574
(1.021)
1.052
(0.157)
0.297
(0.019)

0.944
(0.102)
1.072
(0.135)
0.531
(0.032)
0.278
(0.012)

3.043
(0.095)
3.055
(0.098)
0.863
(0.021)
0.454
(0.006)

8.478
(0.633)
8.558
(0.65)
2.874
(0.224)
1.108
(0.027)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

14

Table 2
Comparison of estimation errors under different setups. The standard deviations are shown
in parentheses. Here `ρ , `F , `1 and `∞ are the spectral, Frobenius, `1 and `∞ matrix norm
resp. The pattern of transition matrix is ‘cluster’.
d=20

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.400
(0.073)
0.383
(0.074)
0.324
(0.032)
0.239
(0.011)

1.093
(0.300)
1.695
(0.506)
0.774
(0.143)
0.291
(0.03)

0.465
(0.084)
0.448
(0.077)
0.352
(0.033)
0.228
(0.011)

1.527
(0.081)
1.527
(0.083)
0.601
(0.025)
0.341
(0.008)

2.736
(0.202)
2.731
(0.251)
1.133
(0.084)
0.562
(0.019)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.557
(0.096)
0.589
(0.132)
0.399
(0.045)
0.241
(0.009)

1.127
(0.285)
2.236
(0.659)
0.848
(0.135)
0.289
(0.027)

0.631
(0.097)
0.664
(0.110)
0.424
(0.038)
0.239
(0.010)

2.095
(0.083)
2.120
(0.093)
0.710
(0.024)
0.393
(0.006)

4.305
(0.308)
4.343
(0.297)
1.607
(0.112)
0.729
(0.020)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.705
(0.115)
0.803
(0.153)
0.460
(0.042)
0.253
(0.011)

1.139
(0.225)
2.986
(0.874)
0.974
(0.15)
0.295
(0.022)

0.783
(0.108)
0.847
(0.133)
0.479
(0.037)
0.258
(0.012)

2.612
(0.106)
2.596
(0.098)
0.794
(0.022)
0.428
(0.006)

6.179
(0.44)
6.171
(0.452)
2.158
(0.156)
0.906
(0.022)

d=30

`∞
`1
`ρ
`F

d=40

`∞
`1
`ρ
`F

d=50

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.841
(0.121)
1.023
(0.177)
0.504
(0.04)
0.267
(0.013)

1.116
(0.198)
3.384
(1.055)
1.015
(0.162)
0.292
(0.019)

0.964
(0.108)
1.070
(0.152)
0.531
(0.036)
0.278
(0.013)

3.052
(0.106)
3.043
(0.09)
0.863
(0.021)
0.454
(0.006)

8.593
(0.703)
8.498
(0.558)
2.874
(0.224)
1.108
(0.027)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

15

Table 3
Comparison of estimation errors under different setups. The standard deviations are shown
in parentheses. Here `ρ , `F , `1 and `∞ are the spectral, Frobenius, `1 and `∞ matrix norm
resp. The pattern of transition matrix is ‘band’.
d=20

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.402
(0.072)
0.385
(0.075)
0.326
(0.033)
0.237
(0.011)

1.048
(0.251)
1.645
(0.476)
0.758
(0.137)
0.289
(0.031)

0.453
(0.068)
0.437
(0.077)
0.347
(0.032)
0.226
(0.01)

1.532
(0.08)
1.531
(0.068)
0.600
(0.021)
0.341
(0.008)

2.738
(0.188)
2.707
(0.183)
1.132
(0.072)
0.563
(0.019)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.574
(0.094)
0.590
(0.131)
0.401
(0.041)
0.240
(0.009)

1.139
(0.319)
2.306
(0.727)
0.863
(0.154)
0.290
(0.028)

0.637
(0.101)
0.654
(0.123)
0.422
(0.038)
0.239
(0.011)

1.537
(0.058)
1.537
(0.057)
0.537
(0.019)
0.315
(0.005)

4.316
(0.314)
4.312
(0.323)
1.607
(0.112)
0.729
(0.02)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.800
(0.134)
0.925
(0.155)
0.489
(0.041)
0.262
(0.013)

1.170
(0.253)
2.772
(0.825)
0.938
(0.138)
0.293
(0.022)

0.886
(0.113)
0.964
(0.127)
0.510
(0.035)
0.270
(0.012)

2.590
(0.098)
2.597
(0.09)
0.793
(0.024)
0.429
(0.006)

6.198
(0.429)
6.221
(0.457)
2.168
(0.142)
0.907
(0.021)

d=30

`∞
`1
`ρ
`F

d=40

`∞
`1
`ρ
`F

d=50

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.839
(0.125)
1.011
(0.175)
0.506
(0.045)
0.267
(0.014)

1.143
(0.222)
3.673
(1.212)
1.047
(0.176)
0.295
(0.02)

0.992
(0.118)
1.105
(0.126)
0.541
(0.034)
0.280
(0.013)

3.056
(0.100)
3.047
(0.099)
0.862
(0.021)
0.454
(0.005)

8.550
(0.579)
8.635
(0.678)
2.872
(0.223)
1.108
(0.027)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

16

Table 4
Comparison of estimation errors under different setups. The standard deviations are shown
in parentheses. Here `ρ , `F , `1 and `∞ are the spectral, Frobenius, `1 and `∞ matrix norm
resp. The pattern of transition matrix is ‘random’.
d=20

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.397
(0.078)
0.379
(0.080)
0.323
(0.036)
0.237
(0.013)

1.060
(0.244)
1.661
(0.574)
0.757
(0.140)
0.286
(0.031)

0.447
(0.07)
0.433
(0.081)
0.344
(0.033)
0.227
(0.012)

1.538
(0.071)
1.532
(0.077)
0.600
(0.021)
0.341
(0.008)

2.704
(0.188)
2.762
(0.235)
1.132
(0.072)
0.563
(0.019)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.567
(0.103)
0.590
(0.132)
0.398
(0.044)
0.240
(0.01)

1.138
(0.264)
2.280
(0.724)
0.865
(0.146)
0.290
(0.026)

0.641
(0.102)
0.658
(0.119)
0.423
(0.041)
0.239
(0.012)

1.541
(0.066)
1.543
(0.062)
0.537
(0.019)
0.315
(0.005)

4.294
(0.292)
4.292
(0.268)
1.607
(0.112)
0.729
(0.020)

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.713
(0.098)
0.813
(0.136)
0.458
(0.036)
0.252
(0.01)

1.210
(0.268)
2.904
(0.851)
0.951
(0.142)
0.296
(0.021)

0.810
(0.104)
0.869
(0.128)
0.482
(0.034)
0.260
(0.012)

2.595
(0.096)
2.600
(0.09)
0.793
(0.024)
0.429
(0.006)

6.148
(0.413)
6.205
(0.432)
2.168
(0.142)
0.907
(0.021)

d=30

`∞
`1
`ρ
`F

d=40

`∞
`1
`ρ
`F

d=50

`∞
`1
`ρ
`F

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.830
(0.109)
1.023
(0.187)
0.504
(0.039)
0.267
(0.013)

1.207
(0.234)
3.529
(1.214)
1.034
(0.180)
0.298
(0.021)

1.002
(0.120)
1.129
(0.151)
0.543
(0.033)
0.282
(0.013)

3.060
(0.108)
3.061
(0.092)
0.865
(0.022)
0.455
(0.005)

8.452
(0.553)
8.457
(0.551)
2.861
(0.198)
1.107
(0.027)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

17

Fig 1: ROC curves under different settings
(a) Pattern is hub

(b) Pattern is cluster

d=20

0.4

0.6

0.6

d=30

d=40

d=40
d=50

0.0

0.0

0.2

d=20

d=30

d=50
0.0

0.4

True Positive Rate (Sensitivity)

0.2

0.6
0.4
0.2

True Positive Rate (Sensitivity)

0.8

cluster

0.8

hub

0.8

0.0

0.2

False Positive Rate (1−Specificity)

0.6

0.8

False Positive Rate (1−Specificity)

(c) Pattern is band

(d) Pattern is random

0.0

0.2

0.4

0.6

False Positive Rate (1−Specificity)

0.8

0.6
0.4

d=20

d=30

d=30

d=40

d=40

d=50
0.0

0.2

d=20

d=50

0.0

0.2

0.4

0.6

True Positive Rate (Sensitivity)

0.8

random

0.8

band

True Positive Rate (Sensitivity)

0.4

0.0

0.2

0.4

0.6

0.8

False Positive Rate (1−Specificity)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

18

S&P 500 index between January 1st, 2003 and January 1st, 2008. We choose
such time range to avoid the effects of the two financial crises in the year of 2001
and 2008, which could make the stock prices non-smooth with sharp drops and
pick-ups. In total, there are 1,258 time points. We first standardize the data to
zero mean and unit variance and detrend the data. We then fit the detrended
data for those stocks which are most smoothed (without obvious change points)
using AR(1) model and perform the Ljung-Box Tests by using Box.test() in
R. We keep those stocks with nonzero coefficients at significant level 0.05. 30
stocks are finally selected and 10 of them are shown in Table 9. The ten selected
stocks in Table 9 are from six sectors including: Consumer Staples, Consumer
Discretionary, Industrials, Financials, Utilities and Energy. The resulting data
matrix is denoted as X̄ ∈ R1258×30 . We use the Priestley-Subba Rao (PSR)
stationarity test (such as stationarity() in R package fractal) to some of
the stocks such as Kellogg Co. and Target Corp. to find that these time series
are not stationary at significant level 0.05. Thus, it is inappropriate to model
the data with the stationary VAR model in [21]. Finally we fit the data X̄
into the sparse TV-VAR model with order k = 1, stationary VAR [21] with
order k = 1, time-varying lasso method [23], time-varying ridge method [20]
and time-varying MLE.
In order to compare the performance of the five methods, we consider the
one-step-ahead prediction for X̄t , t = 1159, ..., 1258, by using X̄Jt ,∗ , where Jt =
{j : t−1158 ≤ j ≤ t−1} as the training set. We estimate the transition matrices
Ât (λ) where λ is the regularization parameter such as the τ in the TV-VAR or
the shrinkage parameter in the ridge method. We use the Epanechnikov kernel
K(v) = 0.75(1 − v 2 )I(|v| ≤ 1) with bandwidth bn = 0.3, which means that only
the last 30% data in the training set are used for prediction for the TV-VAR,
time-varying lasso method, time-varying ridge method and time-varying MLE.
Since stationary VAR [21] model uses all training data, in order to make fair
comparisons, we only treat the last 30% time points (347 observations) in X̄Jt ,∗
as the training set for the stationary VAR. The averaged prediction error for a
specific λ is measured by:
Err(λ) =

1258
1 X
kX̄t,∗ − Ât (λ)X̄t−1,∗ k2 .
100 t=1159

The smallest averaged one-step-ahead prediction errors minλ∗ Err(λ) for five
methods were shown in Table 5. In Figure 2, we show an example for the predicted closing prices from the five methods and the detrended true closing prices
of Target Corp. Clearly, the methods with time-varying structures perform similarly but outperform stationary VAR.
One advantage of the TV-VAR compared with the stationary VAR is that
it can capture the time-varying data structures. If we treat a transition matrix
as an adjacency matrix for an undirected weighted graph, we can visualize the
graph structures at different time points. For instance, Fig 3a, 3b, 3c and 3d are
four graphical representations of the TV-VAR, which clearly show some timeevolving patterns. In the graphs, different vertices represent different stocks and
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

19

the edges represent the cross-sectional dependence between stocks. Bolder lines
indicate stronger dependence and the stocks in the same sector are in the same
color. In order to illustrate the time-varying structure clearly, we only consider
the 10 stocks in Table 9. From these four figures, we observe that stocks in the
same and closely related sectors are often connected. For examples, CME Group
Inc. and Hartford Financial Svc. GP from Financials sector show a dynamic
dependence structure, and Boeing Company in Industrials shows the consistent
correlation with Exxon Mobil Corp. in Energy sector.
Table 5
The prediction errors for the five methods on the stock price data. The standard deviations
are shown in parentheses.

Err(λ∗)
Standard Deviation

TV-VAR

Stat. VAR

Lasso

Ridge

MLE

0.4822
(0.1751)

2.2824
(0.5169)

0.4580
(0.1810)

0.4843
(0.1747)

0.4902
(0.1799)

−0.5

Detrended Close Price

0.0

0.5

Target Corp.

−1.0

True value
TV−VAR
Lasso
Ridge
MLE

−1.5

Stat. VAR

0

20

40

60

80

100

Time

Fig 2: Comparison of the predicted closing prices and the true detrended closing
prices of Target Corp.

5.2. Economic data: exchange rates
We also apply the exchange rates data on international currencies and study how
the exchange rates evolve over time. We use the exchange rates data from the
Federal Reserve Bank at St. Louis and choose the date range from January 1st,
2003 to January 1st, 2008 to avoid the financial crises’ effects on exchange rates.
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

20

Fig 3: The graphs based on adjacency matrices estimated by TV-VAR
(a) Adjacency matrix Â1183

CME (FIN)

(b) Adjacency matrix Â1208

Boeing (IND)

CME (FIN)

Target (CD)

Prudential (FIN)

Kellogg (CS)

Edison (UTI)

Exxon (ENE)

Lockheed (IND)

PepsiCo (CS)

CME (FIN)

CME (FIN)

Exxon (ENE)

PepsiCo (CS)

Hartford (FIN)

Hartford (FIN)

(d) Adjacency matrix Â1258

Kellogg (CS)

Lockheed (IND)

Exxon (ENE)

Lockheed (IND)

PepsiCo (CS)

Target (CD)

Edison (UTI)

Kellogg (CS)

Edison (UTI)

Boeing (IND)

Prudential (FIN)

Target (CD)

Prudential (FIN)

Hartford (FIN)

(c) Adjacency matrix Â1233

Boeing (IND)

Boeing (IND)

Target (CD)

Prudential (FIN)

Kellogg (CS)

Edison (UTI)

Exxon (ENE)

Lockheed (IND)

PepsiCo (CS)

Hartford (FIN)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

21

In total, there are 1,303 time points. We choose 15 major currencies from the
continents including Europe, Australia, North America, Asian, South America
and we normalize each currency by the exchange rate of the U.S. dollar shown in
Table 10. As in Section 5.1, we standardize and detrend the data. The exchange
rates for all currencies to the U.S. dollars are smoothed and in addition the
PSR test rejects the stationarity hypothesis of the exchange rate Euro/U.S. at
the significant level 0.05. Thus, stationary VAR is also not appropriate in this
application, which is also supported by our empirical findings in Table 6 and
Figure 4. We compare the performance of the five methods by the one-stepahead prediction errors for the last 50 data and use the same prediction metric
in Section 5.1. The bandwidth here is bn = 0.3. We do not include the prediction
curve of the stationary VAR in Figure 4 since its errors are too large. Table 6
and Figure 4 show that methods with time-varying structure performs similarly
but much better than the stationary VAR.
Table 6
The prediction errors for the five methods on the exchange rates data. The standard
deviations are shown in parentheses.

Err(λ∗)
Standard Deviation

TV-VAR

Stat. Var

Lasso

Ridge

MLE

0.2990
(0.0924)

8.7637
(5.2835)

0.2554
(0.1037)

0.2634
(0.1032)

0.2638
(0.1032)

0.0
−0.2
−0.4

Detrended Exchange Rate

0.2

0.4

0.6

Exchange Rate: U.K./U.S.

True value
TV−VAR

−0.6

Lasso
Ridge

−0.8

MLE

0

10

20

30

40

50

Time

Fig 4: Compare the predicted and true detrended exchange rate of U.K./U.S.

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

22

6. Proofs
Proof of Lemma 3.1. Since x̃m (t) is a stationary VAR, we have Σ̃(t) = Id×d +
A(t)Σ̃(t)A(t)> . By the extension, Σ(t) = Σ̃(t) for all t ∈ (0, 1). Applying the
chain rule of differentiation w.r.t. t, we get
Σ̇(t) = Ȧ(t)Σ(t)A(t)> + A(t)Σ̇(t)A(t)> + A(t)Σ(t)Ȧ(t)> .
By Lemma 6.1 and the triangle inequality,
|Σ̇(t)|∞
≤|Ȧ(t)|∞ |Σ(t)|`∞ |A(t)|`1 + |A(t)|`1 |Σ̇(t)|∞ |A(t)|`1 + |A(t)|`1 |Σ(t)|`1 |Ȧ(t)|∞ .
Then, (3.4) follows from the assumption that supt∈[0,1] |A(t)|`1 < 1 and that
Σ(t) is symmetric. (3.5) follows from the differentiation of Σ̇(t) w.r.t. t.
Lemma 6.1. Let A, B, C be matrices of compatible dimension for the product
ABC. Then
|ABC|∞

≤

|A|`1 |B|∞ |C|`∞ ,

|ABC|∞

≤

|A|`1 |B|`1 |C|∞ ,

|ABC|∞

≤

|A|∞ |B|`∞ |C|`∞ .

Proof of Lemma 6.1. The first and second inequalities follow from
|ABC|∞

=

max
i,j

X

Aik Bkm Cmj ≤ max
i,j

k,m

X

|Aik |

k

X

|Bkm Cmj |

m

!
≤
≤

max
i

X
k

|Aik | max
j,k

X

|Bkm Cmj |

m

|A|`1 min {|B|∞ |C|`∞ , |B|`1 |C|∞ } .

The third inequality follows from the second one by considering (ABC)> =
C > B > A> and |A> |`1 = |A|`∞ .
Recall the TV-VAR(1) model xi = Ai xi−1 + ei . The following key lemma
presents a large deviation bound for the marginal and lag-one autocovariance
matrices with sub-Gaussian innovations.
Lemma 6.2. Suppose ρ̄ = supi≥0 ρ(Ai ) < 1 and
sup |Σ̇m,jk (t)| ∨ |Σ̈m,jk (t)| ≤ C0 , m = 0, 1,
t∈[0,1]

for some absolute constant C0 < ∞. If ei = (ei1 , · · · , eid ) has iid sub-Gaussian
components, i.e. keij kq ≤ C1 q 1/2 for all q ≥ 1, then, for any fixed t ∈ [bn , 1−bn ],
we have with probability at least 1 − 2d−1 that
!
r
log d
2
|Σ̂t,0 − Σt,0 |∞ ∨ |Σ̂t,1 − Σt,1 |∞ ≤ C bn +
,
(6.1)
nbn
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

23

where the constant C = C(ρ̄, C0 , C1 ) does not depend on n, bn and d.
Proof of Lemma 6.2. First, consider the marginal covariance matrix m = 0.
For each j, k = 1, · · · , p, we write
Σ̂0,jk (t) − Σ0,jk (t) = {Σ̂0,jk (t) − E[Σ̂0,jk (t)]} + {E[Σ̂0,jk (t)] − Σ0,jk (t)} = I + II.
For |ti − t| ≤ bn , by the assumption supt∈[0,1] |Σ̈0,jk (t)| ≤ C0 , Taylor’s expansion
yields that Σ0,jk (ti ) = Σ0,jk (t) + Σ̇0,jk (t)(ti − t) + O(b2n ) for all t ∈ [bn , 1 − bn ].
Note that for ` = 0, 1, 2, we have
sup

n
X


K

t∈[bn ,1−bn ] i=1

ti − t
bn



(ti − t)` − nbn

Z

1

K(v)v ` dv = O(1).

−1

Then, the bias part is controlled by
II

=

n
X

w(t, i)Σ0,jk (ti ) − Σ0,jk (t)

i=1

= Σ̇0,jk (t)

n
X

w(t, i)(ti − t) + O(b2n ) = O((nbn )−1 + b2n ).

(6.2)

i=1

Now, we deal with the stochastic
Pnpart I. Let Yi,jk = Xij Xik − E(Xij Xik ).
Clearly, E(Yi,jk ) = 0 and I =
i=1 w(t, i)Yi,jk . Let Bi,0 = Id and Bi,m =
m
Ai · · · Ai−m+1 for m ≥ 1. Then, supP
and we have the movingi ρ(Bi,m ) ≤ ρ̄
∞
>
>
average (MA) representation of xi = m=0 Bi,m ei−m . Let ξ = (e>
n , en−1 , · · · ) ,
Wt = diag(w(t, n), · · · , w(t, 1)) be n × n diagonal matrix, and
B̃ (j) =

Bn,0,j·

0


0


..

.
0

Bn,1,j·
Bn−1,0,j·
0
..
.

Bn,2,j·
Bn−1,1,j·
Bn−2,0,j·
..
.

···
···
···
..
.

Bn,n−1,j·
Bn−1,n−2,j·
Bn−2,n−3,j·
..
.

Bn−1,n−1,j·
Bn−2,n−2,j·
..
.

···
···
···
..
.

0

0

···

B1,0,j·

B1,1,j·

···

Bn,n,j·





.



Then, we can write I = ξ > (B̃ (j) )> Wt B̃ (k) ξ−tr((B̃ (j) )> Wt B̃ (k) ). By the HansonWright inequality [29, Theorem 1.1], there exists a constant C := C(C1 ) such
that for all x > 0


(∗) := P |ξ > (B̃ (j) )> Wt B̃ (k) ξ − tr((B̃ (j) )> Wt B̃ (k) )| ≥ x
#)
(
"
x
x2
,
.
≤
2 exp −C min
|(B̃ (j) )> Wt B̃ (k) |2F ρ((B̃ (j) )> Wt B̃ (k) )
Observe that (B̃ (j) )> Wt B̃ (k) has the same nonzero eigenvalues as the n × n
1/2
1/2
matrix Wt B̃ (k) (B̃ (j) )> Wt . Therefore, by the Cauchy-Schwartz inequality,
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

24

we get
1/2

|Wt

h
i
1/2
1/2
tr Wt B̃ (j) (B̃ (k) )> Wt B̃ (k) (B̃ (j) )> Wt
h
i
= tr (B̃ (k) )> Wt B̃ (k) (B̃ (j) )> Wt B̃ (j)

1/2 2
|F

B̃ (k) (B̃ (j) )> Wt

=

|(B̃ (k) )> Wt B̃ (k) |F · |(B̃ (j) )> Wt B̃ (j) |2F

≤

and by the definition of matrix spectral norm
1/2

ρ(Wt

1/2

B̃ (k) (B̃ (j) )> Wt

1/2

|) ≤ ρ(Wt
= ρ

1/2

1/2

B̃ (k) )ρ((B̃ (j) )> Wt

((B̃

(k) >

) Wt B̃

Therefore, (*) is bounded by
(
"
2 exp −C min

(k)

)ρ

1/2

)

((B̃ (j) )> Wt B̃ (j) ).

x2
1/2

maxj≤p |Wt

1/2 2
|F

B̃ (j) (B̃ (j) )> Wt

#)

x
1/2

maxj≤p ρ(Wt

,

1/2

B̃ (j) (B̃ (j) )> Wt

.
)

For i = 1, · · · , n and l = 0, · · · , n − i, let
(j)

γi,l =

∞
X

>
[w(t, i)w(t, i + l)]1/2 |Bi,m,j· Bi+l,m+l,j·
|.

m=0

By the Cauchy-Schwartz inequality and the spectral norm bound supi ρ(Bi,m ) ≤
ρ̄m , we have
(j)
γi,l

≤
≤

∞
X

p
X

m=0
∞
X

k=1

[w(t, i)w(t, i + l)]1/2
[w(t, i)w(t, i + l)]1/2

!1/2
2
Bi,m,jk

p
X

!1/2
2
Bi+l,m+l,jk

k=1

ρ̄m ρ̄m+l = [w(t, i)w(t, i + l)]1/2

m=0

ρ̄l
.
1 − ρ̄2

Then, it follows that
1/2

≤

≤

1/2

|Wt B̃ (j) (B̃ (j) )> Wt |2F
Pn
n−1
n−l
2
XX
ρ̄l
1/2
i=1 w(t, i)
+
2
[w(t,
i)w(t,
i
+
l)]
(1 − ρ̄2 )2
(1 − ρ̄2 )2
l=1 i=1
hP
i P

n
n−1 2l
2
2
l=0 ρ̄
i=1 w(t, i)

≤ C(ρ̄)

(1 − ρ̄2 )2
n
X
i=1

w(t, i)2 =

C(ρ̄)
nbn

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

25

and
1/2

ρ(Wt
"
≤

≤

1/2

B̃ (j) (B̃ (j) )> Wt

)

n−1
X ρ̄l
1
2
max
w(t,
i)
+
max [w(t, i)w(t, i + l)]1/2
2
1 − ρ̄ 1≤i≤n
1 − ρ̄2 1≤i≤n−l
l=1
hP
i
n−1 l
2
l=0 ρ̄ [max1≤i≤n w(t, i)]
C(ρ̄)
≤
.
2
1 − ρ̄
nbn

#

Therefore, we have that for any x > 0
sup

P(|

t∈[bn 1−bn ]

n
X

w(t, i)Yi,jk | ≥ x) ≤ 2 exp[−Cnbn min(x2 , x)].

(6.3)

i=1

Now, by (6.2) and using the union bound applied to (6.3), we obtain that there
exists a constant C which only depends on ρ̄, C0 , C1 such that
!
r
log d
2
|Σ̂t,0 − Σt,0 |∞ ≤ C bn +
nbn
holds with probability ≥ 1 − 2d−1 . Similar argument applied to m = 1 shows
that the lag-one autocovariance matrix obeys the same bound in (6.1).
Proof of Lemma 3.4. Denote l = d/2. For 0 < α < 1, note that
max

1≤m≤d

d
X

Aα
mk = 2

k=1

l−1
X

α
Aα
lk + All = 2

k=1

l−1
X

Aα
lk + 1.

k=1

Then, we have
l−1
X

Aα
lk =

k=1

l−1
X
k=1

=

1
(1 + |l − k|2 /d2r )γα

l−1
X
k=1

d2rγα
+ k 2 )γα

(d2r
l

Z
≤

1+
1

=

1 + dr

d2rγα
dx
+ x2 )γα

(d2r
Z

d1−r /2

d−r

≤

r

Z

Cd

1
dx
(1 + x2 )γα

d1−r /2

x−2γα dx,

1

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

26

where C > 0 is a constant depending only on r, γ, α. Then, it follows that

d
O(dr )
if 2γα > 1

X
α
O(dr log d)
if 2γα = 1
max
Amk =
.

1≤m≤d
k=1
O(dr+(1−r)(1−2γα) ) if 1 < 2γα < 1
So, A ∈ Gα− (s, Md ) for our choice of s and Md . Since A is symmetric, A ∈
Gα (s, Md ). Next, we show that A ∈ Gα,β (s, Md , Ld ). Observe that |supp(A)| =
(3d − 2)d/4. By construction, Amk ≤ u is equivalent to π(h◦m , h◦k ) ≥ f −1 (u) =
1
1
1
1
(u− γ − 1) 2 . Therefore, Amk ≤ u if and only if |m − k| ≥ dr (u− γ − 1) 2 . Since
1
1
dr (u− γ − 1) 2 is monotone decreasing with respect to u,
1
1
d
d
∃u ∈ (0, u0 ], s.t. − 2 < dr (u− γ − 1) 2 ≤ − 1.
¯
¯
2
2

So we only need to consider u ∈ [u, u0 ] because if u ∈ (0, u), then |{(m, k) : 0 <
¯
¯
Amk ≤ u}| = |{(m, k) : |m−k| ≥ d/2}| = 0. Therefore, if we choose Ld such that
Ld uβ & 1, then A ∈ Gα,β (s, Md , Ld ). That is, Ld = C u−β = Cd2(1−r)γβ .
¯
¯
p
Proof of Theorem 3.2. Put ε = C(b2n + (log d)/(nbn )). Recall that A>
i =
−1
Σ−1
, that
i−1,0 Σi−1,1 . By Lemma 6.2, we have, with probability at least 1 − 2d
|Σ̂t,0 − Σt,0 |∞ ≤ ε and |Σ̂t,1 − Σt,1 |∞ ≤ ε for any t ∈ [bn , 1 − bn ]. Then, with such
high probability, the true transition matrix Ai is feasible for the minimization
program (2.5) with τ ≥ (1 + Md )ε because
|Σ̂i−1,1 − Σ̂i−1,0 A>
i |∞

= |Σ̂i−1,1 − Σ̂i−1,0 Σ−1
i−1,0 Σi−1,1 |∞
≤

|Σ̂i−1,1 − Σi−1,1 |∞ + |(I − Σ̂i−1,0 Σ−1
i−1,0 )Σi−1,1 |∞

≤ ε + |Σ̂i−1,0 − Σi−1,0 |∞ |A>
i |`1 ≤ τ.
Hence, |Âi |1 ≤ |Ai |1 . But then
>
|Â>
i − Ai | ∞

=

−1
|Â>
i − Σi−1,0 Σi−1,1 |∞

=

>
|Σ−1
i−1,0 (Σi−1,0 Âi − Σi−1,1 )|∞

=

>
>
>
|Σ−1
i−1,0 (Σi−1,0 Âi − Σ̂i−1,0 Âi + Σ̂i−1,0 Âi − Σ̂i−1,1 + Σ̂i−1,1 − Σi−1,1 )|∞

≤

−1
>
>
|Σ−1
i−1,0 (Σ̂i−1,0 − Σ̂i−1,1 )Âi |∞ + |Σi−1,0 (Σ̂i−1,0 Âi − Σ̂i−1,1 )|∞

+|Σ−1
i−1,0 (Σ̂i−1,1 − Σi−1,1 )|∞
≤

|Σ−1
i−1,0 |`1 (|Σ̂i−1,0 − Σ̂i−1,1 |∞ |Âi |`1 + τ + ε)

≤

|Σ−1
i−1,0 |`1 (ε|Ai |`1 + τ + ε)

≤

2τ |Σ−1
i−1,0 |`1 .

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

27

Let u ≥ 0 and Tj = {m : |Ai,jm | ≥ u}, j = 1, · · · , d. Define D(u) =
Pd
maxj≤d m=1 (|Ai,jm | ∧ u). Then, by the triangle inequality and the equivalence between (2.5) and (2.6), we have uniformly in j = 1, · · · , d,
|[Âi − Ai ]j∗ |1

≤

|[Âi ]j,Tjc |1 + |[Ai ]j,Tjc |1 + |[Âi − Ai ]j,Tj |1

=

|[Ai ]j∗ |1 − |[Âi ]j,Tj |1 + |[Ai ]j,Tjc |1 + |[Âi − Ai ]j,Tj |1

≤

2|[Ai ]j,Tjc |1 + 2|[Âi − Ai ]j,Tj |1

≤

2|[Ai ]j,Tjc |1 + 4τ |Σ−1
i−1,0 |`1 |Tj |

≤

2D(u)(1 + 2τ |Σ−1
i−1,0 |`1 /u).

1−α
for Ai ∈ Gα (s, Md ), it follows that
Choose u = 2τ |Σ−1
i,0 |`1 . Since D(u) ≤ su
1−α
|[Âi − Ai ]j∗ |1 ≤ 4D(u) ≤ 4su1−α = C(α)s(|Σ−1
.
i,0 |`1 τ )

The analysis of the other constraint is similar. Now, (3.7) and (3.8) are immediate in view of ρ2 (M ) ≤ |M |`1 |M |`∞ and |M |2F ≤ d|M |∞ |M |`1 for any d × d
matrix M .
Proof of Theorem 3.3. Suppose (m, k) ∈
/ Si , i.e. Ai,mk = 0. On the event
G = {|Âi − Ai |∞ ≤ u] }, we have |Âi,mk | ≤ u] and (m, k) ∈
/ Ŝi . Therefore, by
Theorem 3.2, we have P(Ŝi ⊂ Si ) ≥ 1 − 2d−1 . On the other hand, suppose that
|Ai,mk | > 2u] , then on the event G, it follows from the triangle inequality that
|Ai,mk | − |Âi,mk | ≤ |Âi,mk − Ai,mk | ≤ u] . Therefore, |Âi,mk | > u] and we have
P({(m, k) : |Ai,mk | > 2u] } ⊂ Ŝi ) ≥ 1 − 2d−1 .
Proof of Theorem 3.5. If |Si | = 0, then the first claim is trivial. So we may
assume that Si 6= ∅. Suppose (m, k) ∈
/ Si , then on the event G = {|Âi −
Ai |∞ ≤ u] }, we have |Âi,mk | ≤ u] and (m, k) ∈
/ Ŝi . Therefore, Ŝi ∩ Sic = ∅
−1
on G and P(FPRi = 0) ≥ 1 − 2d
by Theorem 3.2. On the other hand,
let Ni (u) = {(m, k) : 0 < |Ai,mk | ≤ 2u}. For any u > 0, Ni (u) ⊂ Si and
Ŝic ∩Si ⊂ Ni (u)∪(Ŝic ∩Si ∩Ni (u)c ). Hence, |Ŝic ∩Si | ≤ |Ni (u)|+|Ŝic ∩Si ∩Ni (u)c |.
If (m, k) ∈ Ŝic and |Ai,mk | > 2u] , then Ŝic ∩Si ∩Ni (u] )c = ∅ on G. Choose u = u]
and then on the event G, we have
FNRi =

|Ŝic ∩ Si |
|Ni (u] )| |Ŝic ∩ Si ∩ Ni (u] )c |
≤
+
≤ 2β Ld uβ] ,
|Si |
|Si |
|Si |

from which the theorem follows.
Acknowledgements
We thank Fang Han for providing the stock data in Section 5.1. The highperformance computing work was done on the Illinois Campus Cluster Program
(ICCP) at the University of Illinois at Urbana-Champaign.
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

28

References
[1] Lian An and Jian Wang. Exchange rate pass-through: Evidence based
on vector autoregression with sign restrictions. Open Economies Review,
23(2):359–380, 2012.
[2] Andrew Ang and Monika Piazzesi. A no-arbitrage vector autoregression of
term structure dynamics with macroeconomic and latent variables. Journal
of Monetary Economics, 50:745–787, 2003.
[3] David Backus. The canadian–u.s. exchange rate: Evidence from a vector
autoregression. The Review of Economics and Statistics, 68(4):628–637,
1986.
[4] Richard T. Baillie, Robert E. Lippens, and Patrick C. McMahon. Testing rational expectations and efficiency in the foreign exchange market.
Econometrica, 51(3):553–563, 1983.
[5] Sumanta Basu and George Michailidis. Regularized estimation in sparse
high-dimensional time series models. Annals of Statistics, 43(4):1535–1567,
2015.
[6] 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.
[7] Peter J. Bickel and Elizaveta Levina. Covariance Regularization by Thresholding. The Annals of Statistics, 36(6):2577–2604, 2008.
[8] Peter J. Brockwell and Richard A. Davis. Time Series: Theory and Methods. Series in Statistics. Springer, second edition, 1991.
[9] Tony Cai, Weidong Liu, and Xi Luo. A Constrained `1 Minimization Approach to Sparse Precision Matrix Estimation. Journal of American Statistical Association, 106(494):594–607, 2011.
[10] John Y. Campbell and Robert J. Shiller. Stock prices, earnings, and expected dividends. Journal of Finance, 43(3):661–676, 1988.
[11] Xiaohui Chen, Mengyu Xu, and Wei Biao Wu. Covariance and precision
matrix estimation for high-dimensional time series. Annals of Statistics,
41(6):2994–3021, 2013.
[12] Xiaohui Chen, Mengyu Xu, and Wei Biao Wu. Regularized estimation of
linear functionals of precision matrices for high-dimensional time series.
IEEE Transactions on Signal Processing, 64(24):6459–6470, 2016.
[13] Alex Coad. Exploring the processes of firm growth: evidence from a vector
auto-regression. Industrial and Corporate Change, 19(6):1677–1703, 2010.
[14] John H. Cochrane. Permanent and transitory components of gnp and stock
prices. The Quarterly Journal of Economics, 109(1):241–265, 1994.
[15] Rainer Dahlhaus. Fitting Time Series Models to Nonstationary Processes.
The Annals of Statistics, 25(1):1–37, 1997.
[16] Rainer Dahlhaus. Locally stationary processes. arXiv:1109.4174, 2012.
[17] Dana Draghicescu, Serge Guillas, and Wei Biao Wu. Quantile Curve Estimation and Visualization for Nonostationary Time Series. Journal of
Computational and Graphical Statistics, 18(1):1–20, 2009.
[18] Jianqing Fan, Fang Han, and Han Liu. Challenges of big data analysis.
imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

29

National Science Review, (1):293–314,, 2014.
[19] Shaojun Guo, Yazhen Wang, and Qiwei Yao. High dimensional and banded
vector autoregressions. arXiv:1502.07831, 2015.
[20] J.D. Hamilton. Time Series Analysis. Cambridge University Press, 2 edition, 1994.
[21] Fang Han, Huanran Lu, and Han Liu. A direct estimation of high dimensional stationary vector autoregressions. Journal of Machine Learning
Research, 16:3115–3150, 2015.
[22] Fang Han, Sheng Xu, and Han Liu. Rate-optimal estimation of high dimensional time series. Preprint, 2016.
[23] N. J. Hsu, H.L. Hung, and Y.M. Chang. Subset selection for vector autoregressive process using lasso. Computational Statistics and Data Analysis,
52(7):3645–3657, 2008.
[24] Xingguo Li, Tuo Zhao, Xiaoming Yuan, and Han Liu. An R Package flare
for High Dimensional Linear Regression and Precision Matrix Estimation.
Journal of Machine Learning Research, to appear, 2014.
[25] Helmut Lütkepohl. New Introduction to Multiple Time Series Analysis.
Springer, 2007.
[26] Katta G Murty. Linear programming. John Wiley & Sons, New York, 1983.
[27] Giorgio E. Primiceri. Time varying structural vector autoregressions and
monetary policy. Review of Economic Studies, 72(3):821–852, 2005.
[28] Pradeep Ravikumar, Martin J. Wainwright, Garvesh Raskutti, and Bin
Yu. High-dimensional covariance estimation by minimizing `1 -penalized
log-determinant divergence. Electronic Journal of Statistics, 5:935–980,
2011.
[29] Mark Rudelson and Roman Vershynin. Hanson-Wright inequality and subGaussian concentration. Electronic Communication in Probability, 18(7):1–
9, 2013.
[30] Ali Shojaie and George Michailidis. Discovering graphical granger causality
using the truncating lasso penalty. Bioinformatics, 26(18):517–523, 2010.
[31] Christopher A. Sims. Macroeconomics and reality. Econometrica, 48(1):1–
48, 1980.
[32] Christopher A. Sims. Interpreting the macroeconomic time series facts:
The effects of monetary policy. European Economic Review, 36(5):975–
1000, 1992.
[33] Song Song and Peter J. Bickel.
Large vector auto regressions.
arXiv:1106.3915, 2011.
[34] J.H. Stock and M.W. Watson. Implications of dynamic factor models for
var analysis. Working Paper No. 11467, National Bureau of Economic
Research, 2005.
[35] H. Wang, G Li, and C.L. Tsai. Regression coefficient and autoregressive
order shrinkage and selection via the lasso. Journal of the Royal Statistical
Society: Series B, 69(1):63–78, 2007.
[36] Christopher K Wikle and Mevin B Hooten. A general science-based framework for dynamical spatio-temporal models. Test, 19(3):417–451, 2010.
[37] Shuheng Zhou, John Lafferty, and Larry Wasserman. Time Varying Undiimsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

30

rected Graphs. Machine Learning, 80:295–319, 2010.

Appendix A: Detailed setup of the simulation studies
We use sugm.generator() to genereate sparce matrice from multivariate normal distributions with four graph structures- hub, cluster, band and random.
Besides taking the patterns of graphs, the sugm.generator() also takes other
four critical input arguments which are g, prob, v and u shown in Table 7. The
detailed setup of these fours parameters in our simulation is shown in Table 8.
Table 7
The meanings of sugm.generator()’s arguments
Arguments

Meanings
Represents the number of hubs or clusters when the pattern is “hub” or
“cluster”
Represents the probability that an off-diagonal entry will be non-zero.
Assigns the values to the off-diagonal entries in the precision matrix and
controls the magnitude of partial correlations with u.
A positive number which is added to the diagonal entries of the precision
matrix.

g
prob
v
u

Table 8
Detailed setup of sugm.generator()
Hub
d
20
30
40
50

g
8
10
15
20

prob
NULL
NULL
NULL
NULL

Cluster
v
0.001
0.001
0.001
0.001

u
10
10
10
10

d
20
30
40
50

g
8
10
15
20

Band
d
20
30
40
50

g
1
1
1
1

prob
NULL
NULL
NULL
NULL

v
0.001
0.001
0.001
0.001

prob
NULL
NULL
NULL
NULL

v
0.001
0.001
0.001
0.001

u
10
10
10
10

Random
u
10
10
10
10

d
20
30
40
50

g
NULL
NULL
NULL
NULL

prob
0.001
0.001
0.001
0.001

v
0.001
0.001
0.001
0.001

u
10
10
10
10

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

31

Appendix B: Examples for four graph structures

Fig 5: Examples of the baseline coefficient matrices in four graph structures.
Zero entries are colored in grey. Positive entries and negative entries are colored
in red and black respectively.
(a) Hub

(b) Cluster

(c) Band

(d) Random

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

32

Appendix C: ROC curves based on the true thresholding
parameters u]

Fig 6: ROC curves under different settings
(a) Pattern is hub

(b) Pattern is cluster

0.6
0.4

True Positive Rate (Sensitivity)

d=20
0.2

0.6
0.4

d=20
0.2

True Positive Rate (Sensitivity)

0.8

cluster

0.8

hub

d=30

d=30

d=50
0.0

d=40

d=50
0.0

d=40

0.0

0.2

0.4

0.6

0.8

0.0

0.2

False Positive Rate (1−Specificity)

(c) Pattern is band

0.6

0.8

(d) Pattern is random

0.2

d=30

0.6
0.4

True Positive Rate (Sensitivity)

d=20

d=20
0.2

0.4

0.6

0.8

random

0.8

band

True Positive Rate (Sensitivity)

0.4
False Positive Rate (1−Specificity)

d=30

d=50
0.0

d=40

d=50
0.0

d=40

0.0

0.2

0.4
False Positive Rate (1−Specificity)

0.6

0.8

0.0

0.2

0.4

0.6

0.8

False Positive Rate (1−Specificity)

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

X. Ding et al./High-dimensional VAR

33

Appendix D: Selected stocks and currencies in real data analysis

Table 9
Part of the 30 selected stocks

1
2
3
4
5
6
7
8
9
10

Stock (Abb.)

Sector (Abb.)

Kellogg Co. (Kellogg)
Target Corp. (Target)
Boeing Company (Boeing)
CME Group Inc. (CME)
Prudential Financial (Prudential)
Edison Int’l (Edison)
Lockheed Martin Corp. (Lockheed)
PepsiCo Inc. (PepsiCo)
Hartford Financial Svc. GP (Hartford)
Exxon Mobil Corp. (Exxon)

Consumer Staples (CS)
Consumer Discretionary (CD)
Industrial (IND)
Financials (FIN)
Financials (FIN)
Utilities (UTI)
Industrial (IND)
Consumer Staples (CS)
Financials (FIN)
Energy (ENE)

Table 10
Selected currencies
Currencies
1
2
3
4
5
6
7
8

Euro / U.S.
U.K. / U.S.
Swiss / U.S.
Australia / U.S.
New Zealand / U.S.
Canada / U.S.
Singapore / U.S.
Brazil / U.S.

Currencies
9
10
11
12
13
14
15

Mexico / U.S.
South Korea / U.S.
India / U.S.
Thailand / U.S.
South Africa / U.S.
Norway / U.S.
Sweden / U.S.

imsart-ejs ver. 2014/10/16 file: tvvar_Sep-29-2017_EJS1325.tex date: October 3, 2017

