Nonparametric estimation of coefficient function derivatives in varying coefficient models

1.   Introduction

As an important extension of nonparametric models, varying-coefficient models have evolved through three distinct methodological epochs since their conceptual inception by Shumway ^[1]. This framework fundamentally reimagines traditional parametric specifications by enabling covariate effects to fit the effects of predictors that change across levels of other variables, making them particularly useful in scenarios with heterogeneous data or when relationships between variables are not constant. Many statisticians proposed some effective methods to estimate it; see, e.g., Hastie and Tibshirani ^[2] for the smooth spline method, Fan and Zhang ^[3] for the two-step method, Huang and Shen ^[4] for the polynomial spline method, Cai et al. ^[5] for the local kernel smoothing method.

The derivative estimation of the coefficient functions is equally crucial. The first derivative of coefficient functions precisely characterizes the instantaneous rate of change at any given point and facilitates the detection of inflection points of the coefficient function, while enabling the analysis of marginal effect in economic ^[6]. Crucially, the second derivative constitutes the foundation for statistical inference on coefficient functions ^[5], and its proper estimation is essential for constructing valid confidence intervals. However, a large number of researchers have concentrated predominantly on the estimation of the coefficient functions itself, inadvertently causing the derivative functions of the coefficient functions to languish in relative obscurity. Within the framework of the varying-coefficient models, the estimation of derivative functions has been restricted to the local polynomial method ^[7] and the smoothing spline method ^[8]. These methods require the selection of smoothing parameters. However, existing parameter selection approaches yield only satisfactory estimation for the coefficient functions, while producing suboptimal results for derivative function estimation under the same parameters. An intuitive explanation for the failure of smoothing parameters is that coefficient functions and their derivatives often have different smoothness levels. The selected smoothing parameters meet only requirements of the coefficient functions, resulting in poor estimation performance for the derivative functions. Together, these considerations require the development of novel estimation methods specifically designed for derivative functions of coefficient functions.

We now shift our focus to a special case of varying-coefficient models: nonparametric models. Within the framework of nonparametric models, researchers have developed a diverse array of derivative estimation methods that have significantly advanced the field. Consequently, we can draw upon these estimation techniques to devise a novel approach to estimate derivative functions of coefficient functions. In the realm of nonparametric models, local polynomial regression ^[9,10] and smoothing spline ^[11] remain the predominant methods for derivative estimation. Although both methodologies necessitate the selection of a smoothing parameter, the challenge of optimally choosing this parameter remains inadequately resolved. Furthermore, Wang and Zhou ^[12] employed the Bayesian inference framework to estimate derivative functions under the assumption that the unknown function follows a Gaussian process. In addition, Liu and Li ^[13] proposed a plug-in kernel ridge regression estimator to estimate the derivative function.

Unlike these approaches, the difference-quotient method offers an alternative avenue for derivative estimation in nonparametric models. Müller et al. ^[14] introduced derivative estimators based on difference quotients. However, these estimators are characterized by high variance, which complicates the smoothing process. To mitigate this issue, Charnigo et al. ^[15] proposed a linear combination to reduce the variance of symmetric quotients, based on the assumption of an equispaced design between independent variables. However, the random variation inherent in independent variables within random design settings often invalidates theoretical conclusions derived from equispaced design assumptions, limiting their applicability. In response to this, several researchers have developed novel procedures in the equispaced design framework, which have produced promising results ^[16,17,18]. In pursuit of a symmetric difference quotient applicable under random design conditions, Liu and Brabanter ^[19] elucidated the theoretical properties of derivative estimation utilizing weighted symmetric quotients, assuming that independent variables follow a standard uniform distribution. Empirical results demonstrate the superiority of their method over conventional local polynomial and spline approaches. Furthermore, Liu and Kong ^[20] constructed empirical derivatives using weighted symmetric difference quotients to estimate derivatives under the condition of error correlation. In view of the great estimation performance of the weighted symmetric quotients procedure in nonparametric models, we can draw on the weighted symmetric quotient method to estimate the derivatives of coefficient functions.

In this paper, we consider the estimation of the derivatives in a varying-coefficient model where there are differences in the smoothness of the coefficient functions. First, based on the weighted symmetric difference quotients idea, an effective method for estimating the first and second derivatives of the coefficient functions is proposed under the case where the index variables are a uniform random design. Compared to traditional local kernel estimation and spline methods that cannot select appropriate smoothing parameters, we propose an effective approach to choosing the tuning parameters in our estimator, which leads to superior estimation performance over conventional methods. Then, we derive the corresponding theoretical properties and generalize the estimators to the case that index variables are arbitrary distributions.

The paper is organized as follows. In Section 2, we show the first-order derivative estimation of coefficient functions based on variance reducing weighted difference quotients and establish the theoretical properties of conditional bias, variance. In addition, a method for selecting the optimal tuning parameter is given. In Section 3, we extend the procedure to second-order derivatives of the coefficient. In Sections 4 and 5, we perform Monte Carlo experiments and real data analysis to compare the proposed methodology with the local polynomial regression and the two-step method. Some conclusions are presented in Section 6, and the proof details of the theorems are provided in Appendix.

2.   First order derivative estimation

2.1.   Estimation approach of first order derivative

From now on, let $\{(Y_{i}; U_{i}, X_{i1}, \cdots, X_{ip}), i = 1, \cdots, n\}$ be a set of independent random samples from a population $(\textbf{X}, U, Y)$; then, we consider the varying coefficient model

where $Y$ is response variable, $X_{1}, \cdots, X_{p}$ are covariate variables, $U$ is index variable, $a_{j}(\cdot)(j = 1, \cdots, p)$ are unknown coefficient functions, and $\varepsilon$ is the regression error which satisfies $E(\varepsilon|U, X_{1}, \cdots, X_{p}) = 0$ and $Var(\varepsilon|U, X_{1}, \cdots, X_{p}) = \sigma^{2}$.

Now, we consider the special case where the index variable $U$ follows a uniform distribution. Let $U\sim\mathcal{U}(0, 1)$, where $\mathcal{U}(0, 1)$ is the standard uniform distribution between 0 and 1. In addition, let $U_{(i)}, i = 1, 2, \cdots, n$ be the order statistics of $U$ and $\hat{a}_q(\cdot), q = 1, 2, \cdots, p$ be a two-step estimator ^[3] of the coefficient functions which are obtained by using local linear fitting in the second step. Inspired by the success of symmetric difference quotient methods ^[18,19,20] in nonparametric regression for derivative estimation, we can construct weighted symmetric difference quotient estimators to estimate the derivative of $a_j(\cdot)$. Without loss of generality, we consider the derivative estimate of $a_p(\cdot)$. If we want to estimate the derivative of $a_q(\cdot)$, we can swap the position of $a_q(\cdot)$ with $a_p(\cdot)$. Therefore, the proposed weighted symmetric quotient estimators of first order derivative estimator for uniform order statistics is

where $\hat{a}^{(1)}_{p}(U_{(i)})$ denotes the derivative estimator in $U_{(i)}$, and $\omega_{i, j}$ denote weights and satisfy $\omega_{i, j}\geq0$, $\sum_{j = 1}^k\omega_{i, j} = 1$, then $k$ is a given turning parameter.

The idea behind the symmetric quotient estimator in Equation (2.2) mainly originates from the definition of a derivative. Given a point $u_0$, the first-order derivative $a^{(1)}_1(u_0) = \lim_{\delta \to 0}\frac{a_1(u_0+\delta)-a_1(u_0-\delta)}{2\delta}$. Based on the definition, given a point $U_{(i)}$, we can naturally derive the following estimator for the first order derivative

(2.3) represents the special case of (2.2) when $k = 1$. However, in this scenario, the estimator utilizes only information from points $U_{(i-1)}$ and $U_{(i+1)}$, making it highly susceptible to noise interference. To mitigate this, we require additional estimators analogous to the derivative definition, denoted $\frac{\hat{a}_p(U_{(i+j)})-\hat{a}_p(U_{(i-j)})}{U_{(i+j)}-U_{(i-j)}}$, and obtain a robust derivative estimator through weighted averaging, resulting in (2.2). In particular, $k$ in this context plays a role similar to the smoothing parameter in the estimation of the local polynomial: when $k$ is too small, excessive noise is introduced; When $k$ is too large, the estimator becomes overly smoothed. We conducted an experiment to illustrate the impact of varying $k$ on the estimation of the first-order derivative, the results of which are visualized in Figure 1. As shown in Figure 1, when $k = 1$, the estimator shows significant fluctuations due to noise interference. In contrast, when $k = 120$, excessive smoothing leads to poor estimation near the function's peak. Only when $k = 50$ does the estimator achieve satisfactory performance. Therefore, the choice of $k$ is critical in determining the quality of the estimator.

Furthermore, we employ the two-step estimator instead of the local polynomial estimator, motivated by potential differences in the smoothness of the coefficient functions in (2.2). When there is a difference in the smoothness of the coefficient functions, local polynomial estimation may struggle to provide accurate estimates, resulting in a substantial difference between the estimated and true functions. This considerable estimation error can propagate to the estimation of (2.2), thereby hindering the accuracy of the derivative function estimates. However, according to Fan and Zhang ^[3], the two-step method effectively enhances the accuracy of the coefficient function estimation when their smoothness varies. Consequently, we adopt the two-step method to obtain the estimator $\hat{a}_p(\cdot)$.

In (2.2), the weights $\omega_{i, j}$ are unknown. Before delving into the analysis of these weights $\omega_{i, j}$, we first consider the conditional bias and conditional variance order of $\hat{a}^{(1)}_{p}(U_{(i)})$ when $j = 1$. Now, we present the following assumptions.

(A1) $E(X^{2s}_{j}) < \infty$ holds for $s > 2$ and $j = 1, \cdots, p$.

(A2) $nh^4_2 > n^{\rho_1}$, for some $\rho_1 > 0$. $n^{-\rho_2}h^{-1}_2 = o_p(h^2_2)$, for some $\rho_2 > 0$.

(A3) $nh^{\gamma}_1/log(h_1)\rightarrow \infty$, for any $\gamma > s/(s-2)$.

(A4) The function $K(t)$ is a finite symmetric density function with compact support.

(A5) The function $r_{kj}(u) = E(X_{k}X_{j}|U = u)$ has a bounded second derivative in a neighborhood of $u_{0}$ and $r^{(2)}_{kj}(u_{0})\ne0$.

(A6) The coefficient $a_{j}(\cdot)$ has a continuous second derivative.

Assumptions (A1) and (A3)-(A6) are often used in the varying-coefficient model literature ^[3,5,21]. Assumption (A2) imposes specific requirements on the convergence rate of $h_2$, and generally the condition holds as long as $h_2 = o(n^{-0.1})$ is satisfied for some $\rho_1 > 0$ and $\rho_2 > 0$.

Proposition 1. Under the assumptions (A1)-(A6) and $h_1\rightarrow0, h_2\rightarrow0, k\rightarrow \infty$ as $n\rightarrow \infty$, such that $nh_1\rightarrow \infty$, $h_1/h_2\rightarrow0$, second derivative $a^{(2)}_{j}(\cdot)$ is Lipschitz continuous, then let

thus the asymptotic order of condition bias and condition variance of $\hat{a}^{(1)}_{p}(U_{(i)})$ for $j = 1$ are given by

According to Proposition 1, under $j = 1$, as $n$ grows large, the bias term asymptotically converges to zero, while the variance of the estimator grows without bound, resulting in an unreliable estimate. To mitigate this risk of uncontrolled variance, the proposed weighting scheme effectively controls the variance of the parameter $\hat{a}^{(1)}_{p}(U_{(i)})$. The following proposition provides a method to determine the optimal weights $\omega_{i, j}$.

Proposition 2. For $k+1\leq i \leq n-k$, the weights $\omega_{i, j}$ that minimize an upper bound of conditional variance of (2.2), satisfying $\sum^k_{j = 1}\omega_{i, j} = 1$, are given by

Note that optimal weights are obtained by minimizing an upper bound of the conditional variance of (2.2) rather than minimizing the conditional variance. However, in the next section, it is shown that under the optimal weights given by Proposition 2, the conditional variance can be controlled with appropriate parameter settings, and the performance of our proposed estimator is better than the method in simulation.

Remark 1. Our estimator is formulated under the assumption that the index variable follows a standard uniform distribution. In particular, this simplifying assumption aligns with established practices in the literature on varying-coefficient modeling ^[21,22,23], where researchers often impose uniformity on the index variable for empirical implementation. Although this convention facilitates tractable analysis, we explicitly address its limitations in the following sections by deriving derivative estimates when the variable follows nonuniform distributions.

2.2.   Asymptotic properties of the first-order derivative estimator

We establish the upper bounds of asymptotic conditional bias and variance of our proposed estimator (2.2) for the interior points $k+1\leq i \leq n-k$. Now, we will use the following notation. Define

and

for $d = 1, \cdots, p-1$.

Theorem 1. Under assumptions (A1)-(A6) and $h_1\rightarrow0, h_2\rightarrow0, k\rightarrow \infty$ as $n\rightarrow \infty$, such that $nh_1\rightarrow \infty$, $h_1/h_2\rightarrow0$, $k^2/n^2h_2\rightarrow0$, the upper bounds of conditional bias and the conditional variance of $\hat{a}_p(U_{(i)})$ are given by

and

As demonstrated in Theorem 1, we establish asymptotic upper bounds for conditional bias and conditional variance. This corollary provides a precise asymptotic order for conditional bias and conditional variance, thereby characterizing their exact rates of convergence.

Corollary 1. Under assumptions (A1)-(A6) and $h_1\rightarrow0, h_2\rightarrow0, k\rightarrow \infty$ as $n\rightarrow \infty$, such that $nh_1\rightarrow \infty$, $h_1/h_2\rightarrow0$, $k^2/n^2h_2\rightarrow0$, then the exact asymptotic order for conditional bias and conditional variance of $\hat{a}^{(1)}_p(U_{(i)})$ are given by

According to Corollary 1, we can obtain

Now, we compare the MSE of this method with those of the two-step method and the local polynomial method. According to the proof of Fan and Zhang's theorem ^[3], when $a_{j}(\cdot)$ has a continuous second order derivative, we can obtain the convergence rates of the bias and variance of the estimated first-order derivative, where the local linear fitting is employed in the second step, thus $Bias\left[\hat{a}^{(1)}_{tse, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(\text{max}\left\{h^2_{tse, 2}, \, \frac{h^2_{tse, 1}}{h_{tse, 2}}\right\}\right)$ and $Var\left[\hat{a}^{(1)}_{tse, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(\frac{1}{nh^3_{tse, 2}}\right)$, then $MSE\left[\hat{a}^{(1)}_{tse, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(\text{max}\left\{h^4_{tse, 2}, \, \frac{h^4_{tse, 1}}{h^2_{tse, 2}}\, , \frac{1}{nh^3_{tse, 2}}\right\}\right)$. If we select $h_2 = O(n^{-0.13})$, $h_1 = O(n^{-0.14})$, $h_{tse, 2} = O(n^{-0.2})$ and $h_{tse, 1} = O(n^{0.21})$, we have $MSE\left[\hat{a}^{(1)}_{tse, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(n^{-0.4}\right)$. Thus, we need only select $k = O(n^{-l})$ ($0.77\leq l\leq0.8$), then $MSE\left[\hat{a}^{(1)}_p(U_{(i)})|\mathcal{D}\right] = o_p(n^{-0.4})$. Notice that $h_{tse, 2} = O(n^{-0.2})$ is optimal to minimize the mean squared error (MSE) of the coefficient function in the two-step estimation method. Generally, existing bandwidth selection methods aim to maximize the accuracy of the coefficient function estimation. Therefore, our choice of $h_{tse, 2} = O(n^{-0.2})$ is theoretically justified. Therefore, when the parameters $k$, $h_1$, and $h_2$ are appropriately chosen, the convergence rate of $MSE\left[\hat{a}^{(1)}_p(U_{(i)})|\mathcal{D}\right]$ is superior to that of $MSE\left[\hat{a}^{(1)}_{tse, p}(U_{(i)})|\mathcal{D}\right]$.

For local polynomial estimation based on linear expansion, $Bias\left[\hat{a}^{(1)}_{lpe, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(h^2\right)$ and $Var\left[\hat{a}^{(1)}_{lpe, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(\frac{1}{nh^3}\right)$, $MSE\left[\hat{a}^{(1)}_{lpe, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(\text{max}\left\{h^4, \, \frac{1}{nh^3}\right\}\right)$. If we select $h = O(n^{-0.2})$, $h_2 = O(n^{-0.13})$ and $h_1 = O(n^{-0.14})$, we have $MSE\left[\hat{a}^{(1)}_{lpe, p}(U_{(i)})|\mathcal{D}\right] = O_p\left(n^{-0.4}\right)$. Then, by appropriately choosing $k = O(n^{-l})$ ($0.77\leq l\leq0.8$), we can ensure $MSE\left[\hat{a}^{(1)}_p(U_{(i)})|\mathcal{D}\right] = o_p(n^{-0.4})$. The bandwidth $h = O(n^{-0.2})$ also follows the optimal bandwidth that minimizes the mean squared error (MSE) of the coefficient function in the local polynomial framework.

We complete the theoretical analysis of $\hat{a}^{(1)}_p(U_{(i)}$, but the method for obtaining the parameter $k$ remains unknown in (2.2). According to the analysis in Section 2.1, parameter $k$ controls the estimation performance. Notice that the MSE measures the estimation performance of an estimator and controls the bias-variance trade-off of the estimator. Therefore, we can minimize the upper bound of the conditional MSE to choose the optimal turning parameter $k$.

Corollary 2. Under the assumptions of Theorem 1, the optimal $k$ is obtained by minimizing the following formula

Although a fast and easy parameter selection method is provided in Corollary 2, some unknown quantities need to be estimated. The error variance can be estimated by

The $\underset{u\in[0, 1]}{sup}\left|a^{(2)}_p(u)\right|$ can be estimated using the two-step method with a local polynomial expansion of 3 order in the second step. The $E(X_dX_p|U)$ can be obtained based on observation $(U_i, X_{id}X_{ip})$ using the local linear fitting with bandwidth that is selected by cross-validation. Then, $r^-_{pp}, r^+_{dp}, r^*_{dp}$ can be calculated by taking the maximum or minimum of the fitting values of the corresponding $E(X_dX_q|U)$. By adding the estimators above, the optimal value $k_{opt}$ can be obtained by searching the grid over the set of integers $[1, \lfloor\frac{n-1}{2}\rfloor]$, where $\lfloor\cdot \rfloor$ denotes round down.

Due to limitations in these smoothing parameter selection schemes, the first order derivatives of the coefficient functions cannot be reliably estimated for local polynomial estimation, two-step methods, and spline estimation approaches. In Corollary 2, we propose an effective approach for selecting the parameter $k$. Subsequent simulation studies demonstrate the validity of our parameter selection method, showing that it yields derivative estimates with superior accuracy compared to alternative approaches. Furthermore, it should be noted that the estimation of the derivatives also requires the estimation of $\hat{a}_p(U_{i})$, which involves the selection of smoothing parameters $h_1$ and $h_2$. For this purpose, cross-validation methods are sufficient, as they adequately meet the requirement of producing a good estimate of $\hat{a}_p(U_{i})$. Our approach needs only ensure the proper estimation of $\hat{a}_p(U_{i})$, a requirement that standard cross-validation techniques can satisfy fully.

2.3.   Smoothing the estimators and generalizing results and arbitrary distributions

In the previous discussion, we were only able to obtain the first-order derivative estimate at point $U_i$, $i = k+1, \cdots, p-k$, and estimates at other points could not be obtained using the current method. Within the framework of univariate nonparametric estimation, by treating $U_i$ as the covariate variable and $\hat{a}_p(U_i)$ as the response variable, we can apply nonparametric estimation techniques to fit the first-order derivative, thus estimating the first-order derivative at other points. According to the local polynomial estimation for nonparametric regression ^[7], given an arbitrary point $u_0\in[0, 1]$ and kernel function $K(\cdot)$, we can obtain the derivative estimator through the following equation,

where $g$ represents the order of local polynomial fitting, $e_{1, g}$ represents a column vector of length $g$, and the first element is 1, and the remaining elements are 0,

Note that $u_0$ is a predetermined point. Immediately after specifying $u_0$ and the order statistics $U_{(i)}$, we can directly derive $\mathbf{U}(u_0)$ and $\mathbf{W}(u_0)$.

However, in estimating local polynomials, a smoothing parameter $h$ must be determined. Here, we use the plug-in bandwidth selector ^[10], which is a simple and effective bandwidth selection method. In subsequent simulations, we also compute the smoothed derivative estimates. The improvement in the accuracy of the estimation after smoothing further demonstrates the effectiveness of the smoothing method we propose.

In addition, the index variables $U$ do not follow a standard uniform distribution, which means that the above estimation method fails under other distributions. In the frame of nonparametric derivative estimation, to generalize covariate $Z$ with unknown distribution $F$, we can use probability integral transformation (PIT)

where $F(\cdot)$ is continuous, such that the density function $f(x) = F^{(1)}(x)$ exists and let $f$ be bounded away from zero. This means that the index variables $U$ with the distribution function $F$ satisfying the above conditions can be transformed to an index variable following a uniform distribution. Now, the index variable with unknown distribution denotes $Z$. Then, in the varying coefficient models, we use the chain rule

where $a_q(Z) = r_q(F(Z))$. In practice, since the distribution $F(\cdot)$ and density $f$ both are unknown, the kernel density estimator ^[24] can be used to estimate the distribution and density to obtain $\hat{f}(\cdot)$ and $\hat{F}(\cdot)$ based on the plug-in bandwidth $h$ ^[10], namely,

The detailed derivative estimation algorithm in Table 1 is given below.

3.   Second order derivative estimation

3.1.   Estimation approach and asymptotic results of a second derivative

In varying-coefficient models, many researchers are also interested in the second derivative; for example, in analyzing the curvature of certain coefficient functions. However, higher-order derivatives of coefficient functions become more difficult to estimate, and the performance of the traditional estimation method becomes worse and worse with increasing derivative order. In order to obtain a better estimator of the second derivative estimator of coefficient functions, we proposed second order derivative estimator for uniform order statistics, as

where $\hat{a}^{(2)}_{p}(U_{(i)})$ denotes the derivative estimator in $U_{(i)}$, and $\omega_{i, j, 2}$ denote weights and satisfy $\omega_{i, j, 2}\geq0$, $\sum_{j = 1}^k\omega_{i, j, 2} = 1$, then $k_1$, $k_2$ are given turning parameter.

The weight $\boldsymbol{\omega}_{i, 2} = \left\{\omega_{i, 1, 2}, \cdots, \omega_{i, k_2, 2}\right\}^T$ determination method is analogous to the first order derivative estimation, so the optimal weights can be derived by minimizing the conditional variance of $\hat{a}^{(2)}_p(U_{(i)})$. Thus, we obtain the optimal weights $\boldsymbol{\omega}_{i, 2}$ by solving the following optimization problem

Proposition 3. $Var\left[\hat{a}^{(2)}_p(U_{(i)})|\mathcal{D}\right]$ is a quadratic function of $\boldsymbol{\omega}_{i, 2} = \left\{\omega_{i, 1, 2}, \cdots, \omega_{i, k_2, 2}\right\}$, namely,

where $\Sigma_0$ is a $k_2\times k_2$ positive semidefinite matrix, and the specific expression for $\Sigma_0$ can be found in 6.

According to Proposition 3, we can know that the conditional variance of $\hat{a}^{(2)}_p(U_{(i)})$ can be expressed as a quadratic function of $\omega_{i, j, 2}$, which inspires us to use a quadratic programming algorithm to solve (3.2). Consequently, optimal weights $\hat{\boldsymbol{\omega}}_{i, 2}$ can be derived by solving the following quadratic programming problem with the given samples:

There are several well-established algorithms to solve quadratic programming problems; in this paper, we adopt the dual interior point method proposed by ^[25] to optimize and solve the equation (3.4). This method has been implemented in the R package 'LowRankQP' specifically for solving quadratic programming problems, thus we choose to use the R package 'LowRankQP' to solve the optimization problem. In practice, the computational speed depends on the tuning parameter $k_2$. We conduct a simulation experiment to evaluate the time required to obtain $\hat{a}^{(2)}_p(U_{(i)})$, comparing it with the two-step method involving triple expansion in its second step. The result is shown in Table 2. All computations are performed on an Intel(R) Xeon(R) CPU E5-2683 v4 processor (2.10 GHz base frequency). Table 2 reveals that the computational speed of our method is slightly slower than that of the two-step approach, indicating that its performance remains acceptable. Furthermore, employing faster optimization techniques could further reduce computation time.

However, there is a disadvantage that optimal weights are obtained using the weight of the quadratic programming algorithm, that is, we cannot get a theoretical expression for the weights $\omega_{i, j, 2}$, which means that difficulties arise for the theoretical analysis of (3.1). Thus, we assume that the weights $\omega_{i, j, 2}$ are known for the following theorem.

Theorem 2. Under assumptions (A1)–(A5), and $h_1\rightarrow0, h_2\rightarrow0, k_1\rightarrow \infty$ as $n\rightarrow \infty$, such that $nh_1\rightarrow \infty$, $h_1/h_2\rightarrow0$, $k_1^2/n^2h_2\rightarrow0$. Besides, we assume $a_p(\cdot)$ has a third-order derivative $a^{(3)}_p(\cdot)$, and the weights $\omega_{i, j, 2}$ are known. The absolute conditional bias and conditional variance of (3.1) are

and

where $\mathcal{R}_{dp} = max\{0, \underset{U\in[0, 1]}{sup}E(X_dX_p|U)\}, \mathcal{R}^{*}_{dp} = max\{0, \underset{U\in[0, 1]}{sup}e^{T}_{d, p}(E(XX^T|U))^{-1}e_{p, p}\}$, and $\mu_m$, $M$, $r^-_{pp}$, $r^+_{dp}$, $r^*_{dp}$ maintain identical definitions to those specified in Theorem 1.

Corollary 3. Under the assumptions of Theorem 2, and $k_2$ is bounded, then the exact asymptotic order for conditional bias and conditional variance of $\hat{a}^{(2)}_p(U_{(i)})$ are given by

3.2.   Optimal parameter selection for $k_1$ and $k_2$

Similarly to parameter $k$ for the first-order derivative estimation, parameter $k_1, k_2$ should also be selected, which controls the bias-variance trade-off. Based on the upper bound of conditional MSE, the following corollary provides the selection method of $k_1, k_2$.

Corollary 4. Under the assumptions of Theorem 2, we can obtain optimal $k_1, k_2$ that minimizes the upper bound of conditional MSE is

The unknown quantity $\underset{u\in[0, 1]}{sup}\left|a^{(3)}_p(u)\right|$ can be obtained by the two-step method with local polynomial expansion of order $p = 4$ in the second step. Moreover, $\mathcal{R}_{dp}, \mathcal{R}^{*}_{dp}$ are estimated in a way similar to $r^{*}_{dp}$. The estimation methods for the remaining unknown quantity $\sigma^2, r^-_{pp}, r^+_{dp}, r^*_{dp}$ are the same as those in Corollary 2. The optimal value $k_1, k_2$ can then be obtained by a grid search.

3.3.   Smoothing the second order estimators and generalizing it to arbitrary distributions

Similar to (2.13), we can also obtain the smooth estimator provided at arbitrary point $u_0\in[0, 1]$ by local polynomial smoothing, that is,

where $g$ represents the order of local polynomial fitting, $e_{1, g}$ represents a column vector of length $g$, and the first element is 1 and the remaining elements are 0,

The smoothing parameter $h$ is chosen in the same way as in Section 3.3.

To generalize to arbitrary distributions for the index variables $U$, we still use the Probability Integral Transform (PIT) as in (2.14) to transform the index variables $Z$ to $U$. Then, let $a_q(Z) = r_q(F(Z))$, the second derivative of coefficient functions with respect to $Z$ is

where the estimation of unknown $f^{(1)}(x)$ can be estimated using the kernel density derivative estimator. To take a similar procedure to Algorithm 1, we can obtain

4.   Simulation studies

4.1.   Simulation 1

In this section, we perform simulations to demonstrate the performance of the proposed methods. We consider the following varying coefficient model

For the above model, two examples are considered. The true coefficient functions in two examples are given in Table 3 based on the smoothness difference. In Example 1, the smoothness of $a_2(U)$ and $a_3(U)$ is the same, whereas in Example 2, the smoothness of all the coefficient functions differs.

In two examples, $U$ follows a uniform distribution on $\left [0, 1\right]$. For $X_{1}$ and $X_{2}$, we set $(X_{1}, X_{2}, X_{3}, X_{4})^{T}\sim N(\mu, \Sigma)$, where $\mu = (0, 0, 0, 0)^{T}$, $\Sigma = (\sigma_{ij})_{4\times4}$ and $\sigma_{ij} = \rho^{|i-j|}$, $\rho = 0.1\, \, \text{or}\, \, 0.5$.

Here, we set $(X_{1}, X_{2}, X_{3}, X_{4})$, $U$, and $\varepsilon$ to be independent of each other. Let the random error term $\varepsilon$ follow a normal distribution with mean 0 and variance $\sigma^2$, and $\sigma^{2}$ is chosen so that the signal-to-noise ratio is 5:1, namely,

For each example, we perform 200 simulations with sample sizes $n = 400,800$ in each case. To evaluate the performance of the symmetric difference quotients method (SQE), the performance of the symmetric difference quotients smoothing method (S-SQE) of order $g = 3$, the local polynomial estimator (LPE) of order $g = 3$, and the two-step estimator (TSE), where we use local polynomial fitting of order $g = 3$ in the second step, which are also computed. Additionally, we incorporate a smoothing spline estimator (SPE) based on B-spline, where the derivative estimate can be obtained directly by differentiating the B-spline basis functions. For SQE, optimal $k$ and $k_1, k_2$ are selected based on Corollary 2 and 4 in a set $\{1, 2, \ldots, \lfloor n/2\rfloor\}$ and a search space $\{1, 2, \ldots, \lfloor n/4\rfloor\}\otimes\{1, 2, \ldots, \lfloor n/4\rfloor\}$. For LPE, the optimal bandwidth $h_{LLE}$ is selected by the cross-validation method. In the two-step estimation, the initial bandwidth $h_0 = 0.5h_{LLE}$ is taken in the first step, and then we use cross-validation to select the bandwidth in the second step. Furthermore, for SPE, we employ equidistant knots, with the number of knots selected via cross-validation. For each method, the Epanechnikov kernel function is used, i.e., $K(t) = 0.75(1-t^2)_+$. All the aforementioned numerical simulations are implemented in the R software.

To systematically compare the estimation performance between methodologies, we calculate the mean absolute error (MAE) of the derivative estimates for each coefficient function. A critical challenge arises from the inherent scale incompatibility between coefficient functions. As their derivatives exhibit distinct dynamic ranges (varying upper and lower bounds), direct comparison of raw MAE values becomes dimensionally inconsistent. To resolve this, we implement a range-based standardization procedure that normalizes MAE values relative to the characteristic scale of each derivative function. Furthermore, we address the boundary effect paradox observed in LPE, TSE, and SPE, where derivative estimation accuracy deteriorates severely near domain boundaries. In addition, SQE provides reliable estimates only within the interior region of the domain. Consequently, to ensure a fair comparison, we exclude boundary points and focus on evaluating performance exclusively within the interior domain. Thus, the performance measure is defined using the adjusted mean absolute error.

and

We also calculate the average adjusted mean absolute error of the coefficient functions, denoted as mean. For each method, the mean adjusted mean absolute error based on 200 simulations is presented in Tables 4 and 5. Moreover, in Figures 2 and 3, we present the graphs of the first- and second-order derivatives of the estimated coefficient functions for Example 2, when $n = 800$ and $\rho = 0.1$. Additionally, Figures 4 and 5 present boxplots of the root mean squared error (RMSE) for the derivative function estimates under the same examples.

The results in Tables 4 and 5 show that, in terms of the mean adjusted mean absolute error, the average error of the SQE method is significantly smaller than that of the TSE, SPE and LPE methods, thus highlighting the effectiveness of the proposed approach. In the estimation of the derivative with the highest smoothness, LPE consistently performs the worst. This is because the local polynomial method can only employ a fixed bandwidth when estimating the derivative, resulting in poor estimates for the derivatives of coefficient functions with higher smoothness. The limited bandwidth leads to an inadequate approximation of the derivative, especially for functions with higher smoothness. Furthermore, SPE consistently yields poor derivative estimates for coefficient functions with the lowest smoothness, demonstrating that SPE should generally be avoided for derivative estimation. In addition, the estimation error of S-SQE relative to SQE is smaller, demonstrating the effectiveness of our smoothing method in further improving estimation accuracy and reducing the estimation error. As seen in Example 1, the estimation of the SQE and S-SQE method in $\rho = 0.1$ is better than the estimation at $\rho = 0.5$. However, the mean of the adjusted mean absolute error of SQE and S-SQE in the case $\rho = 0.5$ does not increase much compared to the case $\rho = 0.1$. However, the LPE, TSE, and LPE methods exhibit a substantial decrease in the accuracy of the estimation as the correlation increases, suggesting that our method is less sensitive to the correlation between variables. Furthermore, it can be seen that the error in the second-order derivative is larger than that in the first-order derivative, but SQE and S-SQE are able to give an estimate with lower error. In addition, as the sample size increases, the adjusted mean absolute error of the four methods becomes smaller.

From Figures 2 and 3, the estimation of the first-order derivatives of four methods is similar. However, the estimation of the fourth function in the second-order derivatives estimated by the other methods is poor, whereas S-SQE estimates the general shape of the functions. In addiction, the LPE and SPE methods show significantly poorer performance in estimating the second-order derivative compared to the TSE and S-SQE methods, which is attributed to the chosen bandwidth not being optimal for second-order derivative estimation. According to 4 and 5, the boxplot width of S-SQE is nearly the smallest among all methods, while maintaining a stable width across all derivative functions without observable outliers. This demonstrates the superior stability of our proposed estimation method.

4.2.   Simulation 2

In this section, we consider a more complicated example: A six-dimensional varying-coefficient model. We consider the following varying coefficient model

In this example, $U$ follows a uniform distribution on $\left [0, 1\right]$. For $X_{1}$ and $X_{2}$, we set $(X_{1}, X_{2}, X_{3}, X_{4}, X_{5}, X_{6})^{T}\sim N(\mu, \Sigma)$, where $\mu = (0, 0, 0, 0, 0, 0)^{T}$, $\Sigma = (\sigma_{ij})_{6\times6}$ and $\sigma_{ij} = 0.5^{|i-j|}$.

Here, we set $(X_{1}, X_{2}, X_{3}, X_{4}, X_{5}, X_{6})$, $U$, and $\varepsilon$ to be independent of each other. Let the random error term $\varepsilon$ follows a normal distribution with mean 0 and variance $\sigma^2$, and $\sigma^{2}$ is chosen so that the signal-to-noise ratio is 5:1, namely,

The mean of the adjusted mean absolute error based on 200 simulations is presented in Tables 6 and 7.

Like in Simulation 1, Tables 6 and 7 show that LPE, TSE, and SPE struggle to accurately estimate all functions among the second-order derivatives, but SQE and S-SQE maintain a relatively low estimation error. However, in this example, the TSE and SPE methods perform the worst, suggesting that as the dimensionality increases, the estimation accuracy of coefficient function derivatives using TSE and SPE declines significantly. In addition, compared to LPE, TSE, and SPE, the SQE and S-SQE methods demonstrate relatively stable estimation accuracy for the derivatives of each coefficient function. With the exception of the second-order derivative of the sixth coefficient function, the estimation errors for all other coefficient functions remain at a consistent level.

Furthermore, to demonstrate the statistically significant improvement in mean estimation accuracy of S-SQE over competing methods, we conduct paired t-tests ^[26] comparing S-SQE's mean term based on 200 simulations against those of other approaches. We obtain the corresponding p-values, all of which are less than $1\times10^{-10}$, demonstrating that S-SQE achieves significantly improved performance compared to other methods. This further validates that our method yields more stable estimates in the presence of varying smoothness across the coefficient functions.

5.   Real data analysis

To further illustrate the proposed methodologies, we choose to implement these approaches to estimate derivative functions on the Boston Housing Dataset in this section. For this dataset, calculating the first order derivatives of the coefficient functions enables the analysis of the rate of change for each coefficient function at different points, while the second order of derivatives can be used to construct confidence intervals for local polynomial estimation of the coefficient functions. For the given dataset, we construct a varying-coefficient model. Consequently, the estimation of derivative functions carries substantial significance.

The dataset encompasses 506 observations spanning 14 distinct variables and is retrievable from the R package'mlbench'. According to Fan and Huang's ^[27] modeling analysis of the data using a varying coefficient model, we designate MEDV, which denotes the median value of owner-occupied homes measured in thousands of dollars, as our response variable, and LSTAT, indicating the proportion of the population with a lower socioeconomic status, as our index variable. In terms of covariate variables, we adopt the selection criteria established by Wang and Xia through the application of the KLASSO method ^[22]. This has led to the inclusion of four predictor variables: INT (the model intercept), CRIM (the crime rate per capita within each town), RM (the average number of rooms per residential unit), and PTRATIO (the ratio of pupils to teachers in each town). Prior to the deployment of these methodologies, the data undergoes a normalization process for the covariate variables, with the exception of INT and the response variable. In addition, the index variable is subjected to a transformation to ensure that its distribution conforms to a uniform distribution between 0 and 1.

We estimate the first and second order derivatives of the coefficient functions using SQE, LPE, TSE, and SPE. The plots of the estimated first and second derivatives are presented in Figures 6 and 7. To assess the accuracy of these estimates, we identify the points corresponding to the index variables at the maxima and minima of the first-order derivative for each coefficient function at the interior points. We then calculate the mean distance from the second-order derivative estimates at these points to zero, and the results are presented in Table 8. In theory, the second-order derivative at the extrema of the first-order derivative should be zero. Therefore, a smaller distance from the second-order derivative to the zero point indicates a more accurate estimate, while a larger distance suggests a poorer estimation.

As shown in Figure 6, when LSTAT exceeds 0.5, the first-order derivatives of all coefficient functions are negative, indicating a declining trend. This aligns with empirical observations: Areas with higher proportions of lower socioeconomic status populations generally exhibit lower median home values. Notably, a local minimum occurs near LSTAT = 0.7 (corresponding to 15%), where the rate of decline peaks. Beyond this threshold, the descent rate moderates, suggesting that housing price depreciation slows when LSTAT surpasses 15%. Furthermore, the first-order derivatives of coefficient functions estimated by SPE on INT are asymptotically zero, indicating that INT's coefficient function is essentially constant. This conclusion sharply contradicts Wang and Xia's finding ^[22], which confirmed the existence of a varying-coefficient intercept term. Such inconsistency significantly undermines the credibility of SPE.

According to Table 8, we observe substantial distance among LPE, TSE, and SPE estimates, indicating their uniformly poor performance in derivative estimation. For the RM variable, all methods yield substantially large distances, likely attributable to the high amplitude of its first-order derivative function. In addition, with the exception of the RM variable, the distances estimated by SQE are very close to zero, which suggests that our method provides the most accurate estimates. In, summary, our method effectively estimates both the first and second derivatives, and accurately identifies the maxima and minima of the derivatives in this dataset.

6.   Conclusions

We consider the first and second order derivative function estimation problem of the varying coefficient model when the coefficient functions may have different degrees of smoothness, and a weighted symmetric difference quotient estimation method is proposed. Along with the proposed method for estimating derivative functions, we also present an effective approach for selecting tuning parameters. Both the simulation results and the real data analysis demonstrate the effectiveness of our method in estimating the derivative functions. Furthermore, under certain conditions, we provide the asymptotic upper bounds for the conditional bias and conditional variance of the estimator and exact order of conditional bias and variance. In addition, based on these bounds, we propose a parameter selection method derived from the asymptotic upper bound of the conditional MSE. Furthermore, relevant researchers may explore alternative methods to obtain a more precise estimate $\hat{a}_p(\cdot)$. Incorporating $\hat{a}_p(\cdot)$ higher-precision estimate of a into the weighted symmetric difference quotient estimator could potentially further enhance the accuracy of derivative estimation.

In future research, we plan to extend our method to time-varying coefficient models in economics. Additionally, we aim to investigate higher-order derivative estimation of coefficient functions, derivative estimation for multivariate covariates, and the extension of our approach to additive models for derivative function estimation. Furthermore, due to the difficulty of the proof, we do not present an exact expression of asymptotic conditional bias and conditional variance. Thus, obtaining precise asymptotic expressions for the conditional bias and variance will be essential, as this could facilitate the selection of more appropriate tuning parameters, ultimately leading to improved estimation results. Besides, when deriving the upper bounds for the asymptotic bias and variance of the second-order derivative estimator, we assume that the weights are known—a relatively stringent condition, and we aim to explore methods for obtaining the exact expression of the second-order derivative weights to circumvent this condition in the future.

Appendix: Proof of Theorem

We use the following notation.

where $I_p$ is a $p\times p$ matrix.

A1.   Proof of Proposition 2

The proof of this proposition is essentially a simplified version of Corollary 1. For detailed proof, please refer to Corollary 1.

A2.   Proof of Proposition 2

where $H_{(b)} = (1, 0)(\textbf{X}^{T}_{2, (b)}W_{2, (b)}\textbf{X}_{2, (b)})^{-1}\textbf{X}^{T}_{2, (b)}W_{2, (b)}(I_{n}-B_n)$, $H_{l, (b)}$ is the $l$-th element of $H_{(b)}$. Therefore, we obtain an upper bound of $Var(\hat{a}^{(1)}_p(U_i)|\mathcal{D})$. Then, the Lagrange multiplier method is used to solve the minimum point of $\omega_{i, j}$ for this upper bound. The optimal solution is

When the asymptotic conditional bias and variance are calculated, the following lemma on the uniform convergence is used.

Lemma 1. Let $(X_1, Y_1), \ldots, (X_n, Y_n)$ be i.i.d random vectors, where the $Y_i's$ are scalar random variables. Assume further that $E|y|^3 < \infty$ and ${sup}_x\int|y|^s f(x, y)dy < \infty$, where $f$ denotes the joint density of $(X, Y)$. Let $K$ be a bounded positive function with a bounded support, satisfying a Lipschitz condition. Then,

provided that $n^{2\epsilon-1}h\rightarrow \infty$ for some $\epsilon < 1-s^{-1}$.

Proof. This follows immediately from the result obtained by Mack and Silverman ^[28].

Lemma 2. Let $(X_1, Y_1), \ldots, (X_n, Y_n)$ be i.i.d random vectors, where the $Y_i's$ are scalar random variables, $X$ is $2$-dimention random vector with the continuous marginal density function $f(x)$, and $f(x)$ has a compact support. Besides, $nh^4 > n^{\rho_1}$, for some $\rho_1 > 0$, $n^{-\rho_2}h^{-1} = o_p(h^2)$, for some $\rho_2 > 0$. Then

completely as $n\rightarrow0$, where $t = (t_1, t_2)$.

Proof. This follows immediately from the result obtained by Cheng and Taylor ^[29].

Lemma 3. $\mathit{\text{Let}} \mathbf{U} = \{U_1, \cdots, U_n\}$ be a sample generated from the random variable $U$, where $U$ follows the standard uniform distribution. Also, let $U_{(1)}, \cdots, U_{(n)}$ be the corresponding order statistics. Then, we have:

Proof. This follows immediately from the result obtained by David and Nagaraja ^[30].

A3.   Proof of Theorem 1

For the initial estimator $\hat{a}_{q, 0}(U_{(b)})$ of the two-step estimator, we have

Note that by Taylor's expansion, we have

where $\eta_{i, (b)}$ is between $U_i$ and $U_{(b)}$ for $i = 1, \ldots n$.

Therefore,

where

According to the Lemma 1 and continuity assumption for $a^{(2)}_j(\cdot)$, we can obtain

By substituting (A.1) and (A.2) into $J_1$, we have

Similarly, we have

By combining (A.3) and (A.4), we have

According to the standard uniform distribution property in Lemma 3, we have

According to the (A.5) and (A.6), we can obtain an upper bound of the absolute conditional bias

Then, we deal with conditional variance. First, for $H^{}_{(b)}H^{T}_{(c)}$, we have

According to Lemma 2, we can obtain

and

where

Notice, $e^{T}_{2q-1, 2p}(\Omega^{-1}_p(U_l)\otimes A^{-1})C_l = (0, o_p(1))$. Therefore, we have

Let

We expand the $\tau_{rs}$, we have

where

Then, using Lemma 1, we have

Substituting (A.9) into $\tau_{rs}$, we have

where

According to Lemma 2 and the properties of the Kronecker product, we have

Substituting (A.7), (A.8), and (A.10) into the expansion of $H_{(b)} H_{(c)}^T$, we have

Now, we provide an upper bound for $Var\left[\hat{a}_p^{(1)} (U_{(i)}) | \mathcal{D} \right]$. First, the $Var \left[\hat{a}_p^{(1)} (U_{(i)}) | \mathcal{D} \right]$ is expanded by

By substituting (A.11) into the above expansion, we can obtain the upper bound for $Var \left[\hat{a}_p^{(1)} (U_{(i)}) | \mathcal{D} \right]$,

A4.   Proof of Corollary 1

According to the upper bound of the absolute conditional bias and (A.12), we have

According to the upper bound of the absolute conditional variance, we have

Thus, Corollary 1 is proven.

A4.1.   Proof of Proposition 3

Let

and $\xi = (\xi_1, \cdots, \xi_{k_2})^T$. We can obtain

Then, let $\Sigma_0 = \xi\xi^{T}$, thus Proposition 3 is proven.

A5.   Proof of Theorem 2

Similar to the proof of Theorem 1, we can obtain the upper bound of the conditional bias of $\hat{a}^{(2)}_p(U_i)$

Similar to the proof of the upper bound of the conditional variance in Theorem 1, we can substitute (A.11) into (A.13), and we can have

Thus, Theorem 2 is proven.

A6.   Proof of Corollary 3

According to the proof of Theorem 2, we can immediately obtain Corollary 3.

Author contributions

Junfeng Huo proposed the methodology, proved the theoretical properties, implemented the code, performed part of the numerical simulation, and was primarily responsible for writing and revising the manuscript. Mingquan Wang contributed to the numerical simulation and participated in the writing and revision parts of the manuscript. Xiuqing Zhou contributed to the formulation of the research topic and provided critical revisions of the manuscript.

Use of Generative-AI tools declaration

The authors declare that generative AI tools were used to refine the language of sentences in the Introduction section of this article. Specifically, ChatGPT was used to improve language clarity and readability. All other content, including research design, data analysis, and conclusions, was independently developed by the authors without AI assistance. The authors bear full responsibility for the integrity and accuracy of the manuscript.

Conflict of interest

The authors declare that they have no conflict of interest.