Spline Single-Index Prediction Model

TL;DR

This paper introduces a robust single-index prediction model using spline estimators, effective for high-dimensional, weakly dependent data, with proven consistency, normality, and fast computation, demonstrated through simulations and real data application.

Contribution

It proposes a spline-based single-index prediction method that is robust, computationally efficient, and applicable to high-dimensional, dependent data, extending existing models.

Findings

Spline estimator is root-n consistent and asymptotically normal.

The iterative routine is fast enough for large, high-dimensional data.

Application to Icelandic river flow data shows superior forecasting performance.

Abstract

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.

Figures

Equations

$$Y_i=m\left({\mathbf{X}}_i\right)+\sigma \left({\mathbf{X}}_i\right){\epsilon }_i,m\left({\mathbf{X}}_i\right)=E\left(Y_i|{\mathbf{X}}_i\right),$$

$$g\left(\nu \right)=E\left[m\left(\mathbf{X}\right)|{\mathbf{X}}^T{\theta }_0=\nu \right],$$

$$m_{\theta }\left(X_{\theta }\right)=E\left(Y|X_{\theta }\right)=E\left\{m\left(\mathbf{X}\right)|X_{\theta }\right\}.$$

$$R\left(\theta \right)=E\left[{\left\{Y-m_{\theta }\left(X_{\theta }\right)\right\}}^2\right]=E{\left\{m\left(\mathbf{X}\right)-m_{\theta }\left(X_{\theta }\right)\right\}}^2+E{\sigma }^2\left(\mathbf{X}\right),$$

$${\theta }_0=\arg \underset{\theta \in S_{+}^{d-1}}{\min }R\left(\theta \right).$$

$$S_c^{d-1}=\left\{\left({\theta }_1,...,{\theta }_d\right)|\sum _{p=1}^d{\theta }_p^2=1,{\theta }_d\ge \sqrt{1-c^2}\right\},c\in \left(0,1\right)$$

$$R^{\ast }\left({\theta }_{-d}\right)=R\left({\theta }_1,{\theta }_2,...,{\theta }_{d-1},\sqrt{1-{\Vert {\theta }_{-d}\Vert }_2^2}\right),$$

$$S^{\ast }\left({\theta }_{-d}\right)=\frac{\partial }{\partial {\theta }_{-d}}{R}^{\ast }\left({\theta }_{-d}\right)\text{, }{H}^{\ast }\left({\theta }_{-d}\right)=\frac{{\partial }^2}{\partial {\theta }_{-d}\partial {\theta }_{-d}^T}{R}^{\ast }\left({\theta }_{-d}\right).$$

$$C^{(k)}\left(B_a^d\right)=\left\{m|\text{the }k\text{th order partial derivatives of }m\text{ are continuous on }{B}_a^d\right\}$$

$$\left\{\begin{array}[]{ll}c_{f}/\text{Vol}_{d}\left(B_{a}^{d}\right)\leq f\left(\mathbf{x}\right)\leq C_{f}/\text{Vol}_{d}\left(B_{a}^{d}\right),&\mathbf{x}\in B_{a}^{d}\\ f\left(\mathbf{x}\right)\equiv 0,&\mathbf{x}\notin B_{a}^{d}\end{array}\right..$$

$$U_{\theta }=F_d\left(X_{\theta }\right),U_{\theta ,i}=F_d\left(X_{\theta ,i}\right),1\le i\le n,$$

$$F_d\left(\nu \right)={\int }_{-1}^{\nu /a}\frac{\Gamma \left(d+1\right)}{\Gamma {\left\{\left(d+1\right)/2\right\}}^2{2}^d}{\left(1-t^2\right)}^{\left(d-1\right)/2}dt,\nu \in \left[-a,a\right].$$

$$f_{\theta }\left(u\right)=\left\{F_d^{{}^{\prime }}\left(v\right)\right\}{f}_{X_{\theta }}\left(v\right),v=F_d^{-1}\left(u\right),$$

$$P\left(\nu \le X_{\theta }\le \nu +△\nu \right)=P\left(\mathbf{X}\in {\mathcal{D}}_{\nu }\right)={\int }_{{\mathcal{D}}_{\nu }}f\left(\mathbf{x}\right)d\mathbf{x}.$$

$$\frac{c_f{\text{Vol}}_d({\mathcal{D}}_{\nu })}{{\text{Vol}}_d\left(B_a^d\right)}\le P\left(\nu \le X_{\theta }\le \nu +△\nu \right)\le \frac{C_f{\text{Vol}}_d({\mathcal{D}}_{\nu })}{{\text{Vol}}_d\left(B_a^d\right)}.$$

$${\text{Vol}}_d({\mathcal{D}}_{\nu })={\text{Vol}}_{d-1}({\mathcal{J}}_{\nu })△\nu +o\left(△\nu \right),$$

$${\text{Vol}}_{d-1}\left({\mathcal{J}}_{\nu }\right)={\pi }^{\left(d-1\right)/2}{\left(a^2-{\nu }^2\right)}^{\left(d-1\right)/2}/\Gamma \left\{(d+1)/2\right\},$$

$$\frac{{\text{Vol}}_{d-1}({\mathcal{J}}_{\nu })}{{\text{Vol}}_d\left(B_a^d\right)}=\frac{1}{a\sqrt{\pi }}\frac{\Gamma \left(d+1\right)}{{\left\{\Gamma \left(\frac{d+1}{2}\right)\right\}}^2{2}^d}{\left\{1-{\left(\frac{\nu }{a}\right)}^2\right\}}^{\left(d-1\right)/2}.$$

$${\gamma }_{\theta }\left(U_{\theta }\right)=E\left\{m\left(\mathbf{X}\right)|U_{\theta }\right\}=E\left\{m\left(\mathbf{X}\right)|X_{\theta }\right\}=m_{\theta }\left(X_{\theta }\right),$$

$$R\left(\theta \right)=E\left[{\left\{Y-{\gamma }_{\theta }\left(U_{\theta }\right)\right\}}^2\right]=E{\left\{m\left(\mathbf{X}\right)-{\gamma }_{\theta }\left(U_{\theta }\right)\right\}}^2+E{\sigma }^2\left(\mathbf{X}\right).$$

$$t_{1-k}=...=t_{-1}=t_0=0<t_1<...<t_N<1=t_{N+1}=...=t_{N+k},$$

$${\hat{\gamma }}_{\theta }\left(\cdot \right)=\arg \underset{\gamma \left(\cdot \right)\in {\Gamma }^{\left(2\right)}\left[0,1\right]}{\min }\sum _{i=1}^n{\left\{Y_i-\gamma \left(U_{\theta ,i}\right)\right\}}^2,\text{mbox}{\hat{m}}_{\theta }\left(\nu \right)={\hat{\gamma }}_{\theta }\left\{F_d\left(\nu \right)\right\}.$$

$$\hat{R}\left(\theta \right)=n^{-1}\sum _{i=1}^n{\left\{Y_i-{\hat{\gamma }}_{\theta }\left(U_{\theta ,i}\right)\right\}}^2=n^{-1}\sum _{i=1}^n{\left\{Y_i-{\hat{m}}_{\theta }\left(X_{\theta ,i}\right)\right\}}^2,$$

$$\hat{\theta }=\arg \underset{\theta \in S_c^{d-1}}{\min }\hat{R}\left(\theta \right),$$

$$\hat{g}\left(\nu \right)=\left\{\arg \underset{\gamma \left(\cdot \right)\in {\Gamma }^{\left(2\right)}\left[0,1\right]}{\min }\sum _{i=1}^n{\left\{Y_i-\gamma \left(U_{\hat{\theta },i}\right)\right\}}^2\right\}\left\{F_d\left(\nu \right)\right\}.$$

$$0<c_{\sigma }\le \underset{\mathbf{x}\in B_a^d}{\inf }\sigma \left(\mathbf{x}\right)\le \underset{\mathbf{x}\in B_a^d}{\sup }\sigma \left(\mathbf{x}\right)\le C_{\sigma }<\infty .$$

$$\alpha \left(k\right)=\underset{B\in \sigma \left\{{\mathbf{Z}}_s,s\le t\right\},C\in \sigma \left\{{\mathbf{Z}}_s,s\ge t+k\right\}}{\sup }\left|P\left(B\cap C\right)-P\left(B\right)P\left(C\right)\right|,\text{mbox}k\ge 1.$$

$${\hat{\theta }}_{-d}{⟶\theta }_{0,-d},a.s..$$

$$\underset{{\Vert {\theta }_{-d}\Vert }_2\le \sqrt{1-c^2}}{\sup }\left|{\hat{R}}^{\ast }\left({\theta }_{-d}\right)-R^{\ast }\left({\theta }_{-d}\right)\right|⟶0,a.s..$$

$$R^{\ast }\left({\hat{\theta }}_{-d}\left(\omega \right),\omega \right)-R^{\ast }\left({\theta }_{0,-d}\right)<\delta /2+{\hat{R}}^{\ast }\left({\hat{\theta }}_{-d}\left(\omega \right),\omega \right)-R^{\ast }\left({\theta }_{0,-d}\right)<\delta /2+\delta /2=\delta .$$

Full text

[TABLE]

Proof. Denote by $\left(\Omega,\mathcal{F},\mathcal{P}\right)$ the probability space on which all $\left\{\left(\mathbf{X}_{i}^{T},Y_{i}\right)\right\}_{i=1}^{\infty}$ are defined. By Proposition 2.2, given at the end of this section

[TABLE]

So for any $\delta>0$ and $\omega\in\Omega$, there exists an integer $n_{0}\left(\omega\right)$, such that when $n>n_{0}\left(\omega\right)$, $\hat{R}^{*}\left(\mathbf{\theta}_{0,-d},\omega\right)-R^{*}\left(\mathbf{\theta}_{0,-d}\right)<$ $\delta/2$. Note that $\hat{\mathbf{\theta}}_{-d}=\hat{\mathbf{\theta}}_{-d}\left(\omega\right)$ is the minimizer of $\hat{R}^{*}\left(\mathbf{\theta}_{-d},\omega\right)$, so $\hat{R}^{*}\left(\hat{\mathbf{\theta}}_{-d}\left(\omega\right),\omega\right)-R^{*}\left(\mathbf{\theta}_{0,-d}\right)<\delta/2$. Using (2.12), there exists $n_{1}\left(\omega\right)$, such that when $n>n_{1}\left(\omega\right)$, $R^{*}\left(\hat{\mathbf{\theta}}_{-d}\left(\omega\right),\omega\right)-\hat{R}^{*}\left(\hat{\mathbf{\theta}}_{-d}\left(\omega\right),\omega\right)<$ $\delta/2$. Thus, when $n>\max\left(n_{0}\left(\omega\right),n_{1}\left(\omega\right)\right)$,

[TABLE]

According to Assumption A1, $R^{*}$ is locally convex at $\mathbf{\theta}_{0,-d}$, so for any $\varepsilon>0$ and any $\omega$, if $R^{*}\left(\hat{\mathbf{\theta}}_{-d}\left(\omega\right),\omega\right)-R^{*}\left(\mathbf{\theta}_{0,-d}\right)<\delta$, then $\left\|\hat{\mathbf{\theta}}_{-d}\left(\omega\right)\mathbf{-\theta}_{0,-d}\right\|<\varepsilon$ for $n$ large enough , which implies the strong consistency.

[TABLE]

where $\Sigma\left(\mathbf{\theta}_{0}\right)=\left\{H^{*}\left(\mathbf{\theta}_{0,-d}\right)\right\}^{-1}\Psi\left(\mathbf{\theta}_{0}\right)\left\{H^{*}\left(\mathbf{\theta}_{0,-d}\right)\right\}^{-1}$, $H^{*}\left(\mathbf{\theta}_{0,-d}\right)=\left\{l_{pq}\right\}_{p,q=1}^{d-1}$ and $\Psi\left(\mathbf{\theta}_{0}\right)=\left\{\psi_{pq}\right\}_{p,q=1}^{d-1}$ with

[TABLE]

[TABLE]

in which $\dot{\gamma}_{p}$ and $\ddot{\gamma}_{p,q}$ are the values of $\frac{\partial}{\partial\theta_{p}}\gamma_{\mathbf{\theta}}$, $\frac{\partial^{2}}{\partial\theta_{p}\partial\theta_{q}}\gamma_{\mathbf{\theta}}$ taking at $\mathbf{\theta=\theta}_{0}$, for any $p,q=1,2,...,d-1$ and $\gamma_{\mathbf{\theta}}$ is given in (2.6).

Remark 2.4. Consider the Generalized Linear Model (GLM): $Y=g\left(\mathbf{X}^{T}\mathbf{\theta}_{0}\right)+\sigma\left(\mathbf{X}\right)\varepsilon$, where $g$ is a known link function. Let $\tilde{\mathbf{\theta}}$ be the nonlinear least squared estimator of $\mathbf{\theta}_{0}$ in GLM. Theorem 2 shows that under the assumptions A1-A6, the asymptotic distribution of the $\hat{\mathbf{\theta}}_{-d}$ is the same as that of $\tilde{\mathbf{\theta}}$. This implies that our proposed SIP estimator $\hat{\mathbf{\theta}}_{-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. Proposition 2.1 establishes the uniform convergence rate of the derivatives of $\hat{\gamma}_{\mathbf{\theta}}$ up to order 2 to those of $\gamma_{\mathbf{\theta}}$ in $\mathbf{\theta}$. Proposition 2.2 shows that the derivatives of the risk function up to order 2 are uniformly almost surely approximated by their empirical versions.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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 $\mathbf{\theta}_{0}$ and $g$. We first introduce some new notation. For fixed $\mathbf{\theta}$, write the B-spline matrix as $\mathbf{B}_{\mathbf{\theta}}=\left\{B_{j,4}\left(U_{\mathbf{\theta},i}\right)\right\}_{i=1,j=-3}^{n,\text{ }N}$ and

[TABLE]

as the projection matrix onto the cubic spline space $\Gamma_{n,\mathbf{\theta}}^{\left(2\right)}$. For any $p=1,...,d$, denote

[TABLE]

as the first order partial derivatives of $\mathbf{B}_{\mathbf{\theta}}$ and $\mathbf{P}_{\mathbf{\theta}}$ with respect to $\mathbf{\theta}$.

Let $\hat{S}^{*}(\mathbf{\theta}_{-d})$ be the score vector of $\hat{R}^{*}\left(\mathbf{\theta}_{-d}\right)$, i.e.

[TABLE]

The next lemma provides the exact forms of $\hat{S}^{*}(\mathbf{\theta}_{-d})$.

[TABLE]

where for any $p=1,2,...,d$

[TABLE]

where $\mathbf{\dot{B}}_{p}=\left\{\left\{B_{j,3}\left(U_{\mathbf{\theta},i}\right)-B_{j+1,3}\left(U_{\mathbf{\theta},i}\right)\right\}\dot{F}_{d}\left(\mathbf{X}_{\mathbf{\theta},i}\right)h^{-1}X_{i,p}\right\}_{i=1,j=-3}^{n,\text{ }N}$ with

[TABLE]

Proof. For any $p=1,2,...,d$, the derivatives of B-splines in de Boor (2001) implies

[TABLE]

Next, note that

[TABLE]

Since

[TABLE]

and $\frac{\partial}{\partial\theta_{p}}\left(\mathbf{B}_{\mathbf{\theta}}^{T}\mathbf{B}_{\mathbf{\theta}}\right)=\mathbf{\dot{B}}_{p}^{T}\mathbf{B}_{\mathbf{\theta}}+\mathbf{B}_{\mathbf{\theta}}^{T}\mathbf{\dot{B}}_{p}$, thus

[TABLE]

Hence

[TABLE]

Thus, (3.4) follows immediately.

In practice, the estimation is implemented via the following procedure.

Step 1. Standardize the predictor vectors $\left\{\mathbf{X}_{i}\right\}_{i=1}^{n}$ and for each fixed $\mathbf{\theta}\in S_{c}^{d-1}$ obtain the CDF transformed variables $\left\{U_{\mathbf{\theta},i}\right\}_{i=1}^{n}$ of the SIP variable $\left\{X_{\mathbf{\theta},i}\right\}_{i=1}^{n}$ through formula (2.5), where the radius $a$ is taken to be the 95% percentile of $\left\{\|\mathbf{X}_{i}\|\right\}_{i=1}^{n}$.

Step 2. Compute quadratic and cubic B-spline basis at each value $U_{\mathbf{\theta},i}$, where the number of interior knots $N$ is

[TABLE]

Step 3. Find the estimator $\hat{\mathbf{\theta}}$ of $\mathbf{\theta}_{0}$ by minimizing $\hat{R}^{*}$ through the port optimization routine with $\left(0,0,...,1\right)^{T}$ as the initial value and the empirical score vector $\hat{S}^{*}$ in (3.3). If $d<n$, one can take the simple LSE (without the intercept) for data $\left\{Y_{i},\mathbf{X}_{i}\right\}_{i=1}^{n}$ with its last coordinate set positive.

Step 4. Obtain the spline estimator $\hat{g}$ of $g$ by plugging $\hat{\mathbf{\theta}}$ obtained in Step 3 into (2.10).

Remark 3.1. In (3.5), $c_{1}$ and $c_{2}$ are positive integers and $[\nu]$ denotes the integer part of $\nu$. The choice of the tuning parameter $c_{1}$ 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 $c_{1}=1$ to save computing for massive data sets. The first term ensures Assumption A6. The addition constrain $c_{2}$ can be taken from 5 to 10 for smooth monotonic or smooth unimodel regression and $c_{2}>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 $c_{1}=1,c_{2}=5$. All of our codes have been written in R.

Example 1. Consider the model in Xia, Li, Tong and Zhang (2004)

[TABLE]

where $\mathbf{X}=\left(X_{1},X_{2}\right)^{T}{\sim}N(\mathbf{0},I_{2})$, truncated by $[-2.5,2.5]^{2}$ and

[TABLE]

If $\delta=0$, then the underlying true function $m$ is a single-index function, i.e., $m\left(\mathbf{X}\right)=\sqrt{2}\mathbf{X}^{T}\mathbf{\theta}_{0}+4\exp\left\{-2\left(\mathbf{X}^{T}\mathbf{\theta}_{0}\right)^{2}\right\}$, where $\mathbf{\theta}_{0}^{T}=\left(1,1\right)/\sqrt{2}$. While $\delta\neq 0$, then $m$ is not a genuine single-index function. An impression of the bivariate function $m$ for $\delta=0$ and $\delta=1$ can be gained in Figure 1 (a) and (b), respectively.

For $\delta=0,1$, we draw $100$ random realizations of each sample size $n=50,100,300$ respectively. To demonstrate how close our SIP estimator is to the true index parameter $\mathbf{\theta}_{0}$, Table 1 lists the sample mean (MEAN), bias (BIAS), standard deviation (SD), the mean squared error (MSE) of the estimates of $\mathbf{\theta}_{0}$ and the average MSE of both directions. From this table, we find that the SIP estimators are very accurate for both cases $\delta=0$ and $\delta=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

[TABLE]

in which $\mathbf{X}_{i}=\left\{X_{i,1},...,X_{i,d}\right\}^{T}$ and $\varepsilon_{i}$, $i=1,...,n$, are $\stackrel{{\scriptstyle i.i.d}}{{\sim}}$ $N\left(0,1\right)$, $\sigma_{0}=0.2$. In our simulation, the true parameter $\mathbf{\theta}_{0}^{T}=\left.\left(1,1,0,...,0,1\right)\right/\sqrt{3}$ for different sample size $n$ and dimension $d$. The 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\leq 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 $\hat{\mathbf{\theta}}_{-d}$ is strongly consistent of $\mathbf{\theta}_{0,-d}$. To see the convergence, we run $100$ replications and in each replication, the value of $\|\hat{\mathbf{\theta}}-\mathbf{\theta}_{0}\|/\sqrt{d}$ is computed. Figure 3 plots the kernel density estimations of the $100$ $\|\hat{\mathbf{\theta}}-\mathbf{\theta}_{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 [math], 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 Jökulsá 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 ($Y_{t}$), measured in meter cubed per second of Jökulsá Eystri River. The exogenous variables are temperature ($X_{t}$) in degrees Celsius and daily precipitation ($Z_{t}$) 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 (NAAR$X$) 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 $X_{t}$, $Y_{t}$, $Z_{t}$ to denote the detrended series.

In the analysis, we pre-select all the lagged values in the last 7 days (1 week), i.e., the predictor pool is $\left\{Y_{t-1},...,Y_{t-7},X_{t},X_{t-1},...,X_{t-7},Z_{t},Z_{t-1},...,Z_{t-7},\right\}$. 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 $\left\{Y_{t-1},...,Y_{t-4},X_{t},X_{t-1},X_{t-2},Z_{t},Z_{t-1}\right\}$. Based on this selection, we fit the SIP model again and obtain the estimate of the SIP coefficient $\hat{\mathbf{\theta}}=\left\{-0.877,0.382,-0.208,0.125,-0.046,-0.034,0.004,-0.126,0.079\right\}^{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 NAAR$X$, $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 dependence 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 $\sqrt{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.

[TABLE]

where $\mathbf{\alpha:=}\left(\alpha_{-1,2},\alpha_{0,2},...,\alpha_{N,2},...,\alpha_{N,4}\right)$. In particular, under Assumption A2, for any fixed $\mathbf{\theta}$

[TABLE]

Proof. It follows from the B-spline property on page 96 of de Boor (2001), $\sum_{k=2}^{4}\sum_{j=-k+1}^{N}B_{j,k}\equiv 3$ on $\left[0,1\right]$. So the right inequality follows immediate for $r=\infty$. When $1\leq r<\infty$, we use Hölder’s inequality to find

[TABLE]

Since all the knots are equally spaced, $\int_{-\infty}^{\infty}B_{j,k}\left(u\right)du\leq h$, the right inequality follows from

[TABLE]

When $r<1$, we have

[TABLE]

Since $\int_{-\infty}^{\infty}B_{j,k}^{r}\left(u\right)du\leq t_{j+k}-t_{j}=kh$ and

[TABLE]

the right inequality follows in this case as well. For the left inequalities, we derive from Theorem 5.4.2, DeVore and Lorentz (1993)

[TABLE]

for any $0<r\leq\infty,$ so

[TABLE]

Since each $u\in\left[0,1\right]$ appears in at most $k$ intervals $\left(t_{j,}t_{j+k}\right)$, adding up these inequalities, we obtain that

[TABLE]

The left inequality follows.

For any functions $\phi$ and $\varphi$, define the empirical inner product and the empirical norm as

[TABLE]

In addition, if functions $\phi,\varphi$ are $L_{2}\left[0,1\right]$-integrable, define the theoretical inner product and its corresponding theoretical $L_{2}$ norm as

[TABLE]

[TABLE]

Proof. We only prove the case $k=k^{\prime}=4$, all other cases are similar. Let

[TABLE]

with the second moment

[TABLE]

where $\left\{EB_{j,4}\left(U_{\mathbf{\theta},i}\right)B_{j^{\prime},4}\left(U_{\mathbf{\theta},i}\right)\right\}^{2}\sim N^{-2}$, $E\left[B_{j,4}^{2}\left(U_{\mathbf{\theta},i}\right)B_{j^{\prime},4}^{2}\left(U_{\mathbf{\theta},i}\right)\right]\sim N^{-1}$ by Assumption A2. Hence, $E\zeta_{\mathbf{\theta,}j,j^{\prime},i}^{2}\sim N^{-1}$. The $k$-th moment is given by

[TABLE]

where $\left|EB_{j,4}\left(U_{\mathbf{\theta},i}\right)B_{j^{\prime},4}\left(U_{\mathbf{\theta},i}\right)\right|^{k}\sim N^{-k}$, $E\left|EB_{j,4}\left(U_{\mathbf{\theta},i}\right)B_{j^{\prime},4}\left(U_{\mathbf{\theta},i}\right)\right|^{k}\sim N^{-1}$. Thus, there exists a constant $C>0$ such that $E\left|\zeta_{\mathbf{\theta,}j,j^{\prime},i}\right|^{k}\leq C2^{k-1}k!E\zeta_{j,j^{\prime},i}^{2}$. So the Cramér’s condition is satisfied with Cramér’s constant $c^{*}$. By the Bernstein’s inequality (see Bosq (1998), Theorem 1.4, page 31), we have for $k=3$

[TABLE]

where

[TABLE]

[TABLE]

Observe that $5c\delta_{n}=o(1)$ by Assumption A6, then by taking q\such that $\left[\frac{n}{q+1}\right]\geq c_{0}\log n$, $q\geq c_{1}n/\log n$ for some constants $c_{0},c_{1}$, one has $a_{1}=O(n/q)=O\left(\log n\right)$, $a_{2}\left(3\right)=o\left(n^{2}\right)$ via Assumption A6 again. Assumption A5 yields that

[TABLE]

Thus, for fixed $\mathbf{\theta}\in S_{c}^{d-1}$, when $n$ large enough

[TABLE]

We divide each range of $\theta_{p}$, $p=1,2,...,d-1$, into n^{6/(d-1)}\ equally spaced intervals with disjoint endpoints $-1=\theta_{p,0}<\theta_{p,1}<...<\theta_{p,M_{n}}=1$, for $p=1,...,d-1$. Projecting these small cylinders onto $S_{c}^{d-1}$, the radius of each patch $\Lambda_{r}$, $r=1,...,M_{n}$ is bounded by $cM_{n}^{-1}$. Denote the projection of the $M_{n}$ points as $\mathbf{\theta}_{r}=\left(\mathbf{\theta}_{r,-d},\sqrt{1-\left\|\mathbf{\theta}_{r,-d}\right\|_{2}^{2}}\right)$, $r=0,1,...,M_{n}$. Employing the discretization method, $\sup\limits_{\mathbf{\theta}\in S_{c}^{d-1}}\max\limits_{\begin{subarray}{c}1\leq j,j^{\prime}\leq N\end{subarray}}\left|\zeta_{\mathbf{\theta},j,j^{\prime},i}\right|$ is bounded by

[TABLE]

By (A.1) and Assumption A6, there exists large enough value $\delta>0$ such that

[TABLE]

which implies that

[TABLE]

Thus, Borel-Cantelli Lemma entails that

[TABLE]

Employing Lipschitz continuity of the cubic B-spline, one has with probability 1

[TABLE]

Therefore Assumption A2, (A.2), (A.3) and (A.4) lead to the desired result.

Denote by $\Gamma=\Gamma^{\left(0\right)}\cup\Gamma^{(1)}\cup\Gamma^{(2)}$ the space of all linear, quadratic and cubic spline functions on $\left[0,1\right]$. We establish the uniform rate at which the empirical inner product approximates the theoretical inner product for all B-splines $B_{j,k}$ with $k=2,3,4$.

[TABLE]

Proof. Denote without loss of generality,

[TABLE]

for any two $3\left(N+3\right)$-vectors

[TABLE]

Then for fixed $\mathbf{\theta}$

[TABLE]

[TABLE]

According to Lemma A.1, one has for any $\mathbf{\theta}\in S_{c}^{d-1}$,

[TABLE]

[TABLE]

Hence

[TABLE]

[TABLE]

which, together with Lemma A.2, imply (A.5).

A.2. Proof of Proposition 2.1

For any fixed $\mathbf{\theta}$, we write the response $\mathbf{Y}^{T}=\left(Y_{1},...,Y_{n}\right)$ as the sum of a signal vector $\gamma$θ, a parametric noise vector $\mathbf{E}_{\mathbf{\theta}}$ and a systematic noise vector $\mathbf{E}$, i.e.,

[TABLE]

in which the vectors \gamma$${}_{\mathbf{\theta}}^{T}=\left\{\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta,}1}\right),...,\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta,}n}\right)\right\}, $\mathbf{E}^{T}=\left\{\sigma\left(\mathbf{X}_{1}\right)\varepsilon_{1},...,\sigma\left(\mathbf{X}_{n}\right)\varepsilon_{n}\right\}$ and $\mathbf{E}_{\mathbf{\theta}}^{T}=\left\{m\left(\mathbf{X}_{1}\right)-\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta},1}\right),...,m\left(\mathbf{X}_{n}\right)-\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta},n}\right)\right\}$.

Remark A.1. If $m$ is a genuine single-index function, then $\mathbf{E}_{\mathbf{\theta}_{0}}\equiv 0$, thus the proposed SIP model is exactly the single-index model.

Let $\Gamma_{n,\mathbf{\ \theta}}^{\left(2\right)}$ be the cubic spline space spanned by $\left\{\mathbf{B}_{j,4}\left(U_{\mathbf{\theta},i}\right)\right\}_{i=1}^{n}$, $-3\leq j\leq N$ for fixed $\mathbf{\theta}$. Projecting $\mathbf{Y}$ onto $\Gamma_{n,\mathbf{\ \theta}}^{\left(2\right)}$ yields that

[TABLE]

where $\hat{\gamma}_{\mathbf{\theta}}$ is given in (2.8). We break the cubic spline estimation error $\hat{\gamma}_{\mathbf{\theta}}\left(u_{\mathbf{\theta}}\right)-\gamma_{\mathbf{\theta}}\left(u_{\mathbf{\theta}}\right)$ into a bias term $\tilde{\gamma}_{\mathbf{\theta}}\left(u_{\mathbf{\theta}}\right)-\gamma_{\mathbf{\theta}}\left(u_{\mathbf{\theta}}\right)$ and two noise terms $\tilde{\varepsilon}_{\mathbf{\theta}}\left(u_{\mathbf{\theta}}\right)$ and $\hat{\varepsilon}_{\mathbf{\theta}}\left(u_{\mathbf{\theta}}\right)$

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

In the above, we denote by $\mathbf{V}_{n,\mathbf{\theta}}$ the empirical inner product matrix of the cubic B-spline basis and similarly, the theoretical inner product matrix as $\mathbf{V}_{\mathbf{\theta}}$

[TABLE]

In Lemma A.5, we provide the uniform upper bound of $\left\|\mathbf{V}_{n,\mathbf{\theta}}^{-1}\right\|_{\infty}$ and $\left\|\mathbf{V}_{\mathbf{\theta}}^{-1}\right\|_{\infty}$. Before that, we first describe a special case of Theorem 13.4.3 in DeVore and Lorentz (1993).

[TABLE]

with matrices $\mathbf{V}_{\mathbf{\theta}}$ and $\mathbf{V}_{n,\mathbf{\theta}}$ defined in (A.10). In addition, there exists a constant $C>0$ such that

[TABLE]

Proof. First we compute the lower and upper bounds for the eigenvalues of $\mathbf{V}_{n,\mathbf{\theta}}$. Let $\mathbf{w}$ be any $\left(N+4\right)$-vector and denote $\gamma_{\mathbf{w}}\left(u\right)=\sum_{j=-3}^{N}w_{j}B_{j,4}\left(u\right)$, then $\mathbf{B}_{\mathbf{\theta}}\mathbf{w}=\left\{\gamma_{\mathbf{w}}\left(U_{\mathbf{\theta},1}\right),...,\gamma_{\mathbf{w}}\left(U_{\mathbf{\theta},n}\right)\right\}^{T}$ and the definition of $A_{n}$ in (A.5) from Lemma A.3 entails that

[TABLE]

Using Theorem 5.4.2 of DeVore and Lorentz (1993) and Assumption A2, one obtains that

[TABLE]

which, together with (A.13), yield

[TABLE]

Now the order of $A_{n}$ in (A.5), together with (A.14) and (A.15) implies (A.11), in which $c_{V}=c_{f}C,C_{V}=C_{f}C$. Next, denote by $\lambda_{\max}\left(\mathbf{V}_{n,\mathbf{\theta}}\right)$ and $\lambda_{\min}\left(\mathbf{V}_{n,\mathbf{\theta}}\right)$ the maximum and minimum eigenvalue of $\mathbf{V}_{n,\mathbf{\theta}}$, simple algebra and (A.11) entail that

[TABLE]

thus

[TABLE]

Meanwhile, let $\mathbf{w}_{j}=$ the $\left(N+4\right)$-vector with all zeros except the $j$-th element being $1,j=-3,...,N$. Then clearly

[TABLE]

and in particular

[TABLE]

This, together with (A.5) yields that

[TABLE]

which leads to $\kappa\geq C>1,a.s.$ because the definition of B-spline and Assumption A2 ensure that $\left\|B_{0,4}\right\|_{\mathbf{\theta}}^{2}\geq C_{0}\left\|B_{-3,4}\right\|_{\mathbf{\theta}}^{2}$ for some constant $C_{0}>1$. Next applying Lemma A.4 with $\nu=\left(\kappa^{2}-1\right)^{1/16}\left(\kappa^{2}+1\right)^{-1/16}$ and $c_{0}=\nu^{-8}\left\|\mathbf{V}_{n,\mathbf{\theta}}^{-1}\right\|_{2}$, one gets $\left\|\mathbf{V}_{n,\mathbf{\theta}}^{-1}\right\|_{\infty}\leq 2\nu^{-8}N\left(1-\nu\right)^{-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 $Q_{T}\left(m\right)$ 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.

[TABLE]

[TABLE]

Proof. According to Theorem A.1 of Huang (2003), there exists an absolute constant $C>0$, such that

[TABLE]

which proves (A.16) for the case $k=0$. Applying Lemma A.6, one has for $0\leq k\leq 2$

[TABLE]

As a consequence of (A.17) and (A.18) for the case $k=0$, one has

[TABLE]

which, according to the differentiation of B-spline given in de Boor (2001), entails that

[TABLE]

Combining (A.18) and (A.19) proves (A.16) for $k=1,2.$

[TABLE]

[TABLE]

Proof. According to the definition of $\tilde{\gamma}_{\mathbf{\theta}}$ in (A.7), and the fact that $Q_{T}\left(\gamma_{\mathbf{\theta}}\right)$ is a cubic spline on the knots $T$

[TABLE]

which entails that

[TABLE]

Since

[TABLE]

applying (A.19) to the decomposition above produces (A.20). The proof of (A.21) is similar.

[TABLE]

[TABLE]

Proof. To prove (A.22), observe that for any vector $\mathbf{a}\in R^{n}$, with probability $1$

[TABLE]

[TABLE]

To prove (A.23), one only needs to use (A.12), (A.22) and (3.1).

[TABLE]

[TABLE]

Similarly, under Assumptions A2, A4-A6, with probability 1

[TABLE]

[TABLE]

Proof. We decompose the noise variable $\varepsilon_{i}$ into a truncated part and a tail part $\varepsilon_{i}=\varepsilon_{i,1}^{D_{n}}+\varepsilon_{i,2}^{D_{n}}+m_{i}^{D_{n}}$, where $D_{n}=n^{\eta}\left(1/3<\eta<2/5\right)$, $\varepsilon_{i,1}^{D_{n}}=\varepsilon_{i}I\left\{\left|\varepsilon_{i}\right|>D_{n}\right\}$,

[TABLE]

It is straightforward to verify that the mean of the truncated part is uniformly bounded by $D_{n}^{-2}$, so the boundedness of B spline basis and of the function $\sigma^{2}$ entail that

[TABLE]

The tail part vanishes almost surely

[TABLE]

Borel-Cantelli Lemma implies that

[TABLE]

For the truncated part, using Bernstein’s inequality and discretization as in Lemma A.2

[TABLE]

Therefore (A.24) is established as with probability $1$

[TABLE]

The proofs of (A.25), (A.26) are similar as $E\left\{m\left(\mathbf{X}_{i}\right)-\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta},i}\right)\left|U_{\mathbf{\theta},i}\right.\right\}\equiv 0$, but no truncation is needed for (A.26) as $\sup\limits_{\mathbf{\theta}\in S_{c}^{d-1}}\max\limits_{1\leq i\leq n}\left|m\left(\mathbf{X}_{i}\right)-\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta},i}\right)\right|\leq C<\infty$. Meanwhile, to prove (A.27), we note that for any $p=1,...,d$

[TABLE]

According to (2.6), one has $\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta}}\right)\equiv E\left\{m\left(\mathbf{X}\right)|U_{\mathbf{\theta}}\right\}$, hence

[TABLE]

Applying Assumptions A2 and A3, one can differentiate through the expectation, thus

[TABLE]

which allows one to apply the Bernstein’s inequality to obtain that with probability $1$

[TABLE]

which is (A.27).

[TABLE]

Proof. Denote $\mathbf{\hat{a}}\equiv\left(\hat{a}_{-3},\cdots,\hat{a}_{N}\right)^{T}\mathbf{=}\left(\mathbf{B}_{\mathbf{\theta}}^{T}\mathbf{B}_{\mathbf{\theta}}\right)^{-1}\mathbf{B}_{\mathbf{\theta}}^{T}\mathbf{E}=\mathbf{V}_{n,\mathbf{\theta}}^{-1}\left(n^{-1}\mathbf{\mathbf{B}_{\mathbf{\theta}}^{T}\mathbf{E}}\right)$, then $\hat{\varepsilon}_{\mathbf{\theta}}\left(u\right)=\sum_{j=-3}^{N}\hat{a}_{j}B_{j,4}\left(u\right)$, so the order of $\hat{\varepsilon}_{\mathbf{\theta}}\left(u\right)$ is related to that of $\mathbf{\hat{a}}$. In fact, by Theorem 5.4.2 in DeVore and Lorentz (1993)

[TABLE]

where the last inequality follows from (A.12) of Lemma A.5. Applying (A.24) of Lemma A.10, we have established (A.28).

[TABLE]

The proof is similar to Lemma A.11, thus omitted.

The next result evaluates the uniform size of the noise derivatives.

[TABLE]

Proof. Note that

[TABLE]

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

[TABLE]

in which

[TABLE]

[TABLE]

By (A.24), (A.12), (A.22) and (A.23), one derives

[TABLE]

while (A.27) of Lemma A.10, (A.12) of Lemma A.5

[TABLE]

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)

[TABLE]

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

[TABLE]

It is clear from (A.20), (A.30) and (A.31) that with probability $1$

[TABLE]

[TABLE]

Putting together all the above yields (2.14). The proof of (2.15) is similar.

A.3. Proof of Proposition 2.2

[TABLE]

Proof. For the empirical risk function $\hat{R}\left(\mathbf{\theta}\right)$ in (2.9), one has

[TABLE]

[TABLE]

hence

[TABLE]

[TABLE]

[TABLE]

where $\hat{\gamma}_{\mathbf{\theta}}\left(x\right)$ is defined in (2.8). Using the expression of $R\left(\mathbf{\theta}\right)$ in (2.7), one has

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Bernstein inequality and strong law of large number for $\alpha$ mixing sequence imply that

[TABLE]

Now (2.13) of Proposition 2.1 provides that

[TABLE]

which entail that

[TABLE]

[TABLE]

Hence

[TABLE]

The lemma now follows from (A.37), (A.38) and (A.39) and Assumption A6.

[TABLE]

in which

[TABLE]

Furthermore for $k=1,2$

[TABLE]

Proof. Note that for any $p=1,2,...,d$

[TABLE]

[TABLE]

Thus $E\left(\xi_{\mathbf{\theta},i,p}\right)=2E\left[\left\{\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta},i}\right)-Y_{i}\right\}\frac{\partial}{\partial\theta_{p}}\gamma_{\mathbf{\theta}}\left(U_{\mathbf{\theta},i}\right)\right]-\frac{\partial}{\partial\theta_{p}}R\left(\mathbf{\theta}\right)=0$ and

[TABLE]

with

[TABLE]

Bernstein inequality implies that

[TABLE]

Meanwhile, applying (2.13) and (2.14) of Proposition 2.1, one obtains that

[TABLE]

Note that

[TABLE]

Applying (2.13), one gets

[TABLE]

while (A.24), (A.26) and (A.12) entail that with probability $1$

[TABLE]

[TABLE]

thus

[TABLE]

Lastly

[TABLE]

By applying (A.24), (A.26), and (A.12), it is clear that with probability $1$

[TABLE]

while by applying (A.16) of Lemma A.7, one has

[TABLE]

together, the above entail that

[TABLE]

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 $\hat{R}\left(\mathbf{\theta}\right)$ and $R\left(\mathbf{\theta}\right)$ with respect to $\theta_{p}$, $\theta_{q}$ are

[TABLE]

[TABLE]

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 $\hat{S}_{p}^{*}\left(\mathbf{\theta}_{-d}\right)$ be the $p$-th element of $\hat{S}^{*}\left(\mathbf{\theta}_{-d}\right)$ and for $\gamma_{\mathbf{\theta}}$ in (2.6), denote

[TABLE]

where $\dot{\gamma}_{p}$ is value of $\frac{\partial}{\partial\theta_{p}}\gamma_{\mathbf{\theta}}$ taking at $\mathbf{\theta=\theta}_{0}$, for any $p,q=1,2,...,d-1$.

[TABLE]

Proof. For any $p=1,...,d-1$

[TABLE]

Therefore, according to (A.40), (A.41) and (A.48)

[TABLE]

[TABLE]

Since $S^{*}\left(\mathbf{\theta}_{-d}\right)$ attains its minimum at $\mathbf{\theta}_{0,-d}$, for $p=1,...,d-1$

[TABLE]

which yields (A.49).

Proof. It is easy to show that for any $p,q=1,2,...,d$,

[TABLE]

[TABLE]

Note that

[TABLE]

[TABLE]

Thus

[TABLE]

[TABLE]

[TABLE]

Therefore we obtained the desired result.

Proof of Theorem 2. For any $p=1,2,...,d-1$, let

[TABLE]

then

[TABLE]

Note that $\hat{S}^{*}\left(\mathbf{\theta}_{-d}\right)$ attains its minimum at $\hat{\mathbf{\theta}}_{-d}$, i.e., $\hat{S}_{p}^{*}\left(\hat{\mathbf{\theta}}_{-d}\right)\equiv 0$. Thus, for any $p=1,2,...,d-1$, $t_{p}\in\left[0,1\right]$, one has

[TABLE]

then

[TABLE]

Now (2.11) of Theorem 1 and Proposition 2.2 with $k=2$ imply that uniformly in $p,q=1,2,...,d-1$

[TABLE]

where $l_{p,q}$ is given in Theorem 2. Noting that $\sqrt{n}\left(\hat{\mathbf{\theta}}_{-d}\mathbf{-\theta}_{0,-d}\right)$ is represented as

[TABLE]

where $\hat{S}^{*}\left(\mathbf{\theta}_{0,-d}\right)=\left\{\hat{S}_{p}^{*}\left(\mathbf{\theta}_{0,-d}\right)\right\}_{p=1}^{d-1}$ and according to (A.48) and Lemma A.16

[TABLE]

Let $\Psi\left(\mathbf{\theta}_{0}\right)=\left(\psi_{pq}\right)_{p,q=1}^{d-1}$ be the covariance matrix of $\sqrt{n}\left\{\hat{S}_{p}^{*}\left(\mathbf{\theta}_{0,-d}\right)\right\}_{p=1}^{d-1}$ with $\psi_{pq}$ given in Theorem 2. Cramér-Wold device and central limit theorem for $\alpha$ mixing sequences entail that

[TABLE]

Let $\Sigma\left(\mathbf{\theta}_{0}\right)=\left\{H^{*}\left(\mathbf{\theta}_{0,-d}\right)\right\}^{-1}\Psi\left(\mathbf{\theta}_{0}\right)\left[\left\{H^{*}\left(\mathbf{\theta}_{0,-d}\right)\right\}^{T}\right]^{-1}$, with $H^{*}\left(\mathbf{\theta}_{0,-d}\right)$ being the Hessian matrix defined in (2.3). The above limiting distribution of $\sqrt{n}\hat{S}^{*}\left(\mathbf{\theta}_{0,-d}\right)$, (A.52) and Slutsky’s theorem imply that

[TABLE]

References

Bosq, D. (1998). Nonparametric Statistics for Stochastic Processes. Springer-Verlag, New York.

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.

Härdle, W. (1990). Applied Nonparametric Regression. Cambridge University Press, Cambridge.

Härdle, W. and Hall, P. and Ichimura, H. (1993). Optimal smoothing in single-index models. Ann. Statist. 21 157-178.

Hä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 Hä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.

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.