Higher-Order Delta Method at Stationary Functionals

Why a Different Scale Is Needed

For a smooth scalar functional $g(\theta)$, the ordinary delta method uses the linear change in $g$ near the true parameter $\theta_0$. If the gradient vanishes at $\theta_0$, that linear change is zero. A first-order approximation then does not describe the leading random variation: curvature, rather than slope, becomes decisive. The functional often fluctuates on the scale $n^{-1}$, not $n^{-1/2}$.

This situation occurs naturally when estimating a squared distance from a hypothesized value, a variance-like functional at a symmetry point, or a contrast whose first derivative vanishes at the parameter of interest. The second-order delta method identifies both the appropriate scale and the shape of the limiting distribution.

Definition: Let $g:\Theta\to\mathbb{R}$, where $\Theta$ is an open subset of $\mathbb{R}^p$. A point $\theta_0\in\Theta$ is a stationary point of $g$ if $\nabla g(\theta_0)=0$. Its Hessian is the $p\times p$ matrix $H=\nabla^2 g(\theta_0)$, with entries $H_{jk}=\partial^2g(\theta_0)/(\partial\theta_j\partial\theta_k)$. If $g$ has continuous second partial derivatives near $\theta_0$, then $H$ is symmetric.

Suppose an estimator $\hat{\theta}$ is consistent and satisfies $\sqrt n(\hat{\theta}-\theta_0)\xrightarrow{d}Z$, where $Z\sim N_p(0,V)$ and $V$ is a finite, symmetric positive semidefinite matrix. At a stationary point, the first-order term disappears. The next result makes precise what replaces it.

The Second-Order Delta Method

Theorem: Suppose $\theta_0$ is an interior point of $\Theta$, $g$ has continuous second partial derivatives in a neighborhood of $\theta_0$, $\nabla g(\theta_0)=0$, and $\sqrt n(\hat{\theta}-\theta_0)\xrightarrow{d}Z\sim N_p(0,V)$. Then $$ n\{g(\hat{\theta})-g(\theta_0)\} \xrightarrow{d} \frac{1}{2}Z^{\mathsf T}HZ, $$ where $H=\nabla^2 g(\theta_0)$.

The estimator’s root-$n$ limit is still the input to the calculation. The output scale changes because the functional’s first derivative is zero: a quadratic change in an $O_p(n^{-1/2})$ estimation error is $O_p(n^{-1})$.

Proof: Write $h=\hat{\theta}-\theta_0$. A second-order Taylor expansion with a Peano remainder gives $$ g(\theta_0+h)-g(\theta_0) = \nabla g(\theta_0)^{\mathsf T}h +\frac{1}{2}h^{\mathsf T}Hh +r(h), \qquad \frac{r(h)}{\lVert h\rVert^2}\longrightarrow 0 \quad\text{as }h\longrightarrow 0. $$ The gradient term is zero by stationarity. The assumed convergence in distribution implies $\sqrt n\,h=O_p(1)$, and hence $h=O_p(n^{-1/2})$ and $h\xrightarrow{p}0$. The remainder therefore satisfies $nr(h)=o_p(1)$: indeed, $r(h)=o_p(\lVert h\rVert^2)$ and $n\lVert h\rVert^2=O_p(1)$. Consequently, $$ n\{g(\hat{\theta})-g(\theta_0)\} = \frac{1}{2}(\sqrt n\,h)^{\mathsf T}H(\sqrt n\,h)+o_p(1). $$ The map $z\mapsto \frac12 z^{\mathsf T}Hz$ is continuous. The continuous mapping theorem and Slutsky’s theorem now give the claimed convergence. $\square$

Why the Limit Is a Weighted Chi-Square

The quadratic form in the theorem can be expressed using the eigenvalues of a matrix built from the Hessian and the asymptotic covariance. This representation identifies whether the limit is a nonnegative weighted chi-square or a signed combination of chi-squares.

Theorem: Under the second-order delta-method conditions, let $\lambda_1,\ldots,\lambda_p$ be the eigenvalues, including zeros, of the symmetric matrix $V^{1/2}HV^{1/2}$. Then $$ \frac{1}{2}Z^{\mathsf T}HZ \ \overset{d}{=}\ \frac{1}{2}\sum_{j=1}^p\lambda_j U_j, \qquad U_1,\ldots,U_p\ \text{independent with }U_j\sim\chi^2_1. $$ If every $\lambda_j\geq 0$, this is a nonnegative weighted chi-square distribution. If some eigenvalues are negative, the limit is a signed difference of weighted chi-square variables, not an ordinary chi-square distribution.

Proof: Since $V$ is symmetric positive semidefinite, it has a symmetric positive semidefinite square root $V^{1/2}$. A normal vector with covariance $V$ has the same distribution as $V^{1/2}W$, where $W\sim N_p(0,I_p)$. Thus $$ \frac{1}{2}Z^{\mathsf T}HZ \ \overset{d}{=}\ \frac{1}{2}W^{\mathsf T}(V^{1/2}HV^{1/2})W. $$ The matrix $A=V^{1/2}HV^{1/2}$ is symmetric, so the spectral theorem gives an orthogonal matrix $Q$ such that $A=Q\operatorname{diag}(\lambda_1,\ldots,\lambda_p)Q^{\mathsf T}$. Orthogonal transformations preserve the standard multivariate normal distribution, so $Q^{\mathsf T}W$ has independent standard normal coordinates. Calling these coordinates $W_1,\ldots,W_p$, the quadratic form becomes $\frac12\sum_j\lambda_jW_j^2$. Each $W_j^2$ has a $\chi^2_1$ distribution, and the squares remain independent. This proves the result. $\square$

The covariance matters as much as the curvature. Directions in which the estimator has little or no asymptotic variation contribute little or nothing to the limit, even if the functional curves in those directions. The eigenvalues combine both effects, so inspecting the Hessian alone is not enough to determine the limiting weights.

Worked Example: Squared Error at a Bernoulli Symmetry Point

Let $X_1,\ldots,X_n$ be independent Bernoulli observations with success probability $p$, and estimate $p$ by $\hat p=\bar X$. Consider $g(p)=(p-\tfrac12)^2$ when the true value is $p_0=\tfrac12$.

The derivative is $g'(p)=2(p-\tfrac12)$, so $g'(p_0)=0$; the second derivative is $g''(p)=2$. The central limit theorem gives $\sqrt n(\hat p-\tfrac12)\xrightarrow{d}Z$, where $Z\sim N(0,\tfrac14)$, because $p_0(1-p_0)=\tfrac14$. The second-order limit is therefore $$ n\{g(\hat p)-g(p_0)\} =n(\hat p-\tfrac12)^2 \xrightarrow{d} Z^2 \overset{d}{=}\frac14\chi^2_1. $$ The factor $\tfrac14$ follows by writing $Z=\tfrac12 W$ for $W\sim N(0,1)$. Here $g(p_0)=0$, so the result describes the scaled squared estimation error itself.

Worked Example: Two Curvature Directions with Different Variances

Suppose $Y_1,\ldots,Y_n$ are independent bivariate normal vectors with mean $\theta_0=(0,0)^{\mathsf T}$ and covariance $V=\operatorname{diag}(1,4)$. Use their sample mean $\hat{\theta}$ to estimate $\theta_0$, and consider $g(\theta)=\theta_1^2+\tfrac12\theta_2^2$.

The sample mean satisfies $\sqrt n\,\hat{\theta}\sim N_2(0,V)$ exactly. The gradient of $g$ vanishes at $\theta_0$, and its Hessian is $H=\operatorname{diag}(2,1)$. If $Z=(Z_1,Z_2)^{\mathsf T}\sim N_2(0,V)$, then $Z_1=W_1$ and $Z_2=2W_2$ for independent standard normal variables $W_1,W_2$. Direct substitution gives $$ n\{g(\hat{\theta})-g(\theta_0)\} \xrightarrow{d} \frac12 Z^{\mathsf T}HZ =Z_1^2+\frac12 Z_2^2 \overset{d}{=}W_1^2+2W_2^2. $$ The eigenvalues of $V^{1/2}HV^{1/2}=\operatorname{diag}(2,4)$ are $2$ and $4$, producing weights $1$ and $2$ after multiplication by $\tfrac12$. Unequal variances in the two coordinates lead to unequal limiting weights.

Worked Example: A Stationary Point with Opposite Curvatures

Let $Y_1,\ldots,Y_n$ be independent bivariate normal vectors with mean $(0,0)^{\mathsf T}$ and covariance $I_2$. Take $\hat{\theta}=\bar Y$ and $g(\theta)=\theta_1^2-\theta_2^2$. The gradient vanishes at $\theta_0=(0,0)^{\mathsf T}$, and $H=\operatorname{diag}(2,-2)$.

Here $V=I_2$, so the eigenvalues of $V^{1/2}HV^{1/2}$ are $2$ and $-2$. For independent $W_1,W_2\sim N(0,1)$, the limit is $$ n\{g(\hat{\theta})-g(\theta_0)\} \xrightarrow{d} W_1^2-W_2^2 \overset{d}{=}\chi^2_{1,1}-\chi^2_{1,2}, $$ where the two chi-square variables are independent. The limiting variable can be positive or negative because the functional curves upward in one direction and downward in the other. Calling this limit a chi-square distribution would be incorrect.

Beyond Second Order

If the Hessian contributes no variation, the second-order limit may be identically zero. For example, this happens if $V^{1/2}HV^{1/2}=0$, even when $H$ itself is not zero. In that case the second-order theorem is valid but does not identify a nondegenerate limit; a higher-order expansion may be needed.

More generally, suppose $g$ has continuous derivatives through order $r$ near $\theta_0$, its derivatives of orders $1,\ldots,r-1$ vanish there, and $\sqrt n(\hat{\theta}-\theta_0)\xrightarrow{d}Z$. Taylor expansion then yields $$ n^{r/2}\{g(\hat{\theta})-g(\theta_0)\} \xrightarrow{d} \frac{1}{r!}D^rg(\theta_0)[Z,\ldots,Z], $$ where $D^rg(\theta_0)[Z,\ldots,Z]$ is the $r$-th derivative, viewed as a multilinear form, evaluated $r$ times at $Z$. As in the second-order proof, the Taylor remainder is $o_p(\lVert\hat{\theta}-\theta_0\rVert^r)$, and root-$n$ tightness makes its product with $n^{r/2}$ vanish in probability. For $r=2$, this expression is precisely the Hessian quadratic form.

Common Mistakes and Practical Interpretation

A frequent error is to apply the ordinary delta method mechanically and conclude that a zero gradient means there is no asymptotic variation. It means only that the $n^{-1/2}$-scale linear term vanishes. The next term can still be random, and the correct scale for the functional may be $n$.

A second mistake is to assume every quadratic-form limit is chi-square. A chi-square variable is nonnegative, whereas an indefinite Hessian can produce negative eigenvalues and a signed limit. Even when all weights are nonnegative, unequal weights generally give a weighted chi-square, not a single ordinary chi-square. The eigenvalues of $V^{1/2}HV^{1/2}$, including their signs and zeros, are the relevant diagnostic.

These limits matter for inference because they describe the sampling behavior of functionals for which a first-order normal approximation fails. They also warn against applying a usual normal confidence interval on the wrong scale. Before using the approximation, check stationarity, the estimator’s joint asymptotic covariance, and the curvature in directions where the estimator varies. If the quadratic form vanishes, the result gives no nondegenerate approximation by itself.

Check Your Understanding

Use the conditions and calculations in the tutorial to answer the following.

Why does a zero gradient at $\theta_0$ change the natural scaling of $g(\hat{\theta})-g(\theta_0)$ from $n^{1/2}$ to $n$?
For the second-order limit, which matrix’s eigenvalues determine the chi-square weights, and why does the asymptotic covariance appear in that matrix?
What property of the eigenvalues ensures that the limiting quadratic form is nonnegative?
Can a nonzero Hessian still yield a degenerate second-order limit? State the relevant condition involving $V$ and $H$.
If the first nonzero derivative of $g$ at $\theta_0$ has order $r$, what is the corresponding scale for the functional difference under a root-$n$ estimator?