Edgeworth Coefficients for Standard Multivariate Estimates

1. Introduction and Summary

Suppose that $\widehat{w}$ is a standard estimate, as defined in §2, of an unknown parameter $w\in R^q$ of a statistical model, based on a sample of size n. For example $\widehat{w}$ may be a smooth function of a sample mean, or a smooth functional of an empirical distribution. A smooth function of a standard estimate is also a standard estimate: see Withers (2024). §2 summarises the multivariate Edgeworth expansions of Withers and Nadarajah (2010) for the distribution and density of $X_n=n^{1/2}(\widehat{w}-w)$ in powers of $n^{-1/2}$ about the multivariate normal in terms of the Edgeworth coefficients, the $P_r$-coefficients of (11). For $r\ge 1$, these $P_r$ are needed for the $r+1$st term of the multivariate Edgeworth expansions. They are Bell polynomials in the cumulant coefficients of (5). They are given for $r=1$ by (12) and for $r=2,3$ by (12)–(14) in terms of the symmetrizing operator $\mathcal{S}$, and explicitly in Appendix A. §3 derives expansions on ellipsoidal and hyperrectangular sets.

When $q=2$, §4 simplifies these Edgeworth coefficients using dual notation. Examples include the distribution and density of a sample mean, and of bivariate entangled gamma random variables.

§5 and §6 give conclusions and discussion, and suggests some future directions. Appendix B gives explicitly the bivariate Hermite polynomials needed for bivariate Edgeworth expansions to $O\left(n^{-2}\right)$.

Univariate estimates. Suppose that $\widehat{w}$ is a standard estimate of $w\in R$ with respect to n, (typically the sample size). That is, $\widehat{w}$ is non-lattice,

$E\widehat{w}\to w$ as $n\to \infty$, and its rth cumulant can be expanded as

$$\begin{array}{cc}& {\kappa }_r\left(\widehat{w}\right)\approx \sum _{j=r-1}^{\infty }{n}^{-j}{a}_{rj}\mathrm{for}r\ge 1,\hfill \end{array}$$

(1)

where the cumulant coefficients $a_{rj}$ may depend on n but are bounded as $n\to \infty$, and $a_{21}$ is bounded away from 0. Here and below ≈ indicates an asymptotic expansion that need not converge. So (1) holds in the sense that

$${\kappa }_r\left(\widehat{w}\right)=\sum _{j=r-1}^{I-1}{n}^{-j}{a}_{rj}+O\left(n^{-I}\right)\mathrm{for}I\ge r\ge 1,$$

where $y_n=O\left(x_n\right)$ means that $y_n/x_n$ is bounded in n. Withers (1984) replaced the artifical assumptions of Cornish and Fisher (1937) and Fisher and Cornish (1960) by (1), and gave the distribution, density and quantiles of

$$U_n={(n/a_{21})}^{1/2}(\widehat{w}-w)$$

as asymptotic expansions in powers of $n^{-1/2}$:

$$\begin{array}{cc}& P_n\left(u\right)=Prob.(U_n\le u)\approx \mathsf{\Phi }\left(u\right)-\varphi \left(u\right)\sum _{r=1}^{\infty }{n}^{-r/2}{h}_r\left(u\right),\hfill \end{array}$$

(2)

(3)

(4)

where $Prob.\left(A\right)$ is the probability that A is true,

$$\mathsf{\Phi }\left(u\right)=Prob.(N\le u)={\int }_{-\infty }^u\varphi \left(u\right)du,\varphi \left(u\right)={\left(2\mathsf{\pi }\right)}^{-1/2}{e}^{-u^2/2},N\sim \mathcal{N}(0,1)$$

is a unit normal random variable with even moments $E{N}^{2k}=1.3\dots (2k-1)$, and $h_r\left(u\right),{\overline{h}}_r\left(u\right),f_r\left(u\right),g_r\left(u\right)$ are polynomials in u and also in the standardized cumulant coefficients $A_{ri}=a_{ri}/a_{21}^{r/2}$. For example

$$\begin{array}{cc}& h_1\left(u\right)=f_1\left(u\right)=g_1\left(u\right)=A_{11}+A_{32}{H}_2/6,{\overline{h}}_1\left(u\right)=A_{11}{H}_1+A_{32}{H}_3/6,\hfill \\ & h_2\left(u\right)=(A_{11}^2+A_{22})H_1/2+(A_{11}{A}_{32}+A_{43}/4)H_3/6+A_{32}^2{H}_5/72,\hfill \end{array}$$

(5)

where $H_k$ is the kth Hermite polynomial. By Withers (2000), for $I=\sqrt{-1}$,

$$\begin{array}{cc}& H_k=H_k\left(u\right)=\varphi {\left(u\right)}^{-1}{(-d/du)}^k\varphi \left(u\right)=E{(u+IN)}^k\mathrm{for}k\ge 0:\hfill \\ & H_0=1,H_1=u,H_2=u^2-1,H_3=u^3-3u,H_4=u^4-6u^2+3.\hfill \end{array}$$

(6)

(A1) gives $H_k$ for $k\le 9$. Since $h_r\left(u\right)$ is even/odd for r odd/even,

$$\begin{array}{cc}& Prob.\left(\right|U_n|\le u)\approx \mathsf{\Phi }\left(u\right)-2\varphi \left(u\right)\sum _{r=1}^{\infty }{n}^{-r}{h}_{2r}\left(u\right)\mathrm{for}u>0.\hfill \end{array}$$

(7)

From (4), it follows that

$$\begin{array}{cc}& U_n\approx u+\sum _{r=1}^{\infty }{n}^{-r/2}{g}_r\left(N\right)\mathrm{where}N\sim \mathcal{N}(0,1).\hfill \end{array}$$

(8)

For a discussion of when to truncate (2)–(4), see Withers and Nadarajah (2012a).

NOTE 1.1.

Edgeworth’s expansion was for $\widehat{w}$ the mean of n independent identically distributed random variables on R from a distribution with rth cumulant ${\kappa }_r$. So (1) holds with $a_{r,r-1}={\kappa }_r$, and other $a_{ri}=0$. An explicit formula for its general term was given in Withers and Nadarajah (2009).

For examples of the many applications of Cornish-Fisher expansions and the extensions of Hill and Davis (1968), see Withers and Nadarajah (2011, 2014, 2015). Applications to finance include Simonato (2011) and Zhang et al. (2011). For an application to Rayleigh fading amplitudes, see Withers and Nadarajah (2008). TianLi et al. (2009) used them for GPS accuracy. Winterbottom (1974, 1980, 1984) and Abdel-Wahed and Winterbottom (1983) successfully used them for system reliability, even though binomial random variables fall on a lattice.

Now suppose that ${\widehat{w}}_{*}$ is another standard estimate of w with the same asymptotic variance $a_{21}/n$. Denote its standardized cumulant coefficients by $A_{ri*}$. Suppose that

$$U_{n*}={(n/a_{21})}^{1/2}({\widehat{w}}_{*}-w)\sim {\mathsf{\Phi }}_n\left(u\right)={\int }_{-\infty }^u{\varphi }_n\left(u\right)du\to \mathsf{\Phi }\left(u\right)\mathrm{as}n\to \infty .$$

Then one can expand $P_n\left(u\right)$ about ${\mathsf{\Phi }}_n\left(u\right)$ rather than $\mathsf{\Phi }\left(u\right)$. The above expressions for $P_n\left(u\right),h_k\left(u\right)$ become

$$\begin{array}{cc}& P_n\left(u\right)=Prob.(U_n\le u)\approx {\mathsf{\Phi }}_n\left(u\right)-{\varphi }_n\left(u\right)\sum _{r=1}^{\infty }{n}^{-r/2}{h}_{rn}\left(u\right),\mathrm{where}\hfill \\ & h_{1n}\left(u\right)=A_{110}+A_{320}{H}_{2n}/6,\hfill \\ & h_{2n}\left(u\right)=(A_{110}^2+A_{220})H_{1n}/2+(A_{110}{A}_{320}+A_{430}/4)H_{3n}/6+A_{320}^2{H}_{5n}/72,\hfill \\ & A_{rj0}=A_{rj}-A_{rj*},\mathrm{and}{H}_{kn}=H_{kn}\left(u\right)={\varphi }_n{\left(u\right)}^{-1}{(-d/du)}^k{\varphi }_n\left(u\right).\hfill \end{array}$$

So if for example we choose ${\widehat{w}}_{*}$ so that $A_{110}=A_{320}=0$, then

$$h_{1n}\left(u\right)=0,h_{2n}\left(u\right)=A_{220}{H}_{1n}/2+A_{430}{H}_{3n}/24,$$

and the number of terms in the analogues of $h_r,f_r,g_r$ are greatly reduced. The disadvantage is that $H_{kn}\left(u\right)$ is more complicated than $H_k\left(u\right)$. See Withers and Nadarajah (2010). For expansions about a matching Student’s, gamma, or F distribution, see Withers and Nadarajah (2011, 2012b, 2014).

2. Multivariate Edgeworth Expansions

Ordinary Bell polynomials. For a sequence $e=(e_1,e_2,\dots )$ from R, the partial ordinary Bell polynomial ${\tilde{B}}_{rs}={\tilde{B}}_{rs}\left(e\right)$, is defined by the identity

$$\begin{array}{cc}& \mathrm{for}s=0,1,2,\dots ,\mathrm{and}z\in R,S^s=\sum _{r=s}^{\infty }{z}^r{\tilde{B}}_{rs}\left(e\right)\mathrm{where}S=\sum _{r=1}^{\infty }{z}^r{e}_r.\hfill \end{array}$$

(1)

(2)

where ${\delta }_{00}=1,{\delta }_{r0}=0$ for $r\ne 0.$ They are tabled on p309 of Comtet (1974). (The partial exponential Bell polynomials are not used in this paper.) The complete ordinary Bell polynomial, ${\tilde{B}}_r\left(e\right)$ is defined in terms of S of (1) by

$$\begin{array}{cc}& e^S=\sum _{r=0}^{\infty }{z}^r{\tilde{B}}_r\left(e\right).So{\tilde{B}}_0\left(e\right)=1\mathrm{and}\mathrm{for}r\ge 1,{\tilde{B}}_r\left(e\right)=\sum _{s=1}^r{\tilde{B}}_{rs}\left(e\right)/s!:\hfill \end{array}$$

(3)

(4)

Multivariate estimates. Suppose that $\widehat{w}$ is a standard estimate of

$w\in R^q$ with respect to n. That is, $\widehat{w}$ is non-lattice, $E\widehat{w}\to w$ as $n\to \infty$, and for $r\ge 1$, $1\le i_1,\dots ,i_r\le q,$ the rth order cumulants of $\widehat{w}$ can be expanded as

$$\begin{array}{c}\hfill {\overline{k}}^{1-r}=k^{i_1\dots i_r}=\kappa ({\widehat{w}}^{i_1},\dots ,{\widehat{w}}^{i_r})\approx \sum _{d=r-1}^{\infty }{n}^{-d}{\overline{k}}_d^{1-r}\mathrm{where}{\overline{k}}_d^{1-r}=k_d^{i_1\dots i_r},\end{array}$$

(5)

and the cumulant coefficients ${\overline{k}}_d^{1-r}$ may depend on n but are bounded as $n\to \infty$. So the bar replaces each $i_k$ by k. The use of $i_k$ is reserved for this purpose. For example ${\overline{k}}_0^1=w^{i_1}$ and ${\overline{k}}_1^{12}=k_1^{i_1{i}_2}.$ We use this bar notation repeatedly to avoid double subscripts $i_k$. As $n\to \infty$,

$$\begin{array}{cc}& X_n=n^{1/2}(\widehat{w}-w)\stackrel{\mathcal{L}}{\to }X={\mathcal{N}}_q(0,V)\mathrm{for}V=\left({\overline{k}}_1^{12}\right),q\times q,\hfill \end{array}$$

(6)

the multivariate normal on $R^q$, with density and distribution

$$\begin{array}{cc}& {\varphi }_V\left(x\right)={\left(2\mathsf{\pi }\right)}^{-q/2}{\left(detV\right)}^{-1/2}exp(-x^{\prime }{V}^{-1}x/2),{\mathsf{\Phi }}_V\left(x\right)={\int }_{-\infty }^x{\varphi }_V\left(x\right)dx.\hfill \end{array}$$

(7)

V may depend on n, but we assume that $detV$ is bounded away from 0. Set

(8)

(9)

where here and below, we use the tensor summation convention of implicitly summing each $i_k$ over its range $1,\dots ,q$. For example

$$\begin{array}{cc}& {\overline{b}}_2^1{\overline{t}}_1={\overline{k}}_1^1{\overline{t}}_1=\sum _{i_1=1}^q{k}_1^{i_1}{t}_{i_1}\mathrm{and}{\overline{b}}_2^{12}{\overline{t}}_1{\overline{t}}_2={\overline{k}}_1^{12}{\overline{t}}_1{\overline{t}}_2=\sum _{i_1,i_2=1}^q{k}_1^{i_1{i}_2}{t}_{i_1}{t}_{i_2}.\hfill \end{array}$$

(10)

(11)

where the rth Edgeworth coefficient, ${\overline{P}}_r^{1-k}=P_r^{i_1\dots i_k}$, is a function of the ${\overline{k}}_d^{1-r}$. One could use unsymmetrized ${\overline{P}}_r^{1-k}$. For example by (2) ${\tilde{B}}_2\left(e\left(t\right)\right)$ needs the cross term in $e_1{\left(t\right)}^2/2,{\overline{k}}_1^4{\overline{k}}_2^{1-3}{\overline{t}}_1{\overline{t}}_2{\overline{t}}_3{\overline{t}}_4/6$. So we could use ${\overline{k}}_1^4{\overline{k}}_2^{1-3}/6$ in ${\overline{P}}_2^{1-4}$. But as (3) below illustrates, there is a big advantage in making ${\overline{P}}_r^{1-k}$ symmetric in $i_1,\dots ,i_k$ using the operator $\mathcal{S}$ that symmetrizes over $i_1,\dots ,i_k$. So the ${\overline{P}}_r^{1-k}$ for $r=1,2,3$ that we need are

(12)

(13)

(14)

These formulas are mostly new. The terms involving $\mathcal{S}$ are given explicitly for the first time in Appendix A. By Withers and Nadarajah (2010), the distribution and density of $X_n$ of (6), can be expanded as

(15)

(16)

(17)

(18)

(19)

(20)

is the multivariate Hermite polynomial. By Withers (2020), for $I=\sqrt{-1}$,

(21)

(22)

(23)

(24)

(25)

(26)

for ${\overline{\mu }}^{1-2k}$ of (8). These give $p_1\left(x\right)$ and $p_2\left(x\right)$, and so $p_{X_n}\left(x\right)$ of (15) to $O\left(n^{-3/2}\right)$. For ${\overline{H}}^{1-k}$ for $7\le k\le 9$ needed for $p_3\left(x\right)$ when $q=2$, see Appendix B. $P_{rk}\left(x\right)$ is just ${\tilde{p}}_{rk}$ with ${\overline{H}}^{1-k}$ replaced by ${\overline{H}}_{*}^{1-k}$ of (19). For example,

$$\begin{array}{cc}& {\overline{H}}_{*}^1={\overline{J}}^1,{\overline{H}}_{*}^{12}={\overline{J}}^{12}-{\overline{V}}^{12}{\mathsf{\Phi }}_V\left(x\right),{\overline{H}}_{*}^{1-3}={\overline{J}}^{123}-\sum ^3{\overline{J}}^1{\overline{V}}^{23}\mathrm{where}\hfill \\ & {\overline{J}}^{1-k}={\overline{J}}^{1-k}(x,V)=E{\overline{Y}}_1\dots {\overline{Y}}_{k}I(X\le x)={\overline{V}}^{1,k+1}\dots {\overline{V}}^{k,2k}{\overline{M}}^{k+1-2k},\hfill \\ & \mathrm{and}{\overline{M}}^{a-b}={\overline{M}}^{a-b}(x,V)=E{\overline{X}}_a\dots {\overline{X}}_{b}I(X\le x)={\int }_{-\infty }^x{\overline{x}}_a\dots {\overline{x}}_b{\varphi }_V\left(x\right)dx.\hfill \end{array}$$

We call $\left\{{\overline{M}}^{a-b}(x,V)\right\}$the partial moments of${\mathsf{\Phi }}_V\left(x\right).$ So by (15), the Edgeworth expansions to $O\left(n^{-2}\right)$ for the distribution and density of $X_n$ of (6) about those of $X={\mathcal{N}}_q(0,V)$, are given by

$$\begin{array}{cc}& P_1\left(x\right)=e_1(-\partial /\partial x){\mathsf{\Phi }}_V\left(x\right)=\sum _{r=1}^3{\overline{b}}_{r+1}^{1-r}{\overline{O}}^{1-r}{\mathsf{\Phi }}_V\left(x\right)/r!=\sum _{k=1,3}{P}_{1k}\left(x\right),\hfill \\ & \mathrm{where}{P}_{11}\left(x\right)={\overline{k}}_1^1(-{\overline{\partial }}_1){\mathsf{\Phi }}_V\left(x\right),P_{13}\left(x\right)={\overline{k}}_2^{1-3}{\overline{O}}^{1-3}{\mathsf{\Phi }}_V\left(x\right)/6,\hfill \\ & p_1\left(x\right)={\overline{k}}_1^1(-{\overline{\partial }}_1){\varphi }_V\left(x\right)+{\overline{k}}_2^{1-3}{\overline{O}}^{1-3}{\varphi }_V\left(x\right)/6,\hfill \end{array}$$

(27)

$$\begin{array}{cc}& P_2\left(x\right)=\sum _{k=2,4,6}{P}_{2k}\left(x\right),{\tilde{p}}_2\left(x\right)=\sum _{k=2,4,6}{\tilde{p}}_{2k},\hfill \end{array}$$

(28)

$$\begin{array}{cc}& P_3\left(x\right)=\sum _{k=1,3,5,7,9}{P}_{3k}\left(x\right),{\tilde{p}}_3\left(x\right)=\sum _{k=1,3,5,7,9}{\tilde{p}}_{3k},\hfill \end{array}$$

(29)

for ${\tilde{p}}_{rk}$ and $P_{rk}$ of (18). Each has $q^k$ terms, but many are duplicates, as it is of great advantage to make ${\overline{P}}_r^{1-k}$ symmetric in $i_1,\dots ,i_k$. The density of $X_n$ relative to its asymptotic value is

$$p_{X_n}\left(x\right)/{\varphi }_V\left(x\right)\approx 1+\sum _{r=1}^{\infty }{n}^{-r/2}{\tilde{p}}_r\left(x\right)=1+n^{-1/2}{\tilde{p}}_1\left(x\right)+O\left(n^{-1}\right)\mathrm{for}x\in R^q,$$

for ${\tilde{p}}_r\left(x\right)$ of (17). So $n^{-1/2}{\tilde{p}}_1\left(x\right)$ is a simple measure of the accuracy of the Central Limit Theorem (CLT) approximation. An exception is

Example 1.

If the distribution of $\widehat{w}$ is symmetric about w, then for r odd, $p_r\left(x\right)={\tilde{p}}_r\left(x\right)=P_r\left(x\right)=0$ and by (13), ${\tilde{p}}_{26}\left(x\right)=0$. So,

(30)

(31)

since ${\overline{P}}_2^{12}={\overline{k}}_2^{12}/2,{\overline{P}}_2^{1-4}={\overline{k}}_3^{1-4}/24.$ In this case $n^{-1}{\tilde{p}}_2\left(x\right)$ is a measure of the accuracy of the CLT approximation. Also,

$$\begin{array}{c}\hfill Prob.(X_n\le x)={\mathsf{\Phi }}_V\left(x\right)+n^{-1}{P}_2\left(x\right)+O\left(n^{-2}\right)\mathrm{where}{P}_2=P_{22}+P_{24}.\end{array}$$

(32)

Example 2.

Let $\widehat{w}$ be the sample mean from a distribution with finite cross cumulants ${\overline{\kappa }}^{1-r},r\ge 1$. Then $E\widehat{w}=w$, and only the leading coefficient in (5) is non-zero. ${\overline{k}}_1^1={\overline{k}}_2^{12}={\overline{k}}_2^1={\overline{k}}_3^{123}=0\Rightarrow$

${\tilde{p}}_{11}={\tilde{p}}_{22}={\tilde{p}}_{31}={\tilde{p}}_{33}=0$. For $r=1,2,3$, the non-zero $P_r^{1-k}$ are

$$\begin{array}{cc}& {\overline{P}}_1^{1-3}={\overline{\kappa }}^{1-3}/3!,{\overline{P}}_2^{1-4}={\overline{\kappa }}^{1-4}/4!,{\overline{P}}_2^{1-6}=\mathcal{S}{\overline{\kappa }}^{1-3}{\overline{\kappa }}^{4-6}/72,,\hfill \\ & {\overline{P}}_3^{1-5}={\overline{\kappa }}^{1-5}/5!,{\overline{P}}_3^{1-7}=\mathcal{S}{\overline{\kappa }}^{123}{\overline{\kappa }}^{4-7}/144,{\overline{P}}_3^{1-9}=\mathcal{S}{\overline{\kappa }}^{1-3}{\overline{\kappa }}^{4-6}{\overline{\kappa }}^{7-9}/6^3.\hfill \end{array}$$

By (27)–(29), for $1\le r\le 3,P_r\left(x\right)$ and ${\tilde{p}}_r\left(x\right)$ have only r terms.

NOTE 2.1.

For $H_k\left(x\right)$ the univariate Hermite polynomial of (6), and

$1\le j\le q$, we define the jth marginal Hermite polynomial as

$$\begin{array}{cc}& H^{j^k}=E{(y_j+I{Y}_j)}^k={\tau }_j^k{H}_k\left({\tau }_j^{-1}{y}_j\right)\mathrm{where}I=\sqrt{-1},{\tau }_j={\left(V^{jj}\right)}^{1/2}:\hfill \\ & H^j=y_j,H^{j^2}=y_j^2-V^{jj},H^{j^3}=y_j^3-3V^{jj}{y}_j,H^{j^4}=y_j^4-6V^{jj}{y}_j^2+3{\left(V^{jj}\right)}^2.\hfill \end{array}$$

(33)

(34)

Takemura and Takeuchi (1988) showed how to extend (4) to the multivariate case by replacing (8) by

$$\begin{array}{cc}& X_n\approx X+\sum _{r=1}^{\infty }{n}^{-r/2}{g}_r\left(X\right).\hfill \end{array}$$

(35)

It would very useful to obtain these multivariate $g_r\left(X\right)$ explicitly. How does one extend the univariate formula for $g_r\left(x\right)$ in terms of the $h_r\left(x\right)$? Their §5 shows that the Cornish-Fisher expansions (4) are valid under the usual conditions for validity of the Edgeworth expansion (2).

NOTE 2.2.

Standard estimates have a natural extension to Type b estimates . These are estimates for which the cumulant expansion (5) is replaced by

$${\overline{k}}^{1-r}=k^{i_1\dots i_r}=\kappa ({\widehat{w}}^{i_1},\dots ,{\widehat{w}}^{i_r})\approx \sum _{j=2r-2}^{\infty }{n}^{-j/2}{\overline{b}}_j^{1-r}\mathrm{where}{\overline{b}}_j^{1-r}=b_j^{i_1\dots i_r}.$$

Examples are 1-sided confidence interval limits. These have the form

$$\widehat{w}\approx \sum _{j=0}^{\infty }{n}^{-j/2}{t}_j\left(\widehat{\theta }\right),$$

where $\widehat{\theta }$ is a standard estimate and $t_j$ are smooth functions. See Withers and Nadarajah (2010) for more details.

3. Secondary or Derived Expansions

Let V have Hermitian form $H^{\prime }\mathsf{\Lambda }H$ where $H^{\prime }H=I_q$.

$$\mathrm{Set}S=V^{-1/2}=H^{\prime }{\mathsf{\Lambda }}^{-1/2}H\mathrm{and}T=V^{1/2}=H^{\prime }{\mathsf{\Lambda }}^{1/2}H.$$

As in (18)–(20), we use tensor summation and the bar notation, and

$$\begin{array}{cc}& N\sim {\mathcal{N}}_q(0,I_q),X=TN\sim {\mathcal{N}}_q(0,V),Y=V^{-1}X=SN\sim {\mathcal{N}}_q(0,V^{-1}).\hfill \end{array}$$

From the 2nd Edgeworth expansion in (15), it follows that for $C\subset R^q$,

(1)

(2)

(3)

(4)

(5)

$$\begin{array}{cc}& p_{1C}={\tilde{p}}_{11}\left(C\right)+{\tilde{p}}_{13}\left(C\right),{\tilde{p}}_{r1}\left(C\right)={\overline{P}}_r^1{\overline{H}}^1\left(C\right),{\tilde{p}}_{r3}\left(C\right)={\overline{P}}_r^{1-3}{\overline{H}}^{1-3}\left(C\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}& p_{2C}=\sum _{k=2,4,6}{\tilde{p}}_{2k}\left(C\right)\mathrm{where}{\tilde{p}}_{r2}\left(C\right)={\overline{P}}_r^{12}{\overline{H}}^{12}\left(C\right),\hfill \\ & {\tilde{p}}_{r4}\left(C\right)={\overline{P}}_r^{1-4}{\overline{H}}^{1-4}\left(C\right),{\tilde{p}}_{r6}\left(C\right)={\overline{P}}_r^{1-6}{\overline{H}}^{1-6}\left(C\right).\hfill \end{array}$$

(7)

As ${\overline{H}}^{1-k}$ is a linear combination of ${\overline{y}}_1\dots {\overline{y}}_s$, we can write ${\overline{H}}^{1-k}\left(C\right)$ in terms of

$$\begin{array}{cc}& {\overline{Q}}^{1-s}=E{\overline{Y}}_1\dots {\overline{Y}}_{s}I(VY\in C)={\int }_C{\overline{y}}_1\dots {\overline{y}}_s{\varphi }_V\left(x\right)dx:\hfill \\ & {\overline{H}}^1\left(C\right)={\overline{Q}}^1,{\overline{H}}^{12}\left(C\right)={\overline{Q}}^{12}-p_{0C}{\overline{V}}^{12},{\overline{H}}^{1-3}\left(C\right)={\overline{Q}}^{1-3}-\sum ^3{\overline{Q}}^1{\overline{V}}^{23},\hfill \\ & {\overline{H}}^{1-4}\left(C\right)={\overline{Q}}^{1-4}-\sum ^6{\overline{Q}}^{12}{\overline{V}}^{34}+p_{0C}{\overline{\mu }}^{1-4},\hfill \\ & {\overline{H}}^{1-6}\left(C\right)={\overline{Q}}^{1-6}-\sum ^{15}{\overline{Q}}^{1-4}{\overline{V}}^{56}+\sum ^{15}{\overline{Q}}^{12}{\overline{\mu }}^{3-6}-p_{0C}{\overline{\mu }}^{1-6}.\hfill \end{array}$$

(8)

This gives $p_{1C}$ and $p_{2C}$ of (1). So it gives $Prob.(X_n\in C)$ to $O\left(n^{-3/2}\right)$. For $\widehat{w}$ a sample mean, ${\tilde{p}}_{22}\left(C\right)=0$ and for $k=4,6,{\overline{P}}_2^{1-k}$ needed for ${\tilde{p}}_{2k}\left(C\right)$, are given by Example 2.2. If $-C=C$, then for r odd, ${\overline{Q}}^{1-r}={\tilde{p}}_{rk}\left(C\right)=p_{rC}=0,$ so that

$$\begin{array}{cc}& Prob.(X_n\in C)\approx \sum _{r=0}^{\infty }{n}^{-r}{p}_{2rC}={\mathsf{\Phi }}_V\left(C\right)+n^{-1}{p}_{2C}+O\left(n^{-2}\right).\hfill \end{array}$$

(9)

We now consider three such C. Ellipsoidal C. Take $C=\{x:x^{\prime }{V}^{-1}x\le u\}$ for some $u>0$.

$$\mathrm{Then}-C=C,p_{0C}=Prob.(N^{\prime }N<u)=Prob.({\chi }_q^2<u).$$

(We might choose u such that $p_{0C}=.5$ or $.9$.) ${\overline{Q}}^{1-s}$ of (8) is given by

(10)

(11)

(12)

(13)

(14)

(15)

Similarly ${\overline{R}}^{1-k}=0$ unless k is even and $\{i_1,\dots ,i_k\}=\{1^{J_1},\dots ,q^{J_q}\}$ where $J_1,\dots ,J_q$ are even integers and $j^k$ is a string of js of length k.

For $Z_1,Z_2$ of (13), set

$$\begin{array}{cc}& R_{2k}=E{N}_1^{2k}I(N^{\prime }N<u)=E{Z}_1^{k}I(Z_1+Z_2<u).\hfill \\ & \mathrm{Set}{R}_{2k_1,2k_2}=E{N}_1^{2k_1}{N}_2^{2k_2}I(N^{\prime }N<u)=E{Z}_1^{k_1}{Z}_2^{k_2}I(Z_1+Z_2+Z_3<u),\hfill \end{array}$$

where now for $Z_1=N_1^2,Z_2=N_2^2$ and $Z_3=N^{\prime }N-Z_1-Z_2\sim {\chi }_{q-2}^2,$ so that $Z_1,Z_2,Z_3$ are independent. Set

$$\begin{array}{cc}& So,{\overline{S}}_2^{1-4}={\overline{\mu }}^{1-4}-3{\overline{S}}_1^{1-4},\hfill \end{array}$$

$$\begin{array}{cc}& {\overline{Q}}^{1-4}={\overline{S}}_1^{1-4}{r}_4+{\overline{\mu }}^{1-4}{R}_{22},\mathrm{where}{r}_4=R_4-3R_{22},\hfill \end{array}$$

(16)

$$\begin{array}{cc}& {\tilde{p}}_{r4}\left(C\right)={\overline{P}}_r^{1-4}[{\overline{S}}_1^{1-4}{r}_4+3{\overline{V}}^{12}{\overline{V}}^{34}{r}_5]\mathrm{where}{r}_5=R_{22}-2R_2+p_{0C}.\hfill \end{array}$$

(17)

Similarly, using ${\sum }_{j,k,l}^{\prime }$ to mean a sum over distinct $j,k,l$,

$$\begin{array}{cc}& {\overline{Q}}^{1-6}={\overline{S}}_{17}\dots {\overline{S}}_{6,12}[R_{6}I(i_7=\dots =i_{12})+R_{42}{I}_2+R_{222}{I}_3]\hfill \\ & \mathrm{where}{R}_{222}=E{N}_1^2{N}_2^2{N}_3^{2}I(N^{\prime }N<u),I_2=\sum ^{15}{\overline{S}}^{1-4:56},I_3=\sum ^{90}{\overline{S}}^{12:34:56},\hfill \\ & {\overline{S}}^{1-4:56}=\sum _{j\ne k}{S}_{i_{1}j}\dots S_{i_{4}j}{S}_{i_{5}k}{S}_{i_{6}k},{\overline{S}}^{12:34:56}=\sum _{j,k,l}^{\prime }{S}_{i_{1}j}{S}_{i_{2}j}{S}_{i_{3}k}{S}_{i_{4}k}{S}_{i_{5}l}{S}_{i_{6}l}.\hfill \\ & So,{\overline{P}}_r^{1-6}{\overline{Q}}^{1-6}={\overline{P}}_r^{1-6}[{\overline{S}}_1^{1-6}{R}_6+15{\overline{S}}^{1-4:56}{R}_{42}+90{\overline{S}}^{12:34:56}{R}_{222}],\hfill \\ & {\tilde{p}}_{r6}\left(C\right)={\overline{P}}_r^{1-6}[{\overline{Q}}^{1-6}-15{\overline{Q}}^{1-4}{\overline{V}}^{56}+15{\overline{Q}}^{12}{\overline{\mu }}^{3-6}-p_{0C}{\overline{\mu }}^{1-6}]\hfill \\ & ={\overline{P}}_r^{1-6}[{\overline{S}}_1^{1-6}{R}_6+15{\overline{S}}^{1-4:56}{R}_{42}+90{\overline{S}}^{12:34:56}{R}_{222}-15{\overline{S}}_1^{1-4}{\overline{V}}^{56}{r}_4\hfill \end{array}$$

$$\begin{array}{c}\hfill +3{\overline{V}}^{12}{\overline{V}}^{34}{\overline{V}}^{56}{r}_6]\mathrm{where}{r}_6=-15R_{42}+15R_2-p_{0C}.\end{array}$$

(18)

$p_{2C}$ of (9) is now given by (15), (17), (18). Note that

$$\begin{array}{cc}& {\overline{V}}^{12}{\overline{V}}^{34}{\overline{V}}^{56}={\overline{S}}_1^{1-6}+\sum ^3{\overline{S}}_1^{1-4:56}+{\overline{S}}_1^{12:34:56}.\hfill \end{array}$$

We can write $R_{2k_1,\dots }$ in terms of $F_k\left(u\right)=Prob.({\chi }_k^2\le u):$

$$\begin{array}{cc}& R_{2k}=E{Z}_1^k{F}_{q-1}(u-Z_1)={\int }_0^u{z}_1^k{F}_{q-1}(u-z_1)d{F}_1\left(z_1\right)=G_{q-1}(u:k)\mathrm{say},\hfill \\ & R_{2k_1,2k_2}=E{Z}_2^{k_2}{G}_{q-2}(u-Z_2:k_1)={\int }_0^u{z}_2^{k_2}{G}_{q-2}(u-z_2:k_1)d{F}_1\left(z_2\right)\hfill \\ & =G_{q-2}(u:k_1{k}_2)\mathrm{say},\hfill \\ & R_{222}=E{Z}_3{G}_{q-3}(u-Z_3:11)={\int }_0^u{z}_3{G}_{q-3}(u-z_3:11)d{F}_1\left(z_3\right).\hfill \end{array}$$

For $X,X_n$ of (6), set

$$Q=X^{\prime }{V}^{-1}X\sim {\chi }_q^2,Q_n=X_n^{\prime }{V}^{-1}{X}_n,{\widehat{Q}}_n=X_n^{\prime }{\widehat{V}}^{-1}{X}_n=n{(\widehat{w}-w)}^{\prime }{\widehat{V}}^{-1}(\widehat{w}-w),$$

where $\widehat{V}$ is an estimate $\widehat{V}$, (empirical, semi-parametric, or parametric depending on the model used). So

$$Prob.(Q_n\le u)=Prob.({\chi }_q^2\le u)+n^{-1}{p}_{2C}+O\left(n^{-2}\right)=1-\alpha +O\left(n^{-1}\right)\mathrm{if}u={\chi }_{q,1-\alpha }^2,$$

the $1-\alpha$ quantile of ${\chi }_q^2$ . It is often true that

$$Prob.({\widehat{Q}}_n\le u)=Prob.({\chi }_q^2\le u)+O\left(n^{-1}\right).$$

(An exact confidence region for w is only possible if $\widehat{w}$ is normal: see 5.3.2 of Anderson (1958).) ${\chi }_q^2/q=G_{\gamma }/\gamma$ where $\gamma =q/2$ and $G_{\gamma }$ is a gamma random variable with mean $\gamma$. So if $q=2$,

$$p_{0C}=Prob.(G_1<u/2)=1-e^{-u/2}=1-\alpha \mathrm{if}0<u=-2ln\alpha ,$$

$Q_n<-2ln\alpha$ with probability $1-\alpha +O\left(n^{-1}\right)$, and typically,

$${\widehat{Q}}_n\le -2ln\alpha \mathrm{wuth}\mathrm{probability}1-\alpha +O\left(n^{-1}\right).$$

By way of illustration, Figure 3.1 plots the elliptical contours $x=X_n$ when $Q_n=-2ln\alpha$, for $1-\alpha =.5,.9,.99,$ when

$$\begin{array}{cc}& V=\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right),V^{-1}=\left(\begin{array}{cc}2& -1\\ -1& 2\end{array}\right)/3.\hfill \end{array}$$

(19)

So the asymptotic correlation of ${\widehat{w}}_1$ and ${\widehat{w}}_2$ is 1/2.

Figure 3.1: $x=X_n$ when $Q_n=-2ln\alpha$ for $1-\alpha =.5,.9,.99$ courtesy of Dr Paul Teal:

One can do similarly for $q>2$, using ${\chi }_{q,1-\alpha }^2$ and 2-dimensional slices of these ellipsoids. For related references on confidence regions, see Withers and Nadarajah (2012a).

The R functions. If $q=1,$ then $R_{2k_1,2k_2}=R_{222}=0$ for $k_2>0$, and

$$R_{2k}=2g_{2k}\left(u^{1/2}\right)\mathrm{where}{g}_k\left(u\right)=E{N}^{k}I(0<N<u)={\int }_0^u{n}^k\varphi \left(n\right)dn$$

is given in Example 3.2 using a recurrence formula. (A simpler way is just to use (7).)

Now take $q=2$. Then

(20)

To get a recurrence formula for $g_{k_1{k}_2}=g_{k_1{k}_2}\left(v\right)$, set

$$\begin{array}{cc}& f_k\left(u\right)=g_k\left(u^{1/2}\right),A=\varphi \left(n_1\right)n_1^{k_1-1}{f}_{k_2}(v^2-n_1^2).\hfill \\ & \mathrm{Then}{g}_{k_1{k}_2}=-{\int }_0^{v}Ad\varphi \left(n_1\right)={\left[A\right]}_v^0+{\overline{g}}_{k_{1}2k_2}\mathrm{where}\hfill \\ & {\overline{g}}_{k_{1}2k_2}={\int }_0^v\varphi \left(n_1\right)dA=(k_1-1)A_1+A_2,\hfill \\ & A_1={\int }_0^{v}d{n}_1\varphi \left(n_1\right)n_1^{k_1-2}{f}_{k_2}(v^2-n_1^2)=g_{k_1-2,k_2},\hfill \\ & A_2={\int }_0^{v}d{n}_1\varphi \left(n_1\right)n_1^{k_1-1}(-2n_1){\stackrel{.}{f}}_{k_2}(v^2-n_1^2),{\stackrel{.}{f}}_k\left(u\right)=d{f}_k\left(u\right)/du.\hfill \\ & f_k\left(v^2\right)=g_k\left(v\right)\Rightarrow 2v{\stackrel{.}{f}}_{k_2}\left(v^2\right)={\stackrel{.}{g}}_k\left(v\right)=v^k\varphi \left(v\right).\hfill \\ & So{\stackrel{.}{f}}_{k_2}(v^2-n_1^2)={\left[x^{k-1}\varphi \left(x\right)\right]}_{x=v^2-n_1^2}={(v^2-n_1^2)}^{(k-1)/2}\varphi \left(v\right)e^{n_1^2/2}.\hfill \\ & So{A}_2=-\varphi \left(0\right)\varphi \left(v\right)A_3\mathrm{where}\mathrm{for}{B}_{12}=B(k_1+1)/2,k_2+1)/2),\hfill \\ & A_3={\int }_0^{v}d{n}_1{n}_1^{k_1}{v^2-n_1^2}^{(k_2-1)/2}=v^{k_1+k_2}{B}_{12}/2.\hfill \\ & So{A}_2=-\varphi \left(0\right)\varphi \left(v\right)v^{k_1+k_2}{B}_{12}/2=-A_{k_1{k}_2}\mathrm{say}.\hfill \\ & So\mathrm{for}{k}_1\ge 1,g_{k_1{k}_2}=(k_1-1)g_{k_1-2,k_2}+G_{k_1{k}_2}\hfill \\ & \mathrm{where}{G}_{k_1{k}_2}=I(k_1=1)g_{k_2}\left(v\right)+A_{k_1{k}_2}.\hfill \end{array}$$

(21)

This recurrence formula for $g_{k_1{k}_2}$ of (20) gives

$$\begin{array}{cc}& R_{22}=4g_{22}\left(u^{1/2}\right)\mathrm{where}{g}_{22}\left(v\right)=g_{22}=g_{02}+A_{22},\hfill \\ & \mathrm{and}{A}_{22}=v^4{e}^{-v^2}/32\mathrm{as}{B}_{12}=\mathsf{\pi }/8,\hfill \\ & R_{42}=4g_{42}\left(u^{1/2}\right)\mathrm{where}{g}_{42}\left(v\right)=g_{42}=g_{02}+A_{42},\hfill \\ & \mathrm{and}{A}_{42}=v^6{e}^{-v^2}/64\mathrm{as}{B}_{12}=\mathsf{\pi }/16.\hfill \end{array}$$

So we need $g_{02}$. By Example 3.2,

$$g_2\left(u\right)={\mathsf{\Phi }}_0\left(u\right)-u\varphi \left(u\right)\mathrm{where}{\mathsf{\Phi }}_0\left(u\right)=\mathsf{\Phi }\left(u\right)-1/2.$$

Transforming to $u^2=v^2-n_1^2,$

$$\begin{array}{cc}& g_{02}={\int }_0^{v}d{n}_1\varphi \left(n_1\right)g_2\left(u\right)=e^{v^2/2}(A-B),\hfill \\ & \mathrm{where}A={\int }_0^v{(v^2-u^2)}^{1/2}{\mathsf{\Phi }}_0\left(u\right)u\varphi \left(u\right)du,\hfill \\ & B=\varphi \left(0\right){\int }_0^v{u}^2{(v^2-u^2)}^{1/2}\varphi \left(u\right)du={\left(2\mathsf{\pi }\right)}^{-1}{v}^{2}I(-v^2),\hfill \\ & \mathrm{and}I\left(t\right)={\int }_0^1{x}^{1/2}{(1-x)}^{-1/2}{e}^{tx}dx=\sum _{i=0}^{\infty }B(i+3/2,1/2)t^i/i!.\hfill \\ & B(i+3/2,1/2)={\mathsf{\pi }}^{1/2}\mathsf{\Gamma }(i+3/2)/i!.\hfill \\ & \mathsf{\Gamma }(i+3/2)=\mathsf{\Gamma }(3/2){[3/2]}_i\mathrm{where}{\left[t\right]}_i=t(t+1)\dots (t+i-1).\hfill \\ & SoB=(v^2/8)1F_1(3/2,1:-v^2)\hfill \end{array}$$

where $1F_1(3/2,1:t)$ is the confluent or degenerate hypergeometric function: see 9.21 of Gradshteyn and Ryzhik (2000) and p504 of Abramowitz and Stegun (1964). Now suppose that $q>2$. Then

$$\begin{array}{cc}& R_{2k_1,2k_2}=4E{N}_1^{2k_1}{N}_2^{2k_2}I(N_1^2+N_2^2<u-Z,N_1>0,N_2>0)\hfill \\ & =4E{g}_{2k_1,2k_2}\left({(u-Z)}^{1/2}\right)\mathrm{where}Z=\sum _{i=3}^q\sim {\chi }_{q-2}^2,\hfill \\ & \mathrm{and}{R}_{222}=8E{N}_1^2{N}_2^2{N}_3^{2}I(N^{\prime }N<u,N_1>0,N_2>0,N_3>0)\hfill \\ & =8E{N}_3^2{R}_{22}(u-N_3^2-Z)I(N_3>0)\hfill \end{array}$$

where $R_{22}\left(u\right)=R_{22}$ of (20) and $Z={\sum }_{i=4}^q{N}_i^2\sim {\chi }_{q-3}^2$ is independent of $N_3$.

By Theorem 2.2 of Withers and Nadarajah (2012a), for $r\ge 1,$

$$\begin{array}{cc}& p_{rC}=\sum _{j=1}^{3r}{b}_{2r,j}{(U-1)}^j{F}_q\left(x\right)\mathrm{where}{F}_q\left(x\right)=Prob.({\chi }_q^2<x),U^j{F}_q\left(x\right)=F_{q+2j}\left(x\right),\hfill \\ \hfill & b_{2r,j}={\overline{b}}_{2r}^{1-2j}{\overline{\mu }}^{1-2j},{\overline{\mu }}^{1-2j}=E{\overline{X}}_1\dots {\overline{X}}_{2j}.\hfill \\ & \mathrm{For}\mathrm{example,},p_{rC}\mathrm{is}\mathrm{given}\mathrm{by}{b}_{21}=({\overline{k}}_2^{12}+{\overline{k}}_1^1{\overline{k}}_1^2){\overline{V}}^{12},\hfill \\ & b_{22}={\overline{k}}_3^{1-4}{\overline{V}}^{12}{\overline{V}}^{34}/4+{\overline{k}}_1^1{\overline{k}}_2^{34}{\overline{V}}^{12}{\overline{V}}^{34},\hfill \\ \hfill & b_{23}={\overline{k}}_2^{123}{\overline{k}}_2^{456}({\overline{V}}^{12}{\overline{V}}^{34}{\overline{V}}^{56}/4+{\overline{V}}^{14}{\overline{V}}^{25}{\overline{V}}^{36}/6).\hfill \end{array}$$

Its extension to ${\widehat{Q}}_n$ is given in §3 there for parametric and non-parametric models. Take $C=\{x:|{\left(V^{-1/2}x\right)}_j|\le u_j,j=1,\dots ,q\}$. So $-C=C$. For $Y=V^{-1/2}N$ and $N\sim {\mathcal{N}}_q(0,I_q)$,

$$p_{0C}={\mathsf{\Pi }}_{j=1}^q{G}_0\left(u_j\right)\mathrm{where}{G}_0\left(u\right)=Prob.\left(\right|N_1|<u)=2\mathsf{\Phi }\left(u\right)-1.$$

(We might choose $u_j\equiv u$ such that $p_{0C}=.5$ or $.9$.) For $S=V^{-1/2}$,

$$\begin{array}{cc}& {\overline{Q}}^{1-k}=E{\overline{Y}}_1\dots {\overline{Y}}_k{I\left(\right|N|}_j<u_j,j=1,\dots ,q)={\overline{S}}_{1,k+1}\dots {\overline{S}}_{k,2k}{\overline{R}}^{k+1-2k}\hfill \\ & \mathrm{where}\mathrm{now}{\overline{R}}^{1-k}=E{\overline{N}}_1\dots {\overline{N}}_k{I\left(\right|N|}_j<u_j,j=1,\dots ,q)={\overline{F}}^{1-k}(u_1\dots u_q)\mathrm{say}.\hfill \\ & \mathrm{Set}{R}_{2k_1,\dots ,2k_s}=E{N}_1^{2k_1}\dots N_s^{2k_s}{I\left(\right|N|}_j<u_j,j=1,\dots ,q)\hfill \\ & =\left[{\mathsf{\Pi }}_{j=1}^s{G}_{k_j}\left(u_j\right)\right]\left[{\mathsf{\Pi }}_{j=s+1}^q{G}_0\left(u_j\right)\right]\mathrm{where}\hfill \\ & G_k\left(u\right)=E{N}_1^{2k}I\left(\right|N_1|<u)=2g_{2k}\left(u\right)\mathrm{where}\hfill \\ & g_k\left(u\right)=g_k={\int }_0^u{n}^k\varphi \left(n\right)dn=-{\int }_0^u{n}^{k-1}d\varphi \left(n\right)={\left[n^{k-1}\varphi \left(n\right)\right]}_u^0+(k-1)g_{k-2}:\hfill \\ & g_0=\mathsf{\Phi }\left(u\right)-1/2,g_1=\varphi \left(0\right)-\varphi \left(u\right),g_2=g_0-u\varphi \left(u\right),\hfill \\ & g_3=2\varphi \left(0\right)-(u^2+2)\varphi \left(u\right),g_4=3g_0-(u^3+3u)\varphi \left(u\right),\hfill \\ & g_5=2.4\varphi \left(0\right)-(u^4+4u^2+4.2)\varphi \left(u\right),g_6=3.5g_0-(u^5+5u^3+5.3u)\varphi \left(u\right),\hfill \\ & g_7=2.4.6\varphi \left(0\right)-(u^6+6u^4+66.4u^2+6.4.2)\varphi \left(u\right),\hfill \\ & g_8=3.5.7g_0-(u^7+7u^5+7.5u^3+7.5.3u)\varphi \left(u\right).\hfill \\ & \mathrm{As}{R}^{12}=0,{\overline{Q}}^{12}={\overline{V}}^{12}{R}_2\mathrm{by}\left(14\right)\mathrm{where}{R}_{2k}=2^q{g}_{2k}\left(u_1\right){\mathsf{\Pi }}_{j=2}^q{g}_0\left(u_j\right).\hfill \\ & \mathrm{Set}{R}_{2k_1,2k_2}=2^q\left[{\mathsf{\Pi }}_{j=1}^2{g}_{2k_j}\left(u_j\right)\right]\left[{\mathsf{\Pi }}_{j=3}^q{g}_0\left(u_j\right)\right],\hfill \\ & R_{2k_1,2k_2,2k_3}=2^q\left[{\mathsf{\Pi }}_{j=1}^3{g}_{2k_j}\left(u_j\right)\right]\left[{\mathsf{\Pi }}_{j=4}^q{g}_0\left(u_j\right)\right],\hfill \end{array}$$

So now ${\overline{Q}}^{1-4}$ is given by (16) and (17) in terms of these new $R_{2k_1,\dots }$. Similarly ${\overline{Q}}^{1-6}$ and ${\tilde{p}}_{rk}\left(C\right)$ are given by Example 3.1 with these new $R_{2k_1,\dots }$. So now we have $p_{2C}$ of (9). If we choose $u_j\equiv u$, then

$$\begin{array}{cc}& R_{2k}=2^q{g}_{2k}{g}_0^{q-1},R_{2k_1,2k_2}=2^q\left[{\mathsf{\Pi }}_{j=1}^2{g}_{2k_j}\right]g_0^{q-2},\hfill \\ & R_{2k_1,2k_2,2k_3}=2^q\left[{\mathsf{\Pi }}_{j=1}^3{g}_{2k_j}\right]g_0^{q-3}.\hfill \end{array}$$

Take $C=\{x:|x_j|\le u_j,j=1,\dots ,q\}$. So $-C=C$. This choice gives a set of q simultaneous intervals. For $Y=SN$ and $T=V^{1/2}$

$$\begin{array}{cc}& p_{0C}=Q={Prob.\left(\right|\left(VY\right)}_j|<u_j{,j=1,\dots ,q)=Prob.(|\left(TN\right)}_j|<u_j,j=1,\dots ,q),\hfill \\ & {\overline{Q}}^{1-k}=E{\overline{Y}}_1\dots {\overline{Y}}_k{I\left(\right|\left(VY\right)}_j|<u_j,j=1,\dots ,q)={\overline{S}}_{1,k+1}\dots {\overline{S}}_{k,2k}{\overline{R}}^{k+1-2k}\hfill \\ & \mathrm{where}\mathrm{now}{\overline{R}}^{1-k}=E{\overline{N}}_1\dots {\overline{N}}_k{I\left(\right|\left(TN\right)}_j|<u_j,j=1,\dots ,q).\hfill \end{array}$$

So ${\overline{R}}^{1\dots k}=0$ for k odd. ($R^{12}$ is no longer 0.) But each of those non-zero require q numerical integrations. Consider the case $q=2$.

$$\begin{array}{cc}& I\left(\right|X_j|<u_j,j=1,2)=I(X_j<u_j,j=1,2)-I(X_1<-u_1,X_2<u_2)\hfill \\ & -I(X_1<u_1,X_2<-u_2)+I(X_1<-u_1,X_2<-u_2)\Rightarrow \hfill \\ & R^{jk}=F_{jk}(u_1,u_2)-F_{jk}(u_1,-u_2)-F_{jk}(-u_1,u_2)+F_{jk}(-u_1,-u_2)\hfill \\ & \mathrm{where}{F}_{jk}(u_1,u_2)=E{N}_j{N}_{k}I(X_1<u_1,X_2<u_2)\hfill \\ & =\int d{n}_1\int d{n}_2\varphi \left(n_1\right)\varphi \left(n_2\right)I(\sum _{k=1}^2{T}_{jk}{n}_k<u_j,j=1,2).\hfill \end{array}$$

For more on these types of expansions and their inverses, see Hill and Davis (1968) and their simplifications given in §4–§5 of Withers and Nadarajah (2015). However in these papers the need to express results in terms of $S=V^{-1/2}$ did not arise there. For the case of a sample mean with estimated covariance, it is possible to avoid the use of $S=V^{-1/2}$: see Xu and Gupta (2006).

4. The Distribution of $X_n=n^{1/2}(w-w)$ for $q=2$

As above, we set $y=V^{-1}x,Y=V^{-1}x$. We now switch notation to

(1)

(2)

So $H_{k0}$ and $H_{0k}$ are given by (33) with $j=1$ and 2: $H_{10}=y_1,H_{01}=y_2,$

$$\begin{array}{cc}& H_{11}=y_1{y}_2-{\mu }_{11},H_{21}=y_1^2{y}_2-{\mu }_{20}{y}_2-2{\mu }_{11}{y}_1,H_{12}=y_1{y}_2^2-{\mu }_{02}{y}_1-2{\mu }_{11}{y}_2,\hfill \end{array}$$

and so on. The ${\mu }_{ab}$ and $H_{ab}$ needed here are given in Appendix B. Let us write the cumulant expansion (5) as

$$\begin{array}{cc}& {\kappa }_{ab}({\widehat{w}}_1,{\widehat{w}}_2)\approx \sum _{d=a+b-1}^{\infty }{n}^{-d}{k}_{abd}\mathrm{for}a+b\ge 1\mathrm{where}{k}_{abd}=k_d^{1^a{2}^b}.\hfill \end{array}$$

So $k_{100}=w_1,k_{010}=w_2.$ The density of $X_n=n^{1/2}(\widehat{w}-w)$, $p_{X_n}\left(x\right)$, is given by (15) and (18) in terms of ${\tilde{p}}_r\left(x\right)$ and ${\tilde{p}}_{rk}$ of (18). As $bar{P}_r^{1-k}$ is symmetric,

$$\begin{array}{cc}& {\tilde{p}}_{rk}=\sum _{b=0}^k{P}_r(k-b,b)H_{k-b,b}\mathrm{where}{P}_r\left(ab\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{a+b}{a}}\right)P_r^{1^a{2}^b},\hfill \end{array}$$

(3)

for ${\overline{P}}_r^{1-k}$ of (12)–(14). So $P_r\left(ba\right)$ is just $P_r\left(ab\right)$ with 1 and 2 reversed.

$$\begin{array}{cc}& \mathrm{Set}\sum ^2{P}_r\left(ab\right)H_{ab}=P_r\left(ab\right)H_{ab}+P_r\left(ba\right)H_{ba}.\hfill \end{array}$$

Then ${\tilde{p}}_r\left(x\right)$ is given for $r=1,2,3$ by (27)–(29) in terms of

$$\begin{array}{cc}& {\tilde{p}}_{r1}=\sum ^2{P}_r\left(10\right)H_{10},{\tilde{p}}_{r3}=\sum ^2[P_r\left(30\right)H_{30}+P_r\left(21\right)H_{21}],r\mathrm{odd},\hfill \\ & {\tilde{p}}_{22}=\sum ^2{P}_2\left(20\right)H_{20}+P_2\left(11\right)H_{11},\hfill \\ & {\tilde{p}}_{24}=\sum ^2[P_2\left(40\right)H_{40}+P_2\left(31\right)H_{31}]+P_2\left(22\right)H_{22},\hfill \\ & {\tilde{p}}_{26}=\sum ^2[P_2\left(60\right)H_{60}+P_2\left(51\right)H_{51}+P_2\left(42\right)H_{42}]+P_2\left(33\right)H_{33},\hfill \\ & {\tilde{p}}_{35}=\sum ^2[P_3\left(50\right)H_{50}+P_3\left(41\right)H_{41}+P_3\left(32\right)H_{32}],\hfill \\ & {\tilde{p}}_{37}=\sum ^2[P_3\left(70\right)H_{70}+P_3\left(61\right)H_{61}+P_3\left(52\right)H_{52}+P_3\left(43\right)H_{43}],\hfill \\ & {\tilde{p}}_{39}=\sum ^2[P_3\left(90\right)H_{90}+P_3\left(81\right)H_{81}+P_3\left(72\right)H_{72}+P_3\left(63\right)H_{63}+P_3\left(54\right)H_{54}].\hfill \end{array}$$

The $P_r\left(ab\right)$ needed in (3) for ${\tilde{p}}_{rk},r\le 3,$ are as follows.

$$\begin{array}{cc}& \mathrm{For}{\tilde{p}}_{11},P_1\left(10\right)=k_1^1,so{P}_1\left(01\right)=k_1^2.\hfill \\ & \mathrm{For}{\tilde{p}}_{13},P_1\left(30\right)=k_2^{111}/6,so{P}_1\left(03\right)=k_2^{222}/6,P_1\left(21\right)=k_2^{112}/2.\hfill \end{array}$$

$$\begin{array}{cc}& \mathrm{For}{\tilde{p}}_{22},P_2\left(20\right)=k_2^{11}/2+{\left(k_1^1\right)}^2/2,P_2\left(11\right)=k_2^{12}+k_1^1{k}_1^2.\hfill \\ & \mathrm{For}{\tilde{p}}_{24},P_2\left(40\right)=k_3^{1^4}/24+k_1^1{k}_2^{1^3}/6,P_2\left(31\right)=k_3^{1112}/6+k_1^1{k}_2^{112}/3+k_1^2{k}_2^{111}/3,\hfill \\ & P_2\left(22\right)=k_3^{1122}/4+k_1^1{k}_2^{112}/2+k_1^2{k}_2^{122}/2,\hfill \\ & \mathrm{For}{\tilde{p}}_{26},P_2\left(60\right)={\left(k_2^{111}\right)}^2/72,P_2\left(51\right)=k_2^{111}{k}_2^{112}/12,\hfill \\ & P_2\left(42\right)=k_2^{111}{k}_2^{122}/12+{\left(k_2^{112}\right)}^2/8,\hfill \\ & P_2\left(33\right)=k_2^{111}{k}_2^{222}/36+k_2^{112}{k}_2^{122}/4,\hfill \\ & \mathrm{For}{\tilde{p}}_{31},P_3\left(10\right)=k_2^1.\mathrm{For}{\tilde{p}}_{33},P_3\left(30\right)=k_3^{111}/6+k_2^{11}{k}_1^1/2+{\left(k_1^1\right)}^3/6,\hfill \\ & P_3\left(21\right)=[k_3^{112}+2k_1^1{k}_2^{12}+k_1^2{k}_2^{11}+{\left(k_1^1\right)}^2{k}_1^2]/6.\hfill \\ & \mathrm{For}{\tilde{p}}_{35},P_3\left(50\right)=k_4^{1^5}/120+k_3^{1^4}{k}_1^1/24+k_2^{111}[k_2^{11}+{\left(k_1^1\right)}^2]/12,\hfill \\ & P_3\left(41\right)/5=P_3^{1^{4}2}=k_4^{1^{4}2}/120+S_1/24+S_2/12+S_3/12\mathrm{where}\hfill \\ & S_1=(4k_1^1{k}_3^{1112}+k_1^2{k}_3^{1111})/5,S_2=(3k_2^{11}{k}_2^{112}+2k_2^{12}{k}_2^{111})/5,\hfill \\ & S_3=[2k_1^1{k}_1^2{k}_2^{111}+3{\left(k_1^1\right)}^2{k}_2^{112}]/5,\hfill \\ & P_3\left(32\right)=10P_3^{1112}=k_4^{1112^2}/12+k_3^{1112}{k}_1^2/3+k_3^{1122}{k}_1^1/4+k_3^{111}{\left(k_1^2\right)}^2/12\hfill \\ & +k_3^{112}{k}_1^2{k}_1^2/2+k_3^{122}{\left(k_1^1\right)}^2/4.\hfill \\ & \mathrm{For}{\tilde{p}}_{37},P_3\left(70\right)=k_3^{1111}{k}_2^{111}/144+{\left(k_2^{1^3}\right)}^2{k}_1^1/72,\hfill \\ & P_3\left(61\right)=k_2^{112}{k}_3^{1111}/48+k_2^{111}{k}_3^{1112}/36+k_1^2{\left(k_2^{111}\right)}^2/72+k_1^1{k}_2^{111}{k}_2^{112}/12.\hfill \\ & P_3\left(52\right)=k_3^{1111}{k}_2^{122}/48+k_3^{1^{3}2}{k}_2^{112}/12+k_3^{1122}{k}_2^{111}/24+5k_2^{111}{k}_2^{112}{k}_1^2/24\hfill \\ & +k_2^{111}{k}_2^{122}{k}_1^1/12+{\left(k_2^{112}\right)}^2{k}_1^1/12,\hfill \\ & P_3\left(43\right)=35P_3^{1^4{2}^3}=A/144+B/72\mathrm{for}A=k_3^{1^4}{k}_2^{222}+13k_3^{1^{3}2}{k}_2^{122}+17k_3^{1122}{k}_2^{112}\hfill \\ & +k_3^{1222}{k}_2^{111},B=(k_2^{111}{k}_2^{222}+15k_2^{112}{k}_2^{122})k_1^1+(9k_2^{111}{k}_2^{122}+10k_2^{112}{k}_2^{112})k_1^2.\hfill \\ & \mathrm{For}{\tilde{p}}_{39},P_3\left(90\right)={(k_2^{1^3}/6)}^3,P_3\left(81\right)=9{\left(k_2^{111}\right)}^2{k}_2^{112}/6^3,\hfill \\ & \mathrm{and}\mathrm{for}{\overline{a}}^{123}={\overline{k}}_2^{123},P_3\left(72\right)=[9{\left(a^{111}\right)}^2{a}^{122}+3a^{111}{\left(a^{112}\right)}^2]/6^3,\hfill \\ & P_3\left(63\right)=3[{\left(a^{111}\right)}^2{a}^{222}+18a^{111}{a}^{112}{a}^{122}+9{\left(a^{112}\right)}^3]/6^3,\hfill \\ & P_3\left(54\right)=9[a^{111}{a}^{112}{a}^{222}+15a^{111}{\left(a^{122}\right)}^2+27{\left(a^{112}\right)}^2{a}^{122}]/6^3.\hfill \end{array}$$

(4)

(15) and (18) now give the density of $X_n=n^{1/2}(\widehat{w}-w)$ to $O\left(n^{-2}\right)$. Set

$$\begin{array}{cc}& H_{ab}^{*}={(-{\partial }_1)}^a{(-{\partial }_2)}^b{\mathsf{\Phi }}_V\left(x\right)={\int }_{-\infty }^x{H}_{ab}{\varphi }_V\left(x\right)dx.\hfill \\ & So{H}_{ab}^{*}=H_{a-1,b-1}{\varphi }_V\left(x\right)\mathrm{if}a\ge 2,b\ge 1,\hfill \end{array}$$

(5)

$$\begin{array}{cc}& H_{10}^{*}={\int }_{-\infty }^{x_2}{\varphi }_V\left(x\right)d{x}_2={\partial }_1{\mathsf{\Phi }}_V\left(x\right),H_{a0}^{*}={(-{\partial }_1)}^{a-1}{H}_{10}^{*}\mathrm{if}a\ge 2,\hfill \\ & H_{01}^{*}={\int }_{-\infty }^{x_1}{\varphi }_V\left(x\right)d{x}_1={\partial }_2{\mathsf{\Phi }}_V\left(x\right),H_{0b}^{*}={(-{\partial }_2)}^{b-1}{H}_{01}^{*}\mathrm{if}b\ge 1.\hfill \end{array}$$

The distribution of $X_n=n^{1/2}(\widehat{w}-w)$ is given by (15) and (19) in terms of $P_{rk}\left(x\right)$. $P_{rk}\left(x\right)$ is just ${\tilde{p}}_{rk}$ above with $H_{ab}$ replaced by $H_{ab}^{*}$ of (5). That is,

$$\begin{array}{cc}& P_{rk}\left(x\right)=\sum _{b=0}^k{P}_r(k-b,b)H_{k-b,b}^{*}.\hfill \end{array}$$

(6)

(15) and (19) now give $Prob.(X_n\le x)$ to $O\left(n^{-2}\right)$ for $X_n=n^{1/2}(\widehat{w}-w)$. For example for $x={(1,1)}^{\prime }$ and $x={(2,2)}^{\prime }$, the values of $H_{ab}$ are given in Example B.1. For $r=1,2,3$, $p_{rC}$ of (1) and (9) is given by replacing ${\tilde{p}}_{rk},H_{ab}$ above by ${\tilde{p}}_{rk}\left(C\right)$ of (5) and $H_{abC}=H^{1^a{2}^b}\left(C\right)$ of (Section 3). By §3,

$$\begin{array}{cc}& H_{10C}=Q_{10},H_{20C}=Q_{20}-p_{0C}{\mu }_{20},H_{11C}=Q_{11}-p_{0C}{\mu }_{11},\hfill \\ & H_{30C}=Q_{30}-3Q_{10}{\mu }_{20},H_{21C}=Q_{21}-2Q_{10}{\mu }_{11}-Q_{01}{\mu }_{20},\hfill \\ & H_{40C}=Q_{40}-6Q_{20}{\mu }_{20}+p_{0C}{\mu }_{40},H_{31C}=Q_{31}-3Q_{20}{\mu }_{11}-3Q_{11}{\mu }_{20}+p_{0C}{\mu }_{31},\hfill \\ & H_{22C}=Q_{22}-Q_{20}{\mu }_{02}-Q_{02}{\mu }_{20}-4Q_{11}{\mu }_{11}+p_{0C}{\mu }_{22},\hfill \\ & H_{60C}=Q_{60}-15Q_{40}{\mu }_{20}+15Q_{20}{\mu }_{40}-p_{0C}{\mu }_{60},\hfill \\ & H_{51C}=Q_{51}-10Q_{31}{\mu }_{20}+5Q_{11}{\mu }_{40}-5Q_{40}{\mu }_{11}+10Q_{20}{\mu }_{31}{p}_{0C}{\mu }_{51},\hfill \\ & H_{42C}=Q_{42}-6Q_{22}{\mu }_{20}+Q_{02}{\mu }_{40}-8Q_{31}{\mu }_{11}+8Q_{11}{\mu }_{31}-Q_{40}{\mu }_{02}+6Q_{20}{\mu }_{22}\hfill \\ & -p_{0C}{\mu }_{42},\hfill \\ & H_{33C}=Q_{33}+3\sum ^2(Q_{20}{\mu }_{13}-Q_{13}{\mu }_{20})-9Q_{22}{\mu }_{11}+9Q_{11}{\mu }_{22}-p_{0C}{\mu }_{33},\hfill \\ & \mathrm{where}{Q}_{ab}=E{Y}_1^a{Y}_2^{b}I(VY\in C).\hfill \end{array}$$

These expressions can be read off those for $H_{ab}$ in Appendix B. If $-C=C$, then $0=Q_{ab}=H_{abC}$ for $a+b$ odd. For Examples 3.1 and 3.2, $Q_{11}={\mu }_{11}{R}_2,$ and ${\tilde{p}}_{r2}\left(C\right)$ of (4) is given by (15) in terms of ${\overline{P}}_r^{12}{\overline{V}}^{12}={\sum }^2{P}_r\left(20\right){\mu }_{20}+2P_r\left(11\right){\mu }_{11}$. For $r=2$ this is given by (4). Let $\widehat{w}$ be a sample mean of a distribution with cumulants ${\kappa }_{ab}$. By Example 2.2, for $\left(rk\right)=\left(11\right),\left(22\right),\left(31\right),\left(33\right),$ ${\tilde{p}}_{rk}=0$, and we need the following Edgeworth coefficients.

$$\begin{array}{cc}& \mathrm{For}{\tilde{p}}_{13}\mathrm{and}{P}_{13}\left(x\right):P_1\left(30\right)={\kappa }_{30}/3!,P_1\left(21\right)={\kappa }_{21}/2,P_1\left(12\right)={\kappa }_{12}/2.\hfill \\ & \mathrm{For}{\tilde{p}}_{24},P_{24}\left(x\right):P_2\left(40\right)={\kappa }_{40}/4!,P_2\left(31\right)={\kappa }_{31}/6,P_2\left(22\right)={\kappa }_{22}/4.\hfill \\ & \mathrm{For}{\tilde{p}}_{26},P_{26}\left(x\right):P_2\left(60\right)={\kappa }_{30}^2/72,P_2\left(51\right)={\kappa }_{30}{\kappa }_{21}/12,\hfill \\ & P_2\left(42\right)=(2{\kappa }_{30}{\kappa }_{12}+3{\kappa }_{21}^2)/24,P_2\left(33\right)=({\kappa }_{30}{\kappa }_{03}+9{\kappa }_{21}{\kappa }_{12})/6.\hfill \\ & \mathrm{For}{\tilde{p}}_{35},P_{35}\left(x\right):P_3\left(50\right)={\kappa }_{50}/5!,P_3\left(41\right)/5={\kappa }_{41}/5!,P_3\left(32\right)={\kappa }_{32}/12.\hfill \\ & \mathrm{For}{\tilde{p}}_{37},P_{37}\left(x\right):P_3\left(70\right)={\kappa }_{40}{\kappa }_{30}/144,P_3\left(61\right)={\kappa }_{21}{\kappa }_{40}/48+{\kappa }_{30}{\kappa }_{31}/36,\hfill \\ & P_3\left(52\right)={\kappa }_{40}{\kappa }_{12}/48+{\kappa }_{31}{\kappa }_{21}/12+{\kappa }_{22}{\kappa }_{30}/24,\hfill \\ & P_3\left(43\right)=A/144\mathrm{for}A={\kappa }_{40}{\kappa }_{03}+13{\kappa }_{31}{\kappa }_{12}+17{\kappa }_{22}{\kappa }_{21}+{\kappa }_{13}{\kappa }_{30}.\hfill \\ & \mathrm{For}{\tilde{p}}_{39},P_{39}\left(x\right):P_3\left(90\right)={\kappa }_{30}^3/6^3,P_3\left(81\right)={\kappa }_{30}^2{\kappa }_{21}/6^3,\hfill \\ & P_3\left(72\right)=({\kappa }_{30}^2{\kappa }_{12}+3{\kappa }_{30}{\kappa }_{21}^2)/4.6^3,P_3\left(63\right)=({\kappa }_{30}^2{\kappa }_{03}+12{\kappa }_{30}{\kappa }_{21}{\kappa }_{12}\hfill \\ & +15{\kappa }_{21}^3)/28.6^3,P_3\left(54\right)=(5{\kappa }_{30}{\kappa }_{21}{\kappa }_{03}+9{\kappa }_{30}{\kappa }_{12}^2+21{\kappa }_{21}^2{\kappa }_{12})/35.6^3.\hfill \end{array}$$

For $a<b,P_r\left(ab\right)$ is $P_r\left(ba\right)$ with superscripts 1 and 2 reversed. The $P_2\left(ab\right)$ needed for $p_{2C}$ of (1) and (9) are those given above for ${\tilde{p}}_{24},{\tilde{p}}_{26}$. For a specific case of a bivariate sample mean, consider An entangled gamma model. Let $G_0,G_1,G_2$ be independent gamma random variables with means $\gamma ={\gamma }_0,{\gamma }_1,{\gamma }_2$. For $i=1,2$, set $X_i=G_0+G_i,w_i=E{X}_i=\gamma +{\gamma }_i$, and let $\widehat{w}$ be the mean of a random sample of size n distributed as $(X_1,X_2)$. So, $E\widehat{w}=w,$ and $n\widehat{w}\stackrel{\mathcal{L}}{=}{(G_{n0}+G_{n1},G_{n0}+G_{n2})}^{\prime }$ where $G_{n0},G_{n1},G_{n2}$ are independent gamma random variables with means $n\gamma ,n{\gamma }_1,n{\gamma }_2$. The rth order cumulants of $X={(X_1,X_2)}^{\prime }$ are ${\kappa }^{i^r}=(r-1)!w_i,$ and otherwise $(r-1)!\gamma$. For example ${\kappa }_{20}={\kappa }^{11}=w_1,{\kappa }_{02}={\kappa }^{22}=w_2,{\kappa }_{11}={\kappa }^{12}=\gamma ,{\kappa }_{21}={\kappa }^{112}={\kappa }_{12}={\kappa }^{122}=2!\gamma$, and

$$\begin{array}{cc}& V=\left(\begin{array}{cc}{w}_1& \gamma \\ \gamma & w_2\end{array}\right),V^{-1}=\left(\begin{array}{cc}{w}_2& -\gamma \\ -\gamma & w_1\end{array}\right)/D\mathrm{where}D=detV=w_1{w}_2-{\gamma }^2.\hfill \end{array}$$

So, $y=V^{-1}x={(w_2{x}_1-\gamma x_2,-\gamma x_1+w_1{x}_2)}^{\prime }/D.$ Set ${\nu }_i=w_i/\gamma .$ Then V has correlation ${\left({\nu }_1{\nu }_2\right)}^{-1/2}$. This ranges from 0 at $\gamma =0$ to 1 at $\gamma =\infty .$ So ${\widehat{w}}_1$ and ${\widehat{w}}_2$ are positively entangled. (For a negatively entangled example, replace $G_{n0}+G_{n2}$ by $-G_{n0}+G_{n2}$: the correlation is then $-{\left({\nu }_1{\nu }_2\right)}^{-1/2}$. An extension to $R^q$ is $n{\widehat{w}}_i={\lambda }_i{G}_{n0}+G_{ni},i=1,\dots q.$) For $c=0,1,\dots ,$ set

$$\sum ^2{w}_1^c{H}_{ab}=w_1^c{H}_{ab}+w_2^c{H}_{ba}.$$

By Example 3.1, $P_1\left(30\right)=w_1/3,$ (so $P_1\left(03\right)=w_2/3,)P_1\left(22\right)=2\gamma ,$

$$\begin{array}{cc}& {\tilde{p}}_1\left(x\right)={\tilde{p}}_{13}=\sum ^2[w_1{H}_{30}/3+\gamma H_{21}],\hfill \\ & {\tilde{p}}_2\left(x\right)={\tilde{p}}_{24}+{\tilde{p}}_{26}\mathrm{where}{\tilde{p}}_{24}=\sum ^2[w_1{H}_{40}/4+\gamma H_{31}]+\gamma H_{22},\hfill \\ & {\tilde{p}}_{26}=\sum ^2[w_1^2{H}_{60}/18+w_1\gamma (H_{51}+H_{42})/3]+(w_1{w}_2/9+{\gamma }^2)H_{33},\hfill \\ & {\tilde{p}}_3\left(x\right)=\sum _{k=5,7,9}{\tilde{p}}_{3k}\mathrm{where}5{\tilde{p}}_{35}=\sum ^2[w_1{H}_{50}+\gamma (H_{41}+H_{32})],\hfill \\ & 12{\tilde{p}}_{37}=\sum ^2[w_1^2{H}_{70}+7w_1\gamma H_{61}+6(w_1\gamma +2{\gamma }^2)H_{52}+(w_1{w}_2+w_1\gamma +30{\gamma }^2)H_{43}],\hfill \\ & 27{\tilde{p}}_{39}=\sum ^2[w_1^3{H}_{90}+9w_1^2\gamma H_{81}+(9w_1^2\gamma +3w_1{\gamma }^2)H_{72}\hfill \\ & +3(w_1^2{w}_2+18w_1{\gamma }^2+9{\gamma }^3)H_{63}+9(w_1{w}_2\gamma +15w_1{\gamma }^2+27{\gamma }^3)H_{54}].\hfill \end{array}$$

Now consider the entangled exponential case ${\gamma }_i\equiv 1$. So $w_i\equiv 2$, V and $V^{-1}$ are given by (19),

$$\begin{array}{cc}& {\tilde{p}}_1\left(x\right)={\tilde{p}}_{13}=\sum ^2[2H_{30}/3+H_{21}]\hfill \\ & {\tilde{p}}_{24}=\sum ^2[H_{40}/2+H_{31}]+H_{22},9{\tilde{p}}_{26}=2\sum ^2[H_{60}+3H_{51}+3H_{42}]+13H_{33},\hfill \\ & 5{\tilde{p}}_{35}=\sum ^2[2H_{50}+H_{41}+H_{32}],6{\tilde{p}}_{37}=\sum ^2[2H_{70}+7H_{61}+12H_{52}+18H_{43}],\hfill \end{array}$$

(7)

For $x={(1,1)}^{\prime }$ and $x={(2,2)}^{\prime }$, the values of $H_{ab}$ are given in Example B.1.

So if $x={(1,1)}^{\prime }$, then $y={(1,1)}^{\prime }/3,{\tilde{p}}_1\left(x\right)=-62/81\approx -.765,$ so that our measure of the accuracy of the CLT is $n^{-1/2}{\tilde{p}}_1\left(x\right)=-31/81\approx -.383$ for $n=4$ and $-62/243\approx -.255$

$$\begin{array}{cc}& \mathrm{For}x={(1,1)}^{\prime },{\tilde{p}}_{24}=4/81\approx .049,{\tilde{p}}_{26}=11552/3^8\approx 1.761,\hfill \\ & {\tilde{p}}_2\left(x\right)=11876/3^8\approx 1.810,{\tilde{p}}_{35}=352/3^5\approx 1.449,{\tilde{p}}_{37}=65139/3^8\approx 9.928,\hfill \\ & {\tilde{p}}_{39}=5530972/3^{12}\approx 10.408,{\tilde{p}}_3\left(x\right)=11577055/3^{12}\approx 21.784,\hfill \\ & p_{X_n}\left(x\right)/{\varphi }_V\left(x\right)\approx 1-.765n^{-1/2}+1.810n^{-1}+21.784n^{-3/2}+O\left(n^{-2}\right)\hfill \\ & \approx 1-.191+.113+.340\mathrm{if}n=16\text{so that only 3 terms can be used}.\hfill \end{array}$$

If $x={(2,2)}^{\prime }$, then $y={(2,2)}^{\prime }/3,{\tilde{p}}_1\left(x\right)=-64/81\approx -.790$, and our measure of the accuracy of the CLT is $n^{-1/2}{\tilde{p}}_1\left(x\right)=-32/81\approx -.395$ for $n=4$, and $-16/81\approx -.197$ for $n=16$.

$$\begin{array}{cc}& \mathrm{For}x={(2,2)}^{\prime },{\tilde{p}}_{24}=-58/81\approx -.716,{\tilde{p}}_{26}=22910/3^8\approx 3.492,{\tilde{p}}_2\left(x\right)\approx 2.776,\hfill \\ & {\tilde{p}}_{35}=392/3^5\approx 1.613,{\tilde{p}}_{37}=8000/3^7\approx 3.658,{\tilde{p}}_{39}=2528960/3^{12}\approx 4.759,\hfill \\ & {\tilde{p}}_3\left(x\right)=5330264/3^{12}\approx 10.030,\hfill \\ & p_{X_n}\left(x\right)/{\varphi }_V\left(x\right)\approx 1-.790n^{-1/2}+2.776n^{-1}+10.030n^{-3/2}+O\left(n^{-2}\right)\hfill \\ & \approx 1-.198+.174+.157\mathrm{if}n=16.\hfill \end{array}$$

As in Example 2.1, suppose that the distribution of $\widehat{w}$ is symmetric about w. By (31),

$$\begin{array}{cc}& 2{\tilde{p}}_{22}=k_2^{11}{H}_{20}+2k_2^{12}{H}_{11}+k_2^{22}{H}_{02},\hfill \\ & 24{\tilde{p}}_{24}=k_3^{1111}{H}_{40}+4k_3^{1112}{H}_{31}+6k_3^{1122}{H}_{22}+4k_3^{1222}{H}_{13}+k_3^{2222}{H}_{04},\hfill \end{array}$$

and $P_{2k}\left(x\right)$ needed for (32) is ${\tilde{p}}_{2k}$ with $H_{ab}$ replaced by $H_{ab}^{*}$ of (5). Now suppose that $\widehat{w}={\widehat{W}}_1-{\widehat{W}}_2$ where ${\widehat{W}}_1$ and ${\widehat{W}}_2$ are independent copies of a random vector $\widehat{W}.$ Then the cumulants of $\widehat{w}$ of odd order are zero, and the cumulants of even order are twice those of $\widehat{W}.$

Consider the case $\widehat{W}=\widehat{w}$, the bivariate entangled gamma of Example 5.2. So $w={(0,0)}^{\prime }$, the odd cumulants of $\widehat{w}$ are zero, and the odd ${\tilde{p}}_r\left(x\right),P_r\left(x\right)$ are zero. Denote $V,x,y,X,Y$ of Example 3.2 as $V_0,x_0,y_0,X_0,Y_0$. Then

$$\begin{array}{cc}& V=2V_0,X=2^{1/2}{X}_0,Y=2^{-1/2}{Y}_0,\hfill \\ & H_{ab}=2^{-(a+b)/2}{H}_{ab}(x_0,V_0)\mathrm{where}x=2^{1/2}{x}_0,y=2^{-1/2}{y}_0.\hfill \end{array}$$

By Example 2.2, ${\tilde{p}}_{22}=0,P_2\left(40\right)=k_3^{1^4}/24,P_2\left(31\right)={\kappa }_{40}/6,P_2\left(22\right)={\kappa }_{22}/4$, where ${\kappa }_{40}=12({\gamma }_1+{\gamma }_3),{\kappa }_{31}={\kappa }_{22}=12{\gamma }_3.$

$$\begin{array}{cc}& So,{\tilde{p}}_2\left(x\right)=\sum ^2[({\gamma }_1+{\gamma }_3)H_{40}+4{\gamma }_3{H}_{31}]+6{\gamma }_3{H}_{22}.\hfill \end{array}$$

Now consider the exponential case ${\gamma }_i\equiv 1$. So

$$\begin{array}{cc}& {\tilde{p}}_2\left(x\right)={\tilde{p}}_{24}=\sum ^2[H_{40}+2H_{31}]+2H_{22}\hfill \\ & =7/9\approx .778\mathrm{if}x=(1,1)\mathrm{or}-14/9\approx -1.556\mathrm{if}x=(2,2).\hfill \end{array}$$

So for $n=16,n^{-1}{\tilde{p}}_2\left(x\right)\approx .049\mathrm{if}x=(1,1),\mathrm{or}-.194\mathrm{if}x=(2,2).$ Here we used the values of $H_{ab}$ given in the example of Teal (2024). Takeuchi (1978) gave explicit expressions for the Cornish-Fisher expansion when $q=2$.

5. Conclusions

Most estimates of interest are standard estimates, including functions of sample moments, like the sample correlation, and any smooth multivariate function of Fisher’s k-statistics: see Stuart and Ord (1991) for these. In §2 we gave the density and distribution of $X_n=n^{1/2}(\widehat{w}-w)$ to $O\left(n^{-2}\right)$, for $\widehat{w}$ any standard estimate, in terms of functions of the cumulants coefficients ${\overline{k}}_j^{1-r}$ of (5), the coefficients ${\overline{P}}_r^{1-k}$ of (11)–(14). §3 gave explicit detail when $q=2$ using the dual notation $P_r\left(ab\right)$. §3 gave as examples Edgeworth expansions for the probability of $X_n$ lying in an ellipsoidal or hyperrectangular set.

6. Discussion

A good approximation for the distribution of an estimate, is vital for accurate inference. It enables one to explore the distribution’s dependence on underlying parameters. The analytic results given here avoid the need for simulation or jack-knife or bootstrap methods while providing greater accuracy than them. Hall (1992) used the Edgeworth expansion to show that the bootstrap gives accuracy to $O\left(n^{-1}\right)$. Hall (1988) said that “2nd order correctness usually cannot be bettered”. But this is not true using these analytic results. Simulation, while popular, can at best shine a light on behaviour when there is only a small number of parameters.

Estimates based on a sample of independent but not identically distributed random vectors, are also generally standard estimates. For example for a univariate sample mean $\overline{w}=n^{-1}{\sum }_{j=1}^n{X}_{jn}$ where $X_{jn}$ has rth cumulant ${\kappa }_{rjn}$, then ${\kappa }_r\left(\overline{w}\right)=n^{1-r}{\kappa }_r$ where ${\kappa }_r=n^{-1}{\sum }_{j=1}^n{\kappa }_{rjn}$ is the average rth cumulant. For some examples, see Skovgaard (1981a, 1981b) and Withers and Nadarajah (2010, 2020). The last is for a function of a weighted mean of complex random matrices. It gives an application in electrical engineering to channel capacity for multiple arrays.

For conditions for the validity of multivariate Edgeworth expansions, see Skovgaard (1986) and its references.

Cornish and Fisher (1937) and Fisher and Cornish (1960) did not deal with the question of when these expansions diverge. We showed how to confront this in the numerical examples in Example 4.2 and Withers (1984).

Lastly we discuss numerical computation. We have used Teal (2024) for the numerical calculations in Example 4.2. One could also download R-4.4.1 for Windows and google for the routines needed. dmvnorm computes the density function of the multivariate normal specified by mean and co-variance matrix. Use mvtnorm for the multivariate normal. See also qmvnorm for quantiles. rmvnorm generates multivariate normal variables. Googling ’bivariate hermite polynomials r’ gives

https://cran.r-project.org/web/packages/calculus/vignettes/hermite.html

One can then use install.packages("calculus"). However we have not used this route.

Some future directions.

1. Hill and Davis (1968) showed how to generalise the expansions of Cornish and Fisher about $N(0,1)$ to expansions about an arbitrary continuous distribution. Their results are cumbersome as they involve partition theory. In Withers and Nadarajah (2015) we overcame this using Bell polynomials. It would be straightforward to apply these to expansions about ${\chi }_q^2$ in Example 3.1, to obtain the percentiles of $n{(\widehat{w}-w)}^{\prime }{V}^{-1}(\widehat{w}-w)$ and $n{(\widehat{w}-w)}^{\prime }{\widehat{V}}^{-1}(\widehat{w}-w)$. However in the latter case we first need to derive the cumulant coefficients of ${\widehat{V}}^{-1/2}(\widehat{w}-w)$ from those of $\widehat{w}$. This can be done by applying Withers (2024).

2. It would very useful to obtain the multivariate $g_r\left(X\right)$ of (35) explicitly.

3. The multivariate expansions considered here have been about the multivariate normal. However as noted at the end of §1, expansions about other distributions can greatly reduce the number of terms in each $P_r\left(x\right),{\tilde{p}}_r\left(x\right)$ and $p_{rC}$ by matching bias and/or skewness. While this has been done for $q=1$ by Withers and Nadarajah (2011, 2012a, 2014) about Student’s distribution , the F-distribution, and the gamma, to date this has yet to be done for multivariate expansions about, for example, a multivariate gamma distribution.

4. The results given here are the first step for constructing confidence intervals and confidence regions. See Withers (1989).

5. The results here can be extended to tilted (saddlepoint) expansions by applying the results of Withers and Nadarajah (2010). These are very useful where convergence fails, that is, where the CLT cannot be improved upon, typically due to $\widehat{w}$ being in a tail. The tilted version of the multivariate distribution and density of a standard estimate are given by Corollaries 3, 4 there. Tilting was 1st used in statistics by Daniels (1954). He gave an approximation to the density of a sample mean. Jensen (1988) gave a univariate extension to $S_N$ where $S_n$ is the sum of i.i.d. observations and N is Poisson. The extension of the present results from ${\widehat{w}}_n=\widehat{w}$ to ${\widehat{w}}_N$ would be useful for both univariate and multivariate observations. For a review of references on tilting, see Withers and Nadarajah (2010).

6. Teal (2024) wrote a python program to obtain both analytic and numerical values of multivariate normal moments and multivariate Hermite polynomials when $q=2$. It would be useful to have these extended to $q=3$ and 4. (The dual notation for ${\overline{\mu }}^{1-k}$ and ${\overline{H}}^{1-k}$ when $q=3$ or 4 is straightforward.)