Testing for normality in any dimension based on a partial differential equation involving the moment generating function

TL;DR

This paper introduces a new class of affine invariant tests for multivariate normality based on a PDE involving the moment generating function, demonstrating strong power and consistency.

Contribution

It develops a novel testing method using PDEs for the moment generating function, applicable in any dimension, with proven null distribution and power against alternatives.

Findings

Tests are affine invariant and consistent.

Strong power against heavy-tailed distributions.

Connection with multivariate skewness measures.

Abstract

We use a system of first-order partial differential equations that characterize the moment generating function of the $d$-variate standard normal distribution to construct a class of affine invariant tests for normality in any dimension. We derive the limit null distribution of the resulting test statistics, and we prove consistency of the tests against general alternatives. In the case $d > 1$, a certain limit of these tests is connected with two measures of multivariate skewness. The new tests show strong power performance when compared to well-known competitors, especially against heavy-tailed distributions, and they are illustrated by means of a real data set.

Tables

$\gamma$	2.5	3	4	5	7
$𝔼\left[T_{\infty ,\gamma }\right]$	2.6013	0.7787	0.2022	0.0861	0.0277
$\text{Var}\left[T_{\infty ,\gamma }\right]$	4.7153	0.5430	0.0458	0.0094	0.0011

				$\gamma$
$n$	$2.5$	$3$	$4$	$5$	$7$	$10$	$\infty$
				$d=1$
20	$120.42$	$105.16$	$88.41$	$79.48$	$70.66$	$64.74$	$265.14$
50	$219.25$	$173.63$	$130.90$	$111.24$	$93.12$	$82.02$	$125.20$
100	$294.01$	$218.62$	$154.20$	$126.76$	$102.94$	$89.04$	$66.13$
200	$361.62$	$254.88$	$170.48$	$136.65$	$108.60$	$92.88$	$33.89$
300	$395.67$	$271.05$	$176.56$	$140.20$	$110.50$	$94.19$	$22.79$
				$d=2$
20	$535.50$	$413.65$	$300.96$	$249.60$	$203.09$	$174.73$	$628.97$
50	$1,086.28$	$737.69$	$464.16$	$356.56$	$268.57$	$220.27$	$291.96$
100	$1,516.50$	$947.61$	$546.91$	$402.84$	$292.63$	$235.24$	$152.02$
200	$1,867.44$	$1,089.76$	$585.61$	$419.94$	$299.04$	$238.36$	$77.02$
300	$2,035.48$	$1,141.98$	$595.36$	$422.25$	$299.14$	$238.42$	$51.52$
				$d=3$
20	$1,460.39$	$1,044.12$	$695.34$	$549.27$	$423.00$	$350.28$	$1,157.20$
50	$3,444.99$	$2,095.29$	$1,162.22$	$831.36$	$580.61$	$451.38$	$537.68$
100	$5,054.15$	$2,781.23$	$1,384.12$	$941.78$	$628.33$	$477.34$	$278.23$
200	$6,463.29$	$3,267.95$	$1,495.75$	$980.63$	$638.18$	$481.29$	$140.44$
300	$7,108.14$	$3,439.05$	$1,508.96$	$977.42$	$633.30$	$478.42$	$93.78$
				$d=5$
20	$6,346.44$	$4,065.35$	$2,389.07$	$1,759.71$	$1,257.17$	$986.77$	$2,903.55$
50	$20,164.36$	$10,268.49$	$4,655.52$	$2,988.80$	$1,862.85$	$1,340.01$	$1,361.04$
100	$34,187.51$	$15,193.42$	$5,934.77$	$3,545.86$	$2,070.71$	$1,439.25$	$705.30$
200	$47,436.90$	$18,844.25$	$6,578.37$	$3,746.81$	$2,114.14$	$1,450.44$	$355.91$
300	$54,128.44$	$20,321.98$	$6,715.41$	$3,749.27$	$2,091.10$	$1,436.11$	$237.36$

	$CvM$	$AD$	$SW$	$JB$	$Z_3$	$Z_{15}$	$T_{2.5}$	$T_5$	$T_{10}$	$T_{\infty }$
$\text{N}(0,1)$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$
NMIX1	$18$	$20$	$21$	$28$	$24$	$20$	$24$	$23$	$21$	$18$
NMIX2	$19$	$22$	$28$	$37$	$34$	$28$	$34$	$32$	$30$	$26$
$t\left(3\right)$	$58$	$61$	$63$	$69$	$65$	$57$	$65$	$63$	$60$	$52$
$t\left(5\right)$	$28$	$31$	$37$	$44$	$41$	$35$	$41$	$39$	$37$	$32$
$t\left(10\right)$	$12$	$13$	$16$	$21$	$20$	$17$	$20$	$20$	$19$	$17$
$LN(0,\frac{1}{2})$	$83$	$87$	$93$	$85$	$80$	$89$	$76$	$85$	$88$	$91$
$LN(0,\frac{1}{4})$	$31$	$35$	$44$	$39$	$37$	$45$	$34$	$40$	$44$	$47$
${\chi }^2\left(5\right)$	$74$	$81$	$89$	$75$	$69$	$82$	$62$	$74$	$80$	$83$
${\chi }^2\left(15\right)$	$30$	$34$	$43$	$36$	$34$	$43$	$31$	$38$	$41$	$45$
$Logistic(0,1)$	$14$	$16$	$19$	$26$	$24$	$20$	$24$	$23$	$21$	$19$
$Weibull\left(10\right)$	$25$	$28$	$34$	$30$	$28$	$36$	$25$	$31$	$35$	$37$
$Weibull\left(20\right)$	$39$	$44$	$53$	$46$	$44$	$53$	$40$	$48$	$52$	$55$
$P_{VII}\left(5\right)$	$27$	$30$	$36$	$43$	$40$	$35$	$41$	$39$	$37$	$32$
$P_{VII}\left(10\right)$	$10$	$12$	$16$	$21$	$20$	$17$	$20$	$19$	$18$	$16$
$SN\left(3\right)$	$31$	$34$	$40$	$32$	$30$	$39$	$25$	$33$	$37$	$41$
$SN\left(5\right)$	$53$	$59$	$67$	$49$	$43$	$58$	$36$	$49$	$55$	$61$

	$MS$	$MK$	$HZ$	$EN$	$HM$	$HJ$	$T_{2.5}$	$T_5$	$T_{10}$	$T_{\infty }$
$\text{N}(0,1)$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$
NMIX1	$85$	$34$	$75$	$82$	$57$	$73$	$48$	$69$	$80$	$86$
NMIX2	$44$	$48$	$29$	$38$	$57$	$53$	$55$	$54$	$52$	$44$
${\text{t}}_5(0,{\text{I}}_2)$	$53$	$62$	$42$	$51$	$67$	$60$	$60$	$60$	$58$	$53$
${\text{t}}_{10}(0,{\text{I}}_2)$	$24$	$26$	$14$	$19$	$32$	$29$	$29$	$29$	$28$	$25$
${\left({\chi }^2\left(15\right)\right)}^2$	$49$	$19$	$34$	$42$	$26$	$41$	$30$	$39$	$45$	$52$
${\left({\chi }^2\left(20\right)\right)}^2$	$40$	$16$	$27$	$33$	$24$	$34$	$25$	$32$	$37$	$42$
$\text{Logistic}{(0,1)}^2$	$24$	$27$	$15$	$19$	$33$	$28$	$28$	$29$	$28$	$25$
$\text{Gamma}{(5,1)}^2$	$67$	$27$	$52$	$61$	$38$	$57$	$41$	$54$	$62$	$70$
$\text{Gamma}{(4,2)}^2$	$76$	$32$	$64$	$72$	$42$	$66$	$48$	$62$	$71$	$78$
$P_{VII}{\left(10\right)}^2$	$20$	$21$	$11$	$14$	$27$	$23$	$24$	$24$	$23$	$20$
$P_{VII}{\left(20\right)}^2$	$11$	$10$	$7$	$8$	$14$	$12$	$13$	$13$	$12$	$12$
$\text{N}(0,1)\otimes$$t\left(3\right)$	$47$	$52$	$42$	$49$	$61$	$55$	$56$	$56$	$54$	$47$
$\text{N}(0,1)\otimes$${\chi }^2\left(5\right)$	$63$	$25$	$52$	$60$	$36$	$52$	$39$	$49$	$57$	$65$
$\text{N}(0,1)\otimes$${\chi }^2\left(10\right)$	$38$	$15$	$26$	$32$	$21$	$32$	$24$	$30$	$35$	$40$
${𝒮}^2\left(LN(0,\frac{1}{2})\right)$	$26$	$25$	$15$	$21$	$29$	$30$	$31$	$31$	$29$	$26$
${\text{NM}}_2\left(\rho =0.2\right)$	$6$	$6$	$5$	$6$	$6$	$6$	$6$	$6$	$6$	$6$

	$MS$	$MK$	$HZ$	$EN$	$HM$	$HJ$	$T_{2.5}$	$T_5$	$T_{10}$	$T_{\infty }$
$\text{N}{(0,1)}^3$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$
NMIX1	$89$	$36$	$81$	$91$	$59$	$72$	$43$	$66$	$82$	$91$
NMIX2	$71$	$76$	$49$	$66$	$79$	$79$	$79$	$80$	$78$	$72$
${\text{t}}_5(0,{\text{I}}_3)$	$68$	$78$	$55$	$68$	$77$	$73$	$71$	$73$	$73$	$69$
${\text{t}}_{10}(0,{\text{I}}_3)$	$34$	$38$	$18$	$27$	$35$	$38$	$36$	$38$	$38$	$34$
${\left({\chi }^2\left(15\right)\right)}^3$	$52$	$21$	$35$	$49$	$27$	$42$	$31$	$39$	$47$	$55$
${\left({\chi }^2\left(20\right)\right)}^3$	$40$	$16$	$26$	$37$	$21$	$33$	$24$	$30$	$36$	$44$
$\text{Logistic}{(0,1)}^3$	$28$	$30$	$15$	$22$	$33$	$31$	$30$	$31$	$31$	$28$
$\text{Gamma}{(5,1)}^3$	$72$	$30$	$53$	$69$	$39$	$58$	$41$	$53$	$65$	$75$
$\text{Gamma}{(4,2)}^3$	$80$	$36$	$65$	$79$	$46$	$66$	$47$	$61$	$73$	$83$
$P_{VII}{\left(10\right)}^3$	$22$	$22$	$10$	$16$	$24$	$25$	$25$	$26$	$25$	$23$
$P_{VII}{\left(20\right)}^3$	$12$	$10$	$6$	$8$	$14$	$13$	$13$	$13$	$13$	$12$
$\text{N}{(0,1)}^2\otimes$$t\left(3\right)$	$42$	$43$	$29$	$40$	$54$	$48$	$49$	$49$	$48$	$43$
$\text{N}{(0,1)}^2\otimes$${\chi }^2\left(5\right)$	$47$	$18$	$33$	$46$	$28$	$39$	$29$	$36$	$43$	$51$
$\text{N}{(0,1)}^2\otimes$${\chi }^2\left(10\right)$	$26$	$12$	$17$	$24$	$16$	$22$	$17$	$21$	$24$	$28$
${𝒮}^3\left(LN(0,\frac{1}{2})\right)$	$53$	$58$	$18$	$43$	$62$	$58$	$57$	$58$	$58$	$54$
${\text{NM}}_3\left(\rho =0.2\right)$	$8$	$7$	$5$	$6$	$7$	$8$	$8$	$8$	$8$	$8$

	$MS$	$MK$	$HZ$	$EN$	$HM$	$HJ$	$T_{2.5}$	$T_5$	$T_{10}$	$T_{\infty }$
$\text{N}{(0,1)}^5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$	$5$
NMIX1	$82$	$33$	$74$	$94$	$43$	$58$	$34$	$51$	$68$	$86$
NMIX2	$94$	$94$	$68$	$89$	$95$	$95$	$95$	$96$	$95$	$94$
${\text{t}}_5(0,{\text{I}}_5)$	$88$	$94$	$72$	$88$	$89$	$90$	$86$	$89$	$90$	$89$
${\text{t}}_{10}(0,{\text{I}}_5)$	$54$	$58$	$23$	$45$	$51$	$55$	$51$	$55$	$57$	$55$
${\chi }^2{\left(15\right)}^5$	$51$	$22$	$30$	$52$	$26$	$39$	$29$	$36$	$44$	$56$
${\chi }^2{\left(20\right)}^5$	$39$	$16$	$22$	$39$	$20$	$30$	$23$	$28$	$33$	$42$
$\text{Logistic}{(0,1)}^5$	$33$	$34$	$13$	$25$	$31$	$34$	$31$	$34$	$36$	$33$
$\text{Gamma}{(5,1)}^5$	$72$	$33$	$49$	$74$	$37$	$55$	$40$	$51$	$63$	$76$
$\text{Gamma}{(4,2)}^5$	$81$	$40$	$60$	$84$	$40$	$64$	$47$	$59$	$72$	$85$
$P_{VII}{\left(10\right)}^5$	$27$	$25$	$9$	$19$	$26$	$28$	$26$	$28$	$29$	$27$
$P_{VII}{\left(20\right)}^5$	$12$	$9$	$6$	$9$	$12$	$12$	$11$	$12$	$12$	$12$
$\text{N}{(0,1)}^4\otimes$$t\left(3\right)$	$35$	$32$	$16$	$30$	$42$	$39$	$38$	$39$	$39$	$35$
$\text{N}{(0,1)}^4\otimes$${\chi }^2\left(5\right)$	$28$	$13$	$16$	$28$	$19$	$23$	$19$	$22$	$25$	$31$
$\text{N}{(0,1)}^4\otimes$${\chi }^2\left(10\right)$	$16$	$8$	$10$	$15$	$13$	$14$	$12$	$13$	$15$	$18$
${𝒮}^5\left(LN(0,\frac{1}{2})\right)$	$89$	$95$	$77$	$90$	$90$	$90$	$86$	$90$	$91$	$89$
${\text{NM}}_5\left(\rho =0.2\right)$	$12$	$9$	$5$	$11$	$12$	$13$	$12$	$13$	$13$	$12$

$\gamma =2.5$	$\gamma =3$	$\gamma =4$	$\gamma =5$	$\gamma =7$	$\gamma =10$	$\gamma =\infty$
$0.0072$	$0.0061$	$0.0044$	$0.0033$	$0.00210$	$0.0013$	$0.0002$

Equations

$$H_0:\mathrm{The}\mathrm{law}\mathrm{of}X\mathrm{belongs}\mathrm{to}{\mathcal{N}}_d,$$

$$Y_{n,i}^{\top }{Y}_{n,j}=(X_i-{\overline{X}}_n{)}^{\top }{S}_n^{-1}(X_j-{\overline{X}}_n),i,j\in \{1,\dots ,n\},$$

$$\frac{\partial m(t)}{\partial t_j}=t_{j}m(t),t=(t_1,\dots ,t_d{)}^{\top }\in {\mathbb{R}}^d,j=1,\dots ,d.$$

$$m(t)=c_2(t_2,\dots ,t_d)\cdot {\text{e}}^{t_1^2/2}$$

$$\frac{\partial }{\partial t_2}{c}_2(t_2,\dots ,t_d)={\text{e}}^{-t_1^2/2}\cdot \frac{\partial }{\partial t_2}m(t).$$

$$t_2=\frac{\frac{\partial }{\partial t_2}{c}_2(t_2,\dots ,t_d)}{c_2(t_2,\dots ,t_d)},t\in {\mathbb{R}}^d,$$

$$m(t)=\prod _{j=1}^d{\text{e}}^{t_j^2/2}={\text{e}}^{\Vert t{\Vert }^2/2},t\in {\mathbb{R}}^d,$$

$$M_n(t)=\frac{1}{n}\sum _{j=1}^n{\text{e}}^{t^{\top }{Y}_{n,j}}$$

$$T_{n,\gamma }:=n{\int }_{{\mathbb{R}}^d}\Vert M_n^{\prime }(t)-t{M}_n(t){\Vert }^2{w}_{\gamma }(t)\text{d}t,$$

$$w_{\gamma }(t)=\exp \left(-\gamma \Vert t{\Vert }^2\right),t\in {\mathbb{R}}^d.$$

$${\int }_{{\mathbb{R}}^d}{\text{e}}^{t^{\top }\alpha }{w}_{\gamma }(t)\text{d}t$$

$${\int }_{{\mathbb{R}}^d}{\text{e}}^{t^{\top }\alpha }{t}^{\top }\alpha w_{\gamma }(t)\text{d}t$$

$${\int }_{{\mathbb{R}}^d}{\text{e}}^{t^{\top }\alpha }\Vert t{\Vert }^2{w}_{\gamma }(t)\text{d}t$$

$$T_{n,\gamma }=\frac{1}{n}{\left(\frac{\pi }{\gamma }\right)}^{d/2}\sum _{j,k=1}^n\exp \left(\frac{\Vert Y_{n,j,k}^{+}{\Vert }^2}{4\gamma }\right)\left(Y_{n,j}^{\top }{Y}_{n,k}-\frac{\Vert Y_{n,j,k}^{+}{\Vert }^2}{2\gamma }+\frac{d}{2\gamma }+\frac{\Vert Y_{n,j,k}^{+}{\Vert }^2}{4{\gamma }^2}\right),$$

$$W_n(t):=\sqrt{n}\left(M_n^{\prime }(t)-t{M}_n(t)\right)=\frac{1}{\sqrt{n}}\sum _{j=1}^n{\text{e}}^{t^{\top }{Y}_{n,j}}\left(Y_{n,j}-t\right),t\in {\mathbb{R}}^d,$$

$$⟨f,g⟩={\int }_{{\mathbb{R}}^d}f(t)g(t)w_{\gamma }(t)\text{d}t,\Vert f{\Vert }_{{\text{L}}^2}=\sqrt{⟨f,f⟩},$$

$$⟨f,g{⟩}_{\oplus }:=⟨f_1,g_1{⟩}_{{\text{L}}^2}+\dots +⟨f_d,g_d{⟩}_{{\text{L}}^2},$$

$$\Vert f{\Vert }_{\oplus }^2=\sum _{i=1}^d\Vert f_i{\Vert }_{{\text{L}}^2}^2={\int }_{{\mathbb{R}}^d}\Vert f(t){\Vert }^2{w}_{\gamma }(t)\text{d}t,f=(f_1,\dots ,f_d)\in \mathcal{H}.$$

$$K(s,t)={\mathrm{e}}^{(\Vert s{\Vert }^2+\Vert t{\Vert }^2)/2}\left({\mathrm{e}}^{s^{\top }t}\left(t{s}^{\top }+{\mathrm{I}}_{\mathrm{d}}\right)-t{s}^{\top }-(1+s^{\top }t){\mathrm{I}}_d\right),$$

$$\Vert W{\Vert }_{\oplus }^2={\int }_{{\mathbb{R}}^d}\Vert W(t){\Vert }^2{w}_{\gamma }(t)\mathrm{d}t.$$

$$Z_j(t):={\text{e}}^{t^{\top }{Y}_{n,j}}\left(Y_{n,j}-t\right),Z_j^{\circ }(t):={\text{e}}^{t^{\top }{X}_j}\left(X_j-t\right),t\in {\mathbb{R}}^d,$$

$$W_n^{\circ }(t):=\frac{1}{\sqrt{n}}\sum _{j=1}^n{Z}_j^{\circ }(t).$$

$${\Delta }_{n,j}=Y_{n,j}-X_j=\left(S_n^{-1/2}-{\text{I}}_d\right)X_j-S_n^{-1/2}{\overline{X}}_n,j=1,\dots ,n.$$

$${\text{e}}^{t^{\top }{\Delta }_{n,j}}=1+t^{\top }{\Delta }_{n,j}+\frac{1}{2}{\left(t^{\top }{\Delta }_{n,j}\right)}^2\exp \left({\Theta }_{n,j}{t}^{\top }{\Delta }_{n,j}\right),$$

$$W_n(t)=\frac{1}{\sqrt{n}}\sum _{j=1}^n{\text{e}}^{t^{\top }{X}_j}(X_j-t)+A_n(t)+B_n(t)+C_n(t),$$

$$A_n(t)$$

$$B_n(t)$$

$$C_n(t)$$

$$B_n(t)=-\frac{1}{\sqrt{n}}\exp \left(\frac{\Vert t{\Vert }^2}{2}\right)\sum _{j=1}^n\left(\frac{X_j{X}_j^{\top }t-t}{2}+X_j\right)+o_{\mathbb{P}}(1).$$

$${\tilde{Z}}_j(t)={\text{e}}^{t^{\top }{X}_j}(X_j-t)-{\text{e}}^{\Vert t{\Vert }^2/2}\left(X_j{X}_j^{\top }t-t+X_j\right).$$

$$\mathbb{E}\left[T_{\infty ,\gamma }\right]={\left(\frac{\pi }{\gamma -2}\right)}^{d/2}\left(d+\frac{d}{2(\gamma -2)}\right)-\left(\frac{d(d+1)}{2(\gamma -1)}+d\right){\left(\frac{\pi }{\gamma -1}\right)}^{d/2}.$$

$$\mathbb{E}\left[T_{\infty ,\gamma }\right]={\int }_{{\mathbb{R}}^d}\mathbb{E}\Vert W(t){\Vert }^2{w}_{\gamma }(t)\text{d}t.$$

Full text

Testing for normality in any dimension based on a partial differential equation involving the moment generating function

Norbert Henze1, Jaco Visagie2

*1**Institute of Stochastics, Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany

henze@kit.edu*

*2**Department of Statistics, University of Pretoria, South Africa

jaco.visagie@up.ac.za

**Abstract. We use a system of first-order partial differential equations that characterize the moment generating function of the $d$-variate standard normal distribution to construct a class of affine invariant tests for normality in any dimension. We derive the limit null distribution of the resulting test statistics, and we prove consistency of the tests against general alternatives. In the case $d>1$, a certain limit of these tests is connected with two measures of multivariate skewness. The new tests show strong power performance when compared to well-known competitors, especially against heavy-tailed distributions, and they are illustrated by means of a real data set.

Keywords. Moment generating function; test for multivariate normality; direct sum of Hilbert spaces; multivariate skewness; weighted $L^{2}$-statistic

2010 AMS classification: 62H15, 62G20.**

1 Introduction

It is often of interest to check whether an observed $d$-dimensional dataset is compatible with the assumption of coming from a multivariate normal distribution. Such a model check is of practical interest to researchers, as many multivariate techniques rely on the assumption of multivariate normality (for short: multinormality). As a consequence, it is not surprising that there is ongoing interest in testing for multinormality, as can be witnessed by various recent papers on the subject, see for example [28], [33], [10] and [23] as well as the references therein. Research into the practical implementation of these tests is also undertaken regularly, see, e.g., [25] and [23] regarding the implementation in the statistical software package R.

The purpose of this paper is not to review the multitude of tests that hitherto has been proposed (for an account of the state of the art regarding affine invariant procedures before 2002, see [16]), but to introduce and study a new class of tests that is based on a certain partial differential equation. To be specific, let $X_{1},X_{2},\ldots,X_{n},\ldots$ be independent and identically distributed (i.i.d.) copies of a $d$-dimensional random (column) vector $X$, the distribution of which is assumed to be absolutely continuous with respect to $d$-dimensional Lebesgue measure. All random vectors are defined on a common probability space $(\Omega,{\cal A},\mathbb{P})$.

Writing N${}_{d}(\mu,\Sigma)$ for the $d$-dimensional normal distribution with mean vector $\mu$ and non-degenerate covariance matrix $\Sigma$ and ${\cal N}_{d}$ for the class of all non-degenerate $d$-variate normal distributions, a classical problem is to test the null hypothesis

[TABLE]

against general alternatives. Since the class ${\cal N}_{d}$ is closed with respect to full rank affine transformations, any genuine test statistic $T_{n}=T_{n}(X_{1},\ldots,X_{n})$ based on $X_{1},\ldots,X_{n}$ should also be invariant with respect to such transformations, i.e., we should have $T_{n}(AX_{1}+b,\ldots,AX_{n}+b)=T_{n}(X_{1},\ldots,X_{n})$ for each nonsingular ($d\times d$)-matrix $A$ and each $b\in\mathbb{R}^{d}$, see [16] for an account of the importance of affine invariance in connection with testing for multivariate normality. Writing $\overline{X}_{n}=n^{-1}\sum_{j=1}^{n}X_{j}$ for the sample mean and $S_{n}=n^{-1}\sum_{j=1}^{n}(X_{j}-\overline{X}_{n})(X_{j}-\overline{X}_{n})^{\top}$ for the sample covariance matrix of $X_{1},\ldots,X_{n}$, where the superscript $\top$ denotes the transpose of vectors and matrices, a necessary and sufficient condition for a test statistic $T_{n}$ to be affine invariant is that it is based on the scalar products

[TABLE]

where $Y_{n,j}=S_{n}^{-1/2}(X_{j}-\overline{X}_{n})$, $j=1,\ldots,n$, are the so-called scaled residuals of $X_{1},\ldots,X_{n}$, see [16]. Here, $S_{n}^{-1/2}$ denotes the unique symmetric square root of $S_{n}^{-1}$ which, due to the absolute continuity of the distribution of $X$, exists with probability one if $n\geq d+1$, see [9]. The latter condition is tacitly assumed to hold in what follows.

The novel idea for constructing a test of $H_{0}$ is the following: Suppose $X$ is a $d$-dimensional random vector, and assume that the moment generating function (MGF) $m(t)=\mathbb{E}[\exp(t^{\top}X)]$ exists for each $t\in\mathbb{R}^{d}$ and satisfies the system of partial differential equations

[TABLE]

Writing $y^{\prime}(t)$ for the gradient of a function $y:\mathbb{R}^{d}\to\mathbb{R}$ at $t$, (1.2) may be succinctly written as $m^{\prime}(t)=tm(t)$, $t\in\mathbb{R}^{d}$. Taking into account that $m(0)=1$, the unique solution of this equation in the case $d=1$ is $m(t)=\exp(t^{2}/2)$, which is the MGF of the standard normal distribution. If $d>1$, and we fix $t_{2},\ldots,t_{d}$, the solution of (1.2) for $j=1$ is

[TABLE]

for some function $c_{2}:\mathbb{R}^{d-1}\to\mathbb{R}$. We thus have $c_{2}(t_{2},\ldots,t_{d})=\textrm{e}^{-t_{1}^{2}/2}\cdot m(t)$, $t\in\mathbb{R}^{d}$, which shows that $c_{2}$ is differentiable. Moreover,

[TABLE]

Using (1.2) with $j=2$ then gives

[TABLE]

the solution of which is $c_{2}(t_{2},\ldots,t_{d})=c_{3}(t_{3},\ldots,t_{d})\cdot\exp(t_{2}^{2}/2)$ for some function $c_{3}:\mathbb{R}^{d-2}\to\mathbb{R}$. Inserting this expression into (1.3) and continuing this way, we finally obtain

[TABLE]

where $\|\cdot\|$ denotes the Euclidean norm in $\mathbb{R}^{d}$. Notice that this unique solution of (1.2) is the MGF of the standard normal distribution N${}_{d}(0,\textrm{I}_{d})$, where $\textrm{I}_{d}$ is the unit matrix of order $d$.

If $X$ has some non-degenerate normal distribution, the scaled residuals $Y_{n,1},\ldots,Y_{n,n}$ should be approximately independent, with a distribution close to N${}_{d}(0,\textrm{I}_{d})$, at least for large $n$. Hence, a natural approach for testing $H_{0}$ is to consider the empirical MGF

[TABLE]

of $Y_{n,1},\ldots,Y_{n,n}$, and to employ the weighted $L^{2}$-statistic

[TABLE]

where

[TABLE]

Rejection of $H_{0}$ is for large values of $T_{n,\gamma}$. The role of $\gamma>0$ will be discussed later.

Using the relations

[TABLE]

and putting $Y_{n,j,k}^{+}=Y_{n,j}+Y_{n,k}$, the test statistic defined in (1.5) takes the form

[TABLE]

which is amenable to computational purposes. Notice that $T_{n,\gamma}$ is affine invariant.

The remainder of this paper is structured as follows: Section 2 deals with the convergence in distribution of $T_{n,\gamma}$ under $H_{0}$, and Section 3 is devoted to the problem of consistency of the new test (which seems to be new even in the univariate case). In Section 4 we let $\gamma$ tend to infinity while keeping $n$ fixed. Under this setting, $T_{n,\gamma}$ converges to a certain linear combination of two well-known measures of multivariate skewness. In Section 5 we compare the finite-sample power behavior of the new test to that of several classical and recent tests for both univariate and multivariate normality. Section 6 illustrates the procedures by means of a real data set. Section 7 presents some conclusions, while Section 8 contains several technical proofs which do not form part of the main text.

2 The limit null distribution of $T_{n,\gamma}$

In this section, we derive the limit distribution of the test statistic $T_{n,\gamma}$ defined in (1.5) under the null hypothesis (1.1). In view of affine invariance, we will assume $\mathbb{E}\left(X\right)=0$ and $\mathbb{E}\left(XX^{\top}\right)=\textrm{I}_{d}$. Put

[TABLE]

and let $\textrm{L}^{2}:=\textrm{L}^{2}(\mathbb{R}^{d},{\cal B}^{d},w_{\gamma}(t)\textrm{d}t)$ denote the separable Hilbert space of (equivalence classes of) measurable functions $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that are square integrable with respect to the finite measure on the $\sigma$-field ${\cal B}^{d}$ of Borel sets of $\mathbb{R}^{d}$, given by the weight function $w_{\gamma}$. The inner product and the resulting norm on $\textrm{L}^{2}$ will be denoted by

[TABLE]

respectively. Since $W_{n}(t)$ in (2.1) is $\mathbb{R}^{d}$-valued, we consider the Hilbert space ${\cal H}$, which is the $d$-fold (orthogonal) direct sum ${\cal H}:={\textrm{L}}^{2}\oplus\cdots\oplus{\textrm{L}}^{2}$, consisting of all ordered $d$-tuples $f=(f_{1},\ldots,f_{d})\in{\textrm{L}}^{2}\times\cdots\times{\textrm{L}}^{2}$, equipped with the inner product

[TABLE]

where $f=(f_{1},\ldots,f_{d})$, $g=(g_{1},\ldots,g_{d})\in{\cal H}$. Notice that the norm $\|\cdot\|_{\oplus}$ on ${\cal H}$ satisfies

[TABLE]

Since $W_{n}$ is a random element of ${\cal H}$ and $T_{n}=\|W_{n}\|^{2}_{\oplus}$, the aim is to prove $W_{n}\stackrel{{\scriptstyle{\cal D}}}{{\longrightarrow}}W$ for some centred Gaussian random element of ${\cal H}$. By the Continuous Mapping Theorem (CMT), we would then have $T_{n}\stackrel{{\scriptstyle{\cal D}}}{{\longrightarrow}}\|W\|^{2}_{\oplus}$. Here and in the sequel, $\stackrel{{\scriptstyle{\cal D}}}{{\longrightarrow}}$ denotes convergence in distribution of random elements (especially: of random variables). Moreover, $o_{\mathbb{P}}(1)$ refers to convergence in probability to zero of random elements ${\cal H}$. The main result of this section is as follows.

Theorem 2.1

(Convergence of $W_{n}$)*

Suppose that $X$ has some non-degenerate $d$-variate normal distribution. If the weight function $w_{\gamma}$ in (1.6) satisfies $\gamma>2$, there is some centred Gaussian random element $W$ of ${\cal H}$ having covariance (matrix) kernel K(s,t)=\mathbb{E}\big{[}W(s)W(t)^{\top}\big{]}, where*

[TABLE]

$s,t\in\mathbb{R}^{d}$, so that $W_{n}\stackrel{{\scriptstyle{\cal D}}}{{\longrightarrow}}W$ as $n\to\infty$.

Corollary 2.2

The limit distribution of $T_{n,\gamma}$ as $n\to\infty$ under $H_{0}$ is that of

[TABLE]

Proof of Theorem 2.1. To highlight the main idea of the proof, let

[TABLE]

and put

[TABLE]

Notice that $Z_{1}^{\circ},Z_{2}^{\circ},\ldots$ are i.i.d. (centred) random elements of ${\cal H}$. Moreover, the condition $\gamma>2$ implies $\mathbb{E}\|Z_{1}^{\circ}\|^{2}_{\oplus}<\infty$. By a Hilbert space central limit theorem (CLT), see e.g. [26], there is some centred Gaussian random element $W^{\circ}$ of ${\cal H}$ such that $W_{n}^{\circ}\stackrel{{\scriptstyle{\cal D}}}{{\longrightarrow}}W^{\circ}$. The idea now is to approximate $W_{n}$ by a suitable random element $\widetilde{W}_{n}$ of ${\cal H}$ so that $\|W_{n}-\widetilde{W}_{n}\|_{\oplus}=o_{\mathbb{P}}(1)$ as $n\to\infty$, and $\widetilde{W}_{n}(t)=n^{-1/2}\sum_{j=1}^{n}\widetilde{Z}_{j}(t)$, where $\widetilde{Z}_{1},\widetilde{Z}_{2},\ldots$ are i.i.d. centred random elements of ${\cal H}$ satisfying $\mathbb{E}\|\widetilde{Z}_{1}\|^{2}_{\oplus}<\infty$. The assertion would then follow from the CLT in Hilbert spaces and Slutzky’s lemma. To this end, put

[TABLE]

A Taylor expansion gives

[TABLE]

where $|\Theta_{n,j}|\leq 1$. By some algebra, it follows that

[TABLE]

where

[TABLE]

Notice that the first term on the right-hand side of (2.6) is $W_{n}^{\circ}(t)$, as given in (2.3). By Proposition 8.1, we have $\|C_{n}\|_{\oplus}=o_{\mathbb{P}}(1)$, and Proposition 8.2 yields

[TABLE]

In view of Proposition 8.3, (2.6) implies $W_{n}(t)=n^{-1/2}\sum_{j=1}^{n}\widetilde{Z}_{j}(t)+o_{\mathbb{P}}(1)$, where

[TABLE]

A straightforward calculation shows $\mathbb{E}\widetilde{Z}_{j}(t)=0$, $t\in\mathbb{R}^{d}$. Moreover, due to the condition $\gamma>2$, we have $\mathbb{E}\|\widetilde{Z}_{j}\|^{2}_{\oplus}<\infty$. Hence, $\widetilde{Z}_{1},\widetilde{Z}_{2},\ldots$ are i.i.d. centred random elements of ${\cal H}$, and the CLT in Hilbert spaces yields the assertion, since the kernel $K$ figuring in (2.2) is given by $K(s,t)=\mathbb{E}[\widetilde{Z}_{1}(s)\widetilde{Z}_{1}(t)^{\top}]$ (for details in computing $K(s,t)$, see Proposition 8.4).

The following result provides some information on the limit null distribution of $T_{n,\gamma}$.

Theorem 2.3

Let $T_{\infty,\gamma}$ be a random variable with the distribution of $\|W\|^{2}_{\oplus}$, where $W$ is given in Theorem 2.1. We then have

[TABLE]

Proof. By Fubini’s theorem, it follows that

[TABLE]

Writing $\textrm{tr}(D)$ for the trace of a square matrix $D$, we have

[TABLE]

Now, some straightforward algebra yields the assertion.

In the univariate case $d=1$, we also obtained an explicit expression for the variance of $T_{\infty,\gamma}$, by using the relation

[TABLE]

By tedious calculations, it follows that

[TABLE]

where $\beta=\gamma-1$, $\delta=(\beta^{2}-1)^{-1/2}$ and $\eta=(4\beta^{2}-1)^{-1/2}$.

Table 1 contains expectation and variance of $T_{\infty,\gamma}$ in the univariate case for various values of $\gamma$.

3 Consistency

In this section, let $X$ have an absolutely continuous distribution, and suppose that $m(t)=\mathbb{E}[\exp(t^{\top}X)]$ exists for each $t\in\mathbb{R}^{d}$. In view of affine invariance, let w.l.o.g. $\mathbb{E}[X]=0$ and $\mathbb{E}[XX^{\top}]=\textrm{I}_{d}$.

Theorem 3.1

We have

[TABLE]

Proof. Fix $K>0$, and put $S(K):=\{t\in\mathbb{R}^{d}:\|t\|\leq K\}$. Let $B(K)$ be the Banach space of continuous functions $f:S(K)\rightarrow\mathbb{R}$, equipped with the norm $\|f\|_{\infty}=\max_{\|t\|\leq K}|f(t)|$. Recall $M_{n}(t)$ from (1.4), and put

[TABLE]

Let $\xi_{n}=\max_{j=1,\ldots,n}\|\Delta_{n,j}\|$, where $\Delta_{n,j}$ is given in (2.4). From (2.4) we obtain

[TABLE]

Since the existence of $m(t)$ implies $\mathbb{E}\|X\|^{4}<\infty$, Theorem 5.2 of [3] yields

[TABLE]

From $S_{n}-\textrm{I}_{d}=n^{-1}\sum_{j=1}^{n}(X_{j}X_{j}^{\top}-\textrm{I}_{d})-\overline{X}_{n}\overline{X}_{n}^{\top}$ and Kolmogorov’s variance criterion for averages (see [24], p. 73), we have $\lim_{n\to\infty}n^{1/2-\varepsilon}\|S_{n}-\textrm{I}_{d}\|_{2}=0$ $\mathbb{P}\textrm{-a.s.}$ for each $\varepsilon>0$. Hence Proposition 8.5 implies

[TABLE]

Using the notation $\|f\|_{\infty}=\max_{\|t\|\leq K}\|f(t)\|$ also for a function $f:\mathbb{R}^{d}\to\mathbb{R}^{d}$, (2.5) gives

[TABLE]

By the strong law of large numbers (SLLN) in $B(K)$, the first factor on the right-hand side converges almost surely to $\|m\|_{\infty}$, and thus (3.4) entails

[TABLE]

Putting

[TABLE]

the triangle inequality gives

[TABLE]

Invoking (3.5), (3.2), (3.3) and Proposition 8.5, we have

[TABLE]

Writing id for the identity function on $\mathbb{R}^{d}$, the triangle inequality yields

[TABLE]

From (3.6), (3.8) and the SLLN in $B(K)$ it follows that

[TABLE]

Likewise, we have

[TABLE]

and thus

[TABLE]

Upon combining (3.9) and (3.10), and using Fatou’s lemma, we obtain

[TABLE]

Since $K$ was arbitrary, the assertion follows.

Remark 3.2

If $X$ has some non-degenerate non-normal distribution with existing moment generating function, then $m^{\prime}(t)\neq tm(t)$ for at least one $t$. Consequently, the right-hand side of (3.1) is strictly positive, and thus

[TABLE]

Hence, the test is consistent against each such alternative. We conjecture that (3.11) holds for any non-normal alternative distribution. A proof of such a result, however, remains an open problem.

4 The limit $\gamma\to\infty$

In this section, we will show that, for fixed $n$, the statistic $T_{n,\gamma}$, after a suitable scaling, approaches a linear combination of two well-known measures of multivariate skewness as $\gamma\to\infty$.

Theorem 4.1

We have (pointwise on the underlying probability space)

[TABLE]

where

[TABLE]

is nonnegative invariant sample skewness in the sense of [27], and

[TABLE]

denotes sample skewness in the sense of Móri, Rohatgi and Székely, as defined in [29].

Proof. We start with (1.7) and use

[TABLE]

as $\gamma\to\infty$. Multiplying this expression with the term within the rightmost bracket of (1.7), and using the relations $\sum_{j=1}^{n}Y_{n,j}=0$, $\sum_{j=1}^{n}\|Y_{n,j}\|^{2}=nd$, $\sum_{j,k=1}^{n}\|Y_{n,j}+Y_{n,k}\|^{2}=2n^{2}d$,

[TABLE]

as well as

[TABLE]

the result follows by tedious but straightforward algebra.

Remark 4.2

It is interesting to note a similarity between Theorem 4.1 and Theorem 2.1 of [15], who showed that the BHEP-statistic for testing for multivariate normality, after suitable rescaling, approaches the linear combination $2b_{1,d}+3\widetilde{b}_{1,d}$ as a smoothing parameter (called $\beta$ in that paper) tends to [math]. Since $\beta$ and $\gamma$ are related by $\beta=\gamma^{-1/2}$, this corresponds to letting $\gamma$ tend to infinity. The same linear combination $2b_{1,d}+3\widetilde{b}_{1,d}$ also showed up as a limit statistic in [18]. Notice that, in the univariate case, the limit statistic figuring in Theorem 4.1 is nothing but three times squared sample skewness. It should be stressed that tests for multivariate normality based on $b_{1,d}$ or $\widetilde{b}_{1,d}$ (or on related measures of multivariate skewness and kurtosis) lack consistency against general alternatives, see, e.g., [4], [5], [13], and [14].

5 Monte Carlo results

In this section we compare the finite-sample power performance of the newly proposed test to those of several competing tests for normality, both for the univariate and the multivariate case. In the case $d=1$ the competing procedures are

1. a)

the Cramér-von Mises ($CvM$) test, 2. b)

the Anderson-Darling ($AD$) test, 3. c)

the Shapiro-Wilk ($SW$) test, 4. d)

the Jarque-Bera ($JB$) test, 5. e)

the Zghoul ($Z$) test.

The R package $nortest$ contains the functions cvm.test and ad.test, which can be used in order to calculate the test statistic and the associated p-value for each of the first two tests mentioned above, see [12]. The Shapiro-Wilk test can be performed using the $Shapiro.test$ function included in the $stats$ package. The R package $tseries$ contains a function called $jarque.bera.test$, which can be used in order to calculate the test statistic and p-value associated with the fourth test mentioned above, see [34]. Each of these tests are well-known and will not be discussed further.

The test of Zghoul (see [36]) is based on the empirical moment generating function. [36] includes a Monte Carlo study which indicates that the finite-sample power performance of the test compares favorably to that of its competitors, especially against symmetric alternatives with kurtosis higher than that of the normal distribution. However, the mentioned paper fails to provide the mathematical theory underlying this test. [19] more recently provided this theory, including a proof that the test is consistent against general alternatives.

The test statistic of Zghoul is a weighted $L^{2}$-statistic, given by

[TABLE]

where $\gamma>2$ is a smoothing parameter. Based on the finite-sample performance reported in [36], the author recommended setting $\gamma$ equal to either $3$ or $15$. The numerical results presented below include the powers obtained using both of these values for the smoothing parameter; the resulting tests are denoted by $Z_{3}$ and $Z_{15}$ respectively. The test statistic $Z_{n}(\gamma)$ can be rewritten in the computationally amenable form

[TABLE]

This test rejects normality for large values of the test statistic.

Turning our attention to the multivariate case, we consider the finite-sample power performance of the newly proposed test to those of some prominent competing tests, and to two very recent tests for multinormality. These procedures are

1. a)

Mardia’s tests based on skewness and kurtosis, 2. b)

the Henze-Zirkler test, 3. c)

the energy test of Székely and Rizzo, 4. d)

a recent test of Henze, Jiménez-Gamero and Meintanis, 5. e)

a recent test of Henze and Jiménez-Gamero.

a): **Mardia’s tests based on skewness and kurtosis

** Mardia’s test for multinormality based on sample skewness rejects $H_{0}$ for large values of $b_{1,d}$, where $b_{1,d}$ is given in (4.1). Notice that $b_{1,d}$ is a consistent estimator of $\beta_{1,d}=E(X_{1}^{\top}X_{2})^{3}$. Under normality we have $\beta_{1,d}=0$, and the limit distribution of $nb_{1,d}$ as $n\rightarrow\infty$ is $6\chi^{2}_{d(d+1)(d+2)/6}$, see [27]. The limit distribution of $nb_{1,d}$ under elliptical symmetriy has been derived by [5].

Sample kurtosis in the sense of Mardia is given by

[TABLE]

which is an estimator of $\beta_{2,d}=\mathbb{E}\|X_{1}\|^{4}$.

Under normality, we have $\beta_{2,d}=d(d+2)$, and the limit null distribution of $\sqrt{n}(b_{2,d}-d(d+2))$ is the normal distribution $\textrm{N}(0,8d(d+2))$, see [27]. The test based on $b_{2,d}$ rejects $H_{0}$ for large and small values of $\beta_{2,d}$.

The R package $QuantPsyc$ contains a function $mult.norm$, which calculates both Mardia’s skewness and kurtosis measures as well as the p-values associated with the corresponding tests from multivariate normality, see [11]. Below we denote the tests based on Mardia’s skewness and kurtosis by $MS_{n}$ and $MK_{n}$, respectively.

We stress that there are several other measures of skewness and kurtosis, see Sections 3 and 4 of [16]. The deficiences of such measures as statistics for purportly “directed tests” for multivariate normality have been addressed in [4], [5], [13] as well as [14].

b): **The Henze-Zirkler test

** Writing $\Psi_{n}(t)=n^{-1}\sum_{k=1}^{n}\exp(\textrm{i}t^{\top}Y_{n,k})$ for the empirical characteristic function of the scaled residuals $Y_{n,1},\ldots,Y_{n,n}$, [21] proposed the test statistic

[TABLE]

for some fixed $\gamma>0$. The test statistic can be rewritten as

[TABLE]

The Henze-Zirkler test is obtained by setting

[TABLE]

This choice of $\gamma$ corresponds to the optimal bandwidth for a multivariate nonparametric density estimator with a Gaussian kernel. The Henze-Zirkler test is included because of its impressive power performance reported in previous studies, see, e.g., [28].

The Henze-Zirkler test (denoted $HZ_{n}$ below) is programmed in the function $hzTest$ in the R package $MVN$, see [25]. This test rejects normality for large values of the test statistic. $HZ_{n}$ is equivalent to a test by [7], as remarked in [16].

c): **The energy test

** [33] proposed the test statistic

[TABLE]

Here, the first expectation is with respect to a random vector $Z$, which is independent of $Y_{n,j}$ and has the distribution $\textrm{N}_{d}(0,\textrm{I}_{d})$, and

[TABLE]

This test, denoted by $EN_{n}$, is known as the energy test. Rejection of $H_{0}$ is for large values of $EN_{n}$. The function $mvnorm.etest$ in the R package $energy$ calculates this test statistic as well as the corresponding p-value, see [31]. The energy test is also reported to have excellent power performance, see [23].

d): **The Henze–Jiménez-Gamero–Meintanis test

** By generalizing a characterization of normality involving both the characteristic and the moment generating function (see [35]) to the multivariate case, [18] proposed to base a test of $H_{0}$ on the weighted $L^{2}$-statistic

[TABLE]

for some parameter $\gamma>1$. Putting $Y_{jk}^{\pm}=Y_{n,j}\pm Y_{n,k}$, $HM_{n}$ can be rewritten as

[TABLE]

The test based on $HM_{n}$ uses an upper rejection region.

A disadvantage of this test is the need to calculate a four-fold sum in order to evaluate the test statistic. Hence, the amount of computer time required in order to perform this test is substantially more than the time required for the other tests under consideration.

e): **The Henze–Jiménez-Gamero test

** [17] present a multivariate generalization of a recent class of tests for univariate normality (see [19]) based on the empirical moment generating function. The test statistic is

[TABLE]

where $\beta>1$ is a fixed parameter. An alternative representation of $HJ_{n}$ is

[TABLE]

which is amenable to computation. This test rejects $H_{0}$ for large values of the test statistic.

The numerical results below were obtained using the software package R, see [30]. Monte Carlo simulation was used in order to estimate the upper percentiles of the distribution of $T_{n,\gamma}$ for the sample sizes $n=20,50,100,200,300$ and dimensions $d=1,2,3,5$. Table 2 provides the empirical $95$-percentiles of $(16\gamma^{2+d/2}/\pi^{d/2})T_{n,\gamma}$ obtained by simulating one million samples from a $\textrm{N}_{d}(0,\textrm{I}_{d})$ distribution using $\gamma=2.5,3,4,5,7,10,\infty$. Notice that, for ease of comparison, the values associated with $\gamma=\infty$ in the rightmost column of Table 2 have been multiplied by 100, since these values are quite close to 0. The tables presented were obtained using the R package $stargazer$, see [22]. From Table 2, we see that the speed of convergence to the asymptotic null distribution is slow and depends on the value of $\gamma$ as well as on the dimension of the data.

5.1 Power results

In the power study presented in this section we simulate from several univariate and multivariate distributions. The critical values presented in Table 2 are used in order to calculate the power results for the newly proposed test. For the remaining tests, the critical values used are also calculated based on $10^{6}$ simulations under $H_{0}$. The power estimates presented below are based on 10 000 random samples. Notice that, due to the computational complexity of $HM_{n}$, the number of simulations used in order to obtain both critical values and power estimates for this test are reduced by a factor of 10. This step is necessary in order to ensure that the numerical results can be obtained within a reasonable time.

Power results are reported for a sample size of $n=50$ and data dimensions $d=1,2,3,5$. A nominal significance level of $\alpha=0.05$ is used throughout. As is pointed out in [28], the maximum possible standard error for each power estimate is $0.005$. Thus, the 99% confidence interval for each reported power estimate is contained in the interval obtained by subtracting 1% from, and adding 1% to, the stated power.

The results for the univariate case are given in Table 3. The entries are percentages of rejection of $H_{0}$, rounded to the nearest integer. The power against the standard normal distribution shows that the nominal level is maintained very closely. As for the alternative distributions, NMIX1 and NMIX2 denote mixtures of the normal distributions $\textrm{N}(0,1)$ and $\textrm{N}(0,4)$. The mixture NMIX1 gives equal weight to these distributions, while NMIX2 is obtained when the probability of sampling from the standard normal distribution is 0.75. The remaining alternative distributions considered are $t$-distributions with 3, 5 and 10 degress of freedom, the lognormal distributions with parameters $(0,1/2)$ and $(0,1/4)$ (denoted by $LN(\cdot)$), the $\chi^{2}$-distributions with 5 and 15 degrees of freedom, the standard logistic distribution, the Weibull distributions with shape parameters 10 and 20, the Pearson type VII distributions with 5 and 10 degrees of freedom (denoted by $P_{VII}(\cdot)$), and the skew-normal law with skewness parameters 3 and 5 (denoted by $SN(\cdot)$), see [1].

The R Package $PearsonDS$ contains the function $rpearsonVII$, which can be employed to simulate random variables from this distribution, see [6]. The R Package $sn$ disposes of the function $rsn$, which can be used to simulate random variates from the skew normal distribution, see [2].

Tables 3,4, 5 and 6 report the powers calculated in the cases where $d$ equals 1, 2, 3 and 5 respectively. Note that the subscript $n$ in the names of the test statistics is omitted in the tables in order to save space.

The results shown in Table 3 indicate that the newly proposed class of tests exhibit substantial power against the distributions considered. $T_{\infty}$ outperforms each of the competing tests in terms of the estimated powers against the $LN(0,1/4)$, $\chi^{2}(15)$ and $SN(3)$ laws, as well as both of the Weibull distributions considered. The newly proposed tests also proves to be serious competitors against each of the remaining distributions considered, especially for small values of $\gamma$ (in which case $T_{n,\gamma}$ is often only outperformed by the Jarque-Bera test).

We now turn our attention to the case $d>1$. As was the case for the univariate tests, the powers against 16 alternative distributions are reported for each of the data dimensions considered. The powers of each of the tests against the standard normal distribution are also included in the relevant tables. The alternative distributions considered include mixtures of normal laws, distributions with independent marginals, distributions where each marginal is normal except for one, a spherically symmetric distribution, and a distribution for which every marginal is normal, but the joint distribution does not follow the normal law.

The parameter combinations used for the mixtures of normal distributions were taken from [33]. Let $p\textrm{N}_{d}(\mu_{1},\Sigma_{1})+(1-p)\textrm{N}_{d}(\mu_{2},\Sigma_{2})$ denote a normal mixture, where the probability of sampling from $\textrm{N}_{d}(\mu_{1},\Sigma_{1})$ is $p$ and the probability of sampling from $\textrm{N}_{d}(\mu_{2},\Sigma_{2})$ is $1-p$. Let $\mu=0$ and $\mu=3$ denote $d$-dimensional column vectors of 0’s and 3’s, respectively, and let $B_{d}$ denote a ($d\times d$)-matrix containing 1’s on the main diagonal and 0.9’s on each off diagonal entry. The normal mixtures are constructed by combining $\textrm{N}_{d}(0,\textrm{I}_{d})$, $\textrm{N}_{d}(3,\textrm{I}_{d})$ and $\textrm{N}_{d}(0,B_{d})$. The first mixture, denoted by NMIX1, is $0.9\textrm{N}_{d}(0,\textrm{I}_{d})+0.1\textrm{N}_{d}(3,\textrm{I}_{d})$. This distribution is skewed with heavy tails. The second mixture, denoted by NMIX2, is $0.9\textrm{N}_{d}(0,B_{d})+0.1\textrm{N}_{d}(0,\textrm{I}_{d})$. This is a symmetric, heavy-tailed distribution. In addition, we included two multivariate $t$-distributions; the $t_{\nu}(0,\textrm{I}_{d})$-distribution for $\nu=5$ and $\nu=10$. Next, we included distributions with independent marginals, the latter being the $\chi^{2}$-distribution with 15 and 20 degrees of freedom respectively, the $logistic(0,1)$ distribution, the gamma distributions with parameters $(5,1)$ and $(4,2)$, as well as the Pearson Type VII distributions with 10 and 20 degrees of freedom.

Three $d$-dimensional distributions are obtained by combining $d-1$ independent standard normal marginals with one non-normal distribution. This distribution is denoted by $\textrm{N}(0,1)^{d-1}\otimes F$, where $F$ denotes the non-normal marginal distribution. The three alternatives considered for $F$ are the $\chi^{2}$-distributions with 5 and 10 degrees of freedom, respectively, as well as the $t$-distribution with 3 degrees of freedom.

Spherically symmetric distributions can be defined in R using the $EllipticalDistribution$ function from the R package $distrEllipse$, see [32]. Tables 4–6 display the estimated powers of the various tests considered against the $d$-dimensional spherically symmetric distribution, where the radius of the distribution follows a lognormal distribution with parameters 0 and 0.5. This distribution is denoted by $\mathcal{S}^{d}(LN(0,1/2))$.

Let $\rho_{d}$ and $\rho_{d}^{\prime}$ denote positive definite ($d\times d$)-matrices with 1’s on the main diagonal, where $\rho_{d}$ has the constant $\rho$ and $\rho_{d}^{\prime}$ the constant $-\rho$ on each off diagonal entry. The final distribution considered is the mixture $0.5\textrm{N}_{d}(0,\rho_{d})+0.5\textrm{N}_{d}(0,\rho_{d}^{\prime})$. This distribution is a non-normal $d$-variate distribution with normal marginals.

As was the case in the univariate setting, the newly proposed test is associated with several high powers reported in Tables 4, 5, and 6. When comparing the results for the distributions with independent marginals, we see that $T_{\infty}$ outperforms each of the competitors against both of the distributions with $\chi^{2}$ marginals as well as both of the distributions with gamma marginals considered. This is also the case aginst the \textrm{N}(0,1)^{d-1}\otimes$$\chi^{2}(5) and \textrm{N}(0,1)^{d-1}\otimes$$\chi^{2}(10) distributions. The mentioned predominance is for each of the data dimensions considered. Furthermore, the newly proposed test statistic shows high power against the remaining distributions for finite values of $\gamma$. Specifically, when $d=2$ the new test outperforms its competitors against the spherically symmetric distribution. In the case $d=3$, none of the competing tests are able to outperform the newly proposed class of tests against the second normal mixture considered, the $t$-distribution with 10 degrees of freedom or the $PVII(10)$ distribution. The corresponding list of distributions in the case where $d=5$ is obtained by substituting the $t$-distribution with 10 degrees of freedom for the $t$-distribution with 5 degrees of freedom and adding the $PVII(20)$ distribution. Finally, none of the competing tests is able to outperform the newly proposed class of tests against the $NM_{d}(\rho=0.2)$ distribution.

In most of the cases considered, the power of the newly proposed class of tests is a monotone function of $\gamma$. Based on the numerical results presented, it is recommeded that $\gamma=5$ be used when performing the test as this value results in reasonably high power against the majority of the alternatve distributions considered.

6 A real data example

The payoff function of certain types of financial derivatives depends on the joint behaviour of multiple stocks or indexes; an important example is the class of basket options. When calculating the price of a basket option, it is often assumed that the log-returns of the stocks or indexes considered are realized from a multivariate normal distribution (this assumption is an extention of the celebrated Black-Merton-Scholes model for options on a single stock or index). As a result, testing the hypothesis that observed financial log-returns follow a multivariate normal law is of interest when pricing basket options. For more details regarding the pricing of these options, the interested reader is referred to [8].

As a practical application, we consider the log-returns associated with three major indexes traded in the financial market of the United States. 50 daily log-returns were calculated for the period ending 29 December 2017, the relevant prices were downloaded from http://finance.yahoo.com. The first index considered is the Dow Jones Industrial Average (DJIA) index; this index is comprised of a price-weighted average of 30 large publicly owned companies. The second is the Standard & Poor 500 (S&P 500); a market-capitalization weighted index comprising 500 large companies. Finally, we consider the log-returns of the National Association of Securities Dealers Automated Quotations (NASDAQ) composite index. We are interested in testing the hypothesis that the log-returns are realized from a multivariate normal distribution using the newly proposed test.

Table 7 shows the estimated p-values associated with the newly proposed tests for various values of $\gamma$, the reported p-values were obtained using one million Monte Carlo simulations in each case. The results indicate that the hypothesis of multivariate normality is rejected at a 1% significance level for each value of $\gamma$ considered.

7 Conclusion

We proposed and studied a new class of affine invariant tests for normality in any dimension that are based on a partial differential equation involving the moment generating function. Some properties of the limit null distribution of the test statistic $T_{n,\gamma}$ have been derived, and the consistency of this class of tests against general alternatives has been proved under some mild conditions. For fixed $n$, the test statistic $T_{n,\gamma}$, after suitable scaling, approaches a linear combination of two measures of multivariate skewness as $\gamma\to\infty$.

A Monte Carlo study investigates the finite-sample performance of $T_{n,\gamma}$ compared to those of competing tests in the univariate and multivariate settings. The competing tests considered for univariate normality comprise four well-known tests, while, in the multivariate case, we include four prominent classical tests for multinormality and two very recent tests. The numerical results indicate that $T_{n,\gamma}$ often exhibits power greater than those associated with several of its competitors, both in univariate and multivariate settings. Based on the numerical results obtained, it is recommended that $\gamma=5$ is used when performing the test.

Acknowledgement: The authors thank Tobias Jahnke for providing the proof of Proposition 8.5. The second author’s work is based on research supported by the National Research Foundation, South Africa (Research chair: Non-parametric, Robust Statistical Inference and Statistical Process Control, Grant number 71199). Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF.

8 Appendix: Technical proofs

Proposition 8.1

We have $\|C_{n}\|_{\oplus}=o_{\mathbb{P}}(1)$, where $C_{n}$ is given in (2.9).

Proof. Putting $\xi_{j}(t):=\exp(t^{\top}X_{j}\left(t^{\top}\Delta_{n,j}\right)^{2}\exp\left(\Theta_{n,j}t^{\top}\Delta_{n,j}\right)/2$, we have

[TABLE]

Recall $F_{n}$ from (3.7) and put $\Lambda_{n}:=\|S_{n}^{-1/2}-\textrm{I}_{d}\|_{2}F_{n}+\|S_{n}^{-1/2}\|_{2}\|\overline{X}_{n}\|$. We have $|\Theta_{n,j}t^{\top}\Delta_{n,j}|\leq\|t\|\,\Lambda_{n}$. Furthermore, using $|t^{\top}(X_{j}+X_{k})|\leq 2\,\|t\|\,F_{n}$ and

[TABLE]

as well as

[TABLE]

we obtain

[TABLE]

where $\Gamma_{n}=F_{n}+\Lambda_{n}$. Expanding the curly brackets, the leading terms are those that do not involve any of $\Delta_{n,j}$. We concentrate on

[TABLE]

which originates from choosing the first term within each of the curly brackets. The other terms are treated similarly. Notice that

[TABLE]

and that $\|S_{n}^{-1/2}-\textrm{I}_{d}\|^{2}_{2}=O_{\mathbb{P}}(n^{-2})$, $n^{-1}\sum_{j=1}^{n}\|X_{j}\|=O_{\mathbb{P}}(1)$. By Proposition 10.2 of [18], the integral is of order $O_{\mathbb{P}}(\Gamma_{n}^{d+3})\exp(\Gamma_{n}^{2}/\gamma)$ (notice that $1+\gamma$ in that paper corresponds to (our) $\gamma$). From display (10.6) and display (10.7) of [18] we have $\Gamma_{n}=O_{\mathbb{P}}(\sqrt{\log n})$ and $\exp\left(\Gamma_{n}^{2}/\gamma\right)=n^{2/\gamma}\cdot\left(\log n\right)^{(d-2)/\gamma}\,O_{\mathbb{P}}(1)$. Since $\gamma>2$, it follows that

[TABLE]

Proposition 8.2

For $B_{n}$ given in (2.8), we have

[TABLE]

Proof. Observe that $B_{n}(t)=n^{-1/2}\sum_{j=1}^{n}\exp(t^{\top}X_{j})\Delta_{n,j}+R_{n}(t)$, where

[TABLE]

Use

[TABLE]

and $\|\Delta_{n,j}\|\leq\|S_{n}^{-1/2}-\textrm{I}_{d}\|_{2}\,\|F_{n}\|+\|S_{n}^{-1/2}\|\,\|\overline{X}_{n}\|$ with $F_{n}$ given in (3.7) together with $F_{n}=O_{\textrm{P}}(\sqrt{\log n})$ (see Prop. 10.1. of [18]) and

[TABLE]

to show $R_{n}=o_{\mathbb{P}}(1)$. Next, $B_{n}(t)-R_{n}(t)=B_{n,1}(t)-B_{n,2}(t)-B_{n,3}(t),$ where

[TABLE]

Since

[TABLE]

and $\|\sqrt{n}\left(S_{n}^{-1/2}-\textrm{I}_{d}\right)\|_{2}^{2}\cdot\|\overline{X}_{n}\|^{2}=O_{\mathbb{P}}(n^{-1})$, one may use (8.1) and Fubini’s theorem to show $B_{n,2}=o_{\mathbb{P}}(1)$. To proceed, rewrite $B_{n,3}(t)$ in the form

[TABLE]

Since replacing the first factor with its expectation $\exp(\|t\|^{2}/2)$ means adding a term that is asymptotically negligible, we have

[TABLE]

To tackle $B_{n,1}(t)$, we rewrite its transpose in the form

[TABLE]

and use display (2.13) of [20], according to which

[TABLE]

where $D_{n}=O_{\mathbb{P}}(n^{-1/2})$. Putting

[TABLE]

and $Y_{n}(t)=n^{-1}\sum_{j=1}^{n}\left(\exp(t^{\top}X_{j})X_{j}-E(t)\right)$, we have

[TABLE]

Abbreviating the matrix within curly brackets by $S_{n}$, we have $\|S_{n}Y_{n}(t)\|^{2}\leq\|S_{n}\|^{2}_{2}\,\|Y_{n}(t)\|^{2}$. Since $\|S_{n}\|^{2}_{2}=O_{\mathbb{P}}(1)$, it follows that

[TABLE]

Now, observe that $Y_{n}(t)$ is a mean of centred random vectors, and invoking Fubini’s theorem the expectation of the integral is seen to be of order $O(n^{-1})$. Thus, $\|S_{n}Y_{n}\|_{\oplus}=o_{\mathbb{P}}(1)$, and hence (since the matrix $D_{n}$ figuring in (8.5) is asymptotically negligible) we have

[TABLE]

Upon combing this result with (8.2) and (8.4), the assertion follows.

Proposition 8.3

For $A_{n}$ given in (2.7), we have

[TABLE]

Proof. Notice that $A_{n}(t)=A_{n,1}(t)-A_{n,2}(t)-A_{n,3}(t)$, where

[TABLE]

To show $A_{n,2}=o_{\mathbb{P}}(1)$, use n\big{\|}S_{n}^{-1/2}-\textrm{I}_{d}\big{\|}^{2}_{2}\,\|\overline{X}_{n}\|^{2}=O_{\mathbb{P}}(n^{-1}) together with

[TABLE]

and Fubini’s theorem, since \mathbb{E}\big{\|}n^{-1}\sum_{j=1}^{n}\exp(t^{\top}X_{j})(X_{j}-t)\big{\|}^{2}=O(n^{-1}), due to the fact that the summands are centred random vectors. Likewise,

[TABLE]

Since each of the summands is a product of centred random vectors or random variables, we have $\mathbb{E}\|A_{n,3}(t)\|^{2}=O(1/n)$, and Fubini’s theorem yields $A_{n,3}=o_{\mathbb{P}}(1)$. To conclude the proof, observe that, by (8.3),

[TABLE]

where $\Delta_{n}(t)=n^{-1}\sum_{i=1}^{n}\exp(t^{\top}X_{i})(X_{i}-t)X_{i}^{\top}$, $V_{j}(t)=X_{j}X_{j}^{\top}t-t$. Notice that $\Delta_{n}(t)$ is a mean of i.i.d. random matrices, and that

[TABLE]

Straightforward calculations show that replacing $\Delta_{n}(t)$ with the right-hand side of the last equation means adding a term that is asymptotically negligible. Hence, the first term on the right-hand side of (8.6) is

[TABLE]

The second term is $o_{\mathbb{P}}(1)$ since $D_{n}=O_{\mathbb{P}}(n^{-1/2})$.

Proposition 8.4

(Calculation of $K(s,t)$)*

Recall $\widetilde{Z}_{1}(t)$ from (2.10). Putting $m(s)=\exp(\|s\|^{2}/2)$ and $X_{1}=X$, we have*

[TABLE]

Writing $\mathbf{0}$ for the zero matrix of order $d$, we have

[TABLE]

The assertion now follows from straightforward calculations.

Proposition 8.5

Let $(A_{n})$ be a sequence of symmetric positive definite ($d\times d$)-matrices and $(b_{n})$ an increasing sequence of positive real numbers satisfying $\lim_{n\to\infty}b_{n}=\infty$ so that

[TABLE]

We then have $\lim_{n\to\infty}b_{n}\|A_{n}^{-1/2}-{\rm{I}}_{d}\|_{2}=0$.

Proof. Let $\Lambda_{n}=\textrm{diag}(\lambda_{1},\ldots,\lambda_{n})$ be the diagonal matrix consisting of the positive eigenvalues of $A_{n}$ so that $\|A_{n}\|_{2}=\max_{i=1,\ldots,n}\lambda_{i}$ and

[TABLE]

Since the assumptions imply $\|A_{n}-{\rm{I}}_{d}\|_{2}\to 0$, choose $n_{0}$ so large that $\|A_{n}-{\rm{I}}_{d}\|_{2}\leq 1/2$ for each $n\geq n_{0}$. Putting $T_{n}=\textrm{I}_{d}-A_{n}$, we have

[TABLE]

and thus $\|A_{n}^{-1/2}\|_{2}\leq\sqrt{2},n\geq n_{0}$. Now, $A_{n}-\textrm{I}_{d}=(A_{n}^{1/2}+\textrm{I}_{d})A_{n}^{1/2}(\textrm{I}_{d}-A_{n}^{-1/2})$ implies $\textrm{I}_{d}-A_{n}^{-1/2}=A_{n}^{-1/2}(A_{n}^{1/2}+\textrm{I}_{d})^{-1}(A_{n}-\textrm{I}_{d})$, whence

[TABLE]

In view of (8.7) and $\|A_{n}^{-1/2}\|_{2}\leq\sqrt{2},n\geq n_{0}$, the assertion follows.