Estimating a Function and Its Derivatives Under a Smoothness Condition

1. Introduction

We are concerned with the problem of estimating an unknown function $f_{*}:{\mathbb{R}}^d\to \mathbb{R}$ and its partial derivatives over a domain ${[a,b]}^d$ of interest when there is no closed-form formula for $f_{*}$; hence, simulation must be conducted to estimate $f_{*}(\mathbf{x})$ at $\mathbf{x}\in {\mathbb{R}}^d$, or only noisy data on $f_{*}$ are available. We make no assumptions about $f_{*}$ except that it has square integrable partial derivatives of order m.

To place the problem in a more rigorous setting, we consider the situation where we wish to estimate the unknown function $f_{*}:{\mathbb{R}}^d\to \mathbb{R}$ and its partial derivatives by taking observations $(X_1,Y_1),\cdots ,(X_n,$ $Y_n)$ satisfying

for $i=1,2,\cdots ,n$, where $((X_i,{\epsilon }_i):1\le i\le n)$ is a sequence of ${[a,b]}^d\times \mathbb{R}$-valued independent and identically distributed (iid) random vectors satisfying $\mathbb{E}({\epsilon }_i|X_i)=0$ and $\mathbb{E}({\epsilon }_i^2|X_i)={\sigma }^2<\infty$ for $i=1,2,\cdots ,n$.

Under the assumption that $f_{*}$ has square integrable partial derivatives of order m with $m\ge 1$, our goal is to estimate $f_{*}$ and its partial derivatives of order up to $m-1$ as functions over ${[a,b]}^d$ from the data set $(X_1,Y_1),\cdots ,(X_n,Y_n)$.

When the only information available on the unknown regression function $f_{*}$ is the fact that it is “smooth” so that it has square integrable partial derivatives of order m, one way to estimate $f_{*}$ is to fit a “smooth” function to the data set, that is, find the solution to the following problem:

over $f\in \mathcal{F}$ for some appropriately chosen set $\mathcal{F}$ of functions. The next question is “Which set is appropriate for $\mathcal{F}$?” One natural candidate for $\mathcal{F}$ is the set ${\mathcal{F}}^{\prime }$ of functions $f:{\mathbb{R}}^d\to \mathbb{R}$ whose partial derivatives of order m are square integrable, that is, ${\mathcal{F}}^{\prime }=\{f:{\mathbb{R}}^d\to \mathbb{R}:J(f)<\infty \}$, where

$\alpha =({\alpha }_1,\cdots ,{\alpha }_d)$ is a d-tuple of nonnegative integers, $|\alpha |$ is the order of $\alpha$ defined by $|\alpha |={\alpha }_1+\cdots +{\alpha }_d$, and $D^{\alpha }$ denotes the partial differential operator defined by

However, a question arises as to whether there exists a solution to (1). To guarantee the existence of a solution to (1), one of the properties that $\mathcal{F}$ must have is the completeness with respect to a suitably chosen metric (or a pseudometric). It turns out that $\mathcal{F}$ must be larger than ${\mathcal{F}}^{\prime }$ to be complete. However, if we let $\mathcal{F}$ be the set of “generalized functions” whose “weak” derivatives of order m are square integrable, then $\mathcal{F}$ contains ${\mathcal{F}}^{\prime }$ and is large enough to be a semi-Hilbert space in the $\Vert ·{\Vert }_m$ seminorm, where $\Vert f{\Vert }_m\triangleq \sqrt{J(f)}$; see, for example, Meinguet [51, theorem 1 on p. 130]. (A weak derivative of a function f is a generalized version of the classical derivative in the sense that whenever the classical derivative $D^{\alpha }f$ exists, it is also a weak derivative of f.) Hence, throughout this paper, we will assume $\mathcal{F}={\mathcal{F}}_m$, where ${\mathcal{F}}_m$ is the space of generalized functions whose weak partial derivatives of order m are square integrable; see Section 2 of this paper and Oden and Reddy [54] for the precise definitions of ${\mathcal{F}}_m$, generalized functions, and weak derivatives. We will further restrict our attention to the functions f whose smoothness $J(f)$ is bounded by a certain constant. Thus, a natural candidate for the estimator of $f_{*}$ is the solution ${\widehat{f}}_n$ to the following problem:

for some upper bound $U_n$ on $J(f)$. As Proposition 1 of this paper suggests, the solution to Problem (A) exists under the assumption that $2m>d$ and $\{X_1,\cdots ,X_n\}$ are ${\mathcal{P}}_{m-1}$-unisolvent. For definitions, see Assumption 1 in Section 3.

The question now becomes how to compute the solution to Problem (A) numerically. As Proposition 6 of this paper states, the solution to Problem (A) can be found by solving a convex programming problem. However, it is not clear how to compute or estimate $U_n$ from the data set. This necessitates another approach that can estimate $f_{*}$ using a parameter that can be directly estimated from the data set. As an alternative to Problem (A), we consider the solution ${\widehat{g}}_n$ to the following problem:

for some upper bound $S_n$ on the average of squared errors. Problem (B) is appealing from a computational point of view; it is intuitively acceptable that $S_n$ should be close to ${\sigma }^2$, and ${\sigma }^2$ can be readily estimated from the data set.

In the context of nonparametric regression techniques in the statistics literature, Problems (A) and (B) are often contrasted with the following formulation:

where $\lambda$ is a constant called the smoothing parameter. The solution to Problem (C) is often referred to as the “penalized least squares estimator”; see Györfi et al. [27], Wahba [71], and the references therein. One weakness of the penalized least squares estimator is that its performance is highly sensitive to the choice of $\lambda$, and the question of how to choose the smoothing parameter $\lambda$ has been a challenging one. Cross validation (Craven and Wahba [15]) has been a popular method for choosing the smoothing parameter, but it does not guarantee the consistency of the estimators of the derivatives as $n\to \infty$. Thus, when the primary purpose of a modeler is to estimate a derivative of $f_{*}$ and he or she does not want to use any parameter that cannot be directly estimated from the data set, Problem (B) is preferred from the computational point of view.

This paper is motivated by a stock trader whose goal is to estimate “gamma,” which is the second derivative of the value, $f_{*}$, of a stock option as a function of the underlying stock price. The only information that he has is the fact that gamma is a smooth function of the stock price. The trader can use Problem (C) to estimate $f_{*}$ and its second derivative, but he needs to estimate the smoothing parameter $\lambda$. Cross validation does not guarantee the consistency of the estimator of the second derivative of $f_{*}$ as $n\to \infty$. Moreover, the trader does not want to use any parameter that cannot be directly estimated from a data set. In this situation, Problem (B) better serves his purpose because its only parameter $S_n$ can be directly estimated from the data set. In fact, Problem (B) has been used widely within the numerical analysis community; see, for example, Cox [13, p. 530].

Despite its popularity in the numerical analysis community, little is known about how to compute ${\widehat{f}}_n$ and ${\widehat{g}}_n$ numerically and what statistical properties ${\widehat{f}}_n$ and ${\widehat{g}}_n$ possess. The contributions of this paper are showing that ${\widehat{f}}_n$ and ${\widehat{g}}_n$ are computable as solutions to convex programs, establishing the consistency of ${\widehat{f}}_n$ and ${\widehat{g}}_n$ as estimators of $f_{*}$, establishing the consistency of $D^{\alpha }{\widehat{f}}_n$ and $D^{\alpha }{\widehat{g}}_n$ as estimators of $D^{\alpha }{f}_{*}$ for $|\alpha |=1,\cdots ,m-1$, and computing the convergence rate of ${\widehat{f}}_n$ when $n\to \infty$. The key contributions of this paper can be summarized as follows.

1. We identify the relationship between Problems (A) and (B). We prove that for $U_n$ within a certain interval, the solution ${\widehat{f}}_n$ to Problem (A) exists uniquely and becomes a unique solution to Problem (B) with $S_n=J({\widehat{f}}_n)$. Conversely, for $S_n$ within a certain interval, the solution ${\widehat{g}}_n$ to Problem (B) exists uniquely and becomes a unique solution to Problem (A) with $U_n=E_n({\widehat{g}}_n)$.

2. We discuss the computational aspects of Problems (A) and (B). We show that the solutions to Problems (A) and (B) can be found by solving convex programming problems. There exist numerous efficient algorithms that can successfully find the solutions to convex programs (see, for example, Boyd and Vandenberghe [6], Zangwill [76]), so this enables us to easily compute the solutions to Problems (A) and (B).

3. We establish the consistency of ${\widehat{f}}_n$ and ${\widehat{g}}_n$ as estimators of $f_{*}$ and the consistency of $D^{\alpha }{\widehat{f}}_n$ and $D^{\alpha }{\widehat{g}}_n$ as estimators of $D^{\alpha }{f}_{*}$ for $|\alpha |=1,\cdots ,m-1$ as $n\to \infty$. Our main results for Problem (A) state that if $U_n\ge J(f_{*})$ for n sufficiently large and ${\text{lim sup}}_{n\to \infty }{U}_n<\infty$, we have

   $$\begin{array}{l}{\sup }_{\mathbf{x}\in {[a,b]}^d}|{\widehat{f}}_n(\mathbf{x})-f_{*}(\mathbf{x})|\to 0\hspace{1em}\text{almost} \text{surely} (\text{a}.\text{s}.),\hfill \end{array}$$
   and we have
   $$\begin{array}{l}\mathbb{E}{({\widehat{f}}_n(X)-f_{*}(X))}^2\to 0\hspace{0.5em}\text{and}\hspace{0.5em}\mathbb{E}{(D^{\alpha }{\widehat{f}}_n(X)-D^{\alpha }{f}_{*}(X))}^2\to 0\hfill \end{array}$$
   for $|\alpha |=1,\cdots ,m-1$ as $n\to \infty$ under modest assumptions. On the other hand, our main results for Problem (B) state that under some conditions on ${\widehat{g}}_n$ and $S_n$, we have
   $$\begin{array}{l}\underset{\mathbf{x}\in {[a,b]}^d}{\sup }|{\widehat{g}}_n(\mathbf{x})-f_{*}(\mathbf{x})|\to 0\hspace{1em}\text{a}.\text{s}.,\hfill \end{array}$$
   $$\begin{array}{l}\mathbb{E}{({\widehat{g}}_n(X)-f_{*}(X))}^2\to 0\hspace{1em} \text{and} \hspace{1em}\mathbb{E}{(D^{\alpha }{\widehat{g}}_n(X)-D^{\alpha }{f}_{*}(X))}^2\to 0\hfill \end{array}$$
   for $|\alpha |=1,\cdots ,m-1$ as $n\to \infty$.

4. We compute the convergence rate of ${\widehat{f}}_n$ as $n\to \infty$. Specifically, when $U_n\ge J(f_{*})$ for n sufficiently large and ${\text{lim sup}}_{n\to \infty }{U}_n<\infty$, we have

   $$\begin{array}{l}{\left\{\frac{1}{n}{\displaystyle \sum _{i=1}^n{({\widehat{f}}_n(X_i)-f_{*}(X_i))}^2}\right\}}^{1/2}={\mathcal{O}}_p(n^{-m/(2m+d)})\hfill \end{array}$$
   under the assumption that the ${\epsilon }_i$’s are uniformly sub-Gaussian.

1.1. Literature Review

This paper is a contribution to the literature of nonparametric regression for smooth functions. Problem (C) has been studied extensively in the literature, whereas relatively less attention has been paid to Problems (A) and (B). For Problem (C), the existence of a solution is established by Duchon [19], and the convergence rates are derived by Cox [14], Rice and Rosenblatt [57], Utreras [67], and van de Geer [69]. Craven and Wahba [15] and Utreras [66] discuss how to choose the smoothing parameter through cross validation. Duchon [19] has shown that the solution to Problem (C) has an explicit representation as a sum of basis functions. However, this approach often requires solving a system of linear equations that involves a dense and highly ill-conditioned matrix (see Dierckx [17, p. 145]). Numerically more efficient algorithms are developed by Hutchinson and de Hoog [33] and Luo and Wahba [47]. Schoenberg [60] has shown that a solution to Problem (A) is also a solution to Problem (C) for some $\lambda$ appropriately chosen, whereas Kersey [35] and Reinsch [56] studied the relationships between Problems (C) and (B). Green et al. [24] computed convergence rates of the Laplacian smoothing estimator, which can be viewed as a discrete approximation of Problem (C). For further references, see, for example, de Boor [16], Dierckx [17], Green and Silverman [23], Györfi et al. [27], Wahba [71], and Wegman and Wright [74]. For Problem (A), when $d=1$, the existence and uniqueness of the solution are established by Schoenberg [60], and the convergence rates are derived by van de Geer [68]. When $d\ge 2$, we could not find any work concerning how to compute the solution to Problem (A) numerically or what statistical properties the solution to Problem (A) possesses. For Problem (B), when $d=1$, the existence and uniqueness of the solution are established by Reinsch [56]. However, we could not find any work concerning how to compute the solution to Problem (B) and what statistical properties the solution to Problem (B) possesses either when $d=1$ or when $d\ge 2$. Therefore, the main contribution of this paper is to establish a theoretical foundation of Problems (A) and (B) for the case when $d\ge 1$ by suggesting a convex programming formulation, establishing consistency, and computing the convergence rate.

This paper is also a contribution to the literature on shape constrained estimation. Many researchers have studied the problem of estimating the unknown regression function $f_{*}$ when $f_{*}$ is known to possess shape properties, such as monotonicity or convexity. One of the most popular approaches is least squares estimation under shape restrictions. For instance, when $f_{*}$ is known to be monotone, one can fit a monotone function to the data set by solving the following quadratic program:

where $X_i\le X_j$ denotes coordinate-wise monotonicity. The solution to (2) is referred to as the isotonic regression estimator and has well-established theory. Brunk [7] and Brunk [8] introduced isotonic regression, established consistency, and computed convergence rates when $d=1$. Mammen [48] considered estimation of a smooth monotone function when $d=1$. Lim [44] established consistency of isotonic regression and computed its convergence rate when $d\ge 1$ and $f_{*}$ is possibly misspecified. On the other hand, when $f_{*}$ is known to be convex, the constraint of (2) can be replaced by the convexity condition, and the solution to the corresponding quadratic program is referred to as the convex regression estimator. Hildreth [30] first introduced convex regression when $d=1$; Hanson and Pledger [28] established consistency when $d=1$; Groeneboom et al. [26] computed the rate of convergence when $d=1$; Kuosmanen [40] formulated the convex regression estimator as the solution to a quadratic program when $d\ge 1$; Lim and Glynn [46] and Seijo and Sen [61] established consistency of the multivariate convex regression estimator; Lim [44] established consistency of convex regression and computed its convergence rates when $d\ge 1$ and $f_{*}$ is possibly misspecified; Bertsimas and Mundru [4] and Lee et al. [43] proposed computationally efficient algorithms for estimating a concave monotone function and a convex function, respectively; Lim [44] and Yagi et al. [75] proposed hypothesis tests for detecting monotonicity and convexity; Kuosmanen and Johnson [42] discussed an interesting application of production function estimation and proposed a shape constrained estimator for such an application; and Kuosmanen and Johnson [41] established connection between data envelopment analysis and least squares estimation under shape restrictions. Two drawbacks of the convex regression estimator are that (1) it becomes computationally expensive when n increases and that (2) it tends to overfit the data near the boundary of the domain; see, for example, Lim [45, p. 70]. To overcome these shortcomings, Yagi et al. [75] proposed a local polynomial kernel estimator with shape constraints, and Keshvari [36] fitted a convex piecewise linear function to data where the number of linear segments is prespecified. An advantage of Problem (B) is that one can readily compute an estimate of $S_n$ using sample variances. In the context of optimizing behavior in economics, Varian [70] similarly used sample variances to set bounds for acceptable deviations from the fitted function. For a comprehensive survey, see Johnson and Jiang [34] and the references therein. Although the previous literature has focused on monotonicity and convexity constraints, this paper demonstrates that one can work with the smoothness condition only. In particular, Problem (A) replaces the constraint of (2) by a condition on the degree of smoothness (i.e., $J(f)\le U_n$). Therefore, this paper can be viewed as an addition to the literature on nonparametric regression subject to shape constraints.

This paper can be viewed as a contribution to the literature on metamodel estimation or response surface estimation, which has been an important subject of research in the simulation community. A common goal in this setting is to fit a metamodel to a simulation output, where the simulation runs are time consuming. Therefore, most work has focused on a setting where n is relatively small. One of the earlier works is the response surface methodology that is based on the idea of fitting a polynomial function to the data set locally; see, for example, Myers and Montgomery [52]. When the $Y_i$’s are possibly correlated, the kriging-based methods serve as good alternatives; see, for example, Ankenman et al. [3]. The kriging-based methods, in general, do not assume that the regression function $f_{*}$ is m times differentiable, where $m\ge 2$. Hence, they do not produce any estimator of the derivatives of $f_{*}$. Rather, the kriging estimator of $f_{*}(\mathbf{x})$ at any $\mathbf{x}\in {[a,b]}^d$ is expressed as the weighted sum of the $Y_i$’s. The weights are dependent on ${\mathit{x}}_0$ and calculated in a way in which the mean squared deviations are minimized. In Chen et al. [9], under the assumption that the regression function is differentiable, the gradient estimates are observed along with the $Y_i$’s, and they are used to obtain an estimator of $f_{*}$ that shows superior performance to the stochastic kriging estimator (see Ankenman et al. [3]). In the radial basis function estimators, the regression function is approximated by a linear combination of radially symmetric functions, where the coefficients and other parameters are estimated from the data set; see, for example, Franke [21]. The orthogonal series estimators or wavelet estimators represent the regression function by its Fourier series and estimate the coefficients from the data set; see, for example, Donoho and Johnstone [18]. Nonparametric regression methods, such as the kernel regression method (Nadaraya [53], Watson [73]), the local polynomial regression estimators (Cleveland [12]), the smoothing splines (Reinsch [55]), the weighted least squares regression estimators (Ruppert and Wand [58], Salemi et al. [59]), and the nearest neighbor estimators (Shapiro [62], Stone [64]), have received considerable attention in the statistics literature; see, for example, Eubank [20], Green and Silverman [23], Härdle [29], and Wand and Jones [72] for comprehensive surveys. Most of these estimators invariably require the estimation of smoothing parameters that cannot be estimated directly from the data set. In this regard, Problem (B) has the advantage of estimating only one parameter, $S_n$, which has a direct meaning and can be estimated from the sample variances computed from the data set.

1.2. Organization of This Paper

This paper is organized as follows. Section 2 introduces some definitions and preliminaries. We precisely describe the relationship between Problems (A) and (B) in Section 3, whereas Section 4 is concerned with the convex programming formulations of Problems (A) and (B). Section 5 states the main results regarding consistency and convergence rates, and the numerical results can be found in Section 6. Section 7 provides discussions on future research topics. The proofs of all theorems, corollaries, lemmas, and propositions are provided in the Appendix.

2. Definitions and Preliminaries

For $\mathbf{x}=(x_1,\cdots ,x_d)\in {\mathbb{R}}^d$, its norm is given by $\Vert \mathbf{x}\Vert ={\{{\sum }_{j=1}^n{x}_j^2\}}^{1/2}$. For a sequence of random variables $(Z_n:n\ge 1)$ and a sequence of positive real numbers $({\alpha }_n:n\ge 1)$, we say $Z_n={\mathcal{O}}_p({\alpha }_n)$ as $n\to \infty$ if for any $ϵ>0$, there exist constants C and N such that $\mathbb{P}(|Z_n/{\alpha }_n|>C)<ϵ$ for all $n\ge N$.

When $2m>d$, we define ${\mathcal{F}}_m$ by the space of generalized functions whose weak partial derivatives of order m are square integrable, that is,

where ${\mathcal{D}}^{\prime }({\mathbb{R}}^d)$ is the space of Schwartz distributions. A test function on ${\mathbb{R}}^d$ is a real-valued, infinitely differentiable function with compact support. Let $\mathcal{D}({\mathbb{R}}^d)$ be the set of all test functions on ${\mathbb{R}}^d$. A functional on $\mathcal{D}({\mathbb{R}}^d)$ is a mapping that assigns to each $\varphi \in \mathcal{D}({\mathbb{R}}^d)$ a real number. A (Schwartz) distribution on ${\mathbb{R}}^d$ is a continuous linear functional on $\mathcal{D}({\mathbb{R}}^d)$. See Adams and Fournier [2, p. 20] or Oden and Reddy [54, p. 14] for detailed definitions. Let ${(·,·)}_m$ denote the semi-inner product on ${\mathcal{F}}_m$, defined byfor $f,g\in {\mathcal{F}}_m$. The associated seminorm $\Vert ·{\Vert }_m$ is defined by $\Vert f{\Vert }_m\triangleq {(f,f)}_m^{1/2}$. It should be noted that the kernel of this seminorm is the vector space ${\mathcal{P}}_{m-1}$ of all polynomials of total degree less than or equal to $m-1$ defined on ${\mathbb{R}}^d$, the dimension of ${\mathcal{P}}_{m-1}$ isand ${\mathcal{P}}_{m-1}$ is spanned by the M monomials of the total degree less than or equal to $m-1$. We will denote the M monomials of total degree less than or equal to $m-1$ by $p_1,\cdots ,p_M$.

It should be also noted that $\Vert ·{\Vert }_m$ is a seminorm on ${\mathcal{F}}_m$, but a norm on the equivalence classes in ${\mathcal{F}}_m/{\mathcal{P}}_{m-1}$, and ${\mathcal{F}}_m/{\mathcal{P}}_{m-1}$ is a Hilbert space under ${(·,·)}_m$; see Meinguet [51, theorem 1, p. 130] for details.

For any open-bounded subset $\mathrm{\Omega }\subset {\mathbb{R}}^d$, we define $L_2(\mathrm{\Omega })$ by the set of functions f defined on $\mathrm{\Omega }$ satisfying ${\int }_{\mathrm{\Omega }}{\{f(\mathbf{x})\}}^{2}d\mathbf{x}<\infty$. Its $L_2(\mathrm{\Omega })$ norm is defined as follows:

for $f\in L_2(\mathrm{\Omega })$. We also define $W_m(\mathrm{\Omega })$ by the set of functions whose weak partial derivatives up to order m are in $L_2(\mathrm{\Omega })$, that is,

Its $W_m(\mathrm{\Omega })$ norm is defined as follows:

for $f\in W_m(\mathrm{\Omega })$. $W_m(\mathrm{\Omega })$ is referred to as a Sobolev space of order m on $\mathrm{\Omega }$.

3. Relationships Between Problems (A) and (B)

In this section, we study the relationship between Problems (A) and (B). When $d=1$, Schoenberg [60] established some results similar to Lemma 1, Proposition 3, Proposition 4, and Proposition 5 of this paper. Thus, some of our results in this section can be viewed as extensions of those in Schoenberg [60] to the multidimensional case. Throughout this section, we assume that $(X_1,Y_1),\cdots ,(X_n,Y_n)$ are given and that the following assumption holds.

Assumption 1.

m is a positive integer satisfying $2m>d$, and $\{X_1,\cdots ,X_n\}$ is a set of mutually distinct points containing a ${\mathcal{P}}_{m-1}$-unisolvent set. A set of points $\{{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_L\}\subset {\mathbb{R}}^d$ is called ${\mathcal{P}}_{m-1}$-unisolvent if for any $p\in {\mathcal{P}}_{m-1}$, the condition $p({\mathbf{x}}_j)=0$ for all $1\le j\le L$ implies $p(\mathbf{x})=0$ for all $\mathbf{x}\in {\mathbb{R}}^d$. In particular, when $d=1$, $\{{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_L\}$ is ${\mathcal{P}}_{m-1}$-unisolvent if $L\ge m$ and ${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_L$ are mutually distinct.

Propositions 1 and 2 establish the existence of the solutions to Problems (A) and (B) under Assumption 1.

Proposition 1.

Let $U_n\ge 0$ be given. Under Assumption 1, there exists a solution to Problem (A). Furthermore, if ${\widehat{f}}_n$ and ${\tilde{f}}_n$ are two solutions to Problem (A), then ${\widehat{f}}_n(X_i)={\tilde{f}}_n(X_i)$ for $1\le i\le n$.

Proposition 2.

Let $S_n\ge 0$ be given. Under Assumption 1, there exists a solution to Problem (B). Furthermore, if ${\widehat{g}}_n$ and ${\tilde{g}}_n$ are two solutions to Problem (B), then they differ by a polynomial of total degree less than or equal to $m-1$ (i.e., ${\widehat{g}}_n-{\tilde{g}}_n\in {\mathcal{P}}_{m-1}$).

Our next propositions, Propositions 3, 4, and 5, are concerned with the uniqueness of the solutions to Problems (A) and (B). For our analysis, we need to define the following function: ${\mathrm{\Psi }}_n:[0,\infty )\to \mathbb{R}$. For any nonnegative real number u, let ${\mathrm{\Psi }}_n(u)=E_n(f^u)$, where $f^u$ is a solution to Problem (A) with $U_n=u$. The following lemma characterizes the shape of ${\mathrm{\Psi }}_n$.

Lemma 1.

There exists a minimizer, say $f_P$, of $E_n(f)$ over $f\in {\mathcal{P}}_{m-1}$. Under Assumption 1, there exists a unique solution, say $f_I$, to the following problem:

$$\begin{array}{ll}{\mathit{\text{Minimize}}}_{f\in {\mathcal{F}}_m}\hspace{0.17em}J(f)\hspace{0.17em}\hfill & \mathit{\text{\hspace{0.17em}subject}} to\hspace{0.17em}f(X_i)=Y_i\hspace{1em}1\le i\le n.\hfill \end{array}$$

Under Assumption 1, ${\mathrm{\Psi }}_n$ is nonnegative, convex, and strictly decreasing over $[0,J(f_I)]$ with ${\mathrm{\Psi }}_n(0)=E_n(f_P)$ and ${\mathrm{\Psi }}_n(x)=0$ for all $x\ge J(f_I)$.

The following propositions discuss the uniqueness of the solutions to Problems (A) and (B) and the relationship between them.

Proposition 3.

Assume Assumption 1. Let $0\le U_n<J(f_I)$ be given. Then, there exists a unique solution ${\widehat{f}}_n$ to Problem (A), and ${\widehat{f}}_n$ satisfies $J({\widehat{f}}_n)=U_n$. In other words, if ${\widehat{f}}_n$ and ${\tilde{f}}_n$ are two solutions to Problem (A), then ${\widehat{f}}_n(\mathbf{x})={\tilde{f}}_n(\mathbf{x})$ for all $\mathbf{x}\in {\mathbb{R}}^d$. Furthermore, ${\widehat{f}}_n$ is also a solution to Problem (B) with $S_n=E_n({\widehat{f}}_n)$.

Proposition 4.

Assume Assumption 1. Let $0<S_n\le E_n(f_P)$ be given. Then, there exists a unique solution ${\widehat{g}}_n$ to Problem (B), and ${\widehat{g}}_n$ satisfies $E_n({\widehat{g}}_n)=S_n$. In other words, if ${\widehat{g}}_n$ and ${\tilde{g}}_n$ are two solutions to Problem (B), then ${\widehat{g}}_n(\mathbf{x})={\tilde{g}}_n(\mathbf{x})$ for all $\mathbf{x}\in {\mathbb{R}}^d$. Furthermore, ${\widehat{g}}_n$ is also a solution to Problem (A) with $U_n=J({\widehat{g}}_n)$.

Combining Propositions 3 and 4 yields the following proposition.

Proposition 5.

Assume Assumption 1.

1. Let $0\le U_n\le J(f_I)$ be given. Then, there exists a unique solution ${\widehat{f}}_n$ to Problem (A), and ${\widehat{f}}_n$ satisfies $J({\widehat{f}}_n)=U_n$. Furthermore, ${\widehat{f}}_n$ is also a unique solution to Problem (B) with $S_n=E_n({\widehat{f}}_n)$.

2. Let $0<S_n\le E_n(f_P)$ be given. Then, there exists a unique solution ${\widehat{g}}_n$ to Problem (B), and ${\widehat{g}}_n$ satisfies $E_n({\widehat{g}}_n)=S_n$. Furthermore, ${\widehat{g}}_n$ is also a unique solution to Problem (A) with $U_n=J({\widehat{g}}_n)$.

The proofs of Lemma 1 and Propositions 1–5 are provided in Appendix A.

4. Convex Programming Formulations for Problems (A) and (B)

In this section, we discuss the question of how to compute the solutions to Problems (A) and (B) numerically. The following propositions, Propositions 6 and 7, reveal that we can find the solutions to Problems (A) and (B) by solving convex programs. Their proofs are provided in Appendix A.

We start with some definitions. Let ${\mathbf{s}}_1,\cdots ,{\mathbf{s}}_M$ be any ${\mathcal{P}}_{m-1}$-unisolvent set of M fixed points in ${\mathbb{R}}^d$. (M is given in (3).) Let $q_1,\cdots ,q_M$ be the unique polynomials of total degree less than m satisfying $q_i({\mathbf{s}}_j)=1$ if $i=j$ and $q_i({\mathbf{s}}_j)=0$ if $i\ne j$, that is,

for $\mathbf{x}\in {\mathbb{R}}^d$ and $1\le i\le M$, where $e_i$ is the $M\times 1$ column vector whose ith entry is one and zero elsewhere. Let $R:{\mathbb{R}}^d\times {\mathbb{R}}^d\to \mathbb{R}$ be defined byfor $\mathbf{s},\mathbf{t}\in {\mathbb{R}}^d$, where $K_m:{\mathbb{R}}^d\to \mathbb{R}$ is given byfor $\mathbf{z}\in {\mathbb{R}}^d$,and $\mathrm{\Gamma }(·)$ denotes the gamma function.

Proposition 6.

Let $U_n\ge 0$ be given. Under Assumption 1, a solution to Problem (A) exists by Proposition 1. Furthermore, for any solution ${\tilde{f}}_n$ to Problem (A), there exists a solution ${\widehat{f}}_n$ to Problem (A), satisfying ${\tilde{f}}_n(X_i)={\widehat{f}}_n(X_i)$ for $1\le i\le n$, which can be represented by

$${\widehat{f}}_n(\mathbf{x})={\displaystyle \sum _{i=1}^M{\widehat{c}}_i}{p}_i(\mathbf{x})+{\displaystyle \sum _{i=1}^n{\widehat{d}}_i}R(\mathbf{x},X_i)$$(4)

for $\mathbf{x}\in {\mathbb{R}}^d$, where $({\widehat{c}}_1,\cdots ,{\widehat{c}}_M,$ ${\widehat{d}}_1,\cdots ,{\widehat{d}}_n,$ ${\widehat{f}}_n(X_1),\cdots ,{\widehat{f}}_n(X_n))$ is a solution to the following convex program in the decision variables $c_1,\cdots ,c_M,$ $d_1,\cdots ,d_n$, $y_1,$ $\cdots ,y_n\in \mathbb{R}$.

$$\begin{array}{l}\mathit{\text{Minimize}}\hspace{0.17em}\frac{1}{n}{\displaystyle \sum _{i=1}^n{(Y_i-y_i)}^2}\hfill \\ \mathit{\text{subject}}\hspace{0.17em}to\hspace{0.17em}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{n}R}}(X_i,X_j)d_i{d}_j\le U_n\hfill \\ \hspace{0.17em}{\displaystyle \sum _{i=1}^M{c}_i}{p}_i(X_j)+{\displaystyle \sum _{i=1}^n{d}_i}R(X_j,X_i)=y_j\hspace{1em}1\le j\le n.\hfill \end{array}$$(5)

Conversely, any function of form (4), where $({\widehat{c}}_1,\cdots ,{\widehat{c}}_M,$ ${\widehat{d}}_1,\cdots ,{\widehat{d}}_n,$ ${\widehat{f}}_n(X_1),\cdots ,{\widehat{f}}_n(X_n))$ is any solution to (6), is a solution to Problem (A).

Proposition 7.

Let $S_n\ge 0$ be given. Under Assumption 1, a solution to Problem (B) exists by Proposition 2. Furthermore, ${\widehat{g}}_n$ is a solution to Problem (B) if and only if ${\widehat{g}}_n$ is represented by

$${\widehat{g}}_n(\mathbf{x})={\displaystyle \sum _{i=1}^M{\widehat{c}}_i}{p}_i(\mathbf{x})+{\displaystyle \sum _{i=1}^n{\widehat{d}}_i}R(\mathbf{x},X_i)$$(6)

for $\mathbf{x}\in {\mathbb{R}}^d$, where $({\widehat{c}}_1,\cdots ,{\widehat{c}}_M,$ ${\widehat{d}}_1,\cdots ,{\widehat{d}}_n$, ${\widehat{g}}_n(X_1),\cdots ,{\widehat{g}}_n(X_n))$ is a solution to the following convex program in the decision variables $c_1,\cdots ,c_M,$ $d_1,\cdots ,d_n$, $y_1,$ $\cdots ,y_n\in \mathbb{R}$.

$$\begin{array}{l}\mathit{\text{Minimize}}\hspace{0.17em}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{n}R}}(X_i,X_j)d_i{d}_j\hfill \\ \mathit{\text{subject}}\hspace{0.17em}to\hspace{0.17em}\frac{1}{n}{\displaystyle \sum _{i=1}^n{(Y_i-y_i)}^2}\le S_n\hfill \\ {\displaystyle \sum _{i=1}^M{c}_i}{p}_i(X_j)+{\displaystyle \sum _{i=1}^n{d}_i}R(X_j,X_i)=y_j\hspace{1em}1\le j\le n.\hfill \end{array}$$

Remark 1.

It should be noted that when $d=1$, the representation in (4) coincides with the form of a natural spline; see Greville [25, equation 4.1 on p. 3] for the details.

Remark 2.

Proposition 6 reveals that for any solution ${\tilde{f}}_n$ to Problem (A), there exists a solution ${\widehat{f}}_n$ to Problem (A) satisfying ${\tilde{f}}_n(X_i)={\widehat{f}}_n(X_i)$ for $1\le i\le n$, which is a linear combination of the $p_i$’s and $R(·,X_i)$’s. In addition, Proposition 7 implies that a solution to Problem (B) is always a linear combination of the $p_i$’s and $R(·,X_i)$’s.

5. Consistency and Rates of Convergence

In this section, we focus on the consistency and convergence rates of the solutions to Problems (A) and (B). Toward this goal, we need the following assumptions.

Assumption 2.

1. $2m>d$.

2. $X,X_1,X_2,\cdots$ are iid ${[a,b]}^d$-valued random vectors having a common positive continuous density function $\tau :{[a,b]}^d\to \mathbb{R}$. This implies that there exist ${\tau }_1,{\tau }_2>0$ such that ${\tau }_1\le \tau (x)\le {\tau }_2$ for all $x\in {[a,b]}^d$.

Assumption 3.

1. $f_{*}\in {\mathcal{F}}_m$.

2. $(X,Y),(X_1,Y_1),(X_2,Y_2),\cdots$ is a sequence of iid ${[a,b]}^d\times \mathbb{R}$-valued random vectors satisfying $Y=f_{*}(X)+\epsilon$ and $Y_i=f_{*}(X_i)+{\epsilon }_i$ for $i=1,2,\cdots$ .

3. $\mathbb{E}(\epsilon |X)=\mathbb{E}({\epsilon }_i|X_i)=0$ and $\mathbb{E}({\epsilon }^2|X)=\mathbb{E}({\epsilon }_i^2|X_i)={\sigma }^2<\infty$ for $i=1,2,\cdots$ .

Assumption 4.

$\epsilon ,{\epsilon }_1,{\epsilon }_2,\cdots$ are uniformly sub-Gaussian random variables (i.e., there exist positive real numbers A and B such that $\mathbb{E}(\exp (t\epsilon ))\le A \exp (B{t}^2)$ for any $t\in \mathbb{R}$).

Assumption 5.

1. $\mathbb{E}(f_{*}{(X)}^2)<\infty$.

2. $f_{*}$, restricted to ${[a,b]}^d$, is not a polynomial of degree less than or equal to $m-1$.

3. $f_{*}$ is continuous on ${[a,b]}^d$.

We need Assumption 2(i) for the following reason. When $2m>d$ is not assumed, a function in ${\mathcal{F}}_m$ is not well defined on a set of measure 0 because the functions in ${\mathcal{F}}_m$ are determined only outside a set of measure 0. Under the assumption $2m>d$, it is known that the functions in ${\mathcal{F}}_m$ are continuous, so the point evaluation is well defined; see Meinguet [51, theorem 1 on p. 130].

The next lemma shows that Assumption 2 guarantees a very important property of the $X_i$’s, the so-called ${\mathcal{P}}_{m-1}$-unisolvency. The proof of Lemma 2 is provided in Appendix A.

Lemma 2.

Under Assumption 2, $\{X_1,\cdots ,X_n\}$ is a set of mutually distinct points containing a ${\mathcal{P}}_{m-1}$-unisolvent set for n sufficiently large a.s. In other words, Assumption 2 implies Assumption 1 for n sufficiently large a.s.

Lemma 2 ensures that the solutions to Problems (A) and (B) exist for n sufficiently large a.s. under Assumption 2.

We are now ready to present our main results. Theorems 1 and 2 reveal that ${\widehat{f}}_n$ and ${\widehat{g}}_n$ are consistent estimators of $f_{*}$ and that $D^{\alpha }\widehat{f}$ and $D^{\alpha }{\widehat{g}}_n$ are consistent estimators of $D^{\alpha }{f}_{*}$ for $|\alpha |=1,\cdots ,m-1$. Theorem 3 computes a lower bound on the rate of convergence of ${\widehat{f}}_n$. We discuss the rate of convergence of ${\widehat{g}}_n$ in Remark 4.

Theorems 1 and 2 use the following assumptions on ${\widehat{f}}_n$ and ${\widehat{g}}_n$.

Assumption 6.

There exists a constant $c_A$ such that ${\text{lim sup}}_{n\to \infty }J({\widehat{f}}_n)\le c_A$ for all n sufficiently large a.s.

Assumption 7.

There exists a constant $c_B$ such that ${\text{lim sup}}_{n\to \infty }J({\widehat{g}}_n)\le c_B$ for all n sufficiently large a.s.

Assumption 8.

$$\frac{1}{n}{\displaystyle \sum _{i=1}^n({\widehat{f}}_n(}{X}_i)-f_{*}(X_i){)}^2\le \frac{2}{n}{\displaystyle \sum _{i=1}^n{\epsilon }_i}({\widehat{f}}_n(X_i)-f_{*}(X_i))+A_n$$

for a sequence of random variables $(A_n:n\ge 1)$ satisfying ${\text{lim sup}}_{n\to \infty }{A}_n\le 0$ a.s.

Assumption 9.

$$\frac{1}{n}{\displaystyle \sum _{i=1}^n{({\widehat{g}}_n(X_i)-f_{*}(X_i))}^2}\le \frac{2}{n}{\displaystyle \sum _{i=1}^n{\epsilon }_i}({\widehat{g}}_n(X_i)-f_{*}(X_i))+B_n$$

for a sequence of random variables $(B_n:n\ge 1)$ satisfying ${\text{lim sup}}_{n\to \infty }{B}_n\le 0$ a.s.

Although Assumptions 6–9 are somewhat hard to verify, we can establish them by imposing some conditions on $U_n$ and $S_n$. For example, Assumptions 6, 8, and 9 are implied by Assumptions 10, 11, and 12, respectively.

Assumption 10.

${\text{lim sup}}_{n\to \infty }{U}_n<\infty$.

Assumption 11.

$J(f_{*})\le U_n$ for all n sufficiently large a.s.

Assumption 12.

${\text{lim sup}}_{n\to \infty }{S}_n\le {\sigma }^2$.

We are now ready to present Theorem 1, Theorem 2, Corollary 1, and Corollary 2.

Theorem 1.

Under Assumptions 2, 3, 6, and 8,

$$\underset{\mathbf{x}\in {[a,b]}^d}{\sup }|{\widehat{f}}_n(\mathbf{x})-f_{*}(\mathbf{x})|\to 0\hspace{1em}a.s.,$$(7)

$$\mathbb{E}{({\widehat{f}}_n(X)-f_{*}(X))}^2\to 0\hspace{1em} \mathit{\text{and}} \hspace{1em}\mathbb{E}{(D^{\alpha }{\widehat{f}}_n(X)-D^{\alpha }{f}_{*}(X))}^2\to 0$$(8)

for $|\alpha |=1,\cdots ,m-1$ as $n\to \infty$.

Theorem 2.

Under Assumptions 2, 3, 7, and 9,

$$\underset{\mathbf{x}\in {[a,b]}^d}{\sup }|{\widehat{g}}_n(\mathbf{x})-f_{*}(\mathbf{x})|\to 0\hspace{1em}a.s.,$$(9)

$$\mathbb{E}{({\widehat{g}}_n(X)-f_{*}(X))}^2\to 0\hspace{1em} \mathit{\text{and}} \hspace{1em}\mathbb{E}{(D^{\alpha }{\widehat{g}}_n(X)-D^{\alpha }{f}_{*}(X))}^2\to 0$$(10)

for $|\alpha |=1,\cdots ,m-1$ as $n\to \infty$.

Corollary 1.

Under Assumptions 2, 3, 10, and 11, (7) and (8) hold.

Corollary 2.

Under Assumptions 2, 3, 7, and 12, (9) and (10) hold.

We now turn to the convergence rate of ${\widehat{f}}_n$. We need Assumption 13, which is implied by Assumption 11.

Assumption 13.

$$\frac{1}{n}{\displaystyle \sum _{i=1}^n({\widehat{f}}_n(}{X}_i)-f_{*}(X_i){)}^2\le \frac{2}{n}{\displaystyle \sum _{i=1}^n{\epsilon }_i}({\widehat{f}}_n(X_i)-f_{*}(X_i))$$

for all n sufficiently large a.s.

Under Assumption 13, the convergence rate of ${\widehat{f}}_n$ is established in Theorem 3.

Theorem 3.

Under Assumptions 2, 3, 4, 6, and 13,

$${\left\{\frac{1}{n}{\displaystyle \sum _{i=1}^n{({\widehat{f}}_n(X_i)-f_{*}(X_i))}^2}\right\}}^{1/2}={\mathcal{O}}_p(n^{-m/(2m+d)})$$(11)

as $n\to \infty$.

Corollary 3.

Under Assumptions 2, 3, 4, 10, and 11, (11) holds.

Remark 3.

It should be noted that the rate $-m/(2m+d)$ established in Theorem 3 is identical to the optimal rate of convergence for the estimation of a function of smoothness of order m, which was established by Stone [65]. However, as Corollary 3 states, this rate is achieved under Assumptions 10 and 11, which may be hard to verify if sufficient information on $J(f_{*})$ is not available. One should also note that there are many estimators that achieve this optimal rate. See Kohler et al. [38] for a kernel estimator; Cox [14] for a smoothing spline estimator; Györfi et al. [27, theorems 4.3, 5.2., and 6.2] for partitioning, kernel, and nearest neighbor estimators, respectively; and Kohler et al. [37] for partitioning and nearest neighbor estimators. Furthermore, Stone [65] established that the optimal rate of convergence for the estimation of the rth derivative of a function of smoothness of order m is $-(m-r)/(2m+d)$ for $r=0,1,\cdots ,m-1$, and this rate is achieved by the smoothing spline estimator proposed by Cox [14].

We next discuss a situation where the solution to Problem (B) exists uniquely.

Theorem 4.

Consider the problem:

$$\mathit{\text{Minimize}}\hspace{0.17em}\mathbb{E}{(c_1{p}_1(X)+\cdots +c_M{p}_M(X)-f_{*}(X))}^2$$(12)

over $c_1,\cdots ,c_M\in \mathbb{R}$. Under Assumption 5, (12) has a solution, and the minimizing value ${\nu }^{*}$ of (12) is positive. Furthermore, if ${\text{lim sup}}_{n\to \infty }{S}_n<{\sigma }^2+{\nu }^{*}$ a.s., under Assumption 2, Assumption 3(ii), Assumption 3(iii), and Assumption 5, the solution to Problem (B) exists uniquely for n sufficiently large a.s.

The proofs of Theorems 1–4 and Corollaries 1–3 are provided in Appendix B.

Remark 4.

The same convergence rate as in Theorem 3 can be obtained for ${\widehat{g}}_n$ under the assumption

$$\frac{1}{n}{\displaystyle \sum _{i=1}^n{({\widehat{g}}_n(X_i)-f_{*}(X_i))}^2}\le \frac{2}{n}{\displaystyle \sum _{i=1}^n{\epsilon }_i}({\widehat{g}}_n(X_i)-f_{*}(X_i))$$(13)

for all n sufficiently large a.s. (Arguments similar to those in the proof of Theorem 3 can be applied.) However, we were not able to find a suitable condition on $S_n$ under which (13) holds.

Remark 5.

In the case where the error ${\epsilon }_1$ is dependent on $X_1$, our results in Theorems 1, 2, 3, and 4 hold with light modifications in the proofs.

Remark 6.

In many practical situations, the variance of the $Y_i$’s is heterogeneous (i.e., $\text{Var}({\epsilon }_1)$ is dependent on $X_1$). A natural candidate for $S_n$ in such a case is $S_n^{\prime }\triangleq \mathbb{E}(\text{Var}({\epsilon }_1|X_1))=\mathbb{E}(\text{Var}(Y_1|X_1))={\int }_{{[a,b]}^d}\text{Var}(Y_1|\mathbf{x})\tau (\mathbf{x})d\mathbf{x}$, where $\tau (·)$ is the density function of $X_1$. When the $X_i$’s are uniformly distributed, $S_n^{\prime }$ can be estimated by computing the average of $\text{Var}(Y_i|X_i)$ across $i\in \{1,\cdots ,n\}$.

6. Numerical Results

In this section, we observe the numerical behavior of the solutions to Problems (A), (B), and (C). To describe the detailed procedure, we start by noticing that, in practice, we can often observe $Y=f_{*}(X)+\epsilon$ multiple times at a fixed design point $X\in {[a,b]}^d$. Thus, we assume that we can collect iid observations $((X_i,Y_{ij}):1\le i\le n,1\le j\le r)$, where $Y_{ij}=f_{*}(X_i)+{\epsilon }_{ij}$ for $1\le i\le n$ and $1\le j\le r$ and the ${\epsilon }_{ij}$’s, conditional on $X_i$, are iid copies of $\epsilon$ conditional on X. We then compute the average of the $Y_{ij}$’s across multiple replications (i.e., ${\overline{Y}}_i={\sum }_{j=1}^r{Y}_{ij}/r$ for each $i\in \{1,\cdots ,n\}$). Next, the ${\overline{Y}}_i$’s will be used in place of the $Y_i$’s in Problems (A), (B), and (C).

Next, we provide more detailed descriptions on how the solutions to Problems (A)–(C) are obtained.

6.1. Problem (B)

In order to estimate $S_n$, we observe that Remark 6 suggests using the average of the $\text{Var}(Y_i|X_i)$’s as an estimate of $S_n$ when the $X_i$’s are uniformly distributed. Because we use the ${\overline{Y}}_i$’s in place of the $Y_i$’s, we use the average of $\text{Var}({\overline{Y}}_i|X_i)=(1/r)\text{Var}(Y_i|X_i)$ across $i\in \{1,\cdots ,n\}$ as an estimate of $S_n$. Therefore, our estimate of $S_n$ is

where ${\widehat{S}}_i={\sum }_{j=1}^r{(Y_{ij}-{\overline{Y}}_i)}^2/(r-1)$.

We then solve the convex programming problem in Proposition 7 with the $Y_i$’s replaced by the ${\overline{Y}}_i$’s and $S_n$ replaced by (14) using CVX, a package for solving convex programs developed by Grant and Boyd [22]. The solution obtained this way is denoted by ${\widehat{g}}_n$, and its objective value $J({\widehat{g}}_n)$ is denoted by ${\widehat{U}}_n$. ${\widehat{U}}_n$ is used as an estimate of $U_n$ when solving Problem (A).

6.2. Problem (A)

We next solve the convex programming problem in Proposition 6 by using ${\widehat{U}}_n=J({\widehat{g}}_n)$ as an estimate of $U_n$ and the ${\overline{Y}}_i$’s in place of the $Y_i$’s. CVX is used to solve the convex programming formulation.

6.3. Problem (C)

The key issue here is how to determine the smoothing parameter $\lambda$. We use cross validation to determine the smoothing parameter. We also use a fixed sequence of real numbers $({\lambda }_n:n\ge 1)$ decreasing to zero as $n\to \infty$. The smoothing parameter $\lambda$ chosen by cross validation depends on the observed data set $((X_i,Y_{ij}):1\le i\le n,1\le j\le r)$, whereas ${\lambda }_n$ does not depend on the observed data set and depends on n only.

To define the cross validation function, let $f_{\lambda }^{[k]}$ be the minimizer of

over $f\in {\mathcal{F}}_m$. Then, the cross validation function $V(\lambda )$ is defined byfor $\lambda \ge 0$, and the estimator of the smoothing parameter based on cross validation is the minimizer of $V(\lambda )$ over $\lambda \ge 0$. We call this estimator ${\widehat{\lambda }}_{CV}$. The solution to Problem (C) based on cross validation is computed by solving Problem (C) with $\lambda$ replaced by ${\widehat{\lambda }}_{CV}$ and the $Y_i$’s replaced by the ${\overline{Y}}_i$’s; see Wahba [71, equations (2.4.23) and (2.4.24) on p. 33] for the details on how Problem (C) is solved using systems of linear equations.

We also solve Problem (C) with the $Y_i$’s replaced by the ${\overline{Y}}_i$’s and $\lambda$ replaced by ${\lambda }_n$ for $n\ge 1$. The specific choice of $({\lambda }_n:n\ge 1)$ is given in Sections 6.5 and 6.6.

6.4. How to Select m

One needs a priori knowledge on m to determine the right value of m. It should be noted that the smoothness of a function $g:\mathbb{R}\to \mathbb{R}$ is measured by ${\int }_{-\infty }^{\infty }{\{g^{(2)}(x)\}}^{2}dx$. Thus, if $f_{*}$ itself is believed to be smooth, then we set $m=2$. If $f_{*}$ is believed to have smooth partial derivatives of order 2, then we set $m=4$.

6.5. M/M/1 Queue

We consider the case where $f_{*}(x)$ is the steady-state mean waiting time of a customer in a single-server queue under the first come, first served discipline with infinite capacity buffer, exponential service times, exponential interarrival times, unit arrival rate, and a service rate of $x\in [1.5,2.0]$. $f_{*}(x)$ is known to be $1/x(x-1)$; see Hillier and Lieberman [31, p. 781]. We chose an M/M/1 queue because there exists a closed-form formula for the steady-state waiting time of a customer, so we can compare our estimators with the true values. We set $X_i=1.5+i/(2n)-1/(4n)$ for $1\le i\le n$. For each $i\in \{1,\cdots ,n\}$ and $j\in \{1,\cdots ,r\}$ with $r=100$, we generated $Y_{ij}$ by averaging the waiting times of the first 1,000 customers arriving at the queue when the queue is initialized empty and idle. We next computed the solution ${\widehat{g}}_n$ to Problem (B) and the solution ${\widehat{f}}_n$ to Problem (A) with the $Y_i$’s replaced by the ${\overline{Y}}_i$’s and $m=4$. $S_n$ is estimated from (14), and $U_n$ is estimated from $J({\widehat{g}}_n)$. Both Problems (A) and (B) are solved using CVX.

For Problem (C), we computed ${\widehat{\lambda }}_{CV}$ by minimizing $V(\lambda )$ in (15) over $\lambda \in \{{10}^{-11},5·{10}^{-11},{10}^{-10},5·{10}^{-10},\cdots ,1\}$. In addition to ${\widehat{\lambda }}_{CV}$, we used three sequences of real numbers, $\lambda (1)=({10}^{-6}/n:n\ge 1)$, $\lambda (2)=({10}^{-7}/n:n\ge 1)$, and $\lambda (3)=({10}^{-8}/n:n\ge 1)$, for $\lambda$. To select $\lambda (2)$, we tried several sequences of real numbers and selected the one generating the lowest empirical integrated mean square error (EIMSE) values for $f_{*}$ and $f_{*}^{(2)}$. The rate of decay for $({\lambda }_n:n\ge 1)$ was chosen to be of order $1/n$ based on the suggestion made by Cox [13].

To measure the accuracy of ${\widehat{g}}_n$, the solution to Problem (B), we computed the EIMSE of ${\widehat{g}}_n$ as follows:

We also computed the EIMSE of ${\widehat{g}}_n^{(2)}$ as follows:

where ${\widehat{g}}_n^{(2)}$ and $f_{*}^{(2)}$ denote the second derivatives of ${\widehat{g}}_n$ and $f_{*}$, respectively. The EIMSE values of the solutions to Problems (A) and (C) and those of their second derivatives are computed similarly.

Table 1 reports the 95% confidence intervals of the EIMSE values of the estimators of $f_{*}$ computed from Problems (A)–(C) based on 400 iid replications for each n value. The columns labeled as (A) and (B) report the results obtained from Problems (A) and (B), respectively. The columns labeled as CV, $\lambda (1)$, $\lambda (2)$, and $\lambda (3)$ report the results obtained from Problem (C) when $\lambda$ is set as ${\widehat{\lambda }}_{CV}$, $\lambda (1)$, $\lambda (2)$, and $\lambda (3)$, respectively. The values in Table 1 are measured in ${10}^{-4}$. Table 2 reports the 95% confidence intervals of the EIMSE values of the estimators of $f_{*}^{(2)}$ computed from Problems (A)–(C) based on 400 iid replications for each n value.

Table 1. The 95% confidence intervals of the EIMSE values for the estimators of $f_{*}$ computed from Problems (A)–(C) when $f_{*}$ is the steady-state mean waiting time of a customer in an M/M/1 queue. All values are measured in ${10}^{-4}$.

Table 2. The 95% confidence intervals of the EIMSE values for the estimators of $f_{*}^{(2)}$ computed from Problems (A)–(C) when $f_{*}$ is the steady-state mean waiting time of a customer in an M/M/1 queue.

6.6. Stock Trader’s Problem

We consider the stock trader’s problem that motivated this paper. Consider a trader who trades stock options on a daily basis. Each day, his decision is either buy or sell a stock option, and he bases his decision on gamma, which is the second derivative of the value of the stock option as a function of the underlying stock price. The value of the stock option, $f_{*}$, typically has no closed-form formula, so simulation is required to estimate the value of the stock option and its second derivative with respect to the underlying stock price. Because simulation takes a long time to conduct, the estimation of $f_{*}$ and $f_{*}^{(2)}$ is done a day before any decision is made. Because $f_{*}$ and $f_{*}^{(2)}$ depend on the underlying stock price and because we cannot precisely predict the underlying stock price for the next day, $f_{*}$ and $f_{*}^{(2)}$ are typically estimated over a range of possible values of the underlying stock price. Hence, the trader’s goal is to estimate gamma, $f_{*}^{(2)}$, over a range of possible values of the underlying stock price. Furthermore, in his estimation, he does not want to use any parameter that is not estimated from the data set.

To place this problem in a more specific setting, we consider the case where $f_{*}(x)$ is the price of a European call option on a nondividend-paying stock when the underlying stock price is $x\in [0,2]$. We chose a European call option because there exists a closed-form formula for its price, so we can compare our estimates with the true values. Specifically, we consider the case where the strike price is $K=1.3$, the risk-free annual interest rate is $s=0.03$, the stock price volatility is 0.3, the underlying stock has a drift of 0.03, and the time to maturity of the option is $T=1$ year. It is known that the price of this option is given by

for $x\in [0,2]$, where $N(·)$ is the cumulative probability distribution function of a standard normal random variable; see, for example, Hull [32, p. 295].

To compute our estimators, we set $X_i=2i/n-1/n$ for $1\le i\le n$. For each $X_i$, we assumed that the current stock price is $X_i$ and simulated a sample path $(S_t:0\le t\le T)$ of a geometric Brownian motion up to time T with a drift of 0.03 and a volatility of 0.3. We then computed $Y_{ij}$, the price of the option by computing $\exp (-sT)\max (0,S_T-K)$, where $S_T$ is the value of the geometric Brownian motion we generated at time T. This procedure was repeated $r=5,000$ times, generating $((X_i,Y_{ij}):1\le i\le n,1\le j\le r)$. We next computed the solutions to Problems (A) and (B) with the $Y_i$’s replaced by the ${\overline{Y}}_i$’s and $m=4$. $S_n$ is estimated from (14), and $U_n$ is estimated from $J({\widehat{g}}_n)$. Both Problems (A) and (B) are solved with CVX.

For Problem (C), we computed ${\widehat{\lambda }}_{CV}$ by minimizing $V(\lambda )$ in (15) over $\lambda \in \{{10}^{-11},5·{10}^{-11},{10}^{-10},5·{10}^{-10},\cdots ,1\}$. In addition to ${\widehat{\lambda }}_{CV}$, we used three sequences of real numbers, $\lambda (4)=({10}^{-5}/n:n\ge 1)$, $\lambda (5)=({10}^{-6}/n:n\ge 1)$, and $\lambda (6)=({10}^{-7}/n:n\ge 1)$, for $\lambda$. To select $\lambda (5)$, we tried several sequences of real numbers and selected the one generating the lowest EIMSE values for $f_{*}$ and $f_{*}^{(2)}$. The rate of decay for $({\lambda }_n:n\ge 1)$ was chosen to be of order $1/n$ based on the suggestion made by Cox [13].

Table 3 reports the 95% confidence intervals of the EIMSE values of the estimators of $f_{*}$ computed from Problems (A)–(C) based on 400 iid replications for each n value. The columns labeled as (A) and (B) report the results obtained from Problems (A) and (B), respectively. The columns labeled as CV, $\lambda (4)$, $\lambda (5)$, and $\lambda (6)$ report the results obtained from Problem (C) when $\lambda$ is set as ${\widehat{\lambda }}_{CV}$, $\lambda (4)$, $\lambda (5)$, and $\lambda (6)$, respectively. The values in Table 3 are measured in ${10}^{-5}$.

Table 3. The 95% confidence intervals of the EIMSE values for the estimators of $f_{*}$ computed from Problems (A)–(C) when $f_{*}$ is the price of a European call option. All values are measured in ${10}^{-5}$.

Table 4 reports the 95% confidence intervals of the EIMSE values of the estimators of $f_{*}^{(2)}$ computed from Problems (A)–(C) based on 400 iid replications for each n value. The values in Table 3 are measured in ${10}^{-1}$.

Table 4. The 95% confidence intervals of the EIMSE values for the estimators of $f_{*}^{(2)}$ computed from Problems (A)–(C) when $f_{*}$ is the price of a European call option. All values are measured in ${10}^{-1}$.

6.7. Observations from Numerical Experiments

When we compare Problems (A) and (B) with Problem (C) with the smoothing parameter estimated from cross validation, Problem (B) shows superior performance when estimating the second derivative of $f_{*}$. Problem (C) generates extreme estimates in a few percents of time, thereby degrading the overall performance. This phenomenon has been observed by other researchers; see the discussion in Wahba [71, p. 65].

When we compare Problems (A) and (B) with Problem (C) with a fixed smoothing parameter ${\lambda }_n$ for $n\ge 1$, Problem (C) shows superior performance. However, it should be noted that to select ${\lambda }_n$, we tried several sequences of real numbers and selected the one generating the lowest EIMSE values for $f_{*}$ and $f_{*}^{(2)}$. This was possible because we already knew what $f_{*}$ and $f_{*}^{(2)}$ were in all of our numerical examples. In reality, one does not have a priori knowledge of the exact values of $f_{*}(X_i)$ and $f_{*}^{(2)}(X_i)$ for $1\le i\le n$, so finding $({\lambda }_n:n\ge 1)$ this way is not feasible. However, in our numerical experiments, the results from Problem (C) with a fixed smoothing parameter ${\lambda }_n$ give us a sense of how well Problems (A) and (B) perform compared with the best possible performance of Problem (C). Tables 1–4 report that the performance of Problems (A) and (B) is as good as that of Problem (C) with the best possible choice of $({\lambda }_n:n\ge 1)$.

6.8. When Multiple Observations Cannot Be Made at a Fixed Design Point in ${[a,b]}^d$

Throughout Section 6, we assumed that multiple observations can be made at any fixed design point $\mathbf{x}$ in ${[a,b]}^d$. These multiple observations were used to empirically estimate $S_n$. In many applications, multiple observations at a point $\mathbf{x}\in {[a,b]}^d$ may not be available. In such cases, one needs to use other ways to empirically estimate $S_n$. One possible approach is partitioning ${[a,b]}^d$ into a finite number of hyperrectangles, computing the sample variance of the $Y_i$’s whose $X_i$’s fall within each hyperrectangle, and computing the average of these sample variances (weighted by the density function of X evaluated at each hyperrectangle) as an estimate of $S_n$. This approach is based on the observation that $S_n$ should ideally be $\mathbb{E}(\text{Var}(Y|X))$.

To illustrate this approach numerically, we consider the following example: $f_{*}:[0,1]\to \mathbb{R}$ is given by $f_{*}(x)={(x-1/4)}^2$ for $x\in [0,1]$, $m=4$, $X_i=i/n-1/(2n)$ for $1\le i\le n$, and $Y_i=f_{*}(X_i)+{\epsilon }_i$ for $1\le i\le n$, where the ${\epsilon }_i$’s follow a uniform distribution over $[-0.25,0.25]$. We assumed that only one observation $Y_i=f_{*}(X_i)+{\epsilon }_i$ is made for $1\le i\le n$. To estimate $S_n$, we partitioned [0, 1] into the five intervals $I^1=[0,0.2],I^2=[0.2,0.4],I^3=[0.4,0.6],I^4=[0.6,0.8]$, and $I^5=[0.8,1.0]$ and computed the sample variances, say $S^1,S^2,S^3,S^4,$ and $S^5$, in which $S^k$ is the sample variance of the $Y_i$’s whose $X_i$’s fall in $I^k$ for $1\le k\le 5$. We set $S_n$ equal to the average of $S^1,S^2,S^3,S^4$, and $S^5$. We then solved Problem (B) to obtain ${\widehat{g}}_n$ and ${\widehat{g}}_n^{(2)}$. This procedure was repeated 1,600 times independently. Using these 1,600 replications, we computed the 95% confidence intervals of the EIMSE values of ${\widehat{g}}_n$ and ${\widehat{g}}_n^{(2)}$, respectively. Table 5 reports the 95% confidence intervals for a variety of n values. The EIMSE values in Table 5 decrease as n increases, indicating that the choice of $S_n$ is appropriate.

Table 5. The 95% confidence intervals of the EIMSE values for ${\widehat{g}}_n$ and ${\widehat{g}}_n^{(2)}$ when $S_n$ was set as the average of sample variances across five intervals.

7. Conclusions

We conclude this paper with some discussions on future research topics.

7.1. How to Choose $S_n$ When Multiple Observations Are Not Available

In Section 6.8, we discussed a way to empirically estimate $S_n$ when multiple observations are not available at a design point in ${[a,b]}^d$. In this approach, we divide ${[a,b]}^d$ into smaller hyperrectangles, compute the sample variance of the $Y_i$’s whose $X_i$’s fall within each hyperrectangle, and compute the average of these sample variances (weighted by the density function of X evaluated at each hyperrectangle) as an estimate of $S_n$. We conjecture that this approach will generate an asymptotically consistent estimator if the volume of the hyperrectangles decreases to zero and the number of $X_i$’s within each hyperrectangle increases to infinity as $n\to \infty$. A more rigorous study on this conjecture is a good future research topic.

7.2. Combining the Smoothness Condition with Shape Constraints

Because our proposed formulations, Problems (A) and (B), impose only the smoothness condition on the underlying function $f_{*}$, the proposed estimators do not necessarily produce a function that satisfies additional shape conditions on $f_{*}$, such as convexity or monotonicity. For example, even if $f_{*}$ is convex, our estimators are not necessarily convex. Figure 1 shows the solution to Problem (B) when $f_{*}:[0,1]\to \mathbb{R}$ is given by $f_{*}(x)={(x-1/4)}^2$ for $x\in [0,1]$, $m=4$, $n=10$, $X_i=i/n-1/(2n)$ for $1\le i\le n$, and $Y_i=f_{*}(X_i)+{\epsilon }_i$ for $1\le i\le n$, where the ${\epsilon }_i$’s follow a uniform distribution over $[-0.5,0.5]$ and $S_n=1/12$. Even if $f_{*}$ is convex, ${\widehat{g}}_n$ is not necessarily convex.

One possible fix to this phenomenon is incorporating additional conditions to our formulation. A good topic for future research is combining shape constraints, such as convexity and monotonicity, with the smoothness condition in our proposed formulations.

7.3. Establishing Consistency Under a Condition That Is Easier to Verify

Corollary 2 is established under Assumption 7, which is hard to verify. We conjecture that the conclusions of Corollary 2 hold under a more natural condition that $S_n\to {\sigma }^2$ as $n\to \infty$. This remains as an open question. Another open question is computing the convergence rate of ${\widehat{g}}_n$ under suitable conditions on $S_n$.

7.4. Extending Our Formulations to Shape-Constrained Estimation Problems

Problem (A) can be seen as an estimation problem of a smooth function by minimizing the sum of squared errors under the smoothness condition $J(f)\le U_n$. This type of formulation has been used for shape-constrained estimation problems, such as convex regression. For example, Mazumder et al. [49] proposed estimating an unknown convex function by minimizing the sum of squared errors over convex functions satisfying a Lipschitz condition. An interesting question for future research is how to extend Problem (B) to shape-constrained problems, such as convex and isotonic regression problems.