Estimating Derivatives with Local Polynomials

Why Derivatives Require Rescaling

A local polynomial fit approximates the regression function near a point $x$ using powers of the rescaled distance $u=(X_i-x)/h$. Its fitted coefficients therefore describe powers of $u$, not directly powers of $X_i-x$. Recovering a derivative requires undoing that bandwidth scaling and accounting for the factorial in a Taylor expansion. These factors matter: omitting either changes the estimate, and omitting one from a bias constant gives the wrong leading term.

We consider random-design regression with response $Y_i$, covariate $X_i$, and regression function $m(t)=E[Y_i\mid X_i=t]$. The observations $(X_i,Y_i)$, $i=1,\ldots,n$, are independent and identically distributed. The target is the $\nu$-th derivative $m^{(\nu)}(x)$, where $\nu$ is a nonnegative integer and the fitted polynomial degree $p$ satisfies $0\leq \nu\leq p$.

Definition: Let $r_p(u)=(1,u,\ldots,u^p)^\mathsf{T}$, and let $K$ be a nonnegative kernel supported on a bounded interval. The degree-$p$ local polynomial coefficient vector $\hat\beta(x;h)$ minimizes $$ \sum_{i=1}^n K\!\left(\frac{X_i-x}{h}\right) \left\{Y_i-r_p\!\left(\frac{X_i-x}{h}\right)^\mathsf{T}b\right\}^2 $$ over $b\in\mathbb R^{p+1}$. Write $\hat\beta_\nu(x;h)$ for its coordinate indexed by $\nu$, with coordinates indexed from $0$ through $p$. The corresponding estimate of the $\nu$-th derivative is $$ \hat m^{(\nu)}(x)=\nu!h^{-\nu}\hat\beta_\nu(x;h). $$

The scaling follows directly from Taylor's formula: $$ m(x+hu)=\sum_{j=0}^p \frac{m^{(j)}(x)}{j!}h^j u^j+\text{remainder}. $$ Thus the coefficient of $u^j$ in the local expansion is $h^j m^{(j)}(x)/j!$. Solving for the derivative gives $m^{(j)}(x)=j!h^{-j}$ times that coefficient. In particular, estimating a second derivative requires multiplication by $2h^{-2}$, not merely $h^{-2}$.

Worked Example: Recovering Derivatives From a Quadratic Fit

Suppose the regression function is $m(t)=1+2t+3t^2$. Fit a degree-two local polynomial at $x=1$ using bandwidth $h=0.2$. Because this regression function is exactly quadratic, its local Taylor polynomial has no remainder.

First calculate the Taylor coefficients in the rescaled coordinate $u=(t-1)/0.2$. Since $m(1)=6$, $m'(1)=8$, and $m''(1)=6$, the polynomial is $$ m(1+0.2u)=6+(0.2)(8)u+\frac{(0.2)^2(6)}{2}u^2 =6+1.6u+0.12u^2. $$

The fitted coefficients are therefore $\beta_0=6$, $\beta_1=1.6$, and $\beta_2=0.12$. For the first derivative, rescale by $h^{-1}$: $1.6/0.2=8$. For the second derivative, rescale by $2!h^{-2}$: $2(0.12)/(0.2)^2=0.24/0.04=6$. Both estimates recover the true derivatives exactly.

The factor $2!$ is essential. Using only $h^{-2}\beta_2$ would give $3$, half of the true second derivative.

Bias of a Local Polynomial Derivative

The bias order depends both on the polynomial degree and on which derivative is being estimated. The statement below concerns an interior point $x$, so the kernel window lies inside the covariate support for sufficiently small $h$. It also makes explicit the random-design density factor; the leading factor $f_X(x)$ cancels from the bias expansion.

Theorem: Let $p$ be fixed, let $m$ have a continuous $(p+1)$-st derivative near $x$, and suppose the density $f_X$ is continuous at $x$ with $f_X(x)>0$. Let $K$ be bounded, nonnegative, compactly supported, and such that $$ M=\int K(u)r_p(u)r_p(u)^\mathsf{T}\,du $$ is nonsingular. Define $e_\nu$ as the coordinate vector with a one at index $\nu$ and zeros elsewhere, and set $v_{p+1}=\int K(u)r_p(u)u^{p+1}\,du$. For the population local polynomial fit, the derivative estimate has bias $$ \frac{\nu!m^{(p+1)}(x)}{(p+1)!}h^{p+1-\nu} e_\nu^\mathsf{T}M^{-1}v_{p+1} +o\!\left(h^{p+1-\nu}\right). $$ The displayed term is the leading term when its coefficient is nonzero. If it vanishes, the bias is of smaller order and requires examining further terms.

Here the population fit means the minimizer of the expected weighted squared-error criterion. To derive the expansion, write $A_h=\int K(u)r_p(u)r_p(u)^\mathsf{T}f_X(x+hu)\,du$. Its population coefficient vector $\beta_h$ solves $$ A_h\beta_h=\int K(u)r_p(u)m(x+hu)f_X(x+hu)\,du. $$ Let $a_j=h^jm^{(j)}(x)/j!$, and let $a=(a_0,\ldots,a_p)^\mathsf{T}$. Taylor expansion gives $$ m(x+hu)=r_p(u)^\mathsf{T}a+ \frac{m^{(p+1)}(x)}{(p+1)!}h^{p+1}u^{p+1} +o(h^{p+1}) $$ uniformly over the bounded support of $K$. Subtract $A_ha$ from the normal equation. Continuity of $f_X$ gives $A_h\to f_X(x)M$, while the residual on the right is $$ \frac{m^{(p+1)}(x)}{(p+1)!}h^{p+1}f_X(x)v_{p+1} +o(h^{p+1}). $$ Since $M$ is nonsingular and $f_X(x)>0$, it follows that $$ \beta_h-a= \frac{m^{(p+1)}(x)}{(p+1)!}h^{p+1}M^{-1}v_{p+1} +o(h^{p+1}). $$ The coordinate $a_\nu$ equals $h^\nu m^{(\nu)}(x)/\nu!$. Multiplying the coefficient error in coordinate $\nu$ by $\nu!h^{-\nu}$ therefore gives the theorem’s bias expansion, including its factor $\nu!$. $\square$

When the leading constant is nonzero, the generic bias order is $h^{p+1-\nu}$. Differentiation makes the estimate more sensitive to local variation: each additional derivative reduces the power of $h$ in this bias order by one. The constant can nevertheless be zero because of kernel moments or symmetry. In that case the displayed expansion remains valid, but it does not identify the first nonzero bias term.

Worked Example: A Second-Derivative Bias With a Symmetric Kernel

Take $p=2$, estimate the second derivative at $x=0$, and use the uniform kernel $K(u)=1/2$ for $-1\leq u\leq 1$, zero otherwise. Let $m(t)=t^4$. This example shows why the nominal $h^{p+1-\nu}$ term can vanish and why a later term may determine the actual bias.

In the rescaled coordinate, the response function is $m(hu)=h^4u^4$. By symmetry, the coefficient of $u$ in the population quadratic fit is zero. For the intercept and quadratic coefficients, the normal equations are $$ \beta_0+\frac{1}{3}\beta_2=\frac{h^4}{5}, \qquad \frac{1}{3}\beta_0+\frac{1}{5}\beta_2=\frac{h^4}{7}. $$ Here the moments used are $\int K(u)u^2du=1/3$, $\int K(u)u^4du=1/5$, and $\int K(u)u^6du=1/7$. Subtracting one-third of the first equation from the second gives $$ \left(\frac{1}{5}-\frac{1}{9}\right)\beta_2 =h^4\left(\frac{1}{7}-\frac{1}{15}\right), \qquad \frac{4}{45}\beta_2=\frac{8}{105}h^4. $$ Thus $\beta_2=(6/7)h^4$. The estimated second derivative is $2h^{-2}\beta_2=(12/7)h^2$, while the true second derivative $m''(0)$ is zero. The cubic Taylor term is zero for $t^4$ at zero, and symmetry also cancels its relevant moment; the first nonzero bias here is of order $h^2$.

Variance and Bandwidth Scaling

For derivative estimation, the factor $h^{-\nu}$ also amplifies random variation. Suppose additionally that $\operatorname{Var}(Y_i\mid X_i=t)=\sigma^2(t)$, with $\sigma^2$ continuous at $x$, and that the conditional errors have a bounded $2+\delta$ moment near $x$ for some $\delta>0$. Define $$ R=\int K(u)^2r_p(u)r_p(u)^\mathsf{T}\,du. $$ Under the preceding design and kernel conditions, with $h\to0$ and $nh\to\infty$, the leading variance is $$ \operatorname{Var}\{\hat m^{(\nu)}(x)\} \sim \frac{\nu!^2\sigma^2(x)}{nh^{2\nu+1}f_X(x)} e_\nu^\mathsf{T}M^{-1}RM^{-1}e_\nu. $$ This is a leading asymptotic variance: it describes the first-order contribution from the conditional response noise. Local smoothness and the stated moment conditions make the variation due to estimating the local design matrices lower order.

To see the scaling and constant, the empirical design matrix, divided by $nh$, converges to $f_X(x)M$. The covariance of the noise part of the corresponding coefficient score, divided by $nh$, converges to $\sigma^2(x)f_X(x)R$. Inverting the design matrix on both sides yields coefficient covariance asymptotic to $\sigma^2(x)(nhf_X(x))^{-1}M^{-1}RM^{-1}$. Converting coefficient $\nu$ to a derivative multiplies this covariance by $(\nu!h^{-\nu})^2$, giving the stated expression. In particular, the variance order is $1/(nh^{2\nu+1})$.

Worked Example: Computing a First-Derivative Variance

Use a local linear fit, so $p=1$, with the uniform kernel $K(u)=1/2$ on $[-1,1]$. Suppose $f_X(x)=1$ and $\sigma^2(x)=4$, and take $n=1000$, $h=0.2$. The moments give $$ M= \begin{pmatrix}1&0\\0&1/3\end{pmatrix}, \qquad R= \begin{pmatrix}1/2&0\\0&1/6\end{pmatrix}. $$ For the first derivative, $e_1=(0,1)^\mathsf{T}$ and $\nu!=1$. Therefore $$ e_1^\mathsf{T}M^{-1}RM^{-1}e_1 =3\left(\frac{1}{6}\right)3=\frac{3}{2}. $$ The leading variance is consequently $$ \frac{4}{1000(0.2)^3}\cdot\frac{3}{2} =\frac{6}{0.8}=7.5. $$ This is an asymptotic variance approximation, not an exact finite-sample variance. Its $h^{-3}$ factor illustrates the increased variability of a first-derivative estimate as the bandwidth shrinks.

If the leading bias has nonzero order $h^q$, the leading squared bias is of order $h^{2q}$, while the variance is of order $1/(nh^{2\nu+1})$. Balancing these orders gives $h$ of order $n^{-1/(2q+2\nu+1)}$. For the generic nonzero bias term above, $q=p+1-\nu$, yielding $n^{-1/(2p+3)}$. If a moment cancellation changes the bias order, the balance changes as well; one should not use the generic exponent without checking the applicable bias term.

Common Mistakes

Forgetting the factorial. The coefficient of $u^\nu$ represents $h^\nu m^{(\nu)}(x)/\nu!$. The derivative estimate must therefore use $\nu!h^{-\nu}\hat\beta_\nu$.
Scaling the bias only by $h^{p+1}$. That is the coefficient-error scale. Converting the coefficient to a derivative adds $h^{-\nu}$, so the derivative bias scale is $h^{p+1-\nu}$, with the additional factor $\nu!$ in its constant.
Assuming the generic bias constant is always nonzero. Kernel moments can cancel it. The formula gives a leading term only when its displayed constant is nonzero; a zero constant calls for a higher-order expansion.
Ignoring derivative amplification of variance. Multiplication by $h^{-\nu}$ multiplies variance by $h^{-2\nu}$, in addition to the usual local sample-size factor $1/(nh)$.

Key takeaway: A degree-$p$ local polynomial estimates derivative $\nu$ by $\nu!h^{-\nu}\hat\beta_\nu$. Its generic derivative bias is of order $h^{p+1-\nu}$, and its variance is of order $1/(nh^{2\nu+1})$; check the bias constant for moment cancellations before using that generic order as the actual leading order.

Check Your Understanding

Use the coefficient scaling and the bias and variance expansions to answer these questions.

A local cubic fit uses $h=0.1$ and has fitted coefficient $\hat\beta_2=0.015$. What second-derivative estimate does this coefficient imply?
For degree $p=3$, what is the generic bias order for estimating the first derivative, provided the leading constant is nonzero?
Why does a second-derivative estimate have a variance factor $h^{-5}$ rather than $h^{-3}$?
In the quartic example, which symmetry property makes the cubic-order bias term vanish?
If the coefficient of the nominal leading bias term is zero, what additional step is needed to identify the actual bias order?

Estimating Derivatives with Local Polynomials

What You'll Learn