Local Quadratic Estimation of the Curvature in a Functional Single Index Model

TL;DR

This paper develops a method for estimating the curvature of the link function in a functional single index model, addressing the challenge of unknown coefficient and link functions, with theoretical convergence results and practical considerations.

Contribution

It introduces a nested optimization approach for estimating the second derivative of the link function in a functional single index model, including convergence rates and root-n consistency.

Findings

The second derivative of the link function can be estimated with a specific convergence rate.

The argument of the link function can be estimated with root-n consistency.

Estimates are sensitive to initial values in practical implementation.

Abstract

The nonlinear effects of environmental variability on species abundance plays an important role in the maintenance of ecological diversity. Nonetheless, many common models use parametric nonlinear terms pre-determining ecological conclusions. Motivated by this concern, we study the estimate of the second derivative (curvature) of the link function g in a functional single index model. Since the coefficient function and the link function are both unknown, the estimate is expressed as a nested optimization. For a fixed and unknown coefficient function, the link function and its second derivative are estimated by local quadratic approximation, then the coefficient function is estimated by minimizing the MSE of the model. In this paper, we derive the rate of convergence of the estimation. In addition, we prove that the argument of g, can be estimated root-n consistently. However, practical implementation of the method requires solving a nonlinear optimization problem, and our results show that the estimates of the link function and the coefficient function are quite sensitive to the choices of starting values.

Tables

Table 1. Table 1: Simulation results by GCV using rescaled bandwidth h ​ ‖ 𝒄 ‖ ℎ norm 𝒄 h\left\|\bm{c}\right\| and constraint β ^ ^ 𝛽 \hat{\beta} .

		g1			g2			g3
Initial	n	RSE	RASE	RASE2	RSE	RASE	RASE2	RSE	RASE	RASE2
True	100	0.2887	0.0959	9.4156	0.2785	0.0945	7.2648	0.3442	0.0918	5.0951
	1000	0.0829	0.0336	0.9003	0.0731	0.0334	0.5429	0.1095	0.0276	0.8620
Linear	100	1.1401	0.3401	1.3869	1.3569	0.4166	1.2014	0.7049	0.1393	0.2507
	1000	1.9010	0.2340	0.5205	1.0719	0.3421	1.3075	0.5287	0.0891	0.0729
Equal	100	0.8389	0.2820	1.0933	0.8342	0.2566	1.4349	0.8329	0.1857	0.2589
	1000	0.8502	0.3159	1.0133	0.8437	0.2838	1.4043	0.8508	0.1921	0.0905
Random	100	1.1911	0.1461	0.5597	1.1617	0.1979	0.6186	1.2603	0.0942	0.2043
	1000	1.2704	0.1320	0.4354	1.3021	0.1890	0.4413	1.2300	0.0796	0.0637

Table 2. Table 2: Simulation results of random starting value with rescaled and original bandwidths selected by GCV

	g1		g2		g3
n	original	rescaled	original	rescaled	original	rescaled
100	1.2793	0.5597	2.1497	0.6186	0.2360	0.2043
1000	1.1297	0.4354	2.2081	0.4413	0.0703	0.0637

Table 3. Table 3: Simulation results by 10 10 10 -fold CV using rescaled bandwidth h ​ ‖ 𝒄 ‖ ℎ norm 𝒄 h\left\|\bm{c}\right\| and constraint β ^ ^ 𝛽 \hat{\beta} .

		g1			g2			g3
Initial	n	RSE	RASE	RASE2	RSE	RASE	RASE2	RSE	RASE	RASE2
True	100	0.1921	0.0868	1.2106	0.2052	0.0904	4.5719	0.1982	0.0690	0.7954
	1000	0.0680	0.0290	0.4551	0.0620	0.0271	0.4016	0.0774	0.0243	0.2762
Linear	100	1.1383	0.3294	1.4303	1.3370	0.4245	1.3388	0.6995	0.1331	0.2692
	1000	1.9015	0.2232	0.5060	1.0419	0.3185	1.3652	0.5297	0.0852	0.0786
Equal	100	0.8486	0.2789	1.0731	0.8385	0.2509	1.4728	0.8349	0.1836	0.3349
	1000	0.8442	0.2977	1.0181	0.8500	0.2755	1.5396	0.8447	0.1827	0.1045
Random	100	1.2998	0.1339	0.5191	1.2086	0.1979	0.7160	1.4296	0.0875	0.2362
	1000	1.5607	0.1288	0.4209	1.3540	0.1687	0.4418	1.7622	0.0798	0.0744

Table 4. Table 4: Simulation results of random starting value with rescaled and original bandwidths selected by 10 10 10 -fold cross-validation

	g1		g2		g3
n	original	rescaled	original	rescaled	original	rescaled
100	1.2710	0.5191	2.2005	0.7160	0.2768	0.2362
1000	1.0145	0.4209	2.2812	0.4418	0.0752	0.0744

Table 5. Table 5: Simulation results of true starting value with rescaled and original bandwidths selected by 10 10 10 -fold cross-validation

	g1		g2		g3
n	original	rescaled	original	rescaled	original	rescaled
100	5.9302	1.9266	33.811	7.9015	4.5630	0.9791
1000	2.4004	1.0660	2.2494	1.3416	1.2444	0.3558

Table 6. Table 6: Simulation results of true starting value with rescaled and original bandwidths selected by GCV

	g1		g2		g3
n	original	rescaled	original	rescaled	original	rescaled
100	51.397	7.3786	38.007	6.5712	50.816	5.9508
1000	9.0220	1.5416	8.9590	1.8170	7.8245	1.1493

Table 7. Table 7: CV values

		g1		g2		g3
Initial	n	GCV	10-fold	GCV	10-fold	GCV	10-fold
True	100	0.0398	3.8686	0.0399	3.3823	0.0377	1.4970
	1000	0.0401	0.9468	0.0402	0.7882	0.0400	0.4781
Linear	100	0.1559	5.2259	0.2184	4.7373	0.0498	1.6148
	1000	0.0793	2.0295	0.1272	2.2896	0.0439	0.7265
Equal	100	0.1401	4.7640	0.1046	3.9670	0.0795	1.8530
	1000	0.1358	2.6925	0.1318	2.1995	0.0751	1.1081
Random	100	0.0611	4.2920	0.0957	3.9025	0.0459	1.6237
	1000	0.0591	2.4463	0.0850	2.2702	0.0468	1.0243

Equations

$$\mathbf{G}=g\left(\mathbf{E}\right)+ϵ,$$

$$\mathbf{E}={\int }_0^{1}X\left(t\right){\beta }^0\left(t\right)\mathrm{d}t,$$

$$Y=g\left({\int }_0^{1}X\left(t\right){\beta }^0\left(t\right)\mathrm{d}t\right)+ϵ.$$

$$Y={\int }_0^{1}X\left(t\right){\beta }^0\left(t\right)\mathrm{d}t+ϵ,$$

$$Y=g\left(\mathbf{X}{\mathbf{\beta }}^0\right)+ϵ,$$

$${\hat{g}}_i\left({\mathbf{X}}_i\mathbf{\beta }|\mathbf{\beta },h\right)=\frac{{\sum }_{j\ne i}{Y}_{j}K\left(\frac{{\mathbf{X}}_i\mathbf{\beta }-{\mathbf{X}}_j\mathbf{\beta }}{h}\right)}{{\sum }_{j\ne i}K\left(\frac{{\mathbf{X}}_i\mathbf{\beta }-{\mathbf{X}}_j\mathbf{\beta }}{h}\right)},$$

$$\hat{S}\left(\mathbf{\beta },h\right)=\sum _{i=1}^n{\left[Y_i-{\hat{g}}_i\left({\mathbf{X}}_i\mathbf{\beta }|\mathbf{\beta },h\right)\right]}^2.$$

$$\mathrm{E}\left[\frac{\partial g\left(\mathbf{X\beta }\right)}{\partial \mathbf{X}}\right]\stackrel{z\doteq \mathbf{X}\mathbf{\beta }}{=}\mathrm{E}\left[\frac{\partial g\left(z\right)}{\partial z}\mathbf{\beta }\right]=\mathrm{E}\left[\frac{\partial g\left(z\right)}{\partial z}\right]\mathbf{\beta }\doteq \gamma \mathbf{\beta }.$$

$$\frac{1}{n}\sum _{j=1}^n{\left[\hat{g}\left(\int X_j\hat{\beta }\right)-g\left(\int X_j{\beta }^0\right)\right]}^2=O\left(n^{-c}\right),$$

$$Y=g\left({\int }_0^{1}X\left(t\right){\beta }^0\left(t\right)\mathrm{d}t\right)+ϵ,$$

$$g\left(x\right)\approx g\left(x_0\right)+g^{\prime }\left(x_0\right)\left(x-x_0\right)+g^{\prime \prime }\left(x_0\right)\frac{{\left(x-x_0\right)}^2}{2}.$$

$$\left(\hat{a},\hat{b},\hat{c}\right)=$$

$$\underset{a,b,c}{\inf }\sum _{i=1}^n\left\{{\left[Y_i-a-b\left({\int }_0^1{X}_i\left(t\right)\beta \left(t\right)\mathrm{d}t-u\right)-c\frac{{\left({\int }_0^1{X}_i\left(t\right)\beta \left(t\right)\mathrm{d}t-u\right)}^2}{2}\right]}^{2}K\left(\frac{{\int }_0^1{X}_i\left(t\right)\beta \left(t\right)\mathrm{d}t-u}{h_n}\right)\right\},$$

$$\hat{\beta }=\underset{\beta \in \mathcal{B}}{\inf }\sum _{i=1}^n{\left[Y_i-\hat{g}\left(\int X_i\beta \right)\right]}^2,$$

$$\hat{g}\left(\int X_i\beta \right)=\frac{\sum _{j:j\ne i}{Y}_{j}K\left(\frac{\int X_i\beta -\int X_j\beta }{h_n}\right)}{\sum _{j:j\ne i}K\left(\frac{\int X_i\beta -\int X_j\beta }{h_n}\right)}.$$

$${\left(\hat{g}\left(\int X_j\beta \right),{\hat{g}}^{\prime }\left(\int X_j\beta \right),{\hat{g}}^{\prime \prime }\left(\int X_j\beta \right)\right)}^{\top }={\left({\mathbb{X}}_{\beta ,j}^{\top }{\mathbb{K}}_{\beta ,j}{\mathbb{X}}_{\beta ,j}\right)}^{-1}\left({\mathbb{X}}_{\beta ,j}^{\top }{\mathbb{K}}_{\beta ,j}\right)\mathbf{Y},$$

$${\hat{g}}^{\prime }\left(\int X_j\beta \right)$$

$${\hat{g}}^{\prime \prime }\left(\int X_j\beta \right)$$

$${\mathbb{X}}_{\beta ,j}=\left(1,\int \mathbf{X}\beta -\left(\int X_j\beta \right)1,\frac{{\left(\int \mathbf{X}\beta -\left(\int X_j\beta \right)1\right)}^2}{2}\right),$$

$${\mathbb{K}}_{\beta ,j}=\text{diag}\left(K\left(\frac{\int X_1\beta -\int X_j\beta }{h_n}\right),\cdots ,K\left(\frac{\int X_n\beta -\int X_j\beta }{h_n}\right)\right).$$

$$\left|g^{\left(k\right)}\left(u\right)-g^{\left(k\right)}\left(v\right)\right|\le D_2\left|u-v\right|,$$

$$X_i\left(t\right)=\sum _{j=1}^{\infty }{c}_{ij}{\varphi }_j\left(t\right),$$

$$\beta \left(t\right)=\sum _{j=1}^{\infty }{b}_j{\varphi }_j\left(t\right).$$

$$\hat{\beta }\left(t\right)=\sum _{j=1}^{\infty }{\hat{b}}_j{\varphi }_j\left(t\right),{\beta }^0\left(t\right)=\sum _{j=1}^{\infty }{b}_j^0{\varphi }_j\left(t\right).$$

$$\int X_i\beta =\sum _{j=1}^{\infty }{c}_{ij}{b}_j.$$

$$\sum _{j=p_n+1}^{\infty }{c}_{ij}{b}_j=O\left(p_n^{-\lambda }\right),$$

$${\mathbb{X}}_{\beta }^{\top }{\mathbb{K}}_{\beta }=$$

$$=$$

$$T_j^p\left(\beta \right)=$$

$$T_j^0\left(\beta \right)=$$

$$=$$

Full text

1 Introduction

1.1 Ecological Motivation

Within mathematical ecology, nonlinear responses to environmental variability play an important role in maintaining the diversity of competing species. Species competing for the same resources can nonetheless co-exist by exploiting differing environmental conditions; see Hutchinson (1961), Chesson and Warner (1981) and Ellner (1987). For an individual species, environmental fluctuation can accelerate growth rate (Drake, 2005; Koons et al., 2009) or sometimes decrease long-term population growth rates (Lewontin and Cohen, 1969). The nonlinearity of these responses also plays an important role in forecasting the effect of increased environmental variability under climate change. The motivating data for this study come from long-term observations of communities of prairie plants in which Artemesia Triparta – sage brush – is a dominant species and we wish to understand its responses to climate given by temperature and precipitation.

Traditional statistical models for plant growth make parametric assumptions that imply specific forms of nonlinearity, particularly in the presence of high-dimensional covariates. Instead, we use a nonparametric growth model of an individual plant or animal:

[TABLE]

where $\bm{G}$ and $\bm{E}$ are the growth and environment of a plant, $g$ is a link function to be estimated and $\epsilon$ is the random error. To answer the ecological question, “Would the growth be higher if we just gave the plant a constant environment at the average of $\bm{E}$?”, we need to compare $g\left[\mathrm{E}\left(\bm{E}\right)\right]$ and $\mathrm{E}\left[g\left(\bm{E}\right)\right]$.

If the link function $g$ is convex, $g\left[\mathrm{E}\left(\bm{E}\right)\right]\leq\mathrm{E}\left[g\left(\bm{E}\right)\right]$ by Jensen’s inequality, and the plant grows better in a varying environment. Otherwise, if the link function $g$ is concave, a constant environment is preferred. Assuming a smooth function $g$, convexity is equivalent to $g^{\prime\prime}\left(s\right)>0$, for all $s$ in the domain of $g$. Therefore, in this paper we consider the problem of estimating the curvature of the link function $g$.

To finalize this model, the environment $\bm{E}$ is described by the recent history of temperature and rainfall recorded at up to daily resolution. Since plants may be impacted by climate events over a long period of time (see Dahlgren and Ehrlén, 2011; Clark et al., 2011), we will consider the past two years of data. Following Teller et al. (2016), these are thought of as functional covariates leading to a representation of $\bm{E}$ as a functional linter term:

[TABLE]

where $\beta^{0}\left(t\right)$ is the coefficient function to be estimated, and $X\left(t\right)$ is the covariate function we observed, typically a measurement of climate history.

The growth model of a plant is now given by

[TABLE]

This is Functional Single Index model, introduced in Chen et al. (2011) and Ma (2016).

In functional data analysis, a functional linear model (FLM) is defined as

[TABLE]

which is often used in modeling the relationship between a functional covariate and a scalar response. To assess curvature, we need a more flexible model than the FLM. A generalized functional linear model (GFLM) is proposed in Müller and Stadtmüller (2005), James (2002) and Escabias et al. (2007), which has the same form as the functional single index model but with a known link function $g$. The functional single index model could be considered as an extension to the GFLM, as it is more flexible and could model a variety of real-world data.

Compared to the generalized functional linear model, estimation of the link function $g$ based on a unknown coefficient function $\beta^{0}$ is challenging. Even if $\beta^{0}$ is known, estimating the second derivative of a nonparametric function directly is difficult. In this paper, we prove a theoretical convergence rate for an estimate of $g^{\prime\prime}$ in the functional single index model, even if there are some bias in estimating the coefficient function $\beta^{0}$.

The convergence rates that we derive are based on finding a global solution to a nonlinear optimization problem using a bandwidth that decreases at a known rate with $n$. However, this requires overcoming several practical issues. First, to find an optimum for $\beta^{0}$, we rely on nonlinear optimization methods which require an initial value from which to search for a minimum. Our experiments demonstrate that the performance of the estimate can depend critically on this choice of initial condition and natural choices which provide good estimates of $g$ do not necessarily work well for $g^{\prime\prime}$. Further, the optimal choice of bandwidth can be quite different between estimation targeting $g$ and that targeting $g^{\prime\prime}$ and we provide a heuristic post-cross-validation modification to improve the estimate of bandwidth. We expect similar rates of convergence will hold for alternative non-parametric estimators of $g$, penalized splines, for example, but that the specifics of smoothing parameter selection and nonlinear optimization can be expected to be quite different. A detailed analysis of the optimization problem is beyond the scope of this paper.

1.2 Previous Results

In this section, we will introduce previous theoretical and empirical results for Single Index model and Functional Single Index model.

1.2.1 Single Index Model

There has been considerable research on the single index model, where the coefficient $\bm{\beta}^{0}$ is finite dimensional. The Single Index Model is defined as

[TABLE]

where $\bm{X}$ is the covariate and $\bm{\beta}^{0}$ is the coefficient vector. There are three methods to estimate the link function $g$ and the coefficient vector $\bm{\beta}^{0}$. The Projection Pursuit Regression (PPR) approach introduced in Hardle et al. (1993) is a two-step estimation procedure:

Estimate the link function $g$ by the kernel method

[TABLE]

where $h$ is the bandwidth. 2. 2.

Estimate the coefficient $\bm{\beta}^{0}$ by minimizing the mean squared error

[TABLE]

Hardle et al. (1993) proved that the coefficient vector $\bm{\beta}^{0}$ can be estimated root-$n$ consistently. Ichimura (1993) showed the asymptotic normality of the estimator. The other two approaches provide new methods to estimate the coefficient vector. The Average Derivative approach in Hristache et al. (2001) showed that

[TABLE]

If we could find a consistent estimator of the average derivative $\mathrm{E}\left[\frac{\partial g\left(\bm{X\beta}^{0}\right)}{\partial\bm{X}}\right]$, we can get a consistent estimator of the coefficient $\bm{\beta}^{0}$ up to a scale. Normally, we require the coefficient vector to be norm $1$. Stoker (1986) proposed two consistent estimators of the average derivative.

The sliced inverse regression method in Li (1991) considered the estimation of the coefficient vector as a dimension-reduction problem. Any linear combination of the coefficient vector $\bm{\beta}^{0}$ is assumed to be an effective dimension-reduction (EDR) direction. They conduct a principle component analysis on the inverse regression space $\mathrm{E}\left(\bm{X}\left|Y\right.\right)$, and estimate the coefficient vector $\bm{\beta}^{0}$ by the largest component.

1.2.2 Functional Single Index Model

There are only a few papers in the functional single index model. In Chen et al. (2011), similar to the projection pursuit regression in the single index model, the link function $g$ and the coefficient function $\beta^{0}$ are estimeated by a two-step procedure. The coefficient function $\beta^{0}$ is reduced to a finite dimensional coefficient vector by a spline basis. Under some assumptions, Chen et al. (2011) showed that

[TABLE]

for $c>0$. In Ma (2016), two spline basis were used to represent the coefficient function and the link function, respectively, and the MSE was minimized iteratively until convergence. Ma (2016) constructed a asymptotic simultaneous confidence band for the coefficient function $\beta^{0}$. Our estimates follow Chen et al. (2011) but will examine the properties of $\hat{g}^{\prime\prime}$. By a clever decomposition of squared error, Chen et al. (2011) were able to avoid the need to directly account for the estimate of $\beta^{0}$. Unlike that case, to examine $g^{\prime\prime}$ we will need to obtain the $\sqrt{n}$ convergence rate for $\int X\beta^{0}$ directly, before we can examine our target.

2 Estimation Procedure

Suppose that we observe $n$ environment histories and responses $\left(X_{1}\left(t\right),Y_{1}\right),\cdots,\left(X_{n}\left(t\right),Y_{n}\right)$, independent and identically distributed as $\left(X\left(t\right),Y\right)$, where $t\in\left[0,1\right]$, with

[TABLE]

where $Y$ is the scalar response variable, and $X\left(t\right)$ is the covariate function. For the purpose of simplification, we assume that the predicator $X$ and the coefficient function $\beta^{0}$ are defined in the domain $\left[0,1\right]$, and $\epsilon$ is a Gaussian random error.

To answer our ecological question, we are interested in estimating the second derivative (curvature) of the link function $g$. The estimate of the coefficient function $\beta^{0}$ is denoted as $\hat{\beta}$. Define a Hilbert space $\mathcal{B}$ as the set of the coefficient functions $\beta$, where $\hat{\beta}\left(t\right),\beta^{0}\left(t\right)\in\mathcal{B}$.

To estimate $\beta^{0}$, $g$ and $g^{\prime\prime}$, we use a local quadratic approximation. By Taylor’s expansion, at a fixed point $x$, the link function $g$ can be approximated by

[TABLE]

Fix $u$, where $u$ is in the domain of the link function $g$, the curvature, denoted as $\hat{g}^{\prime\prime}$, is estimated by minimizing the weighted sum of squares

[TABLE]

where $K$ is a kernel function and $h_{n}$ is the bandwidth. The estimators are then $\left(\hat{g}\left(u\right),\hat{g}^{\prime}\left(u\right),\hat{g}^{\prime\prime}\left(u\right)\right)=\left(\hat{a},\hat{b},\hat{c}\right)$.

The coefficient function $\beta^{0}\left(t\right)$ is unknown in the penalized weighted sum of square $\left(\ref{a}\right)$. We estimate it by minimizing the MSE

[TABLE]

where

[TABLE]

Since the kernel function $K$ is only defined in $\left[-1,1\right]$, we constraint the domain of the estimate of $g$ or $g^{\prime\prime}$ to be in $\left[-1,1\right]$ by normalizing the coefficients of $g$ or $g^{\prime\prime}$ under an orthonormal basis to be $1$ after optimization procedure.

Denote a column vector $\bm{Y}=\left(Y_{1},\cdots,Y_{n}\right)^{\top}$. Fix $j\in\left\{1,\cdots,n\right\}$, and $\beta\in\mathcal{B}$, the estimated $\hat{g}\left(\int X_{j}\beta\right)$, $\hat{g}^{\prime}\left(\int X_{j}\beta\right)$ and $\hat{g}^{\prime\prime}\left(\int X_{j}\beta\right)$ can be calculated as

[TABLE]

where

[TABLE]

where $\mathbb{A}_{k\cdot}$ denotes the $k^{\text{th}}$ row of a matrix $\mathbb{A}$, and the $\left(n\times 3\right)$-dimensional matrix $\mathbb{X}_{\beta,j}$ is

[TABLE]

with $\int\bm{X}\beta\doteq\left(\int X_{1}\beta,\cdots,\int X_{n}\beta\right)^{\top}$, $\bm{1}$ is a $n$-dimensional column vector of ones, and the $\left(n\times n\right)$-dimensional matrix $\mathbb{K}_{\beta,j}$ is

[TABLE]

The estimation of the coefficient function $\beta^{0}$ and the link function $g$ is therefore a nested procedure, summarized in $\left(\ref{step1}\right)$, $\left(\ref{step2}\right)$ and $\left(\ref{step3}\right)$. Following Ma (2016), the identifiability of the model is ensured by adding a constraint on the coefficient function, such that $\int_{0}^{1}\beta^{2}\left(t\right)\mathrm{d}t=1.$

3 Assumptions

In deriving a convergence rate for $\frac{1}{n}\sum\limits_{i=1}^{n}\mathrm{E}\left[\hat{g}^{\prime\prime}\left(\int X_{i}\hat{\beta}\right)-g^{\prime\prime}\left(\int X_{i}\beta^{0}\right)\right]^{2}$, we make the following assumptions in the functional single index model.

The observations $\left(X_{i}\left(t\right),Y_{i}\right)$, where $i=1,\cdots,n$, are independent and identically distributed. Each covariate function $X_{i}\left(t\right)$ is a square-integrable random function defined in the interval $\left[0,1\right]$. The random error $\epsilon$ is independent from $X$, and has zero mean and variance $\sigma^{2}$. 2. 2.

The dependent variable $Y$ has the $m$th-order absolute moment, where $m\geq 2$. This is an assumption from Ichimura (1993). The finite moment $m$ is used in establishing the main convergence theorem. 3. 3.

The link function $g$ and the curvature $g^{\prime\prime}$ are bounded and satisfy the Lipschitz condition such that

[TABLE]

for all $u$ and $v$, where $D_{2}>0$ and $k=0,2$. The Lipschitz condition ensures that if $\beta^{0}$ can be estimated root-$n$ consistently, the distance between $g^{\prime\prime}\left(\int X\beta^{0}\right)$ and $g^{\prime\prime}\left(\int X\hat{\beta}\right)$ can be controlled. 4. 4.

The kernel function $K$ is nonnegative and symmetric with support $\left[-1,1\right]$, and $\int_{-1}^{1}K\left(s\right)\mathrm{d}s=1$. Assume that $K$ is three times continuously differentiable, with $\left|K^{\left(3\right)}\left(s\right)\right|\leq D_{3}$, for any $s\in\left[-1,1\right]$ and $D_{3}>0$. Since the kernel $K$ satisfies a Lipschitz condition, the Nadaraya-Watson estimator of the link function $g$ also has a Lipschitz condition. 5. 5.

For some orthonormal basis $\left\{\phi_{k}\left(t\right):k=1,2,\cdots\right\}$, for each $i=1,\cdots,n$, there exists a sequence of random variables $\left\{c_{ij}\right\}_{j=1}^{\infty}$, such that

[TABLE]

and

[TABLE]

Assume that $\mathrm{E}\left(c_{ij}\right)=0$.

In particular, we have

[TABLE]

For any $\beta\in\mathcal{B}$, we can write

[TABLE]

We observe that an orthonormal basis approximation of the covariate function and coefficient function transforms an integration to an infinite sum. In addition, define a sequence $p_{n}$ such that $p_{n}\rightarrow\infty$ as $n\rightarrow\infty$, we require

[TABLE]

where $\lambda>0$, and $p_{n}=o\left(\frac{1}{h_{n}}\right)$. Condition $\left(\ref{dim}\right)$ ensures that the integration $\int X\beta$ can be approximated by a finite sum of coefficients under an orthonormal basis. 6. 6.

Assume that $\sup_{\beta\in\mathcal{B};x}f\left(x\left|\beta\right.\right)<\infty$, where $f\left(x\left|\beta\right.\right)$ is the probability density of $\int X\beta$.

4 Convergence Rates

By the definition of $\mathbb{X}_{\beta}^{\top}$ and $\mathbb{K}_{\beta}$, we can calculate

[TABLE]

Denote

[TABLE]

and

[TABLE]

Denote $u_{j}\doteq\int X_{j}\beta$, we have

[TABLE]

and

[TABLE]

where $\mu_{p}\left(K\right)\doteq\int_{-1}^{1}m^{p}K\left(m\right)\mathrm{d}m$. Since the kernel function $K$ is symmetric, $T_{j}^{p}\left(\beta\right)=0$ if $p$ is an odd number.

We have

[TABLE]

The determinant of the matrix $\mathbb{X}^{\top}_{\beta,j}\mathbb{K}_{\beta,j}\mathbb{X}_{\beta,j}$ is

[TABLE]

By definitions of $S_{1}\left(\beta;j\right)$ and $S_{2}\left(\beta;j\right)$ in $\left(\ref{s1}\right)$ and $\left(\ref{s2}\right)$, we can get

[TABLE]

4.1 Convergence Rate of $\beta$

For any $\beta\in\mathcal{B}$, define

[TABLE]

and

[TABLE]

Denote $\bm{c}_{i}=\left(c_{i1},\cdots\right)^{\top}$ and $\bm{b}=\left(b_{1},\cdots\right)^{\top}$, where $i=1,\cdots,n$. Define the subspaces $C,B\subset R^{\infty}$ such that $\bm{c}_{i}\in C$ and $\bm{b}\in B$.

• Lemma1.

For a sequence of positive numbers $M_{n}$, suppose that $nh_{n}^{2}\epsilon_{0n}^{2}M_{n}^{-2}\rightarrow\infty$, then

[TABLE]

where

[TABLE]

as $n\rightarrow\infty$.

Proof.

We have

[TABLE]

where

[TABLE]

By Assumption $4$, the kernel function $K$ is bounded. We will apply Bernstein’s inequality (see Appendix A) to the above equation with

[TABLE]

[TABLE]

and

[TABLE]

where $K_{1}$, $C_{1}$ and $C_{2}$ are constants. By Bernstein’s inequality, we can get

[TABLE]

The assumption of Lemma $A.5$ in Ichimura (1993) is that the sequence $\left\{M_{n}\right\}_{n=1}^{\infty}$ should satisfy $\epsilon_{0n}h_{n}M_{n}^{m-1}\rightarrow\infty$. Since $h_{n}\rightarrow 0$ and $\epsilon_{0n}\rightarrow 0$, we need $M_{n}\rightarrow\infty$. Therefore, in the denominator, we have $\frac{4}{3}M_{n}K_{1}\epsilon_{0n}=o\left(\frac{2M_{n}^{2}K_{2}}{h_{n}}\right)$. If $nh_{n}^{2}\epsilon_{0n}^{2}M_{n}^{-2}\rightarrow\infty$, then

[TABLE]

and

[TABLE]

∎

Following similar arguments, we can derive Lemmas $2$ and $3$ below.

• Lemma2.

For a sequence of positive numbers $M_{n}$, suppose that $nh_{n}^{4}\epsilon_{1n}^{2}M_{n}^{-2}\rightarrow\infty$, then

[TABLE]

where

[TABLE]

as $n\rightarrow\infty$.

• Lemma3.

For a sequence of positive numbers $M_{n}$, suppose that $nh_{n}^{6}\epsilon_{2n}^{2}M_{n}^{-2}\rightarrow\infty$, then

[TABLE]

where

[TABLE]

as $n\rightarrow\infty$.

Now, we show the root-$n$ consistency of the estimator $\hat{\beta}$.

• Theorem4.

The estimator $\hat{\beta}$ is consistent if Assumptions $1-6$ hold, and the bandwidth sequence satisfies $nh_{n}^{8}\rightarrow 0$ and $nh_{n}^{6}\rightarrow\infty$.

Proof.

Theorem $5.1$ in Ichimura (1993) states the consistency of $\hat{\beta}$, while $\beta^{0}$ is a coefficient vector. Theorem $5.1$ is based on Lemma $5.1$, and Lemma $5.1$ is based on Lemmas $A.2-A.10$. The proof of Lemmas $A.2-A.7$ will be the same whenever $\beta^{0}$ is a vector or a function. Lemmas $1-3$ above are the functional version of Lemmas $A.8-A.10$. We need to figure out the constraint on the smoothing parameter $h_{n}$. The constraints are

• •

Lemma $A.2-A.4$

[TABLE]

• •

Lemma $A.5-A.7$

[TABLE]

• •

Lemma $A.8-A.10$

[TABLE]

We require $nh_{n}^{2}\epsilon_{0n}^{2}M_{n}^{-2}\rightarrow\infty$. Since $\epsilon_{0n}\rightarrow 0$ and $M_{n}^{-2}\rightarrow 0$, we need to have $nh_{n}^{2}\rightarrow\infty$. Following the same argument, we need $nh_{n}^{4}\rightarrow\infty$ and $nh_{n}^{6}\rightarrow\infty$. ∎

To prove the convergence rate of the functional single index model, we need to find a convergence rate for $\hat{\beta}$, $\hat{\bm{b}}$ or $\int X\hat{\beta}$.

• Theorem5.

Suppose that $nh_{n}^{6}\rightarrow\infty$, $nh_{n}^{8}\rightarrow 0$ and $\frac{nh_{n}^{3+\frac{3}{m-1}}}{-\log h_{n}}\rightarrow\infty$, then we have

[TABLE]

for $i=1,\cdots,n$.

Proof.

For each $p_{n}$ such that $p_{n}\rightarrow\infty$ as $n\rightarrow\infty$, define

[TABLE]

Since the kernel function $K$ satisfies the Lipschitz condition, and the estimated link function $g$ is

[TABLE]

the estimated $\hat{g}$ also satisfies the Lipschitz condition. Since $\sum\limits_{j=p_{n}+1}^{\infty}c_{ij}b_{j}=O\left(p_{n}^{-\lambda}\right)$, we have

[TABLE]

as $n\rightarrow\infty$, where $K_{2}$ is a constant. Therefore, $\hat{g}\left(\sum\limits_{j=1}^{p_{n}}c_{ij}b_{j}\right)\rightarrow\hat{g}\left(\sum\limits_{j=1}^{\infty}c_{ij}b_{j}\right)$. By the definition of $J_{n,p_{n}}\left(\beta\right)$, we can get $J_{n,p_{n}}\left(\beta\right)\rightarrow J_{n}\left(\beta\right)$ as $n\rightarrow\infty$, for any $\beta\in\mathcal{B}$.

Hence, by Lemma $5.4$ in Ichimura (1993), we have

[TABLE]

where $\bm{0}_{p_{n}}$ is a $p_{n}$-dimensional mean vector, and $\mathbb{V}^{-1}\mathbb{\Sigma}\mathbb{V}^{-1}$ is a $\left(p_{n}\times p_{n}\right)$-dimensional covariance matrix. Suppose that $\mathcal{X}$ is the $\sigma$-algebra generated by $\left(X_{1},\cdots,X_{n}\right)$, the $\left(k,m\right)$th-term of the matrix $\mathbb{V}$ is

[TABLE]

For any $i=1,\cdots,n$, we have

[TABLE]

where the second equality is the bias property of the kernel density estimate. We can calculate

[TABLE]

where the first term converges to [math] by the equation $\left(\ref{kk}\right)$, since $h_{n}\rightarrow 0$ as $n\rightarrow\infty$, and the second term converges to [math] by $\left(\ref{ye}\right)$. Therefore, we have

[TABLE]

as $n\rightarrow\infty$. By the Slutsky’s Theorem, we can get

[TABLE]

Therefore,

[TABLE]

For any $\beta\in\mathcal{B}$ and any $j\in\left\{1,\cdots,n\right\}$, we have

[TABLE]

where $u=\int X_{j}\beta$.

For any $i=1,\cdots,n$, denote $c_{\cdot\cdot}=c_{ik}$, for any $k\in\left\{1,\cdots,n\right\}$, and $c_{\cdot k}=c_{ik}$ for a fix $k$. For any $k,m\in\left\{1,\cdots,n\right\}$, we have

[TABLE]

since $X_{i}$ are independent and identically distributed with $\mathrm{E}\left(c_{ij}\right)=0$, for any $i=1,\cdots,n$ and $j=1,\cdots$.

The $\left(k,m\right)$th-term of the matrix $\mathbb{\Sigma}$ is

[TABLE]

The diagonal term of the covariance matrix $\mathbb{V}^{-1}\mathbb{\Sigma}\mathbb{V}^{-1}$ is

[TABLE]

Define the truncated version of $X_{i}\left(t\right)$ and $\beta\left(t\right)$ as

[TABLE]

where $i=1,\cdots,n$.

We have for any $i=1,\cdots,n$, define $\bm{c}_{i,p_{n}}=\left(c_{i1},\cdots,c_{ip_{n}}\right)$,

[TABLE]

Therefore, $\sqrt{n}\left(\int X_{i,p_{n}}\hat{\beta}_{p_{n}}-\int X_{i,p_{n}}\beta^{0}_{p_{n}}\right)$ converges to a normal distribution with the covariance matrix $\text{var}\left(c_{\cdot\cdot}\right)\cdot\frac{p_{n}^{2}h_{n}^{2}}{\text{var}\left(c_{\cdot\cdot}\right)}=p_{n}^{2}h_{n}^{2}$. By Assumption $5$, since $p_{n}=o_{p}\left(\frac{1}{h_{n}}\right)$, we can get $\text{var}\left(c_{\cdot\cdot}\right)\cdot\frac{p_{n}^{2}h_{n}^{2}}{\text{var}\left(c_{\cdot\cdot}\right)}=o_{p}\left(1\right)$. When $n\rightarrow\infty$, $p_{n}\rightarrow\infty$, we have

[TABLE]

for $i=1,\cdots,n$. ∎

4.2 Main Theorem

• Theorem6.

If $nh_{n}^{6}\rightarrow\infty$, $nh_{n}^{8}\rightarrow 0$ and $\frac{nh_{n}^{3+\frac{3}{m-1}}}{-\log h_{n}}\rightarrow\infty$, we have

[TABLE]

Proof.

Since $X_{i}$ are independent and identically distributed, for any $j\in\left\{1,\cdots,n\right\}$, we only need to find the convergence rate of $\mathrm{E}\left[\hat{g}^{\prime\prime}\left(\int X_{j}\hat{\beta}\right)-g^{\prime\prime}\left(\int X_{j}\beta^{0}\right)\right]^{2}$. We can decompose it into three terms:

[TABLE]

where

[TABLE]

and

[TABLE]

By Lemma $7$, $8$ and $9$ in the Appendix, we have

[TABLE]

Combining these three terms, we obtain

[TABLE]

∎

5 Practical Implementation

5.1 Initialization

In $\left(\ref{step1}\right)$, the coefficient function $\beta^{0}$ is estimated by minimizing the mean square error of $g$, which is a nonlinear optimization problem. The coefficient function is approximated by a $K$-dimensional Fourier basis, as

[TABLE]

where $\bm{c}$ is a $K$-dimensional column vector, and $t\in\left[0,1\right]$. In order to ensure identifiability, we constraint the coefficient function $\beta^{0}$ to have $\int\beta^{0}=0$ and $\|\beta^{0}\|_{2}=1$. Since the Fourier basis is an orthonormal basis, a constraint on $\beta^{0}$ is equivalent to a constraint on the coefficient vector, such that $\left\|\bm{c}\right\|_{2}=1$. The first constraint can be enforced by dropping any constant terms from the Fourier basis. To compensate for rescaling $\bm{c}$, we also rescale the bandwidth $h$ such that $h=h\left\|\bm{c}\right\|$ in the optimization step.

We use R function optim to minimize the MSE, and an initial value of the coefficient vector $\bm{c}$ is needed, denoted as $\bm{c}_{\text{init}}$. We use three different methods to select the initial value.

Assume each item in $\bm{c}_{\text{init}}$ is equal, such that $\bm{c}_{\text{init}}=\left(\frac{1}{\sqrt{K}},\cdots,\frac{1}{\sqrt{K}}\right)$. We don’t have any previous knowledge about the start point, so for simplicity, choose a vector while all items are equal to each other. 2. 2.

Assume $g\left(s\right)=s$, the coefficient vector $\bm{c}_{\text{init}}$ is estimated by minimizing the ordinary least square, and normalize it to be $\left\|\bm{c}_{\text{init}}\right\|_{2}=1$. To obtain an initial for $\beta^{0}$, we need to specify a structure for the link function $g$. A linear structure of $g$ is obviously the simplest, and could be calculated easily. 3. 3.

Generate $1000$ different standard normal distribution random initial vectors (with length $\bm{c}$), and select the best $10$ initial vectors by penalized mean squared error, when the bandwidth is equal to the mean of bandwidth sequence. For each bandwidth $h$, select the initial that minimizes the penalized MSE. That means that in the cross-validation step, we select different initial for different bandwidth $h$. The selected initial is determined by the bandwidth and corresponding penalized MSE.

5.2 Cross-validation

We need to select the bandwidth $h_{n}$ in the kernel density estimation of $g$ and $g^{\prime\prime}$. We examine two cross-validation methods:

$10$-fold cross-validation. We partition the dataset into $10$ subsamples. Each time, we use $9$ subsamples as the training set and the remaining subsample as the validation set. We will observe that the $10$-fold cross validation method produces similar results with the GCV method introduced below. 2. 2.

Fix $j\in\left\{1,\cdots,n\right\}$, define

[TABLE]

where $S_{0}\left(\beta;j\right)$ is a $n$-dimensional row vector. Denote a $\left(n\times n\right)$-dimensional smoother matrix $\mathbb{S}_{h_{n}}=\left(S_{0}\left(\hat{\beta};1\right),\cdots,S_{0}\left(\hat{\beta};n\right)\right)^{\top}$. We can get $\hat{\bm{Y}}=\mathbb{S}_{h_{n}}\bm{Y}$, where $\hat{\bm{Y}}\doteq\left(\hat{g}\left(\int X_{1}\hat{\beta}\right),\cdots,\hat{g}\left(\int X_{n}\hat{\beta}\right)\right)^{\top}$. The generalized cross-validation criterion is

[TABLE]

where $\mathbb{I}$ is a $n$-dimensional identity matrix. Find a bandwidth $h_{n}$ minimizing the $\text{GCV}\left(h_{n}\right)$. Note that this does not account for estimating $\beta^{0}$.

5.3 Simulation Study

In order to obtain simulated functional data, we defined a $25$-dimensional Fourier basis $\bm{\psi}\left(t\right)$, where $t\in\left[0,1\right]$. The covariate function $X\left(t\right)$ is defined in $\left[0,1\right]$, such that

[TABLE]

where $\eta_{i}=\frac{i-1}{24}\cdot\mathrm{N}\left(0,1\right)$. The coefficient function is

[TABLE]

where $\bm{a}=\left(0,1,1,0.5,0,\cdots,0\right)^{\top}$. We use three different link functions:

$g\left(s\right)=e^{-s}$. 2. 2.

$g\left(s\right)=-s^{2}$. 3. 3.

$g\left(s\right)=s$.

In order to measure the performance of our estimators, we define the MSE of the estimated $\beta^{0}$ and $g^{\left(k\right)}$ to be

[TABLE]

and

[TABLE]

where $\hat{Y}^{\left(k\right)}_{i}=\hat{g}^{\left(k\right)}\left(\int_{0}^{1}X_{i}\left(t\right)\hat{\beta}\left(t\right)\mathrm{d}t\right)$.

Theorem $6$ shows that for consistency, the bandwidth must scale between $O\left(n^{-\frac{1}{6}}\right)$ and $O\left(n^{-\frac{1}{8}}\right)$. Assume that the optimal bandwidth for the curvature $g^{\prime\prime}$ is $n^{-\frac{1}{7}}$, by Chen et al. (2011), we can expect that the optimal bandwidth for the link function $g$ is $n^{-\frac{1}{5}}$. We use either the GCV or the $10$-fold cross-validation to select a bandwidth h, and rescale it to be $h^{\frac{5}{7}}$. As we have discussed before, we constrain the coefficient $\left\|\beta^{0}\right\|_{2}=1$; this was achieved by rescaling the bandwidth in our objective function. After solving the nonlinear optimization problem, we rescale both $\bm{c}$ and the bandwidth. The following table shows the simulation results of the GCV method. The number of data points we use are $n=100$ and $n=1000$. In Table $2$, we also show the RASE$2$ results without rescaling the bandwidth $h$ to be $h\left\|\bm{c}\right\|$ and $h^{\frac{5}{7}}$, if we start from random initials. We can conclude that re-scaling does matter to the final results and it reduces the error of RASE$2$. The $10$-fold cross-validation results in the Appendix C also confirm that by re-scaling the bandwidth, we improve our estimate for $g^{\prime\prime}$.

We observe that RSE and RASE, as expected, achieve best performance when we initialize our optimizer at the true values. However, the more natural Linear initialization strategy does not outperform initializing at Equal coefficients. For both RASE and RASE$2$, a more intensive search over initializations pays off; in the case of RASE$2$ this even outperforms starting from true values. We suspect that this is associated with differing optimal smoothness criteria. Table $2$ compares our results to when we do not use the $h^{\frac{5}{7}}$ re-scaling where applying this has a significant effect; by starting from a position far from the optimum, a large bandwidth may have the effect of smoothing the objective function, indirectly improving our re-estimate of $g"$.

Overall, the random initial produces a much better results compared to other initial strategies. Starting from $1000$ initial vectors, we have a chance to select the best $10$ initial vectors, which increases the probability of selecting a “good” starting points and decreases the chance of converging to a local minimum for the non-linear optimization problem, although this comes at a significant computational prices. In addition, we do observe an improvement of RASE and RASE$2$ when $n=1000$ compared to $n=100$.

In order to provide a visual sense of the performance of our estimate, in Figures $1$ and $2$, we plot the estimates of the link function $g\left(s\right)=e^{-s}$ and $g\left(s\right)=s$, respectively. We observe that the estimate and the true curve of the link function $g$ almost overlap with each other, but the second derivative has significantly larger error relative to the truth.

6 Ecological Data

6.1 Model Formulation

Examining our ecological questions, the purpose of estimating the second derivative of the link function in a functional single index model is to figure out whether the link function $g$ is convex or concave. Then, we can answer the question: in which environment, constant or varying, the plant will grow better. We apply our nested estimation method to plant growth dataset. In this dataset, there are several variables:

$\tt{logarea.t1,logarea.t0:}$ the plant’s logarithm of area at time $t0$ and $t1$, where $t0$ is the observation start time and $t1$ is the end time. A relatively large quantity indicates a high growth rate of the plant at that time. 2. 2.

$\tt{W}$: a measure of plant competition. Taken to be a scalar covariate. 3. 3.

$\tt{p.00-p.36}$: discrete aggregated temporal record of precipitation, denoted as $p\left(s\right)$. 4. 4.

$\tt{t.00-t.36}$: discrete aggregated temporal record of temperature, denoted as $t\left(s\right)$.

The precipitation and temperature histories are modeled as two covariate functions. Assume that the response variable is $logarea.t1-logarea.t0$, a Functional Single Index model is:

[TABLE]

where the coefficient $\alpha$, the functions $g$, $\beta_{1}$ and $\beta_{2}$ need to be estimated.

6.2 Results

We used three different starting values: linear, equal and random. For each starting point, we selected the curve with minimum GCV value. Then, we selected the estimate with the minimum GCV value among all starting points. The random starting point is selected. The plot of the estimated $g$, $g^{\prime\prime}$, and the coefficient functions $\beta_{1}$ and $\beta_{2}$ is in the Figure $3$. Since the estimated $g^{\prime\prime}$ is always negative, the link function $g$ is concave. We could conclude that the species will grow better with a constant environment.

7 Conclusion

To answer the ecological question, we need to figure out the convexity or concavity of the link function $g$, or equivalently, find whether the second derivative is positive or negative. In this paper, we used the local quadratic method to approximate the link function $g$, and estimated the curvature of $g$ and the coefficient function $\beta^{0}$ by a nested optimization procedure. Under some assumptions, we showed that the coefficient function $\beta^{0}$ could be estimated root-$n$ consistently. In addition, the rate of convergence of the curvature $g^{\prime\prime}$ is $\frac{1}{n}\sum\limits_{i=1}^{n}\mathrm{E}\left[\hat{g}^{\prime\prime}\left(\int X_{i}\hat{\beta}\right)-g^{\prime\prime}\left(\int X_{i}\beta^{0}\right)\right]^{2}=O\left(h_{n}^{4}+\frac{1}{nh_{n}^{4}}\right)$.

In the simulation study, we used three different link functions, convex, concave and neither convex nor concave. While we derive convergence rates for the curvature of $g$, our simulation results demonstrate the numerical challenges that accompany Functional Single Index models. We can estimate $g$ fairly well, but our estimates of $g^{\prime\prime}$ are sensitive to the choice of initial condition, requiring considerable care in optimization. We expect that these numerical challenges are specific to the estimators employed, but they suggest that alternative means for the influence of environmental variability on plant growth are warranted.

Appendix A Bernstein’s Inequality

Bernstein’s Inequality. Let $Y_{1n},\cdots,Y_{nn}$ be independent random variables with means [math] and bounded ranges, that is $\left|Y_{in}\right|\leq c_{n}$. Write $\sigma_{in}^{2}$ for the variance of $Y_{in}$. Suppose $V_{n}\geq\sigma_{1n}^{2}+\cdots+\sigma_{nn}^{2}$. Then for each $\eta_{n}>0$,

[TABLE]

Appendix B Lemmas for Theorem $6$

• Lemma7.

If $nh_{n}^{6}\rightarrow\infty$, $nh_{n}^{8}\rightarrow 0$ and $\frac{nh_{n}^{3+\frac{3}{m-1}}}{-\log h_{n}}\rightarrow\infty$, then we have

[TABLE]

Proof.

Denote $\bm{\epsilon}=\left(\epsilon_{1},\cdots,\epsilon_{n}\right)^{\top}$. The term could be expressed as

[TABLE]

We bound the next two terms as:

[TABLE]

and

[TABLE]

where $M_{1}$, $M_{2}$ are constants and $u=\int X_{j}\hat{\beta}$.

Therefore, we have

[TABLE]

∎

• Lemma8.

If $nh_{n}^{6}\rightarrow\infty$, $nh_{n}^{8}\rightarrow 0$ and $\frac{nh_{n}^{3+\frac{3}{m-1}}}{-\log h_{n}}\rightarrow\infty$, then we have

[TABLE]

Proof.

Calculate

[TABLE]

Write

[TABLE]

for $i=1,\cdots,n$, we decompose it into two terms $A$ and $B$, where $A$ is

[TABLE]

and $B$ is

[TABLE]

Denote $K_{ij}\doteq K\left(\frac{\int X_{i}\hat{\beta}-\int X_{j}\hat{\beta}}{h_{n}}\right)$, the term $A$ can be calculated as

[TABLE]

where

[TABLE]

and

[TABLE]

Therefore, the term $A$ is

[TABLE]

The term $B$ is bounded by

[TABLE]

where $E_{1}$ is a constant.

Combining the terms $A$ and $B$, we have

[TABLE]

∎

• Lemma9.

If $nh_{n}^{6}\rightarrow\infty$, $nh_{n}^{8}\rightarrow 0$ and $\frac{nh_{n}^{3+\frac{3}{m-1}}}{-\log h_{n}}\rightarrow\infty$, then we have

[TABLE]

Proof.

By Theorem $5$, we have $\left(\int X_{i}\hat{\beta}-\int X_{i}\beta^{0}\right)=O\left(\frac{1}{\sqrt{n}}\right)$. Since $g^{\prime\prime}$ satisfies the Lipschitz condition, we can get

[TABLE]

where $K_{3}$ is a constant. ∎

Appendix C $10$-fold CV Results

The simulation results using the $10$-fold cross-validation is:

For comparison purpose, we estimate the curvature with or without rescaling the bandwidth, where the results are in Table $4$.

Appendix D True Starting Value

In Table $5$ and $6$, we compare the RASE$2$ results of true starting value with rescaled and original bandwidths.

Appendix E CV Values

The CV values for both GCV and $10$-fold cross-validation are in Table $7$.