Multidimensional Edgeworth Expansions

• Multidimensional Edgeworth expansions are refinements of Gaussian approximations that incorporate higher-order cumulant corrections and Hermite polynomials to improve statistical accuracy.
• They utilize analytic mechanisms such as characteristic functions, Fourier inversion, and Stein’s method to systematically derive correction terms.
• The framework extends to high-dimensional, infinite-dimensional, and mixed-normal settings, offering precise approximations even under complex dependency structures.

Searching arXiv for recent and foundational papers on multidimensional Edgeworth expansions to ground the article in published work. Multidimensional Edgeworth expansions are refinements of Gaussian approximation in which a limiting normal law is supplemented by explicit higher-order correction terms built from cumulants, Hermite polynomials, derivatives of Gaussian densities, or operator-valued analogues. In the strict classical setting, the object is a vector $Z_n\in\mathbb{R}^k$ or a smooth functional of a $k$-dimensional sample mean; in more recent work, the same organizing idea appears for mixed-normal limits, path-space laws on Hilbert spaces, Poisson random measures on $\mathbb{R}^d$, high-dimensional linear statistics, random matrix cocycles, and finite-width neural network outputs (Podolskij et al., 2018, Jirak et al., 2020, Celli, 22 May 2026).

1. Classical form and scope

The classical multidimensional template starts from normalized sums of i.i.d. random vectors and expands their laws or expectations around a Gaussian reference. In the formulation quoted for random vectors, one writes

$$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$

or, for probabilities,

$$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$

where $Z$ is Gaussian with matching covariance and the coefficients depend on cumulants, especially the third and fourth cumulants (Privault, 2018).

A central modern formulation treats a statistic as a smooth function of a multivariate sample mean,

$W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$

with $Z_i\in\mathbb{R}^k$, $\bar Z=n^{-1}\sum_{j=1}^n Z_j$, and $H:\mathbb{R}^k\to\mathbb{R}$ sufficiently smooth. In that setting, the multidimensional content lies in the cumulants of the $k$0-dimensional mean, even when the final statistic is scalar (Wei et al., 3 Nov 2025). This same pattern reappears for many familiar objects: sample variance, covariance, Pearson’s correlation, ratio statistics, and Z-score type functionals are all treated as smooth functions of multivariate means in the cited work (Wei et al., 3 Nov 2025).

A recurrent source of ambiguity is that “multidimensional” does not always mean that the target approximation itself is $k$1-valued. In several important papers, the approximated variable is scalar, while multidimensionality enters through the ambient space, the input array, the parameter block, or the geometry of the underlying random measure. The Poisson–Malliavin setting on $k$2 and panel-data expansions built from $k$3 arrays are explicit examples of this broader usage (Privault, 2018, Gao et al., 2022).

2. Analytic mechanisms

Classical proofs rely on characteristic functions, cumulant expansions, and Fourier inversion. That structure remains visible even in recent work: cumulants are encoded through Taylor expansions of $k$4, correction polynomials are generated from cumulant series, and inversion yields Gaussian densities multiplied by Hermite-type terms (Wei et al., 3 Nov 2025). In weakly dependent stationary processes, the transition from a Berry–Esseen bound to a second-order Edgeworth expansion is described through a Berry–Esseen characteristic $k$5, with a parallel integrated characteristic $k$6 for the Wasserstein metric $k$7 (Jirak et al., 2020).

Stein’s method provides an alternative route. For the one-dimensional normal target, the Stein equation

$k$8

is combined with Stein identities and, in the continuous case, Stein kernels. Repeated use of Stein equations produces the two-term correction

$k$9

with $\mathbb{R}^d$0 remainder for bounded test functions, including indicators of half-lines (Fang et al., 2022). In infinite-dimensional Gaussian settings, Stein’s method is coupled with Malliavin calculus to produce functional Edgeworth expansions on $\mathbb{R}^d$1, with correction terms expressed as contractions of higher Fréchet derivatives against explicit tensors or covariance perturbations (Coutin et al., 2014).

For mixed-normal limits, Malliavin calculus enters differently. In the Euler approximation of a diffusion, the relevant object is the joint vector $\mathbb{R}^d$2, where $\mathbb{R}^d$3 is the principal statistic and $\mathbb{R}^d$4 is an auxiliary vector. The approximating density is built from a conditional Gaussian core $\mathbb{R}^d$5 and random symbols $\mathbb{R}^d$6 and $\mathbb{R}^d$7, which encode adaptive and anticipative corrections through Malliavin derivatives of the random covariance and of the auxiliary variable (Podolskij et al., 2018). In mod-$\mathbb{R}^d$8 theory, by contrast, the central analytic object is an exponentially normalized Laplace transform $\mathbb{R}^d$9 converging to a random analytic limit $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$0, and the Edgeworth coefficients are generated from deterministic cumulants $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$1 and random cumulants $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$2 (Kabluchko et al., 2016).

3. Validity conditions and obstructions

The non-lattice or Cramér condition remains a central theme, but its correct formulation depends strongly on dimension and structure. In higher dimensions, the paper on Bernoulli weighted means defines a random vector $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$3 to be semi-lattice if there exists a nonzero vector $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$4 such that $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$5 is concentrated on a one-dimensional lattice. It then proves the equivalence

$$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$6

thereby identifying the correct geometric obstruction in dimensions $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$7 (Cauvin, 2022). A recurrent misconception is therefore that “not supported on a full lattice” is the relevant condition in higher dimensions; the cited result shows that the obstruction is the existence of a lattice projection.

For smooth functions of sample means, the general partial Cramér’s condition (GPCC) enlarges the admissible class beyond the original partial Cramér condition. Instead of requiring a one-coordinate conditional characteristic function to decay, GPCC permits a block $$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$8 and asks for

$$\mathbb{E}[f(S_n)] = \mathbb{E}[f(Z)] + n^{-1/2}A_1(f) + n^{-1}A_2(f) + \cdots,$$9

where $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$0 is the conditional joint characteristic function of the first $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$1 coordinates given the remaining ones (Wei et al., 3 Nov 2025). This is what allows vectors such as $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$2 and $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$3 to enter Edgeworth theory when standard partial Cramér conditions fail (Wei et al., 3 Nov 2025).

In weakly dependent stationary processes, the cited “sharp transition” results show that a central limit theorem alone is not enough: a second-order Edgeworth expansion with $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$4 Kolmogorov error is valid if and only if the Berry–Esseen characteristic $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$5 is itself $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$6 (Jirak et al., 2020). In that framework, no non-lattice or Cramér condition is imposed, and the characteristic $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$7 or its Wasserstein analogue $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$8 becomes the decisive object (Jirak et al., 2020).

For symmetric statistics, the decisive assumptions are different again. The linear Hoeffding projection must satisfy a localized Cramér condition, while the quadratic component must have a genuinely non-reducible orthogonal part. The resulting two-term Edgeworth expansion then has remainder $$P(S_n\in A)=\int_A \phi_\Sigma(x)\Big(1+n^{-1/2}P_1(x)+n^{-1}P_2(x)+\cdots\Big)\,dx+o(n^{-r/2}),$$9 for a large class of asymptotically linear symmetric statistics (Götze et al., 2021).

4. Dependent data, Poisson geometry, and studentization

One important line of work replaces i.i.d. vectors by structured randomness. For compensated Poisson integrals over $Z$0,

$Z$1

cumulant-type operators

$Z$2

generate an Edgeworth-type expansion in which the third cumulant is the first non-Gaussian correction (Privault, 2018). In the deterministic case $Z$3, the third cumulant is

$Z$4

and the refinement replaces an $Z$5-norm bound by the signed cubic integral. The cited radial example shows that when the third cumulant vanishes, the convergence rate can improve from the Berry–Esseen $Z$6 rate to an $Z$7 rate (Privault, 2018). This is a concrete instance in which a multidimensional spatial base space produces a sharper scalar Gaussian approximation.

In panel data with serial dependence and cross-sectional dependence, the scalar statistic

$Z$8

admits a second-order cdf expansion

$Z$9

with uniform remainder $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$0 under the stated dependence assumptions (Gao et al., 2022). The same paper derives a Berry–Esseen bound of order $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$1 generally and $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$2 under additional weak cross-sectional dependence (Gao et al., 2022).

Studentization under weak dependence produces a distinctly nonclassical structure. For strongly mixing data, the Edgeworth expansion of a studentized statistic is a superposition of three series: one in powers of $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$3, one in powers of $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$4 arising from the standard error of the studentizing factor, and one in powers of the bias of the studentizing factor (Lahiri, 2010). The cited analysis emphasizes that the nuisance dimension $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$5 of lag-covariance estimators grows with $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$6, so the effective multidimensionality resides in the studentizer as much as in the original statistic (Lahiri, 2010).

For symmetric statistics, Hoeffding’s decomposition supplies the expansion mechanism. A statistic is decomposed into linear, quadratic, cubic, and higher-order canonical components, and the resulting Edgeworth coefficients involve mixed cumulant-like quantities built from the kernels $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$7, $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$8, and $W_n=n^{1/2}\bigl(H(\bar Z)-H(\mu)\bigr),$9. Under the stated moment, Cramér, and non-reducibility conditions, the remainder is $Z_i\in\mathbb{R}^k$0 (Götze et al., 2021).

5. Functional, geometric, and profile-level extensions

A major extension of the subject replaces finite-dimensional vectors by function-space laws. In the Stein-based functional theory, Brownian paths are embedded into a fractional Sobolev–Besov space $Z_i\in\mathbb{R}^k$1, transferred isometrically to $Z_i\in\mathbb{R}^k$2, and compared with an infinite-dimensional Gaussian measure $Z_i\in\mathbb{R}^k$3. For a normalized compensated Poisson process, the law $Z_i\in\mathbb{R}^k$4 satisfies an expansion of the form

$Z_i\in\mathbb{R}^k$5

with remainder $Z_i\in\mathbb{R}^k$6 for $Z_i\in\mathbb{R}^k$7 (Coutin et al., 2014). For the linear interpolation of Brownian motion, the correction is instead expressed through powers of the covariance difference $Z_i\in\mathbb{R}^k$8, and the two expansions are explicitly described as “rather different” (Coutin et al., 2014). This shows that “multidimensional” can mean genuinely infinite-dimensional.

The Euler approximation of a diffusion provides a mixed-normal functional analogue. Here the finite-dimensional vector of errors at times $Z_i\in\mathbb{R}^k$9 is embedded into a joint vector $\bar Z=n^{-1}\sum_{j=1}^n Z_j$0, where $\bar Z=n^{-1}\sum_{j=1}^n Z_j$1 contains random scaling variables such as $\bar Z=n^{-1}\sum_{j=1}^n Z_j$2 and possibly other functionals. The approximation is not a fixed Gaussian density but a mixture $\bar Z=n^{-1}\sum_{j=1}^n Z_j$3 corrected by differential operators weighted by the random symbols $\bar Z=n^{-1}\sum_{j=1}^n Z_j$4 and $\bar Z=n^{-1}\sum_{j=1}^n Z_j$5 (Podolskij et al., 2018). A plausible implication is that, in mixed-normal settings, the role played by cumulants in classical Edgeworth theory is assumed by random symbols adapted to the limiting conditional covariance.

In products of positive random matrices, the expanded variable is scalar,

$\bar Z=n^{-1}\sum_{j=1}^n Z_j$6

but the underlying system is a multidimensional random matrix cocycle acting on the positive cone and projective space. Under assumptions A2, A3, and a third logarithmic moment, the expansion

$\bar Z=n^{-1}\sum_{j=1}^n Z_j$7

holds uniformly on compact sets (Xiao et al., 2022). Here the third cumulant $\bar Z=n^{-1}\sum_{j=1}^n Z_j$8 appears together with boundary corrections $\bar Z=n^{-1}\sum_{j=1}^n Z_j$9 and $H:\mathbb{R}^k\to\mathbb{R}$0, which are specific to the cocycle geometry (Xiao et al., 2022).

Mod-$H:\mathbb{R}^k\to\mathbb{R}$1 theory supplies yet another direction. For profiles $H:\mathbb{R}^k\to\mathbb{R}$2, exponentially normalized Laplace transforms $H:\mathbb{R}^k\to\mathbb{R}$3 converge almost surely to a random analytic function $H:\mathbb{R}^k\to\mathbb{R}$4, and one obtains an almost-sure Edgeworth expansion uniform in both $H:\mathbb{R}^k\to\mathbb{R}$5 and $H:\mathbb{R}^k\to\mathbb{R}$6 on compact intervals (Kabluchko et al., 2016). The corresponding applications to one-split branching random walks and random trees yield refined asymptotics for profile, mode, width, and occupation numbers (Kabluchko et al., 2016).

6. High-dimensional regimes, misconceptions, and open directions

High-dimensional linear statistics in Hilbert spaces fit naturally into projection-based Edgeworth theory. For a stationary $H:\mathbb{R}^k\to\mathbb{R}$7-valued process and linear statistics $H:\mathbb{R}^k\to\mathbb{R}$8, the one-dimensional expansions hold uniformly over directions $H:\mathbb{R}^k\to\mathbb{R}$9 satisfying the stated tail conditions, and no non-lattice or Cramér condition is required (Jirak et al., 2020). This suggests that genuinely high-dimensional approximations can sometimes be reduced to a family of one-dimensional projections, provided the dependence and coefficient tails are controlled.

The most explicit recent multivariate result in the supplied material concerns finite-width neural networks. For network outputs evaluated at a finite set of inputs, the output vector $k$00 is conditionally Gaussian, and its Edgeworth approximation is written directly as a multivariate signed density involving products of one-dimensional Hermite polynomials and expectations of products of entries of the relative covariance matrix $k$01 (Celli, 22 May 2026). The paper proves multidimensional Edgeworth expansions of arbitrary order $k$02 with total variation error of order $k$03, together with matching lower bounds. It further shows that replacing the prior by its Edgeworth approximation perturbs the Bayesian posterior by at most $k$04 in total variation (Celli, 22 May 2026). In this literature, high dimensionality is therefore not merely an asymptotic nuisance; it is the native state space of the approximation.

A different high-dimensional application appears in sphericity testing under two-step monotone incomplete data. There, the scalar likelihood-ratio statistic has a standardized Edgeworth approximation

$k$05

and the paper derives computable bounds for

$k$06

through an explicit decomposition of the Fourier integral into $k$07, $k$08, and $k$09 terms (Sato et al., 31 Mar 2026). The simulations reported there show that the proposed Edgeworth expansion is more accurate than the existing asymptotic chi-square approximations in high-dimensional settings (Sato et al., 31 Mar 2026).

Several open directions recur across the cited literature. A fully vector-valued Poisson Stein–Malliavin Edgeworth theory is identified as a natural extension of the scalar Poisson results (Privault, 2018). In panel models, the Edgeworth expansion itself is stated only for scalar statistics, even though the final parameter is multivariate normal, and the multivariate extension is left implicit (Gao et al., 2022). For Bernoulli weighted means, the main expansion remains one-dimensional, while the mixed semi-lattice versus non-semi-lattice geometry is presented as a first step toward a genuine multidimensional theory for mixed discrete/continuous vectors (Cauvin, 2022). For Euler schemes, the paper is scalar in the diffusion state variable but already fully multidimensional in the statistic vector $k$10, suggesting that higher-dimensional SDEs would require stronger non-degeneracy arguments rather than a different conceptual framework (Podolskij et al., 2018).

Taken together, these works show that multidimensional Edgeworth expansions are not a single theorem but a family of higher-order Gaussian approximation principles whose concrete form depends on what counts as “dimension”: vector-valued sums, blocks of sample means, geometric base spaces, path spaces, projective state spaces, or high-dimensional output layers. The shared structure is the same throughout: a Gaussian leading term, correction terms organized by cumulants or their analogues, and a remainder whose validity depends on explicit analytic, geometric, or probabilistic conditions (Privault, 2018, Jirak et al., 2020, Celli, 22 May 2026).