SPLINE SINGLE-INDEX PREDICTION MODEL
Li Wang and Lijian Yang
University of Georgia and Michigan State University
Abstract: For the past two decades, single-index model, a special case of projection pursuit regres-
sion, has proven to be an efficient way of coping with the high dimensional problem in nonparamet-
ric regression. In this paper, based on weakly dependent sample, we investigate the single-index
prediction (SIP) model which is robust against deviation from the single-index model. The single-
index is identified by the best approximation to the multivariate prediction function of the response
variable, regardless of whether the prediction function is a genuine single-index function. A poly-
nomial spline estimator is proposed for the single-index prediction coefficients, and is shown to be
root-n consistent and asymptotically normal. An iterative optimization routine is used which is
sufficiently fast for the user to analyze large data of high dimension within seconds. Simulation
experiments have provided strong evidence that corroborates with the asymptotic theory. Appli-
cation of the proposed procedure to the rive flow data of Iceland has yielded superior out-of-sample
rolling forecasts.
Key words and phrases: B-spline, geometric mixing, knots, nonparametric regression, root-n rate,
strong consistency.
1. Introduction
XTi , Yi
= {Xi,1, ...,Xi,d, Yi}ni=1 be a length n realization of a (d+ 1)-dimensional
strictly stationary process following the heteroscedastic model
Yi = m (Xi) + σ (Xi) εi,m (Xi) = E (Yi|Xi) , (1.1)
in which E (εi |Xi ) = 0, E
ε2i |Xi
= 1, 1 ≤ i ≤ n. The d-variate functions m, σ are the
unknown mean and standard deviation of the response Yi conditional on the predictor vector
Xi, often estimated nonparametrically. In what follows, we let
XT , Y, ε
have the stationary
distribution of
XTi , Yi, εi
. When the dimension of X is high, one unavoidable issue is the
“curse of dimensionality”, which refers to the poor convergence rate of nonparametric esti-
mation of general multivariate function. Much effort has been devoted to the circumventing
of this difficulty. In the words of Xia, Tong, Li and Zhu (2002), there are essentially two
approaches: function approximation and dimension reduction. A favorite function approxima-
tion technique is the generalized additive model advocated by Hastie and Tibshirani (1990),
Address for correspondence: Lijian Yang, Department of Statistics and Probability, Michigan State Univer-
sity, East Lansing, MI 48824, USA. E-mail: yang@stt.msu.edu
http://arxiv.org/abs/0704.0302v2
2 LI WANG AND LIJIAN YANG
see also, for example, Mammen, Linton and Nielsen (1999), Huang and Yang (2004), Xue
and Yang (2006 a, b), Wang and Yang (2007). An attractive dimension reduction method is
the single-index model, similar to the first step of projection pursuit regression, see Friedman
and Stuetzle (1981), Hall (1989), Huber (1985), Chen (1991). The basic appeal of single-index
model is its simplicity: the d-variate function m (x) = m (x1, ..., xd) is expressed as a univariate
function of xT θ0 =
p=1 xpθ0,p. Over the last two decades, many authors had devised various
intelligent estimators of the single-index coefficient vector θ0 = (θ0,1, ..., θ0,d)
, for instance,
Powell, Stock and Stoker (1989), Härdle and Stoker (1989), Ichimura (1993), Klein and Spady
(1993), Härdle, Hall and Ichimura (1993), Horowitz and Härdle (1996), Carroll, Fan, Gijbels
and Wand (1997), Xia and Li (1999), Hristache, Juditski and Spokoiny (2001). More recently,
Xia, Tong, Li and Zhu (2002) proposed the minimum average variance estimation (MAVE) for
several index vectors.
All the aforementioned methods assume that the d-variate regression function m (x) is
exactly a univariate function of some xT θ0 and obtain a root-n consistent estimator of θ0. If
this model is misspecified (m is not a genuine single-index function), however, a goodness-of-fit
test then becomes necessary and the estimation of θ0 must be redefined, see Xia, Li, Tong
and Zhang (2004). In this paper, instead of presuming that underlying true function m is
a single-index function, we estimate a univariate function g that optimally approximates the
multivariate function m in the sense of
g (ν) = E
m (X)|XT θ0 = ν
, (1.2)
where the unknown parameter θ0 is called the SIP coefficient, used for simple interpretation
once estimated; XT θ0 is the latent SIP variable; and g is a smooth but unknown function used
for further data summary, called the link prediction function. Our method therefore is clearly
interpretable regardless of the goodness-of-fit of the single-index model, making it much more
relevant in applications.
We propose estimators of θ0 and g based on weakly dependent sample, which includes
many existing nonparametric time series models, that are (i) computationally expedient and
(ii) theoretically reliable. Estimation of both θ0 and g has been done via the kernel smoothing
techniques in existing literature, while we use polynomial spline smoothing. The greatest
advantages of spline smoothing, as pointed out in Huang and Yang (2004), Xue and Yang
(2006 b) are its simplicity and fast computation. Our proposed procedure involves two stages:
estimation of θ0 by some
n-consistent θ̂, minimizing an empirical version of the mean squared
error, R(θ) = E{Y − E(Y |XT θ)}2; spline smoothing of Y on XT θ̂ to obtain a cubic spline
estimator ĝ of g. The best single-index approximation to m(x) is then m̂(x) = ĝ
xT θ̂
Under geometrically strong mixing condition, strong consistency and
n-rate asymptotic
SINGLE-INDEX PREDICTION MODEL 3
normality of the estimator θ̂ of the SIP coefficient θ0 in (1.2) are obtained. Proposition 2.2 is
the key in understanding the efficiency of the proposed estimator. It shows that the derivatives
of the risk function up to order 2 are uniformly almost surely approximated by their empirical
versions.
Practical performance of the SIP estimators is examined via Monte Carlo examples. The
estimator of the SIP coefficient performs very well for data of both moderate and high dimension
d, of sample size n from small to large, see Tables 1 and 2, Figures 1 and 2. By taking advantages
of the spline smoothing and the iterative optimization routines, one reduces the computation
burden immensely for massive data sets. Table 2 reports the computing time of one simulation
example on an ordinary PC, which shows that for massive data sets, the SIP method is much
faster than the MAVE method. For instance, the SIP estimation of a 200-dimensional θ0 from
a data of size 1000 takes on average mere 2.84 seconds, while the MAVE method needs to spend
2432.56 seconds on average to obtain a comparable estimates. Hence on account of criteria (i)
and (ii), our method is indeed appealing. Applying the proposed SIP procedure to the rive flow
data of Iceland, we have obtained superior forecasts, based on a 9-dimensional index selected
by BIC, see Figure 5.
The rest of the paper is organized as follows. Section 2 gives details of the model spec-
ification, proposed methods of estimation and main results. Section 3 describes the actual
procedure to implement the estimation method. Section 4 reports our findings in an extensive
simulation study. The proposed SIP model and the estimation procedure are applied in Section
5 to the rive flow data of Iceland. Most of the technical proofs are contained in the Appendix.
2. The Method and Main Results
2.1. Identifiability and definition of the index coefficient
It is obvious that without constraints, the SIP coefficient vector θ0 = (θ0,1, ..., θ0,d)
identified only up to a constant factor. Typically, one requires that ‖θ0‖ = 1 which entails
that at least one of the coordinates θ0,1, ..., θ0,d is nonzero. One could assume without loss of
generality that θ0,d > 0, and the candidate θ0 would then belong to the upper unit hemisphere
Sd−1+ =
(θ1, ..., θd) |
p=1 θ
p = 1, θd > 0
For a fixed θ = (θ1, ..., θd)
, denote Xθ = X
T θ, Xθ,i = X
i θ, 1 ≤ i ≤ n. Let
mθ (Xθ) = E (Y |Xθ) = E {m (X) |Xθ} . (2.1)
Define the risk function of θ as
R (θ) = E
{Y −mθ (Xθ)}2
= E {m (X)−mθ (Xθ)}2 + Eσ2 (X) , (2.2)
4 LI WANG AND LIJIAN YANG
which is uniquely minimized at θ0 ∈ Sd−1+ , i.e.
θ0 = arg min
θ∈Sd−1+
R (θ) .
Remark 2.1. Note that Sd−1+ is not a compact set, so we introduce a cap shape subset of
Sd−1+
Sd−1c =
(θ1, ..., θd) |
θ2p = 1, θd ≥
1− c2
, c ∈ (0, 1)
Clearly, for an appropriate choice of c, θ0 ∈ Sd−1c , which we assume in the rest of the paper.
Denote θ−d = (θ1, ..., θd−1)
, since for fixed θ ∈ Sd−1+ , the risk function R (θ) depends only
on the first d− 1 values in θ, so R (θ) is a function of θ−d
R∗ (θ−d) = R
θ1, θ2, ..., θd−1,
1− ‖θ−d‖22
with well-defined score and Hessian matrices
S∗ (θ−d) =
R∗ (θ−d) , H
∗ (θ−d) =
∂θ−d∂θ
R∗ (θ−d) . (2.3)
Assumption A1: The Hessian matrix H∗ (θ0,−d) is positive definite and the risk function R
is locally convex at θ0,−d, i.e., for any ε > 0, there exists δ > 0 such that R
∗ (θ−d)−R∗ (θ0,−d) <
δ implies ‖θ−d − θ0,−d‖2 < ε.
2.2. Variable transformation
Throughout this paper, we denote by Bda =
x ∈ Rd |‖x‖ ≤ a
the d-dimensional ball with
radius a and center 0 and
∣the kth order partial derivatives of m are continuous on B
the space of k-th order smooth functions.
Assumption A2: The density function of X, f (x) ∈ C(4)
, and there are constants
0 < cf ≤ Cf such that
cf/Vold
≤ f (x) ≤ Cf/Vold
, x ∈ Bda
f (x) ≡ 0, x /∈ Bda
For a fixed θ, define the transformed variables of the SIP variable Xθ
Uθ = Fd (Xθ) , Uθ,i = Fd (Xθ,i) , 1 ≤ i ≤ n, (2.4)
SINGLE-INDEX PREDICTION MODEL 5
in which Fd is the a rescaled centered Beta {(d+ 1) /2, (d+ 1) /2} cumulative distribution
function, i.e.
Fd (ν) =
∫ ν/a
Γ (d+ 1)
Γ {(d+ 1) /2}2 2d
1− t2
)(d−1)/2
dt, ν ∈ [−a, a] . (2.5)
Remark 2.2. For any fixed θ, the transformed variable Uθ in (2.4) has a quasi-uniform [0, 1]
distribution. Let fθ (u) be the probability density function of Uθ, then for any u ∈ [0, 1]
fθ (u) =
d (v)
fXθ (v) , v = F
d (u) ,
in which fXθ (v) = lim△ν→0 P (ν ≤ Xθ ≤ ν +△ν). Noting that xθ is exactly the projection of
x on θ, let Dν = {x|ν ≤ xθ ≤ ν +△ν} ∩Bda, then one has
P (ν ≤ Xθ ≤ ν +△ν) = P (X ∈ Dν) =
f (x) dx.
According to Assumption A2
cfVold(Dν)
Vold (B
≤ P (ν ≤ Xθ ≤ ν +△ν) ≤
CfVold(Dν)
Vold (B
On the other hand
Vold(Dν) = Vold−1(Jν)△ν + o (△ν) ,
where Jν = {x|xθ = v} ∩Bda. Note that the volume of Bda is πd/2ad/Γ (d/2 + 1) and
Vold−1 (Jν) = π(d−1)/2
a2 − ν2
)(d−1)/2
Γ {(d+ 1)/2} ,
Vold−1(Jν)
Vold (B
Γ (d+ 1)
}(d−1)/2
Therefore 0 < cf ≤ fθ (u) ≤ Cf <∞, for any fixed θ and u ∈ [0, 1].
In terms of the transformed SIP variable Uθ in (2.4), we can rewrite the regression function
mθ in (2.1) for fixed θ
γθ (Uθ) = E {m (X) |Uθ} = E {m (X) |Xθ} = mθ (Xθ) , (2.6)
then the risk function R (θ) in (2.2) can be expressed as
R (θ) = E
{Y − γθ (Uθ)}2
= E {m (X)− γθ (Uθ)}2 + Eσ2 (X) . (2.7)
2.3. Estimation Method
6 LI WANG AND LIJIAN YANG
Estimation of both θ0 and g requires a degree of statistical smoothing, and all estimation
here is carried out via cubic spline. In the following, we define the estimator θ̂ of θ0 and the
estimator ĝ of g.
To introduce the space of splines, we pre-select an integer n1/6 ≪ N = Nn ≪ n1/5 (log n)−2/5,
see Assumption A6 below. Divide [0, 1] into (N + 1) subintervals Jj = [tj, tj+1), j = 0, ..., N −
1, JN = [tN , 1], where T := {tj}Nj=1 is a sequence of equally-spaced points, called interior knots,
given as
t1−k = ... = t−1 = t0 = 0 < t1 < ... < tN < 1 = tN+1 = ... = tN+k,
in which tj = jh, j = 0, 1, ..., N +1, h = 1/ (N + 1) is the distance between neighboring knots.
The j-th B-spline of order k for the knot sequence T denoted by Bj,k is recursively defined by
de Boor (2001).
Denote by Γ(k−2) = Γ(k−2) [0, 1] the space of all C(k−2) [0, 1] functions that are polynomials
of degree k−1 on each interval. For fixed θ, the cubic spline estimator γ̂θ of γθ and the related
estimator m̂θ of mθ are defined as
γ̂θ (·) = arg min
γ(·)∈Γ(2)[0,1]
{Yi − γ (Uθ,i)}2 , m̂θ (ν) = γ̂θ {Fd (ν)} . (2.8)
Define the empirical risk function of θ
R̂ (θ) = n−1
{Yi − γ̂θ (Uθ,i)}2 = n−1
{Yi − m̂θ (Xθ,i)}2 , (2.9)
then the spline estimator of the SIP coefficient θ0 is defined as
θ̂ = arg min
θ∈Sd−1c
R̂ (θ) ,
and the cubic spline estimator of g is m̂θ with θ replaced by θ̂, i.e.
ĝ (ν) =
arg min
γ(·)∈Γ(2)[0,1]
Yi − γ
{Fd (ν)} . (2.10)
2.4. Asymptotic results
Before giving the main theorems, we state some other assumptions.
Assumption A3: The regression function m ∈ C(4)
for some a > 0.
Assumption A4: The noise ε satisfies E (ε |X) = 0, E
ε2 |X
= 1 and there exists a positive
constant M such that sup
|ε|3 |X = x
< M . The standard deviation function σ (x) is
continuous on Bda,
0 < cσ ≤ inf
x∈Bda
σ (x) ≤ sup
x∈Bda
σ (x) ≤ Cσ <∞.
SINGLE-INDEX PREDICTION MODEL 7
Assumption A5: There exist positive constants K0 and λ0 such that α (n) ≤ K0e−λ0n holds
for all n, with the α-mixing coefficient for
XTi , εi
defined as
α (k) = sup
B∈σ{Zs,s≤t},C∈σ{Zs,s≥t+k}
|P (B ∩ C)− P (B)P (C)| , k ≥ 1.
Assumption A6: The number of interior knots N satisfies: n1/6 ≪ N ≪ n1/5 (log n)−2/5.
Remark 2.3. Assumptions A3 and A4 are typical in the nonparametric smoothing literature,
see for instance, Härdle (1990), Fan and Gijbels (1996), Xia, Tong Li and Zhu (2002). By
the result of Pham (1986), a geometrically ergodic time series is a strongly mixing sequence.
Therefore, Assumption A5 is suitable for (1.1) as a time series model under aforementioned
assumptions.
We now state our main results in the next two theorems.
Theorem 1. Under Assumptions A1-A6, one has
θ̂−d−→ θ0,−d, a.s.. (2.11)
Proof. Denote by (Ω,F ,P) the probability space on which all
XTi , Yi
are defined. By
Proposition 2.2, given at the end of this section
‖θ−d‖2≤
R̂∗ (θ−d)−R∗ (θ−d)
−→ 0, a.s.. (2.12)
So for any δ > 0 and ω ∈ Ω, there exists an integer n0 (ω), such that when n > n0 (ω),
R̂∗ (θ0,−d, ω) − R∗ (θ0,−d) < δ/2. Note that θ̂−d = θ̂−d (ω) is the minimizer of R̂∗ (θ−d, ω),
so R̂∗
θ̂−d (ω) , ω
− R∗ (θ0,−d) < δ/2. Using (2.12), there exists n1 (ω), such that when
n > n1 (ω), R
θ̂−d (ω) , ω
− R̂∗
θ̂−d (ω) , ω
< δ/2. Thus, when n > max (n0 (ω) , n1 (ω)),
θ̂−d (ω) , ω
−R∗ (θ0,−d) < δ/2 + R̂∗
θ̂−d (ω) , ω
−R∗ (θ0,−d) < δ/2 + δ/2 = δ.
According to Assumption A1, R∗ is locally convex at θ0,−d, so for any ε > 0 and any ω, if
θ̂−d (ω) , ω
−R∗ (θ0,−d) < δ, then
∥θ̂−d (ω)−θ0,−d
∥ < ε for n large enough , which implies
the strong consistency.
Theorem 2. Under Assumptions A1-A6, one has
θ̂−d−θ0,−d
d−→ N {0,Σ (θ0)} ,
where Σ (θ0) = {H∗ (θ0,−d)}−1Ψ(θ0) {H∗ (θ0,−d)}−1, H∗ (θ0,−d) = {lpq}d−1p,q=1 and Ψ(θ0) =
{ψpq}d−1p,q=1 with
lp,q = −2E [{γ̇pγ̇q + γθ0 γ̈p,q} (Uθ0)] + 2θ0,qθ
E [{γ̇pγ̇d (Uθ0) + γθ0 γ̈p,d} (Uθ0)]
+2θ−30,dE [(γθ0 γ̇d) (Uθ0)]
θ20,d + θ
I{p=q} + θ0,pθ0,qI{p 6=q}
+2θ0,pθ
0,dE [{γ̇pγ̇q + γθ0 γ̈p,q} (Uθ0)]− 2θ0,pθ0,qθ
γ̇2d + γθ0 γ̈d,d
(Uθ0)
8 LI WANG AND LIJIAN YANG
ψpq = 4E
γ̇p − θ0,pθ−10,dγ̇d
γ̇q − θ0,qθ−10,dγ̇d
(Uθ0) {γθ0 (Uθ0)− Y }
in which γ̇p and γ̈p,q are the values of
∂θp∂θq
γθ taking at θ = θ0, for any p, q = 1, 2, ..., d−1
and γθ is given in (2.6).
Remark 2.4. Consider the Generalized Linear Model (GLM): Y = g
XT θ0
+σ (X) ε, where
g is a known link function. Let θ̃ be the nonlinear least squared estimator of θ0 in GLM.
Theorem 2 shows that under the assumptions A1-A6, the asymptotic distribution of the θ̂−d
is the same as that of θ̃. This implies that our proposed SIP estimator θ̂−d is as efficient as if
the true link function g is known.
The next two propositions play an important role in our proof of the main results. Propo-
sition 2.1 establishes the uniform convergence rate of the derivatives of γ̂θ up to order 2 to
those of γθ in θ. Proposition 2.2 shows that the derivatives of the risk function up to order 2
are uniformly almost surely approximated by their empirical versions.
Proposition 2.1. Under Assumptions A2-A6, with probability 1
θ∈Sd−1c
u∈[0,1]
|γ̂θ (u)− γθ (u)| = O
log n+ h4
, (2.13)
1≤p≤d
θ∈Sd−1c
1≤i≤n
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
log n√
, (2.14)
1≤p,q≤d
θ∈Sd−1c
1≤i≤n
∂θp∂θq
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
log n√
. (2.15)
Proposition 2.2. Under Assumptions A2-A6, one has for k = 0, 1, 2
‖θ−d‖≤
∂kθ−d
R̂∗ (θ−d)−R∗ (θ−d)
= o(1), a.s..
Proofs of Theorem 2, Propositions 2.1 and 2.2 are given in Appendix.
3. Implementation
In this section, we will describe the actual procedure to implement the estimation of θ0
and g. We first introduce some new notation. For fixed θ, write the B-spline matrix as
Bθ = {Bj,4 (Uθ,i)}n, Ni=1,j=−3 and
Pθ = Bθ
θ (3.1)
as the projection matrix onto the cubic spline space Γ
n,θ. For any p = 1, ..., d, denote
Ḃp =
Bθ, Ṗp =
SINGLE-INDEX PREDICTION MODEL 9
as the first order partial derivatives of Bθ and Pθ with respect to θ.
Let Ŝ∗(θ−d) be the score vector of R̂
∗ (θ−d), i.e.
Ŝ∗(θ−d) =
R̂∗ (θ−d) . (3.2)
The next lemma provides the exact forms of Ŝ∗(θ−d).
Lemma 3.1. For the score vector of R̂∗ (θ−d) defined in (3.2), one has
Ŝ∗ (θ−d) = −n−1
ṖpY − θpθ−1d Y
, (3.3)
where for any p = 1, 2, ..., d
ṖpY = 2Y
T (I−Pθ) Ḃp
θ Y, (3.4)
where Ḃp =
{Bj,3 (Uθ,i)−Bj+1,3 (Uθ,i)} Ḟd (Xθ,i)h−1Xi,p
}n, N
i=1,j=−3
Ḟd (x) =
Γ (d+ 1)
aΓ {(d+ 1) /2}2 2d
I (|x| ≤ a) .
Proof. For any p = 1, 2, ..., d, the derivatives of B-splines in de Boor (2001) implies
Ḃp =
Bj,4 (Uθ,i)
}n, N
i=1,j=−3
Bj,4 (Uθ,i)
}n, N
i=1,j=−3
Bj,3 (Uθ,i)
tj+3 − tj
Bj+1,3 (Uθ,i)
tj+4 − tj+1
Ḟd (Xθ,i)Xi,p
}n, N
i=1,j=−3
{Bj,3 (Uθ,i)−Bj+1,3 (Uθ,i)} Ḟd (Xθ,i) h−1Xi,p
}n, N
i=1,j=−3
Next, note that
Ṗp = Ḃp
θ +Bθ
= Ḃp
θ +Bθ
θ +Bθ
Since
BTθ Bθ
BTθ Bθ
BTθ Bθ
θ Bθ +
)−1 ∂
BTθ Bθ
and ∂
BTθ Bθ
= ḂTp Bθ +B
θ Ḃp, thus
)−1 (
p Bθ +B
θ Ḃp
10 LI WANG AND LIJIAN YANG
Hence
Ṗp = (I−Pθ) Ḃp
θ +Bθ
p (I−Pθ) .
Thus, (3.4) follows immediately.
In practice, the estimation is implemented via the following procedure.
Step 1. Standardize the predictor vectors {Xi}ni=1 and for each fixed θ ∈ Sd−1c obtain the
CDF transformed variables {Uθ,i}ni=1 of the SIP variable {Xθ,i}
through formula (2.5), where
the radius a is taken to be the 95% percentile of {‖Xi‖}ni=1.
Step 2. Compute quadratic and cubic B-spline basis at each value Uθ,i, where the number
of interior knots N is
N = min
n1/5.5
, (3.5)
Step 3. Find the estimator θ̂ of θ0 by minimizing R̂
∗ through the port optimization routine
with (0, 0, ..., 1)
as the initial value and the empirical score vector Ŝ∗ in (3.3). If d < n, one
can take the simple LSE (without the intercept) for data {Yi,Xi}ni=1 with its last coordinate set
positive.
Step 4. Obtain the spline estimator ĝ of g by plugging θ̂ obtained in Step 3 into (2.10).
Remark 3.1. In (3.5), c1 and c2 are positive integers and [ν] denotes the integer part of ν. The
choice of the tuning parameter c1 makes little difference for a large sample and according to our
asymptotic theory there is no optimal way to set these constants. We recommend using c1 = 1
to save computing for massive data sets. The first term ensures Assumption A6. The addition
constrain c2 can be taken from 5 to 10 for smooth monotonic or smooth unimodel regression
and c2 > 10 if has many local minima and maxima, which is very unlikely in application.
4. Simulations
In this section, we carry out two simulations to illustrate the finite-sample behavior of our
SIP estimation method. The number of interior knots N is computed according to (3.5) with
c1 = 1, c2 = 5. All of our codes have been written in R.
Example 1. Consider the model in Xia, Li, Tong and Zhang (2004)
Y = m (X) + σ0ε, σ0 = 0.3, 0.5, ε
i.i.d∼ N(0, 1)
where X = (X1,X2)
T ∼N(0, I2), truncated by [−2.5, 2.5]2 and
m (x) = x1 + x2 + 4exp
− (x1 + x2)2
x21 + x
. (4.1)
If δ = 0, then the underlying true functionm is a single-index function, i.e., m (X) =
2XT θ0+
XT θ0
, where θT0 = (1, 1) /
2. While δ 6= 0, then m is not a genuine single-index
SINGLE-INDEX PREDICTION MODEL 11
function. An impression of the bivariate function m for δ = 0 and δ = 1 can be gained in
Figure 1 (a) and (b), respectively.
Table 1: Report of Example 1 (Values out/in parentheses: δ = 0/δ = 1)
σ0 n θ0 BIAS SD MSE Average MSE
5e− 04 0.00825 7e− 05
(−0.00236) (0.02093) (0.00044) 7e− 05
−6e− 04 0.00826 7e− 05 (0.00043)
(0.00174) (0.02083) (0.00043)
−0.00124 0.00383 2e− 05
(−0.00129) (0.01172) (0.00014) 2e− 05
−0.00124 0.00383 2e− 05 (0.00014)
(0.00110) (0.01160) (0.00013)
0.00121 0.01346 0.00018
(−0.00137) (0.02257) (0.00051) 0.00018
−0.00147 0.01349 0.00018 (0.00051)
(0.00062) (0.02309) (0.00052)
−0.00204 0.00639 4e− 05
(−0.00229) (0.01205) (0.00015) 4e− 05
0.00197 0.00637 4e− 05 (0.00015)
(0.00208) (0.01190) (0.00014)
For δ = 0, 1, we draw 100 random realizations of each sample size n = 50, 100, 300 respec-
tively. To demonstrate how close our SIP estimator is to the true index parameter θ0, Table 1
lists the sample mean (MEAN), bias (BIAS), standard deviation (SD), the mean squared error
(MSE) of the estimates of θ0 and the average MSE of both directions. From this table, we find
that the SIP estimators are very accurate for both cases δ = 0 and δ = 1, which shows that
our proposed method is robust against the deviation from single-index model. As we expected,
when the sample size increases, the SIP coefficient is more accurately estimated. Moreover, for
n = 100, 300, the total average is inversely proportional to n.
Example 2. Consider the heteroscedastic regression model (1.1) with
m (X) = sin
, σ (X) = σ0
5− exp
5 + exp
) , (4.2)
in which Xi = {Xi,1, ...,Xi,d}T and εi, i = 1, ..., n, are
i.i.d∼ N (0, 1), σ0 = 0.2. In our simulation,
the true parameter θT0 = (1, 1, 0, ..., 0, 1)/
3 for different sample size n and dimension d. The
12 LI WANG AND LIJIAN YANG
superior performance of SIP estimators is borne out in comparison with MAVE of Xia, Tong,
Li and Zhu (2002). We also investigate the behavior of SIP estimators in the previously
unemployed cases that sample size n is smaller than or equal to d, for instance, n = 100, d =
100, 200 and n = 200, d = 200, 400. The average MSEs of the d dimensions are listed in Table
2, from which we see that the performance of the SIP estimators are quite reasonable and in
most of the scenarios n ≤ d, the SIP estimators still work astonishingly well where the MAVEs
become unreliable. For n = 100, d = 10, 50, 100, 200, the estimates of the link prediction
function g from model (4.2) are plotted in Figure 2, which is rather satisfactory even when
dimension d exceeds the sample size n.
Theorem 1 indicates that θ̂−d is strongly consistent of θ0,−d. To see the convergence, we
run 100 replications and in each replication, the value of ‖θ̂ − θ0‖/
d is computed. Figure
3 plots the kernel density estimations of the 100 ‖θ̂ − θ0‖ in Example 2, in which dimension
d = 10, 50, 100, 200. There are four types of line characteristics which correspond to the two
sample sizes, the dotted-dashed line (n = 100), dotted line (n = 200), dashed line (500) and
solid line (n = 1000). As sample sizes increasing, the squared errors are becoming closer to 0,
with narrower spread out, confirmative to the conclusions of Theorem 1.
Lastly, we report the average computing time of Example 2 to generate one sample of size
n and perform the SIP or MAVE procedure done on the same ordinary Pentium IV PC in
Table 2. From Table 2, one sees that our proposed SIP estimator is much faster than the
MAVE. The computing time for MAVE is extremely sensitive to sample size as we expected.
For very large d, MAVE becomes unstable to the point of the breaking down in four cases.
5. An application
In this section we demonstrate the proposed SIP model through the river flow data of
Jökulsá Eystri River of Iceland, from January 1, 1972 to December 31, 1974. There are 1096
observations, see Tong (1990). The response variables are the daily river flow (Yt), measured in
meter cubed per second of Jökulsá Eystri River. The exogenous variables are temperature (Xt)
in degrees Celsius and daily precipitation (Zt) in millimeters collected at the meteorological
station at Hveravellir.
This data set was analyzed earlier through threshold autoregressive (TAR) models by
Tong, Thanoon and Gudmundsson (1985), Tong (1990), and nonlinear additive autoregressive
(NAARX) models by Chen and Tsay (1993). Figure 4 shows the plots of the three time series,
from which some nonlinear and non-stationary features of the river flow series are evident. To
make these series stationary, we remove the trend by a simple quadratic spline regression and
these trends (dashed lines) are shown in Figure 4. By an abuse of notation, we shall continue
to use Xt, Yt, Zt to denote the detrended series.
SINGLE-INDEX PREDICTION MODEL 13
In the analysis, we pre-select all the lagged values in the last 7 days (1 week), i.e., the
predictor pool is {Yt−1, ..., Yt−7,Xt,Xt−1, ...,Xt−7, Zt, Zt−1, ..., Zt−7, }. Using BIC similar to
Huang and Yang (2004) for our proposed spline SIP model with 3 interior knots, the following
9 explanatory variables are selected from the above set {Yt−1, ..., Yt−4,Xt,Xt−1,Xt−2, Zt, Zt−1}.
Based on this selection, we fit the SIP model again and obtain the estimate of the SIP coefficient
θ̂ = {−0.877, 0.382,−0.208, 0.125,−0.046,−0.034, 0.004,−0.126, 0.079}T . Figure 5 (a) and (b)
display the fitted river flow series and the residuals against time.
Next we examine the forecasting performance of the SIP method. We start with estimating
the SIP estimator using only observations of the first two years, then we perform the out-of-
sample rolling forecast of the entire third year. The observed values of the exogenous variables
are used in the forecast. Figure 5 (c) shows this SIP out-of-sample rolling forecasts. For the
purpose of comparison, we also try the MAVE method, in which the same predictor vector is
selected by using BIC. The mean squared prediction error is 60.52 for the SIP model, 61.25 for
MAVE, 65.62 for NAARX, 66.67 for TAR and 81.99 for the linear regression model, see Chen
and Tsay (1993). Among the above five models, the SIP model produces the best forecasts.
6. Conclusion
In this paper we propose a robust SIP model for stochastic regression under weak depen-
dence regardless if the underlying function is exactly a single-index function. The proposed
spline estimator of the index coefficient possesses not only the usual strong consistency and
n-rate asymptotically normal distribution, but also is as efficient as if the true link function
g is known. By taking advantage of the spline smoothing method and the iterative method,
the proposed procedure is much faster than the MAVE method. This procedure is especially
powerful for large sample size n and high dimension d and unlike the MAVE method, the
performance of the SIP remains satisfying in the case d > n.
Acknowledgment
This work is part of the first author’s dissertation under the supervision of the second
author, and has been supported in part by NSF award DMS 0405330.
Appendix
A.1. Preliminaries
In this section, we introduce some properties of the B-spline.
Lemma A.1. There exist constants c > 0 such that for
j=−k+1 αj,kBj,k up to order k = 4
ch1/r ‖α‖r ≤
j=−k+1 αj,kBj,k
3r−1h
)1/r ‖α‖r , 1 ≤ r ≤ ∞
ch1/r ‖α‖r ≤
j=−k+1 αj,kBj,k
≤ (3h)1/r ‖α‖r , 0 < r < 1
14 LI WANG AND LIJIAN YANG
where α := (α−1,2, α0,2, ..., αN,2, ..., αN,4). In particular, under Assumption A2, for any fixed θ
ch1/2 ‖α‖2 ≤
j=−k+1
αj,kBj,k
≤ Ch1/2 ‖α‖2 .
Proof. It follows from the B-spline property on page 96 of de Boor (2001),
j=−k+1Bj,k ≡
3 on [0, 1]. So the right inequality follows immediate for r = ∞. When 1 ≤ r < ∞, we use
Hölder’s inequality to find
j=−k+1
αj,kBj,k
j=−k+1
|αj,k|r Bj,k
j=−k+1
1−1/r
= 31−1/r
j=−k+1
|αj,k|r Bj,k
Since all the knots are equally spaced,
−∞Bj,k (u) du ≤ h, the right inequality follows from
j=−k+1
αj,kBj,k (u)
du ≤ 3r−1h ‖α‖rr .
When r < 1, we have
j=−k+1
αj,kBj,k
j=−k+1
|αj,k|r Brj,k.
Since
j,k (u) du ≤ tj+k − tj = kh and
j=−k+1
αj,kBj,k (u)
du ≤ ‖α‖rr
Brj,k (u) du ≤ 3h ‖α‖
the right inequality follows in this case as well. For the left inequalities, we derive from Theorem
5.4.2, DeVore and Lorentz (1993)
|αj,k| ≤ C1h−1/r
∫ tj+1
j=−k+1
αj,kBj,k (u)
for any 0 < r ≤ ∞, so
|αj,k|r ≤ Cr1h−1
∫ tj+1
j=−k+1
αj,kBj,k (u)
SINGLE-INDEX PREDICTION MODEL 15
Since each u ∈ [0, 1] appears in at most k intervals (tj,tj+k), adding up these inequalities, we
obtain that
‖α‖rr ≤ C1h
∫ tj+k
j=−k+1
αj,kBj,k (u)
du ≤ 3Ch−1
j=−k+1
αj,kBj,k
The left inequality follows.
For any functions φ and ϕ, define the empirical inner product and the empirical norm as
〈φ,ϕ〉θ =
φ (u)ϕ (u) fθ (u) du, ‖φ‖22,n,θ = n
φ2 (Uθ,i) .
In addition, if functions φ,ϕ are L2 [0, 1]-integrable, define the theoretical inner product and
its corresponding theoretical L2 norm as
‖φ‖22,θ =
φ2 (u) fθ (u) du, 〈φ,ϕ〉n,θ = n
φ (Uθ,i)ϕ (Uθ,i) .
Lemma A.2. Under Assumptions A2, A5 and A6, with probability 1,
θ∈Sd−1c
k,k′=2,3,4
1≤j,j′≤N
Bj,k, Bj′,k′
Bj,k, Bj′,k′
log n
Proof. We only prove the case k = k′ = 4, all other cases are similar. Let
ζθ,j,j′,i = Bj,4 (Uθ,i)Bj′,4 (Uθ,i)− EBj,4 (Uθ,i)Bj′,4 (Uθ,i) ,
with the second moment
Eζ2θ,j,j′,i = E
B2j,4 (Uθ,i)B
j′,4 (Uθ,i)
EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
where
EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
}2 ∼ N−2, E
B2j,4 (Uθ,i)B
j′,4 (Uθ,i)
∼ N−1 by Assumption A2.
Hence, Eζ2θ,j,j′,i ∼ N−1. The k-th moment is given by
∣ζθ,j,j′,i
∣Bj,4 (Uθ,i)Bj′,4 (Uθ,i)− EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
≤ 2k−1
∣Bj,4 (Uθ,i)Bj′,4 (Uθ,i)
∣EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
where
∣EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
k ∼ N−k, E
∣EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
k ∼ N−1. Thus, there exists
a constant C > 0 such that E
∣ζθ,j,j′,i
k ≤ C2k−1k!Eζ2j,j′,i. So the Cramér’s condition is satisfied
with Cramér’s constant c∗. By the Bernstein’s inequality (see Bosq (1998), Theorem 1.4, page
31), we have for k = 3
ζθ,j,j′,i
≤ a1 exp
25m22 + 5c
+ a2 (k)α
q + 1
])6/7
16 LI WANG AND LIJIAN YANG
where
δn = δ
log n√
, a1 = 2
δ2 (nN)
log2 n
25m22 + 5c
, m22 ∼ N−1,
a2 (3) = 11n
, m3 = max
1≤i≤n
∥ζθ,j,j′,i
≤ cN1/3.
Observe that 5cδn = o(1) by Assumption A6, then by taking q such that
≥ c0 log n,
q ≥ c1n/ log n for some constants c0, c1, one has a1 = O(n/q) = O (log n), a2 (3) = o
Assumption A6 again. Assumption A5 yields that
q + 1
])6/7
K0 exp
q + 1
])}6/7
≤ Cn−6λ0c0/7.
Thus, for fixed θ ∈ Sd−1c , when n large enough
ζθ,j,j′,i
≤ c log n exp
−c2δ2 log n
+ Cn2−6λ0c0/7. (A.1)
We divide each range of θp, p = 1, 2, ..., d − 1, into n6/(d−1) equally spaced intervals with
disjoint endpoints −1 = θp,0 < θp,1 < ... < θp,Mn = 1, for p = 1, ..., d − 1. Projecting these
small cylinders onto Sd−1c , the radius of each patch Λr, r = 1, ...,Mn is bounded by cM
Denote the projection of the Mn points as θr =
θr,−d,
1− ‖θr,−d‖22
, r = 0, 1, ...,Mn.
Employing the discretization method, sup
θ∈Sd−1c
1≤j,j′≤N
∣ζθ,j,j′,i
∣ is bounded by
0≤r≤Mn
1≤j,j′≤N
∣ζθr,j,j′,i
∣+ sup
0≤r≤Mn
1≤j,j′≤N
∣ζθ,j,j′,i − ζθr,j,j′,i
∣ . (A.2)
By (A.1) and Assumption A6, there exists large enough value δ > 0 such that
ζθr,j,j′,i
≤ n−10,
which implies that
1≤j,j′≤N
ζθr ,j,j′,i
N2Mnn
−10 ≤ C
n−3 <∞.
Thus, Borel-Cantelli Lemma entails that
0≤r≤Mn
1≤j,j′≤N
ζθr ,j,j′,i
log n√
, a.s.. (A.3)
SINGLE-INDEX PREDICTION MODEL 17
Employing Lipschitz continuity of the cubic B-spline, one has with probability 1
0≤r≤Mn
1≤j,j′≤N
ζθ,j,j′,i − ζθr,j,j′,i
M−1n h
−6) . (A.4)
Therefore Assumption A2, (A.2), (A.3) and (A.4) lead to the desired result.
Denote by Γ = Γ(0)∪Γ(1)∪Γ(2) the space of all linear, quadratic and cubic spline functions
on [0, 1]. We establish the uniform rate at which the empirical inner product approximates the
theoretical inner product for all B-splines Bj,k with k = 2, 3, 4.
Lemma A.3. Under Assumptions A2, A5 and A6, one has
An = sup
θ∈Sd−1c
γ1,γ2∈Γ
〈γ1, γ2〉n,θ − 〈γ1, γ2〉θ
‖γ1‖2,θ ‖γ2‖2,θ
log n
, a.s.. (A.5)
Proof. Denote without loss of generality,
j=−k+1
αjkBj,k, γ2 =
j=−k+1
βjkBj,k,
for any two 3 (N + 3)-vectors
α =(α−1,2, α0,2, ..., αN,2, ..., αN,4) , β =(β−1,2, β0,2, ..., βN,2, ..., βN,4) .
Then for fixed θ
〈γ1, γ2〉n,θ =
j=−k+1
αj,kBj,k (Uθ,i)
j=−k+1
βj,kBj,k (Uθ,i)
j=−k+1
j′=−k+1
αj,kβj′,k′
Bj,k, Bj′,k′
‖γ1‖22,θ =
j=−k+1
j′=−k+1
αj,kαj′,k′
Bj,k, Bj′ ,k′
‖γ2‖22,θ =
j=−k+1
j′=−k+1
βj,kβj′,k′
Bj,k, Bj′ ,k′
According to Lemma A.1, one has for any θ ∈ Sd−1c ,
c1h ‖α‖22 ≤ ‖γ1‖
2,θ ≤ c2h ‖α‖
2 , c1h ‖β‖
2 ≤ ‖γ2‖
2,θ ≤ c2h ‖β‖
c1h ‖α‖2 ‖β‖2 ≤ ‖γ1‖2,θ ‖γ2‖2,θ ≤ c2h ‖α‖2 ‖β‖2 .
18 LI WANG AND LIJIAN YANG
Hence
An = sup
θ∈Sd−1c
γ1∈γ,γ2∈Γ
〈γ1, γ2〉n,θ − 〈γ1, γ2〉θ
‖γ1‖2,θ ‖γ2‖2,θ
‖α‖∞ ‖β‖∞
c1h ‖α‖2 ‖β‖2
× sup
θ∈Sd−1c
k,k′=2,3,4
1≤j,j′≤N
Bj,k, Bj′ ,k′
Bj,k, Bj′ ,k′
An ≤ c0h−1 sup
θ∈Sd−1c
k,k′=2,3,4
1≤j,j′≤N
Bj,k, Bj′ ,k′
Bj,k, Bj′ ,k′
which, together with Lemma A.2, imply (A.5).
A.2. Proof of Proposition 2.1
For any fixed θ, we write the response YT = (Y1, ..., Yn) as the sum of a signal vector γθ,
a parametric noise vector Eθ and a systematic noise vector E, i.e.,
Y = γθ +Eθ +E,
in which the vectors γTθ = {γθ (Uθ,1) , ..., γθ (Uθ,n)}, ET = {σ (X1) ε1, ..., σ (Xn) εn} and ETθ =
{m (X1)− γθ (Uθ,1) , ...,m (Xn)− γθ (Uθ,n)}.
Remark A.1. If m is a genuine single-index function, then Eθ0 ≡ 0, thus the proposed SIP
model is exactly the single-index model.
Let Γ
n, θ be the cubic spline space spanned by {Bj,4 (Uθ,i)}
, −3 ≤ j ≤ N for fixed θ.
Projecting Y onto Γ
n, θ yields that
γ̂θ = {γ̂θ (Uθ,1) , ..., γ̂θ (Uθ,n)}T = ProjΓ(2)
γθ + ProjΓ(2)
Eθ + ProjΓ(2)
where γ̂θ is given in (2.8). We break the cubic spline estimation error γ̂θ (uθ) − γθ (uθ) into a
bias term γ̃θ (uθ)− γθ (uθ) and two noise terms ε̃θ (uθ) and ε̂θ (uθ)
γ̂θ (uθ)− γθ (uθ) = {γ̃θ (uθ)− γθ (uθ)}+ ε̃θ (uθ) + ε̂θ (uθ) , (A.6)
where
γ̃θ (u) = {Bj,4 (u)}T−3≤j≤N V
〈γθ, Bj,4〉n,θ
, (A.7)
ε̃θ (u) = {Bj,4 (u)}T−3≤j≤N V
〈Eθ, Bj,4〉n,θ
, (A.8)
ε̂θ (u) = {Bj,4 (u)}T−3≤j≤N V
〈E, Bj,4〉n,θ
. (A.9)
SINGLE-INDEX PREDICTION MODEL 19
In the above, we denote by Vn,θ the empirical inner product matrix of the cubic B-spline basis
and similarly, the theoretical inner product matrix as Vθ
Vn,θ =
θ Bθ =
Bj′,4, Bj,4
j,j′=−3
,Vθ =
Bj′,4, Bj,4
j,j′=−3 . (A.10)
In Lemma A.5, we provide the uniform upper bound of
∞. Before
that, we first describe a special case of Theorem 13.4.3 in DeVore and Lorentz (1993).
Lemma A.4. If a bi-infinite matrix with bandwidth r has a bounded inverse A−1 on l2 and
κ = κ (A) := ‖A‖2
is the condition number of A, then
∞ ≤ 2c0 (1− ν)
, with
c0 = ν
−2r ∥
, ν =
κ2 − 1
)1/4r (
κ2 + 1
)−1/4r
Lemma A.5. Under Assumptions A2, A5 and A6, there exist constants 0 < cV < CV such
that cVN
−1 ‖w‖22 ≤ wTVθw ≤ CVN−1 ‖w‖
2 and
−1 ‖w‖22≤ w
Vn,θw ≤ CVN−1 ‖w‖22 , a.s., (A.11)
with matrices Vθ and Vn,θ defined in (A.10). In addition, there exists a constant C > 0 such
θ∈Sd−1c
≤ CN, a.s., sup
θ∈Sd−1c
∞ ≤ CN. (A.12)
Proof. First we compute the lower and upper bounds for the eigenvalues of Vn,θ. Let w be any
(N + 4)-vector and denote γw (u) =
j=−3wjBj,4 (u), then Bθw = {γw (Uθ,1) , ..., γw (Uθ,n)}
and the definition of An in (A.5) from Lemma A.3 entails that
‖γw‖22,θ (1−An) ≤ w
Vn,θw = ‖γw‖22,n,θ ≤ ‖γw‖
2,θ (1 +An) . (A.13)
Using Theorem 5.4.2 of DeVore and Lorentz (1993) and Assumption A2, one obtains that
‖w‖22 ≤ ‖γw‖
2,θ = w
Vθw =
wjBj,4
‖w‖22 , (A.14)
which, together with (A.13), yield
−1 ‖w‖22 (1−An) ≤ w
Vn,θw ≤ CfCN−1 ‖w‖22 (1 +An) . (A.15)
Now the order of An in (A.5), together with (A.14) and (A.15) implies (A.11), in which cV =
cfC,CV = CfC. Next, denote by λmax (Vn,θ) and λmin (Vn,θ) the maximum and minimum
eigenvalue of Vn,θ, simple algebra and (A.11) entail that
−1 ≥ ‖Vn,θ‖2 = λmax (Vn,θ) ,
= λ−1min (Vn,θ) ≤ c
V N, a.s.,
20 LI WANG AND LIJIAN YANG
κ := ‖Vn,θ‖2
= λmax (Vn,θ)λ
min (Vn,θ) ≤ CV c
V <∞, a.s..
Meanwhile, let wj = the (N + 4)-vector with all zeros except the j-th element being 1, j =
−3, ..., N . Then clearly
j Vn,θwj =
B2j,4 (Uθ,i) = ‖Bj,4‖
, ‖wj‖2 = 1,−3 ≤ j ≤ N
and in particular
0 Vn,θw0 ≤ λmax (Vn,θ) ‖w0‖2 = λmax (Vn,θ) ,
−3Vn,θw−3 ≥ λmin (Vn,θ) ‖w−3‖2 = λmin (Vn,θ) .
This, together with (A.5) yields that
κ = λmax (Vn,θ)λ
min (Vn,θ) ≥
wT0 Vn,θw0
wT−3Vn,θw−3
‖B0,4‖2n,θ
‖B−3,4‖2n,θ
‖B0,4‖2θ
‖B−3,4‖2θ
1 +An
which leads to κ ≥ C > 1, a.s. because the definition of B-spline and Assumption A2 ensure
that ‖B0,4‖2θ ≥ C0 ‖B−3,4‖
for some constant C0 > 1. Next applying Lemma A.4 with
κ2 − 1
)1/16 (
κ2 + 1
)−1/16
and c0 = ν
, one gets
≤ 2ν−8N (1− ν)−1 =
CN, a.s.. Hence part one of (A.12) follows. Part two of (A.12) is proved in the same fashion.
In the following, we denote by QT (m) the 4-th order quasi-interpolant of m corresponding
to the knots T , see equation (4.12), page 146 of DeVore and Lorentz (1993). According to
Theorem 7.7.4, DeVore and Lorentz (1993), the following lemma holds.
Lemma A.6. There exists a constant C > 0, such that for 0 ≤ k ≤ 2 and γ ∈ C(4) [0, 1]
(γ −QT (γ))(k)
h4−k,
Lemma A.7. Under Assumptions A2, A3, A5 and A6, there exists an absolute constant C > 0,
such that for function γ̃θ (u) in (A.7)
θ∈Sd−1c
(γ̃θ − γθ)
h4−k, a.s., 0 ≤ k ≤ 2, (A.16)
Proof. According to Theorem A.1 of Huang (2003), there exists an absolute constant C > 0,
such that
θ∈Sd−1c
‖γ̃θ − γθ‖∞ ≤ C sup
θ∈Sd−1c
γ∈Γ(2)
‖γ − γθ‖∞ ≤ C
h4, a.s., (A.17)
SINGLE-INDEX PREDICTION MODEL 21
which proves (A.16) for the case k = 0. Applying Lemma A.6, one has for 0 ≤ k ≤ 2
θ∈Sd−1c
{QT (γθ)− γθ}
≤ C sup
θ∈Sd−1c
h4−k ≤ C
h4−k, (A.18)
As a consequence of (A.17) and (A.18) for the case k = 0, one has
θ∈Sd−1c
‖QT (γθ)− γ̃θ‖∞ ≤ C
h4, a.s.,
which, according to the differentiation of B-spline given in de Boor (2001), entails that
θ∈Sd−1c
{QT (γθ)− γ̃θ}
h4−k, a.s., 0 ≤ k ≤ 2. (A.19)
Combining (A.18) and (A.19) proves (A.16) for k = 1, 2.
Lemma A.8. Under Assumptions A1, A2, A4 and A5, there exists an absolute constant C > 0,
such that
1≤p≤d
θ∈Sd−1c
{γ̃θ (Uθ,i)− γθ (Uθ,i)}ni=1
h3, a.s., (A.20)
1≤p,q≤d
θ∈Sd−1c
∂θp∂θq
{γ̃θ (Uθ,i)− γθ (Uθ,i)}ni=1
h2, a.s.. (A.21)
Proof. According to the definition of γ̃θ in (A.7), and the fact that QT (γθ) is a cubic spline
on the knots T
{{QT (γθ)− γ̃θ} (Uθ,i)}ni=1 = Pθ {{QT (γθ)− γθ} (Uθ,i)}
which entails that
{{QT (γθ)− γ̃θ} (Uθ,i)}ni=1 =
Pθ {{QT (γθ)− γθ} (Uθ,i)}ni=1
= Ṗp {{QT (γθ)− γθ} (Uθ,i)}ni=1 +Pθ
{{QT (γθ)− γθ} (Uθ,i)}ni=1 .
Since
{{QT (γθ)− γθ} (Uθ,i)}ni=1 =
(Uθ,i)
{QT (γθ)− γθ} (Uθ,i)Xip
applying (A.19) to the decomposition above produces (A.20). The proof of (A.21) is similar.
22 LI WANG AND LIJIAN YANG
Lemma A.9. Under Assumptions A2, A5 and A6, there exists a constant C > 0 such that
θ∈Sd−1c
∥n−1BTθ
∞ ≤ Ch, a.s., sup
1≤p≤d
θ∈Sd−1c
n−1ḂTp
≤ C, a.s., (A.22)
θ∈Sd−1c
‖Pθ‖∞ ≤ C, a.s., sup
1≤p≤d
θ∈Sd−1c
≤ Ch−1, a.s.. (A.23)
Proof. To prove (A.22), observe that for any vector a ∈ Rn, with probability 1
∥n−1BTθ a
∞ ≤ ‖a‖∞ max−3≤j≤N
Bj,4 (Uθ,i)
≤ Ch ‖a‖∞ ,
n−1ḂTp a
≤ ‖a‖∞ max−3≤j≤N
{(Bj,3 −Bj+1,3) (Uθ,i)} Ḟd (Xθ,i)Xi,p
≤ C ‖a‖∞ .
To prove (A.23), one only needs to use (A.12), (A.22) and (3.1).
Lemma A.10. Under Assumptions A2 and A4-A6, one has with probability 1
θ∈Sd−1c
BTθ E
= max
−3≤j≤N
Bj,4 (Uθ,i)σ (Xi) εi
log n√
, (A.24)
1≤p≤d
θ∈Sd−1c
BTθ E
= sup
1≤p≤d
θ∈Sd−1c
ḂTpE
log n√
. (A.25)
Similarly, under Assumptions A2, A4-A6, with probability 1
θ∈Sd−1c
BTθ Eθ
= sup
θ∈Sd−1c
−3≤j≤N
Bj,4 (Uθ,i) {m (Xi)− γθ (Uθ,i)}
log n√
(A.26)
1≤p≤d
θ∈Sd−1c
BTθ Eθ
log n√
, a.s.. (A.27)
Proof. We decompose the noise variable εi into a truncated part and a tail part εi = ε
i,1 +
εDni,2 +m
i , where Dn = n
η (1/3 < η < 2/5), εDni,1 = εiI {|εi| > Dn},
εDni,2 = εiI {|εi| ≤ Dn} −m
i = E [εiI {|εi| ≤ Dn} |Xi] .
It is straightforward to verify that the mean of the truncated part is uniformly bounded by
D−2n , so the boundedness of B spline basis and of the function σ
2 entail that
θ∈Sd−1c
Bj,4 (Uθ,i) σ (Xi)m
n−2/3
SINGLE-INDEX PREDICTION MODEL 23
The tail part vanishes almost surely
P {|εn| > Dn} ≤
D−3n <∞.
Borel-Cantelli Lemma implies that
Bj,4 (Uθ,i) σ (Xi) ε
, for any k > 0.
For the truncated part, using Bernstein’s inequality and discretization as in Lemma A.2
θ∈Sd−1c
1≤j≤N
Bj,4 (Uθ,i) σ (Xi) ε
log n/
, a.s..
Therefore (A.24) is established as with probability 1
θ∈Sd−1c
n−2/3
log n/
log n/
The proofs of (A.25), (A.26) are similar as E {m (Xi)− γθ (Uθ,i) |Uθ,i } ≡ 0, but no truncation
is needed for (A.26) as sup
θ∈Sd−1c
1≤i≤n
|m (Xi)− γθ (Uθ,i)| ≤ C <∞. Meanwhile, to prove (A.27),
we note that for any p = 1, ..., d
[Bj,4 (Uθ,i) {m (Xi)− γθ (Uθ,i)}]
According to (2.6), one has γθ (Uθ) ≡ E {m (X) |Uθ}, hence
E [Bj,4 (Uθ) {m (X)− γθ (Uθ)}] ≡ 0,−3 ≤ j ≤ N, θ ∈ Sd−1c .
Applying Assumptions A2 and A3, one can differentiate through the expectation, thus
[Bj,4 (Uθ) {m (X)− γθ (Uθ)}]
≡ 0, 1 ≤ p ≤ d,−3 ≤ j ≤ N, θ ∈ Sd−1c ,
which allows one to apply the Bernstein’s inequality to obtain that with probability 1
[Bj,4 (Uθ,i) {m (Xi)− γθ (Uθ,i)}]
(nh)−1/2 log n
which is (A.27).
Lemma A.11. Under Assumptions A2 and A4-A6, for ε̂θ (u) in (A.9), one has
θ∈Sd−1c
u∈[0,1]
|ε̂θ (u)| = O
log n
, a.s.. (A.28)
24 LI WANG AND LIJIAN YANG
Proof. Denote â ≡ (â−3, · · · , âN )T =
BTθ Bθ
BTθ E = V
n−1BTθ E
, then ε̂θ (u) =
j=−3 âjBj,4 (u), so the order of ε̂θ (u) is related to that of â. In fact, by Theorem 5.4.2
in DeVore and Lorentz (1993)
θ∈Sd−1c
u∈[0,1]
|ε̂θ (u)| ≤ sup
θ∈Sd−1c
‖â‖∞ =
θ∈Sd−1c
n−1BTθ E
≤ CN sup
θ∈Sd−1c
∥n−1BTθ E
∞ , a.s.,
where the last inequality follows from (A.12) of Lemma A.5. Applying (A.24) of Lemma A.10,
we have established (A.28).
Lemma A.12. Under Assumptions A2 and A4-A6, for ε̃θ (u) in (A.8), one has
θ∈Sd−1c
u∈[0,1]
|ε̃θ (u)| = O
log n
, a.s.. (A.29)
The proof is similar to Lemma A.11, thus omitted.
The next result evaluates the uniform size of the noise derivatives.
Lemma A.13. Under Assumptions A2-A6, one has with probability 1
1≤p≤d
θ∈Sd−1c
1≤i≤n
ε̂θ (Uθ,i)
(nh3)−1/2 log n
, (A.30)
1≤p≤d
θ∈Sd−1c
1≤i≤n
ε̃θ (Uθ,i)
(nh3)−1/2 log n
, (A.31)
1≤p,q≤d
θ∈Sd−1c
1≤i≤n
∂θp∂θq
ε̂θ (Uθ,i)
(nh5)−1/2 log n
, (A.32)
1≤p,q≤d
θ∈Sd−1c
1≤i≤n
∂θp∂θq
ε̃θ (Uθ,i)
(nh5)−1/2 log n
. (A.33)
Proof. Note that
ε̂θ (Uθ,i)
= (I−Pθ) Ḃp
θ E+Bθ
p (I−Pθ)E.
Applying (A.24) and (A.25) of Lemma A.10, (A.12) of Lemma A.5, (A.22) and (A.23) of
Lemma A.9, one derives (A.30). To prove (A.31), note that
ε̃θ (Uθ,i)
{PθEθ} = ṖpEθ +Pθ
Eθ = T1 + T2, (A.34)
in which
(I−Pθ) Ḃp −Bθ
(I−Pθ) Ḃp −Bθ
BTθ Bθ
ḂTpBθ
BTθ Bθ
BTθ Eθ
SINGLE-INDEX PREDICTION MODEL 25
T2 = Bθ
BTθ Bθ
BTθ Eθ
By (A.24), (A.12), (A.22) and (A.23), one derives
θ∈Sd−1c
‖T1‖∞ = O
n−1/2N3/2 log n
, a.s., (A.35)
while (A.27) of Lemma A.10, (A.12) of Lemma A.5
θ∈Sd−1c
‖T2‖∞ = N ×O
n−1/2h−1/2 log n
n−1/2h−3/2 log n
, a.s.. (A.36)
Now, putting together (A.34), (A.35) and (A.36), we have established (A.31). The proof for
(A.32) and (A.33) are similar.
Proof of Proposition 2.1. According to the decomposition (A.6)
|γ̂θ (u)− γθ (u)| = |{γ̃θ (u)− γθ (u)}+ ε̃θ (u) + ε̂θ (u)| .
Then (2.13) follows directly from (A.16) of Lemma A.7, (A.28) of Lemma A.11 and (A.29) of
Lemma A.12. Again by definitions (A.8) and (A.9), we write
{(γ̂θ − γθ) (Uθ,i)} =
(γ̃θ − γθ) (Uθ,i) +
γ̃θ (Uθ,i) +
ε̂θ (Uθ,i) .
It is clear from (A.20), (A.30) and (A.31) that with probability 1
1≤p≤d
θ∈Sd−1c
1≤i≤n
(γ̃θ − γθ) (Uθ,i)
1≤p≤d
θ∈Sd−1c
1≤i≤n
ε̃θ (Uθ,i)
ε̂θ (Uθ,i)
)−1/2
log n
Putting together all the above yields (2.14). The proof of (2.15) is similar.
A.3. Proof of Proposition 2.2
Lemma A.14. Under Assumptions A2-A6, one has
θ∈Sd−1c
R̂ (θ)−R (θ)
= o(1), a.s..
Proof. For the empirical risk function R̂ (θ) in (2.9), one has
R̂ (θ) = n−1
{γ̂θ (Uθ,i)−m (Xi)− σ (Xi) εi}2
26 LI WANG AND LIJIAN YANG
= n−1
{γ̂θ (Uθ,i)− γθ (Uθ,i) + γθ (Uθ,i)−m (Xi)− σ (Xi) εi}2 ,
hence
R̂ (θ) = n−1
{γ̂θ (Uθ,i)− γθ (Uθ,i)}2 + n−1
σ2 (Xi) ε
+2n−1
{γ̂θ (Uθ,i)− γθ (Uθ,i)} {γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
{γθ (Uθ,i)−m (Xi)}2 + 2n−1
{γθ (Uθ,i)−m (Xi)}σ (Xi) εi,
where γ̂θ (x) is defined in (2.8). Using the expression of R (θ) in (2.7), one has
θ∈Sd−1c
∣R̂ (θ)−R (θ)
∣ ≤ I1 + I2 + I3 + I4,
I1 = sup
θ∈Sd−1c
{γ̂θ (Uθ,i)− γθ (Uθ,i)}2
I2 = sup
θ∈Sd−1c
{γ̂θ (Uθ,i)− γθ (Uθ,i)} {γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
I3 = sup
θ∈Sd−1c
{γθ (Uθ,i)−m (Xi)}2 − E {γθ (Uθ)−m (X)}2
I4 = sup
θ∈Sd−1c
σ2 (Xi) ε
i − Eσ2 (X)
{γθ (Uθ,i)−m (Xi)}σ (Xi) εi
Bernstein inequality and strong law of large number for α mixing sequence imply that
I3 + I4 = o(1), a.s.. (A.37)
Now (2.13) of Proposition 2.1 provides that
θ∈Sd−1c
u∈[0,1]
|γ̂θ (u)− γθ (u)| = O
n−1/2h−1/2 log n+ h4
, a.s.,
which entail that
I1 = O
n−1/2h−1/2 log n
, a.s., (A.38)
I2 ≤ O
(nh)−1/2 log n+ h4
× sup
θ∈Sd−1c
|γθ (Uθ,i)−m (Xi)− σ (Xi) εi| .
Hence
I2 ≤ O
n−1/2h−1/2 log n+ h4
, a.s.. (A.39)
The lemma now follows from (A.37), (A.38) and (A.39) and Assumption A6.
SINGLE-INDEX PREDICTION MODEL 27
Lemma A.15. Under Assumptions A2 - A6, one has
θ∈Sd−1c
1≤p≤d
R̂ (θ)−R (θ)
− n−1
ξθ,i,p
n−1/2
, a.s., (A.40)
in which
ξθ,i,p = 2 {γθ (Uθ,i)− Yi}
γθ (Uθ,i)−
R (θ) , E (ξθ,i,p) = 0. (A.41)
Furthermore for k = 1, 2
θ∈Sd−1c
R̂ (θ)−R (θ)
n−1/2h−1/2−k log n+ h4−k
, a.s.. (A.42)
Proof. Note that for any p = 1, 2, ..., d
R̂ (θ) = n−1
{γ̂θ (Uθ,i)− Yi}
γ̂θ (Uθ,i) ,
R (θ) = E
{γθ (Uθ)−m (X)}
γθ (Uθ)
{γθ (Uθ)−m (X)− σ (X) ε}
γθ (Uθ)
Thus E (ξθ,i,p) = 2E
{γθ (Uθ,i)− Yi} ∂∂θpγθ (Uθ,i)
R (θ) = 0 and
R̂ (θ)−R (θ)
= (2n)
ξθ,i,p + J1,θ,p + J2,θ,p + J3,θ,p, (A.43)
J1,θ,p = n
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
(γ̂θ − γθ) (Uθ,i) ,
J2,θ,p = n
{γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
(γ̂θ − γθ) (Uθ,i) ,
J3,θ,p = n
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
γθ (Uθ,i) .
Bernstein inequality implies that
θ∈Sd−1c
1≤p≤d
ξθ,i,p
n−1/2 log n
, a.s.. (A.44)
28 LI WANG AND LIJIAN YANG
Meanwhile, applying (2.13) and (2.14) of Proposition 2.1, one obtains that
θ∈Sd−1c
1≤p≤d
|J1,θ,p| = O
log n+ h4
)−1/2
log n+ h3
n−1h−2 log2 n+ h7
, a.s.. (A.45)
Note that
J2,θ,p = n
{γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
(γ̃θ − γθ) (Uθ,i)
−n−1 (E+Eθ)T
{Pθ (E+Eθ)} .
Applying (2.13), one gets
θ∈Sd−1c
1≤p≤d
J2,θ,p + n
−1 (E+Eθ)
{Pθ (E+Eθ)}
, a.s.,
while (A.24), (A.26) and (A.12) entail that with probability 1
θ∈Sd−1c
1≤p≤d
n−1 (E+Eθ)
{Pθ (E+Eθ)}
log n
×N ×N ×O
log n
n−1N log2 n
θ∈Sd−1c
1≤p≤d
|J2,θ,p| = O
h3 + n−1N log2 n
, a.s.. (A.46)
Lastly
J3,θ,p − n−1
(γ̃θ − γθ)
γθ (Uθ,i) = n
−1 (E+Eθ)
BTθ Bθ
By applying (A.24), (A.26), and (A.12), it is clear that with probability 1
θ∈Sd−1c
1≤p≤d
n−1BTθ E+n
BTθ Bθ
log n
×N ×O
h+ (nN)
log n
n−1 log2 n+ (nN)−1/2 log n
while by applying (A.16) of Lemma A.7, one has
θ∈Sd−1c
1≤p≤d
(γ̃θ − γθ)
γθ (Uθ,i)
, a.s.,
SINGLE-INDEX PREDICTION MODEL 29
together, the above entail that
θ∈Sd−1c
1≤p≤d
|J3,θ,p| = O
h4 + n−1 log2 n+ (nN)−1/2 log n
, a.s.. (A.47)
Therefore, (A.43), (A.45), (A.46), (A.47) and Assumption A6 lead to (A.40), which, together
with (A.44), establish (A.42) for k = 1.
Note that the second order derivative of R̂ (θ) and R (θ) with respect to θp, θq are
{γ̂θ (Uθ,i)− Yi}
∂θp∂θq
γ̂θ (Uθ,i) +
γ̂θ (Uθ,i)
γ̂θ (Uθ,i)
E {γθ (Uθ)−m (X)}
∂θp∂θq
γθ (Uθ) + E
γθ (Uθ)
γθ (Uθ)
The proof of (A.42) for k = 2 follows from (2.13), (2.14) and (2.15).
Proof of Proposition 2.2. The result follows from Lemma A.14, Lemma A.15, equations
(A.50) and (A.51).
A.4. Proof of the Theorem 2
Let Ŝ∗p (θ−d) be the p-th element of Ŝ
∗ (θ−d) and for γθ in (2.6), denote
ηi,p := 2
γ̇p − θ0,pθ−10,dγ̇d
(Uθ0,i) {γθ0 (Uθ0,i)− Yi} , (A.48)
where γ̇p is value of
γθ taking at θ = θ0, for any p, q = 1, 2, ..., d − 1.
Lemma A.16. Under Assumptions A2-A6, one has
1≤p≤d−1
Ŝ∗p (θ0,−d)− n−1
n−1/2
, a.s.. (A.49)
Proof. For any p = 1, ..., d − 1
Ŝ∗p (θ−d)− S∗p (θ−d) =
− θpθ−1d
R̂ (θ)−R (θ)
Therefore, according to (A.40), (A.41) and (A.48)
ηi,p = n
ξθ0,i,p − θ0,pθ
ξθ0,i,d, E (ηi,p) = 0,
1≤p≤d−1
Ŝ∗p (θ0,−d)− S∗p (θ0,−d)− n−1
n−1/2
, a.s..
Since S∗ (θ−d) attains its minimum at θ0,−d, for p = 1, ..., d − 1
S∗p (θ0,−d) ≡
− θpθ−1d
R (θ)
which yields (A.49).
30 LI WANG AND LIJIAN YANG
Lemma A.17. The (p, q)-th entry of the Hessian matrix H∗ (θ0,−d) equals lp,q given in Theo-
rem 2.
Proof. It is easy to show that for any p, q = 1, 2, ..., d,
R (θ) =
E {m (X)− γθ (Uθ)}2 = −2E
γθ (Uθ)
γθ (Uθ)
∂θp∂θq
R (θ) = −2E
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θq
γθ (Uθ)
Note that
R∗ (θ−d) =
R (θ)−
R (θ) , (A.50)
∂θp∂θq
R∗ (θ−d) =
∂θp∂θq
R (θ)−
∂θp∂θd
R (θ)−
∂θd∂θq
R (θ)
1− ‖θ−d‖22
R (θ) +
∂θd∂θd
R (θ) . (A.51)
R∗ (θ−d) = −2E
γθ (Uθ)
γθ (Uθ)
+ 2θ−1
γθ (Uθ)
γθ (Uθ)
∂θp∂θq
R∗ (θ−d) = −2E
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θq
γθ (Uθ)
+2θqθ
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θd
γθ (Uθ)
1− ‖θ−d‖22
γθ (Uθ)
γθ (Uθ)
+2θpθ
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θq
γθ (Uθ)
−2θpθqθ−2d E
γθ (Uθ)
+ γθ (Uθ)
∂θd∂θd
γθ (Uθ)
Therefore we obtained the desired result.
Proof of Theorem 2. For any p = 1, 2, ..., d − 1, let
fp (t) = Ŝ
tθ̂−d + (1− t) θ0,−d
, t ∈ [0, 1],
fp (t) =
tθ̂−d+(1− t) θ0,−d
θ̂q − θ0,q
SINGLE-INDEX PREDICTION MODEL 31
Note that Ŝ∗ (θ−d) attains its minimum at θ̂−d, i.e., Ŝ
≡ 0. Thus, for any p = 1, 2, ..., d−
1, tp ∈ [0, 1], one has
−Ŝ∗p (θ0,−d) = Ŝ∗p
− Ŝ∗p (θ0,−d) = fp (1)− fp (0)
∂θqθp
tpθ̂−d + (1− tp) θ0,−d
q=1,...,d−1
θ̂−d−θ0,−d
−Ŝ∗ (θ0,−d) =
∂θq∂θp
tpθ̂−d + (1− tp) θ0,−d
p,q=1,...,d−1
θ̂−d − θ0,−d
Now (2.11) of Theorem 1 and Proposition 2.2 with k = 2 imply that uniformly in p, q =
1, 2, ..., d − 1
∂θq∂θp
tpθ̂−d+(1− tp) θ0,−d
−→ lq,p, a.s., (A.52)
where lp,q is given in Theorem 2. Noting that
θ̂−d−θ0,−d
is represented as
∂θq∂θp
tpθ̂−d + (1− tp) θ0,−d
p,q=1,...,d−1
nŜ∗ (θ0,−d) ,
where Ŝ∗ (θ0,−d) =
Ŝ∗p (θ0,−d)
and according to (A.48) and Lemma A.16
Ŝ∗p (θ0,−d) = n
ηp,i + o
n−1/2
, a.s., E (ηp,i) = 0.
Let Ψ (θ0) = (ψpq)
p,q=1
be the covariance matrix of
Ŝ∗p (θ0,−d)
with ψpq given in
Theorem 2. Cramér-Wold device and central limit theorem for α mixing sequences entail that
nŜ∗ (θ0,−d)
d−→ N {0,Ψ(θ0)} .
Let Σ (θ0) = {H∗ (θ0,−d)}−1Ψ(θ0)
{H∗ (θ0,−d)}T
, withH∗ (θ0,−d) being the Hessian matrix
defined in (2.3). The above limiting distribution of
nŜ∗ (θ0,−d), (A.52) and Slutsky’s theorem
imply that
θ̂−d−θ0,−d
d−→ N {0,Σ (θ0)} .
References
Bosq, D. (1998). Nonparametric Statistics for Stochastic Processes. Springer-Verlag, New
York.
32 LI WANG AND LIJIAN YANG
Carroll, R., Fan, J., Gijbles, I. and Wand, M. P. (1997). Generalized partially linear single-
index models. J. Amer. Statist. Assoc. 92 477-489.
Chen, H. (1991). Estimation of a projection -persuit type regression model. Ann. Statist. 19
142-157.
de Boor, C. (2001). A Practical Guide to Splines. Springer-Verlag, New York.
DeVore, R. A. and Lorentz, G. G. (1993). Constructive Approximation: Polynomials and
Splines Approximation. Springer-Verlag, Berlin.
Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Applications. Chapman
and Hall, London.
Friedman, J. H. and Stuetzle, W. (1981). Projection pursuit regression. J. Amer. Statist.
Assoc. 76 817-823.
Härdle, W. (1990). Applied Nonparametric Regression. Cambridge University Press, Cam-
bridge.
Härdle, W. and Hall, P. and Ichimura, H. (1993). Optimal smoothing in single-index models.
Ann. Statist. 21 157-178.
Härdle, W. and Stoker, T. M. (1989). Investigating smooth multiple regression by the method
of average derivatives. J. Amer. Statist. Assoc. 84 986-995.
Hall, P. (1989). On projection pursuit regression. Ann. Statist. 17 573-588.
Hastie, T. J. and Tibshirani, R. J. (1990). Generalized Additive Models. Chapman and Hall,
London.
Horowitz, J. L. and Härdle, W. (1996). Direct semiparametric estimation of single-index
models with discrete covariates. J. Amer. Statist. Assoc. 91 1632-1640.
Hristache, M., Juditski, A. and Spokoiny, V. (2001). Direct estimation of the index coefficients
in a single-index model. Ann. Statist. 29 595-623.
Huang, J. Z. (2003). Local asymptotics for polynomial spline regression. Ann. Statist. 31
1600-1635.
Huang, J. and Yang, L. (2004). Identification of nonlinear additive autoregressive models. J.
R. Stat. Soc. Ser. B Stat. Methodol. 66 463-477.
SINGLE-INDEX PREDICTION MODEL 33
Huber, P. J. (1985). Projection pursuit (with discussion). Ann. Statist. 13 435-525.
Ichimura, H. (1993). Semiparametric least squares (SLS) and weighted SLS estimation of
single-index models Journal of Econometrics 58 71-120.
Klein, R. W. and Spady. R. H. (1993). An efficient semiparametric estimator for binary
response models. Econometrica 61 387-421.
Mammen, E., Linton, O. and Nielsen, J. (1999). The existence and asymptotic properties of
a backfitting projection algorithm under weak conditions. Ann. Statist. 27 1443-1490.
Pham, D. T. (1986). The mixing properties of bilinear and generalized random coefficient
autoregressive models. Stochastic Anal. Appl. 23 291-300.
Powell, J. L., Stock, J. H. and Stoker, T. M. (1989). Semiparametric estimation of index
coefficients. Econometrica. 57 1403-1430.
Tong, H. (1990) Nonlinear Time Series: A Dynamical System Approach. Oxford, U.K.:
Oxford University Press.
Tong, H., Thanoon, B. and Gudmundsson, G. (1985) Threshold time series modeling of two
icelandic riverflow systems. Time Series Analysis in Water Resources. ed. K. W. Hipel,
American Water Research Association.
Wang, L. and Yang, L. (2007). Spline-backfitted kernel smoothing of nonlinear additive
autoregression model. Ann. Statist. Forthcoming.
Xia, Y. and Li, W. K. (1999). On single-index coefficient regression models. J. Amer. Statist.
Assoc. 94 1275-1285.
Xia, Y., Li, W. K., Tong, H. and Zhang, D. (2004). A goodness-of-fit test for single-index
models. Statist. Sinica. 14 1-39.
Xia, Y., Tong, H., Li, W. K. and Zhu, L. (2002). An adaptive estimation of dimension
reduction space. J. R. Stat. Soc. Ser. B Stat. Methodol. 64 363-410.
Xue, L. and Yang, L. (2006 a). Estimation of semiparametric additive coefficient model. J.
Statist. Plann. Inference 136, 2506-2534.
Xue, L. and Yang, L. (2006 b). Additive coefficient modeling via polynomial spline. Statistica
Sinica 16 1423-1446.
Table 2: Report of Example 2
Sample Size n Dimension d
Average MSE Time
MAVE SIP MAVE SIP
4 0.00020 0.00018 1.91 0.19
10 0.00031 0.00043 2.17 0.10
30 0.00106 0.00285 2.77 0.13
50 0.00031 0.00043 3.29 0.10
100 0.00681 0.00620 5.94 0.31
200 0.00529 0.00407 27.90 0.49
4 0.00008 0.00008 3.28 0.09
10 0.00012 0.00017 3.93 0.13
30 0.00017 0.00058 5.41 0.15
50 0.00032 0.00127 8.48 0.16
100 — 0.00395 — 0.44
200 — 0.00324 — 0.73
4 0.00004 0.00003 5.32 0.17
10 0.00005 0.00007 7.49 0.24
30 0.00006 0.00017 10.08 0.26
50 0.00007 0.00030 15.42 0.24
100 0.00015 0.00061 40.81 0.54
200 — 0.00197 — 1.44
4 0.00002 0.00001 14.44 0.76
10 0.00002 0.00003 24.54 0.79
30 0.00002 0.00008 32.51 0.83
50 0.00002 0.00010 52.93 0.89
100 0.00003 0.00012 143.07 0.99
200 0.00004 0.00020 386.80 1.96
400 — 0.00054 — 4.98
4 0.00001 0.00001 33.57 1.95
10 0.00001 0.00001 62.54 3.64
30 0.00001 0.00002 92.41 1.95
50 0.00001 0.00003 155.38 2.72
100 0.00001 0.00005 275.73 1.81
200 0.00008 0.00006 2432.56 2.84
400 — 0.00010 — 9.35
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
−1 0 1 2
−1 0 1 2
Figure 1: Example 1. (a) and (b) Plots of the actual surface m in model (4.1) with re-
spect to δ = 0, 1; (c) and (d) Plots of various univariate functions with respect to δ = 0, 1:
XTi θ̂, Yi
, 1 ≤ i ≤ 50 (dots); the univariate function g (solid line); the estimated function of
g by plugging in the true index coefficient θ0 (dotted line); the estimated function of g by
plugging in the estimated index coefficient (dashed line) θ̂ = (0.69016, 0.72365)T for δ = 0 and
(0.72186, 0.69204)T for δ = 1.
−2 −1 0 1 2
n= 100 ,  d= 10
−2 −1 0 1 2
n= 100 ,  d= 50
−2 −1 0 1 2
n= 100 ,  d= 100
−3 −2 −1 0 1 2 3
n= 100 ,  d= 200
Figure 2: Example 2. Plots of the spline estimator of g with the estimated index parameter θ̂
(dotted curve), cubic spline estimator of g with the true index parameter θ0 (dashed curves),
the true function m (x) in (4.2) (solid curve), and the data scatter plots (dots).
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=10
n=100
n=200
n=500
n=1000
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=50
n=100
n=200
n=500
n=1000
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=100
n=100
n=200
n=500
n=1000
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=200
n=100
n=200
n=500
n=1000
Figure 3: Example 2. Kernel density estimators of the 100 ‖θ̂ − θ0‖/
0 200 400 600 800 1000
0 200 400 600 800 1000
0 200 400 600 800 1000
Figure 4: Time plots of the daily Jökulsá Eystri River data (a) river flow Yt (solid line) with its
trend (dashed line) (b) temperature Xt (solid line) with its trend (dashed line) (c) precipitation
Zt (solid line) with its trend (dashed line).
+++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++
++++++
+++++
+++++
++++++++++
++++++++++
+++++++
+++++++
++++++++++++
+++++
++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++
++++++++++++++++++
+++++
+++++
++++++
+++++++
+++++
+++++
+++++
+++++
++++++++
+++++
++++++++++++++
+++++++
+++++++++++++++++++++++++
++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++
++++++++++
+++++++++++++++++++++
++++++
++++++
++++++
++++++++++
+++++++++
+++++
++++++++++++++++++++++++++++++
+++++++
++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0 200 400 600 800 1000
0 200 400 600 800 1000
800 900 1000 1100
++++++++++++++++++++++++++++++++++++++++++++++++++
++++++
++++++++++
+++++++++++++++++++++
++++++
++++++
++++++
++++++++++
+++++++++
+++++
++++++++++++++++++++++++++++++
+++++++
++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Figure 5: (a) The scatter plot of the river flow (“+”) and the fitted plot of the river flow (line)
and (b) Residuals of the fitted SIP model (c) Out-of-sample rolling forecasts (line) of the river
flow for the entire third year (“+”) based on the first two years’ river flow.

  For the past two decades, single-index model, a special case of projection
pursuit regression, has proven to be an efficient way of coping with the high
dimensional problem in nonparametric regression. In this paper, based on weakly
dependent sample, we investigate the single-index prediction (SIP) model which
is robust against deviation from the single-index model. The single-index is
identified by the best approximation to the multivariate prediction function of
the response variable, regardless of whether the prediction function is a
genuine single-index function. A polynomial spline estimator is proposed for
the single-index prediction coefficients, and is shown to be root-n consistent
and asymptotically normal. An iterative optimization routine is used which is
sufficiently fast for the user to analyze large data of high dimension within
seconds. Simulation experiments have provided strong evidence that corroborates
with the asymptotic theory. Application of the proposed procedure to the rive
flow data of Iceland has yielded superior out-of-sample rolling forecasts.

Introduction
XTi , Yi
= {Xi,1, ...,Xi,d, Yi}ni=1 be a length n realization of a (d+ 1)-dimensional
strictly stationary process following the heteroscedastic model
Yi = m (Xi) + σ (Xi) εi,m (Xi) = E (Yi|Xi) , (1.1)
in which E (εi |Xi ) = 0, E
ε2i |Xi
= 1, 1 ≤ i ≤ n. The d-variate functions m, σ are the
unknown mean and standard deviation of the response Yi conditional on the predictor vector
Xi, often estimated nonparametrically. In what follows, we let
XT , Y, ε
have the stationary
distribution of
XTi , Yi, εi
. When the dimension of X is high, one unavoidable issue is the
“curse of dimensionality”, which refers to the poor convergence rate of nonparametric esti-
mation of general multivariate function. Much effort has been devoted to the circumventing
of this difficulty. In the words of Xia, Tong, Li and Zhu (2002), there are essentially two
approaches: function approximation and dimension reduction. A favorite function approxima-
tion technique is the generalized additive model advocated by Hastie and Tibshirani (1990),
Address for correspondence: Lijian Yang, Department of Statistics and Probability, Michigan State Univer-
sity, East Lansing, MI 48824, USA. E-mail: yang@stt.msu.edu
http://arxiv.org/abs/0704.0302v2
2 LI WANG AND LIJIAN YANG
see also, for example, Mammen, Linton and Nielsen (1999), Huang and Yang (2004), Xue
and Yang (2006 a, b), Wang and Yang (2007). An attractive dimension reduction method is
the single-index model, similar to the first step of projection pursuit regression, see Friedman
and Stuetzle (1981), Hall (1989), Huber (1985), Chen (1991). The basic appeal of single-index
model is its simplicity: the d-variate function m (x) = m (x1, ..., xd) is expressed as a univariate
function of xT θ0 =
p=1 xpθ0,p. Over the last two decades, many authors had devised various
intelligent estimators of the single-index coefficient vector θ0 = (θ0,1, ..., θ0,d)
, for instance,
Powell, Stock and Stoker (1989), Härdle and Stoker (1989), Ichimura (1993), Klein and Spady
(1993), Härdle, Hall and Ichimura (1993), Horowitz and Härdle (1996), Carroll, Fan, Gijbels
and Wand (1997), Xia and Li (1999), Hristache, Juditski and Spokoiny (2001). More recently,
Xia, Tong, Li and Zhu (2002) proposed the minimum average variance estimation (MAVE) for
several index vectors.
All the aforementioned methods assume that the d-variate regression function m (x) is
exactly a univariate function of some xT θ0 and obtain a root-n consistent estimator of θ0. If
this model is misspecified (m is not a genuine single-index function), however, a goodness-of-fit
test then becomes necessary and the estimation of θ0 must be redefined, see Xia, Li, Tong
and Zhang (2004). In this paper, instead of presuming that underlying true function m is
a single-index function, we estimate a univariate function g that optimally approximates the
multivariate function m in the sense of
g (ν) = E
m (X)|XT θ0 = ν
, (1.2)
where the unknown parameter θ0 is called the SIP coefficient, used for simple interpretation
once estimated; XT θ0 is the latent SIP variable; and g is a smooth but unknown function used
for further data summary, called the link prediction function. Our method therefore is clearly
interpretable regardless of the goodness-of-fit of the single-index model, making it much more
relevant in applications.
We propose estimators of θ0 and g based on weakly dependent sample, which includes
many existing nonparametric time series models, that are (i) computationally expedient and
(ii) theoretically reliable. Estimation of both θ0 and g has been done via the kernel smoothing
techniques in existing literature, while we use polynomial spline smoothing. The greatest
advantages of spline smoothing, as pointed out in Huang and Yang (2004), Xue and Yang
(2006 b) are its simplicity and fast computation. Our proposed procedure involves two stages:
estimation of θ0 by some
n-consistent θ̂, minimizing an empirical version of the mean squared
error, R(θ) = E{Y − E(Y |XT θ)}2; spline smoothing of Y on XT θ̂ to obtain a cubic spline
estimator ĝ of g. The best single-index approximation to m(x) is then m̂(x) = ĝ
xT θ̂
Under geometrically strong mixing condition, strong consistency and
n-rate asymptotic
SINGLE-INDEX PREDICTION MODEL 3
normality of the estimator θ̂ of the SIP coefficient θ0 in (1.2) are obtained. Proposition 2.2 is
the key in understanding the efficiency of the proposed estimator. It shows that the derivatives
of the risk function up to order 2 are uniformly almost surely approximated by their empirical
versions.
Practical performance of the SIP estimators is examined via Monte Carlo examples. The
estimator of the SIP coefficient performs very well for data of both moderate and high dimension
d, of sample size n from small to large, see Tables 1 and 2, Figures 1 and 2. By taking advantages
of the spline smoothing and the iterative optimization routines, one reduces the computation
burden immensely for massive data sets. Table 2 reports the computing time of one simulation
example on an ordinary PC, which shows that for massive data sets, the SIP method is much
faster than the MAVE method. For instance, the SIP estimation of a 200-dimensional θ0 from
a data of size 1000 takes on average mere 2.84 seconds, while the MAVE method needs to spend
2432.56 seconds on average to obtain a comparable estimates. Hence on account of criteria (i)
and (ii), our method is indeed appealing. Applying the proposed SIP procedure to the rive flow
data of Iceland, we have obtained superior forecasts, based on a 9-dimensional index selected
by BIC, see Figure 5.
The rest of the paper is organized as follows. Section 2 gives details of the model spec-
ification, proposed methods of estimation and main results. Section 3 describes the actual
procedure to implement the estimation method. Section 4 reports our findings in an extensive
simulation study. The proposed SIP model and the estimation procedure are applied in Section
5 to the rive flow data of Iceland. Most of the technical proofs are contained in the Appendix.
2. The Method and Main Results
2.1. Identifiability and definition of the index coefficient
It is obvious that without constraints, the SIP coefficient vector θ0 = (θ0,1, ..., θ0,d)
identified only up to a constant factor. Typically, one requires that ‖θ0‖ = 1 which entails
that at least one of the coordinates θ0,1, ..., θ0,d is nonzero. One could assume without loss of
generality that θ0,d > 0, and the candidate θ0 would then belong to the upper unit hemisphere
Sd−1+ =
(θ1, ..., θd) |
p=1 θ
p = 1, θd > 0
For a fixed θ = (θ1, ..., θd)
, denote Xθ = X
T θ, Xθ,i = X
i θ, 1 ≤ i ≤ n. Let
mθ (Xθ) = E (Y |Xθ) = E {m (X) |Xθ} . (2.1)
Define the risk function of θ as
R (θ) = E
{Y −mθ (Xθ)}2
= E {m (X)−mθ (Xθ)}2 + Eσ2 (X) , (2.2)
4 LI WANG AND LIJIAN YANG
which is uniquely minimized at θ0 ∈ Sd−1+ , i.e.
θ0 = arg min
θ∈Sd−1+
R (θ) .
Remark 2.1. Note that Sd−1+ is not a compact set, so we introduce a cap shape subset of
Sd−1+
Sd−1c =
(θ1, ..., θd) |
θ2p = 1, θd ≥
1− c2
, c ∈ (0, 1)
Clearly, for an appropriate choice of c, θ0 ∈ Sd−1c , which we assume in the rest of the paper.
Denote θ−d = (θ1, ..., θd−1)
, since for fixed θ ∈ Sd−1+ , the risk function R (θ) depends only
on the first d− 1 values in θ, so R (θ) is a function of θ−d
R∗ (θ−d) = R
θ1, θ2, ..., θd−1,
1− ‖θ−d‖22
with well-defined score and Hessian matrices
S∗ (θ−d) =
R∗ (θ−d) , H
∗ (θ−d) =
∂θ−d∂θ
R∗ (θ−d) . (2.3)
Assumption A1: The Hessian matrix H∗ (θ0,−d) is positive definite and the risk function R
is locally convex at θ0,−d, i.e., for any ε > 0, there exists δ > 0 such that R
∗ (θ−d)−R∗ (θ0,−d) <
δ implies ‖θ−d − θ0,−d‖2 < ε.
2.2. Variable transformation
Throughout this paper, we denote by Bda =
x ∈ Rd |‖x‖ ≤ a
the d-dimensional ball with
radius a and center 0 and
∣the kth order partial derivatives of m are continuous on B
the space of k-th order smooth functions.
Assumption A2: The density function of X, f (x) ∈ C(4)
, and there are constants
0 < cf ≤ Cf such that
cf/Vold
≤ f (x) ≤ Cf/Vold
, x ∈ Bda
f (x) ≡ 0, x /∈ Bda
For a fixed θ, define the transformed variables of the SIP variable Xθ
Uθ = Fd (Xθ) , Uθ,i = Fd (Xθ,i) , 1 ≤ i ≤ n, (2.4)
SINGLE-INDEX PREDICTION MODEL 5
in which Fd is the a rescaled centered Beta {(d+ 1) /2, (d+ 1) /2} cumulative distribution
function, i.e.
Fd (ν) =
∫ ν/a
Γ (d+ 1)
Γ {(d+ 1) /2}2 2d
1− t2
)(d−1)/2
dt, ν ∈ [−a, a] . (2.5)
Remark 2.2. For any fixed θ, the transformed variable Uθ in (2.4) has a quasi-uniform [0, 1]
distribution. Let fθ (u) be the probability density function of Uθ, then for any u ∈ [0, 1]
fθ (u) =
d (v)
fXθ (v) , v = F
d (u) ,
in which fXθ (v) = lim△ν→0 P (ν ≤ Xθ ≤ ν +△ν). Noting that xθ is exactly the projection of
x on θ, let Dν = {x|ν ≤ xθ ≤ ν +△ν} ∩Bda, then one has
P (ν ≤ Xθ ≤ ν +△ν) = P (X ∈ Dν) =
f (x) dx.
According to Assumption A2
cfVold(Dν)
Vold (B
≤ P (ν ≤ Xθ ≤ ν +△ν) ≤
CfVold(Dν)
Vold (B
On the other hand
Vold(Dν) = Vold−1(Jν)△ν + o (△ν) ,
where Jν = {x|xθ = v} ∩Bda. Note that the volume of Bda is πd/2ad/Γ (d/2 + 1) and
Vold−1 (Jν) = π(d−1)/2
a2 − ν2
)(d−1)/2
Γ {(d+ 1)/2} ,
Vold−1(Jν)
Vold (B
Γ (d+ 1)
}(d−1)/2
Therefore 0 < cf ≤ fθ (u) ≤ Cf <∞, for any fixed θ and u ∈ [0, 1].
In terms of the transformed SIP variable Uθ in (2.4), we can rewrite the regression function
mθ in (2.1) for fixed θ
γθ (Uθ) = E {m (X) |Uθ} = E {m (X) |Xθ} = mθ (Xθ) , (2.6)
then the risk function R (θ) in (2.2) can be expressed as
R (θ) = E
{Y − γθ (Uθ)}2
= E {m (X)− γθ (Uθ)}2 + Eσ2 (X) . (2.7)
2.3. Estimation Method
6 LI WANG AND LIJIAN YANG
Estimation of both θ0 and g requires a degree of statistical smoothing, and all estimation
here is carried out via cubic spline. In the following, we define the estimator θ̂ of θ0 and the
estimator ĝ of g.
To introduce the space of splines, we pre-select an integer n1/6 ≪ N = Nn ≪ n1/5 (log n)−2/5,
see Assumption A6 below. Divide [0, 1] into (N + 1) subintervals Jj = [tj, tj+1), j = 0, ..., N −
1, JN = [tN , 1], where T := {tj}Nj=1 is a sequence of equally-spaced points, called interior knots,
given as
t1−k = ... = t−1 = t0 = 0 < t1 < ... < tN < 1 = tN+1 = ... = tN+k,
in which tj = jh, j = 0, 1, ..., N +1, h = 1/ (N + 1) is the distance between neighboring knots.
The j-th B-spline of order k for the knot sequence T denoted by Bj,k is recursively defined by
de Boor (2001).
Denote by Γ(k−2) = Γ(k−2) [0, 1] the space of all C(k−2) [0, 1] functions that are polynomials
of degree k−1 on each interval. For fixed θ, the cubic spline estimator γ̂θ of γθ and the related
estimator m̂θ of mθ are defined as
γ̂θ (·) = arg min
γ(·)∈Γ(2)[0,1]
{Yi − γ (Uθ,i)}2 , m̂θ (ν) = γ̂θ {Fd (ν)} . (2.8)
Define the empirical risk function of θ
R̂ (θ) = n−1
{Yi − γ̂θ (Uθ,i)}2 = n−1
{Yi − m̂θ (Xθ,i)}2 , (2.9)
then the spline estimator of the SIP coefficient θ0 is defined as
θ̂ = arg min
θ∈Sd−1c
R̂ (θ) ,
and the cubic spline estimator of g is m̂θ with θ replaced by θ̂, i.e.
ĝ (ν) =
arg min
γ(·)∈Γ(2)[0,1]
Yi − γ
{Fd (ν)} . (2.10)
2.4. Asymptotic results
Before giving the main theorems, we state some other assumptions.
Assumption A3: The regression function m ∈ C(4)
for some a > 0.
Assumption A4: The noise ε satisfies E (ε |X) = 0, E
ε2 |X
= 1 and there exists a positive
constant M such that sup
|ε|3 |X = x
< M . The standard deviation function σ (x) is
continuous on Bda,
0 < cσ ≤ inf
x∈Bda
σ (x) ≤ sup
x∈Bda
σ (x) ≤ Cσ <∞.
SINGLE-INDEX PREDICTION MODEL 7
Assumption A5: There exist positive constants K0 and λ0 such that α (n) ≤ K0e−λ0n holds
for all n, with the α-mixing coefficient for
XTi , εi
defined as
α (k) = sup
B∈σ{Zs,s≤t},C∈σ{Zs,s≥t+k}
|P (B ∩ C)− P (B)P (C)| , k ≥ 1.
Assumption A6: The number of interior knots N satisfies: n1/6 ≪ N ≪ n1/5 (log n)−2/5.
Remark 2.3. Assumptions A3 and A4 are typical in the nonparametric smoothing literature,
see for instance, Härdle (1990), Fan and Gijbels (1996), Xia, Tong Li and Zhu (2002). By
the result of Pham (1986), a geometrically ergodic time series is a strongly mixing sequence.
Therefore, Assumption A5 is suitable for (1.1) as a time series model under aforementioned
assumptions.
We now state our main results in the next two theorems.
Theorem 1. Under Assumptions A1-A6, one has
θ̂−d−→ θ0,−d, a.s.. (2.11)
Proof. Denote by (Ω,F ,P) the probability space on which all
XTi , Yi
are defined. By
Proposition 2.2, given at the end of this section
‖θ−d‖2≤
R̂∗ (θ−d)−R∗ (θ−d)
−→ 0, a.s.. (2.12)
So for any δ > 0 and ω ∈ Ω, there exists an integer n0 (ω), such that when n > n0 (ω),
R̂∗ (θ0,−d, ω) − R∗ (θ0,−d) < δ/2. Note that θ̂−d = θ̂−d (ω) is the minimizer of R̂∗ (θ−d, ω),
so R̂∗
θ̂−d (ω) , ω
− R∗ (θ0,−d) < δ/2. Using (2.12), there exists n1 (ω), such that when
n > n1 (ω), R
θ̂−d (ω) , ω
− R̂∗
θ̂−d (ω) , ω
< δ/2. Thus, when n > max (n0 (ω) , n1 (ω)),
θ̂−d (ω) , ω
−R∗ (θ0,−d) < δ/2 + R̂∗
θ̂−d (ω) , ω
−R∗ (θ0,−d) < δ/2 + δ/2 = δ.
According to Assumption A1, R∗ is locally convex at θ0,−d, so for any ε > 0 and any ω, if
θ̂−d (ω) , ω
−R∗ (θ0,−d) < δ, then
∥θ̂−d (ω)−θ0,−d
∥ < ε for n large enough , which implies
the strong consistency.
Theorem 2. Under Assumptions A1-A6, one has
θ̂−d−θ0,−d
d−→ N {0,Σ (θ0)} ,
where Σ (θ0) = {H∗ (θ0,−d)}−1Ψ(θ0) {H∗ (θ0,−d)}−1, H∗ (θ0,−d) = {lpq}d−1p,q=1 and Ψ(θ0) =
{ψpq}d−1p,q=1 with
lp,q = −2E [{γ̇pγ̇q + γθ0 γ̈p,q} (Uθ0)] + 2θ0,qθ
E [{γ̇pγ̇d (Uθ0) + γθ0 γ̈p,d} (Uθ0)]
+2θ−30,dE [(γθ0 γ̇d) (Uθ0)]
θ20,d + θ
I{p=q} + θ0,pθ0,qI{p 6=q}
+2θ0,pθ
0,dE [{γ̇pγ̇q + γθ0 γ̈p,q} (Uθ0)]− 2θ0,pθ0,qθ
γ̇2d + γθ0 γ̈d,d
(Uθ0)
8 LI WANG AND LIJIAN YANG
ψpq = 4E
γ̇p − θ0,pθ−10,dγ̇d
γ̇q − θ0,qθ−10,dγ̇d
(Uθ0) {γθ0 (Uθ0)− Y }
in which γ̇p and γ̈p,q are the values of
∂θp∂θq
γθ taking at θ = θ0, for any p, q = 1, 2, ..., d−1
and γθ is given in (2.6).
Remark 2.4. Consider the Generalized Linear Model (GLM): Y = g
XT θ0
+σ (X) ε, where
g is a known link function. Let θ̃ be the nonlinear least squared estimator of θ0 in GLM.
Theorem 2 shows that under the assumptions A1-A6, the asymptotic distribution of the θ̂−d
is the same as that of θ̃. This implies that our proposed SIP estimator θ̂−d is as efficient as if
the true link function g is known.
The next two propositions play an important role in our proof of the main results. Propo-
sition 2.1 establishes the uniform convergence rate of the derivatives of γ̂θ up to order 2 to
those of γθ in θ. Proposition 2.2 shows that the derivatives of the risk function up to order 2
are uniformly almost surely approximated by their empirical versions.
Proposition 2.1. Under Assumptions A2-A6, with probability 1
θ∈Sd−1c
u∈[0,1]
|γ̂θ (u)− γθ (u)| = O
log n+ h4
, (2.13)
1≤p≤d
θ∈Sd−1c
1≤i≤n
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
log n√
, (2.14)
1≤p,q≤d
θ∈Sd−1c
1≤i≤n
∂θp∂θq
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
log n√
. (2.15)
Proposition 2.2. Under Assumptions A2-A6, one has for k = 0, 1, 2
‖θ−d‖≤
∂kθ−d
R̂∗ (θ−d)−R∗ (θ−d)
= o(1), a.s..
Proofs of Theorem 2, Propositions 2.1 and 2.2 are given in Appendix.
3. Implementation
In this section, we will describe the actual procedure to implement the estimation of θ0
and g. We first introduce some new notation. For fixed θ, write the B-spline matrix as
Bθ = {Bj,4 (Uθ,i)}n, Ni=1,j=−3 and
Pθ = Bθ
θ (3.1)
as the projection matrix onto the cubic spline space Γ
n,θ. For any p = 1, ..., d, denote
Ḃp =
Bθ, Ṗp =
SINGLE-INDEX PREDICTION MODEL 9
as the first order partial derivatives of Bθ and Pθ with respect to θ.
Let Ŝ∗(θ−d) be the score vector of R̂
∗ (θ−d), i.e.
Ŝ∗(θ−d) =
R̂∗ (θ−d) . (3.2)
The next lemma provides the exact forms of Ŝ∗(θ−d).
Lemma 3.1. For the score vector of R̂∗ (θ−d) defined in (3.2), one has
Ŝ∗ (θ−d) = −n−1
ṖpY − θpθ−1d Y
, (3.3)
where for any p = 1, 2, ..., d
ṖpY = 2Y
T (I−Pθ) Ḃp
θ Y, (3.4)
where Ḃp =
{Bj,3 (Uθ,i)−Bj+1,3 (Uθ,i)} Ḟd (Xθ,i)h−1Xi,p
}n, N
i=1,j=−3
Ḟd (x) =
Γ (d+ 1)
aΓ {(d+ 1) /2}2 2d
I (|x| ≤ a) .
Proof. For any p = 1, 2, ..., d, the derivatives of B-splines in de Boor (2001) implies
Ḃp =
Bj,4 (Uθ,i)
}n, N
i=1,j=−3
Bj,4 (Uθ,i)
}n, N
i=1,j=−3
Bj,3 (Uθ,i)
tj+3 − tj
Bj+1,3 (Uθ,i)
tj+4 − tj+1
Ḟd (Xθ,i)Xi,p
}n, N
i=1,j=−3
{Bj,3 (Uθ,i)−Bj+1,3 (Uθ,i)} Ḟd (Xθ,i) h−1Xi,p
}n, N
i=1,j=−3
Next, note that
Ṗp = Ḃp
θ +Bθ
= Ḃp
θ +Bθ
θ +Bθ
Since
BTθ Bθ
BTθ Bθ
BTθ Bθ
θ Bθ +
)−1 ∂
BTθ Bθ
and ∂
BTθ Bθ
= ḂTp Bθ +B
θ Ḃp, thus
)−1 (
p Bθ +B
θ Ḃp
10 LI WANG AND LIJIAN YANG
Hence
Ṗp = (I−Pθ) Ḃp
θ +Bθ
p (I−Pθ) .
Thus, (3.4) follows immediately.
In practice, the estimation is implemented via the following procedure.
Step 1. Standardize the predictor vectors {Xi}ni=1 and for each fixed θ ∈ Sd−1c obtain the
CDF transformed variables {Uθ,i}ni=1 of the SIP variable {Xθ,i}
through formula (2.5), where
the radius a is taken to be the 95% percentile of {‖Xi‖}ni=1.
Step 2. Compute quadratic and cubic B-spline basis at each value Uθ,i, where the number
of interior knots N is
N = min
n1/5.5
, (3.5)
Step 3. Find the estimator θ̂ of θ0 by minimizing R̂
∗ through the port optimization routine
with (0, 0, ..., 1)
as the initial value and the empirical score vector Ŝ∗ in (3.3). If d < n, one
can take the simple LSE (without the intercept) for data {Yi,Xi}ni=1 with its last coordinate set
positive.
Step 4. Obtain the spline estimator ĝ of g by plugging θ̂ obtained in Step 3 into (2.10).
Remark 3.1. In (3.5), c1 and c2 are positive integers and [ν] denotes the integer part of ν. The
choice of the tuning parameter c1 makes little difference for a large sample and according to our
asymptotic theory there is no optimal way to set these constants. We recommend using c1 = 1
to save computing for massive data sets. The first term ensures Assumption A6. The addition
constrain c2 can be taken from 5 to 10 for smooth monotonic or smooth unimodel regression
and c2 > 10 if has many local minima and maxima, which is very unlikely in application.
4. Simulations
In this section, we carry out two simulations to illustrate the finite-sample behavior of our
SIP estimation method. The number of interior knots N is computed according to (3.5) with
c1 = 1, c2 = 5. All of our codes have been written in R.
Example 1. Consider the model in Xia, Li, Tong and Zhang (2004)
Y = m (X) + σ0ε, σ0 = 0.3, 0.5, ε
i.i.d∼ N(0, 1)
where X = (X1,X2)
T ∼N(0, I2), truncated by [−2.5, 2.5]2 and
m (x) = x1 + x2 + 4exp
− (x1 + x2)2
x21 + x
. (4.1)
If δ = 0, then the underlying true functionm is a single-index function, i.e., m (X) =
2XT θ0+
XT θ0
, where θT0 = (1, 1) /
2. While δ 6= 0, then m is not a genuine single-index
SINGLE-INDEX PREDICTION MODEL 11
function. An impression of the bivariate function m for δ = 0 and δ = 1 can be gained in
Figure 1 (a) and (b), respectively.
Table 1: Report of Example 1 (Values out/in parentheses: δ = 0/δ = 1)
σ0 n θ0 BIAS SD MSE Average MSE
5e− 04 0.00825 7e− 05
(−0.00236) (0.02093) (0.00044) 7e− 05
−6e− 04 0.00826 7e− 05 (0.00043)
(0.00174) (0.02083) (0.00043)
−0.00124 0.00383 2e− 05
(−0.00129) (0.01172) (0.00014) 2e− 05
−0.00124 0.00383 2e− 05 (0.00014)
(0.00110) (0.01160) (0.00013)
0.00121 0.01346 0.00018
(−0.00137) (0.02257) (0.00051) 0.00018
−0.00147 0.01349 0.00018 (0.00051)
(0.00062) (0.02309) (0.00052)
−0.00204 0.00639 4e− 05
(−0.00229) (0.01205) (0.00015) 4e− 05
0.00197 0.00637 4e− 05 (0.00015)
(0.00208) (0.01190) (0.00014)
For δ = 0, 1, we draw 100 random realizations of each sample size n = 50, 100, 300 respec-
tively. To demonstrate how close our SIP estimator is to the true index parameter θ0, Table 1
lists the sample mean (MEAN), bias (BIAS), standard deviation (SD), the mean squared error
(MSE) of the estimates of θ0 and the average MSE of both directions. From this table, we find
that the SIP estimators are very accurate for both cases δ = 0 and δ = 1, which shows that
our proposed method is robust against the deviation from single-index model. As we expected,
when the sample size increases, the SIP coefficient is more accurately estimated. Moreover, for
n = 100, 300, the total average is inversely proportional to n.
Example 2. Consider the heteroscedastic regression model (1.1) with
m (X) = sin
, σ (X) = σ0
5− exp
5 + exp
) , (4.2)
in which Xi = {Xi,1, ...,Xi,d}T and εi, i = 1, ..., n, are
i.i.d∼ N (0, 1), σ0 = 0.2. In our simulation,
the true parameter θT0 = (1, 1, 0, ..., 0, 1)/
3 for different sample size n and dimension d. The
12 LI WANG AND LIJIAN YANG
superior performance of SIP estimators is borne out in comparison with MAVE of Xia, Tong,
Li and Zhu (2002). We also investigate the behavior of SIP estimators in the previously
unemployed cases that sample size n is smaller than or equal to d, for instance, n = 100, d =
100, 200 and n = 200, d = 200, 400. The average MSEs of the d dimensions are listed in Table
2, from which we see that the performance of the SIP estimators are quite reasonable and in
most of the scenarios n ≤ d, the SIP estimators still work astonishingly well where the MAVEs
become unreliable. For n = 100, d = 10, 50, 100, 200, the estimates of the link prediction
function g from model (4.2) are plotted in Figure 2, which is rather satisfactory even when
dimension d exceeds the sample size n.
Theorem 1 indicates that θ̂−d is strongly consistent of θ0,−d. To see the convergence, we
run 100 replications and in each replication, the value of ‖θ̂ − θ0‖/
d is computed. Figure
3 plots the kernel density estimations of the 100 ‖θ̂ − θ0‖ in Example 2, in which dimension
d = 10, 50, 100, 200. There are four types of line characteristics which correspond to the two
sample sizes, the dotted-dashed line (n = 100), dotted line (n = 200), dashed line (500) and
solid line (n = 1000). As sample sizes increasing, the squared errors are becoming closer to 0,
with narrower spread out, confirmative to the conclusions of Theorem 1.
Lastly, we report the average computing time of Example 2 to generate one sample of size
n and perform the SIP or MAVE procedure done on the same ordinary Pentium IV PC in
Table 2. From Table 2, one sees that our proposed SIP estimator is much faster than the
MAVE. The computing time for MAVE is extremely sensitive to sample size as we expected.
For very large d, MAVE becomes unstable to the point of the breaking down in four cases.
5. An application
In this section we demonstrate the proposed SIP model through the river flow data of
Jökulsá Eystri River of Iceland, from January 1, 1972 to December 31, 1974. There are 1096
observations, see Tong (1990). The response variables are the daily river flow (Yt), measured in
meter cubed per second of Jökulsá Eystri River. The exogenous variables are temperature (Xt)
in degrees Celsius and daily precipitation (Zt) in millimeters collected at the meteorological
station at Hveravellir.
This data set was analyzed earlier through threshold autoregressive (TAR) models by
Tong, Thanoon and Gudmundsson (1985), Tong (1990), and nonlinear additive autoregressive
(NAARX) models by Chen and Tsay (1993). Figure 4 shows the plots of the three time series,
from which some nonlinear and non-stationary features of the river flow series are evident. To
make these series stationary, we remove the trend by a simple quadratic spline regression and
these trends (dashed lines) are shown in Figure 4. By an abuse of notation, we shall continue
to use Xt, Yt, Zt to denote the detrended series.
SINGLE-INDEX PREDICTION MODEL 13
In the analysis, we pre-select all the lagged values in the last 7 days (1 week), i.e., the
predictor pool is {Yt−1, ..., Yt−7,Xt,Xt−1, ...,Xt−7, Zt, Zt−1, ..., Zt−7, }. Using BIC similar to
Huang and Yang (2004) for our proposed spline SIP model with 3 interior knots, the following
9 explanatory variables are selected from the above set {Yt−1, ..., Yt−4,Xt,Xt−1,Xt−2, Zt, Zt−1}.
Based on this selection, we fit the SIP model again and obtain the estimate of the SIP coefficient
θ̂ = {−0.877, 0.382,−0.208, 0.125,−0.046,−0.034, 0.004,−0.126, 0.079}T . Figure 5 (a) and (b)
display the fitted river flow series and the residuals against time.
Next we examine the forecasting performance of the SIP method. We start with estimating
the SIP estimator using only observations of the first two years, then we perform the out-of-
sample rolling forecast of the entire third year. The observed values of the exogenous variables
are used in the forecast. Figure 5 (c) shows this SIP out-of-sample rolling forecasts. For the
purpose of comparison, we also try the MAVE method, in which the same predictor vector is
selected by using BIC. The mean squared prediction error is 60.52 for the SIP model, 61.25 for
MAVE, 65.62 for NAARX, 66.67 for TAR and 81.99 for the linear regression model, see Chen
and Tsay (1993). Among the above five models, the SIP model produces the best forecasts.
6. Conclusion
In this paper we propose a robust SIP model for stochastic regression under weak depen-
dence regardless if the underlying function is exactly a single-index function. The proposed
spline estimator of the index coefficient possesses not only the usual strong consistency and
n-rate asymptotically normal distribution, but also is as efficient as if the true link function
g is known. By taking advantage of the spline smoothing method and the iterative method,
the proposed procedure is much faster than the MAVE method. This procedure is especially
powerful for large sample size n and high dimension d and unlike the MAVE method, the
performance of the SIP remains satisfying in the case d > n.
Acknowledgment
This work is part of the first author’s dissertation under the supervision of the second
author, and has been supported in part by NSF award DMS 0405330.
Appendix
A.1. Preliminaries
In this section, we introduce some properties of the B-spline.
Lemma A.1. There exist constants c > 0 such that for
j=−k+1 αj,kBj,k up to order k = 4
ch1/r ‖α‖r ≤
j=−k+1 αj,kBj,k
3r−1h
)1/r ‖α‖r , 1 ≤ r ≤ ∞
ch1/r ‖α‖r ≤
j=−k+1 αj,kBj,k
≤ (3h)1/r ‖α‖r , 0 < r < 1
14 LI WANG AND LIJIAN YANG
where α := (α−1,2, α0,2, ..., αN,2, ..., αN,4). In particular, under Assumption A2, for any fixed θ
ch1/2 ‖α‖2 ≤
j=−k+1
αj,kBj,k
≤ Ch1/2 ‖α‖2 .
Proof. It follows from the B-spline property on page 96 of de Boor (2001),
j=−k+1Bj,k ≡
3 on [0, 1]. So the right inequality follows immediate for r = ∞. When 1 ≤ r < ∞, we use
Hölder’s inequality to find
j=−k+1
αj,kBj,k
j=−k+1
|αj,k|r Bj,k
j=−k+1
1−1/r
= 31−1/r
j=−k+1
|αj,k|r Bj,k
Since all the knots are equally spaced,
−∞Bj,k (u) du ≤ h, the right inequality follows from
j=−k+1
αj,kBj,k (u)
du ≤ 3r−1h ‖α‖rr .
When r < 1, we have
j=−k+1
αj,kBj,k
j=−k+1
|αj,k|r Brj,k.
Since
j,k (u) du ≤ tj+k − tj = kh and
j=−k+1
αj,kBj,k (u)
du ≤ ‖α‖rr
Brj,k (u) du ≤ 3h ‖α‖
the right inequality follows in this case as well. For the left inequalities, we derive from Theorem
5.4.2, DeVore and Lorentz (1993)
|αj,k| ≤ C1h−1/r
∫ tj+1
j=−k+1
αj,kBj,k (u)
for any 0 < r ≤ ∞, so
|αj,k|r ≤ Cr1h−1
∫ tj+1
j=−k+1
αj,kBj,k (u)
SINGLE-INDEX PREDICTION MODEL 15
Since each u ∈ [0, 1] appears in at most k intervals (tj,tj+k), adding up these inequalities, we
obtain that
‖α‖rr ≤ C1h
∫ tj+k
j=−k+1
αj,kBj,k (u)
du ≤ 3Ch−1
j=−k+1
αj,kBj,k
The left inequality follows.
For any functions φ and ϕ, define the empirical inner product and the empirical norm as
〈φ,ϕ〉θ =
φ (u)ϕ (u) fθ (u) du, ‖φ‖22,n,θ = n
φ2 (Uθ,i) .
In addition, if functions φ,ϕ are L2 [0, 1]-integrable, define the theoretical inner product and
its corresponding theoretical L2 norm as
‖φ‖22,θ =
φ2 (u) fθ (u) du, 〈φ,ϕ〉n,θ = n
φ (Uθ,i)ϕ (Uθ,i) .
Lemma A.2. Under Assumptions A2, A5 and A6, with probability 1,
θ∈Sd−1c
k,k′=2,3,4
1≤j,j′≤N
Bj,k, Bj′,k′
Bj,k, Bj′,k′
log n
Proof. We only prove the case k = k′ = 4, all other cases are similar. Let
ζθ,j,j′,i = Bj,4 (Uθ,i)Bj′,4 (Uθ,i)− EBj,4 (Uθ,i)Bj′,4 (Uθ,i) ,
with the second moment
Eζ2θ,j,j′,i = E
B2j,4 (Uθ,i)B
j′,4 (Uθ,i)
EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
where
EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
}2 ∼ N−2, E
B2j,4 (Uθ,i)B
j′,4 (Uθ,i)
∼ N−1 by Assumption A2.
Hence, Eζ2θ,j,j′,i ∼ N−1. The k-th moment is given by
∣ζθ,j,j′,i
∣Bj,4 (Uθ,i)Bj′,4 (Uθ,i)− EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
≤ 2k−1
∣Bj,4 (Uθ,i)Bj′,4 (Uθ,i)
∣EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
where
∣EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
k ∼ N−k, E
∣EBj,4 (Uθ,i)Bj′,4 (Uθ,i)
k ∼ N−1. Thus, there exists
a constant C > 0 such that E
∣ζθ,j,j′,i
k ≤ C2k−1k!Eζ2j,j′,i. So the Cramér’s condition is satisfied
with Cramér’s constant c∗. By the Bernstein’s inequality (see Bosq (1998), Theorem 1.4, page
31), we have for k = 3
ζθ,j,j′,i
≤ a1 exp
25m22 + 5c
+ a2 (k)α
q + 1
])6/7
16 LI WANG AND LIJIAN YANG
where
δn = δ
log n√
, a1 = 2
δ2 (nN)
log2 n
25m22 + 5c
, m22 ∼ N−1,
a2 (3) = 11n
, m3 = max
1≤i≤n
∥ζθ,j,j′,i
≤ cN1/3.
Observe that 5cδn = o(1) by Assumption A6, then by taking q such that
≥ c0 log n,
q ≥ c1n/ log n for some constants c0, c1, one has a1 = O(n/q) = O (log n), a2 (3) = o
Assumption A6 again. Assumption A5 yields that
q + 1
])6/7
K0 exp
q + 1
])}6/7
≤ Cn−6λ0c0/7.
Thus, for fixed θ ∈ Sd−1c , when n large enough
ζθ,j,j′,i
≤ c log n exp
−c2δ2 log n
+ Cn2−6λ0c0/7. (A.1)
We divide each range of θp, p = 1, 2, ..., d − 1, into n6/(d−1) equally spaced intervals with
disjoint endpoints −1 = θp,0 < θp,1 < ... < θp,Mn = 1, for p = 1, ..., d − 1. Projecting these
small cylinders onto Sd−1c , the radius of each patch Λr, r = 1, ...,Mn is bounded by cM
Denote the projection of the Mn points as θr =
θr,−d,
1− ‖θr,−d‖22
, r = 0, 1, ...,Mn.
Employing the discretization method, sup
θ∈Sd−1c
1≤j,j′≤N
∣ζθ,j,j′,i
∣ is bounded by
0≤r≤Mn
1≤j,j′≤N
∣ζθr,j,j′,i
∣+ sup
0≤r≤Mn
1≤j,j′≤N
∣ζθ,j,j′,i − ζθr,j,j′,i
∣ . (A.2)
By (A.1) and Assumption A6, there exists large enough value δ > 0 such that
ζθr,j,j′,i
≤ n−10,
which implies that
1≤j,j′≤N
ζθr ,j,j′,i
N2Mnn
−10 ≤ C
n−3 <∞.
Thus, Borel-Cantelli Lemma entails that
0≤r≤Mn
1≤j,j′≤N
ζθr ,j,j′,i
log n√
, a.s.. (A.3)
SINGLE-INDEX PREDICTION MODEL 17
Employing Lipschitz continuity of the cubic B-spline, one has with probability 1
0≤r≤Mn
1≤j,j′≤N
ζθ,j,j′,i − ζθr,j,j′,i
M−1n h
−6) . (A.4)
Therefore Assumption A2, (A.2), (A.3) and (A.4) lead to the desired result.
Denote by Γ = Γ(0)∪Γ(1)∪Γ(2) the space of all linear, quadratic and cubic spline functions
on [0, 1]. We establish the uniform rate at which the empirical inner product approximates the
theoretical inner product for all B-splines Bj,k with k = 2, 3, 4.
Lemma A.3. Under Assumptions A2, A5 and A6, one has
An = sup
θ∈Sd−1c
γ1,γ2∈Γ
〈γ1, γ2〉n,θ − 〈γ1, γ2〉θ
‖γ1‖2,θ ‖γ2‖2,θ
log n
, a.s.. (A.5)
Proof. Denote without loss of generality,
j=−k+1
αjkBj,k, γ2 =
j=−k+1
βjkBj,k,
for any two 3 (N + 3)-vectors
α =(α−1,2, α0,2, ..., αN,2, ..., αN,4) , β =(β−1,2, β0,2, ..., βN,2, ..., βN,4) .
Then for fixed θ
〈γ1, γ2〉n,θ =
j=−k+1
αj,kBj,k (Uθ,i)
j=−k+1
βj,kBj,k (Uθ,i)
j=−k+1
j′=−k+1
αj,kβj′,k′
Bj,k, Bj′,k′
‖γ1‖22,θ =
j=−k+1
j′=−k+1
αj,kαj′,k′
Bj,k, Bj′ ,k′
‖γ2‖22,θ =
j=−k+1
j′=−k+1
βj,kβj′,k′
Bj,k, Bj′ ,k′
According to Lemma A.1, one has for any θ ∈ Sd−1c ,
c1h ‖α‖22 ≤ ‖γ1‖
2,θ ≤ c2h ‖α‖
2 , c1h ‖β‖
2 ≤ ‖γ2‖
2,θ ≤ c2h ‖β‖
c1h ‖α‖2 ‖β‖2 ≤ ‖γ1‖2,θ ‖γ2‖2,θ ≤ c2h ‖α‖2 ‖β‖2 .
18 LI WANG AND LIJIAN YANG
Hence
An = sup
θ∈Sd−1c
γ1∈γ,γ2∈Γ
〈γ1, γ2〉n,θ − 〈γ1, γ2〉θ
‖γ1‖2,θ ‖γ2‖2,θ
‖α‖∞ ‖β‖∞
c1h ‖α‖2 ‖β‖2
× sup
θ∈Sd−1c
k,k′=2,3,4
1≤j,j′≤N
Bj,k, Bj′ ,k′
Bj,k, Bj′ ,k′
An ≤ c0h−1 sup
θ∈Sd−1c
k,k′=2,3,4
1≤j,j′≤N
Bj,k, Bj′ ,k′
Bj,k, Bj′ ,k′
which, together with Lemma A.2, imply (A.5).
A.2. Proof of Proposition 2.1
For any fixed θ, we write the response YT = (Y1, ..., Yn) as the sum of a signal vector γθ,
a parametric noise vector Eθ and a systematic noise vector E, i.e.,
Y = γθ +Eθ +E,
in which the vectors γTθ = {γθ (Uθ,1) , ..., γθ (Uθ,n)}, ET = {σ (X1) ε1, ..., σ (Xn) εn} and ETθ =
{m (X1)− γθ (Uθ,1) , ...,m (Xn)− γθ (Uθ,n)}.
Remark A.1. If m is a genuine single-index function, then Eθ0 ≡ 0, thus the proposed SIP
model is exactly the single-index model.
Let Γ
n, θ be the cubic spline space spanned by {Bj,4 (Uθ,i)}
, −3 ≤ j ≤ N for fixed θ.
Projecting Y onto Γ
n, θ yields that
γ̂θ = {γ̂θ (Uθ,1) , ..., γ̂θ (Uθ,n)}T = ProjΓ(2)
γθ + ProjΓ(2)
Eθ + ProjΓ(2)
where γ̂θ is given in (2.8). We break the cubic spline estimation error γ̂θ (uθ) − γθ (uθ) into a
bias term γ̃θ (uθ)− γθ (uθ) and two noise terms ε̃θ (uθ) and ε̂θ (uθ)
γ̂θ (uθ)− γθ (uθ) = {γ̃θ (uθ)− γθ (uθ)}+ ε̃θ (uθ) + ε̂θ (uθ) , (A.6)
where
γ̃θ (u) = {Bj,4 (u)}T−3≤j≤N V
〈γθ, Bj,4〉n,θ
, (A.7)
ε̃θ (u) = {Bj,4 (u)}T−3≤j≤N V
〈Eθ, Bj,4〉n,θ
, (A.8)
ε̂θ (u) = {Bj,4 (u)}T−3≤j≤N V
〈E, Bj,4〉n,θ
. (A.9)
SINGLE-INDEX PREDICTION MODEL 19
In the above, we denote by Vn,θ the empirical inner product matrix of the cubic B-spline basis
and similarly, the theoretical inner product matrix as Vθ
Vn,θ =
θ Bθ =
Bj′,4, Bj,4
j,j′=−3
,Vθ =
Bj′,4, Bj,4
j,j′=−3 . (A.10)
In Lemma A.5, we provide the uniform upper bound of
∞. Before
that, we first describe a special case of Theorem 13.4.3 in DeVore and Lorentz (1993).
Lemma A.4. If a bi-infinite matrix with bandwidth r has a bounded inverse A−1 on l2 and
κ = κ (A) := ‖A‖2
is the condition number of A, then
∞ ≤ 2c0 (1− ν)
, with
c0 = ν
−2r ∥
, ν =
κ2 − 1
)1/4r (
κ2 + 1
)−1/4r
Lemma A.5. Under Assumptions A2, A5 and A6, there exist constants 0 < cV < CV such
that cVN
−1 ‖w‖22 ≤ wTVθw ≤ CVN−1 ‖w‖
2 and
−1 ‖w‖22≤ w
Vn,θw ≤ CVN−1 ‖w‖22 , a.s., (A.11)
with matrices Vθ and Vn,θ defined in (A.10). In addition, there exists a constant C > 0 such
θ∈Sd−1c
≤ CN, a.s., sup
θ∈Sd−1c
∞ ≤ CN. (A.12)
Proof. First we compute the lower and upper bounds for the eigenvalues of Vn,θ. Let w be any
(N + 4)-vector and denote γw (u) =
j=−3wjBj,4 (u), then Bθw = {γw (Uθ,1) , ..., γw (Uθ,n)}
and the definition of An in (A.5) from Lemma A.3 entails that
‖γw‖22,θ (1−An) ≤ w
Vn,θw = ‖γw‖22,n,θ ≤ ‖γw‖
2,θ (1 +An) . (A.13)
Using Theorem 5.4.2 of DeVore and Lorentz (1993) and Assumption A2, one obtains that
‖w‖22 ≤ ‖γw‖
2,θ = w
Vθw =
wjBj,4
‖w‖22 , (A.14)
which, together with (A.13), yield
−1 ‖w‖22 (1−An) ≤ w
Vn,θw ≤ CfCN−1 ‖w‖22 (1 +An) . (A.15)
Now the order of An in (A.5), together with (A.14) and (A.15) implies (A.11), in which cV =
cfC,CV = CfC. Next, denote by λmax (Vn,θ) and λmin (Vn,θ) the maximum and minimum
eigenvalue of Vn,θ, simple algebra and (A.11) entail that
−1 ≥ ‖Vn,θ‖2 = λmax (Vn,θ) ,
= λ−1min (Vn,θ) ≤ c
V N, a.s.,
20 LI WANG AND LIJIAN YANG
κ := ‖Vn,θ‖2
= λmax (Vn,θ)λ
min (Vn,θ) ≤ CV c
V <∞, a.s..
Meanwhile, let wj = the (N + 4)-vector with all zeros except the j-th element being 1, j =
−3, ..., N . Then clearly
j Vn,θwj =
B2j,4 (Uθ,i) = ‖Bj,4‖
, ‖wj‖2 = 1,−3 ≤ j ≤ N
and in particular
0 Vn,θw0 ≤ λmax (Vn,θ) ‖w0‖2 = λmax (Vn,θ) ,
−3Vn,θw−3 ≥ λmin (Vn,θ) ‖w−3‖2 = λmin (Vn,θ) .
This, together with (A.5) yields that
κ = λmax (Vn,θ)λ
min (Vn,θ) ≥
wT0 Vn,θw0
wT−3Vn,θw−3
‖B0,4‖2n,θ
‖B−3,4‖2n,θ
‖B0,4‖2θ
‖B−3,4‖2θ
1 +An
which leads to κ ≥ C > 1, a.s. because the definition of B-spline and Assumption A2 ensure
that ‖B0,4‖2θ ≥ C0 ‖B−3,4‖
for some constant C0 > 1. Next applying Lemma A.4 with
κ2 − 1
)1/16 (
κ2 + 1
)−1/16
and c0 = ν
, one gets
≤ 2ν−8N (1− ν)−1 =
CN, a.s.. Hence part one of (A.12) follows. Part two of (A.12) is proved in the same fashion.
In the following, we denote by QT (m) the 4-th order quasi-interpolant of m corresponding
to the knots T , see equation (4.12), page 146 of DeVore and Lorentz (1993). According to
Theorem 7.7.4, DeVore and Lorentz (1993), the following lemma holds.
Lemma A.6. There exists a constant C > 0, such that for 0 ≤ k ≤ 2 and γ ∈ C(4) [0, 1]
(γ −QT (γ))(k)
h4−k,
Lemma A.7. Under Assumptions A2, A3, A5 and A6, there exists an absolute constant C > 0,
such that for function γ̃θ (u) in (A.7)
θ∈Sd−1c
(γ̃θ − γθ)
h4−k, a.s., 0 ≤ k ≤ 2, (A.16)
Proof. According to Theorem A.1 of Huang (2003), there exists an absolute constant C > 0,
such that
θ∈Sd−1c
‖γ̃θ − γθ‖∞ ≤ C sup
θ∈Sd−1c
γ∈Γ(2)
‖γ − γθ‖∞ ≤ C
h4, a.s., (A.17)
SINGLE-INDEX PREDICTION MODEL 21
which proves (A.16) for the case k = 0. Applying Lemma A.6, one has for 0 ≤ k ≤ 2
θ∈Sd−1c
{QT (γθ)− γθ}
≤ C sup
θ∈Sd−1c
h4−k ≤ C
h4−k, (A.18)
As a consequence of (A.17) and (A.18) for the case k = 0, one has
θ∈Sd−1c
‖QT (γθ)− γ̃θ‖∞ ≤ C
h4, a.s.,
which, according to the differentiation of B-spline given in de Boor (2001), entails that
θ∈Sd−1c
{QT (γθ)− γ̃θ}
h4−k, a.s., 0 ≤ k ≤ 2. (A.19)
Combining (A.18) and (A.19) proves (A.16) for k = 1, 2.
Lemma A.8. Under Assumptions A1, A2, A4 and A5, there exists an absolute constant C > 0,
such that
1≤p≤d
θ∈Sd−1c
{γ̃θ (Uθ,i)− γθ (Uθ,i)}ni=1
h3, a.s., (A.20)
1≤p,q≤d
θ∈Sd−1c
∂θp∂θq
{γ̃θ (Uθ,i)− γθ (Uθ,i)}ni=1
h2, a.s.. (A.21)
Proof. According to the definition of γ̃θ in (A.7), and the fact that QT (γθ) is a cubic spline
on the knots T
{{QT (γθ)− γ̃θ} (Uθ,i)}ni=1 = Pθ {{QT (γθ)− γθ} (Uθ,i)}
which entails that
{{QT (γθ)− γ̃θ} (Uθ,i)}ni=1 =
Pθ {{QT (γθ)− γθ} (Uθ,i)}ni=1
= Ṗp {{QT (γθ)− γθ} (Uθ,i)}ni=1 +Pθ
{{QT (γθ)− γθ} (Uθ,i)}ni=1 .
Since
{{QT (γθ)− γθ} (Uθ,i)}ni=1 =
(Uθ,i)
{QT (γθ)− γθ} (Uθ,i)Xip
applying (A.19) to the decomposition above produces (A.20). The proof of (A.21) is similar.
22 LI WANG AND LIJIAN YANG
Lemma A.9. Under Assumptions A2, A5 and A6, there exists a constant C > 0 such that
θ∈Sd−1c
∥n−1BTθ
∞ ≤ Ch, a.s., sup
1≤p≤d
θ∈Sd−1c
n−1ḂTp
≤ C, a.s., (A.22)
θ∈Sd−1c
‖Pθ‖∞ ≤ C, a.s., sup
1≤p≤d
θ∈Sd−1c
≤ Ch−1, a.s.. (A.23)
Proof. To prove (A.22), observe that for any vector a ∈ Rn, with probability 1
∥n−1BTθ a
∞ ≤ ‖a‖∞ max−3≤j≤N
Bj,4 (Uθ,i)
≤ Ch ‖a‖∞ ,
n−1ḂTp a
≤ ‖a‖∞ max−3≤j≤N
{(Bj,3 −Bj+1,3) (Uθ,i)} Ḟd (Xθ,i)Xi,p
≤ C ‖a‖∞ .
To prove (A.23), one only needs to use (A.12), (A.22) and (3.1).
Lemma A.10. Under Assumptions A2 and A4-A6, one has with probability 1
θ∈Sd−1c
BTθ E
= max
−3≤j≤N
Bj,4 (Uθ,i)σ (Xi) εi
log n√
, (A.24)
1≤p≤d
θ∈Sd−1c
BTθ E
= sup
1≤p≤d
θ∈Sd−1c
ḂTpE
log n√
. (A.25)
Similarly, under Assumptions A2, A4-A6, with probability 1
θ∈Sd−1c
BTθ Eθ
= sup
θ∈Sd−1c
−3≤j≤N
Bj,4 (Uθ,i) {m (Xi)− γθ (Uθ,i)}
log n√
(A.26)
1≤p≤d
θ∈Sd−1c
BTθ Eθ
log n√
, a.s.. (A.27)
Proof. We decompose the noise variable εi into a truncated part and a tail part εi = ε
i,1 +
εDni,2 +m
i , where Dn = n
η (1/3 < η < 2/5), εDni,1 = εiI {|εi| > Dn},
εDni,2 = εiI {|εi| ≤ Dn} −m
i = E [εiI {|εi| ≤ Dn} |Xi] .
It is straightforward to verify that the mean of the truncated part is uniformly bounded by
D−2n , so the boundedness of B spline basis and of the function σ
2 entail that
θ∈Sd−1c
Bj,4 (Uθ,i) σ (Xi)m
n−2/3
SINGLE-INDEX PREDICTION MODEL 23
The tail part vanishes almost surely
P {|εn| > Dn} ≤
D−3n <∞.
Borel-Cantelli Lemma implies that
Bj,4 (Uθ,i) σ (Xi) ε
, for any k > 0.
For the truncated part, using Bernstein’s inequality and discretization as in Lemma A.2
θ∈Sd−1c
1≤j≤N
Bj,4 (Uθ,i) σ (Xi) ε
log n/
, a.s..
Therefore (A.24) is established as with probability 1
θ∈Sd−1c
n−2/3
log n/
log n/
The proofs of (A.25), (A.26) are similar as E {m (Xi)− γθ (Uθ,i) |Uθ,i } ≡ 0, but no truncation
is needed for (A.26) as sup
θ∈Sd−1c
1≤i≤n
|m (Xi)− γθ (Uθ,i)| ≤ C <∞. Meanwhile, to prove (A.27),
we note that for any p = 1, ..., d
[Bj,4 (Uθ,i) {m (Xi)− γθ (Uθ,i)}]
According to (2.6), one has γθ (Uθ) ≡ E {m (X) |Uθ}, hence
E [Bj,4 (Uθ) {m (X)− γθ (Uθ)}] ≡ 0,−3 ≤ j ≤ N, θ ∈ Sd−1c .
Applying Assumptions A2 and A3, one can differentiate through the expectation, thus
[Bj,4 (Uθ) {m (X)− γθ (Uθ)}]
≡ 0, 1 ≤ p ≤ d,−3 ≤ j ≤ N, θ ∈ Sd−1c ,
which allows one to apply the Bernstein’s inequality to obtain that with probability 1
[Bj,4 (Uθ,i) {m (Xi)− γθ (Uθ,i)}]
(nh)−1/2 log n
which is (A.27).
Lemma A.11. Under Assumptions A2 and A4-A6, for ε̂θ (u) in (A.9), one has
θ∈Sd−1c
u∈[0,1]
|ε̂θ (u)| = O
log n
, a.s.. (A.28)
24 LI WANG AND LIJIAN YANG
Proof. Denote â ≡ (â−3, · · · , âN )T =
BTθ Bθ
BTθ E = V
n−1BTθ E
, then ε̂θ (u) =
j=−3 âjBj,4 (u), so the order of ε̂θ (u) is related to that of â. In fact, by Theorem 5.4.2
in DeVore and Lorentz (1993)
θ∈Sd−1c
u∈[0,1]
|ε̂θ (u)| ≤ sup
θ∈Sd−1c
‖â‖∞ =
θ∈Sd−1c
n−1BTθ E
≤ CN sup
θ∈Sd−1c
∥n−1BTθ E
∞ , a.s.,
where the last inequality follows from (A.12) of Lemma A.5. Applying (A.24) of Lemma A.10,
we have established (A.28).
Lemma A.12. Under Assumptions A2 and A4-A6, for ε̃θ (u) in (A.8), one has
θ∈Sd−1c
u∈[0,1]
|ε̃θ (u)| = O
log n
, a.s.. (A.29)
The proof is similar to Lemma A.11, thus omitted.
The next result evaluates the uniform size of the noise derivatives.
Lemma A.13. Under Assumptions A2-A6, one has with probability 1
1≤p≤d
θ∈Sd−1c
1≤i≤n
ε̂θ (Uθ,i)
(nh3)−1/2 log n
, (A.30)
1≤p≤d
θ∈Sd−1c
1≤i≤n
ε̃θ (Uθ,i)
(nh3)−1/2 log n
, (A.31)
1≤p,q≤d
θ∈Sd−1c
1≤i≤n
∂θp∂θq
ε̂θ (Uθ,i)
(nh5)−1/2 log n
, (A.32)
1≤p,q≤d
θ∈Sd−1c
1≤i≤n
∂θp∂θq
ε̃θ (Uθ,i)
(nh5)−1/2 log n
. (A.33)
Proof. Note that
ε̂θ (Uθ,i)
= (I−Pθ) Ḃp
θ E+Bθ
p (I−Pθ)E.
Applying (A.24) and (A.25) of Lemma A.10, (A.12) of Lemma A.5, (A.22) and (A.23) of
Lemma A.9, one derives (A.30). To prove (A.31), note that
ε̃θ (Uθ,i)
{PθEθ} = ṖpEθ +Pθ
Eθ = T1 + T2, (A.34)
in which
(I−Pθ) Ḃp −Bθ
(I−Pθ) Ḃp −Bθ
BTθ Bθ
ḂTpBθ
BTθ Bθ
BTθ Eθ
SINGLE-INDEX PREDICTION MODEL 25
T2 = Bθ
BTθ Bθ
BTθ Eθ
By (A.24), (A.12), (A.22) and (A.23), one derives
θ∈Sd−1c
‖T1‖∞ = O
n−1/2N3/2 log n
, a.s., (A.35)
while (A.27) of Lemma A.10, (A.12) of Lemma A.5
θ∈Sd−1c
‖T2‖∞ = N ×O
n−1/2h−1/2 log n
n−1/2h−3/2 log n
, a.s.. (A.36)
Now, putting together (A.34), (A.35) and (A.36), we have established (A.31). The proof for
(A.32) and (A.33) are similar.
Proof of Proposition 2.1. According to the decomposition (A.6)
|γ̂θ (u)− γθ (u)| = |{γ̃θ (u)− γθ (u)}+ ε̃θ (u) + ε̂θ (u)| .
Then (2.13) follows directly from (A.16) of Lemma A.7, (A.28) of Lemma A.11 and (A.29) of
Lemma A.12. Again by definitions (A.8) and (A.9), we write
{(γ̂θ − γθ) (Uθ,i)} =
(γ̃θ − γθ) (Uθ,i) +
γ̃θ (Uθ,i) +
ε̂θ (Uθ,i) .
It is clear from (A.20), (A.30) and (A.31) that with probability 1
1≤p≤d
θ∈Sd−1c
1≤i≤n
(γ̃θ − γθ) (Uθ,i)
1≤p≤d
θ∈Sd−1c
1≤i≤n
ε̃θ (Uθ,i)
ε̂θ (Uθ,i)
)−1/2
log n
Putting together all the above yields (2.14). The proof of (2.15) is similar.
A.3. Proof of Proposition 2.2
Lemma A.14. Under Assumptions A2-A6, one has
θ∈Sd−1c
R̂ (θ)−R (θ)
= o(1), a.s..
Proof. For the empirical risk function R̂ (θ) in (2.9), one has
R̂ (θ) = n−1
{γ̂θ (Uθ,i)−m (Xi)− σ (Xi) εi}2
26 LI WANG AND LIJIAN YANG
= n−1
{γ̂θ (Uθ,i)− γθ (Uθ,i) + γθ (Uθ,i)−m (Xi)− σ (Xi) εi}2 ,
hence
R̂ (θ) = n−1
{γ̂θ (Uθ,i)− γθ (Uθ,i)}2 + n−1
σ2 (Xi) ε
+2n−1
{γ̂θ (Uθ,i)− γθ (Uθ,i)} {γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
{γθ (Uθ,i)−m (Xi)}2 + 2n−1
{γθ (Uθ,i)−m (Xi)}σ (Xi) εi,
where γ̂θ (x) is defined in (2.8). Using the expression of R (θ) in (2.7), one has
θ∈Sd−1c
∣R̂ (θ)−R (θ)
∣ ≤ I1 + I2 + I3 + I4,
I1 = sup
θ∈Sd−1c
{γ̂θ (Uθ,i)− γθ (Uθ,i)}2
I2 = sup
θ∈Sd−1c
{γ̂θ (Uθ,i)− γθ (Uθ,i)} {γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
I3 = sup
θ∈Sd−1c
{γθ (Uθ,i)−m (Xi)}2 − E {γθ (Uθ)−m (X)}2
I4 = sup
θ∈Sd−1c
σ2 (Xi) ε
i − Eσ2 (X)
{γθ (Uθ,i)−m (Xi)}σ (Xi) εi
Bernstein inequality and strong law of large number for α mixing sequence imply that
I3 + I4 = o(1), a.s.. (A.37)
Now (2.13) of Proposition 2.1 provides that
θ∈Sd−1c
u∈[0,1]
|γ̂θ (u)− γθ (u)| = O
n−1/2h−1/2 log n+ h4
, a.s.,
which entail that
I1 = O
n−1/2h−1/2 log n
, a.s., (A.38)
I2 ≤ O
(nh)−1/2 log n+ h4
× sup
θ∈Sd−1c
|γθ (Uθ,i)−m (Xi)− σ (Xi) εi| .
Hence
I2 ≤ O
n−1/2h−1/2 log n+ h4
, a.s.. (A.39)
The lemma now follows from (A.37), (A.38) and (A.39) and Assumption A6.
SINGLE-INDEX PREDICTION MODEL 27
Lemma A.15. Under Assumptions A2 - A6, one has
θ∈Sd−1c
1≤p≤d
R̂ (θ)−R (θ)
− n−1
ξθ,i,p
n−1/2
, a.s., (A.40)
in which
ξθ,i,p = 2 {γθ (Uθ,i)− Yi}
γθ (Uθ,i)−
R (θ) , E (ξθ,i,p) = 0. (A.41)
Furthermore for k = 1, 2
θ∈Sd−1c
R̂ (θ)−R (θ)
n−1/2h−1/2−k log n+ h4−k
, a.s.. (A.42)
Proof. Note that for any p = 1, 2, ..., d
R̂ (θ) = n−1
{γ̂θ (Uθ,i)− Yi}
γ̂θ (Uθ,i) ,
R (θ) = E
{γθ (Uθ)−m (X)}
γθ (Uθ)
{γθ (Uθ)−m (X)− σ (X) ε}
γθ (Uθ)
Thus E (ξθ,i,p) = 2E
{γθ (Uθ,i)− Yi} ∂∂θpγθ (Uθ,i)
R (θ) = 0 and
R̂ (θ)−R (θ)
= (2n)
ξθ,i,p + J1,θ,p + J2,θ,p + J3,θ,p, (A.43)
J1,θ,p = n
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
(γ̂θ − γθ) (Uθ,i) ,
J2,θ,p = n
{γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
(γ̂θ − γθ) (Uθ,i) ,
J3,θ,p = n
{γ̂θ (Uθ,i)− γθ (Uθ,i)}
γθ (Uθ,i) .
Bernstein inequality implies that
θ∈Sd−1c
1≤p≤d
ξθ,i,p
n−1/2 log n
, a.s.. (A.44)
28 LI WANG AND LIJIAN YANG
Meanwhile, applying (2.13) and (2.14) of Proposition 2.1, one obtains that
θ∈Sd−1c
1≤p≤d
|J1,θ,p| = O
log n+ h4
)−1/2
log n+ h3
n−1h−2 log2 n+ h7
, a.s.. (A.45)
Note that
J2,θ,p = n
{γθ (Uθ,i)−m (Xi)− σ (Xi) εi}
(γ̃θ − γθ) (Uθ,i)
−n−1 (E+Eθ)T
{Pθ (E+Eθ)} .
Applying (2.13), one gets
θ∈Sd−1c
1≤p≤d
J2,θ,p + n
−1 (E+Eθ)
{Pθ (E+Eθ)}
, a.s.,
while (A.24), (A.26) and (A.12) entail that with probability 1
θ∈Sd−1c
1≤p≤d
n−1 (E+Eθ)
{Pθ (E+Eθ)}
log n
×N ×N ×O
log n
n−1N log2 n
θ∈Sd−1c
1≤p≤d
|J2,θ,p| = O
h3 + n−1N log2 n
, a.s.. (A.46)
Lastly
J3,θ,p − n−1
(γ̃θ − γθ)
γθ (Uθ,i) = n
−1 (E+Eθ)
BTθ Bθ
By applying (A.24), (A.26), and (A.12), it is clear that with probability 1
θ∈Sd−1c
1≤p≤d
n−1BTθ E+n
BTθ Bθ
log n
×N ×O
h+ (nN)
log n
n−1 log2 n+ (nN)−1/2 log n
while by applying (A.16) of Lemma A.7, one has
θ∈Sd−1c
1≤p≤d
(γ̃θ − γθ)
γθ (Uθ,i)
, a.s.,
SINGLE-INDEX PREDICTION MODEL 29
together, the above entail that
θ∈Sd−1c
1≤p≤d
|J3,θ,p| = O
h4 + n−1 log2 n+ (nN)−1/2 log n
, a.s.. (A.47)
Therefore, (A.43), (A.45), (A.46), (A.47) and Assumption A6 lead to (A.40), which, together
with (A.44), establish (A.42) for k = 1.
Note that the second order derivative of R̂ (θ) and R (θ) with respect to θp, θq are
{γ̂θ (Uθ,i)− Yi}
∂θp∂θq
γ̂θ (Uθ,i) +
γ̂θ (Uθ,i)
γ̂θ (Uθ,i)
E {γθ (Uθ)−m (X)}
∂θp∂θq
γθ (Uθ) + E
γθ (Uθ)
γθ (Uθ)
The proof of (A.42) for k = 2 follows from (2.13), (2.14) and (2.15).
Proof of Proposition 2.2. The result follows from Lemma A.14, Lemma A.15, equations
(A.50) and (A.51).
A.4. Proof of the Theorem 2
Let Ŝ∗p (θ−d) be the p-th element of Ŝ
∗ (θ−d) and for γθ in (2.6), denote
ηi,p := 2
γ̇p − θ0,pθ−10,dγ̇d
(Uθ0,i) {γθ0 (Uθ0,i)− Yi} , (A.48)
where γ̇p is value of
γθ taking at θ = θ0, for any p, q = 1, 2, ..., d − 1.
Lemma A.16. Under Assumptions A2-A6, one has
1≤p≤d−1
Ŝ∗p (θ0,−d)− n−1
n−1/2
, a.s.. (A.49)
Proof. For any p = 1, ..., d − 1
Ŝ∗p (θ−d)− S∗p (θ−d) =
− θpθ−1d
R̂ (θ)−R (θ)
Therefore, according to (A.40), (A.41) and (A.48)
ηi,p = n
ξθ0,i,p − θ0,pθ
ξθ0,i,d, E (ηi,p) = 0,
1≤p≤d−1
Ŝ∗p (θ0,−d)− S∗p (θ0,−d)− n−1
n−1/2
, a.s..
Since S∗ (θ−d) attains its minimum at θ0,−d, for p = 1, ..., d − 1
S∗p (θ0,−d) ≡
− θpθ−1d
R (θ)
which yields (A.49).
30 LI WANG AND LIJIAN YANG
Lemma A.17. The (p, q)-th entry of the Hessian matrix H∗ (θ0,−d) equals lp,q given in Theo-
rem 2.
Proof. It is easy to show that for any p, q = 1, 2, ..., d,
R (θ) =
E {m (X)− γθ (Uθ)}2 = −2E
γθ (Uθ)
γθ (Uθ)
∂θp∂θq
R (θ) = −2E
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θq
γθ (Uθ)
Note that
R∗ (θ−d) =
R (θ)−
R (θ) , (A.50)
∂θp∂θq
R∗ (θ−d) =
∂θp∂θq
R (θ)−
∂θp∂θd
R (θ)−
∂θd∂θq
R (θ)
1− ‖θ−d‖22
R (θ) +
∂θd∂θd
R (θ) . (A.51)
R∗ (θ−d) = −2E
γθ (Uθ)
γθ (Uθ)
+ 2θ−1
γθ (Uθ)
γθ (Uθ)
∂θp∂θq
R∗ (θ−d) = −2E
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θq
γθ (Uθ)
+2θqθ
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θd
γθ (Uθ)
1− ‖θ−d‖22
γθ (Uθ)
γθ (Uθ)
+2θpθ
γθ (Uθ)
γθ (Uθ) + γθ (Uθ)
∂θp∂θq
γθ (Uθ)
−2θpθqθ−2d E
γθ (Uθ)
+ γθ (Uθ)
∂θd∂θd
γθ (Uθ)
Therefore we obtained the desired result.
Proof of Theorem 2. For any p = 1, 2, ..., d − 1, let
fp (t) = Ŝ
tθ̂−d + (1− t) θ0,−d
, t ∈ [0, 1],
fp (t) =
tθ̂−d+(1− t) θ0,−d
θ̂q − θ0,q
SINGLE-INDEX PREDICTION MODEL 31
Note that Ŝ∗ (θ−d) attains its minimum at θ̂−d, i.e., Ŝ
≡ 0. Thus, for any p = 1, 2, ..., d−
1, tp ∈ [0, 1], one has
−Ŝ∗p (θ0,−d) = Ŝ∗p
− Ŝ∗p (θ0,−d) = fp (1)− fp (0)
∂θqθp
tpθ̂−d + (1− tp) θ0,−d
q=1,...,d−1
θ̂−d−θ0,−d
−Ŝ∗ (θ0,−d) =
∂θq∂θp
tpθ̂−d + (1− tp) θ0,−d
p,q=1,...,d−1
θ̂−d − θ0,−d
Now (2.11) of Theorem 1 and Proposition 2.2 with k = 2 imply that uniformly in p, q =
1, 2, ..., d − 1
∂θq∂θp
tpθ̂−d+(1− tp) θ0,−d
−→ lq,p, a.s., (A.52)
where lp,q is given in Theorem 2. Noting that
θ̂−d−θ0,−d
is represented as
∂θq∂θp
tpθ̂−d + (1− tp) θ0,−d
p,q=1,...,d−1
nŜ∗ (θ0,−d) ,
where Ŝ∗ (θ0,−d) =
Ŝ∗p (θ0,−d)
and according to (A.48) and Lemma A.16
Ŝ∗p (θ0,−d) = n
ηp,i + o
n−1/2
, a.s., E (ηp,i) = 0.
Let Ψ (θ0) = (ψpq)
p,q=1
be the covariance matrix of
Ŝ∗p (θ0,−d)
with ψpq given in
Theorem 2. Cramér-Wold device and central limit theorem for α mixing sequences entail that
nŜ∗ (θ0,−d)
d−→ N {0,Ψ(θ0)} .
Let Σ (θ0) = {H∗ (θ0,−d)}−1Ψ(θ0)
{H∗ (θ0,−d)}T
, withH∗ (θ0,−d) being the Hessian matrix
defined in (2.3). The above limiting distribution of
nŜ∗ (θ0,−d), (A.52) and Slutsky’s theorem
imply that
θ̂−d−θ0,−d
d−→ N {0,Σ (θ0)} .
References
Bosq, D. (1998). Nonparametric Statistics for Stochastic Processes. Springer-Verlag, New
York.
32 LI WANG AND LIJIAN YANG
Carroll, R., Fan, J., Gijbles, I. and Wand, M. P. (1997). Generalized partially linear single-
index models. J. Amer. Statist. Assoc. 92 477-489.
Chen, H. (1991). Estimation of a projection -persuit type regression model. Ann. Statist. 19
142-157.
de Boor, C. (2001). A Practical Guide to Splines. Springer-Verlag, New York.
DeVore, R. A. and Lorentz, G. G. (1993). Constructive Approximation: Polynomials and
Splines Approximation. Springer-Verlag, Berlin.
Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Applications. Chapman
and Hall, London.
Friedman, J. H. and Stuetzle, W. (1981). Projection pursuit regression. J. Amer. Statist.
Assoc. 76 817-823.
Härdle, W. (1990). Applied Nonparametric Regression. Cambridge University Press, Cam-
bridge.
Härdle, W. and Hall, P. and Ichimura, H. (1993). Optimal smoothing in single-index models.
Ann. Statist. 21 157-178.
Härdle, W. and Stoker, T. M. (1989). Investigating smooth multiple regression by the method
of average derivatives. J. Amer. Statist. Assoc. 84 986-995.
Hall, P. (1989). On projection pursuit regression. Ann. Statist. 17 573-588.
Hastie, T. J. and Tibshirani, R. J. (1990). Generalized Additive Models. Chapman and Hall,
London.
Horowitz, J. L. and Härdle, W. (1996). Direct semiparametric estimation of single-index
models with discrete covariates. J. Amer. Statist. Assoc. 91 1632-1640.
Hristache, M., Juditski, A. and Spokoiny, V. (2001). Direct estimation of the index coefficients
in a single-index model. Ann. Statist. 29 595-623.
Huang, J. Z. (2003). Local asymptotics for polynomial spline regression. Ann. Statist. 31
1600-1635.
Huang, J. and Yang, L. (2004). Identification of nonlinear additive autoregressive models. J.
R. Stat. Soc. Ser. B Stat. Methodol. 66 463-477.
SINGLE-INDEX PREDICTION MODEL 33
Huber, P. J. (1985). Projection pursuit (with discussion). Ann. Statist. 13 435-525.
Ichimura, H. (1993). Semiparametric least squares (SLS) and weighted SLS estimation of
single-index models Journal of Econometrics 58 71-120.
Klein, R. W. and Spady. R. H. (1993). An efficient semiparametric estimator for binary
response models. Econometrica 61 387-421.
Mammen, E., Linton, O. and Nielsen, J. (1999). The existence and asymptotic properties of
a backfitting projection algorithm under weak conditions. Ann. Statist. 27 1443-1490.
Pham, D. T. (1986). The mixing properties of bilinear and generalized random coefficient
autoregressive models. Stochastic Anal. Appl. 23 291-300.
Powell, J. L., Stock, J. H. and Stoker, T. M. (1989). Semiparametric estimation of index
coefficients. Econometrica. 57 1403-1430.
Tong, H. (1990) Nonlinear Time Series: A Dynamical System Approach. Oxford, U.K.:
Oxford University Press.
Tong, H., Thanoon, B. and Gudmundsson, G. (1985) Threshold time series modeling of two
icelandic riverflow systems. Time Series Analysis in Water Resources. ed. K. W. Hipel,
American Water Research Association.
Wang, L. and Yang, L. (2007). Spline-backfitted kernel smoothing of nonlinear additive
autoregression model. Ann. Statist. Forthcoming.
Xia, Y. and Li, W. K. (1999). On single-index coefficient regression models. J. Amer. Statist.
Assoc. 94 1275-1285.
Xia, Y., Li, W. K., Tong, H. and Zhang, D. (2004). A goodness-of-fit test for single-index
models. Statist. Sinica. 14 1-39.
Xia, Y., Tong, H., Li, W. K. and Zhu, L. (2002). An adaptive estimation of dimension
reduction space. J. R. Stat. Soc. Ser. B Stat. Methodol. 64 363-410.
Xue, L. and Yang, L. (2006 a). Estimation of semiparametric additive coefficient model. J.
Statist. Plann. Inference 136, 2506-2534.
Xue, L. and Yang, L. (2006 b). Additive coefficient modeling via polynomial spline. Statistica
Sinica 16 1423-1446.
Table 2: Report of Example 2
Sample Size n Dimension d
Average MSE Time
MAVE SIP MAVE SIP
4 0.00020 0.00018 1.91 0.19
10 0.00031 0.00043 2.17 0.10
30 0.00106 0.00285 2.77 0.13
50 0.00031 0.00043 3.29 0.10
100 0.00681 0.00620 5.94 0.31
200 0.00529 0.00407 27.90 0.49
4 0.00008 0.00008 3.28 0.09
10 0.00012 0.00017 3.93 0.13
30 0.00017 0.00058 5.41 0.15
50 0.00032 0.00127 8.48 0.16
100 — 0.00395 — 0.44
200 — 0.00324 — 0.73
4 0.00004 0.00003 5.32 0.17
10 0.00005 0.00007 7.49 0.24
30 0.00006 0.00017 10.08 0.26
50 0.00007 0.00030 15.42 0.24
100 0.00015 0.00061 40.81 0.54
200 — 0.00197 — 1.44
4 0.00002 0.00001 14.44 0.76
10 0.00002 0.00003 24.54 0.79
30 0.00002 0.00008 32.51 0.83
50 0.00002 0.00010 52.93 0.89
100 0.00003 0.00012 143.07 0.99
200 0.00004 0.00020 386.80 1.96
400 — 0.00054 — 4.98
4 0.00001 0.00001 33.57 1.95
10 0.00001 0.00001 62.54 3.64
30 0.00001 0.00002 92.41 1.95
50 0.00001 0.00003 155.38 2.72
100 0.00001 0.00005 275.73 1.81
200 0.00008 0.00006 2432.56 2.84
400 — 0.00010 — 9.35
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
−1 0 1 2
−1 0 1 2
Figure 1: Example 1. (a) and (b) Plots of the actual surface m in model (4.1) with re-
spect to δ = 0, 1; (c) and (d) Plots of various univariate functions with respect to δ = 0, 1:
XTi θ̂, Yi
, 1 ≤ i ≤ 50 (dots); the univariate function g (solid line); the estimated function of
g by plugging in the true index coefficient θ0 (dotted line); the estimated function of g by
plugging in the estimated index coefficient (dashed line) θ̂ = (0.69016, 0.72365)T for δ = 0 and
(0.72186, 0.69204)T for δ = 1.
−2 −1 0 1 2
n= 100 ,  d= 10
−2 −1 0 1 2
n= 100 ,  d= 50
−2 −1 0 1 2
n= 100 ,  d= 100
−3 −2 −1 0 1 2 3
n= 100 ,  d= 200
Figure 2: Example 2. Plots of the spline estimator of g with the estimated index parameter θ̂
(dotted curve), cubic spline estimator of g with the true index parameter θ0 (dashed curves),
the true function m (x) in (4.2) (solid curve), and the data scatter plots (dots).
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=10
n=100
n=200
n=500
n=1000
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=50
n=100
n=200
n=500
n=1000
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=100
n=100
n=200
n=500
n=1000
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Density Estimation, d=200
n=100
n=200
n=500
n=1000
Figure 3: Example 2. Kernel density estimators of the 100 ‖θ̂ − θ0‖/
0 200 400 600 800 1000
0 200 400 600 800 1000
0 200 400 600 800 1000
Figure 4: Time plots of the daily Jökulsá Eystri River data (a) river flow Yt (solid line) with its
trend (dashed line) (b) temperature Xt (solid line) with its trend (dashed line) (c) precipitation
Zt (solid line) with its trend (dashed line).
+++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++
++++++
+++++
+++++
++++++++++
++++++++++
+++++++
+++++++
++++++++++++
+++++
++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++
++++++++++++++++++
+++++
+++++
++++++
+++++++
+++++
+++++
+++++
+++++
++++++++
+++++
++++++++++++++
+++++++
+++++++++++++++++++++++++
++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++
++++++++++
+++++++++++++++++++++
++++++
++++++
++++++
++++++++++
+++++++++
+++++
++++++++++++++++++++++++++++++
+++++++
++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0 200 400 600 800 1000
0 200 400 600 800 1000
800 900 1000 1100
++++++++++++++++++++++++++++++++++++++++++++++++++
++++++
++++++++++
+++++++++++++++++++++
++++++
++++++
++++++
++++++++++
+++++++++
+++++
++++++++++++++++++++++++++++++
+++++++
++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Figure 5: (a) The scatter plot of the river flow (“+”) and the fitted plot of the river flow (line)
and (b) Residuals of the fitted SIP model (c) Out-of-sample rolling forecasts (line) of the river
flow for the entire third year (“+”) based on the first two years’ river flow.
