Multivariate Edgeworth Density Expansions

Why the Multivariate Correction Has a Tensor Form

For a vector-valued statistic, departures from normality can involve more than the skewness of each coordinate considered separately. Dependence can produce mixed third-order effects, such as a nonzero joint cumulant involving two coordinates and a third coordinate. The first multivariate Edgeworth correction collects all of these effects in a third-order cumulant tensor.

We develop the formula for sums of independent, identically distributed vectors with mean zero and identity covariance. Whitening extends the result to any positive definite covariance matrix. The approximation describes the density of a standardized sum; under suitable regularity conditions, its error is uniformly of order \(n^{-1}\).

Definitions and the First-Order Formula

Write \(\phi_d(x)=(2\pi)^{-d/2}\exp(-x^\mathsf{T}x/2)\) for the standard \(d\)-variate normal density. The multivariate Hermite polynomial of order three is \(H_{ijk}(x)=(-1)^3\phi_d(x)^{-1}\partial_i\partial_j\partial_k\phi_d(x)\), where \(\partial_i\) denotes differentiation with respect to \(x_i\). Direct differentiation gives

where \(\delta_{ij}=1\) if \(i=j\), and is zero otherwise. The first correction is a sum over every ordered triple of coordinate indices, including mixed-index terms.

The factor \(1/6\) is the reciprocal of \(3!\), arising from the cubic term in the characteristic-function expansion. The factor \(n^{-1/2}\) reflects the scaling of third cumulants for a standardized sum.

Derivation from the Characteristic Function

The derivation makes both the tensor contraction and the Hermite-polynomial sign explicit. Let \(\psi(u)=E[\exp(iu^\mathsf{T}Z)]\). For \(u\) near zero, finite third moments and the centering and covariance assumptions give

At the argument \(u=t/\sqrt n\), the characteristic function of \(W_n\) is \(\psi(t/\sqrt n)^n\). Keeping the normal term and the first cubic term yields

Under the inverse Fourier convention with kernel \(e^{-it^\mathsf{T}x}\), the inverse transform of \((it_i)(it_j)(it_k)e^{-t^\mathsf{T}t/2}\) is \((-1)^3\partial_i\partial_j\partial_k\phi_d(x)\). By the definition of \(H_{ijk}\), this equals \(H_{ijk}(x)\phi_d(x)\). Thus inversion gives the stated positive-sign correction.

A Uniform Density Result and Its Justification

A formal characteristic-function calculation alone does not ensure a density expansion. The following sufficient conditions make the inversion argument precise. They are stronger than needed in some settings, but permit a short proof with an explicit uniform error.

Proof. Finite fourth moments allow a Taylor expansion at zero with a fourth-order remainder. Taking the logarithm near \(\psi(0)=1\) gives \(\log\psi(u)=-\|u\|^2/2+C_3(u)+O(\|u\|^4)\), where \(C_3(u)=\frac16\sum_{i,j,k}\kappa_{ijk}(iu_i)(iu_j)(iu_k)\). Also, for some \(c>0\) and sufficiently small \(\delta>0\), \(|\psi(u)|\leq \exp(-c\|u\|^2)\) whenever \(\|u\|\leq\delta\). This follows by expanding \(|\psi(u)|^2\) at zero and using the positive definite covariance \(I_d\).

On \(\|t\|\leq\delta\sqrt n\), the logarithmic expansion at \(t/\sqrt n\) has exponent \(-\|t\|^2/2+n^{-1/2}C_3(t)+O(\|t\|^4/n)\). Shrink \(\delta\), if necessary, so the cubic term in absolute value is bounded by a sufficiently small constant times \(\|t\|^2\) throughout this region. Taylor's formula for the exponential, together with the Gaussian decay just established, then bounds the difference between \(\psi(t/\sqrt n)^n\) and \(e^{-\|t\|^2/2}[1+n^{-1/2}C_3(t)]\) by \[ \frac{C}{n}(\|t\|^4+\|t\|^6)e^{-c'\|t\|^2} \] for constants \(C,c'>0\). Its integral over \(\mathbb R^d\) is \(O(n^{-1})\).

For the remaining region, the assumptions imply that \(\sup_{\|u\|\geq\delta}|\psi(u)|\leq q<1\). Indeed, the limit condition bounds the characteristic function away from one outside a sufficiently large ball; on the intervening compact annulus continuity and the same condition give a maximum strictly below one. The latter follows because \(|\psi(u)|=1\) at a nonzero \(u\) would imply \(|\psi(ku)|=1\) for every positive integer \(k\), contradicting the limit condition. Consequently, for \(n\geq m\), \[ \int_{\|t\|>\delta\sqrt n}|\psi(t/\sqrt n)|^n\,dt \leq n^{d/2}q^{\,n-m}\int_{\mathbb R^d}|\psi(u)|^m\,du, \] which is exponentially small. The Gaussian comparison term also has an exponentially small integral on this region. Thus the full characteristic-function difference has \(L^1\) norm \(O(n^{-1})\).

The integrability assumption makes \(\psi(t/\sqrt n)^n\) integrable for \(n\geq m\). Fourier inversion therefore gives a continuous density for \(W_n\); the uniform bound on the difference of the inverse transforms is at most \((2\pi)^{-d}\) times the \(L^1\) difference just bounded. This proves the uniform expansion for that continuous version. Since two densities of the same distribution agree almost everywhere, the expansion holds almost everywhere for any other density version, but need not hold at every point for an arbitrarily modified version. \(\square\)

Worked Examples

Using the Formula Beyond Identity Covariance

Suppose \(Y\) has mean \(\mu\) and positive definite covariance matrix \(\Sigma\). Define the whitened vector \(Z=\Sigma^{-1/2}(Y-\mu)\), where \(\Sigma^{-1/2}\) is the symmetric inverse square root. Apply the formula to independent copies of \(Z\). For \(T_n=n^{-1/2}\sum_{r=1}^n(Y_r-\mu)\), the whitened statistic is \(\Sigma^{-1/2}T_n\). If \(T_n\) has a density, transform back by the Jacobian: its approximation at \(t\) is \(|\Sigma|^{-1/2}g_n(\Sigma^{-1/2}t)\), with the cumulant tensor in \(g_n\) calculated for \(Z\).

The correction matters when approximating probabilities or densities for vector statistics whose finite-sample distributions are not close enough to normal for the Gaussian term alone. The tensor form identifies which joint asymmetries contribute, while whitening keeps covariance effects distinct from higher-order departures from normality.

Common Pitfalls

• Keeping only marginal skewnesses. This drops mixed cumulants such as \(\kappa_{112}\). Include the full index sum unless independence or another model feature proves those entries vanish.
• Changing the Hermite sign. With \(H_{ijk}=(-1)^3\phi_d^{-1}\partial_i\partial_j\partial_k\phi_d\), the correction has a plus sign. The sign follows from the inverse Fourier transform of \((it_i)(it_j)(it_k)\), not from an informal analogy with a univariate formula.
• Treating the approximation as an exact density. The bracket in \(g_n\) can be negative in the tails. An Edgeworth approximation is an asymptotic expansion, not necessarily a nonnegative density.
• Ignoring density versions. A density can be changed on a set of measure zero without changing its distribution. A pointwise uniform theorem applies to the continuous version obtained by Fourier inversion; for another version the conclusion is only almost everywhere.