Inference of high-dimensional linear models with
time-varying coefficients
Xiaohui Chen, Yifeng He
University of Illinois at Urbana-Champaign
Email: xhchen@illinois.edu, he3@illinois.edu

arXiv:1506.03909v3 [stat.ME] 16 Mar 2017

March 17, 2017
Abstract
We propose a pointwise inference algorithm for high-dimensional linear models with time-varying coefficients. The method is based on a novel combination
of the nonparametric kernel smoothing technique and a Lasso bias-corrected ridge
regression estimator. Due to the non-stationarity feature of the model, dynamic
bias-variance decomposition of the estimator is obtained. With a bias-correction
procedure, the local null distribution of the estimator of the time-varying coefficient
vector is characterized for iid Gaussian and heavy-tailed errors. The limiting null
distribution is also established for Gaussian process errors, and we show that the
asymptotic properties differ between short-range and long-range dependent errors.
Here, p-values are adjusted by a Bonferroni-type correction procedure to control the
familywise error rate (FWER) in the asymptotic sense at each time point. The finite
sample size performance of the proposed inference algorithm is illustrated with synthetic data and an application to learn brain connectivity by using the resting-state
fMRI data for Parkinson’s disease.

1

Introduction

We consider the time-varying coefficient models (TVCM)
y(t) = x(t)> β(t) + e(t)

(1)

where t ∈ [0, 1] is the time index, y(·) the response process, x(·) the p × 1 deterministic
predictor process, β(·) the p × 1 time varying coefficient vector, and e(·) the mean zero
stationary error process. The response and predictors are observed at ti = i/n, i = 1, ..., n,
i.e. yi = y(ti ), xi = x(ti ), and ei = e(ti ) with a known covariance matrix Σe = Cov(e) where
1

e = (e1 , · · · , en )> . TVCM is useful for capturing the dynamic associations in the regression
models and longitudinal data analysis [16], and it has broad applications in biomedical
engineering, environmental science, and econometrics. In this paper, we focus on the
pointwise inference for the time-varying coefficient vector β(t) in the high-dimensional
double asymptotic framework min(p, n) → ∞.
Nonparametric estimation and inference of the TVCM in fixed dimension has been
extensively studied, see e.g. [27, 16, 12, 15, 39, 25, 8, 38, 41]. In high dimensions, variable
selection and estimation of varying-coefficient models using basis expansions have been
studied in [35, 36, 31]. Our primary objective is not to estimate β(t), but rather to
perform statistical inference on the coefficients. In particular, for any t ∈ (0, 1), we wish
to test the local hypothesis, for j = 1, · · · , p,
H0,j,t : βj (t) = 0 VS H1,j,t : βj (t) 6= 0.

(2)

By assigning p-values at each time point, we construct a sequence of estimators of the
coefficient vectors that allows us to assess the uncertainty of the dynamic patterns in such
as brain connectivity networks. Confidence intervals and hypothesis testing problems of
lower-dimensional functionals of the high-dimensional constant coefficient vector β(t) ≡
β, ∀t ∈ [0, 1], have been studied in [5, 37, 19]. To the best of our knowledge, little has been
done for inference of the high-dimensional TVCM and our goal is to fill the inference gap
between the classical TVCM and the high-dimensional linear model.
While the existing literature on high-dimensional linear models is based on iid errors,
[5, 37, 19], we provide an asymptotic theory for answering the question that to which
extent the statistical validity of inferences based on iid errors can hold for dependent
errors. Allowing temporal dependence is of the practical interest as many datasets such as
fMRI data are spatio-temporal and the errors are naturally correlated in the time domain.
Theoretical analysis has revealed that the temporal dependence has delicate impact on the
asymptotic rates for estimating the covariance structures, [10, 11]. Therefore, it is useful
to build an inference procedure that is also robust in the time series context. The error
process ei is modelled as a stationary linear process
ei =

∞
X

am ξi−m ,

m=0

2

(3)

where a0 = 1 and ξi are iid mean-zero random variables (a.k.a. innovations) with variance
σ 2 . When the ξi are normal, the linear processes of form (3) are Gaussian processes that
cover the autoregressive and moving-average (ARMA) models with iid Gaussian innovations as special cases. For the linear process, we deal with both weak and strong temporal
dependences. In particular, if am = O((m + 1)−% ) and % > 1/2, then ei is well-defined
and has (i) short-range dependence (SRD) if % > 1, (ii) long-range dependence (LRD) if
P
1/2 < % < 1. For the SRD processes, it is clear that ∞
m=0 |am | < ∞ and therefore the
long-run variance is finite.
The paper is organized as follows. In Section 2, we describe our method in details.
Asymptotic theory is presented in Section 3. Section 4 presents some simulation results
and Section 5 demonstrates an application to an fMRI dataset. The paper concludes in
Section 6 with a discussion of some future work. Proofs and some implementation issues
are available in the Appendix.

2
2.1

Method
Notations and Preliminary

Let K be a non-negative symmetric function with bounded support in [−1, 1],
−1

1, and let bn be a bandwidth parameter satisfying bn = o(1) and n

R1
−1

K(x)dx =

= o(bn ). For each

time point t ∈ $ = [bn , 1 − bn ], the Nadaraya-Waston smoothing weight is defined as

 Pn Kbn (ti −t)
if |ti − t| ≤ bn
m=1 Kbn (tm −t)
w(i, t) =
,
(4)

0
otherwise
where Kb (·) = K(·/b). Let Nt = {i : |ti − t| ≤ bn } be the bn -neighborhood of time t,
|Nt | be the cardinality of the discrete set Nt , Wt = diag(w(i, t)i∈Nt ) be the |Nt | × |Nt |
diagonal matrix with w(i, t), i ∈ Nt on the diagonal, and let Rt = span(xi : i ∈ Nt ) be the
subspace in Rp spanned by xi , the rows of design matrix X in the Nt neighborhood. Let
1/2
1/2
Xt = (w(i, t)1/2 xi )>
yi )>
ei )>
i∈Nt , Yt = (w(i, t)
i∈Nt , and Et = (w(i, t)
i∈Nt . Denote Ip as the

p × p identity matrix. We write the singular value decomposition (SVD) of Xt as
Xt = P DQ>
3

(5)

where P and Q are |Nt | × r, and p × r matrices such that P > P = Q> Q = Ir , and
D = diag(d1 , · · · , dr ) is a diagonal matrix containing the r nonzero singular values of Xt .
Now let PRt be the projection matrix onto Rt ,
PRt = Xt> (Xt Xt> )− Xt = QQ> ,

(6)

where (Xt Xt> )− = P D−2 P > is the pseudo-inverse matrix of Xt Xt> . Let θ(t) = PRt β(t) be
the projection of β(t) onto Rt such that B(t) = θ(t) − β(t) is the projection bias. Let
1/2

1/2

Ω(λ) = (Xt> Xt + λIp )−1 Xt> Wt Σe,t Wt Xt (Xt> Xt + λIp )−1

(7)

be the covariance matrix of the time-varying ridge estimator defined in (9), where Σe,t =
Cov((ei )i∈Nt ) and λ > 0 is the shrinkage parameter of the ridge estimator. Let Ωmin (λ) =
minj≤p Ωjj (λ) be the smallest diagonal entry of Ω(λ). For a generic vector b ∈ Rp , we write
P
P
|b|q = ( pj=1 |bj |q )1/q if q > 0, and |b|0 = pj=1 1(bj 6= 0) if q = 0. Let wt = inf i∈Nt w(i, t)
and wt = supi∈Nt w(i, t). For an n × n square symmetric matrix M and an n × m rectangle
matrix R, we use ρi (M ) and σi (R) to denote the i-th largest eigenvalues of M and singular
values of R, respectively. If k = rank(R), then σ1 (R) ≥ σ2 (R) ≥ · · · ≥ σk (R) > 0 =
σk+1 (R) = · · · = σmax(m,n) (R), zeros being padded to the last max(m, n) − k singular
values. We take ρmax (M ), ρmin (M ) and ρmin6=0 (M ) as the maximum, minimum and nonzero
minimum eigenvalues of M , respectively, and |M |∞ = max1≤j,k≤p |Mjk |. Let
ρmax (M, s) =

b> M b
.
|b|0 ≤s,b6=0 b> b
max

If M is nonnegative definite, then ρmax (M, s) is the restricted maximum eigenvalues of M
at most s columns and rows.
The p-dimensional coefficient vector β(t) is decomposed into two parts via projecting
onto the |Nt |-dimensional linear subspace spanned by the rows of Xt and its orthogonal
complement; see Figure 1(a). A key advantage of this decomposition is that the projected
part can be conveniently estimated in closed-form, for example, by the ridge estimator since
it lies in the row space of Xt and thus is amenable for the subsequent inferential analysis.
In the high-dimensional situation, this projection introduces a non-negligible shrinkage
bias in estimating β(t) and therefore we may lose information because p  |Nt |. On the
other hand, the shrinkage bias can be corrected by a consistent estimator of β(t). As a
4

(a) Bias correction by projection to the (b) Smoothly time-varying row space
row space of Xt .

of Xt .

Figure 1: Intuition of the proposed algorithm in Section 2.2.
particular example, we use the Lasso estimator, though any sparsity-promoting estimator
attaining the same convergence rate as the Lasso should work. Because of the time-varying
nature of the nonzero functional β(t), the smoothness on the row space of Xt along the
time index t is necessary to apply nonparametric smoothing technique; see Fig. 1(b). As
a special case, when the nonzero components β(t) ≡ β are constant functions and the
error process is iid Gaussian, our algorithm is the same as that of [5]. Here, we emphasize
that (i) coefficient vectors are time-varying (i.e. non-constant), (ii) errors are allowed to
have heavy-tails by assuming milder polynomial moment conditions and to have temporal
dependence, including both SRD and LRD processes. There are other inferential methods
for high-dimensional linear models such as [37, 19]. We do not explore specific choices
here since our contribution is a general framework of combining nonparametric smoothing
and bias-correction methods to make inference for high-dimensional TVCM. However, we
expect that a non-stationary generalization would be feasible for those methods as well.
Some simulation comparisons are provided for time-varying versions of the bias-correction
methods in Section 4.

5

2.2

Inference algorithm

First, we estimate the projection bias B(t) by B̃(t) = (PRt − Ip )β̃(t), where β̃(t) is the
time-varying Lasso (tv-Lasso) estimator
X
2
β̃(t) = arg minp
w(i, t)(yi − x>
i b) + λ1 |b|1
b∈R

(8)

i∈Nt

= arg minp |Yt − Xt b|22 + λ1 |b|1 .
b∈R

Next, we estimate θ(t) = PRt β(t) using the time-varying ridge (tv-ridge) estimator
X
2
2
θ̃(t) = arg minp
w(i, t)(yi − x>
i b) + λ2 |b|2
b∈R

i∈Nt

= (Xt> Xt + λ2 Ip )−1 Xt> Yt .

(9)

We defer the discussion of tuning parameters choice λ1 and λ2 to Section 3. Our tv-Lasso
bias-corrected tv-ridge regression estimator for β(t) is
β̂(t) = θ̃(t) − B̃(t).

(10)

Based on β̂(t) = (β̂1 (t), · · · , β̂p (t))> , we calculate the raw two-sided p-values for individual
coefficients
"
P̃j = 2 1 − Φ

|β̂j (t)| −

λ11−ξ maxk6=j
1/2
Ωjj (λ2 )

|(PRt )jk |

!#
,

j = 1, · · · , p,

(11)

where ξ ∈ [0, 1) is user pre-specified number that depends on the number of nonzero β(t).
In particular, if |supp(β(t))| is bounded, then we can choose ξ = 0. Generally, following
[5], we use ξ = 0.05 in our numeric examples to allow the number of nonzero components
in β(t) to diverge at proper rates. Let v(t) = (V1 (t), · · · , Vp (t))> ∼ N (0, Ω(λ2 )) and define
the distribution function


h

i
−1/2
F (z) = P min 2 1 − Φ Ωjj (λ2 )|Vj (t)| ≤ z .
j≤p

(12)

We adjust the P̃j for multiplicity by Pj = F (P̃j + ζ), where ζ is another pre-defined small
number [5] that accommodates asymptotic approximation errors. Our decision rule is
defined as: reject H0,j,t if Pj ≤ α for α ∈ (0, 1). For iid errors, since Σe = σ 2 In and
Ω(λ2 ) = σ 2 (Xt> Xt + λ2 Ip )−1 Xt> Wt Xt (Xt> Xt + λ2 Ip )−1 ,
6

we see that F (·) is independent of σ. Therefore, F (·) can be easily estimated by repeatedly
sampling from the multivariate Gaussian distribution N (0, Ω(λ2 )). A similar observation
has been made in [5].

3

Asymptotic results

In this section, we present the asymptotic theory of the inference algorithm in Section 2.2.
First, we state the main assumptions for iid Gaussian errors.
1. Error. The errors ei ∼ N (0, σ 2 ) are independent and identically distributed (iid).
2. Sparsity. β(·) is uniformly s-sparse, i.e. supt∈[0,1] |St∗ | ≤ s, where St∗ = {j : βj (t) 6=
0} is the support set.
3. Smoothness.
(a) β(·) is twice differentiable with bounded and continuous first and second derivatives in the coordinatewise sense, i.e. βj (·) ∈ C 2 ([0, 1], C0 ) for each j = 1, · · · , p
and C0 is an upper bound for the partial derivatives.
(b) The bn -neighborhood covariance matrix Σ̂t = |Nt |−1

P

i∈Nt

> 
xi x>
i := Xt Xt

satisfies
ρmax (Σ̂t , s) ≤ ε−2
0 < ∞.

(13)

4. Non-degeneracy.
lim inf Ωmin (λ) > 0.

(14)

λ↓0

5. Identifiability.
(a) The minimum nonzero eigenvalue condition
ρmin6=0 (Σ̂t ) ≥ ε20 > 0.
(b) The restricted eigenvalue condition
(
φ0 = inf

φ > 0 : min

|S|=s

inf

|bS c |1 ≤3|bS |1

b> Σ̂t b
φ2
≥
|bS |21
s

(15)

)
holds for all t ∈ [0, 1]

> 0,
(16)

where Σ̂t = Xt> Xt is the kernel smoothed covariance matrix of the predictors.
7

6. Kernel. The kernel function K(·) is nonnegative, symmetric around 0 with bounded
support in [−1, 1].
Here, we comment the assumptions and their implications. Assumption 1 and 6 are standard. The Gaussian distribution is non-essential and can be relaxed to sub-Gaussian and
heavier tailed distributions (Theorem 3.4). Assumption 2 is a sparsity condition for the
nonzero functional components and allows that s → ∞ slower than min(p, n). It is a key
condition for maintaining the low-dimensional structure when the dimension p grows with
the sample size n. By the argument of Theorem 5 in [40], it implies that the number of the
first and second non-vanishing derivatives of β(t) is bounded by s almost surely on [0, 1].
Assumption 3 ensures the smoothness of the time-varying coefficient vectors and the design
matrix so that nonparametric smoothing techniques are applicable. Examples of Assumption 3(a) include the quadratic functions β(t) = β + αt + ξt2 /2 and the periodic functions
β(t) = β + α sin(t) + ξ cos(t) with |α|∞ + |ξ|∞ ≤ C0 . Assamption 3(b) can be viewed as
Lipschitz continuity on the local design matrix that is smoothly evolving, [41]. It is weaker
than the condition that ρmax (Σ̂t ) ≤ ε−2
0 because the latter may grow to infinity much faster
than the restricted form (13). Assumption 4 is required for a non-degenerated stochastic
component of the proposed estimator which is used for the inference purpose. Assumption
5(a) and 5(b) together impose the identifiability conditions for recovering the coefficient
vectors. The analogous condition of the time-invariant version has been extensively used
in literature to derive theoretical properties of the Lasso model; see e.g. [2, 33].
For the tv-lasso bias-corrected tv-ridge estimator (10), we have the following representation theorem.
Theorem 3.1 (Representation). Fix t ∈ $ and let
"
#1/2
X
` 2
Lt,` = max
w(i, t) Xij
, ` = 1, 2, · · · ,
j≤p

λ0 = 4σLt,2

p
log p,

(17)

i∈Nt

and λ1 ≥ 2(λ0 + 2C0 Lt,1 bn (s|Nt |wt )1/2 ε−1
0 ). If λ2 = o(1), Assumptions 1-6 hold, and
C ≤ |Nt |wt ≤ |Nt |wt ≤ C −1 for some C ∈ (0, 1), then β̂(t) admits the decomposition
β̂(t) = β(t) + z(t) + γ(t),

(18)

z(t) ∼ N (0, Ω(λ2 )),
λ2 |θ(t)|2 + 2C0 s1/2 bn 4λ1 s
+ 2 |PRt − Ip |∞ ,
|γj (t)| ≤
Cε20
φ0
8

(19)
j = 1, · · · , p,

(20)

with probability tending to one. If βj (t) = 0, then we have
−1/2

Ωjj
where

(λ2 )(β̂j (t) − γj (t)) ∼ N (0, 1),

λ2 |θ(t)|2 + 2C0 s1/2 bn 4λ1 s
|γj (t)| ≤
+ 2 max |(PRt )jk |
Cε20
φ0 k6=j

(21)

(22)

with probability tending to one.
Remark 1. The decomposition (18) can be viewed as a local version of the one proposed
in [5] (Proposition 2). However, due to the time-varying nature of the nonzero coefficient
vectors, both the stochastic component z(t) in (19) and the bias component γ(t) in (20)
differ from [5]. First, our bound (20) for bias has three terms arising from: ridge shrinkage, non-stationarity and Lasso correction, and each has localized features depending on
the bandwidth bn of the sliding window and the smoothness parameter C0 . Second, the
stochastic part (19) also has time-dependent features in the covariance matrix (i.e. Ω(λ2 )
implicitly depends on t though Xt ) and the scale of normal random vector is different
from [5]. Delicate balance among them allows us to perform valid statistical inference such
as hypothesis testing and confidence interval construction for the coefficients and, more
broadly, their lower-dimensional linear functionals.
Example 3.1. Consider the uniform kernel K(x) = 0.5I(|x| ≤ 1) as an important special
case, the kernel used for our numeric experiments in Section 4. In this case, wt = (2nbn )−1
and |Nt |wt = |Nt |wt = 1. It is easily verified that under the local null hypothesis H0,j,t ,
(22) can be simplified to


1/2
γj (t) = O λ2 |θ(t)|2 + s bn + λ1 s max |(PRt )jk | .
k6=j

From this, it is clear that the three terms correspond to bias of ridge-shrinkage, nonstationarity and Lasso-correction. The first and last components have dynamic features
and the non-stationary bias is controlled by the bandwidth and sparsity parameters. The
condition C ≤ |Nt |wt ≤ |Nt |wt ≤ C −1 in Theorem 3.1 rules out the case that the kernel
does not use the boundary rows in the localized window and therefore avoids any jump in
the time-dependent row subspaces.

9

Remark 2. In Theorem 3.1, the penalty level for the tv-Lasso estimator λ1 can be chosen as
√
O(σLt,2 log p + Lt,1 s1/2 bn ). The second term in the penalty is due to the non-stationarity
of β(t) and the factor s1/2 arises from the weak coordinatewise smoothness requirement
on its derivatives (Assumption 3(a)). In the Lasso case with β(t) ≡ β and w(i, t) ≡ n−1 ,
P
an ideal order of the penalty level λ1 is σn−1 maxj≤p ( ni=1 Xij2 )1/2 (log p)1/2 see e.g. [2].
P
In the standardized design case n−1 ni=1 Xij2 = 1 so that Lt,1 = 1 and Lt,2 = n−1/2 , the
Lasso penalty is O(σ(n−1 log p)1/2 ), while the tv-Lasso has an additional term s1/2 bn that
may cause a larger bias. In our case, we estimate the time-varying coefficient vectors by
smoothing the data points in the localized window. Thus, it is unnatural to standardize the
reweighted local design matrix to have unit `2 length and the additional bias O(s1/2 bn ) is
due to non-stationarity. If the Xij are iid Gaussian random variables without standardization and we interpret the linear model as conditional on X, then, under the uniform kernel,
we have L2t,2 = OP (log p/|Nt |) and, in the Lasso case, the penalty level is O(σ|Nt |−1/2 log p).
If s = O(log p) and bn = O((log p/n)1/3 ), then the choice in Theorem 3.1 has the same
order as the Lasso with constant coefficient vector.
Based on Theorem 3.1, we can prove that the inference algorithm in Section 2.2 asymptotically controls the familywise error rate (FWER). Let α ∈ (0, 1) and FPα (t) be the
number of false rejections of H0,j,t based on the adjusted p-values. In the asymptotic
statement, p := p(n) is a function of n such that p → ∞ as n → ∞.
Theorem 3.2 (Pointwise inference: multiple testing). If the conditions of Theorem 3.1
hold and
λ2 |θ(t)|2 + s1/2 bn = o(Ωmin (λ2 )1/2 ),

(23)

then we have for each fixed t ∈ $
lim sup P(FPα (t) > 0) ≤ α.

(24)

n→∞

The proof of Theorem 3.2 is standard by combining the argument of Theorem 2 in [5]
and Theorem 3.1. Therefore, we omit the proof. Condition (23) essentially requires that
the shrinkage and non-stationarity biases of the tv-ridge estimator together are dominated
by the variance; see also the representation (18), (19), (20), and (21). This is mild condition
for two reasons. First, in view that variance of the tv-ridge estimator is lower bounded
when λ2 is small enough; c.f. (14), the first term is quite weak in the sense that the
10

tv-ridge estimator acts on a much smaller subspace with dimension |Nt | than the original
p-dimensional vector space. Second, for the choice of penalty parameter of λ1 in Theorem
3.1, the term s1/2 bn in (23) is at most λ1 . Hence, the bias correction (including the
projection and non-stationary parts) in the inference algorithm (11) has a dominating
effect on the second term of (23). Consequently, provided λ2 is small enough, the bias
correction step in computing the raw p-value asymptotically approximates the stochastic
component in the tv-ridge estimator.
Remark 3. The Bonferroni correction (12) for the raw p-values is often conservative and
thus it may be sub-optimal in power. In our simulation studies, it seems that detection
power is reasonable while the FWER is controlled at 0.05; c.f. Table 1 and 2. To improve
the power, one can consider the control of the false discovery proportion (FDP) by the
principal factor approximation (PFA) method proposed in [14, 13]. By Theorem 3.1, under
the global null hypothesis H0,t : β1 (t) = · · · = βp (t) = 0, we have β̂(t)−γ(t) ∼ N (0, Ω(λ2 ))
with a known covariance matrix Ω(λ2 ). Therefore, our test statistic is jointly normal and
the PFA can be applied to control the FDP if Ω(λ2 ) can be well approximated by the
covariance matrix of a factor model plus a weakly dependent component. [14] provided a
practical procedure to estimate the FDP.
Next, we relax the iid assumption on the errors to allow temporal dependence.
Theorem 3.3 (Gaussian process errors). Suppose that the error process ei is a meanzero stationary Gaussian process of form (3) such that |am | ≤ K(m + 1)−% for some
% ∈ (1/2, 1) ∪ (1, ∞) and finite constant K > 0. Under Assumptions 2-6 and the notation
of Theorem 3.1 with
(
λ0 =

√
4σLt,2 |a|1 log p
√
C%,K σLt,2 n1−% log p

if % > 1
if 1 > % > 1/2

,

(25)

where a = (a0 , a1 , · · · )> , we have the representation of β̂(t) in (18)–(22) with probability
tending to one.
From Theorem 3.3, the temporal dependence strength has a dichotomous effect on the
choice of λ0 , and therefore on the asymptotic properties of β̂(t). For ei with SRD, we have
√
|a|1 < ∞ and λ0  σLt,2 log p. Therefore, the bias-correction part γ(t) of estimating
β(t) has the same rate of convergence as the iid error case. The temporal effect only plays
11

a role in the long-run covariance matrix of the stochastic part z(t). If ei has LRD, then
the temporal dependence has impact on both γ(t) and z(t). In addition, the choice of the
bandwidth parameter bn is different from the SRD and iid cases. In particular, the optimal
bandwidth for % ∈ (1/2, 1) is O((log p/n% )1/3 ) which is much larger than O((log p/n)1/3 )
in the iid and SRD cases, assuming s is bounded. The boundary case % = 1 can also be
characterized; details are omitted.
We also relax the moment condition on the errors that, in the iid error case, are assumed
to be zero-mean Gaussian. First, it is easy to relax this assumption to distributions with
sub-Gaussian tails (see Definition 7.1 in the Supplementary Material) and Theorem 3.1
and 3.2 continue to hold, in view that the large deviation inequality and the Gaussian
approximation for a weighted partial sum of the error process only depend on the tail
behavior and therefore on moments of ei . Second and more importantly, the sub-Gaussian
assumption may be knocked down to allow iid noise processes with algebraic tails, or
equivalently ei with moments up to a finite order. The consequence of this relaxation
is that a larger penalty parameter for the tv-Lasso is needed for errors with polynomial
moments. Let Ξ be the square root matrix of Ω(λ2 )/σ 2 (i.e. Ω(λ2 ) = σ 2 ΞΞ> ) and ξ j be
the j-th row of Ξ.
Theorem 3.4 (Heavy-tailed errors). Under the conditions of Theorem 3.1 with E|ei |q <
∞, q > 2, choose

λ0 = Cq max (pµn,q )1/q , σLt,2 (log p)1/2 ,
for large enough Cq > 0,
(26)
P
where µn,q = i∈Nt |w(i, t)Xij |q . If |ξ j |q = o(|ξ j |2 ) for all j = 1, · · · , p, then (18) holds
with probability tending to one and Theorem 3.2 holds.
The assumption |ξ j |q = o(|ξ j |2 ) is needed to ensure the asymptotic validity of the
Gaussian approximation of the ridge component (19) for non-Gaussian data.

4
4.1

Simulation studies
Simulation setup

In the simulation studies, we generated the n × p design matrix with iid rows sampled
from N (0, ΣX ) for n = 300 and p = 300. We considered two covariance structures on the
12

design matrix: (i) ΣX = Ip ; (ii) ΣX = T , where T = (tjk )pj,k=1 and tjk = 0.5|j−k| . The
time-varying coefficient vectors β(t) had s = 3 non-zero elements and p − 3 zeros for all
t ∈ [0, 1]. The non-zero elements in β(t) were generated by sampling nodes from a uniform
distribution U (−2.5, 2.5) at regular time points and smoothly interpolating on the interval
[0, 1] using the cubic splines. We simulated the following stationary error processes.
1. The ei are iid N (0, 1).
2. ei = ϕei−1 + ξi is an AR(1) process where ϕ ∈ {0.2, 0.5} and ξi are iid N (0, 1).
√
3. The ei are iid Student’s t(3)/ 3.
4. ei is a long-memory process ei =

P∞

m=0 (m

+ 1)−% ξi−m , where % = 0.75 and ξi are iid

Gaussian with mean zero.
We compared the performance of the proposed method with the following.
1. (TV-Lasso) - The time-varying Lasso, the kernel smoothed time-varying LASSO
defined in (8), where λ1 is selected by the cross-validation (CV).
2. (FP-Lasso) - The false-positive Lasso, where λ1 is tuned to match the FWER of the
proposed method. This allowed us to compare the power at similar levels of FWER.
3. (TV-LDPE) - An adaptation of the de-biased LASSO inference procedure by [37] to
the kernel smoothed, time-varying setting.
4. (TV-SDL) - An adaptation of the SDL test of [19] to the kernel smoothed, timevarying setting.
5. (Non-TV) - The original non-time-varying method of [5] that ignores the dynamic
structures. The penalty parameter λ1 in the Lasso is set to the scaled Lasso parameter
p
2 log(p)/n.
In all time-varying models, we used the kernel bandwidth bn = 0.1. For the proposed
method, we used λ2 = 1/n and ζ = 0. We let Pj,t,m be the multiplicity-adjusted p-value
for testing the hypothesis H0,j,t : βj (t) = 0 for t ∈ $ = [bn , 1 − bn ] in the m-th Monte
Carlo simulation for m = 1, · · · , M . For TV-LDPE, TV-SDL, the proposed method and
its non-tv version, we adopted the following performance measures.
13

1. The (averaged) false positive rate (FPR) over the interval $,
M
XXX
1
1(Pj,t,m ≤ α).
n(1 − 2bn )(p − s)M t∈$ j∈S c m=1

2. The (averaged) false negative rate (FNR) over the interval $,
M
XXX
1
1(Pj,t,m > α).
n(1 − 2bn )sM t∈$ j∈S m=1

3. The (averaged) FWER over the interval $,
M
XX
1
FWER =
1(min Pj,t,m ≤ α).
n(1 − 2bn )M t∈$ m=1 j∈S c

For the Lasso-based methods (TV-Lasso and FP-Lasso), the probabilities are replaced
by the corresponding indicators of whether or not the estimated coefficients are zero.

4.2

Empirical results

For each simulation setup, we report the FPR, RNR, FWER, and the root mean square
errors (RMSE) of the estimates. The results are shown in Tables 1 and 2, from which
we make several observations. First, TV-Lasso and the method of [5] do not control the
FWER, while the proposed method can control the FWER at the nominal level α = 0.05
in all setups. Second, the proposed method has uniformly higher power than FP-Lasso,
the TV-Lasso tuned to match the FWER with our method. This is probably explained
by the bias of the `1 regularization in the TV-Lasso. Third, for the design matrix with iid
Gaussian entries, TV-LDPE and TV-SDL have comparable performance as our proposed
method in terms of the power, while all three methods have the FWER controlled below
0.05. TV-LDPE and TV-SDL are more sensitive to the design matrix than the proposed
method; for the Toeplitz design matrix T , TV-LDPE and TV-SDL seem to lose control on
the FWER. Moreover, the FWER, FNR, and RMSE are larger as the dependence level of
the error process grows and as the tail of the error distribution becomes thicker. Finally,
the proposed method is computationally more economical than the competing methods
TV-LDPE and TV-SDL. Table 3 shows the runtimes on an Intel i5-4790K with the Intel
MKL linear algebra libraries (software and platform: R 3.2.2 for Windows).
14

Table 1: Simulation results. n = 300, p = 300, s = 3.
iid

iid

xi ∼ N (0, Ip ), e ∼ N (0, In )
Method

FPR

FNR

FWER

TV-Lasso

7.51 ×

10−2

0.0551

1

FP-Lasso

1.50 × 10−4

0.2352

0.0344

Proposed

1.50 × 10−4

0.1889

TV-LDPE

1.29 × 10−4

0.1981

TV-SDL

1.53 × 10−4

Non-TV

4.89 × 10−1

xi ∼ N (0, T ), e ∼ N (0, In )
RMSE

FWER

RMSE

0.0537

1.55 ×

0.0684

1

0.0520

0.1124

2.81 × 10−5

0.5170

0.0063

0.1318

0.0346

0.1838

2.81 × 10−5

0.4072

0.0063

0.1984

0.0254

0.2652

3.77 × 10−4

0.3615

0.0743

0.2727

0.1848

0.0357

0.1316

1.08 × 10−3

0.3006

0.2119

0.1409

0.5100

0.4600

0.5762

2.58 × 10−1

0.7100

0.2400

0.6084

iid

iid

xi ∼ N (0, T ), e ∼ AR(1) with ϕ = 0.2

FPR

FNR

FWER

RMSE

FPR

FNR

TV-Lasso

8.93 × 10−2

0.0563

1

0.0533

1.51 × 10−1

0.0629

1

0.0519

FP-Lasso

1.78 ×

10−4

0.2363

0.0384

0.1099

7.21 × 10−5

0.4709

0.0173

0.1302

Proposed

1.78 × 10−4

0.1891

0.0376

0.1836

7.21 × 10−5

0.3434

0.0173

0.1920

TV-LDPE

9.85 ×

10−5

0.1995

0.0215

0.2799

3.70 × 10−4

0.3620

0.0843

0.2725

TV-SDL

1.97 × 10−4

0.1883

0.0419

0.1374

1.07 × 10−3

0.3032

0.2165

0.1404

Non-TV

4.29 × 10−1

0.5633

0.4100

0.5557

2.58 × 10−1

0.6933

0.2200

0.6088

iid

xi ∼ N (0, Ip ), e ∼ AR(1) with ϕ = 0.5
Method

5

FNR

10−1

xi ∼ N (0, Ip ), e ∼ AR(1) with ϕ = 0.2
Method

FPR

FWER

RMSE

iid

xi ∼ N (0, T ), e ∼ AR(1) with ϕ = 0.5

FPR

FNR

FWER

RMSE

FPR

FNR

FWER

RMSE

TV-Lasso

7.55 × 10−2

0.0544

1

0.0537

1.51 × 10−1

0.0611

1

0.0518

FP-Lasso

1.93 ×

10−4

0.2402

0.0431

0.1124

9.81 × 10−5

0.4805

0.0222

0.1303

Proposed

1.93 × 10−4

0.1809

0.0422

0.1836

9.81 × 10−5

0.3347

0.0235

0.1918

TV-LDPE

9.82 × 10−5

0.2028

0.0229

0.2756

3.71 × 10−4

0.3618

0.0849

0.2659

TV-SDL

1.94 × 10−4

0.1898

0.0425

0.1374

1.01 × 10−3

0.2993

0.2555

0.1439

Non-TV

5.49 × 10−1

0.4500

0.5200

0.5314

1.70 × 10−1

0.8300

0.1500

0.6004

Data example: learning brain connectivity

We illustrate our proposed method in an application to model the functional brain connectivity in a Parkinson’s disease study. The problem is to construct brain connectivity
networks from the resting-state functional magnetic resonance imaging (fMRI) data, where
slowly time-varying graphs have implications in modeling brain connectivity networks.
Traditional correlation analysis of the resting state blood-oxygen-level-dependent (BOLD)
signals of the brain showed considerable temporal variation on small time scales, [9, 18].
In view of the high spatial resolution of fMRI data, brain networks of subjects at rest are
believed to be structurally homogeneous with subtle fluctuations in some, but a small number of, connectivity edges, [20, 17]. A popular approach to learn brain connectivity is the
neighborhood selection procedure, [23]. Therefore, high-dimensional TVCM with a small
15

Table 2: Simulation results (continued). n = 300, p = 300, s = 3.
√
iid
iid
xi ∼ N (0, Ip ), ei ∼ t(3)/ 3

Method

FPR

FNR

TV-Lasso

9.14 × 10−2

0.0518

FP-Lasso

1.96 × 10−4

0.2193

Proposed

1.96 × 10−4

0.1708

TV-LDPE

1.26 ×

10−4

TV-SDL
Non-TV

FWER

√
iid
iid
xi ∼ N (0, T ), ei ∼ t(3)/ 3

RMSE

FPR

FNR

1

0.0539

1.57 × 10−1

0.0605

1

0.0506

0.0398

0.1120

4.20 × 10−5

0.4547

0.0124

0.1209

0.0439

0.1834

4.20 × 10−5

0.3129

0.0125

0.1885

0.2041

0.0275

0.2659

3.62 × 10−4

0.3043

0.0810

0.2719

1.90 × 10−4

0.1903

0.0403

0.1321

9.54 × 10−4

0.2814

0.1960

0.1381

5.16 × 10−1

0.4800

0.4800

0.5289

2.24 × 10−1

0.7700

0.1500

0.5962

iid

xi ∼ N (0, Ip ), e ∼ LRD with % = 0.75
Method

FWER

RMSE

iid

xi ∼ N (0, T ), e ∼ LRD with % = 0.75

FP(%)

FN(%)

FWER

RMSE

FPR

FNR

FWER

RMSE

TV-Lasso

7.25 × 10−2

0.0652

1

0.0499

1.77 × 10−1

0.0805

1

0.0556

FP-Lasso

1.60 ×

10−4

0.2496

0.0450

0.1091

1.37 × 10−4

0.5138

0.0229

0.1381

Proposed

1.60 × 10−4

0.1783

0.0433

0.1806

1.37 × 10−4

0.3376

0.0243

0.1953

TV-LDPE

1.08 ×

10−4

0.2067

0.0238

0.2653

3.68 × 10−4

0.3648

0.0859

0.2691

TV-SDL

2.10 × 10−4

0.1924

0.0501

0.1386

9.29 × 10−4

0.3014

0.0186

0.1448

Non-TV

5.19 × 10−1

0.4800

0.4800

0.5304

5.49 × 10−1

0.4500

0.4100

0.5280

Table 3: Runtime per 10 simulations.
Method

Runtime (in minutes)

TV-Lasso

0.5

FP-Lasso

9.5

Proposed

13

TV-LDPE

> 1000

TV-SDL

> 1000

Non-TV

< 0.5

number of nonzero components is a natural approach to study the time-evolving sparse
brain connectivity networks, in which the time-varying coefficients reflect the dynamic features of the corresponding edges in the networks. The neighborhood selection approach
is an approximation to the full multivariate distributions while ignoring the correlation
among the node-wise responses and this may cause certain power loss in finite samples.
[23] showed that in terms of variable selection, these two approaches are asymptotically
equivalent.
Our data example uses fMRI data collected from a study of patients with Parkinson’s
disease (PD) and their respective normal controls. PD is typically characterized by deviations in functional connectivity between various regions of the brain. Additionally, the
16

resting state functional connectivity has been shown as a candidate biomarker for PD
progression and treatment, where more advanced stages or manifestations of PD are associated with greater deviations from normal connectivity. Each resting state data matrix
in our example contains 240 time points and 52 brain regions of interest (ROI). The time
points are evenly sampled and the time indices are normalized to [0,1].
The brain connectivity network was constructed using the neighborhood selection procedure. In essence, it is a sequence of time-varying linear regressions enumerating each
ROI as the response variable and sparsely regressing on all the other ROIs. Figure 2 and 3
show the estimated graphs of a normal subject and a PD subject at three sequential time
points around t = 0.25 based on the proposed method. Red nodes are ROIs known to be
associated with motor control and blue nodes are ROIs either known to be unrelated to
motor control or whose functions in humans are not well understood. Different patterns of
connectivity in the networks can be found by comparing normal and PD subjects. From
the graphs, there are slow changes in the networks over time: most edges are preserved
on a small time scale, but there are a few edges evolving over time. For instance, in a PD
subject, ROI 1 and ROI 40 are unconnected in the first network but they are connected
in the second network and remain connected in the third network.
We also plot the connectivity graphs generated by the competing methods TV-LDPE
and TV-SDL; see Figure 4 and 5. These graphs are denser than with the proposed method
and are harder to interpret. Besides, as commented in the simulation studies, TV-LDPE
and TV-SDL are more computationally expensive. The method of [5] cannot be used to
capture the dynamic features of brain connectivity networks.

6

Discussions

This paper presents a pointwise inference algorithm for high-dimensional TVCM that can
asymptotically control the FWER. Based on the current work, an interesting improvement
would be to study simultaneous inference. Construction of the simultaneous confidence
band (SCB) for the time-varying coefficients is useful for testing their parametric forms in
high dimensions. This is a more challenging topic, which requires substantial additional
work and probability tools that are beyond the scope of this paper.
Our brain connectivity application is a subject-by-subject analysis. To perform the
17

Figure 2: Connectivity networks in control subject around t = 0.25 based on the proposed
method.

Figure 3: Connectivity networks in Parkinson’s Disease subject around t = 0.25 based on
the proposed method.

group analysis on the population level, a hierarchical linear model is more appropriate.
When p is fixed, the linear mixed-effects model is widely used in performing the multilevel group analysis in fMRI, [1, 21, 30]. The reason is that the generalized least squares
(GLS) estimator of a two-level model is inferentially equivalent to the GLS estimator of
the corresponding single-level model, provided that the second-level covariance is the sum
of the group covariance and the covariance of the first-level estimate, [1]. Extension of
the testing problem on the population parameters based on the ridge and Lasso estimates
when p → ∞ is another interesting future research topic.

18

Figure 4: A connectivity network based on TV-LDPE.

Figure 5: A connectivity network based on TV-SDL.

7

Proof

Proof of Theorem 3.1. Observe that Xt β(t) = Xt θ(t) since θ(t) = PRt β(t). Using the
closed-form formulae for the tv-ridge estimator (9) and by (36), we have
bias(θ̃(t)) = E(θ̃(t)) − θ(t) = (Xt> Xt + λ2 Ip )−1 Xt> [Xt θ(t) + Mt Xt β 0 (t) + Xt ξ] − θ(t), (27)
where |ξ|∞ ≤ C0 b2n /2 and |ξ|0 ≤ s almost surely, t ∈ $. First, we bound the shrinkage
bias of the tv-ridge estimator. By the argument in Section 3 of [28], we can show that
2
−1 >
(Xt> Xt + λ2 Ip )−1 Xt> Xt θ(t) − θ(t) = −Q(λ−1
2 D + Ir ) Q θ(t).

19

It follows from Lemma 7.2 that
|θ(t)|2
(28)
2
ρmin (λ−1
2 D + Ir )


λ2
λ2 |θ(t)|2
λ2 |θ(t)|2
,
=
|θ(t)|2 ≤
≤
2
λ2 + minj≤r dj
|Nt |wt ε20
ρmin6=0 (Σ̂t )

2
−1 >
|Q(λ−1
2 D + Ir ) Q θ(t)|2 ≤

where d2j = ρj (Σ̂t ), j = 1, · · · , r. Next, we deal with the non-stationary bias of the tvridge estimator (27) by a similar argument for (28). Indeed, let Q⊥ be the orthogonal
>
complement of Q such that Q>
⊥ Q⊥ = Ip−r and Q⊥ Q = 0(p−r)×r . Denote Γ = [Q; Q⊥ ].

Then, ΓΓ> = Γ> Γ = Ip . By the SVD of Xt , (5), we have
−1 > >
Γ Xt Mt Xt β 0 (t)
(Xt> Xt + λ2 Ip )−1 Xt> Mt Xt β 0 (t) = Γ Γ> (QD2 Q> + λ2 Ip )Γ
"
#
!−1 "
#
>
>
Q
Q
= [Q; Q⊥ ]
(QD2 Q> + λ2 Ip )[Q; Q⊥ ]
QDP > Mt Xt β 0 (t)
>
>
Q⊥
Q⊥
#
!"
0
2
−1
>
(D
+
λ
I
)
0
DP
M
X
β
(t)
2
r
t
t
= [Q; Q⊥ ]
0
λ−1
0
2 Ip−r
= Q(D + λ2 D−1 )−1 P > Mt Xt β 0 (t).
Hence, by Lemma 7.2 we have
|(Xt> Xt

−1

+ λ2 Ip )

Xt> Mt Xt β 0 (t)|2

bn |Xt β 0 (t)|2
C0 bn (s|Nt |wt )1/2 ε−1
0
≤
≤
,
ρmin (D + λ2 D−1 )
minj≤r (dj + λ2 /dj )

where wt = supi∈Nt w(i, t). Since λ2 = o(1) and dj ≥ (|Nt |wt )1/2 ε0 , the denominator of
last expression is lower bounded by [(|Nt |wt )1/2 ε0 + λ2 /((|Nt |wt )1/2 ε0 )] for large enough n.
Therefore, we have
|(Xt> Xt + λ2 Ip )−1 Xt> Mt Xt β 0 (t)|2 ≤

C0 bn s1/2
C0 bn (s|Nt |wt )1/2
≤
.
(|Nt |wt )1/2 ε20
Cε20

(29)

Similarly, an upper bound for the remainder term of (27) can be established. We have
|(Xt> Xt + λ2 Ip )−1 Xt> Mt Xt ξ|2 ≤

C0 b2n s1/2
,
2Cε20

for almost surely t ∈ $.

(30)

In addition, θ̃(t) − E[θ̃(t)] = (Xt> Xt + λ2 Ip )−1 Xt> Et is the stochastic part of the tv-ridge
estimator. Since the ei ∼ N (0, σ 2 In ) are iid, Et ∼ N (0, σ 2 Wt ). Hence, θ̃(t) − E[θ̃(t)] ∼
N (0, Ω(λ2 )), where Ω(λ) is defined in (7), and thus
Var(θ̃j (t)) = Ωjj (λ2 ) ≥ Ωmin (λ2 ).
20

(31)

Now, we consider the initial tv-lasso estimator. By Lemma 7.3,
|β̃(t) − β(t)|1 ≤ 4φ−2
0 λ1 s.

(32)

Then, (18), (19), and (20) follow by assembling (28), (29), (30), and (32) into (10),
β̂(t) = β(t) + bias(θ̃(t)) + {θ̃(t) − E[θ̃(t)]} − {(PRt − Ip )[β̃(t) − β(t)]}.
The marginal representation (21) and (22) follow from similar arguments by noting that
P
Bj (t) = k6=j (PRt )jk βk (t) under H0,j,t .
Proof of Theorem 3.3. The proof of Theorem 3.3 is similar to that of Theorem 3.1 so we
1/2

1/2

only highlight the difference involving the error process. First, Cov(Et ) = Wt Σe,t Wt .
Second, instead of using (37) in proving Lemma 7.3, we use Lemma 7.4 to get for all λ > 0
!


X
λ2
if % > 1,
P max
w(i, t)Xij ei ≥ λ ≤ 2p exp − 2
2 2
j≤p
σ
|a|
2L
1
t,2
i∈N
t

!
P max
j≤p

X

w(i, t)Xij ei ≥ λ


≤ 2p exp −

i∈Nt

C% λ2
L2t,2 n2(1−%) σ 2 K 2


if % ∈ (1/2, 1).

Proof of Theorem 3.4. The proof essentially follows the lines of that of Theorem 3.1, but
with differences in requiring a larger penalty parameter λ1 of the tv-Lasso. First, by the
Nagaev inequality [24], we have for any ε > 0,
!
X

pµn,q
2
P max
w(i, t)Xij ei ≥ σLt,2 ε ≤ (1 + 2/q)q κq
+
2p
exp
−c
ε
,
q
q
j≤p
(σL
t,2 ε)
i∈N
t

where cq = 2e−q (q + 2)−2 and κq is the q-th absolute moment of e1 . Then, choosing


(pµn,q )1/q
1/2
ε = Cq max
, (log p)
for large enough Cq > 0,
σLt,2
we have maxj≤p |

P

i∈Nt

1/2

w(i, t)Xij ei | = OP (λ0 ). Second, let Ξ = (Xt> Xt + λ2 Ip )−1 Xt> Wt


and Et = (ei )>
i∈Nt . Recall that θ̃(t) − E[θ̃(t)] = ΞEt . By the Gaussian approximation

21

2
) defined on a
[29, Theorem B], there exist iid Gaussian random variables gi ∼ N (0, σ 2 ξji

possibly richer probability space such that for every t > 0
!
P
q
X
i∈Nt E|ξji ei |
q
P θ̃j (t) − E[θ̃j (t)] −
gi ≥ t ≤ (Cq)
.
q
t
i∈N
t

Thus, it follows that θ̃j (t) − E[θ̃j (t)] = N (0, Ωjj (λ2 )) + OP (|ξ j |q ). As Ωjj (λ2 ) = σ 2 |ξ j |22 , the
proof is complete for by assumption, |ξ j |q = o(|ξ j |2 ).

Supplementary Material
The supplementary material contains additional technical lemmas and discusses some
implementation issues.

Acknowledgements
We would like to thank an anonymous referee, an Associate Editor, and the Co-Editor
for their many helpful comments that led to improvements of this paper. This research was
partially supported by NSF DMS-1404891 and UIUC Research Board Award RB15004. We
thank Aiping Liu (University of British Columbia) for providing the Parkinson’s Disease
data. The high-performance computing of this work was completed on the Illinois Campus
Cluster Program at UIUC.

References
[1] Christian F. Beckman, Mark Jenkinson, and Stephen M. Smith. General multilevel
linear modeling for group analysis in fmri. NeuroImage, 20:1052–1063, 2003.
[2] Peter Bickel, Ya’acov Ritov, and Alexandre Tsybakov. Simultaneous analysis of lasso
and dantzig selector. Annals of Statistics, 37(4):1705–1732, 2009.
[3] Peter J. Bickel and Elizaveta Levina. Regularized Estimation of Large Covariance
Matrices. The Annals of Statistics, 36(1):199–227, 2008.
[4] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration Inequalities:
A Nonasymptotic Theory of Independence. Oxford, 2013.
22

[5] Peter Bühlmann. Statistical significance in high-dimensional linear models. Bernoulli,
19(4):1212–1242, 2013.
[6] Peter Bühlmann and Sara van de Geer. Statistics for High-Dimensional Data: Methods, Theory and Applications. Springer Series in Statistics, 2011.
[7] Tony Cai, Cun-Hui Zhang, and Harrison Zhou. Optimal Rates of Convergence for
Covariance Matrix Estimation. The Annals of Statistics, 38(4):2118–2144, 2010.
[8] Zongwu Cai. Trending time-varying coefficient time series models with serially correlated errors. Journal of Econometrics, 136(1):163–188, 2007.
[9] Catie Chang and Gary H. Glover. Time-frequency dynamics of resting-state brain
connectivity measured with fMRI. NeuroImage, 50:81–98, 2010.
[10] 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.
[11] 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.
[12] W. S. Cleveland, E. Grosse, and W. M. Shyu. Local regression models. In Statistical
Models in S (Chambers, J. M. and Hastie, T. J., eds) Wadsworth & Brooks, Pacific
Grove., pages 309–376, 1991.
[13] Jianqing Fan and Xu Han. Estimation of the false discovery proportion with unknown
dependence. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2016.
[14] Jianqing Fan, Xu Han, and Weijie Gu. Estimating false discovery proportion under arbitrary covariance dependence. Journal of American Statistical Association,
107(499):1019–1035, 2012.
[15] Jianqing Fan and Wenyang Zhang. Statistical estimation in varying coefficient models.
Annals of Statistics, 27(5):1491–1518, 1999.
23

[16] Donald R. Hoover, John A. Rice, Colin O Wu, and Li-Ping Yang. Nonparametric smoothing estimates of time-varying coefficient models with longitudinal data.
Biometrika, 85(4):809–822, 1998.
[17] M. Hutchison, Mirsattari S., C. Jones, J. Gati, and L. Leugn. Functional networks in
the anesthetized rat brain revealed by independent component analysis of resting-state
fMRI. J. Neurophysiol, 103:3398–3406, 2010.
[18] M. Hutchison, T. Womelsdorf, J. Gati, S. Everling, and R. Menon. Resting-state
networks show dynamic functional connectivity in awake humans and anesthetized
macaques. Human Brain Mapping, 34:2154–2177, 2013.
[19] Adel Javanmard and Andrea Montanari. Confidence inter hypothesis testing in highdimensional regression under the gaussian random design model: Asymptotic theory.
IEEE Transactions on Information Theory, 60(10):6522–6554, 2014.
[20] V. Kiviniemi, H. Haanpää, J-H Kantola, Jauhiainen J., Vainionpää V., S. Alahuhta,
and O. Tervonen. Midazolam sedation increases fluctuation and synchrony of the
resting brain BOLD signal. Magn Reson Imaging, 23:531–537, 2005.
[21] Martin Lindquist. The statistical analysis of fmri data. Statistical Science, 23(4):439–
464, 2008.
[22] T. L. McMurry and D. N. Politis. Banded and tapered estimates for autocovariance
matrices and the linear process bootstrap. J. Time Ser. Anal., 31:471–482, 2010.
[23] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable
selection with the Lasso. Annals of Statistics, 34(3):1049–1579, 2006.
[24] S.V. Nagaev. Large deviations of sums of independent random variables. Annals of
Probability, 7(5):745–789, 1979.
[25] Susan Orbe, Eva Ferreira, and Juan Rodriguez-Poo. Nonparametric estimation of time
varying parameters under shape restrictions. Journal of Econometrics, 126(1):53–77,
2005.

24

[26] Sidney I. Resnick. Extreme Values, Regular Variation, and Point Processes. Applied
Probability. Springer-Verlag, 1987.
[27] Peter M. Robinson. Nonparametric estimation of time-varying parameters. In Statistical Analysis and Forecasting of Economic Structural Change (P. Hackl, ed.) Berlin:
Springer., 1989.
[28] Jun Shao and Xinwei Deng. Estimation in high-dimensional linear models with deterministic design matrices. Annals of Statistics, 40(2):812–831, 2012.
[29] Qi-Man Shao. Strong approximation theorems for independent random variables and
their applications. Journal of Multivariate Analysis, 52:107–130, 1995.
[30] Martha Skup. Longitudinal fmri analysis: a review of methods. Statistics and Its
Interface, 3(2):235–252, 2010.
[31] Rui Song, Feng Yi, and Hui Zou. On varying-coefficient independence screening for
high-dimensional varying-coefficient models. Statistica Sinica, 24(4):1735–1752, 2014.
[32] Tingni Sun and Chun-Hui Zhang. Scaled sparse linear regression. Biometrika, 99:879–
898, 2012.
[33] Sara van de Geer and Peter Bühlmann. On the conditions used to prove oracle results
for the Lasso. Electronic Journal of Statistics, 3:1360–1392, 2009.
[34] Boying Wang and Fuzhen Zhang. Some inequalities for the eigenvalues of the product
of positive semidefinite Hermitian matrices. Linear Algebra and Its Applications,
160:113–118, 1992.
[35] Fengrong Wei, Jian Huang, and Hongzhe Li. Variable selection and estimation in
high-dimensional varying-coefficient models. Statistica Sinica, 21:1515–1540, 2011.
[36] Lan Xue and Annie Qu. Variable selection in high-dimensional varying coefficient
models with global optimality. Journal of Machine Learning Research, 13:1973–1998,
2012.

25

[37] Cun-Hui Zhang and Stephanie S. Zhang. Confidence intervals for low dimensional
parameters in high dimensional linear models. Journal of the Royal Statistical Society:
Series B (Statistical Methodology), 71(1):217–242, 2013.
[38] Ting Zhang and Wei Biao Wu. Inference of time-varying regression models. Annals
of Statistics, 40(3):1376–1402, 2012.
[39] Wenyang Zhang, Sik-Yum Lee, and Xinyuan Song. Local polynomial fitting in semivarying coefficient model. Journal of Multivariate Analysis, 82:166–188, 2002.
[40] Shuheng Zhou, John Lafferty, and Larry Wasserman. Time-varying undirected graphs.
Machine Learning, 80(2-3):295–319, 2010.
[41] Zhou Zhou and Wei Biao Wu. Simultaneous inference of linear models with time
varying coefficients. Journal of the Royal Statistical Society, 72(4):513–531, 2010.

Department of Statistics, 725 S Wright Street, Champaign, IL USA 61820
E-mail: xhchen@illinois.edu, he3@illinois.edu

Supplementary Material
Additional technical lemmas
Lemma 7.1. Let X be an n × p matrix and D = diag(d1 , · · · , dn ) with |di | ≤ b and b ≥ 0.
Then
ρmax (X > DX, s) ≤ 2bρmax (X > X, s).
If di ∈ [0, b], then ρmax (X > DX, s) ≤ bρmax (X > X, s).
−
+
Proof. Let As = {a ∈ Rp : |a|2 ≤ 1, |a|0 ≤ s}. Write di = d+
i − di , where di = max(di , 0)

and d−
i = max(−di , 0) are the positive and negative parts, respectively. By definition
ρmax (X > DX, s) = max |a> X > DXa| = max |tr(D(Xaa> X > ))|
a∈As

= max
a∈As

a∈As

n
X

(d+
i

−

> >
d−
i )(Xaa X )ii

n
X
≤ 2b max
(Xaa> X > )ii
a∈As

i=1
>

>

a∈As

26

>

i=1

= 2b max tr(Xaa X ) = 2b max a X Xa = 2bρmax (X > X, s),
a∈As

>

because X > X is nonnegative definite. The second claim follows from the same lines with
d−
i = 0.
Lemma 7.2. Let t ∈ $ and Σ̂t be the kernel smoothed sample covariance at time t and
Σ̂t = Xt > Xt . Suppose that Xt has full row rank. Assume further (15), (13) and assumption 6 hold, then we have
ρmin6=0 (Σ̂t ) ≥ |Nt |wt ε20 ,

(33)

ρmax (Σ̂t , s) ≤ |Nt |wt ε−2
0 .

(34)

Proof. Since Xt is of full row rank, r = |Nt |. Note that Xt = (|Nt |Wt )1/2 Xt , ρi (Σ̂t ) =
σi2 (Xt ) and ρi (Σ̂t ) = σi2 (Xt ). By the generalized Marshall-Olkin inequality, see e.g. [34,
Theorem 4], assumption 6 and (15), we have
1/2

1/2

ρmin6=0 (Σ̂t ) = ρmin (Xt Xt> ) = |Nt |ρmin (Wt Xt Xt > Wt )
= |Nt |ρmin (Xt Xt > Wt ) ≥ |Nt |ρmin (Wt )ρmin (Xt Xt > ) ≥ |Nt |wt ε20 .
The second inequality (34) follows from assumption 3(b) and Lemma 7.1 applying to
Σ̂t = |Nt |Xt > Wt Xt and Wt ≥ 0.
Lemma 7.3. Suppose assumption 1, 2, 3 and 5(a) hold. Let t ∈ $ be fixed and λ0 be
defined in (17). Then, for λ1 ≥ 2(λ0 + 2C0 Lt,1 s1/2 ε−1
0 bn |Nt |w t ) where λ0 is defined in (17),
we have, with probability 1 − 2p−1 ,
|Xt [β̃(t) − β(t)]|22 + λ1 |β̃(t) − β(t)|1 ≤ 4λ21

s
.
φ20

(35)

Proof. By definition (8),
|Yt − Xt β̃(t)|22 + λ1 |β̃(t)|1 ≤ |Yt − Xt β(t)|22 + λ1 |β(t)|1 ,
which implies that
D
E
|Xt [β̃(t) − β(t)]|22 + λ1 |β̃(t)|1 ≤ λ1 |β(t)|1 + 2 Yt − Xt β(t), Xt [β̃(t) − β(t)] .
By assumption 2 and Taylor’s expansion in the bn -neighborhood of t, we see that
Yt − Xt β(t) = Et + Mt Xt β 0 (t) + Xt ξ,
27

(36)

where Mt = diag((ti − t)i∈Nt ) and ξ is a vector such that |ξ|∞ ≤ C0 b2n /2 and |ξ|0 ≤ s.
P
Let J = {2|Et> Xt |∞ ≤ λ0 }. Observe that |Et> Xt |∞ = maxj≤p | i∈Nt w(i, t)Xij ei | and, by
assumption 1,
!
X

0, σ 2

w(i, t)Xij ei ∼ N

X

w(i, t)2 Xij2

.

(37)

i∈Nt

i∈Nt

Then, by the standard Gaussian tail bound and the union bound, we obtain that
P


 2
p
p
ε
i∈Nt w(i, t)Xij ei
2
2
P max
≥ ε + 2 log p ≤ P(max |Zj | ≥ ε + 2 log p) ≤ 2 exp −
j≤p
j≤p
σLt,2
2
for all ε > 0, where Zj ∼ N (0, 1). Now, choose ε = (2 log p)1/2 and λ0 = 4σLt,2 (log p)1/2 ,
we have P(J ) ≥ 1 − 2p−1 . Further, we have
|β 0 (t)> Xt> Mt Xt [β̃(t) − β(t)]| ≤ |β̃(t) − β(t)|1 |Xt> Mt Xt β 0 (t)|∞
!1/2
X
 0 > > 2
1/2
w(i, t)Xij2
≤ |β̃(t) − β(t)|1 max
β (t) Xt Mt Xt β 0 (t)
j≤p

(Cauchy-Schwarz)

i∈Nt

q
≤ |β̃(t) − β(t)|1 Lt,1 ρmax (Xt> Mt2 Xt , s)|β 0 (t)|2 (assumption 2)
q
≤ |β̃(t) − β(t)|1 Lt,1 C0 s1/2 bn ρmax (Xt> Xt , s) (Lemma 7.1, assumption 2 and 3)
≤ |β̃(t) − β(t)|1 Lt,1 C0 (|Nt |wt s)1/2 bn ε−1
0

(Lemma 7.2, equation (34)).

Similarly, we can show that |ξ > Xt> Xt [β̃(t)−β(t)]| = O(Lt,1 (|Nt |wt s)1/2 b2n ε−1
0 |β̃(t)−β(t)|1 ).
Therefore, it follows that, with probability at least (1 − 2p−1 ),
D
E


Yt − Xt β(t), Xt [β̃(t) − β(t)] ≤ λ0 + 2Lt,1 C0 (|Nt |wt s)1/2 bn ε−1
0 (1 + o(1)) |β̃(t) − β(t)|1 .
Now, choose λ1 ≥ 2(λ0 + 2Lt,1 C0 (|Nt |wt s)1/2 bn ε−1
0 ), we get
2|Xt [β̃(t) − β(t)]|22 + 2λ1 |β̃(t)|1 ≤ λ1 |β̃(t) − β(t)|1 + 2λ1 |β(t)|1 .
Denote S0 := S0 (t) = supp(β(t)). By the same argument as [6, Lemma 6.3], it is easy to
see that, on J ,
2|Xt [β̃(t) − β(t)]|22 + λ1 |β̃ S0c (t)|1 ≤ 3λ1 |β̃ S0 (t) − β S0 (t)|1 .
But then, (35) follows from the restricted eigenvalue condition (assumption 5) with the
elementary inequality 4ab ≤ a2 + 4b2 that
2|Xt [β̃(t) − β(t)]|22 + λ1 |β̃(t) − β(t)|1 ≤ 4λ1 |β̃ S0 (t) − β S0 (t)|1 ≤ |Xt [β̃(t) − β(t)]|22 + 4λ21 s/φ20 .

28

Definition 7.1. A mean zero random variable is said to be sub-Gaussian with variance
factor σ 2 if
log E(eλX ) ≤ λ2 σ 2 /2

for all λ ∈ R.

Lemma 7.4. Let ξi be iid sub-Gaussian random variables with mean zero and variance
P∞
factor σ 2 , and ei =
m=0 am ξi−m be a linear process. Let w = (w1 , · · · , wn ) be a real
Pn
vector and Sn = i=1 wi ei be the weighted partial sum of ei .
1. (Short-range dependence). If |a|1 =

P∞

|ai | < ∞, then for all x > 0 we have


x2
P(|Sn | ≥ x) ≤ 2 exp −
.
(38)
2|w|22 |a|21 σ 2
i=0

2. (Long-range dependence). Suppose K = supm≥0 |am |(m + 1)% < ∞, where 1/2 < % <
1. Then, there exists a constant C% that only depends on % such that


C % x2
P(|Sn | ≥ x) ≤ 2 exp − 2 2(1−%) 2 2 .
|w|2 n
σ K
Proof. Put am = 0 if m < 0 and we may write Sn =

P

m∈Z bm ξm ,

where bm =

(39)
Pn

i=1

wi ai−m .

By the Cauchy-Schwarz inequality,
X

b2m ≤

m∈Z

X

n
X

m∈Z

i=1

!

n
X

wi2 |ai−m |

!
≤ |w|22 |a|21 .

|ai−m |

i=1

Then, (38) follows from the Cramér-Chernoff bound [4]. Let ām = maxl≥m |al | and Am =
Pn
Pm
−%
≤ C% K(n + 1)1−% , where C% = (1 − %)−1 . Then,
l=0 (l + 1)
l=0 |al |. Note that An ≤ K
we have

n
X

b2m ≤

m=1−n

n
X

n
X

m=1−n

i=1

!
wi2 |ai−m |

n
X

!
|ai−m |

≤ |w|22 A22n .

i=1

If m ≤ −n, then |bm | ≤ |w|1 ā1−m and therefore
X
m≤−n

b2m ≤ |w|21

X

ā21−m ≤ C% n|w|22 K 2 n1−2% ,

m≤−n

where the last inequality follows from Karamata’s theorem; see e.g. [26]. Hence, the proof is
complete by invoking the Cramér-Chernoff bound for sub-Gaussian random variables.

29

Some implementation issues
We assumed that the noise variance-covariance matrix Σe is known. In the iid error case
Σe = σ 2 In , we have seen that the distribution F (·) is independent of σ 2 and therefore
its value does not affect the inference procedure. The noise variance only impacts the
tuning parameter of the initial Lasso estimator. In practice, we can use the scaled Lasso
to estimate σ 2 in our numeric studies. Given that |σ̂/σ − 1| = oP (1) [32], the theoretical
properties of our estimator (10) remains the same if we plug in the scaled Lasso variance
output to our method. For temporally dependent stationary error process, estimation of Σe
becomes more subtle since it involves n autocovariance parameters. We propose a heuristic
strategy: first, run the tv-Lasso estimator and get the residuals; then calculate the sample
autocovariance matrix and apply a banding or tapering operation Bh (Σ) = {σjk 1(|j − k| ≤
h)}pj,k=1 [3, 7, 22].
We provide some justification on the heuristic strategy for SRD time series models. To
simplify explanation, we consider the uniform kernel and the bandwidth bn = 1. Suppose
we have an oracle where β(t) is known and we have access to the error process e(t). Let
Σ∗e be the oracle sample covariance matrix of ei with the Toeplitz structure i.e. the h-th
P
∗
subdiagonal of Σ∗e is σe,h
= n−1 n−h
i=1 ei ei+h . We first compare the oracle estimator and
the true error covariance matrix Σe . Let α > 0 and define
(
)
p
X
T (α, C1 , C2 ) = M ∈ ST p×p :
|mk | ≤ C1 h−α , ρj (M ) ∈ [C2 , C2−1 ], ∀j = 1, · · · , p ,
k=h+1

where ST p×p is the set of all p × p symmetric Toeplitz matrices. If ei has SRD, then
Σe ∈ T (% − 1, C1 , C2 ). By the argument in [3] and Lemma 7.4, we can show that
ρmax (Bh (Σ∗e ) − Σe )

≤
.P

ρmax (Bh (Σ∗e ) − Bh (Σe )) + ρmax (Bh (Σe ) − Σe )
r
log h
h
+ h−(%−1) .
n

Choosing h∗  (n/ log n)1/(2%) , we get
ρmax (Bh (Σ∗e ) − Σe ) = OP



log n
n

!
 %−1
2%
.

This oracle rate is sharper than the one established in [3] for regularizing more general
bandable matrices if n = o(p). Here, the improved rate is due to the Toeplitz structure
30

in Σe . Since Σe has uniformly bounded eigenvalues from zero and infinity, the banded
oracle estimator Bh (Σ∗e ) can be used as a benchmark to assess the tv-Lasso residuals
Ẽt = Yt − Xt β̃(t).
Proposition 7.5. Suppose Σe ∈ T (%−1, C1 , C2 ) and conditions of Lemma 7.3 are satisfied
except that (ei ) is an SRD stationary Gaussian process with % > 1. Then
ρmax (Bh (Σ̂e ) − Bh (Σ∗e )) = OP (hλ1 s1/2 ).

(40)

With the choice h∗  (n0 / log n0 )1/2% where n0 = |Nt |, we have

ρmax (Bh (Σ̂e ) − Σe )) = OP

log n
n0

0

 %−1
2%


+

0

n
log n0

 2%1

r

!!

s log p
+ sbn )
n0

.

(41)

It is interesting to note that the price we pay to choose h for not knowing the error
process is the second term in (41). Bandwidth selection for the smoothing parameter bn is
a theoretically challenging task in the high dimension. Asymptotic optimal order for the
parameter is available up to some unknown constants depending on the data generation
parameters. We shall use the cross-validation (CV) in our simulation studies and real data
analysis.
Proof of Proposition 7.5. Since we consider the uniform kernel, we may assume bn =
1, |Nt | = n and then rescale. Observe that
2
max |σ̂e,k
|k|≤h

−

∗2
σe,k
|

n−k
1 X
= max
(êi êi+k − ei ei+k )
|k|≤h n
i=1
n−k
n−k
X
1 X
êi (êi+k − ei+k ) +
ei+k (êi − ei )
|k|≤h n
i=1
i=1
!1/2 n−k
!1/2
n−k
X
1 X 2
≤ max
êi
(êi+k − ei+k )2
|k|≤h n
i=1
i=1
!1/2 n−k
!1/2
n−k
X
1 X 2
+ max
ei+k
(êi − ei )2
|k|≤h n
i=1
i=1

!1/2
!1/2 
!1/2
n
n
n
X
X
X
1
1
1

≤ 
ê2
+
e2
(êi − ei )2
.
n i=1 i
n i=1 i
n i=1

≤ max

31

By Lemma 7.3,
n

1X
(êi − ei )2 = |Ẽt − Et |22 = |Xt [β̃(t) − β(t)]|22 = OP (λ21 s).
n i=1
Then, it follows from the last expression and n−1

Pn

2
i=1 ei

= OP (1) that

2
∗2
max |σ̂e,k
− σe,k
| = OP (λ1 s1/2 ).
|k|≤h

Therefore
∗2
2
| = OP (hλ1 s1/2 ).
− σe,k
ρmax (Bh (Σ̂e ) − Bh (Σ∗e )) . h max |σ̂e,k
|k|≤h

32

