Testing multivariate normality by zeros of the harmonic oscillator in characteristic function spaces
Philip Dörr Bruno Ebner Norbert Henze






White dwarfs as a probe of dark energy
    B. Novosyadlyj^1,2
    December 30, 2023
======================================



 MSC 2010 subject
classifications. Primary 62E10 Secondary 60E10, 62G10

Key words and phrases test for multivariate normality; affine invariance; consistency; empirical characteristic function; harmonic oscillator; neighborhood-of-model validation



We study a novel class of affine invariant and consistent tests for normality	in any dimension.
The tests are based on a characterization of the standard d-variate normal distribution as the unique solution
of an initial value problem of a partial differential equation motivated by the harmonic oscillator,
which is a special case of a Schrödinger operator.  We derive the asymptotic distribution of the
test statistics under the hypothesis of normality as well as under fixed and contiguous alternatives.
The tests are consistent against general alternatives, exhibit strong power performance for finite samples,
and they are applied to a classical data set due to  R.A. Fisher. The results can also be used for a neighborhood-of-model validation
procedure.




§ INTRODUCTION

The multivariate normal distribution plays a key role in classical and hence widely used procedures.
Serious statistical inference that involves the assumption of multivariate normality should therefore start with
 a test of fit to this model. There is a continuing interest in this testing problem, as evidenced by a multitude of papers.
 The proposed tests may be roughly classified as follows: <cit.> consider
 tests based on the empirical characteristic function, while <cit.> employ the
 empirical moment generating function.
A classical (and still popular) approach is to consider measures of multivariate skewness and kurtosis
(see, e.g., <cit.>), as supposedly diagnostic tools
with regard to the kind of deviation from normality when this hypothesis has been rejected, but the deficiencies of those measures in this regard have been clearly demonstrated
(see, e.g., <cit.>). Other approaches involve generalizations of tests for univariate normality
<cit.>, the examination of nonlinearity of dependence <cit.>, canonical correlations <cit.>
and the notion of energy, see <cit.>. For a survey of affine invariant tests for multivariate normality, see  <cit.>.
 Monte Carlo studies can be found in  <cit.>.

To be specific, let X,X_1,…,X_n,… be a sequence of independent identically distributed (i.i.d.) d-dimensional random (column) vectors,
which are defined on some common probability space (Ω, A,ℙ). Here, d≥1 is a fixed integer, which means that
the univariate case is deliberately not excluded.  We write ℙ^X for the distribution of X. The d-variate normal distribution with expectation
μ and nonsingular covariance matrix Σ will be denoted by N_d(μ,Σ). Furthermore,

    𝒩_d={N_d(μ,Σ):μ∈ℝ^d, Σ∈ℝ^d× d }

stands for the family of nondegenerate d-variate normal distributions. A check of the assumption of multivariate normality means to test the hypothesis

    H_0: ℙ^X∈𝒩_d,

against general alternatives.

Writing I_d for the unit matrix of order d, our novel idea for testing H_0 is to use a characterization of the Fourier transform
 of N_d(0,I_d) as the unique solution of an initial value problem of a partial differential equation
motivated by the harmonic oscillator, which is a special case of a Schrödinger operator.
The proposed test statistic is based on the squared norm of a functional of the empirical characteristic function in a suitably weighted L^2-space.
This statistic is close to zero under the hypothesis (<ref>), and rejection will be for large values of the test statistic.

Let L^2(^d) be the space of square integrable functions, equipped with the usual norm and scalar product ⟨·,·⟩.
Consider for sufficiently regular f ∈ L^2(^d) the partial differential equation

    {[ Δ f(x) = (x^2-d) f(x),   x∈^d,;                        f(0)=1. ].

Here, Δ stands for the Laplace operator, and · denotes the Euclidean norm. Notice that we can rewrite (<ref>) as

(-Δ + x^2 - d) f(x)= 0

or, equivalently, as

    ∑_j=1^d(-∂^2/∂ x_j^2 + x_j^2 - 1) f(x) = 0,   x=(x_1,…,x_d) ∈^d.

  The operator -Δ + x^2 - d is known as the harmonic oscillator, see <cit.>. In the univariate case,
   (<ref>) reduces to a fixed point problem (or, equivalently, to the problem of finding the eigenfunction that corresponds
   to the eigenvalue 1) of the Hermite operator, see equation (1.1.9) in <cit.>.
   The solution of this problem is the 0^th Hermite function, which coincides with the solution given in the following theorem.



The characteristic function

    ψ(t)=exp(-t^2/2),   t∈^d,

of the d-variate standard normal distribution N_d(0,I_d)
is the unique solution of (<ref>).




Let f be an arbitrary solution of (<ref>). Writing i for the imaginary unit, we introduce the creation and annihilation operators
a_j=x_j+ip_j and a_j^⋆=x_j-ip_j, j=1,…,d,
where p_j=-i∂/∂ x_j, j=1,…,d. For each j ∈{1,…,d} we have

    a_j^⋆ a_j     =    (x_j-∂/∂ x_j)(x_j+∂/∂ x_j) =
     x_j^2+x_j∂/∂ x_j-∂/∂ x_jx_j-∂^2/∂ x_j^2
        =     -∂^2/∂ x_j^2+x_j^2-1.

So we can rewrite (<ref>) as
∑_j=1^da_j^⋆ a_j f = 0,
which implies

    ⟨ f,∑_j=1^da_j^⋆ a_j f⟩ 
    = ∑_j=1^d⟨ f,a_j^⋆ a_j f ⟩ =∑_j=1^d⟨ a_jf, a_j f⟩ =∑_j=1^da_j f^2=0

and thus a_jf=0 for each j ∈{1,…,d}. We therefore have (x+∇)f=0, where ∇ denotes the gradient operator.
By the last statement and the product rule, it follows that

    ∇(exp(x^2/2) f)=xexp(x^2/2)f-xexp(x^2/2)f=0

which, in view of the condition f(0)=1, completes the proof.





  The operator H=-Δ + x^2 is the Hermite operator in ^d, and ψ is the product of the one-dimensional 0^th Hermite functions.
  Therefore, since Hψ=dψ (as we have shown in Theorem <ref>), ψ is the eigenfunction associated with the eigenvalue d.
  For details on the Hermite operator in ^d and corresponding eigenfunctions, see p. 5 of <cit.>.



In this paper, we study a family of affine invariant test statistics for H_0 that is based on the characterization given in Theorem <ref>.
Since the class  𝒩_d is closed under full rank affine transformations, Theorem <ref> does not restrict the scope of the testing problem.
We make the tacit standing assumption that ℙ^X is absolutely continuous with respect to Lebesgue measure, and that
 n ≥ d+1.  Let X_n= n^-1∑_j=1^nX_j denote the sample mean and S_n= n^-1∑_j=1^n(X_j-X_n)(X_j-X_n)^⊤
the sample covariance matrix of X_1,…,X_n, where x^⊤ means transposition of a column vector x.
The assumptions made above guarantee that S_n is invertible almost surely, see <cit.>.
The test statistic will be based on the so-called scaled residuals

    Y_n,j=S_n^-1/2(X_j-X_n),    j =1,…,n,

which represent an empirical standardization of the data. Here, S_n^-1/2 is the unique symmetric
positive definite square root of S_n^-1. The test statistic will be based on the
empirical characteristic function

    ψ_n(t)=1/n∑_j=1^nexp(it^⊤ Y_n,j),   t∈^d,

of Y_n,1,…,Y_n,n.  Notice that an application of the Laplace operator Δ to ψ_n yields

    Δψ_n(t)=-1/n∑_j=1^nY_n,j^2exp(it^⊤ Y_n,j),   t∈^d.



Motivated by (<ref>) and Theorem <ref>, we propose the weighted L^2-statistic

    T_n,a   =   n∫|Δψ_n(t)-Δψ(t) |^2w_a(t) t
       =   n∫|1/n∑_j=1^nY_n,j^2exp(it^⊤ Y_n,j)+(t^2-d)exp(-t^2/2)|^2w_a(t) t,

where

    w_a(t)= exp(-at^2),    t ∈^d,

and a>0 is a fixed constant. Moreover, |z| is the modulus of a complex number z, and integration is, unless otherwise specified, over ^d.
In principle, other weight functions
than  w_a are conceivable in the definition of T_n,a. Since, for c ∈^d (see  <cit.>, p. 3601),

    ∫cos(t^⊤ c)exp(-at^2)t=(π/a)^d/2exp(-c^2/4a),

as well as

    ∫(t^2 - d)cos(t^⊤ c)exp( -(a +1/2)t^2)t  hhhhhhhhhhhhhhhhhhhhhhhhhhh
          hhhhhhhhhhhhhhhh =   -(2π)^d/2(c^2+2da(2a+1))/(2a+1)^2+d/2exp( -c^2/2(2a + 1)),


    ∫(t^2-d)^2exp(-(a+1)t^2)t   =  π^d/2/(a+1)^2+d/2(a(a+1)d^2+d(d+2)/4),

an attractive feature of the choice of w_a is that the test statistic takes the simple form

    T_n,a   =   (π/a)^d/21/n∑_j,k=1^nY_n,j^2Y_n,k^2exp(-1/4aY_n,j-Y_n,k^2)
          - 2(2π)^d/2/(2a+1)^2+d/2∑_j=1^nY_n,j^2(Y_n,j^2+2da(2a+1))exp(-1/2Y_n,j^2/2a+1)
          +  nπ^d/2/(a+1)^2+d/2(a(a+1)d^2+d(d+2)/4),

which is amenable to computational purposes. Moreover, T_n,a depends only on the scalar products Y_n,i^⊤ Y_n,j = (X_i-X_n)^⊤ S_n^-1(X_j-X_n),
where i,j ∈{1,…,n}. This shows that T_n,a is invariant with respect to full rank affine transformations of X_1,…,X_n. Moreover,
not even the square root S_n^-1/2 of S_n^-1 is needed.

The rest of the paper is organized as follows: In Section <ref>, we show that, as the tuning parameter a tends to infinity,
the test statistic T_n,a, after a suitable scaling, converges to a  certain measure of multivariate skewness. On the other hand, a time-honored measure of multivariate
kurtosis emerges as a → 0. Section <ref> presents a basic Hilbert space central limit theorem, which proves beneficial for obtaining the
limit distribution of T_n,a both under H_0 and under fixed alternatives to normality.
In Section <ref>, we derive the limit null distribution of T_n,a as n →∞. Section <ref> considers the behavior of T_n,a with respect to contiguous
alternatives to H_0. In Section <ref>, we show that the test for multivariate normality that rejects H_0 for large values of T_n,a is consistent against general
alternatives. Moreover, the limit distribution of T_n,a under a fixed alternative  distribution is seen to be normal. Since the
variance of this normal distribution can be estimated from the data, there is an asymptotic confidence interval for the  measure of distance from normality
under alternative distributions that is inherent in the procedure. Furthermore, there is the option for a neighborhood-of-model validation procedure.
The results of a simulation study, presented in  Section <ref>, show
that the novel test is strong with respect to prominent competitors.
In Section <ref>, the new tests are applied to the Iris flower data set due to R.A. Fisher.
The paper concludes with some remarks. For the sake of readability, most of the proofs and some
auxiliary results are deferred to Section <ref>. Finally, the following abbreviations, valid for t,x ∈^d, will be used in several sections:

    CS^+(t,x) := cos(t^⊤ x) + sin(t^⊤ x),    CS^-(t,x) = cos(t^⊤ x) - sin(t^⊤ x).











§ THE LIMITS A →∞ AND A → 0

The results of this section show that the class of tests for multivariate normality based on T_n,a is
"closed at the boundaries" a →∞ and a → 0 and thus  shed some light on the tuning parameter a, which figures in the weight function w_a given in (<ref>).
We first consider the case a →∞.

 Elementwise on the underlying probability space, we have

    lim_a→∞2a^d/2+1/nπ^d/2  T_n,a = 1/n^2∑_j,k=1^nY_n,j^2Y_n,k^2Y_n,j^⊤ Y_n,k.



 For short, put Y_j := Y_n,j. From the representation of T_n,a we have

    2a^d/2+1T_n,a/nπ^d/2    =    2a/n^2∑_j,k =1^n Y_j^2Y_k^2exp(-1/4aY_n,j-Y_k^2)
           - (2a/2a + 1)^d/2+12/n(2a + 1)∑_j=1^nY_j^2(Y_j^2+2da(2a + 1))exp(-1/2Y_j^2/2a + 1)
           + (a/a+1)^d/2+12/a+1(a(a+1)d^2+d(d+2)/4)
        =:     U_n,1 - U_n,2 + U_n,3

(say). Since ∑_j=1^n Y_j^2 = nd, an expansion of the exponential function yields

    U_n,1 = 2ad^2 - d/n∑_j=1^n Y_j^4 + 1/n^2∑_j,k=1^n Y_j^2 Y_k^2 Y_j^⊤ Y_k + o(1)

as a →∞. To tackle U_n,2, we use

    (2a/2a + 1)^d/2+1 = ( 1+ 1/2a)^-(d/2+1) = 1 - (d/2+1) 1/2a + O(a^-2)

and, after some algebra, obtain

    U_n,2 = 4ad^2 - (d/2+1)2d^2 - d/n∑_j=1^n Y_j^4 + o(1).

Finally, a binomial expansion of (a/(a+1))^d/2+1 yields

    U_n,3 = 2ad^2 - (d/2+1) 2d^2 + o(1).

Summing up, the assertion follows.





The limit

    b_1,d :=  1/n^2∑_j,k=1^nY_n,j^2Y_n,k^2Y_n,j^⊤ Y_n,k

(say), which figures on the right-hand side of (<ref>),
is a measure of multivariate (sample) skewness,  introduced by  Móri, Rohatgi and Székely, see <cit.>.
A much older time-honored measure of multivariate (sample) skewness is skewness in the sense of Mardia (see <cit.>), which
is given by

    b_1,d := 1/n^2∑_j,k=1^n (Y_n,j^⊤ Y_n,k)^3.

It is interesting to compare Theorem <ref> with similar results found in connection with other weighted L^2-statistics
that have been studied for testing H_0. Thus, by Theorem 2.1 of <cit.>, the time-honored
class of BHEP-statistics for testing for multivariate normality (see <cit.>), after suitable rescaling, approaches the linear combination 2b_1,d + 3 b_1,d,
as a smoothing parameter (called β in that paper) tends to 0. Since β and a are related by β = a^-1/2, this corresponds to letting a tend to infinity.
The same linear combination 2b_1,d + 3 b_1,d also showed up as a limit statistic in <cit.> and <cit.>.
Notice that, in the univariate case, the
limit statistic b_1,d figuring in Theorem <ref> is nothing but three times squared sample skewness.
We stress that tests for multivariate normality based on b_1,d or b_1,d or on related measures of multivariate skewness and kurtosis
lack consistency against general alternatives, see, e.g., <cit.>.


We now consider the case a → 0.  Since, elementwise on the underlying probability space, the expressions in (<ref>) and (<ref>)
have finite limits as a → 0, and since the double sum figuring in (<ref>) converges to ∑_j=1^n Y_n,j^4 as a → 0,
we have the following result.

 Elementwise on the underlying probability space, we have

    lim_a→ 0(a/π)^d/2 T_n,a =  1/n∑_j=1^n Y_n,j^4.




The limit statistic on the right-hand side of (<ref>) is Mardia's celebrated measure
b_2,d of multivariate sample kurtosis, see <cit.>. Together with Theorem <ref>, this result shows that,
just like the class of BHEP tests for multivariate normality (see <cit.>), also the class of tests based on T_n,a is "closed at the boundaries"
a →∞ and a → 0. Notably, Mardia's measure of kurtosis shows up for the first time in connection with limits of
weighted L^2-statistics for testing for multivariate normality. The corresponding limit statistic for the class of BHEP tests is, up to a linear transformation,
n^-1∑_j=1^n exp(-Y_n,j^2/2), see Theorem 3.1 of       <cit.>.











§ A BASIC HILBERT SPACE CENTRAL LIMIT THEOREM

In this chapter, we present a basic Hilbert space central limit theorem. This theorem implies the limit distribution of  T_n,a under the null hypothesis  (<ref>),
but it is also beneficial for proving a limit normal distribution of T_n,a under fixed alternatives to H_0. Throughout this section,
we assume that the underlying distribution satisfies X^4 < ∞. Moreover, we let (X) =0 and (X X^⊤) =  I_d in view of affine invariance.
To motivate the benefit of a Hilbert space setting and for later purposes, it will be convenient to represent T_n,a in a different way.



Recall ψ(t) from (<ref>), and
let

    m(t) := (d-t^2) ψ(t),    t ∈^d,


    Z_n(t) := 1/√(n)∑_j=1^n {Y_n,j^2 (cos (t^⊤ Y_n,j) + sin(t^⊤ Y_n,j)) - m(t) },    t ∈^d.

We then have

    T_n,a = ∫ Z_n^2(t)   w_a(t)   d t.




The proof follows by straightforward algebra using the addition theorems for the sine function and the cosine function and the fact that
∫sin(t^⊤ y) m(t) w_a(t)  d t=0, y ∈^d.


A convenient setting for asymptotics will be the separable Hilbert space  of (equivalence classes of) measurable functions f:^d →
satisfying ∫ f^2(t) w_a(t)   dt < ∞. The scalar product and the  norm in  will be denoted by

    ⟨ f,g ⟩_ = ∫ f(t)g(t)   w_a(t)   dt,   f_ = ⟨ f,f ⟩_^1/2,    f,g ∈,

respectively. Notice that

    T_n,a = Z_n_^2.

Recalling CS^±(t,x) from (<ref>), and putting

    μ(t) := [X^2  CS^+(t,X) ],    t ∈^d,

the main object of this section is the random element V_n of , defined by

    V_n(t) := 1/√(n)∑_j=1^n ( Y_n,j^2  CS^+(t,Y_n,j) - μ(t) ),    t ∈^d.

Observe that V_n = Z_n if the distribution of X is N_d(0, I_d), since then the functions μ and m coincide.
We will show that, as n →∞, V_n converges in distribution to a centred Gaussian random element V of .
The only technical problem in proving such a result is the fact that V_n is based on the scaled residuals Y_n,1,…,Y_n,n and not
on X_1,…,X_n. If V_n^0(t) denotes the modification of V_n(t) that results from replacing Y_n,j with X_j, a
Hilbert space central limit theorem holds for V_n^0, since the summands comprising V_n^0(t) are i.i.d. square-integrable centred
random elements of . Writing  for convergence in distribution of random elements of  and random variables, the basic idea
to prove V_n  V is to find a random element V_n of , such that V_n  V and
V_n - V_n = o_(1) as n →∞. In what follows, the stochastic Landau symbol o_(1) refers to
convergence to zero in probability in , i.e., we have to show

    V_n - V_n^2_ = ∫( V_n(t) - V_n(t))^2 w_a(t)   dt  = o_(1)   as  n →∞.

To state the main result of this section, let ψ_X(t) =  [ exp( it^⊤ X)], t ∈^d, denote the characteristic function of X, and put

    ψ_X^+(t) :=  Re ψ_X(t)  +  Im ψ_X(t),   ψ_X^-(t) :=  Re ψ_X(t)  -  Im ψ_X(t),

where Rew and Imw stand for the real and the imaginary part of a complex number w, respectively.
For a twice continuously differentiable function f: ^d →, let Hf(t) denote the Hessian matrix
of f, evaluated at t. Furthermore, recall the gradient operator ∇ and the Laplace operator Δ from
Section <ref>.


Let

    V_n(t) := 1/√(n)∑_j=1^n v(t,X_j),    t ∈^d,

where

    v(t,x)     =     v_1(t,x) + v_2(t,x) + v_3(t,x) + v_4(t,x),
    
    v_1(t,x)     =    x^2  CS^+(t,x),  v_2(t,x) = 1/2 t^⊤(xx^⊤ -  I_d ) ∇Δψ_X^+(t),
    
    v_3(t,x)     =    (2 ∇ψ_X^-(t) + Δψ_X^-(t)t )^⊤ x,  v_4(t,x) = x^⊤ Hψ_X^+(t)x.

We then have (<ref>).


The proof of Proposition <ref> is given in Section <ref>.


Since (X) =0, (XX^⊤) =  I_d and [X^2  CS^+(t,X)] = - Δψ_X^+(t), we have (writing tr for trace)
v(t,X) =  - Δψ_X^+(t) +  tr (H ψ_X^+(t))  = 0.
Thus, v(·,X_1), …, v(·,X_n) are i.i.d. centred square-integrable random elements of , and the central limit theorem in
Hilbert spaces gives V_n  V for some centred  Gaussian element V of .  In view of (<ref>) and Slutsky's lemma,
we therefore can state the main result of this section.


Let X,X_1,X_2, … be i.i.d. random vectors satisfying X^4 < ∞, (X) =0 and (XX^⊤ ) =  I_d.
For the sequence of random elements V_n defined in (<ref>) we have

    V_n  V   as  n →∞,

where V is a centred Gaussian element of  having  covariance kernel

    L(s,t) =  [ v(s,X)v(t,X)],    s,t ∈^d,

where v(t,x) is given in (<ref>).











§ THE LIMIT NULL DISTRIBUTION OF T_N,A


In this section we derive the limit distribution of T_n,a under the null hypothesis  (<ref>). In view of affine invariance, we assume
that X has a d-variate standard normal distribution. Since the process Z_n(t) given in (<ref>) is nothing  but V_n, as  defined in (<ref>),
in this special case, we have the following result.


 Suppose that X has some non-degenerate normal distribution. Putting

    d_2 := d+2,    d_4 := d+4,

we have the following:


  a) There is a centred Gaussian random element Z of  with covariance kernel

    K(s,t)   =   ψ(s - t)((s - t^2 - d_2)^2  -  2d_2) +ψ(s)ψ(t){-(s^⊤ t)^2/2(s^2 - d_4)(t^2 - d_4)
           +2d_2(s^2 +t^2)  -s^4  -t^4 -s^2t^2  -  s^⊤ t (s^2 - d_2)(t^2 - d_2) - dd_2},

s,t ∈^d, such that, with Z_n defined in (<ref>), we have
Z_n  Z  in  as  n →∞.

  b)
We have

    T_n,a∫ Z^2(t)   w_a(t)   dt.




Notice that b) follows from a) and the continuous mapping theorem.
Part a) follows from Theorem <ref>. For the special case X ∼ N_d(0, I_d), we have
ψ_X^+(t) = ψ_X^-(t) = exp(-t^2/2) = ψ(t), which entails ∇ψ(t)  =  - t ψ(t),
Δψ(t)  =  (t^2-d ) ψ(t),
Hψ(t)  =  ( tt^⊤ -  I_d ) ψ(t),  and
∇Δψ(t)  =  t ψ(t) ( 2+d-t^2 ).
Thus, the function v(t,x) figuring in the statement of Proposition <ref> takes the special form

    h(x,t)     =    x^2  CS^+(t,x) - (2ψ(t) + m(t)) t^⊤ x - ψ(t)x^2 
           + (2 ψ(t) + m(t)/2) (t^⊤ x)^2 - ( ψ(t) + m(t)/2) t^2.

Long but straightforward computations, using symmetry arguments and the identities

    [X^2(s^⊤ X)^2]   =   (d+2)s^2,
    [X^4cos(t^⊤ X)]   =   ((d+2-t^2)(d-t^2)-2t^2)ψ(t),
    [X^2s^⊤ X sin(t^⊤ X)]   =   s^⊤ t(d+2-t^2)ψ(t),
    [X^2(s^⊤ X)^2 cos(t^⊤ X)]   =   (d+2-t^2)(s^2-(s^⊤ t)^2)ψ(t)-2(s^⊤ t)^2ψ(t),

s,t∈^d, show that the covariance kernel K(s,t) = [h(s,X)h(t,X)] takes the form given above.





Let T_∞,a be a random variable with the limit null distribution of T_n,a, i.e., with the distribution of
∫ Z^2(t) w_a(t)   dt. Since
(T_∞,a) = ∫ K(t,t) w_a(t)   dt,
the following result may be obtained by straightforward but tedious manipulations of integrals.


Putting

    c_j(a,d) := π^d/2 d/(a+1)^d/2 +j,    j =1,2,3,4,

we have

    (T_∞,a)     =     d(d+2)((π/a)^d/2 -  (π/a+1)^d/2)
           - c_4(a,d) (d+2)(d+4)(d+6)/32+c_3(a,d)(d+2)(d+3)(d+4)/8
           - c_2(a,d) (d+2)(d^2+4d+14)/8 + c_1(a,d) (d-2)(d+2)/2.



The quantiles of the distribution of T_∞,a can be approximated by a Monte Carlo method, see Section <ref>.








§ CONTIGUOUS ALTERNATIVES

In this section, we consider a triangular array (X_n1, …, X_nn), n ≥ d+1, of
rowwise i.i.d. random vectors with Lebesgue density

    f_n(x) = φ(x) ( 1 + g(x)/√(n)),    x ∈^d,

where φ(x) = (2π)^-d/2exp(-x^2/2), x ∈^d, is the density of the
standard normal distribution N_d(0, I_d), and g is some bounded measurable function  satisfying
∫ g(x) φ(x)   d x = 0.
We assume that n is sufficiently large to render g nonnegative. Recall Z_n(t) from (<ref>).

 Under the triangular X_n,1, …, X_n,n given above, we have

    Z_n  Z + c    as  n →∞,

where Z is the centred random element of  figuring in Theorem <ref>, and

    c(t) = ∫ h(x,t) g(x) φ(x)   dx,

with h(x,t) given in (<ref>).



Since the reasoning uses standard LeCam theory on contiguous probability measures and closely parallels  that given in Section 3 of <cit.>,
it will be omitted.


 Under the conditions of Theorem <ref>, we have

    T_n,a∫(Z(t) + c(t))^2 w_a(t)   d t.



From Theorem <ref> and the above corollary, we conclude that the test for
multivariate normality based on T_n,a is able to detect alternatives which converge to the normal distribution at
the rate n^-1/2, irrespective of the underlying dimension d. The test of Bowman and Foster (see <cit.>) is
a prominent example of an affine invariant tests for multivariate normality
that is consistent against each fixed non-normal alternative distribution but nevertheless lacks this property of n^-1/2-consistency,
see Section 7 of <cit.>.








§ FIXED ALTERNATIVES AND CONSISTENCY

In this section we assume that X,X_1,X_2, … are i.i.d. with a distribution that is absolutely continuous with respect to Lebesgue measure, and that
X^4 < ∞. In view of affine invariance, we make the additional assumptions (X) =0 and (XX^⊤) =  I_d. Recall m(t) from
(<ref>) and μ(t) from (<ref>).  Writing  for -almost sure convergence, our first result is a strong limit for T_n,a/n.


If X^4 < ∞, we have

    T_n,a/nΔ_a= ∫(μ(t) - m(t))^2 w_a(t)   dt    as  n →∞.



The proof of Theorem <ref> is given in Section <ref>.


Let ψ_X(t) =  [exp( it^⊤ X)], t ∈^d, denote the characteristic function of X. We have
Δψ_X(t) = - [X^2 exp( it^⊤ X)]. Since Δψ(t) = (t^2-d)ψ(t), some algebra gives

    Δ_a = ∫|Δψ_X(t) - Δψ(t) |^2 w_a(t)   d t.

By Theorem <ref>, we have Δ_a=0 if and only if X has the normal distribution N_d(0, I_d).
In view of Theorem  <ref>, T_n,a→∞ -almost  surely under any alternative distribution
satisfying X^4 < ∞. Since, according to Theorem <ref>,  the sequence of critical values of a level-α-test of H_0 that rejects H_0 for
large values of T_n,a stays bounded, we have the following result.


 The test for multivariate normality based on T_n,a is consistent against any fixed alternative
distribution satisfying 𝔼X^4 < ∞.




By analogy with  Theorem <ref>, the next result shows that  the (population) measure of multivariate skewness in the sense  of Móri, Rohatgi and Székely (see <cit.>)
emerges as the limit of Δ_a, after a suitable normalization, as a →∞.

  If X^6 < ∞, then

    lim_a →∞2 a^d/2+1/π^d/2Δ_a =  (X^2 X)^2.



The proof of Theorem <ref> is given in Section <ref>.



We now show that the limit distribution of √(n)(T_n,a/n- Δ_a) as n →∞ is centred normal.
This fact is essentially a consequence of Theorem 1 in <cit.>. The reasoning is as follows:
By (<ref>), we have T_n,a = Z_n^2_,
where Z_n is given in (<ref>). Putting z(t) := μ(t)-m(t), t ∈^d,  display (<ref>) shows that
Δ_a = z^2_.
Now, defining Z_n^*(t) := n^-1/2 Z_n(t), it follows that

    √(n)( T_n,a/n - Δ_a )     =    √(n)( Z_n^*^2_ - z^2_) = √(n)⟨ Z_n^* -z,Z_n^*+z ⟩_
        =    √(n)⟨ Z_n^* -z,2z + Z_n^*-z ⟩_
        =     2 ⟨√(n)(Z_n^* -z),z ⟩_ + 1/√(n)√(n)(Z_n^* - z)^2_.

A little thought shows that √(n)(Z_n^*(t) - z(t)) = V_n(t), where V_n is given in (<ref>).
By Theorem <ref>, V_n  V for a centred Gaussian element V of .
 As a consequence, the  second summand in (<ref>) is o_(1) as n →∞, and, by the continuous mapping theorem,
 the first summand
converges in distribution to 2⟨ V,z⟩_. The latter random variable has a centred normal distribution with variance σ_a^2 := 4  [⟨ V,z ⟩_ ^2].
The following theorem summarizes our findings.


For a fixed alternative distribution satisfying X^4 < ∞, (X) =0 and (XX^⊤) =  I_d, we have

    √(n)( T_n,a/n - Δ_a )  N(0,σ_a^2)     as  n →∞,

where

    σ_a^2 = 4 ∬ L(s,t) z(s) z(t)   w_a(s) w_a(t)   ds dt

and L(s,t) is given in (<ref>).



To complete the proof, notice that, by Fubini's theorem,

    σ_a^2       =      4  [⟨ V,z ⟩_ ^2]  =  4 [(∫ V(s)   z(s)   w_a(s)   d s)  (∫ V(t)   z(t)   w_a(t)   d t) ]
        =      4 ∬[V(s)V(t)]   z(s)   z(t)   w_a(s) w_a(t)   ds dt.



Under slightly stronger conditions on ℙ^X, there is a consistent estimator of σ_a^2.
To obtain such an estimator, we replace L(s,t)  as well as z(s) and z(t) figuring in (<ref>)
with suitable empirical counterparts. To this end, notice that, by (<ref>), we have

    L(s,t) = ∑_i=1^4 ∑_j=1^4 L^i,j(s,t),

where

    L^i,j(s,t) = [ v_i(s,X)v_j(t,X)]

and v_j(t,x), j ∈{1,2,3,4}, are given in (<ref>), (<ref>).
  Since ∇Δψ_X^±(t) =  ∓[ CS^∓(t,X) X^2X], ∇ψ_X^±(t) = ±[ CS^∓(t,X)X],
Δψ_X^±(t) = - [ CS^±(t,X)X^2] and Hψ_X^+(t) = - [ CS^+(t)(t,X) X X^⊤], parts a) – d) of the following lemma show that the unknown quantities ∇Δψ_X^+(t), ∇ψ_X^-(t),
Δψ_X^-(t) and Hψ_X^+(t) that figure in the expressions of v_2(t,x), v_3(t,x) and v_4(t,x) can be replaced with consistent
estimators that are based on the scaled residuals Y_n,1, …,Y_n,n defined in (<ref>).

 
If X_1^6 < ∞, we have




  a) Ψ_1,n(t) := n^-1∑_j=1^n  CS^+(t, Y_n,j)Y_n,j -∇ψ_X^-(t),

  b) Ψ_2,n(t) := n^-1∑_j=1^n  CS^+(t, Y_n,j)Y_n,jY_n,j^⊤ -Hψ_X^+(t),

  c) Ψ_3,n^±(t) := n^-1∑_j=1^n  CS^±(t, Y_n,j)Y_n,j^2  -Δψ_X^±(t),

  d) Ψ_4,n^±(t) := n^-1∑_j=1^n  CS^±(t, Y_n,j)Y_n,j^2Y_n,j±∇Δψ_X^∓(t),

  e) Ψ_5,n(t) := n^-1∑_j=1^n  CS^+(t, Y_n,j)Y_n,j^2Y_n,jY_n,j^⊤HΔψ_X^+(t).



The proof of Lemma <ref> is given in Section 9.



In view of Lemma <ref>, a suitable estimator of L(s,t) defined in (<ref>) is

    L_n(s,t) = ∑_i=1^4 ∑_j=1^4 L_n^i,j(s,t),

where

    L_n^i,j(s,t) = 1/n∑_k=1^n v_n,i(s,Y_n,k)v_n,j(t,Y_n,k),

and

    v_n,1(s,Y_n,k)     =    Y_n,k^2  CS^+(s,Y_n,k),    v_n,2(s,Y_n,k) = - 1/2 s^⊤ (Y_n,kY_n,k^⊤ -  I_d) Ψ_4,n^-(s),
    
    v_n,3(s,Y_n,k)     =     - (2Ψ_1,n(s) + Ψ_3,n^-(s)s)^⊤ Y_n,k,    v_n,4(s,Y_n,k) = -Y_n,k^⊤Ψ_2,n(s) Y_n,k.


By straightforward algebra we have

    L^1,1_n(s, t)    =   n^-1∑_j=1^n Y_n,j^4cos((t-s)^⊤Y_n,j) + n^-1∑_j=1^n Y_n,j^4sin((t+s)^⊤Y_n,j),
    
        L^1,2_n(s, t)    =   -1/2Ψ_4,n^-(s)^⊤Ψ_5,n(t)s + 1/2Ψ_3,n^+(t)Ψ_4,n^-(s)^⊤s,
    
        L^1,3_n(s, t)    =   (-2Ψ_1,n(s) - Ψ_3,n^-(s)s)^⊤Ψ_4,n^+(t)
    
        L^1,4_n(s, t)    =   -n^-1∑_j=1^n Y_n,j^2 CS^+(t, Y_n,j)Y_n,j^⊤Ψ_2,n(s)Y_n,j
    
        L^2,2_n(s, t)    =   1/4Ψ_4,n^-(t)^⊤n^-1∑_j=1^n Y_n,jY_n,j^⊤tΨ_4,n^-(s)^⊤(Y_n,jY_n,j^⊤ - I_d)s,
    
        L^2,3_n(s, t)    =   Ψ_4,n^-(t)^⊤n^-1∑_j=1^n Y_n,jY_n,j^⊤t(Ψ_1,n(s) + 1/2Ψ_3,n^-(s)s)^⊤Y_n,j,
    
        L^2,4_n(s, t)    =   1/2Ψ_4,n^-(t)^⊤n^-1∑_j=1^n (Y_n,jY_n,j^⊤ - I_d)tY_n,j^⊤Ψ_2,n(s)Y_n,j,
    
        L^3,3_n(s, t)    =   (2Ψ_1,n(t) + Ψ_3,n^-(t)t)^⊤(2Ψ_1,n(s) + Ψ_3,n^-(s)s),
    
        L^3,4_n(s, t)    =   n^-1∑_j=1^n (2Ψ_1,n(s) + Ψ_3,n^-(s)s)^⊤Y_n,jY_n,j^⊤Ψ_2,n(t)Y_n,j,
    
        L^4,4_n(s, t)    =   n^-1∑_j=1^n Y_n,j^⊤Ψ_2,n(t)Y_n,jY_n,j^⊤Ψ_2,n(s)Y_n,j.

Notice that, by symmetry, L_n^i,j = L_n^j,i if i ≠ j.
Since z(s) = μ(s)-m(s) = [X^2  CS^+(s,X)] - m(s), a natural estimator of z(s) is

    z_n(s) = 1/n∑_k=1^n  CS^+(s,Y_n,k) Y_n,k^2 -m(s).

Writing ℙ⟶ for convergence in probability, we then have the following result.


Suppose 𝔼X^6 < ∞, (X) =0 and (XX^⊤) =  I_d. Let

    σ_n,a^2 = 4 ∬ L_n(s,t)   z_n(s)   z_n(t)   w_a(s) w_a(t)   ds dt,

where L_n(s,t) and z_n(s) are defined in (<ref>) and (<ref>), respectively. We then have

    σ_n,a^2 ℙ⟶σ_a^2.



The proof of Theorem <ref> is given in Section <ref>.


Under the conditions of Theorem <ref>, Theorem <ref> and Sluzki's lemma yield

    √(n)/σ_n,a( T_n,a/n - Δ_a )  N(0,1)     as  n →∞,

provided that σ_a^2 >0. We thus obtain the following asymptotic confidence interval for Δ_a.


Let α∈ (0, 1), and write Φ_1-α/2 for the (1-α/2)-quantile of the standard normal law. Then

    I_n,a,α := [T_n,a/n - σ_n,a/√(n)Φ_1-α/2, T_n,a/n + σ_n,a/√(n)Φ_1-α/2]

is an asymptotic confidence interval for Δ_a.




We consider the case d=1, a=0.1 and the following standardized symmetric alternative distributions:
        the uniform distribution U(-√(3), √(3)), the Laplace distribution
         L(0, 1/√(2)), and the logistic distribution Lo(0, √(3)/π).
         The corresponding characteristic functions and their second derivatives are given by

    φ_ U(t)    = sin(√(3)t)/√(3)t,    φ_ U”(t)    = √(3)(2 - 3t^2)sin(√(3)t) - 6tcos(√(3)t)/3t^3, 
    φ_ L(t)    = 2/2 + t^2,    φ_ L”(t)    = 12t^2 - 8/(2 + t^2)^3, 
    φ_ Lo(t)    = √(3)t/sinh(√(3)t),    φ_ Lo”(t)    = 3/23√(3)t - 2sinh(2√(3)t) + √(3)tcosh(2√(3)t)/sinh(√(3)t)^3.

The pertaining values of Δ_0.1 are Δ_0.1, U≈0.3322, Δ_0.1, L≈0.127 and Δ_0.1, Lo≈0.033.
Table <ref> shows the percentages of coverage of the confidence intervals I_n,0.1,α
for α = 0.05 and several sample sizes. Each entry is based on 5000 Monte-Carlo-replications.
The results are quite satisfactory for n ≥ 50.




 A further application of (<ref>) addresses a fundamental drawback inherent in any goodness of fit test, see <cit.>.
If a level-α-test of H_0 does not reject H_0, the hypothesis H_0 is by no means validated or confirmed, since each model is wrong, and
there is probably only not enough evidence to reject H_0. However, if we define a certain "essential distance"
 δ_0 >0, we can consider the inverse testing problem

    H_δ_0: Δ_a(F) ≥δ_0   against  K_δ_0: Δ_a(F) < δ_0.

Here, the dependence of Δ_a on the underlying distribution has been made explicit.

From (<ref>), we obtain an asymptotic level-α-neighborhood-of-model-validation
 test of H_δ_0 against K_δ_0, which rejects H_δ_0 if and only if

    T_n.a/n ≤ δ_0 - σ_n,a/√(n)Φ^-1(1-α).

Indeed, from  (<ref>) we have for each F ∈ H_δ_0

    lim sup_n→∞_F (T_n.a/n≤δ_0 - σ_n,a/√(n)Φ^-1(1 -α))
         =     lim sup_n→∞_F (√(n)/σ_n,a(T_n,a/n  - δ_0) ≤ -Φ^-1(1 -α) )
       ≤   α .

Moreover, it follows that

    lim_n→∞_F (T_n,a/n≤δ_0 - σ_n,a/√(n)Φ^-1(1-α))  =  α

for each F such that Δ_a(F) = δ_0. Finally, we have

    lim_n→∞_F(T_n,a/n≤δ_0 - σ_n,a/√(n)Φ^-1(1-α) ) = 1

if Δ_a(F) < δ_0. Thus, the test is consistent against each alternative.




The double integral figuring (<ref>) may be calculated explicitly. To this end, notice that
σ_n,a^2 = ∑_i,j=1^4 σ_n,a^i,j,
  where

    σ_n,a^i,j = 4 ∬ L_n^i,j(s,t) z_n(s)z_n(t) w_a(s)w_a(t)   ds dt

and σ_n,a^i,j = σ_n,a^j,i, i,j ∈{1,…,4}, by symmetry.
We put

    q_1,a(y)    := ∫ m(t)  CS^+(t, y)w_a(t) dt = (2π)^d/2/(2a+1)^d/2+2(y^2 + 2da(2a+1))exp(-1/2y^2/2a+1), 
    
        p_1,a(y,z)    := ∫ CS^+(t, y) CS^+(t, z)w_a(t) dt = (π/a)^d/2exp(-y - z^2/4a), 
    
        p_2,a(y,z)    := ∫ CS^+(t, y)  CS^-(t,z)tw_a(t) dt = (π/a)^d/21/2aexp(-y - z^2/4a)(y-z), 
    
        q_2,a(y)    := ∫ m(t) CS^-(t, y)tw_a(t) dt 
       = (2π)^d/2/(2a+1)^d/2+3(2(2a+1)(1-ad) - y^2)exp(-1/2y^2/2a+1)y

for y,z ∈ℝ^d and

    P_n^1,a,1    :=    1/n^2∑_j,k=1^n Y_j^2 Y_k^2 p_1,a(Y_j, Y_k) - 1/n∑_j=1^n Y_j^2 q_1,a(Y_j), 
    P_n^1,a,1    :=    1/n^2∑_j,k=1^n Y_j^2Y_k^2p_1,a(Y_j,Y_k)Y_k - 1/n∑_j=1^n Y_j^2q_1,a(Y_j)Y_j, 
    P_n^1,a,2    :=    1/n^2∑_j,k=1^n Y_j^2p_1,a(Y_j, Y_k)Y_k - 1/n∑_j=1^n q_1,a(Y_j)Y_j, 
    P_n^1,a    :=    1/n^2∑_j,k=1^n Y_j^2Y_k^2p_1,a(Y_j, Y_k)Y_kY_k^⊤ - 1/n∑_j=1^n Y_j^2Y_jY_j^⊤q_1,a(Y_j), 
    
        P_n^1,a,2(Y_j)     :=    1/n^2∑_k,ℓ =1^n   Y_k^2(Y_j^⊤Y_ℓ)^2p_1,a(Y_k, Y_l) - 1/n∑_k=1^n (Y_j^⊤Y_k)^2q_1,a(Y_k), 
    
        P_n^1,a,3(Y_j)     :=    1/n∑_k=1^n Y_k^2p_1,a(Y_j, Y_k) - q_1,a(Y_j), 
    
        P_n^2,a,1    :=    1/n^2∑_j,k=1^n Y_j^2Y_k^2Y_k^⊤p_2,a(Y_j, Y_k) - 1/n∑_j=1^n Y_j^2Y_j^⊤q_2,a(Y_j), 
    
        P_n^2,a,2    :=    1/n^2∑_j,k=1^n Y_j^2Y_k^2Y_k^⊤P_1,a p_2,a(Y_j, Y_k) - 1/n∑_j=1^n Y_j^2Y_j^⊤P_1,aq_2,a(Y_j), 
    P_n^2,a    :=    1/n^2∑_j,k=1^n Y_j^2Y_k^2p_2,a(Y_j, Y_k) - 1/n∑_j=1^n Y_j^2q_2,a(Y_j), 
    
        P_n^2,a,3(Y_j)     :=    1/n^2∑_k, ℓ =1^n  Y_k^2Y_ℓ^2Y_k^⊤Y_jY_j^⊤p_2,a(Y_k, Y_ℓ) - 1/n∑_k=1^n Y_k^2Y_k^⊤Y_jY_j^⊤q_2,a(Y_k),

where Y_j is shorthand for Y_n,j.
Notice that P_n^1,a,1, P_n^1,a,2 and P_n^2,a are vectors and
 P_n^1,a is a matrix. By straightforward but tedious manipulations of the integrals in (<ref>), each σ_n^i,j
 is seen to be an  arithmetic mean of functions of the scaled residuals, namely:

    σ^1,1_n,a    =    4/n∑_j=1^n Y_j^4P_n^1,a,3(Y_j)^2,   σ^1,2_n,a = 2P_n^1,a,1P_n^2,a,1 - 2P_n^2,a,2, 
    σ^1,3_n,a    =     -4(2P_n^1,a,2 + P_n^2,a)^⊤P_n^1,a,1,   σ^1,4_n,a = -4/n∑_j=1^n Y_j^2P_n^1,a,2(Y_j)P_n^1,a,3(Y_j),    
    σ^2,2_n,a    =    1/n∑_j=1^n P_2,a,3(Y_j)^2 - (P_n^2,a,1)^2,   σ^2,3_n,a = 4/n∑_j=1^n P_2,a,3(Y_j)(P_n^1,a,2 + 1/2P_n^2,a)^⊤Y_j,

    
    σ^2,4_n,a    =    2/n∑_j=1^n (P_n^2,a,3(Y_j) - P_n^2,a,1)P_n^1,a,2(Y_j),   σ^3,3_n,a = 4 ·2P_n^1,a,2 + P_n^2,a^2, 
    σ^3,4_n,a    =    4/n∑_j=1^n P_n^1,a,2(Y_j)(2P_n^1,a,2 + P_n^2,a)^⊤Y_j,   σ^4,4_n,a = 4/n∑_j=1^n P_n^1,a,2(Y_j)^2.












§ SIMULATIONS

In this section, we present the results of a Monte Carlo simulation study on the finite-sample power performance of the test based on T_n,a with that of
 several competitors. Since different procedures have been used for univariate and multivariate data, we distinguish the cases d=1 and d ≥ 2.
 In the univariate case, the sample sizes are n∈{20,50,100}, and in the multivariate case
  we restrict the simulations to n∈{20,50}. The nominal level of significance is fixed throughout all simulations to 0.05.
  Critical values for T_n,a  (in fact, for a scaled version of T_n,a in order to obtain values of a similar magnitude across the range of values for d and a considered)
  have been simulated under H_0 with 100 000 replications, see Table <ref>. The critical values in the rows in Table <ref> denoted by "∞" represent approximations of quantiles of the distribution of the limit random element T_∞,a=∫ Z^2(t)   w_a(t)   dt in Theorem <ref> b). To obtain these values we simulate (say) m i.i.d. random supporting points U_1,…,U_m where U_1∼_d(0,(2a)^-1 I_d), which are adapted to the integration over ^d with respect to the weight function w_a(t)=exp(-at^2)  and a large number (say) ℓ of random variables Z_j=1/d^2mX_j^2, j=1,…,ℓ, with i.i.d. X_j∼_m(0,Σ_K), where the covariance matrix Σ_k is given by Σ_K=(K(U_k_1,U_k_2))_k_1,k_2∈{1,…,m} and K is the covariance kernel in Theorem <ref> a). Next, we calculate the empirical 95% quantile of Z_1,…,Z_ℓ. Each approximation was simulated with ℓ=100,000 and m=1,000 for d∈{2,3,5,10} and each entry in Tables <ref> - <ref>, which exhibit percentages of rejection of H_0 of the tests under consideration against various alternative distributions,  is based on 10 000 replications. Entries that are marked with ∗ denote 100%.
  The simulations have been performed using the statistical computing environment R, see <cit.>.





 §.§ Testing univariate normality

For testing the hypothesis H_0 that the distribution of X is univariate normal (i.e., ℙ^X ∈𝒩_1),
we considered the following competitors to the new test statistic:


  * the Shapiro–Wilk test (SW), see <cit.>,

  * the Shapiro–Francia test (SF), see <cit.>,

  * the Anderson–Darling test (AD), see <cit.>,

  * the Baringhaus–Henze–Epps–Pulley test (BHEP), see <cit.>,

  * the del Barrio–Cuesta-Albertos–Mátran–Rodríguez-Rodríguez test (BCMR), see <cit.>,

  * the Betsch–Ebner test (BE), see <cit.>.

For the implementation of the first three tests in R
we refer to the package nortest by <cit.>. Tests based on the empirical characteristic function
 are represented by the BHEP-test with tuning parameter a > 0. The statistic is given in (<ref>),
 a = 1 is fixed, and the critical values are taken from <cit.>.

The BCMR-test is based on the L^2-Wasserstein distance, see section 3.3 in <cit.>, and is given by

    BCMR =
    	n ( 1 - 1/S_n^2( ∑_k = 1^n X_(k)∫_k - 1/n^k/nΦ^-1(t)  dt )^2 ) - ∫_1/n + 1^n/n + 1t (1 - t)/( φ( Φ^-1(t) ) )^2 dt.

Here, X_(k) is the k-th order statistic of X_1, …, X_n, S_n^2 is the sample variance,
 and Φ^-1 is the quantile function of the standard normal law. Simulated critical values can be found in <cit.>, or in Table 4 of <cit.>.

The BE-test is based on a L^2-distance between the empirical zero-bias transformation and the empirical distribution.
Since the zero-bias transformation has a unique fixed point under normality, this distance is minimal under H_0. The statistic is given by

    BE_a     =    2/n∑_1 ≤ j < k ≤ n{( 1 - Φ( Y_(k)√(a)) ) ( (Y_(j)^2 - 1)(Y_(k)^2 - 1) + a Y_(j) Y_(k)) .
          . + a/√(2 π a) exp(- Y_(k)^22 a) ( - Y_(j)^2 Y_(k) + Y_(k) + Y_(j))
    	}
           + 1/n∑_j=1^n{( 1 - Φ(Y_j√(a)) ) ( Y_j^4 + (a - 2) Y_j^2 + 1 ) . 
          . + a/√(2 π a) exp(- Y_j^22 a) ( 2 Y_j - Y_j^3 ) },

where Y_1, …, Y_n is shorthand for the scaled residuals Y_n, 1, …, Y_n,n, Y_(1)≤…≤ Y_(n) are the order statistics
of Y_1,…,Y_n, and Φ stands for the distribution function of the standard normal law.
  The implementation employs a bootstrap procedure to find a data driven version of the tuning parameter a, see <cit.>.
  We used B=400 bootstrap replications and the same grid of tuning parameters as in <cit.>, p. 19.

The alternative distributions are chosen to fit the extensive power study of univariate normality tests by <cit.>,
 in order to ease the comparison with a wide variety of other existing tests. As representatives of symmetric distributions
 we simulate the Student t_ν-distribution with ν∈{3, 5, 10} degrees of freedom, as well as the uniform distribution U(-√(3), √(3)).
 The asymmetric distributions are represented by the χ^2_ν-distribution with ν∈{5, 15} degrees of freedom,
 the Beta distributions B(1, 4) and B(2, 5), the Gamma distributions Γ(1, 5) and Γ(5, 1), parametrized by their shape and rate parameter,
 the Gumbel distribution Gum(1, 2) with location parameter 1 and scale parameter 2, the lognormal distribution LN(0, 1)
 as well as the Weibull distribution W(1, 0.5) with scale parameter 1 and shape parameter 0.5.
 As representatives of bimodal distributions we take the mixture of normal distributions NMix(p, μ, σ^2), where the random variables are generated by

    (1 - p)   N(0, 1) + p   N(μ, σ^2),    p ∈ (0, 1),  μ∈,  σ > 0.


Table <ref> shows that the empirical power estimates of the new test T_n,a outperform the other strong
procedures for the symmetric t-distribution,  and they can compete for most of the other alternatives.
Interestingly, the power does not differ too much when varying the tuning parameter a,
although an effect is clearly visible, especially for the uniform distribution.
A data driven choice as in <cit.> might be of benefit also in connection with the new testing procedure.




 §.§ Testing multivariate normality

In this subsection we consider testing the hypothesis H_0 that the distribution of X is multivariate normal (i.e., belongs to 𝒩_d),
for the dimensions d ∈{2,3,5}. As competitors to the new test statistic we chose


  * the Henze–Visagie test (HV), see <cit.>,

  * the Henze–Jiménez-Gamero test (HJG), see <cit.>,

  * the Baringhaus–Henze–Epps–Pulley test (BHEP), see <cit.>.

The HV-test uses a weighted L^2-type statistic based on a characterization of the moment generating function
that employs a system of first-order partial differential equations. The statistic is defined by

    HV_γ= 1/n(π/γ)^d/2∑_j,k = 1^n exp( Y_n,j + Y_n,k^2/4 γ)(Y_n,j^⊤ Y_n,k+Y_n,j + Y_n,k^2(1/4 γ^2-1/2 γ)+d/2γ),

where γ>2. We followed the recommendation of the authors in <cit.> and fixed γ=5.
Since the limiting statistic HV_∞ for γ→∞ is a linear combination of sample skewness
in the sense of Mardia and that of Móri, Rohatgi and Székely, we also included HV_∞.

The HJG-test uses a weighted L^2-distance between the empirical moment generating function of the standardized sample
 and the moment generating function of the standard normal distribution. The test statistic is given by

    HJG_β= 1/nβ^d/2∑_j,k = 1^n exp( Y_n,j + Y_n,k^2/4 β) - 2/√(β - 1/2)∑_j = 1^n exp( Y_n,j^2/4 β - 2)+ n/(β - 1)^d/2,

with β > 0. In our simulation we fix β=1.5. For each of the tests based on HV_5, HV_∞ and HJG_1.5,
 critical values were simulated with 100 000 replications.

Finally, the now classical BHEP-test examines the weighted L^2-distance between the empirical characteristic function
 of the standardized data and the characteristic function of the d-variate standard normal distribution. The statistic has the simple form

    BHEP_a   =   1/n^2∑_j,k = 1^n exp( -a^2/2Y_n,j - Y_n,k^2 )
           - 2(1+a^2)^-d/21/n∑_j = 1^n exp( -a^2Y_n,j^2/2(1+a^2))+ (1+2a^2))^-d/2,

with a tuning parameter a>0. A variety of values of a, i.e., a∈{0.1,0.25,0.5,0.75,1,2,3,5,10}, has been considered.
 Critical values can be found in Tables I - III of <cit.>, whereas missing critical values have been simulated separately with 100 000 replications.

The alternative distributions are chosen to fit the simulation study in <cit.> and are defined as follows. We denote by NMix(p,μ,Σ) the normal mixture distributions generated by

    (1 - p)   N_d(0,  I_d) + p   N_d(μ, Σ),    p ∈ (0, 1),  μ∈^d,  Σ > 0,

where Σ > 0 stands for a positive definite matrix. We write in the notation of above μ=3 for a d-variate
vector of 3's and Σ= B_d for a (d × d)-matrix containing 1's on the main diagonal and 0.9's on each off-diagonal entry.
 We denote by t_ν(0, I_d) the multivariate t-distribution with ν degrees of freedom, see <cit.>.
 By DIST^d(ϑ) we denote the d-variate random vector generated by independently simulated components of the distribution
 DIST with parameter vector (ϑ), where DIST is taken to be the Cauchy distribution C, the logistic distribution L,
  the Gamma distribution Γ as well as the Pearson Type VII  distribution P_VII, with ϑ denoting in this specific
   case the degrees of freedom. The spherical symmetric distributions where simulated using the R package distrEllipse,
  see <cit.>, and are denoted by 𝒮^d(), where DIST stands for the distribution of the radii
  and was chosen to be the exponential, the beta and the χ^2-distribution.

From Tables <ref> to <ref>, it is obvious that T_n,a outperforms the competing tests for most
of the alternatives considered, again showing that the tuning parameter has –compared to the BHEP test– little effect on the power.
 As in the univariate case T_n,a has very strong power against the multivariate t-distribution.
 If the radial distribution of the spherical symmetric alternatives has compact support, the BHEP test exhibits
  a better performance than T_n,a. The HV_5- and the HJG_1.5-test have a good power,  but they are
   mostly dominated by the BHEP-test and the T_n,a-test. Again, a data driven choice of the tuning parameter a
   would be beneficial for the test, but to the best of our knowledge no reliable multivariate method is yet available.
   We suggest to choose a small tuning parameter like a=0.25.










§ REAL DATA EXAMPLE: THE IRIS DATA SET

In 1936 R.A. Fisher presented the now classical data set called Iris Flower, see Table I in <cit.>.
The data consist of the four variables sepal, length, sepal width, petal length, and petal width, measured
on fifty specimens of each of three types or iris, namely Iris setosa, Iris versicolor, and Iris virginica.
This data set is included in the statistical language R, and it can be downloaded from the UCI Machine Learning Repository,
see <cit.>. That reference provides a list of articles that use this specific data set to validate clustering methods,
neural networks or learning algorithms, and it presents a typical test case for statistical classification techniques in machine learning,
such as support vector machines. A visualization of two-dimensional projections of the data set is given in Figure <ref>.




In Table <ref> we present empirical p-values simulated with 10 000 replications.
As can be seen, the test does not reject the hypothesis of normality on a small significance
level (like α=0.01) for the different species for each of the tuning parameters considered.
For the Iris setosa data, however, an increase of the significance level to 0.05
results in a rejection of the hypothesis for a=2 and a=3. For the whole data set,
 we observe a small p-value due to the mixture of the three populations and consequently reject the hypothesis H_0.



§ CONCLUDING REMARKS

We proved consistency of the test for multivariate normality based on T_n,a against each alternative distribution that
 satisfies the moment condition X^4 < ∞. Intuitively, the  test should be "all the more consistent" if X^4 = ∞.
 In fact, we conjecture consistency of the new test against any non-normal alternative distribution.

The limiting random element T_∞,a=Z^2_ℍ from Theorem <ref> b)
has the same distribution as ∑_j=1^∞λ_jN_j^2, where the N_j are i.i.d. standard normal random variables,
 and the λ_j are the positive eigenvalues corresponding to eigenfunctions of the linear integral operator
 Kf(s)=∫_^dK(s,t)f(t)  w_a(t)   dt associated with the covariance kernel K from Theorem <ref> a), i.e.,
  we have λ f(s)=∫_^dK(s,t)f(t)  w_a(t)   dt. These positive eigenvalues clearly depend on K
  and the weight function w_a. It is hardly possible to obtain a closed form expression for λ_1,λ_2,…. It would be interesting to approximate the eigenvalues either numerically or by some Monte Carlo method, since the largest eigenvalue has a crucial influence on the approximate Bahadur efficiency, see <cit.> and the monograph <cit.>, as well as <cit.> for results on distribution-free L_p-type statistics.



§ ACKNOWLEDGEMENT

The authors thank Ioannis Anapolitanos for his help in proving Theorem <ref> and for sharing his knowledge of Schrödinger operators.











§ PROOFS AND AUXILIARY RESULTS

This section contains the proofs of some of the theorems as well as some auxiliary results that are used in the main text.
Recall our standing assumption that the distribution of X is absolutely continuous. In what follows, let

    Δ_n,j = (S_n^-1/2-I_d)X_j - S_n^-1/2X_n,    j=1,…,n.

Notice that Δ_n,j = Y_n,j-X_j, where the scaled residuals Y_n,j are given in (<ref>).



Let X, X_1,X_2, … be i.i.d. random vectors satisfying X^4 < ∞,  (X) =0 and (X X^⊤) =  I_d.
We then have:


  a)∑_j=1^n Δ_n,j^2 = O_(1).

  b) 1/n∑_j=1^n Δ_n,j^2  0.

  c) max_j=1,…,nΔ_n,j = o_(n^-1/4).



 a) Notice that

    Δ_n,j^2 = X_j^⊤(S_n^-1/2-I_d)^2X_j - 2 X_n^⊤ S_n^-1/2 (S_n^-1/2-I_d)X_j + X_nS_n^-1X_n

and thus, using tr(AB) =  tr(BA), where tr(C) denotes the trace of a square matrix C,

    ∑_j=1^n Δ_n,j^2  =  tr((S_n^-1/2 -I_d)^2 ∑_j=1^n X_jX_j^⊤) - 2nX_n^⊤ S_n^-1/2(S_n^-1/2 -I_d)X_n
    + nX_n S_n^-1X_n .

Due to X_1^4 < ∞, the central limit theorem  implies

    √(n)(S_n^-1 -  I_d) = - S_n^-1√(n)( S_n -  I_d) = O_(1).

Since √(n)(S_n^-1/2-  I_d)(S_n^-1/2+  I_d) = √(n)(S_n^-1 -  I_d), it follows that √(n)(S_n^-1/2-  I_d) = O_(1).
In view of ∑_j=1^n X_jX_j^⊤ = O_(n) and √(n)X_n = O_(1), we are done.

b) After dividing both sides of (<ref>) by n, the first summand on the right hand side converges to 0 -almost surely because
n^-1∑_j=1^n X_jX_j^⊤ I_d,  and S_n^-1/2 I_d. The same holds for the other two terms,
since X_n  0.


c) Let A_ sp be the spectral norm of a matrix A. Then

    Δ_n,j≤S_n^-1/2 -  I_d_ sp X_j + S_n^-1/2_ spX_n

and thus

    max_j=1,…,nΔ_n,j≤S_n^-1/2 -  I_d_ sp max_j=1,…,nX_j + S_n^-1/2_ spX_n.

From Theorem 5.2 of <cit.> we have max_j=1,…,nX_j/n^1/4 0.
Since √(n)S_n^-1/2 -  I_d_ sp = O_(1), S_n^-1/2_ sp = O_(1) and √(n)X_n = O_(1), we are done.
 



Let X, X_1,X_2, … be i.i.d. random vectors satisfying X^6 < ∞,  (X) =0 and (X X^⊤) =  I_d.
We then have

    1/n∑_j=1^n Δ_n,j^k X_j^ℓ 0

as n →∞ for each k ∈ℕ and ℓ∈ℕ_0 such that k+ ℓ≤ 6.



From (<ref>) we have

    Δ_n,j^k  X_j^ℓ≤ 2^k S_n^-1/2 -  I_d_ sp^k  X_j^k+ℓ + 2^k S_n^-1/2_ sp^k X_n^k X_j^ℓ.

Since S_n^-1/2 -  I_d_ sp 0, S_n^-1/2 I_d and  X_n 0, the assertion follows from
the strong law of large numbers.




Proof of Proposition <ref>

In what follows, we put Y_j = Y_n,j and Δ_j = Δ_n,j for the sake of brevity.
Notice that
cos(t^⊤ Y_j) = cos(t^⊤ X_j) - sin(Θ_j) t^⊤Δ_j, sin(t^⊤ Y_j) = sin(t^⊤ X_j) + cos(Γ_j)t^⊤Δ_j,
where Θ_j, Γ_j depend on X_1,…,X_n and t and satisfy

    |Θ_j - t^⊤ X_j| ≤ |t^⊤Δ_j|,      |Γ_j - t^⊤ X_j| ≤ |t^⊤Δ_j|.

Since Y_j^2 = X_j^2 + Δ_j^2 + 2 X_j^⊤Δ_j, it follows that
V_n(t) = ∑_k=1^6 V_n,k(t),  where

    V_n,1(t)     =    1/√(n)∑_j=1^n {X_j^2  CS^+(t,X_j) - μ(t) },
    
    V_n,2(t)     =    1/√(n)∑_j=1^n X_j^2 t^⊤Δ_j (cos(Γ_j) - sin(Θ_j)),
    
    V_n,3(t)     =    2/√(n)∑_j=1^n X_j^⊤Δ_j  CS^+(t,X_j), 
    
    V_n,4(t)     =    2/√(n)∑_j=1^n X_j^⊤Δ_j t^⊤Δ_j (cos(Γ_j) - sin(Θ_j)),
    
    V_n,5(t)     =    1/√(n)∑_j=1^n Δ_j^2  CS^+(t,X_j),
    
    V_n,6(t)     =    1/√(n)∑_j=1^n Δ_j^2 t^⊤Δ_j (cos(Γ_j) - sin(Θ_j)).

We first show that V_n,ℓ = o_(1) if ℓ∈{4,5,6}. As for V_n,4, notice that, by the Cauchy–Schwarz inequality,
|V_n,4(t)| ≤ 4t n^-1/2max_i=1,…,nX_i∑_j=1^n Δ_j^2.
Since X^4 < ∞, Theorem 5.2 of <cit.> yields max_i=1,…,nX_i = o_(n^1/4). In view of
Proposition <ref> a), we have V_n,4 = o_(1).  The same proposition immediately also gives
V_n,5 = o_(1). Since
|V_n,6(t)| ≤ 2t n^-1/2max_i=1,…,nΔ_i∑_j=1^n Δ_j^2,
we have V_n,6 = o_(1) in view of Proposition <ref> a) and Proposition <ref> c).

We now consider V_n,2(t). Since |V_n,2(t)| ≤ 2t n^-1∑_j=1^n X_j^2 n^1/4max_i=1,…,nΔ_i n^^1/4,
the law of large numbers and Proposition  <ref> c) show that V_n,2 = o_(n^1/4). In view of
 (<ref>) and Proposition <ref> c), the error is thus o_(1) if we replace both  Γ_j and Θ_j
 with t^⊤ X_j. Moreover, plugging Δ_j = (S_n^-1/2- I_d)X_j - S_n^-1/2X_n into the definition of V_n,2(t),
 the error is o_(1) if we replace S_n^-1/2X_n with X_n. Recalling (<ref>), we thus obtain
 
    V_n,2(t) = 1/√(n)∑_j=1^n X_j^2 t^⊤{(S_n^-1/2 -  I_d) X_j - X_n } CS^-(t,X_j) + o_(1).

We now use display (2.13) of <cit.>, according to which

    2 √(n)(S_n^-1/2-  I_d) = - 1/√(n)∑_j=1^n (X_jX_j^⊤ -  I_d ) + O_(n^-1/2).

Plugging this expression into (<ref>) we obtain

    V_n,2(t)     =     - 1/2·1/n∑_k=1^n X_k^2  CS^-(t,X_k) t^⊤1/√(n)∑_j=1^n (X_jX_j^⊤ -  I_d) X_k 
           - 1/n∑_k=1^n X_k^2  CS^-(t,X_k)  1/√(n)∑_j=1^n t^⊤ X_j + o_(1).

In (<ref>) we now use the fact that tr(AB) =  tr(BA), where tr denotes trace and
 AB is a square matrix. Furthermore, the error is o_(1) if we replace n^-1∑_j=1^n X_j^2 CS^-(t,X_j) and n^-1∑_j=1^n X_j^2 X_j  CS^-(t,X_j) with
 their almost sure limits
 [X^2  CS^-(t,X)] = - Δψ_X^-(t)
 and
 [X^2 X  CS^-(t,X)] = - ∇Δψ_X^+(t),
  respectively. We therefore obtain
 
    V_n,2(t) = 1/√(n)∑_j=1^n {1/2∇Δψ_X^+(t) (X_jX_j^⊤ -  I_d ) t + Δψ_X^-(t) t^⊤ X_j } + o_(1).

In the same way, we proceed  with V_n,3(t) and, using [X  CS^+(t,X) X^⊤] = - Hψ_X^+(t) as well as [X CS^+(t,X)] = - ∇ψ_X^-(t),
finally arrive at

    V_n,3(t) = 1/√(n)∑_j=1^n { 2 X_j^⊤∇ψ_X^-(t) + X_j^⊤ H ψ_X^-(t) X_j + μ(t)} + o_(1).

By adding (<ref>), (<ref>) and (<ref>), we have V_n = V_n + o_(1),
where V_n is given in (<ref>).□








Proof of Theorem <ref>

By analogy with Z_n(t), as  defined in (<ref>), let

    Z^0_n(t) := 1/√(n)∑_j=1^n { CS^+(t,X_j) - m(t) },    t ∈^d.

A straightforward calculation gives n^-1/2(Z_n(t) - Z^0_n(t)) = ∑_j=1^3 U_n,j(t),
where

    U_n,1(t)     =    1/n∑_j=1^n X_j^2 (  CS^+(t,Y_n,j) -  CS^+(t,X_j)  ),
    
    U_n,2(t)     =    2/n∑_j=1^n X_j^⊤Δ_n,j CS^+(t,Y_n,j),   
    U_n,3(t)  =  1/n∑_j=1^n Δ_n,j^2  CS^+(t,Y_n,j).

Since |cos(t^⊤ Y_n,j) - cos(t^⊤ X_j)| ≤t Δ_n,j, |sin(t^⊤ Y_n,j) - sin(t^⊤ X_j)| ≤t Δ_n,j,
the Cauchy–Schwarz inequality gives |U_n,1(t)| ≤ 2 t( n^-1∑_j=1^n X_j^4)^1/2(n^-1∑_j=1^n Δ_n,j^2 )^1/2.



By the strong law of large numbers, we have n^-1∑_j=1^n X_j^4 X^4. In view of Proposition
<ref> b), we thus obtain U_n,1^2_ 0. Next, notice that, again by the Cauchy–Schwarz inequality,
|U_n,2(t)| ≤ 4 ( n^-1∑_j=1^n X_j^2)^1/2(n^-1∑_j=1^n Δ_n,j^2 )^1/2.



Hence, we have U_n,2^2_ 0. Finally, |U_n,3(t)|  ≤ 2 n^-1∑_j=1^n Δ_n,j^2 which, in view of Proposition
<ref> b), shows that U_n,3^2_ 0. Summarizing, we have

    n^-1/2(Z_n(· ) - Z^0_n(· ))_ 0.

By the strong law of large numbers in Banach spaces, it follows that n^-1/2 Z^0_n(·) μ(·) - m(·) as n →∞
in . In view of (<ref>), we thus obtain

    T_n,a/n =  n^-1/2 Z_n(· )^2_μ(·) - m(·)^2_ = ∫(μ(t)-m(t))^2 w_a(t)   d t.

□


Proof of Theorem <ref>

Starting with (<ref>), we have

    Δ_a     =    ∫(  (X^2  CS^+(t,X)) )^2 exp(-at^2)   d t 
           - 2 ∫ (d-t)^2  (X^2  CS^+(t,X)) exp( - 2a+1/2t^2 )  d t 
           + ∫ (d-t^2)^2 exp(-(a+1)t^2)  d t
        =:     I_1,a -  I_2,a + I_3,a,

say. Letting X_1,X_2 be independent copies of X, Fubini's theorem, the addition theorems for the sine function and the cosine function,
considerations of symmetry and (<ref>) yield

    I_1,a = ( π/a)^d/2[X_1^2 X_2^2 exp( - X_1-X_2^2/4a)].

From (<ref>) and (<ref>), we have

    I_2,a    =    2 (2π)^d/2/(2a+1)^d/2+2 [X^2 ( X^2 + 2da(2a+1)) exp(- X^2/2(2a+1)) ],
    
    I_3,a    =    π^d/2/(a+1)^d/2+2( a(a+1)d^2 + d(d+2)/4)

and thus

    2a (a/π)^d/2Δ_a     =      2a  [X_1^2 X_2^2 exp( - X_1-X_2^2/4a)]
           -  2 (2a)^d/2+1/(2a+1)^d/2+2[X^2 ( X^2 + 2da(2a+1)) exp(- X^2/2(2a+1)) ]
           + 2 a^d/2+1/(a+1)^d/2+2( a(a+1)d^2 + d(d+2)/4)
        =:     J_1,a - J_2,a + J_3,a,

say. An expansion of the exponential terms, dominated convergence (notice that exp(-u) ≤ 1-u + u^2 if u ≥ 0) and a binomial expansion gives

    J_1,a    =     2ad^2 - d X_1^4 + ( X_1^2 X_2^2 X_1^⊤ X_2 ) + O(a^-1),
    
    J_2,a    =      4 a d^2 - d^3 - 2 d^2 - d X_1^4 + O(a^-1),
    
    J_3,a    =     2ad^2 - d^3 - 2d^2 + O(a^-1)

and thus

    lim_a →∞ 2a (a/π)^d/2Δ_a     =    (X_1^2 X_2^2 X_1^⊤ X_2 ) =
    (X_1^2 X_1)^⊤(X_2^2 X_2)
        =     (X^2 X)^2.

Notice that the condition 𝔼X^6 < ∞ is not only needed for the existence of the final limit, but it also
occurs when dealing with J_2,a.□




Proof of Lemma <ref>

Putting Y_j = Y_n,j and Δ_j = Y_j - X_j, we have

    cos(t^⊤Y_j)    = cos(t^⊤X_j) - t^⊤Δ_jsin(t^⊤X_j) + ε_j(t), 
    sin(t^⊤Y_j)    = sin(t^⊤X_j) + t^⊤Δ_jcos(t^⊤X_j) + η_j(t),

where |ε_j(t)|, |η_j(t)| ≤t^2Δ_j^2 (see <cit.>, p. 8) and thus

    CS^±(t,Y_j) = CS^±(t, X_j) ± t^⊤Δ_j CS^∓(t, X_j) + ε_j(t) + η_j(t).

To prove a), notice that Ψ_1,n(t) = n^-1∑_j=1^n CS^+(t, X_j)X_j + R_1,n(t), where

    R_1,n(t)    = 1/n∑_j=1^n (t^⊤Δ_jCS^-(t, X_j) + ε_j(t) + η_j(t))X_j 
       + 1/n∑_j=1^n (CS^+(t, X_j) + t^⊤Δ_jCS^-(t, X_j) + ε_j(t) + η_j(t))Δ_j.

Since |CS^±(t, X_j)| ≤ 2 and |ε_j(t) + η_j(t)| ≤ 2 t^2 Δ_j^2, the Cauchy–Schwarz inequality yields

    |R_1,n(t)|    ≤2t/n∑_j=1^n Δ_jX_j + 2t^2/n∑_j=1^n Δ_j^2X_j + 2/n∑_j=1^n Δ_j
       + 2t/n∑_j=1^n Δ_j^2 + 2t^2/n∑_j=1^n Δ_j^3.

In view of Proposition <ref>, each summand converges to zero almost surely, which proves a).
The remaining assertions b), …, e) are treated similarly. To tackle b) and e), one can show negligibility of terms
by using the fact that the spectral norm ·_ sp of a matrix satisfies xy^⊤_ sp = x y, if x,y are
  (column) vectors in ℝ^d. The condition X_1^6 < ∞ is needed for part e), since

    Y_j^2Y_jY_j^⊤_sp   ≤(X_j^2 + 2X_jΔ_j + Δ_j^2)^2 
       = X_j^4 + 4X_j^3Δ_j + 6X_j^2Δ_j^2 + 4X_jΔ_j^3 + Δ_j^4,

and multiplication with Δ_j^2 (which arises from an expansion of CS^+(t, Y_n,j)) gives
 monomials of order 6.□



Proof of Theorem <ref>

The proof is similar to the proof of Theorem 5 of <cit.> and will therefore only be sketched.
From (<ref>),  (<ref>) and  (<ref>), we have σ_a^2 = ∑_i,j=1^4 σ_a^i,j, where

    σ_a^i,j = 4 ∬ L^i,j(s,t)z(s)z(t) w_a(s)w_a(t)   ds  dt,

and L^i,j(s,t) is given in (<ref>). It thus suffices to prove
σ_n,a^i,jℙ⟶σ_a^i,j for each pair (i,j), where
σ_n,a^i,j is given in (<ref>). The first step of the proof is to replace L_n^i,j(s,t) with
L_n,0^i,j(s,t), which arises from L_n^i,j(s,t) given in (<ref>) by throughout replacing Y_n,k with X_k in the functions v_n,j(s,Y_n,k), j ∈{1,…,4}.
Notice that this replacement also refers to the quantities Ψ^-_4,n(s), Ψ_1,n(s),  Ψ^-_3,n(s) and Ψ_2,n(s) that figure in the
definition of v_n,2, v_n,3 and v_n,4. Moreover, we replace z_n(s) with z_n,0(s) = n^-1∑_j=1^n  CS^+(s,X_j)X_j^2 - m(s).
Putting

    σ_n,0,a^i,j = 4 ∬ L_n,0^i,j(s,t)z_n,0(s)z_n,0(t) w_a(s)w_a(t)   ds  dt,

Fubini's theorem shows that σ_n,0,a^i,jℙ⟶σ_a^i,j. It thus remains to prove

    σ_n,a^i,j - σ_n,0,a^i,j  = o_ℙ(1)

as n →∞. To tackle σ_n,a^i,j - σ_n,0,a^i,j, we put

    Λ_n(s,t) = L_n^i,j(s,t)z_n(s)z_n(t) - L_n,0^i,j(s,t)z_n,0(s)z_n,0(t)

and notice that

    Λ_n(s,t)     =    (L_n^i,j(s,t) - L_n,0^i,j(s,t)) z_n(s)z_n(t)
             + L_n,0^i,j(s,t)(z_n(s)z_n(t)- z_n,0(s)z_n,0(t) ), 
    
     z_n(s)z_n(t)-z_n,0(s) z_n,0(t)     =    (z_n(s)-z_n,0(s)) (z_n(t)-z_n,0(t))
             + z_n,0(s) (z_n(t)-z_n,0(t)) + z_n,0(t)(z_n(s)-z_n,0(s)).

From (<ref>), we have |z_n(s)| ≤ 2 d+ m(s). Moreover, |z_n,0(s)| ≤ 2n^-1∑_j=1^n X_j^2 + m(s).
A Taylor expansion shows that |z_n(s) - z_n,0(s)| is bounded from above by a finite sum of terms of the type
s^ℓ n^-1∑_j=1^n X_j^βΔ_n,j^γ, where s ∈{0,1}, γ≥ 1 and β + γ≤ 3.
Next, by straightforward calculations we obtain that |L_n,0^i,j(s,t)| is bounded from above by s^ℓt^m times finitely many products
of the type n^-1∑_j=1^n X_j^β, where ℓ, m ∈{0,1} and β∈{1,2,3,4}.
It now follows from Proposition <ref> that the term figuring in (<ref>), multiplied with w_a(s)w_a(t) and integrated over ℝ^d ×ℝ^d,
converges to zero in probability.

As for the term figuring in (<ref>), we have |z_n(s)z_n(t)| ≤ (2d+m(s))(2d+m(t)). To tackle
L_n^i,j(s,t) - L_n,0^i,j(s,t), we confine ourselves to the case
i=1, j=2, since the other cases can be treated in a similar way. From (<ref>), we  have

    L^1,2_n(s, t)     =      -1/2{1/n∑_j=1^n  CS^-(s,Y_j)Y_j^2Y_j}^⊤{1/n∑_j=1^n  CS^+(t, Y_j)Y_j^2Y_jY_j^⊤ s}
            + 1/2{1/n∑_j=1^n  CS^+(t, Y_j)Y_j^2}{1/n∑_j=1^n  CS^-(s,Y_j)Y_j^2Y_j^⊤ s}.

Moreover,

    L_n,0^1,2(s,t)     =     -1/2{1/n∑_j=1^n  CS^-(s, X_j)X_j^2X_j}^⊤{1/n∑_j=1^n  CS^+(t, X_j)X_j^2X_jX_j^⊤ s}
            + 1/2{1/n∑_j=1^n  CS^+(t, X_j)X_j^2}{1/n∑_j=1^n  CS^-(s,X_j)X_j^2X_j^⊤ s}.

  Using (<ref>), some calculations show that, for each i ∈{1,2}, the norm of the difference of the i^ th curly bracket in (<ref>) and the   i^ th curly bracket in (<ref>)
  is bounded from above by  finite sum of terms of the type
  
    s^ℓt^m  1/n∑_j=1^n X_j^βΔ_n,j^γ,

   where ℓ, m ∈{0,1,2},
  γ≥ 1, and β + γ≤ 5. Since the same holds for the norm of the difference of the i^ th curly bracket in (<ref>) and the   i^ th curly bracket in (<ref>), i ∈{1,2}, it follows that |L^1,2_n(s, t) - L_n,0^1,2(s,t)| is bounded from above by a finite sum of terms, which are either products of two terms of type (<ref>), or a product of a term of type
  (<ref>) and one of the terms inside one of the curly brackets in (<ref>) or (<ref>). From Proposition <ref>, it follows that
  the term figuring in (<ref>), multiplied with w_a(s)w_a(t) and integrated over ℝ^d ×ℝ^d,
converges to zero in probability.□












unsrt





10

AS:2015
J. S. Allison and L. Santana.
On a data-dependent choice of the tuning parameter appearing in
  certain goodness-of-fit tests.
Journal of Statistical Computation and Simulation,
  85(16):3276–3288, 2015.

AE:2009
J. A. V. Alva and E. G. Estrada.
A generalization of Shapiro–Wilk's test for multivariate
  normality.
Communications in Statistics - Theory and Methods,
  38(11):1870–1883, 2009.

AD:1952
T. W. Anderson and D. A. Darling.
Asymptotic theory of certain "goodness of fit" criteria based on
  stochastic processes.
Annals of Mathematical Statistics, 23:193–212, 1952.

A:2007
M. A. Arcones.
Two tests for multivariate normality based on the characteristic
  function.
Mathematical Methods of Statistics, 16(3):177–201, 2007.

B:1960
R. R. Bahadur.
Stochastic comparison of tests.
The Annals of Mathematical Statistics, 31(2):276–295, 06 1960.

BEH:2017
L. Baringhaus, B. Ebner, and N. Henze.
The limit distribution of weighted L^2-Goodness-of-Fit
  statistics under fixed alternatives, with applications.
The Annals of the Institute of Statistical Mathematics,
  69(5):969–995, 2017.

BH:1988
L. Baringhaus and N. Henze.
A consistent test for multivariate normality based on the empirical
  characteristic function.
Metrika, 35(1):339–348, 1988.

BH:1991
L. Baringhaus and N. Henze.
Limit distributions for measures of multivariate skewness and
  kurtosis based on projections.
Journal of Multivariate Analysis, 38(1):51–69, 1991.

BH:1992
L. Baringhaus and N. Henze.
Limit distributions for Mardia's measure of multivariate skewness.
The Annals of Statistics, 20(4):1889–1902, 1992.

BN:1963
O. Barndorff-Nielsen.
On the limit behavior of extreme order statistics.
The Annals of Mathematical Statistics, 34(3):992–1002, 1963.

BE:2019
S. Betsch and B. Ebner.
Testing normality via a distributional fixed point property in the
  stein characterization.
TEST, Feb 2019.

BF:1993
A. Bowman and P. Foster.
Adaptive smoothing and density-based tests of multivariate normality.
Journal of the American Statistical Association, 88:529–537,
  1993.

CS:1978
D. Cox and N. Small.
Testing multivariate normality.
Biometrika, 65:263–272, 1978.

delBarrio2000
E. del Barrio, J. A. Cuesta-Albertos, C. Matrán, S. Csörgö, C. M.
  Cuadras, T. de Wet, E. Giné, R. Lockhart, A. Munk, and W. Stute.
Contributions of empirical and quantile processes to the asymptotic
  theory of goodness-of-fit tests.
TEST, 9(1):1–96, 2000.

DBCMRR:1999
E. del Barrio, J. A. Cuesta-Albertos, C. Matran, and J. M. Rodriguez-Rodriguez.
Tests of goodness of fit based on the L_2-Wasserstein distance.
The Annals of Statistics, 27(4):1230–1239, 1999.

DH:2008
J. A. Doornik and H. Hansen.
An omnibus test for univariate and multivariate normality.
Oxford Bulletin of Economics and Statistics, 70:927–939, 2008.

DG:2019
D. Dua and C. Graff.
UCI machine learning repository, 2017.

EP:1973
M. L. Eaton and M. D. Perlman.
The non-singularity of generalized sample covariance matrices.
The Annals of Statistics, 1(4):710–717, 1973.

E:2012
B. Ebner.
Asymptotic theory for the test for multivariate normality by Cox
  and Small.
Journal of Multivariate Analysis, 111:368–379, 2012.

FSN:2007
P. J. Farrell, M. Salibian-Barrera, and K. Naczk.
On tests for multivariate normality and associated simulation
  studies.
Journal of Statistical Computation and Simulation,
  77(12):1065–1080, 2007.

F:1936
R. A. Fisher.
The use of multiple measurements in taxonomic problems.
Annals of Eugenics, 7(2):179–188, 1936.

AG:2009
A. Genz and F. Bretz.
Computation of Multivariate Normal and t Probabilities.
Lecture Notes in Statistics, 195. Springer, Dordrecht, 2009.

GL:2015
J. Gross and U. Ligges.
nortest: Tests for Normality, 2015.
R package version 1.0-4.

G:2011
S. J. Gustafson.
Mathematical Concepts of Quantum Mechanics.
Springer, Heidelberg, 2nd edition, 2011.

Hen:1990
N. Henze.
An approximation to the limit distribution of the Epps–Pulley
  test statistic for normality.
Metrika, 37:7–18, 1990.

HEN:1994b
N. Henze.
The asymptotic behaviour of a variant of multivariate kurtosis.
Communications in Statistics - Theory and Methods,
  23(4):1047–1061, 1994.

HEN:1994a
N. Henze.
On Mardia's kurtosis test for multivariate normality.
Communications in Statistics - Theory and Methods,
  23(4):1031–1045, 1994.

HN:1997
N. Henze.
Extreme smoothing and testing for multivariate normality.
Statistics & Probability Letters, 35(3):203–213, 1997.

H:1997
N. Henze.
Limit laws for multivariate skewness in the sense of Móri,
  Rohatgi and Székely.
Statistics & Probability Letters, 33(3):299–307, 1997.

H:2002
N. Henze.
Invariant tests for multivariate normality: A critical review.
Statistical Papers, 43:467–506, 2002.

HJG:2019
N. Henze and M. D. Jiménez-Gamero.
A class of tests for multinormality with i.i.d. and GARCH data
  based on the empirical moment generating function.
TEST, 28(2):499––521, 2019.

HJM:2019
N. Henze, M. D. Jiménez-Gamero, and S. G. Meintanis.
Characterizations of multinormality and corresponding tests of fit,
  including for GARCH models.
Econometric Theory, 35(3):510––546, 2019.

HM:2019
N. Henze and C. Mayer.
More good news on the hkm test for multivariate reflected symmetry
  about an unknown centre.
Annals of the Institute of Statistical Mathematics, pages
  1–30, 2019.

HNE:2009
N. Henze, Y. Nikitin, and B. Ebner.
Integral distribution-free statistics of lp-type and their asymptotic
  comparison.
Computational Statistics & Data Analysis, 53(9):3426 – 3438,
  2009.

HV:2019
N. Henze and J. Visagie.
Testing for normality in any dimension based on a partial
  differential equation involving the moment generating function.
The Annals of the Institute of Statistical Mathematics,
  https://doi.org/10.1007/s10463-019-00720-8:1––29.

HW:1997
N. Henze and T. Wagner.
A new approach to the BHEP tests for multivariate normality.
Journal of Multivariate Analysis, 61(1):1–23, 1997.

HZ:1990
N. Henze and B. Zirkler.
A class of invariant and consistent tests for multivariate normality.
Communications in Statistics - Theory and Methods,
  19(10):3595–3617, 1990.

KTO:2007
K. Kankainen, S. Taskinen, and H. Oja.
Tests of multinormality based on location vectors and scatter
  matrices.
Statistical Methods & Applications, 16:357–359, 2007.

KP:2018
I. Kim and S. Park.
Likelihood ratio tests for multivariate normality.
Communications in Statistics - Theory and Methods,
  47(8):1923–1934, 2018.

K:2009
É. Krauczi.
A study of the quantile correlation test for normality.
TEST, 18(1):156–165, 2009.

MA:1973
J. Malkovich and A. Afifi.
On tests of multivariate normality.
Journal of the American Statistical Association,
  68(341):176–179, 1973.

MAR:1970
K. V. Mardia.
Measures of multivariate skewness and kurtosis with applications.
Biometrika, 57(3):519–530, 1970.

MAR:1974
K. V. Mardia.
Applications of some measures of multivariate skewness and kurtosis
  in testing normality and robustness studies.
Sankhya Series B, 36(2):115–128, 1974.

MM:2005
C. J. Mecklin and D. J. Mundfrom.
A Monte Carlo comparison of the type I and type II error
  rates of tests of multivariate normality.
Journal of Statistical Computation and Simulation,
  75(2):93–107, 2005.

MRS:1993
T. F. Móri, V. K. Rohatgi, and G. J. Székely.
On multivariate skewness and kurtosis.
Teor. Veroyatn. Primen., 38(3):675–679, 1993.

N:1995
Y. Nikitin.
Asymptotic Efficiency of Nonparametric Tests.
Cambridge University Press, 1995.

P:2005
J. Pudelko.
On a new affine invariant and consistent test for multivariate
  normality.
Probability and Mathematical Statistics, 25(1):43–54, 2005.

R:2019
R Core Team.
R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria, 2019.

RDC:2010
X. Romão, R. Delgado, and A. Costa.
An empirical power comparison of univariate goodness-of-fit tests for
  normality.
Journal of Statistical Computation and Simulation,
  80(5):545–591, 2010.

RKSC:2006
P. Ruckdeschel, M. Kohl, T. Stabla, and F. Camphausen.
S4 classes for distributions.
R News, 6(2):2–6, 2006.

SF:1972
S. S. Shapiro and R. S. Francia.
An approximate analysis of variance test for normality.
Journal of the American Statistical Association,
  67(337):215–216, 1972.

SW:1965
S. S. Shapiro and M. B. Wilk.
An analysis of variance test for normality (complete samples).
Biometrika, 52(3/4):591–611, 1965.

S:2006
B. Sürücü.
Goodness-of-fit tests for multivariate distributions.
Communications in Statistics - Theory and Methods,
  35(7):1319–1331, 2006.

SR:2005
G. J. Székely and M. L. Rizzo.
A new test for multivariate normality.
Journal of Multivariate Analysis, 93:58–80, 2005.

T:2009
C. Tenreiro.
On the choice of the smoothing parameter for the BHEP
  goodness-of-fit test.
Computational Statistics & Data Analysis, 53:1038–1053, 2009.

T:1993
S. Thangavelu.
Lectures on Hermite and Laguerre expansions.
Mathematical Notes ; 42. Princeton Univ. Press, Princeton, NJ, 1993.

T:2014
M. Thulin.
Tests for multivariate normality based on canonical correlations.
Statistical Methods & Applications, 23(2):189–208, 2014.

VPMV:2016
V. Voinov, N. Pya, R. Makarov, and Y. Voinov.
New invariant and consistent chi-squared type goodness-of-fit tests
  for multivariate normality and a related comparative simulation study.
Communications in Statistics - Theory and Methods,
  45(11):3249–3263, 2016.




