The Multivariate Delta Method

Why a Joint Limit Needs a Multivariate Transformation Rule

Many estimators are built from several quantities at once. A statistic may be a product of two sample means, a log ratio of two estimated probabilities, or a vector of summaries computed from a common data set. If those quantities have a joint limiting distribution, a differentiable function of them usually has a limit determined by its local linear approximation.

The multivariate delta method makes that idea precise. The derivative is now a matrix, called the Jacobian, and it transforms the entire limiting random vector. The joint limit matters: dependence among the coordinates can affect the transformed limit even when each coordinate’s marginal distribution is known.

Jacobian and First-Order Approximation

Definition: Let $U$ be an open subset of $\mathbb{R}^k$, and let $g:U\to\mathbb{R}^m$. The function $g$ is differentiable at $\theta\in U$ if there is an $m\times k$ matrix $J_g(\theta)$ such that $$ g(\theta+h)=g(\theta)+J_g(\theta)h+r(h), \qquad \frac{\|r(h)\|}{\|h\|}\longrightarrow 0 \quad \text{as } h\to 0. $$ The matrix $J_g(\theta)$ is the Jacobian of $g$ at $\theta$. Its entry in row $i$, column $j$ is the partial derivative of output coordinate $g_i$ with respect to input coordinate $j$, when those partial derivatives exist and give the derivative.

The condition says that, near $\theta$, the change in $g$ is approximately the matrix $J_g(\theta)$ applied to the change in its input. The remainder is small relative to the size of that input change. For a limit theorem, the essential issue is whether this remainder remains negligible after multiplication by the rate.

The Multivariate First-Order Delta Method

Theorem (Multivariate first-order delta method): Let $X_n$ be random vectors in $\mathbb{R}^k$, let $U\subseteq\mathbb{R}^k$ be open, and let $\theta\in U$. Suppose $X_n\in U$ almost surely for every $n$, $g:U\to\mathbb{R}^m$ is differentiable at $\theta$, $a_n>0$ is deterministic with $a_n\to\infty$, and $$ a_n(X_n-\theta)\xrightarrow{d}Z $$ for a random vector $Z$ in $\mathbb{R}^k$. Then $$ a_n\{g(X_n)-g(\theta)\}\xrightarrow{d}J_g(\theta)Z. $$

The almost-sure domain assumption is part of the theorem: it ensures $g(X_n)$ is defined for every $n$. The differentiability assumption is local, but the random inputs must approach the point where the derivative is taken. That follows from the rate condition, as the proof shows.

The scaled deviations are bounded in probability.
Convergence in distribution implies $a_n(X_n-\theta)=O_p(1)$. Since $a_n\to\infty$, this also implies $X_n\xrightarrow{p}\theta$.

Expand the transformation.
Differentiability writes the change in $g$ as its Jacobian term plus a remainder that is small relative to $\|X_n-\theta\|$.

Scale and take the limit.
The scaled remainder converges to zero in probability. The Jacobian term converges in distribution to $J_g(\theta)Z$, so Slutsky’s theorem gives the result.

Proof. Write $V_n=a_n(X_n-\theta)$. Since $V_n\xrightarrow{d}Z$, the sequence $V_n$ is tight, so $V_n=O_p(1)$. For every $\varepsilon>0$, $$ P(\|X_n-\theta\|>\varepsilon) =P(\|V_n\|>a_n\varepsilon)\longrightarrow 0, $$ because $a_n\varepsilon\to\infty$ and a tight sequence has arbitrarily small tail probabilities beyond sufficiently large bounds. Thus $X_n\xrightarrow{p}\theta$.

Differentiability gives a remainder function $r$ such that, almost surely, $$ a_n\{g(X_n)-g(\theta)\} =J_g(\theta)V_n+a_n r(X_n-\theta). $$ We verify that the second term is $o_p(1)$. For any $\eta>0$, differentiability gives a $\delta>0$ such that $0<\|h\|<\delta$ implies $\|r(h)\|\leq\eta\|h\|$; also $r(0)=0$. Given any $\epsilon>0$, choose $M$ so that $P(\|V_n\|>M)<\epsilon$ for all sufficiently large $n$, using tightness. On the event $\|V_n\|\leq M$, we have $\|X_n-\theta\|\leq M/a_n<\delta$ for sufficiently large $n$, and therefore $$ a_n\|r(X_n-\theta)\|\leq\eta a_n\|X_n-\theta\| =\eta\|V_n\|\leq\eta M. $$ For any target bound $b>0$, choose $\eta<b/M$ (with any positive $M$ satisfying the tightness bound). Then the probability that the scaled remainder exceeds $b$ is at most $\epsilon$ for all sufficiently large $n$. Since $\epsilon$ is arbitrary, $a_n r(X_n-\theta)=o_p(1)$. Finally, the continuous mapping theorem gives $J_g(\theta)V_n\xrightarrow{d}J_g(\theta)Z$, and Slutsky’s theorem allows addition of the $o_p(1)$ remainder. This proves the claim. $\square$

Worked Applications

Worked Example: The Product of Two Sample Means

Let $\bar A_n$ and $\bar B_n$ be sample means from independent samples of size $n$, with Bernoulli success probabilities $p=0.4$ and $q=0.7$, respectively. Set $X_n=(\bar A_n,\bar B_n)$, $\theta=(p,q)$, and $g(u,v)=uv$. The joint central limit theorem gives $$ \sqrt n(X_n-\theta)\xrightarrow{d}Z, \qquad Z\sim N\left((0,0), \begin{pmatrix} 0.24&0\\ 0&0.21 \end{pmatrix}\right), $$ since $p(1-p)=0.24$, $q(1-q)=0.21$, and the samples are independent.

The domain is $U=\mathbb{R}^2$, so $g(X_n)$ is defined almost surely. The Jacobian is a $1\times2$ row matrix: $$ J_g(u,v)=\begin{pmatrix}v&u\end{pmatrix}, \qquad J_g(p,q)=\begin{pmatrix}0.7&0.4\end{pmatrix}. $$ Thus $$ \sqrt n(\bar A_n\bar B_n-pq)\xrightarrow{d}0.7Z_1+0.4Z_2. $$ This is a centered normal variable. Its variance is $$ 0.7^2(0.24)+0.4^2(0.21) =0.49(0.24)+0.16(0.21) =0.1176+0.0336 =0.1512. $$ Therefore the limit is $N(0,0.1512)$. The independence assumption is why the off-diagonal entries of the joint limiting covariance are zero.

Worked Example: A Log Ratio of Multinomial Probabilities

In $n$ independent trials, suppose the probabilities of three categories are $0.2$, $0.3$, and $0.5$, and let $N_1,N_2$ be the first two category counts. To keep the log ratio defined for every sample, define $$ X_n=\left(\frac{N_1+1}{n+3},\frac{N_2+1}{n+3}\right), \qquad \theta=(0.2,0.3). $$ Both coordinates of $X_n$ are strictly positive, so it lies in $U=(0,\infty)^2$ almost surely. The multinomial central limit theorem, together with the fact that this smoothing changes each raw proportion by $O(1/n)$, gives $$ \sqrt n(X_n-\theta)\xrightarrow{d}Z, \qquad Z\sim N\left((0,0), \begin{pmatrix} 0.16&-0.06\\ -0.06&0.21 \end{pmatrix}\right). $$ The diagonal terms are $0.2(0.8)=0.16$ and $0.3(0.7)=0.21$; the covariance is $-0.2(0.3)=-0.06$.

For $g(u,v)=\log(u/v)$, its Jacobian is $$ J_g(u,v)=\begin{pmatrix}1/u&-1/v\end{pmatrix}, \qquad J_g(\theta)=\begin{pmatrix}5&-10/3\end{pmatrix}. $$ Hence the limiting variance is $$ 25(0.16)+\left(\frac{100}{9}\right)(0.21) +2(5)\left(-\frac{10}{3}\right)(-0.06) =4+\frac{7}{3}+2 =\frac{25}{3}. $$ Consequently, $$ \sqrt n\left\{\log\left(\frac{N_1+1}{N_2+1}\right)-\log\left(\frac{0.2}{0.3}\right)\right\} \xrightarrow{d}N\left(0,\frac{25}{3}\right). $$ The covariance term contributes $2$; omitting it would give the wrong variance.

Worked Example: A Vector-Valued Transformation

Let $Z_n$ have distribution $N(0,\Sigma)$ for every $n$, where $$ \Sigma=\begin{pmatrix}1&0.5\\0.5&4\end{pmatrix}, \qquad X_n=\begin{pmatrix}2\\-1\end{pmatrix}+\frac{1}{\sqrt n}Z_n. $$ Take $g(u,v)=(u+v,uv)$, a map from $\mathbb{R}^2$ to $\mathbb{R}^2$. Here $\theta=(2,-1)$, $\sqrt n(X_n-\theta)=Z_n\xrightarrow{d}Z\sim N(0,\Sigma)$, and $g(\theta)=(1,-2)$.

The Jacobian and its value at $\theta$ are $$ J_g(u,v)= \begin{pmatrix} 1&1\\ v&u \end{pmatrix}, \qquad J_g(\theta)= \begin{pmatrix} 1&1\\ -1&2 \end{pmatrix}. $$ The theorem gives $\sqrt n\{g(X_n)-g(\theta)\}\xrightarrow{d}J_g(\theta)Z$, a centered bivariate normal vector. Its covariance, calculated by multiplying the matrices, is $$ J_g(\theta)\Sigma J_g(\theta)^{\mathsf T} = \begin{pmatrix}1&1\\-1&2\end{pmatrix} \begin{pmatrix}1&0.5\\0.5&4\end{pmatrix} \begin{pmatrix}1&-1\\1&2\end{pmatrix} = \begin{pmatrix}6&7.5\\7.5&15\end{pmatrix}. $$ For instance, the first diagonal entry is $1+4+2(0.5)=6$; the second is $1+16-2=15$. This example illustrates that the method applies when the transformed limit itself has several coordinates.

Normal Limits and the Role of Dependence

When $Z\sim N(0,\Sigma)$, the transformed limit $J_g(\theta)Z$ is also multivariate normal, possibly degenerate, with covariance $J_g(\theta)\Sigma J_g(\theta)^{\mathsf T}$. Indeed, every linear combination $c^{\mathsf T}J_g(\theta)Z$ is the linear combination $(J_g(\theta)^{\mathsf T}c)^{\mathsf T}Z$, hence is normal with mean zero and variance $$ c^{\mathsf T}J_g(\theta)\Sigma J_g(\theta)^{\mathsf T}c. $$ This verifies the stated covariance and normal form. It also explains why a joint limit, not just coordinate-by-coordinate limits, is needed: the off-diagonal elements of $\Sigma$ contribute to the transformed covariance.

Domain Conditions and a Common Pitfall

A frequent mistake is to apply $g$ to $X_n$ without checking that the transformation is defined there. For example, convergence $X_n\xrightarrow{p}\theta$ with $\theta$ inside an open domain $U$ implies $P(X_n\in U)\to1$, but it does not imply $P(X_n\in U)=1$ for every $n$. Thus, if $g$ is only defined on $U$, $g(X_n)$ may be undefined on an event of positive probability.

There are two standard ways to handle this. The theorem above assumes $X_n\in U$ almost surely for every $n$, as in the smoothed multinomial example. Alternatively, define the transformed statistic arbitrarily on $X_n\notin U$ when $P(X_n\notin U)\to0$. Such a change on events whose probabilities tend to zero does not alter a distributional limit: for any bounded continuous test function, the difference between its expectations under the two versions is bounded by twice its supremum norm times $P(X_n\notin U)$, which tends to zero. This qualification makes the transformed random variable well-defined while preserving the asymptotic conclusion.

The Jacobian must also be evaluated at the limiting point $\theta$, not at a random estimate or at a point chosen for algebraic convenience. The first-order result requires differentiability there; if the derivative vanishes or the map is not differentiable, a different approximation may be needed.

Key takeaway: If the scaled input vector has a joint limit and the transformation is differentiable at its limiting point, the scaled transformed difference converges to the Jacobian applied to that joint limit. Verify the rate, differentiability, and domain before using the formula.

Check Your Understanding

Try these questions without looking back at the examples.

In the multivariate delta method, what dimensions must the Jacobian have if the input has $k$ coordinates and the output has $m$ coordinates?
Why does $a_n(X_n-\theta)\xrightarrow{d}Z$ with $a_n\to\infty$ imply $X_n\xrightarrow{p}\theta$?
Suppose the limiting covariance of a two-coordinate input has a nonzero off-diagonal entry. Where does that dependence enter the variance of a scalar transformed limit?
Why is convergence of $X_n$ to an interior point not, by itself, enough to guarantee that $g(X_n)$ is defined for every $n$?
For $g(u,v)=(u+v,uv)$, write the Jacobian and evaluate it at $(3,2)$.

What You'll Learn