# Sub-Dimensional Mardia Measures of Multivariate Skewness and Kurtosis

Joydeep Chowdhury  Subhajit Dutta  Reinaldo B. Arellano-Valle  Marc G. Genton
Statistics Program, King Abdullah University of Science and Technology,
Thuwal, Saudi Arabia  Department of Mathematics and Statistics, Indian
Institute of Technology, Kanpur 208016, India  Department of Statistics,
Pontificia Universida Católica de Chile, Santiago 22, Chile

###### Abstract

Mardia’s measures of multivariate skewness and kurtosis summarize the
respective characteristics of a multivariate distribution with two numbers.
However, these measures do not reflect the sub-dimensional features of the
distribution. Consequently, testing procedures based on these measures may
fail to detect skewness or kurtosis present in a sub-dimension of the
multivariate distribution. We introduce sub-dimensional Mardia measures of
multivariate skewness and kurtosis, and investigate the information they
convey about all sub-dimensional distributions of some symmetric and skewed
families of multivariate distributions. The maxima of the sub-dimensional
Mardia measures of multivariate skewness and kurtosis are considered, as these
reflect the maximum skewness and kurtosis present in the distribution, and
also allow us to identify the sub-dimension bearing the highest skewness and
kurtosis. Asymptotic distributions of the vectors of sub-dimensional Mardia
measures of multivariate skewness and kurtosis are derived, based on which
testing procedures for the presence of skewness and of deviation from Gaussian
kurtosis are developed. The performances of these tests are compared with some
existing tests in the literature on simulated and real datasets.

###### keywords:

asymptotic distribution, measures of multivariate skewness and kurtosis,
multivariate normality test , skew-normal distribution , skew-$t$ distribution
, symmetric distribution

††journal: review

corollarytheorem corollary lemmatheorem lemma propositiontheorem proposition

##  1 Introduction

Consider a $p$-variate random vector $\mathbf{X}=(X_{1},\ldots,X_{p})^{\top}$
from a multivariate distribution with mean vector
$\boldsymbol{\mu}\in\mathbb{R}^{p}$ and $p\times p$ covariance matrix
$\boldsymbol{\Sigma}$. Mardia [1] defined the measures of multivariate
skewness and kurtosis:

$\displaystyle\beta_{1,p}=\text{E}\left[\\{(\mathbf{X}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{Y}-\boldsymbol{\mu})\\}^{3}\right]\in\mathbb{R}_{+},$
(1.1)
$\displaystyle\beta_{2,p}=\text{E}\left[\\{(\mathbf{X}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{X}-\boldsymbol{\mu})\\}^{2}\right]\in\mathbb{R}_{+},$
(1.2)

respectively, where $\mathbf{X}$ and $\mathbf{Y}$ are independent and
identically distributed. These measures are invariant under affine
transformations and each provide a single number to summarize the skewness and
the kurtosis of that $p$-dimensional distribution. Their sample counterparts
have known asymptotic distribution and can be used to test for normality [1,
2]. For the multivariate normal distribution, it is well known that
$\beta_{1,p}=0$ and $\beta_{2,p}=p(p+2)$.

One drawback of Mardia’s measures of multivariate skewness and kurtosis is
that they summarize the information about skewness and kurtosis too much. Sub-
dimensional distributions may exhibit evidence of skewness or kurtosis, which
may not be reflected in the overall Mardia measures of multivariate skewness
or kurtosis. For example, consider the celebrated Fisher’s iris data for the
species ‘iris setosa’ [3]. In Table 1, the p-values of Mardia’s test of
skewness are presented for four sub-dimensions along with the whole dataset,
that is, the sub-dimension $(1,2,3,4)$. Here, the variables 1, 2, 3 and 4
correspond to sepal length, sepal width, petal length and petal width,
respectively.

Table 1: p-values of Mardia skewness test in some sub-dimensions of Fisher’s iris setosa data. Sub-dimensions | (4) | (1, 4) | (2, 4) | (3, 4) | (1, 2, 3, 4)
---|---|---|---|---|---
p-values | 0.001 | 0.012 | 0.019 | 0.018 | 0.236

One can see that in the complete dataset, there is no significant evidence of
skewness, while evidence of skewness in the reported sub-dimensions of
dimension one (consisting of the fourth variable, petal width) and of
dimension two is quite strong at the 5% level as reflected by the
corresponding p-values. Therefore, this motivates the investigation of the
Mardia measures of multivariate skewness and kurtosis on sub-dimensional
marginals.

Let $\mathbf{X}_{qi}$ denote a subvector of dimension $q$ for $1\leq q\leq p$
from the random vector $\mathbf{X}$, and let $\boldsymbol{\mu}_{qi}$ and
$\boldsymbol{\Sigma}_{qi}$ be the corresponding entries of $\boldsymbol{\mu}$
and $\boldsymbol{\Sigma}$, for $i=1,\ldots,Q_{q}$ with $Q_{q}={p\choose q}$.
We define the following sub-dimensional Mardia measures of multivariate
skewness and kurtosis:

$\displaystyle\beta_{1,q,i}=\text{E}\left[\\{(\mathbf{X}_{qi}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{Y}_{qi}-\boldsymbol{\mu}_{qi})\\}^{3}\right]\in\mathbb{R}_{+},$
(1.3)
$\displaystyle\beta_{2,q,i}=\text{E}\left[\\{(\mathbf{X}_{qi}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{X}_{qi}-\boldsymbol{\mu}_{qi})\\}^{2}\right]\in\mathbb{R}_{+},$
(1.4)

for $i=1,\ldots,Q_{q}$, where $\mathbf{X}_{qi}$ and $\mathbf{Y}_{qi}$ are
independent and identically distributed. When $q=p$, we get back
$\beta_{1,p,1}\equiv\beta_{1,p}$ and $\beta_{2,p,1}\equiv\beta_{2,p}$. For
$q=1,\ldots,p$, we collect these measures in the following vectors:

$\displaystyle\mathbf{M}_{1,q}=\left(\beta_{1,q,1},\ldots,\beta_{1,q,Q_{q}}\right)^{\top}\in\mathbb{R}_{+}^{Q_{q}},$
(1.5)
$\displaystyle\mathbf{M}_{2,q}=\left(\beta_{2,q,1},\ldots,\beta_{2,q,Q_{q}}\right)^{\top}\in\mathbb{R}^{Q_{q}}_{+},$
(1.6)

each of dimension $Q_{q}={p\choose q}$. We call (1.5) and (1.6) the $q$-th
vectors of sub-dimensional Mardia measures of multivariate skewness and
kurtosis.

For the multivariate normal distribution, ${\cal
N}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, we have:

$\displaystyle\mathbf{M}_{1,q}^{\cal
N}=\mathbf{0}_{Q_{q}}\quad\text{and}\quad\mathbf{M}_{2,q}^{\cal
N}=q(q+2)\mathbf{1}_{Q_{q}},$ (1.7)

for all $1\leq q\leq p$, where $\mathbf{0}_{Q_{q}}$ and $\mathbf{1}_{Q_{q}}$
are $Q_{q}$-dimensional vectors of zeros and of ones, respectively.

We further define:

$\displaystyle\mathbf{M}_{1}^{*}=\left(\mathbf{M}_{1,1}^{\top},\mathbf{M}_{1,2}^{\top},\ldots,\mathbf{M}_{1,p}^{\top}\right)^{\top}=\left(\beta_{1,1,1},\ldots,\beta_{1,p,1}\right)^{\top}\in\mathbb{R}_{+}^{2^{p}-1},$
(1.8)
$\displaystyle\mathbf{M}_{2}^{*}=\left(\mathbf{M}_{2,1}^{\top},\mathbf{M}_{2,2}^{\top},\ldots,\mathbf{M}_{2,p}^{\top}\right)^{\top}=\left(\beta_{2,1,1},\ldots,\beta_{2,p,1}\right)^{\top}\in\mathbb{R}^{2^{p}-1}_{+}.$
(1.9)

Here, $\mathbf{M}_{1}^{*}$ and $\mathbf{M}_{2}^{*}$ collect all the sub-
dimensional Mardia measures of multivariate skewness and kurtosis.

Figure 1: (A) Contour plot of ${\cal
SN}_{2}(\mathbf{0},\boldsymbol{\Omega},\boldsymbol{\alpha})$, which has
marginal skewness $\beta_{1,1,1}=\beta_{1,1,2}=0.130$. (B) Contour plot of
${\cal SN}_{2}(\mathbf{0},\mathbf{I}_{2},\boldsymbol{\alpha}^{*})$, which has
marginal skewness $\beta_{1,1,1}=0.889$ and $\beta_{1,1,2}=0$. Both skew-
normal distributions have the same Mardia measure of multivariate skewness
$\beta_{1,2}=0.889$.

Due to affine invariance, Mardia’s measures of multivariate skewness and
kurtosis cannot detect any difference between distributions which are affine
transformations of one another, but whose marginal skewness and/or kurtosis
values are very different. In Figure 1, we depict such an example for the
skew-normal distribution [4] with $p=2$. Let
$\boldsymbol{\Omega}=(\omega_{ij})$, where
$\omega_{ij}=0.5+0.5\mathbb{I}(i=j)$ and $\boldsymbol{\alpha}=(5,5)^{\top}$.
Consider the skew-normal distribution ${\cal
SN}_{2}(\mathbf{0},\boldsymbol{\Omega},\boldsymbol{\alpha})$ and its
‘canonical form’ ${\cal
SN}_{2}(\mathbf{0},\mathbf{I}_{2},\boldsymbol{\alpha}^{*})$; see [5,
Proposition 4] and [6, subsection 5.1.8]. The contour plots of the densities
of ${\cal SN}_{2}(\mathbf{0},\boldsymbol{\Omega},\boldsymbol{\alpha})$ and
${\cal SN}_{2}(\mathbf{0},\mathbf{I}_{2},\boldsymbol{\alpha}^{*})$ are
presented in Figure 1. From Propositions 3 and 4 in [5], it follows that the
skew-normal distributions ${\cal
SN}_{2}(\mathbf{0},\boldsymbol{\Omega},\boldsymbol{\alpha})$ and ${\cal
SN}_{2}(\mathbf{0},\mathbf{I}_{2},\boldsymbol{\alpha}^{*})$ have the same
values of Mardia’s measure of multivariate skewness, $\beta_{1,2}=0.889$.
However, the distribution ${\cal
SN}_{2}(\mathbf{0},\mathbf{I}_{2},\boldsymbol{\alpha}^{*})$ has all its
skewness in its first component $X_{1}$ with the marginal distribution of
$X_{2}$ being symmetric, while both the components of the distribution ${\cal
SN}_{2}(\mathbf{0},\boldsymbol{\Omega},\boldsymbol{\alpha})$ have the same
marginal skewness. This further supports the study of Mardia’s measures of
multivariate skewness and kurtosis on sub-dimensional marginals.

In this paper, we start by studying the forms of $\mathbf{M}_{1,q}$ and
$\mathbf{M}_{2,q}$ in (1.5) and (1.6) for some parametric classes of non-
Gaussian distributions, including, for example, the multivariate Student’s $t$
distribution and the multivariate skew-normal and skew-$t$ distributions. Then
we propose tests of normality for each dimension $q$ based on the maximum
entry of $\mathbf{M}_{1,q}$ and of $\mathbf{M}_{2,q}$, as well as global tests
of multivariate normality based on the maximum entry of $\mathbf{M}_{1}^{*}$
and of $\mathbf{M}_{2}^{*}$ in (1.8) and (1.9). One important advantage of our
approach is the ability of the test, when it rejects multivariate normality,
to identify the dimension $q$ and the associated sub-dimensions for which the
rejection occurs.

This paper is organized as follows. The vectors of sub-dimensional Mardia
measures of multivariate skewness and kurtosis in the case of some symmetric
distributions are investigated in section 2, whereas in the case of some
skewed distributions are considered in section 3. The invariance of these sub-
dimensional measures under location-scale transformations is studied in
section 4. The new hypothesis tests are introduced in section 5 and their
asymptotic distributions are established in section 6. The results of a Monte
Carlo simulation study of the empirical sizes and powers of the new tests, as
well as of the sub-dimensional detection, are reported in section 7. Sub-
dimensional data analyses of Fisher’s iris data and of wind speed data near a
wind farm in Saudi Arabia are presented in section 8. The paper ends with a
discussion in section 9.

##  2 Sub-Dimensional Mardia Measures for Some Symmetric Distributions

First, let $\mathbf{X}\buildrel
d\over{=}\mathbf{S}\circ\mathbf{T}\in\mathbb{R}^{p}$ be a standardized
symmetric random vector, in which $\circ$ represent the componentwise
(Hadamard) product,
$\mathbf{S}=\mbox{sign}(\mathbf{X})=(\mbox{sign}(X_{1}),\ldots,\mbox{sign}(X_{p}))^{\top}\sim{\cal
U}\\{-1,+1\\}^{p}$ discrete uniform, and $\mathbf{S}$ is independent of
$\mathbf{T}=|\mathbf{X}|=(|X_{1}|,\ldots,|X_{p}|)^{\top}$. Let
$\mathbf{X}_{\pi}=\mathbf{S}_{\pi}\circ\mathbf{X}=\mathbf{S}_{\pi}\circ\mathbf{S}\circ\mathbf{T}$,
where $\mathbf{S}_{\pi}$ is formed by a permutation
$\pi=(\pi_{1},\ldots,\pi_{p})$ of the components of $\mathbf{S}$. Note that
$\mathbf{s}\circ\mathbf{t}=\mathbf{D}(\mathbf{s})\mathbf{t}=\mathbf{D}(\mathbf{t})\mathbf{s}$,
where $\mathbf{D}(\mathbf{a})=\mbox{diag}(a_{1},\ldots,a_{p})$. Suppose that
$\mathbf{X}$ has finite fourth moment, with E$(\mathbf{X})={\bf 0}_{p}$ and
Var$(\mathbf{X})=\mathbf{I}_{p}$. This implies
$\text{E}(\|\mathbf{X}\|^{2})=p$, where $\|\cdot\|$ denotes the Euclidean
norm. Also, let $\mathbf{X}_{\pi_{*}}=\mathbf{S}_{\pi_{*}}\circ\mathbf{X}$,
where $\pi_{*}$ represents a permutation where no component remains in its
original position. For example, for $p=3$,
$\pi\in\\{(1,2,3),(2,1,3),(2,3,1),(3,2,1),(3,1,2),(1,3,2)\\}$ and
$\pi_{*}\in\\{(2,3,1),(3,1,2)\\}$. Since $\mathbf{S}_{\pi_{*}}\circ\mathbf{S}$
has mean vector ${\bf 0}_{p}$ and covariance matrix $\mathbf{I}_{p}$, we then
have
$\text{E}(\mathbf{X}_{\pi_{*}})=\text{E}(\mathbf{S}_{\pi_{*}}\circ\mathbf{S})\circ\text{E}(|\mathbf{X}|)={\bf
0}_{p}$ and Var$(\mathbf{X}_{\pi_{*}})=\mathbf{I}_{p}$. Also, obviously
$\|\mathbf{X}_{\pi_{*}}\|\buildrel d\over{=}\|\mathbf{X}_{\pi}\|\buildrel
d\over{=}\|\mathbf{X}\|$, hence
$\text{E}(\|\mathbf{X}_{\pi_{*}}\|^{k})=\text{E}(\|\mathbf{X}_{\pi}\|^{k})=\text{E}(\|\mathbf{X}\|^{k})$,
$k\geq 1$. So, for the random vector $\mathbf{X}_{\pi_{*}}$, we have
$\beta_{1,p}=0$ and $\beta_{2,p}=\text{E}(\|\mathbf{X}\|^{4})$. In particular:

  1. 1.

If $\mathbf{X}\sim{\cal N}_{p}({\bf 0}_{p},\mathbf{I}_{p})$, with probability
density function (pdf) given by
$f(\mathbf{x})=(2\pi)^{-p/2}\exp(-\|\mathbf{x}\|^{2})$,
$\mathbf{x}\in\mathbb{R}^{p}$, then $\|\mathbf{X}_{\pi_{*}}\|^{2}\buildrel
d\over{=}\|\mathbf{X}\|^{2}\sim\chi_{p}^{2}$ and $\beta_{2,p}=p(p+2)$.
Therefore, for any $q$-dimensional subvector $\mathbf{X}_{q}$ of $\mathbf{X}$,
we have $\mathbf{X}_{q}\sim{\cal N}_{q}({\bf 0}_{q},\mathbf{I}_{q})$ and so
$\beta_{1,q}=\beta_{1,p}=0$ and $\beta_{2,q}=q(q+2)$. Therefore, the $q$-th
vectors of sub-dimensional Mardia measures of multivariate skewness and
kurtosis are as in (1.7).

  2. 2.

If $\mathbf{X}$ is spherically distributed with pdf
$f(\mathbf{x})=h(\|\mathbf{x}\|^{2})$, $\mathbf{x}\in\mathbb{R}^{p}$, for some
density generator function $h$, i.e., $h(u)>0$ for $u>0$ and
$\int_{0}^{\infty}u^{p/2-1}h(u)\mbox{d}u=\Gamma(p/2)/\pi^{p/2}$, and
Var$(\mathbf{X})=\sigma^{2}\mathbf{I}_{p}$, then $\mathbf{X}\buildrel
d\over{=}R\mathbf{U}^{(p)}$, where $R\buildrel d\over{=}\|\mathbf{X}\|$, with
$\text{E}(R^{2})=p\sigma^{2}$, and $R$ is independent of
$\mathbf{U}^{(p)}\buildrel d\over{=}\mathbf{X}/\|\mathbf{X}\|$, the uniform
vector on the $p$-dimensional unit sphere. Therefore, $\beta_{1,p}=0$ by
symmetry and

$\beta_{2,p}=\frac{p^{2}\text{E}(R^{4})}{\\{\text{E}(R^{2})\\}^{2}}=p(p+2)(\kappa+1),$

where $\kappa=\gamma_{2}=\\{\text{E}(X_{1}^{4})-3\\}/3$ is the excess of
kurtosis in the corresponding spherical univariate distribution, which can be
computed from the relation:

$\kappa+1=\frac{\beta_{2,p}}{p(p+2)}=\frac{p}{p+2}\frac{\text{E}(R^{4})}{\\{\text{E}(R^{2})\\}^{2}},$

with the assumption E$(R^{4})<\infty$. Moreover, any $q$-dimensional subvector
$\mathbf{X}_{q}$ of $\mathbf{X}$, is also spherically distributed with
stochastic representation $\mathbf{X}_{q}=R_{q}\mathbf{U}^{(q)}$, where
$R_{q}=\sqrt{B_{q}}R$ with Beta $B_{q}\sim{\cal B}(q/2,(p-q)/2)$ and $R$,
$B_{q}$ and $\mathbf{U}^{(q)}$ are mutually independent. Using the fact that
$\text{E}(B_{q}^{s})=\\{\Gamma(p/2)\Gamma(q/2+s)\\}/\\{\Gamma(q/2)\Gamma(p/2+s)\\}$,
we find again that $\beta_{1,q}=0$ and

$\beta_{2,q}=\frac{q^{2}\text{E}(B_{q}^{2})\text{E}(R^{4})}{\\{\text{E}(B_{q})\text{E}(R^{2})\\}^{2}}=\frac{q(q+2)}{p(p+2)}\beta_{2,p}=q(q+2)(\kappa+1).$

Hence, for spherical distributions:

$\mathbf{M}_{1,q}^{\cal
SPH}=\mathbf{0}_{Q_{q}}\quad\mbox{and}\quad\mathbf{M}_{2,q}^{\cal
SPH}=q(q+2)(\kappa+1)\mathbf{1}_{Q_{q}},$ (2.1)

for all $1\leq q\leq p.$ For example:
a) If $\mathbf{X}\sim{\cal N}_{p}({\bf 0}_{p},\mathbf{I}_{p})$, then
$R^{2}\sim\chi_{p}^{2}$, with $\text{E}(R^{2})=p$, $\text{E}(R^{4})=p(p+2)$
and so

$\kappa=\frac{p}{p+2}\frac{\text{E}(R^{4})}{\\{\text{E}(R^{2})\\}^{2}}-1=0.$

b) If $\mathbf{X}\sim t_{p}({\bf 0}_{p},\mathbf{I}_{p},\nu)$, where
$t_{p}(\boldsymbol{\xi},\boldsymbol{\Omega},\nu)$ denotes the multivariate
Student’s $t$ distribution with location vector $\boldsymbol{\xi}$, dispersion
matrix $\boldsymbol{\Omega}$, $\nu$ degrees of freedom with $\nu>4$ and pdf
$c_{p}(\nu)|\boldsymbol{\Omega}|^{-p/2}\\{1+(\mathbf{x}-\boldsymbol{\xi})^{\top}\boldsymbol{\Omega}^{-1}(\mathbf{x}-\boldsymbol{\xi})/\nu\\}^{-(\nu+p)/2}$,
$\mathbf{x}\in\mathbb{R}^{p}$, where
$c_{p}(\nu)=\Gamma\\{(\nu+p)/2\\}/\\{\Gamma(\nu/2)(\nu\pi)^{p/2}\\}$, then
$R^{2}\sim pF_{p,\nu}$ with

$\text{E}(R^{2k})=p^{k}\left(\frac{\nu}{p}\right)^{k}\frac{\Gamma(p/2+k)\Gamma(\nu/2-k)}{\Gamma(p/2)\Gamma(\nu/2)},~{}~{}\nu\geq
2k,$

thus $\kappa=2/(\nu-4)$.
c) If $\mathbf{X}\sim{\cal EP}_{p}({\bf 0}_{p},\mathbf{I}_{p},\nu)$, where
${\cal EP}_{p}(\boldsymbol{\xi},\boldsymbol{\Omega},\nu)$ denotes the
multivariate exponential power distribution, with location vector
$\boldsymbol{\xi}$, dispersion matrix $\boldsymbol{\Omega}$ and kurtosis
parameter $\nu>0$, and pdf
$c_{p}(\nu)|\boldsymbol{\Omega}|^{-p/2}\exp[-\\{(\mathbf{x}-\boldsymbol{\xi})^{\top}\boldsymbol{\Omega}^{-1}(\mathbf{x}-\boldsymbol{\xi})\\}^{\nu}/2\\}]$,
$\mathbf{x}\in\mathbb{R}^{p}$, where
$c_{p}(\nu)=\\{p\Gamma(p/2)\\}/\\{\Gamma(p/2\nu+1)2^{p/2\nu+1}\pi^{p/2}\\}$,
then $R^{2}\buildrel d\over{=}V^{1/\nu}$ with Gamma $V\sim{\cal
G}(p/2\nu,1/2)$. Thus, we find that

$\text{E}(R^{2k})=\frac{2^{k/\nu}\Gamma\\{(p+2k)/(2\nu)\\}}{\Gamma\\{p/(2\nu)\\}},~{}~{}k\geq
1,$

so

$\kappa=\frac{p}{p+2}\frac{\Gamma\\{(p+4)/(2\nu)\\}\Gamma\\{p/(2\nu)\\}}{[\Gamma\\{(p+2)/(2\nu)\\}]^{2}}-1.$

In this case, note that $\kappa=\kappa^{(p)}$ depends on the dimension $p$.
This is due to the fact that each marginal distribution of a $p$-variate
exponential power distribution depends on $p$. In particular, the univariate
marginal distributions are not equivalent to the univariate one obtained by
putting $p=1$. Another relevant characteristic of this distribution is that it
allows lighter $(\nu>1)$ and heavier $(0<\nu\leq 1)$ tails than the normal
distribution.

  3. 3.

The class of random vectors $\mathbf{X}_{\pi-\pi_{*}}$ is also very
interesting because E$(\mathbf{X}_{\pi-\pi_{*}})\neq{\bf 0}_{p}$, but some of
its marginal distributions have zero mean vector. For instance, for $p=3$, a
random vector in this class is given by
$(S_{11}|X_{1}|,S_{23}|X_{2}|,S_{23}|X_{3}|)^{\top},$ where
$S_{ij}=S_{i}S_{j}$, in which the first component has mean
$\text{E}(|X_{1}|)>0$, while the remaining components have zero mean.

##  3 Sub-Dimensional Mardia Measures for Some Skewed Distributions

Consider a generalized skew-normal [7] random vector $\mathbf{X}\sim{\cal
GSN}_{p}({\bf 0}_{p},\boldsymbol{\Omega},\boldsymbol{\lambda})$ with
probability density function
$f_{\mathbf{X}}(\mathbf{x})=2\phi_{p}(\mathbf{x};\boldsymbol{\Omega})G(\boldsymbol{\lambda}^{\top}\mathbf{x})$,
$\mathbf{x}\in\mathbb{R}^{p}$ and some absolutely continuous symmetric
cumulative distribution function $G$ defined on $\mathbb{R}$. Since
$\beta_{1,p}$ and $\beta_{2,p}$ are invariant with respect to nonsingular
linear transformations, they can be computed by using the canonical
representations $\mathbf{Z}=\boldsymbol{\Gamma}\mathbf{X}$ and
$\mathbf{Z}^{\prime}=\boldsymbol{\Gamma}\mathbf{Y}$ of $\mathbf{X}$ and
$\mathbf{Y}$, respectively, in which $\boldsymbol{\Gamma}$ is a $p\times p$
orthonormal matrix such that
$\boldsymbol{\Gamma}\boldsymbol{\Omega}\boldsymbol{\Gamma}^{\top}=\mathbf{I}_{p}$
and $\boldsymbol{\Gamma}\boldsymbol{\lambda}=\lambda_{*}{\bf e}_{1:p}$, where
$\lambda_{*}=\sqrt{\boldsymbol{\lambda}^{\top}\boldsymbol{\Omega}\boldsymbol{\lambda}}$
and ${\bf e}_{1}$ is the first $p$-dimensional unit vector. Then
$Z_{1}\sim{\cal GSN}_{1}(0,1,\lambda_{*})$, which has mean $\mu_{*}$ and
variance $1-\mu_{*}^{2}$, and $Z_{1}$ is independent of
$\mathbf{Z}_{2}=(Z_{2},\ldots,Z_{p})^{\top}\sim{\cal N}_{p-1}({\bf
0}_{p-1},\mathbf{I}_{p-1})$, where
$\mu_{*}=2\lambda_{*}\text{E}\\{G^{\prime}(\lambda_{*}Y)\\}=2\frac{d\mbox{\scriptsize
E}\\{G(\lambda_{*}Y)\\}}{d\lambda_{*}},$ with $Y\sim{\cal HN}_{1}(0,1)$ half-
normal. The same holds for $Z_{1}^{\prime}$ and
$\mathbf{Z}_{2}^{\prime}=(Z_{2}^{\prime},\ldots,Z_{p}^{\prime})^{\top}$. Note
that $\boldsymbol{\mu}_{\mathbf{Z}}=\mu_{*}{\bf e}_{1}$ and
$\boldsymbol{\Sigma}_{\mathbf{Z}}=\mathbf{I}_{p}-\mu_{*}^{2}{\bf e}_{1}{\bf
e}_{1}^{\top}.$ Since
$(\mathbf{X}-\boldsymbol{\mu}_{\mathbf{X}})^{\top}\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}(\mathbf{X}^{\prime}-\boldsymbol{\mu}_{\mathbf{X}})=(\mathbf{Z}-\boldsymbol{\mu}_{\mathbf{Z}})^{\top}\boldsymbol{\Sigma}_{\mathbf{Z}}^{-1}(\mathbf{Z}^{\prime}-\boldsymbol{\mu}_{\mathbf{Z}})=Z_{01}Z_{01}^{\prime}+\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2}^{\prime},$
where $Z_{01}=(Z_{1}-\mu_{*})/\sqrt{1-\mu_{*}^{2}}$ and
$Z_{01}^{\prime}=(Z_{1}^{\prime}-\mu_{*})/\sqrt{1-\mu_{*}^{2}}$ are
standardized and independent random variables:

$\displaystyle\beta_{1,p}$
$\displaystyle=\text{E}\\{(Z_{01}Z_{01}^{\prime}+\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2}^{\prime})^{3}\\}=\text{E}\\{(Z_{01}Z_{01}^{\prime})^{3}\\}+3\text{E}\\{(Z_{01}Z_{01}^{\prime})^{2}(\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2}^{\prime})\\}+3\text{E}\\{(Z_{01}Z_{01}^{\prime})(\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2}^{\prime})^{2}\\}+\text{E}\\{(\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2}^{\prime})^{3}\\}$
$\displaystyle=\text{E}\\{(Z_{01}Z_{01}^{\prime})^{3}\\}=\gamma_{1,1}^{2},$

where $\gamma_{1,1}\equiv\gamma_{1}=\sqrt{\beta_{1,1}}=\text{E}(Z_{01}^{3})$.
Similarly, we have
$(\mathbf{X}-\boldsymbol{\mu}_{\mathbf{X}})^{\top}\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}(\mathbf{X}-\boldsymbol{\mu}_{\mathbf{X}})=Z_{01}^{2}+\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2}$
and thus:

$\displaystyle\beta_{2,p}$
$\displaystyle=\text{E}\\{(Z_{01}^{2}+\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2})^{2}\\}=\text{E}(Z_{01}^{4})+2\text{E}(Z_{01}^{2})\text{E}\\{(\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2})\\}+\text{E}\\{(\mathbf{Z}_{2}^{\top}\mathbf{Z}_{2})^{2}\\}=\text{E}(Z_{01}^{4})=\gamma_{2,1}+p(p+2),$

with $\gamma_{2,1}\equiv\gamma_{2}=\text{E}(Z_{01}^{4})-3=\beta_{2,1}-3$ being
the excess of kurtosis. In this case, we can also consider the Mardia
multivariate excess of kurtosis index $\gamma_{2,p}=\beta_{2,p}-p(p+2)$ which
equals $\gamma_{2,1}$.

In the skew-normal case with $G(x)=\Phi(x)$, we have
$\gamma_{1}=b(2b^{2}-1)\gamma_{*}^{3/2}$ and
$\gamma_{2}=2b^{2}(2-3b^{2})\gamma_{*}^{2}$, where
$\gamma_{*}=\frac{\lambda_{*}^{2}}{1+(1-b^{2})\lambda_{*}^{2}}$,
$b=\sqrt{2/\pi}$ and
$\lambda_{*}=\sqrt{\boldsymbol{\lambda}^{\top}\boldsymbol{\Omega}\boldsymbol{\lambda}}$
as defined above. However, the indices $\gamma_{1}$ and $\gamma_{2}$
corresponding to the subvectors of $\mathbf{X}$ are different through
subvectors with different dimensions. As was shown in [5], if
$\mathbf{X}\sim{\cal
SN}_{p}(\mathbf{0}_{p},\boldsymbol{\Omega},\boldsymbol{\lambda})$ and
$\mathbf{A}$ is a $q\times p$ fixed matrix, then
$\mathbf{A}\mathbf{X}\sim{\cal
SN}_{q}\left(\mathbf{0}_{q},\boldsymbol{\Omega}_{A},\boldsymbol{\lambda}_{A}\right),$
where $\boldsymbol{\Omega}_{A}=\mathbf{A}\boldsymbol{\Omega}\mathbf{A}^{\top}$
and
$\boldsymbol{\lambda}_{A}=\frac{(\mathbf{A}\boldsymbol{\Omega}\mathbf{A}^{\top})^{-1}\mathbf{A}\boldsymbol{\Omega}\boldsymbol{\lambda}}{\sqrt{1+\boldsymbol{\lambda}^{\top}(\boldsymbol{\Omega}-\boldsymbol{\Omega}\mathbf{A}^{\top}\boldsymbol{\Omega}_{A}^{-1}\mathbf{A}\boldsymbol{\Omega})\boldsymbol{\lambda}}}.$
In particular, if we consider the partition
$\mathbf{X}=(\mathbf{X}_{q}^{\top},\mathbf{X}_{p-q}^{\top})^{\top}$ with the
corresponding partition for the scale matrix
$\boldsymbol{\Omega}=(\boldsymbol{\Omega}_{rs})_{r,s=q,p-q}$ and the skewness
vector
$\boldsymbol{\lambda}=(\boldsymbol{\lambda}_{q}^{\top},\boldsymbol{\lambda}_{p-q}^{\top})^{\top}$,
we then have:

$\mathbf{X}_{q}=\mathbf{A}\mathbf{X}\sim{\cal
SN}_{q}\left(\mathbf{0}_{q},\boldsymbol{\Omega}_{qq},\boldsymbol{\lambda}^{(q)}\right),\quad\boldsymbol{\lambda}^{(q)}=\frac{\boldsymbol{\lambda}_{q}+\boldsymbol{\Omega}_{qq}^{-1}\boldsymbol{\Omega}_{q,p-q}\boldsymbol{\lambda}_{p-q}}{\sqrt{1+\boldsymbol{\lambda}_{p-q}^{\top}\boldsymbol{\Omega}_{p-q,p-q:q}\boldsymbol{\lambda}_{p-q}}},$

where $\mathbf{A}\mathbf{A}^{\top}=\mathbf{I}_{q}$ and
$\boldsymbol{\Omega}_{p-q,p-q:q}=\boldsymbol{\Omega}_{p-q,p-q}-\boldsymbol{\Omega}_{p-q,q}\boldsymbol{\Omega}_{qq}^{-1}\boldsymbol{\Omega}_{q,p-q}$.
The associated canonical transformation to $\mathbf{X}_{q}$ has summary
skewness parameter:

$\lambda_{*}^{(q)}=\sqrt{\boldsymbol{\lambda}^{(q)\top}\boldsymbol{\Omega}_{qq}\boldsymbol{\lambda}^{(q)}}=\sqrt{\frac{(\boldsymbol{\lambda}_{q}+\boldsymbol{\Omega}_{qq}^{-1}\boldsymbol{\Omega}_{q,p-q}\boldsymbol{\lambda}_{p-q})^{\top}\boldsymbol{\Omega}_{qq}(\boldsymbol{\lambda}_{q}+\boldsymbol{\Omega}_{qq}^{-1}\boldsymbol{\Omega}_{q,p-q}\boldsymbol{\lambda}_{p-q})}{1+\boldsymbol{\lambda}_{p-q}^{\top}\boldsymbol{\Omega}_{p-q,p-q:q}\boldsymbol{\lambda}_{p-q}}}.$

Therefore, we have $\beta_{1,q,i}=\gamma_{1}^{2(q)}$ and
$\beta_{2,q,i}=\gamma_{2}^{(q)}+q(q+2)$, where for $k=1,2$,
$\gamma_{k}^{(q)}=\gamma_{k}$ for $q=p$, and for $q<p$ it must be computed as
$\gamma_{k}$ but with $\lambda_{*}$ replaced by $\lambda_{*}^{(q)}$. In
particular, if $\boldsymbol{\Omega}=\mathbf{I}_{p}$, then
$\lambda_{*}^{(q)}=\sqrt{\frac{\boldsymbol{\lambda}_{q}^{\top}\boldsymbol{\lambda}_{q}}{1+\boldsymbol{\lambda}_{p-q}^{\top}\boldsymbol{\lambda}_{p-q}}}$,
with $\lambda_{*}^{(q)}=\lambda_{*}$ if $q=p$, and
$\lambda_{*}^{(q)}=\lambda_{i}/\sqrt{1+\sum_{j\neq i}\lambda_{j}^{2}}$ for the
$i$-th marginal component if $q=1$.

Next, we consider the multivariate skew-$t$ distribution as described in [6,
section 6.2]. Let $\mathbf{X}\sim{\cal
ST}_{p}(\boldsymbol{\xi},\boldsymbol{\Omega},\boldsymbol{\alpha},\nu)$, where
${\cal ST}_{p}$ denotes a $p$-dimensional skew-$t$ distribution with degree of
freedom $\nu$, location vector $\boldsymbol{\xi}$, scale matrix
$\boldsymbol{\Omega}$ and skewness vector $\boldsymbol{\alpha}$. Define:

$\displaystyle\delta_{*}=\left(\frac{\boldsymbol{\alpha}^{\top}\boldsymbol{\Omega}\boldsymbol{\alpha}}{1+\boldsymbol{\alpha}^{\top}\boldsymbol{\Omega}\boldsymbol{\alpha}}\right)^{1/2}\quad\text{and}\quad
b_{\nu}=\frac{\sqrt{\nu}\Gamma\\{(\nu-1)/2\\}}{\sqrt{\pi}\Gamma(\nu/2)}.$

Based on $\delta_{*}$ and $b_{\nu}$, we set $\mu_{*}=b_{\nu}\delta_{*}$ and
$\sigma_{*}^{2}=\\{\nu/(\nu-2)\\}-\mu_{*}^{2}$. Mardia’s measures of
multivariate skewness and kurtosis for $\mathbf{X}$ are [6]:

$\displaystyle\beta_{1,p}=\beta_{1}^{*}+3(p-1)\frac{\mu_{*}^{2}}{(\nu-3)\sigma_{*}^{2}}\;\text{if
}\nu>3,$ (3.1)
$\displaystyle\beta_{2,p}=\beta_{2}^{*}+(p^{2}-1)\frac{\nu-2}{\nu-4}+\frac{2(p-1)}{\sigma_{*}^{2}}\left\\{\frac{\nu}{\nu-4}-\frac{(\nu-1)\mu_{*}^{2}}{\nu-3}\right\\}-p(p+2)\;\text{if
}\nu>4,$ (3.2)

where

$\displaystyle\beta_{1}^{*}=\frac{\mu_{*}^{2}}{\sigma_{*}^{3}}\left\\{\frac{\nu(3-\delta_{*}^{2})}{\nu-3}-\frac{3\nu}{\nu-2}+2\mu_{*}^{2}\right\\}^{2},$
$\displaystyle\beta_{2}^{*}=\frac{1}{\sigma_{*}^{4}}\left\\{\frac{3\nu^{2}}{(\nu-2)(\nu-4)}-\frac{4\mu_{*}^{2}\nu(3-\delta_{*}^{2})}{\nu-3}+\frac{6\mu_{*}^{2}\nu}{\nu-2}-3\mu_{*}^{4}\right\\}.$

Next, denote the subvector corresponding to $\beta_{1,q,i}$ and
$\beta_{2,q,i}$ as $\mathbf{X}_{q}$, and consider the partition
$\mathbf{X}=(\mathbf{X}_{q}^{\top},\mathbf{X}_{p-q}^{\top})^{\top}$ with the
corresponding partitions for the location vector
$\boldsymbol{\xi}=(\boldsymbol{\xi}_{q}^{\top},\boldsymbol{\xi}_{p-q}^{\top})^{\top}$,
scale matrix $\boldsymbol{\Omega}=(\boldsymbol{\Omega}_{rs})_{r,s=q,p-q}$ and
the skewness vector
$\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{q}^{\top},\boldsymbol{\alpha}_{p-q}^{\top})^{\top}$
as in the earlier part of this section. Then:

$\displaystyle\mathbf{X}_{q}\sim{\cal
ST}_{q}\left(\boldsymbol{\xi}_{q},\boldsymbol{\Omega}_{qq},\boldsymbol{\alpha}^{(q)},\nu\right),\quad\boldsymbol{\alpha}^{(q)}=\frac{\boldsymbol{\alpha}_{q}+\boldsymbol{\Omega}_{qq}^{-1}\boldsymbol{\Omega}_{q,p-q}\boldsymbol{\alpha}_{p-q}}{\sqrt{1+\boldsymbol{\alpha}_{p-q}^{\top}\boldsymbol{\Omega}_{p-q,p-q:q}\boldsymbol{\alpha}_{p-q}}},$

where
$\boldsymbol{\Omega}_{p-q,p-q:q}=\boldsymbol{\Omega}_{p-q,p-q}-\boldsymbol{\Omega}_{p-q,q}\boldsymbol{\Omega}_{qq}^{-1}\boldsymbol{\Omega}_{q,p-q}$.
Now, $\beta_{1,q,i}$ and $\beta_{2,q,i}$ are obtained by replacing
$\boldsymbol{\alpha}$ and $\boldsymbol{\Omega}$ with
$\boldsymbol{\alpha}^{(q)}$ and $\boldsymbol{\Omega}_{qq}$, respectively, in
the formulae given in (3.1) and (3.2).

##  4 Invariance Under Location-Scale Transformations

It is well known that Mardia’s measures of multivariate skewness and kurtosis
(1.1) and (1.2) are invariant under affine transformations. We show next that
the sub-dimensional Mardia measures of multivariate skewness and kurtosis
(1.3) and (1.4) are only invariant under location and scale transformation,
unless the multivariate distribution is spherically invariant.

Let ${\rm
MD}^{2}(\mathbf{X})=(\mathbf{X}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{X}-\boldsymbol{\mu})$
be the squared Mahalanobis distance. Consider the affine transformation
$\mathbf{Y}=\mathbf{A}\mathbf{X}+\mathbf{b}$,
$\mathbf{A}\in\mathbb{R}^{p\times p}$, with $|\mathbf{A}|\neq 0$,
$\mathbf{b}\in\mathbb{R}^{p}$. Then, it is immediate that ${\rm
MD}^{2}(\mathbf{Y})={\rm MD}^{2}(\mathbf{A}\mathbf{X}+\mathbf{b})={\rm
MD}^{2}(\mathbf{X})={\rm MD}^{2}(\mathbf{Z})$, where
$\mathbf{Z}=\boldsymbol{\Sigma}^{-1/2}(\mathbf{X}-\boldsymbol{\mu})$.

Next, partition
$\mathbf{X}=(\mathbf{X}_{q}^{\top},\mathbf{X}_{p-q}^{\top})^{\top}$ and
similarly for $\mathbf{Y}$ and $\mathbf{Z}$. Then,
$\mathbf{X}_{q}=\mathbf{B}_{q}\mathbf{X}$, where
$\mathbf{B}_{q}=(\mathbf{I}_{q},\,{\bf O})\in\mathbb{R}^{q\times p}$ with
rank$(\mathbf{B}_{q})=q$. Similarly, $\mathbf{Y}_{q}=\mathbf{B}_{q}\mathbf{Y}$
and $\mathbf{Z}_{q}=\mathbf{B}_{q}\mathbf{Z}$. Therefore:

$\displaystyle{\rm MD}^{2}(\mathbf{X}_{q})$ $\displaystyle=$
$\displaystyle(\mathbf{B}_{q}\mathbf{X}-\mathbf{B}_{q}\boldsymbol{\mu})^{\top}(\mathbf{B}_{q}\boldsymbol{\Sigma}\mathbf{B}_{q}^{\top})^{-1}(\mathbf{B}_{q}\mathbf{X}-\mathbf{B}_{q}\boldsymbol{\mu})=(\mathbf{X}-\boldsymbol{\mu})^{\top}\mathbf{B}_{q}^{\top}(\mathbf{B}_{q}\boldsymbol{\Sigma}\mathbf{B}_{q}^{\top})^{-1}\mathbf{B}_{q}(\mathbf{X}-\boldsymbol{\mu})$
$\displaystyle=$
$\displaystyle\mathbf{Z}^{\top}\boldsymbol{\Sigma}^{1/2}\mathbf{B}_{q}^{\top}(\mathbf{B}_{q}\boldsymbol{\Sigma}\mathbf{B}_{q}^{\top})^{-1}\mathbf{B}_{q}\boldsymbol{\Sigma}^{1/2}\mathbf{Z}=\mathbf{Z}^{\top}\mathbf{P}\mathbf{Z},$

where $\mathbf{P}\geq 0$ is a $p\times p$ orthogonal projection matrix, that
is, symmetric with $\mathbf{P}^{2}=\mathbf{P}$ and rank$(\mathbf{P})=q$.
Similarly, ${\rm
MD}^{2}(\mathbf{Y}_{q})=\mathbf{Z}^{\top}\mathbf{P}_{\mathbf{A}}\mathbf{Z}$
with
$\mathbf{P}_{\mathbf{A}}=\boldsymbol{\Sigma}^{1/2}\mathbf{A}^{\top}\mathbf{B}_{q}^{\top}(\mathbf{B}_{q}\mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^{\top}\mathbf{B}_{q}^{\top})^{-1}\mathbf{B}_{q}\mathbf{A}\boldsymbol{\Sigma}^{1/2}$,
where $\mathbf{P}_{\mathbf{A}}\geq 0$ is a $p\times p$ orthogonal projection
matrix, that is, symmetric with
$\mathbf{P}_{\mathbf{A}}^{2}=\mathbf{P}_{\mathbf{A}}$ and
rank$(\mathbf{P}_{\mathbf{A}})=q$. Therefore, if
$\mathbf{P}_{\mathbf{A}}=\mathbf{P}$ then ${\rm MD}^{2}(\mathbf{X}_{q})={\rm
MD}^{2}(\mathbf{Y}_{q})={\rm
MD}^{2}(\mathbf{Z}_{q})=\mathbf{Z}_{q}^{\top}\mathbf{Z}_{q}$. In particular,
this holds when the matrix $\mathbf{A}$ is diagonal, which proves invariance
under location and scale transformations.

More generally, simple calculations show that
$\mathbf{P}_{\mathbf{A}}=\mathbf{P}$ if and only if:

$\mathbf{A}_{12}\boldsymbol{\Sigma}_{21}\mathbf{A}_{11}^{\top}+\mathbf{A}_{11}\boldsymbol{\Sigma}_{12}\mathbf{A}_{12}^{\top}+\mathbf{A}_{12}\boldsymbol{\Sigma}_{22}\mathbf{A}_{12}^{\top}={\bf
O}.$ (4.1)

In addition to the trivial case $\mathbf{A}_{12}={\bf O}$, Equation (4.1)
holds, for instance, when
$\mathbf{A}_{12}\perp(\boldsymbol{\Sigma}_{21},\boldsymbol{\Sigma}_{22})$,
where $\perp$ means orthogonal.

Moreover, if $\mathbf{Z}$ is spherically distributed, that is,
$\boldsymbol{\Gamma}\mathbf{X}\buildrel d\over{=}\mathbf{X}$ for any
orthogonal matrix $\boldsymbol{\Gamma}$, then the invariance holds for any
affine transformation with $|\mathbf{A}|\neq 0$. Indeed, in this case
$\mathbf{P}_{\mathbf{A}}=\boldsymbol{\Gamma}_{\mathbf{A}}^{\top}\mathbf{D}_{\mathbf{A}}\boldsymbol{\Gamma}_{\mathbf{A}}$
with $\mathbf{D}_{\mathbf{A}}=\mathbf{B}_{q}^{\top}\mathbf{B}_{q}$ for some
orthogonal matrix $\boldsymbol{\Gamma}_{\mathbf{A}}$ and therefore:

$\mathbf{Z}^{\top}\mathbf{P}_{\mathbf{A}}\mathbf{Z}=\mathbf{Z}^{\top}\boldsymbol{\Gamma}^{\top}_{\mathbf{A}}\mathbf{B}_{q}^{\top}\mathbf{B}_{q}\boldsymbol{\Gamma}_{\mathbf{A}}\mathbf{Z}=(\boldsymbol{\Gamma}_{\mathbf{A}}\mathbf{Z})^{\top}\mathbf{B}_{q}^{\top}\mathbf{B}_{q}(\boldsymbol{\Gamma}_{\mathbf{A}}\mathbf{Z})\buildrel
d\over{=}\mathbf{Z}^{\top}\mathbf{B}_{q}^{\top}\mathbf{B}_{q}\mathbf{Z}=\mathbf{Z}_{q}^{\top}\mathbf{Z}_{q},$

and similarly when $\mathbf{A}=\mathbf{I}_{p}$.

In summary, although Mardia’s measures of multivariate skewness and kurtosis
are invariant under affine transformations, this is generally not the case for
the sub-dimensional measures as shown in this section and illustrated in the
simple case of a bivariate skew-normal distribution in Figure 1. Therefore,
the sub-dimensional Mardia measures of multivariate skewness and kurtosis are
informative for testing normality in sub-dimensions as proposed in the next
section.

##  5 Sub-Dimensional Estimation and Testing of Hypotheses

Let $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ be a random sample from a
distribution $F$ in $\mathbb{R}^{p}$. Let $\mathbf{X}_{qij}$ denote a sub-
vector of dimension $q$ obtained from $\mathbf{X}_{j}$, $j=1,\ldots,n$. The
dimension of $\mathbf{X}_{qij}$ corresponds to that of $\mathbf{X}_{qi}$
defined in section 1. More precisely, suppose
$\mathbf{X}_{qi}=\mathbf{P}_{qi}\mathbf{X}$, where $\mathbf{P}$ is a diagonal
matrix whose $i$th diagonal element is either 1 or 0 (depending on whether the
$i$th coordinate of $\mathbf{X}$ is included in $\mathbf{X}_{qi}$ or not).
Then, $\mathbf{X}_{qij}=\mathbf{P}_{qi}\mathbf{X}_{j}$ for all $j$. Let
$\bar{\mathbf{X}}_{qi}$ and $\mathbf{S}_{qi}$ be the sample mean and the
sample covariance matrix of the observations
$\mathbf{X}_{qi1},\ldots,\mathbf{X}_{qin}$, respectively. The sample
estimators of $\beta_{1,q,i}$ and $\beta_{2,q,i}$ defined in (1.3) and (1.4)
are:

$\displaystyle
b_{1,q,i}=\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}\\{(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\mathbf{S}_{qi}^{-1}(\mathbf{X}_{qik}-\bar{\mathbf{X}}_{qi})\\}^{3},$
$\displaystyle
b_{2,q,i}=\frac{1}{n}\sum_{j=1}^{n}\\{(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\mathbf{S}_{qi}^{-1}(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})\\}^{2},$

respectively. Consequently, the estimates of $\mathbf{M}_{1,q}$ and
$\mathbf{M}_{2,q}$ defined in (1.5) and (1.6) are
$\mathbf{m}_{1,q}=(b_{1,q,1},\ldots,b_{1,q,Q_{q}})^{\top}$ and
$\mathbf{m}_{2,q}=(b_{2,q,1},\ldots,b_{2,q,Q_{q}})^{\top}$, respectively.
Similarly, the estimates of $\mathbf{M}_{1}^{*}$ and $\mathbf{M}_{2}^{*}$
defined in (1.8) and (1.9) are
$\mathbf{m}_{1}^{*}=(b_{1,1,1},\ldots,b_{1,p,1})^{\top}$ and
$\mathbf{m}_{2}^{*}=(b_{2,1,1},\ldots,b_{2,p,1})^{\top}$, respectively.

The quantities $\mathbf{m}_{1,q}$ and $\mathbf{m}_{2,q}$ provide information
about the skewness and kurtosis present in all the $q$-dimensional sub-vectors
constructed from the sample. On the other hand, $\mathbf{m}_{1}^{*}$ and
$\mathbf{m}_{2}^{*}$ reflect the skewness and kurtosis present in all possible
sub-dimensions of the sample. Based on these quantities, tests of skewness and
kurtosis can be constructed.

###  5.1 Testing Skewness

From (2.4) in [1], it follows that $\beta_{1,q,i}\geq 0$ for all $q,i$.
However, for symmetric distributions, $\beta_{1,q,i}=0$ for all $q,i$ as
discussed in section 2,. So, if any of the sub-dimensions bears skewness, we
shall have $\beta_{1,q,i}>0$ for the $(q,i)$-pair corresponding to that sub-
dimension, which implies that
$\max\\{\beta_{1,q,i}\;|\;q=1,\ldots,p;i=1,\ldots,Q_{q}\\}>0$. So, a
hypothesis, which tests for skewness in all sub-dimensions of the
distribution, can be formulated as follows:

$\displaystyle\text{H}_{0}^{(s)}:\max_{q,i}\beta_{1,q,i}=0\;\;\text{and}\;\;\text{H}_{\text{A}}^{(s)}:\max_{q,i}\beta_{1,q,i}>0.$
(5.1)

The above hypothesis is equivalent to:

$\displaystyle\tilde{\text{H}}_{0}^{(s)}:\beta_{1,q,i}=0\text{ for all
}q,i\;\;\text{and}\;\;\tilde{\text{H}}_{\text{A}}^{(s)}:\beta_{1,q,i}>0\text{
for some }q,i.$

The usual Mardia’s skewness test [1] only tests for:

$\displaystyle\bar{\text{H}}_{0}^{(s)}:\beta_{1,p}=0\;\;\text{and}\;\;\bar{\text{H}}_{\text{A}}^{(s)}:\beta_{1,p}>0,$

and thus it may be less efficient in providing information about skewness
supported on a smaller sub-dimension. In this aspect, testing for (5.1) can be
expected to be more efficient than the usual Mardia’s skewness test.

The null hypothesis $\text{H}_{0}^{(s)}$ in (5.1) should be rejected when the
maximum of $b_{1,q,i}$ is large. However, directly comparing the quantities
$b_{1,q,i}$ is not proper, because they have different means and standard
deviations under the null hypothesis. When the underlying distribution $F$ is
Gaussian, it follows from Equation (2.26) in [1] that the asymptotic
expectation and asymptotic standard deviation of $nb_{1,q,i}$ is $q(q+1)(q+2)$
and $\sqrt{12q(q+1)(q+2)}$, respectively. So, while comparing the quantities
$b_{1,q,i}$, it is appropriate to center and scale them first. Because our aim
is to detect non-Gaussian features in the sample, we center and scale
$b_{1,q,i}$ using its asymptotic expectation and standard deviation under
Gaussianity, and consider:

$\displaystyle\tilde{b}_{1,q,i}=\frac{nb_{1,q,i}-q(q+1)(q+2)}{\sqrt{12q(q+1)(q+2)}},\quad
q=1,\ldots,p;\;i=1,\ldots,Q_{q}.$

We reject $\text{H}_{0}^{(s)}$ in (5.1) when $\max_{q,i}\tilde{b}_{1,q,i}$ is
large, and we shall denote this test as the MaxS test. Let
$\tilde{\mathbf{m}}_{1,q}$ and $\tilde{\mathbf{m}}_{1}^{*}$ be the centered
and scaled analogues of $\mathbf{m}_{1,q}$ and $\mathbf{m}_{1}^{*}$, formed by
replacing $b_{1,q,i}$ by $\tilde{b}_{1,q,i}$. Note that
$\max_{q,i}\tilde{b}_{1,q,i}$ is the maximum of the coordinates of
$\tilde{\mathbf{m}}_{1}^{*}$. To find the p-value of the MaxS test or to
construct the cutoff for a given level, we need the asymptotic distribution of
$\max_{q,i}\tilde{b}_{1,q,i}$, which is derived in section 6.

Sometimes, it may be the case that there is information about possible
presence of skewness in the sub-vectors of a fixed dimension, say, $q_{0}$,
where $q_{0}<p$, but it is unknown exactly which subvector has a skewed
distribution. Then, it would be judicious to construct the test based on
$\max_{i}\tilde{b}_{1,q_{0},i}$ only. In such a situation, the hypothesis is:

$\displaystyle\text{H}_{0}^{q_{0},(s)}:\max_{i}\beta_{1,q_{0},i}=0\;\;\text{and}\;\;\text{H}_{\text{A}}^{q_{0},(s)}:\max_{i}\beta_{1,q_{0},i}>0.$
(5.2)

Here, the null hypothesis $\text{H}_{0}^{q_{0},(s)}$ in (5.2) is rejected for
large values of $\max_{i}\tilde{b}_{1,q_{0},i}$. We denote this test as the
$\text{MaxS}_{q_{0}}$ test. To find the p-value (or the cutoff) for rejection
at a given level, we derive the asymptotic distribution of
$\max_{i}\tilde{b}_{1,q_{0},i}$ in section 6.

###  5.2 Testing Kurtosis

Similar to tests of skewness, tests of kurtosis can be also constructed based
on the $b_{2,q,i}$ quantities. For a $p$-dimensional Gaussian distribution,
the value of Mardia’s measure of kurtosis $\beta_{2,p}$ is $p(p+2)$; see [1].
The quantity $\beta_{2,p}$ measures how heavy-tailed the distribution is. For
distributions with tails heavier than the Gaussian distribution, e.g., the
multivariate Student’s $t$ distributions, we have $\beta_{2,p}>p(p+2)$.
Similarly, for distributions with lighter tails than the Gaussian
distribution, we have $\beta_{2,p}<p(p+2)$. However, Mardia’s test of kurtosis
[1] again checks the overall kurtosis of all the dimensions, and would not be
efficient in checking whether the distribution in a particular sub-dimension
deviates from Gaussianity in terms of kurtosis.

Suppose one wants to test whether the kurtosis of the distribution of any sub-
dimension deviates from the Gaussian distribution. If this is the case, then
we have $\left|\beta_{2,q,i}-q(q+2)\right|>0$, where the pair $(q,i)$
corresponds to that sub-dimension. For this scenario, the suitable hypotheses
are:

$\displaystyle\text{H}_{0}^{(k)}:\max_{q,i}\left|\beta_{2,q,i}-q(q+2)\right|=0\;\;\text{and}\;\;\tilde{\text{H}}_{\text{A}}^{(k)}:\max_{q,i}\left|\beta_{2,q,i}-q(q+2)\right|>0.$
(5.3)

Like in the case of $b_{1,q,i}$, $b_{2,q,i}$ has different expectation and
standard deviation for different $q$. When the underlying distribution $F$ is
Gaussian, from the derivations in subsection 3.2 in [1], the asymptotic
expectation and the asymptotic standard deviation of $b_{2,q,i}$ are $q(q+2)$
and $\sqrt{\\{8q(q+2)\\}/n}$. So, we center and scale $b_{2,q,i}$, and
consider:

$\displaystyle\tilde{b}_{2,q,i}=\frac{b_{2,q,i}-q(q+2)}{\sqrt{\\{8q(q+2)\\}/n}},\quad
q=1,\ldots,p;\;i=1,\ldots,Q_{q}.$

The null hypothesis $\text{H}_{0}^{(k)}$ in (5.3) is rejected if
$\max_{q,i}|\tilde{b}_{2,q,i}|$ is large. We denote this test as the MaxK
test. To find the p-value (or the cutoff) of the test at a given level, we use
the asymptotic distribution of $\max_{q,i}|\tilde{b}_{2,q,i}|$ derived in
section 6.

Now, suppose one knows that there is possible deviation from the Gaussian
distribution in terms of kurtosis in some sub-dimension with dimension
$q_{0}<p$, but the exact sub-dimension is unknown. In such a case, the
appropriate hypotheses would be:

$\displaystyle\text{H}_{0}^{q_{0},(k)}:\max_{i}\left|\beta_{2,q_{0},i}-q_{0}(q_{0}+2)\right|=0\;\;\text{and}\;\;\text{H}_{\text{A}}^{q_{0},(k)}:\max_{i}\left|\beta_{2,q_{0},i}-q_{0}(q_{0}+2)\right|>0.$
(5.4)

Here, the null hypothesis $\text{H}_{0}^{q_{0},(k)}$ in (5.4) is rejected if
$\max_{i}|\tilde{b}_{2,q_{0},i}|$ is large. This test is denoted as the
$\text{MaxK}_{q_{0}}$ test. We employ the asymptotic distribution of
$\max_{i}|\tilde{b}_{2,q_{0},i}|$ derived in section 6 to find the p-value (or
the cutoff) of the test at a given level.

##  6 Asymptotic Distributions and Implementation of Tests

In this section, the asymptotic distributions of the quantities
$\max_{q,i}\tilde{b}_{1,q,i}$, $\max_{i}\tilde{b}_{1,q_{0},i}$,
$\max_{q,i}|\tilde{b}_{2,q,i}|$ and $\max_{i}|\tilde{b}_{2,q_{0},i}|$
introduced in section 5 are derived. Based on the respective asymptotic
distributions, the implementations of the MaxS, $\text{MaxS}_{q_{0}}$, MaxK
and $\text{MaxK}_{q_{0}}$ tests corresponding to (5.1), (5.2), (5.3) and
(5.4), respectively, are described. To derive the aforementioned asymptotic
distributions, we shall use certain linearizations corresponding to the terms
$b_{1,q,i}$ and $b_{2,q,i}$.

###  6.1 Skewness

Given any $q,i$, define:

$\displaystyle h_{qi}(\mathbf{x},\mathbf{y})$
$\displaystyle=\left\\{(\mathbf{x}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{y}-\boldsymbol{\mu}_{qi})\right\\}^{3}-3(\mathbf{x}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{x}-\boldsymbol{\mu}_{qi})(\mathbf{x}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{y}-\boldsymbol{\mu}_{qi})$
$\displaystyle\quad-3(\mathbf{y}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{y}-\boldsymbol{\mu}_{qi})(\mathbf{x}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{y}-\boldsymbol{\mu}_{qi})+3(q+2)(\mathbf{x}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{y}-\boldsymbol{\mu}_{qi}).$
(6.1)

Define the integral operator $\mathbf{h}_{qi}$ by:

$\displaystyle\\{\mathbf{h}_{qi}(g)\\}(\mathbf{x})=\text{E}\left\\{h_{qi}(\mathbf{x},\mathbf{X}_{qi})g(\mathbf{X}_{qi})\right\\}.$
(6.2)

If $\text{E}(\|\mathbf{X}\|^{3})<\infty$, then the integral operator
$\mathbf{h}_{qi}$ is well defined on $L_{2}[\mathbf{X}_{qi}]$, the space of
measurable functions $g$ which are square integrable with respect to the
distribution of $\mathbf{X}_{qi}$. Under appropriate assumptions, it can be
established that $\mathbf{h}_{qi}$ has only finitely many nonzero eigenvalues,
where the number $K(q)$ of nonzero eigenvalues depends on $q$ (see [8]). Let
$\lambda_{qik}$ be the eigenvalues and $f_{qik}(\cdot)$ be the corresponding
eigenfunctions of $\mathbf{h}_{qi}$ for $k=1,2,\ldots,K(q)$. Then, we have:

$\displaystyle
h_{qi}(\mathbf{x},\mathbf{y})=\sum_{k=1}^{K(q)}\lambda_{qik}f_{qik}(\mathbf{x})f_{qik}(\mathbf{y})$
(6.3)

for all $\mathbf{x},\mathbf{y}$. The function $h_{qi}(\mathbf{x},\mathbf{y})$
is closely related to the quantity $b_{1,q,i}$. The eigenvalues and
eigenfunctions of $h_{qi}(\mathbf{x},\mathbf{y})$ are used to establish a
relationship between $b_{1,q,i}$ and the average of independent random vectors
in the following theorem. This linearization will be used to derive the
asymptotic distributions of the quantities $\max_{q,i}\tilde{b}_{1,q,i}$ and
$\max_{i}\tilde{b}_{1,q_{0},i}$. Define:

$\displaystyle\mathbf{m}_{4,qi}$
$\displaystyle=\text{E}\left[\left\\{(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})\right\\}^{2}\right],$
$\displaystyle\mathbf{m}_{6,qi}$
$\displaystyle=\text{E}\left[\left\\{(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})\right\\}^{3}\right].$

We have the following theorem on $b_{1,q,i}$.

######  Theorem 6.1.

Suppose that $\text{E}(\|\mathbf{X}\|^{6})<\infty$ and the distribution of
$\mathbf{X}$ is elliptical. Let
$\sigma_{qi}^{2}=\text{E}\\{(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})^{2}\\}$
if $q=1$. For $j=1,\ldots,n$, define $\mathbf{u}_{qij}$ as:

$\displaystyle\mathbf{u}_{qij}=\begin{cases}&\sigma_{qi}^{-3}(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})\\{(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})^{2}-3\sigma_{qi}^{2}\\}\text{
if }q=1,\\\
&\left(\sqrt{\lambda_{qi1}}f_{qi1}(\mathbf{X}_{qij}),\ldots,\sqrt{\lambda_{qiK(q)}}f_{qiK(q)}(\mathbf{X}_{qij})\right)^{\top}\text{
otherwise}.\end{cases}$

Let $\bar{\mathbf{u}}_{qi}=n^{-1}\sum_{j=1}^{n}\mathbf{u}_{qij}$. Then:

$\displaystyle nb_{1,q,i}$
$\displaystyle=\left\|\sqrt{n}\bar{\mathbf{u}}_{qi}\right\|_{2}^{2}+o_{P}(1)$

as $n\to\infty$. Further, for $q>1$, $K(q)=q+q(q-1)(q+4)/6$:

$\displaystyle\lambda_{qik}=\begin{cases}&(3/q)\left\\{\mathbf{m}_{6,qi}/(q+2)-2\mathbf{m}_{4,qi}+(q+2)q\right\\}\text{
for }k=1,\ldots,q,\\\ &6\mathbf{m}_{6,qi}/\left\\{q(q+2)(q+4)\right\\}\text{
for }k=(q+1),\ldots,K(q),\end{cases}$

and $\text{E}\\{f_{qik}(\mathbf{X}_{qij})\\}=0$ for all $q,i,j$ and $k$.

###### Proof.

When $q=1$, from the arguments in the proof of Theorem 1 in [9], it follows
that, as $n\to\infty$:

$\displaystyle nb_{1,qi}$
$\displaystyle=\left[n^{-1/2}\sum_{j=1}^{n}\sigma_{qi}^{-3}(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})\\{(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})^{2}-3\sigma_{qi}^{2}\\}\right]^{2}+o_{P}(1)=\left(\sqrt{n}\bar{\mathbf{u}}_{qi}\right)^{2}+o_{P}(1).$

For $q>1$, let
$\mathbf{Y}_{qij}=\boldsymbol{\Sigma}_{qi}^{-1/2}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)$
for all $i,j$. Then $\mathbf{Y}_{qij}$ has an elliptical distribution, which
is identical over $j$ for fixed $i$, with
$\text{E}(\mathbf{Y}_{qij})=\mathbf{0}$ and
$\text{E}(\mathbf{Y}_{qij}\mathbf{Y}_{qij}^{\top})=\mathbf{I}_{q}$. Also,
$\text{E}(\|\mathbf{Y}_{qij}\|^{6})<\infty$ for all $i,j$, and hence
$\mathbf{Y}_{qij}$s satisfy the conditions of Lemma 2.1 in [8]. From an
application of this lemma:

$\displaystyle
nb_{1,qi}=n^{-1}\sum_{j=1}^{n}\sum_{k=1}^{n}h_{qi}(\mathbf{X}_{qij},\mathbf{X}_{qik})+o_{P}(1)$
(6.4)

as $n\to\infty$. Next, from the arguments in the proof of Theorem 2.2 in [8],
it follows that the integral operator $\mathbf{h}_{qi}$ defined in Equation
(6.2) has only two non-zero distinct eigenvalues, which are:

$\displaystyle\gamma_{qi1}=(3/q)\left\\{\mathbf{m}_{6,qi}/(q+2)-2\mathbf{m}_{4,qi}+(q+2)q\right\\},$
$\displaystyle\gamma_{qi2}=6\mathbf{m}_{6,qi}/\left\\{q(q+2)(q+4)\right\\},$

with associated multiplicities $\nu_{qi1}=q$ and $\nu_{qi2}=q(q-1)(q+4)/6$,
respectively. So, we can take $\lambda_{qik}=\gamma_{qi1}$ for
$k=1,\ldots,\nu_{qi1}$ and $\lambda_{qik}=\gamma_{qi2}$ for
$k=(\nu_{qi1}+1),\ldots,(\nu_{qi1}+\nu_{qi2})$. Consequently,
$K(q)=q+q(q-1)(q+4)/6$. From (6.3) and (6.4), we get:

$\displaystyle nb_{1,qi}$
$\displaystyle=n^{-1}\sum_{j=1}^{n}\sum_{l=1}^{n}\sum_{k=1}^{K(q)}\lambda_{qik}f_{qik}(\mathbf{X}_{qij})f_{qik}(\mathbf{X}_{qil})+o_{P}(1)=\sum_{k=1}^{K(q)}\left\\{n^{-1/2}\sum_{j=1}^{n}\sqrt{\lambda_{qik}}f_{qik}(\mathbf{X}_{qij})\right\\}^{2}+o_{P}(1)$
$\displaystyle=\left\|\sqrt{n}\bar{\mathbf{u}}_{qi}\right\|_{2}^{2}+o_{P}(1)$

as $n\to\infty$. Finally, from the arguments in the proof of Theorem 2.1 in
[10], it follows that $\text{E}\\{f_{qik}(\mathbf{X}_{qij})\\}=0$ for all
$q,i,j$ and $k$. ∎

Define
$\mathbf{U}_{j}=\left(\mathbf{u}_{1,1,j}^{\top},\ldots,\mathbf{u}_{q,i,j}^{\top},\ldots,\mathbf{u}_{p,1,j}^{\top}\right)^{\top}$
and
$\mathbf{U}_{qj}=\left(\mathbf{u}_{q,1,j}^{\top},\ldots,\mathbf{u}_{q,Q_{q},j}^{\top}\right)^{\top}$.
Let $G(\cdot)$ and $G_{q}(\cdot)$ be such that:

$\displaystyle
G(\mathbf{U}_{j})=\max\left\\{\frac{\left\|\mathbf{u}_{1,1,j}\right\|_{2}^{2}-6}{\sqrt{72}},\ldots,\frac{\left\|\mathbf{u}_{q,i,j}\right\|_{2}^{2}-q(q+1)(q+2)}{\sqrt{12q(q+1)(q+2)}},\ldots,\frac{\left\|\mathbf{u}_{p,1,j}\right\|_{2}^{2}-p(p+1)(p+2)}{\sqrt{12p(p+1)(p+2)}}\right\\},$

and

$\displaystyle
G_{q}(\mathbf{U}_{qj})=\max\left\\{\frac{\left\|\mathbf{u}_{q,1,j}\right\|_{2}^{2}-q(q+1)(q+2)}{\sqrt{12q(q+1)(q+2)}},\ldots,\frac{\left\|\mathbf{u}_{q,1,j}\right\|_{2}^{2}-q(q+1)(q+2)}{\sqrt{12q(q+1)(q+2)}}\right\\}.$

Clearly, $G(\cdot)$ and $G_{q}(\cdot)$ are continuous functions. From these
observations, we derive the asymptotic null distributions of the test
statistics $\max_{q,i}\tilde{b}_{1,q,i}$ for the MaxS test and
$\max_{i}\tilde{b}_{1,q_{0},i}$ for the $\text{MaxS}_{q_{0}}$ test in the next
theorem.

######  Theorem 6.2.

Let $\boldsymbol{\Omega}$ and $\boldsymbol{\Omega}_{q}$ be the dispersion
matrices of $\mathbf{U}_{1}$ and $\mathbf{U}_{q1}$, respectively. Let
$\mathbf{W}$ and $\mathbf{W}_{q_{0}}$ be zero-mean Gaussian random vectors
with dispersion matrices $\boldsymbol{\Omega}$ and
$\boldsymbol{\Omega}_{q_{0}}$, respectively. Assume
$\text{E}(\|\mathbf{X}\|^{6})<\infty$ and the distribution of $\mathbf{X}$ is
elliptical. Then, as $n\to\infty$:
$\max_{q,i}n\tilde{b}_{1,q,i}\stackrel{{\scriptstyle
d}}{{\rightarrow}}G(\mathbf{W})$ and
$\max_{i}n\tilde{b}_{1,q_{0},i}\stackrel{{\scriptstyle
d}}{{\rightarrow}}G_{q_{0}}(\mathbf{W}_{q_{0}})$.

###### Proof.

Since we have $\text{E}\\{f_{qik}(\mathbf{X}_{qij})\\}=0$ for all $q,i,j$ and
$k$ from Theorem 6.1, it follows that $\text{E}(\mathbf{U}_{qj})=\bf 0$ and
$\text{E}(\mathbf{U}_{j})=\bf 0$ for all $j$. Further, the distributions of
$\mathbf{U}_{qj}$ are independent and identical for all $j$, and the same is
true for the distributions of $\mathbf{U}_{j}$. Define
$\bar{\mathbf{U}}=n^{-1}\sum_{j=1}^{n}\mathbf{U}_{j}$ and
$\bar{\mathbf{U}}_{q}=n^{-1}\sum_{j=1}^{n}\mathbf{U}_{qj}$. It follows from
the multivariate central limit theorem that
$\sqrt{n}\bar{\mathbf{U}}\stackrel{{\scriptstyle d}}{{\rightarrow}}\mathbf{W}$
and $\sqrt{n}\bar{\mathbf{U}}_{q_{0}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathbf{W}_{q_{0}}$ as $n\to\infty$. Now, note that
$\max_{q,i}n\tilde{b}_{1,q,i}=G(\sqrt{n}\bar{\mathbf{U}})$ and
$\max_{i}n\tilde{b}_{1,q_{0},i}=\sqrt{n}\bar{\mathbf{U}}_{q_{0}}$. Since
$G(\cdot)$ and $G_{q}(\cdot)$ are continuous, the proof of the theorem is
completed from an application of the continuous mapping theorem. ∎

We implement the tests of skewness under Gaussianity of the null hypotheses.
To derive the p-values of the tests of the hypotheses described in (5.1) and
(5.2), we need to first estimate $\boldsymbol{\Omega}$ and
$\boldsymbol{\Omega}_{q_{0}}$, the dispersion matrices of $\mathbf{U}_{1}$ and
$\mathbf{U}_{q_{0}}$, respectively. However, for that, we need to estimate the
random vectors $\mathbf{U}_{j}$ and $\mathbf{U}_{q_{0}j}$, which involve
unknown population quantities $\boldsymbol{\Sigma}_{qi}$,
$\boldsymbol{\mu}_{qi}$, $\lambda_{qik}$ and $f_{qik}(\mathbf{X}_{qij})$. Here
$\boldsymbol{\Sigma}_{qi}$ is estimated by the sample covariance matrix
$\widehat{\boldsymbol{\Sigma}}_{qi}=(n-1)^{-1}\sum_{j=1}^{n}(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}$
and $\boldsymbol{\mu}_{qi}$ is estimated by the corresponding sample mean
$\widehat{\boldsymbol{\mu}}_{qi}=n^{-1}\sum_{j=1}^{n}\mathbf{X}_{qij}$. Since
we are testing for Gaussianity, $\lambda_{qik}$ is derived considering the
underlying null distribution to be Gaussian, and it can be verified that
$\lambda_{qik}=6$ for all $q,i,k$ in this case. It can also be verified that
under Gaussianity, when $q=1$, $\text{Var}(\mathbf{u}_{qij})=6$ for all $i,j$.
So, while estimating $\mathbf{u}_{qij}$ for $q=1$, first we compute
$\tilde{\mathbf{u}}_{qij}=\widehat{\sigma}_{qi}^{-3}(\mathbf{X}_{qij}-\widehat{\boldsymbol{\mu}}_{qi})\\{(\mathbf{X}_{qij}-\widehat{\boldsymbol{\mu}}_{qi})^{2}-3\widehat{\sigma}_{qi}^{2}\\}$,
where $\widehat{\sigma}_{qi}^{2}$ is the appropriate component of
$\widehat{\boldsymbol{\Sigma}}_{qi}$. Then $\mathbf{u}_{qij}$ is estimated by
$\widehat{\mathbf{u}}_{qij}=\left\\{\sqrt{6}/\sqrt{\text{Var}_{n}(\tilde{\mathbf{u}}_{qij})}\right\\}\tilde{\mathbf{u}}_{qij}$,
where $\text{Var}_{n}(\tilde{\mathbf{u}}_{qij})$ denotes the sample variance
of $\tilde{\mathbf{u}}_{qij}$. In this way, it is ensured that
$\text{Var}_{n}(\tilde{\mathbf{u}}_{qij})=6$. Next, recall that the function
$f_{qik}(\cdot)$ is the eigenfunction of the integral operator
$\mathbf{h}_{qi}$ corresponding to the eigenvalue $\lambda_{qik}$. From the
arguments in Theorem 2.2 in [8], it follows that $f_{qik}(\cdot)$ are
spherical harmonic functions (see [11]). Explicit expressions of these
spherical harmonic functions can be derived (see [12], [13]). However, their
numerical approximation may be unstable due to the involvement of
hypergeometric functions (see [12, p. 1554]). For this reason and ease of
computation, the random variables $f_{qik}(\mathbf{X}_{qij})$ are estimated
through the eigendecomposition of the $n\times n$ matrix
$\widehat{H}_{qi}=(\widehat{h}_{qi}(\mathbf{X}_{qij},\mathbf{X}_{qil}))$,
where:

$\displaystyle\widehat{h}_{qi}(\mathbf{X}_{qij},\mathbf{X}_{qil})$
$\displaystyle=\left\\{(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qil}-\bar{\mathbf{X}}_{qi})\right\\}^{3}+3(q+2)\left\\{(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qil}-\bar{\mathbf{X}}_{qi})\right\\}$
$\displaystyle\quad-3(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qil}-\bar{\mathbf{X}}_{qi})$
$\displaystyle\quad-3(\mathbf{X}_{qil}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qil}-\bar{\mathbf{X}}_{qi})(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qil}-\bar{\mathbf{X}}_{qi}).$

We compute the $K(q)$ eigenvectors of the matrix $\widehat{H}_{qi}$
corresponding to the $K(q)$ eigenvalues with the largest magnitudes, and
arrange them by descending order of the magnitudes of their corresponding
eigenvalues. The eigenvalues are then discarded, and we use $\lambda_{qik}=6$
assuming Gaussianity. Each computed eigenvector is multiplied by $\sqrt{n}$ to
maintain its correspondence with the integral operator $\mathbf{h}_{qi}$. Each
such eigenvector, denoted as $\widehat{f}_{qik}$, is the estimated value of
the vector
$(f_{qik}(\mathbf{X}_{qi1}),\ldots,f_{qik}(\mathbf{X}_{qin}))^{\top}$. From
these estimated values, we estimate the random vectors $\mathbf{U}_{j}$ and
$\mathbf{U}_{q_{0}j}$, and the corresponding estimates are denoted as
$\widehat{\mathbf{U}}_{j}$ and $\widehat{\mathbf{U}}_{q_{0}j}$, respectively.
The dispersion matrices $\boldsymbol{\Omega}$ and
$\boldsymbol{\Omega}_{q_{0}}$ are estimated by:

$\displaystyle\widehat{\boldsymbol{\Omega}}$
$\displaystyle=(n-1)^{-1}\sum_{j=1}^{n}\left(\widehat{\mathbf{U}}_{j}-n^{-1}\sum_{l=1}^{n}\widehat{\mathbf{U}}_{l}\right)\left(\widehat{\mathbf{U}}_{j}-n^{-1}\sum_{l=1}^{n}\widehat{\mathbf{U}}_{l}\right)^{\top},$
$\displaystyle\widehat{\boldsymbol{\Omega}}_{q_{0}}$
$\displaystyle=(n-1)^{-1}\sum_{j=1}^{n}\left(\widehat{\mathbf{U}}_{q_{0}j}-n^{-1}\sum_{l=1}^{n}\widehat{\mathbf{U}}_{q_{0}l}\right)\left(\widehat{\mathbf{U}}_{q_{0}j}-n^{-1}\sum_{l=1}^{n}\widehat{\mathbf{U}}_{q_{0}l}\right)^{\top}.$

Now, to compute the p-values of the MaxS test for (5.1), we generate 1000
independent zero-mean Gaussian random vectors
$\tilde{W}_{1},\ldots,\tilde{W}_{1000}$ with the dispersion matrix
$\widehat{\boldsymbol{\Omega}}$. The proportion of the values
$G(\tilde{W}_{1}),\ldots,G(\tilde{W}_{1000})$ larger than
$\max_{q,i}n\tilde{b}_{1,q,i}$ is taken as the p-value of the null hypothesis
in (5.1). The p-value of the test $\text{MaxS}_{q_{0}}$ in (5.2) is derived
similarly.

###  6.2 Kurtosis

Next, we derive the asymptotic null distributions of
$\max_{q,i}|\tilde{b}_{2,q,i}|$ and of $\max_{i}|\tilde{b}_{2,q_{0},i}|$. Here
also, we first derive a linearization of $b_{2,q,i}$.

######  Theorem 6.3.

Let $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ be independent and identically
distributed random vectors in $\mathbb{R}^{p}$ with
$\text{E}(\|\mathbf{X}\|^{8})<\infty$. Let
$\mathbf{A}_{qi}=\boldsymbol{\Sigma}_{qi}^{-1}\text{E}\\{(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi})^{\top}\\}$.
Define:

$\displaystyle\mathbf{Z}_{qij}=\begin{bmatrix}\left\\{\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)\right\\}^{2}-\text{E}\left[\left\\{\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)\right\\}^{2}\right]\\\
\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)-\text{E}\left\\{\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)\right\\}\\\
\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\end{bmatrix},$

and

$\displaystyle\mathbf{a}_{qi}=\left(1,-2,-4\text{E}\left[\left\\{\left(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi}\right)^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi}\right)\right\\}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi}\right)\right]\right)^{\top}.$

Then,
$n^{1/2}(b_{2,q,i}-\beta_{2,q,i})=n^{-1/2}\sum_{j=1}^{n}\mathbf{a}_{qi}^{\top}\mathbf{Z}_{qij}+o_{P}(1)$
as $n\to\infty$.

###### Proof.

Recall that $\beta_{2,q,i}$ and $b_{2,q,i}$ are invariant under location
transformations. So, without loss of generality, we can assume
$\boldsymbol{\mu}=\text{E}(\mathbf{X})=\mathbf{0}$, which means
$\boldsymbol{\mu}_{qi}=\mathbf{0}$ and
$\text{Var}(\mathbf{X})=\text{E}(\mathbf{X}\mathbf{X}^{\top})=\boldsymbol{\Sigma}$.

The arguments are similar to those in the proof of Theorem 2.1 in [14]. In
this proof, $\mathbf{Y}_{n}=O_{P}(b_{n})$ means that the sequence
$\mathbf{Y}_{n}/b_{n}$ is bounded in Euclidean/matrix norm, while
$\mathbf{Y}_{n}=o_{P}(b_{n})$ means that
$b_{n}^{-1}\mathbf{Y}_{n}\stackrel{{\scriptstyle P}}{{\to}}0$. Let:

$\displaystyle\mathbf{B}_{qin}=n^{1/2}\left\\{n^{-1}\sum_{j=1}^{n}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)^{\top}-\boldsymbol{\Sigma}_{qi}\right\\}.$

From the weak law of large numbers, we get $\mathbf{B}_{qi}=O_{P}(1)$ for all
$q,i$. Also:

$\displaystyle
n^{1/2}\left(\mathbf{S}_{qi}-\boldsymbol{\Sigma}_{qi}\right)=\mathbf{B}_{qi}-n^{1/2}\left(\bar{\mathbf{X}}_{qi}-\boldsymbol{\mu}_{qi}\right)\left(\bar{\mathbf{X}}_{qi}-\boldsymbol{\mu}_{qi}\right)^{\top}.$
(6.5)

From the multivariate central limit theorem, we have that
$n^{1/2}\left(\bar{\mathbf{X}}_{qi}-\boldsymbol{\mu}_{qi}\right)\left(\bar{\mathbf{X}}_{qi}-\boldsymbol{\mu}_{qi}\right)^{\top}=O_{P}(n^{-1/2})$,
which, applied to (6.5), yields:

$\displaystyle\mathbf{S}_{qi}=\boldsymbol{\Sigma}_{qi}+n^{-1/2}\mathbf{B}_{qi}+O_{P}(n^{-1}),\text{or,
}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{S}_{qi}=\mathbf{I}_{q}+n^{-1/2}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qi}+O_{P}(n^{-1}).$

Note that
$n^{1/2}(\mathbf{S}_{qi}-\boldsymbol{\Sigma}_{qi})=\mathbf{B}_{qin}-n^{1/2}\bar{\mathbf{X}}_{qi}\bar{\mathbf{X}}_{qi}^{\top},$
where:

$\mathbf{B}_{qin}=n^{1/2}\left(\frac{1}{n}\sum_{j=1}^{n}{\mathbf{X}}_{qij}{\mathbf{X}}_{qij}^{\top}-\boldsymbol{\Sigma}_{qi}\right)=O_{P}(1)$

and $n^{1/2}\bar{\mathbf{X}}_{qi}\bar{\mathbf{X}}_{qi}^{T}=O_{P}(n^{-1/2})$ by
the multivariate central limit theorem. It follows that:

$\mathbf{S}_{qi}=\boldsymbol{\Sigma}_{qi}+n^{-1/2}\mathbf{B}_{qin}+O_{P}(n^{-1}),$

and thus:

$\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{S}_{qi}=\mathbf{I}_{qi}+n^{-1/2}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}+O_{P}(n^{-1}).$

This means $\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{S}_{qi}$ is invertible for
all sufficiently large $n$ with probability approaching 1, and we have:

$(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{S}_{qi})^{-1}=\mathbf{I}_{qi}-n^{-1/2}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}+O_{P}(n^{-1}),$

which implies that:

$\mathbf{S}_{qi}^{-1}=\boldsymbol{\Sigma}_{qi}^{-1}-n^{-1/2}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\boldsymbol{\Sigma}_{qi}^{-1}+O_{P}(n^{-1}).$
(6.6)

Now:

$\displaystyle\\{(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\mathbf{S}_{qi}^{-1}(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})\\}^{2}$
$\displaystyle=(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})^{2}+4(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})^{2}+(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\bar{\mathbf{X}}_{qi})^{2}$
$\displaystyle\quad-4(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij}\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})+2(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij}\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\bar{\mathbf{X}}_{qi})$
$\displaystyle\quad-4(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij}\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\bar{\mathbf{X}}_{qi}).$
(6.7)

Using this and the expression for the inverse in (6.6) above, we get for the
first term in (6.7) to be:

$\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij}=\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}-n^{-1/2}\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}+O_{P}(n^{-1}).$

Squaring both sides, we get:

$\displaystyle(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})^{2}$
$\displaystyle=(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})^{2}-2n^{-1/2}\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}+O_{P}(n^{-1}).$

Using the fact that the trace of a matrix is invariant under cyclic
permutations, we have:

$\displaystyle\text{tr}(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})=\text{tr}\\{\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}\mathbf{X}_{qij}^{\top})^{2}\\}.$

Since
$n^{-1}\sum_{j=1}^{n}(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}\mathbf{X}_{qij}^{\top})^{2}=\mathbf{A}_{qi}+o_{P}(1)$
with
$\mathbf{A}_{qi}=\text{E}\\{(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qi}\mathbf{X}_{qi}^{\top})^{2}\\}$,
we now get:

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})^{2}=\frac{1}{n}\sum_{j=1}^{n}(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})^{2}-2n^{-1/2}\text{tr}(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\mathbf{A}_{qi})+o_{P}(n^{-1}).$

For the second term in (6.7), again using equation (6.6) and the fact that
$n^{-1}\sum_{j=1}^{n}\mathbf{X}_{qij}\mathbf{X}_{qij}^{\top}=O_{P}(1)$ (using
the weak law of large numbers) we obtain:

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})^{2}=O_{P}(n^{-1/2}).$

Furthermore, for the third term in (6.7):
$(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\bar{\mathbf{X}}_{qi})^{2}=O_{P}(n^{-2})$.
Since,
$n^{-1}\sum_{j=1}^{n}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})=\mathbf{A}_{1,q,i}+o_{P}(1)$
with
$\mathbf{A}_{1,q,i}=\text{E}\\{\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qi}(\mathbf{X}_{qi}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qi})\\}$,
it is easy to see that for the fourth term:

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})=\mathbf{A}_{1,q,i}\bar{\mathbf{X}}_{qi}+O_{P}(n^{-1}).$

Finally, for the last two terms:

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}(\mathbf{X}_{qij}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\bar{\mathbf{X}}_{qi})=O_{P}(n^{-1}),\text{and}\,\frac{1}{n}\sum_{j=1}^{n}(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\mathbf{X}_{qij})(\bar{\mathbf{X}}_{qi}^{\top}\mathbf{S}_{qi}^{-1}\bar{\mathbf{X}}_{qi})=O_{P}(n^{-1}).$

Summarizing, we obtain:

$\displaystyle
b_{2,qi}=\frac{1}{n}\sum_{j=1}^{n}(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})^{2}-2n^{-1/2}\text{tr}(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\mathbf{A}_{qi})-4\mathbf{A}_{1,q,i}\bar{\mathbf{X}}_{qi}+O_{P}(n^{-1}).$

Observing that:

$\displaystyle\text{tr}(\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{B}_{qin}\mathbf{A}_{qi})$
$\displaystyle=\text{tr}(\mathbf{B}_{qin}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1})=n^{1/2}\left\\{\frac{1}{n}\sum_{j=1}^{n}\mathbf{X}_{qij}^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}-\text{E}(\mathbf{X}_{qi}^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qi})\right\\}$

we obtain:

$\displaystyle
n^{1/2}(b_{2,qi}-\beta_{2,qi})=n^{-1/2}\sum_{j=1}^{n}\mathbf{a}_{qi}^{\top}\mathbf{Z}_{qij}+o_{P}(1),$

where:

$\displaystyle\mathbf{Z}_{qij}=\begin{bmatrix}\left(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}\right)^{2}-\text{E}\\{(\mathbf{X}_{qij}^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})^{2}\\}\\\
\mathbf{X}_{qij}^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij}-\text{E}(\mathbf{X}_{qij}^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\mathbf{X}_{qij})\\\
\mathbf{X}_{qij}\end{bmatrix}$

is a $(2+q)$-dimensional vector and
$\mathbf{a}_{qi}=(1,-2,-4\mathbf{A}_{1,q,i})^{\top}$. ∎

Tests of kurtosis are usually conducted for testing normality (see [14]). When
the underlying distribution is Gaussian, we can simplify the quantities
$\mathbf{A}_{qi}$ and $\mathbf{a}_{qi}$ described in Theorem 6.3, and
consequently, the linearization becomes simpler.

######  Corollary .

Let $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ be independent and identically
distributed Gaussian random vectors in $\mathbb{R}^{p}$. Define
$\mathbf{Y}_{qij}=\left\\{(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})\right\\}^{2}-2(q+2)(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})^{\top}\boldsymbol{\Sigma}_{qi}^{-1}(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}).$
Then,
$n^{1/2}(b_{2,q,i}-\beta_{2,q,i})=n^{-1/2}\sum_{j=1}^{n}\\{\mathbf{Y}_{qij}-\text{E}(\mathbf{Y}_{qij})\\}+o_{P}(1)$
as $n\to\infty$.

###### Proof.

Define
$\tilde{\mathbf{X}}_{qij}=\boldsymbol{\Sigma}_{qi}^{-1/2}(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi})$.
Under the assumption of the corollary, $\tilde{\mathbf{X}}_{qij}$s are
independent zero-mean Gaussian random vectors with the identity matrix as
their dispersion matrix. It follows that:

$\displaystyle\text{E}\left[\left\\{\left(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi}\right)^{\top}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi}\right)\right\\}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qi1}-\boldsymbol{\mu}_{qi}\right)\right]$
$\displaystyle=\boldsymbol{\Sigma}_{qi}^{-1/2}\text{E}\left\\{\left(\tilde{\mathbf{X}}_{qi1}^{\top}\tilde{\mathbf{X}}_{qi1}\right)\tilde{\mathbf{X}}_{qi1}\right\\}$
$\displaystyle=\boldsymbol{\Sigma}_{qi}^{-1/2}\times\bf 0=\bf 0.$ (6.8)

Next, we have:

$\displaystyle\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)^{\top}\mathbf{A}_{qi}\boldsymbol{\Sigma}_{qi}^{-1}\left(\mathbf{X}_{qij}-\boldsymbol{\mu}_{qi}\right)=\tilde{\mathbf{X}}_{qij}^{\top}\text{E}\left[\tilde{\mathbf{X}}_{qij}\left(\tilde{\mathbf{X}}_{qij}^{\top}\tilde{\mathbf{X}}_{qij}\right)\tilde{\mathbf{X}}_{qij}^{\top}\right]\tilde{\mathbf{X}}_{qij}=(q+2)\tilde{\mathbf{X}}_{qij}^{\top}\tilde{\mathbf{X}}_{qij}.$
(6.9)

The proof follows from (6.8), (6.9) and the linearization in Theorem 6.3. ∎

The asymptotic null distributions of $\max_{q,i}|\tilde{b}_{2,q,i}|$ and
$\max_{i}|\tilde{b}_{2,q_{0},i}|$ are derived from the linearization in
subsection 6.2.

######  Theorem 6.4.

Let $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ be independent and identically
distributed Gaussian random vectors in $\mathbb{R}^{p}$. Define:

$\displaystyle\tilde{\mathbf{Y}}_{j}=\left(\frac{\mathbf{Y}_{1,1,j}}{\sqrt{24}},\ldots,\frac{\mathbf{Y}_{q,i,j}}{\sqrt{8q(q+2)}},\ldots,\frac{\mathbf{Y}_{p,1,j}}{\sqrt{8p(p+2)}}\right)^{\top},$
$\displaystyle\tilde{\mathbf{Y}}_{q,j}=\left(\frac{\mathbf{Y}_{q,1,j}}{\sqrt{8q(q+2)}},\ldots,\frac{\mathbf{Y}_{q,Q_{q},j}}{\sqrt{8q(q+2)}}\right)^{\top}.$

Let $\boldsymbol{\Gamma}$ and $\boldsymbol{\Gamma}_{q}$ be the dispersion
matrices of $\tilde{\mathbf{Y}}_{1}$ and $\tilde{\mathbf{Y}}_{q,1}$,
respectively. Let $\mathbf{W}$ and $\mathbf{W}_{q}$ be zero-mean Gaussian
random vectors with dispersion matrices $\boldsymbol{\Gamma}$ and
$\boldsymbol{\Gamma}_{q}$, respectively. Then,
$\sqrt{n}\max_{q,i}|\tilde{b}_{2,q,i}|\stackrel{{\scriptstyle
d}}{{\rightarrow}}\|\mathbf{W}\|_{\infty}$ and
$\sqrt{n}\max_{i}|\tilde{b}_{2,q_{0},i}|\stackrel{{\scriptstyle
d}}{{\rightarrow}}\|\mathbf{W}_{q}\|_{\infty}$ as $n\to\infty$, respectively.

###### Proof.

Under the assumption of Gaussianity in the theorem, $\beta_{2,q,i}=q(q+2)$ for
all $q$. Since the $L_{\infty}$ norm in the Euclidean space is continuous, the
proof follows from an application of the multivariate central limit theorem in
subsection 6.2 and then applying the continuous mapping theorem. ∎

The tests of kurtosis are implemented under Gaussianity of the null
hypotheses. To compute the p-values for the MaxK test in (5.3) and the
$\text{MaxK}_{q_{0}}$ test in (5.4), we need to first estimate the dispersion
matrices $\boldsymbol{\Gamma}$ and $\boldsymbol{\Gamma}_{q}$, and for that, we
need to estimate the random vectors $\tilde{\mathbf{Y}}_{j}$ and
$\tilde{\mathbf{Y}}_{q,j}$. So, we first estimate $\mathbf{Y}_{qij}$ by:

$\displaystyle\widehat{\mathbf{Y}}_{qij}=\\{(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})\\}^{2}-2(q+2)(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi})^{\top}\widehat{\boldsymbol{\Sigma}}_{qi}^{-1}(\mathbf{X}_{qij}-\bar{\mathbf{X}}_{qi}),$

where $\widehat{\boldsymbol{\Sigma}}_{qi}$ is the sample dispersion matrix of
$\mathbf{X}_{qij}$. The estimates of $\tilde{\mathbf{Y}}_{j}$ and
$\tilde{\mathbf{Y}}_{qj}$ are constructed from these
$\widehat{\mathbf{Y}}_{qij}$s. Now, it can be derived that under Gaussianity,
$\text{Var}(\mathbf{Y}_{qi1})=8q(q+2)$ for all $q,i$. So, the diagonal entries
of the estimates of $\boldsymbol{\Gamma}$ and $\boldsymbol{\Gamma}_{q}$ are
fixed to be 1. Now, the off-diagonal entries are obtained from
$\text{Corr}(\mathbf{Y}_{q_{1},i_{1},1},\mathbf{Y}_{q_{2},i_{2},1})$. Let us
denote the estimated dispersion matrices as $\widehat{\boldsymbol{\Gamma}}$
and $\widehat{\boldsymbol{\Gamma}}_{q}$. Next, we generate 1000 independent
zero mean Gaussian random vectors
$\tilde{\mathbf{W}}_{1},\ldots,\tilde{\mathbf{W}}_{1000}$ with dispersion
matrix $\widehat{\boldsymbol{\Gamma}}$. The p-value of the MaxK test in (5.3)
is the proportion of the values
$\|\tilde{\mathbf{W}}_{1}\|_{\infty},\ldots,\|\tilde{\mathbf{W}}_{1000}\|_{\infty}$
larger than $\max_{q,i}|\tilde{b}_{2,q,i}|$. The p-value of the
$\text{MaxK}_{q_{0}}$ test in (5.4) is computed similarly.

###  6.3 Testing Gaussianity Based on Both Skewness and Kurtosis

Based on the tests of skewness and kurtosis, a test of Gaussianity can be
constructed. Analogous to (5.1) and (5.3), the null hypothesis here is:

$\displaystyle\text{H}_{0}^{(g)}:\text{The underlying distribution is
Gaussian}.$ (6.10)

Similarly, analogous to (5.2) and (5.4), the null hypothesis is:

$\displaystyle\text{H}_{0}^{q_{0},(g)}:\text{All $q_{0}$-dimensional subsets
of the data follow some Gaussian distribution}.$ (6.11)

Here, if at least one of the null hypotheses of the corresponding skewness
test or kurtosis test is rejected, then we reject the null hypothesis of
Gaussianity (after Bonferroni correction). For example, if any of the MaxS
test and the MaxK test rejects their null hypotheses, then (6.10) is also
rejected. The test for (6.10) is denoted as MaxSK test. Similarly, if any of
the $\text{MaxS}_{q_{0}}$ test or the $\text{MaxK}_{q_{0}}$ test rejects their
null hypothesis for a fixed $q_{0}$, then (6.11) is rejected, and we denote
this test by $\text{MaxSK}_{q_{0}}$ test.

##  7 Simulation Study

In this section, the performance of our proposed tests are investigated in
terms of the estimated sizes and powers using some simulated models. The
estimated powers of our tests are compared with the corresponding Mardia’s
tests. We also compare the estimated powers of ours tests with several tests
of Gaussianity.

###  7.1 Estimated Sizes

For estimating the sizes of our tests, we consider
$\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ being a random sample from the
$p$-variate Gaussian distribution ${\cal
N}_{p}(\mathbf{0},\boldsymbol{\Sigma})$, where
$\boldsymbol{\Sigma}=(\sigma_{ij})$ with $\sigma_{ij}=0.5+0.5\mathbb{I}(i=j)$.
We take $p=5$. The number of replicates to estimate the sizes of our tests is
taken as $1000$. The sample size $n$ is varied.

Table 2: Estimated sizes of the tests for nominal level 5% in ${\cal N}_{5}(\mathbf{0},\boldsymbol{\Sigma})$ based on 1000 replicates. Test | $n=100$ | $n=200$ | $n=500$ | $n=1000$
---|---|---|---|---
MaxS | 0.048 | 0.057 | 0.047 | 0.053
MaxK | 0.053 | 0.056 | 0.058 | 0.055
MaxSK | 0.056 | 0.070 | 0.052 | 0.056
$\text{MaxS}_{1}$ | 0.050 | 0.052 | 0.043 | 0.050
$\text{MaxS}_{2}$ | 0.050 | 0.053 | 0.041 | 0.050
$\text{MaxS}_{3}$ | 0.053 | 0.048 | 0.049 | 0.049
$\text{MaxS}_{4}$ | 0.047 | 0.054 | 0.055 | 0.051
$\text{MaxS}_{5}$ | 0.036 | 0.043 | 0.045 | 0.054
$\text{MaxK}_{1}$ | 0.057 | 0.049 | 0.043 | 0.060
$\text{MaxK}_{2}$ | 0.034 | 0.039 | 0.041 | 0.054
$\text{MaxK}_{3}$ | 0.027 | 0.027 | 0.046 | 0.047
$\text{MaxK}_{4}$ | 0.027 | 0.039 | 0.055 | 0.048
$\text{MaxK}_{5}$ | 0.081 | 0.065 | 0.051 | 0.057
$\text{MaxSK}_{1}$ | 0.053 | 0.062 | 0.046 | 0.055
$\text{MaxSK}_{2}$ | 0.041 | 0.055 | 0.054 | 0.046
$\text{MaxSK}_{3}$ | 0.040 | 0.034 | 0.047 | 0.046
$\text{MaxSK}_{4}$ | 0.036 | 0.037 | 0.057 | 0.052
$\text{MaxSK}_{5}$ | 0.049 | 0.057 | 0.038 | 0.056

From Table 2, it can be seen that the tests MaxS and MaxK have estimated sizes
close to the nominal level 5% for $n=100$ and up. The estimated sizes of
$\text{MaxS}_{5}$ and $\text{MaxK}_{2}$, $\text{MaxK}_{3}$, $\text{MaxK}_{4}$,
$\text{MaxK}_{5}$ and $\text{MaxSK}_{4}$ deviate slightly from the nominal
level 5% for $n=100$, and require higher sample sizes to converge to the
nominal level.

###  7.2 Estimated Powers

We now compare the powers of the tests with the usual Mardia skewness (MS) and
kurtosis (MK) tests along with several tests of normality.

The following tests of Gaussianity are considered for the comparison of
performances. The test by Henze and Zirkler [15] is denoted as the HZ test.
The test of normality developed by Royston [16, 17, 18, 19] is denoted as the
R test. The testing procedure described by Doornik and Hansen [20] is denoted
as the DH test. The skewness-based test of normality described by Kankainen,
Taskinen and Oja [21] is denoted as the $\text{KTO}_{\text{S}}$ test, and the
kurtosis-based test of normality described in the same paper is denoted as the
$\text{KTO}_{\text{K}}$ test. The test developed by Bowman and Shenton [22] is
denoted as the BS test. The testing procedure studied by Villasenor Alva and
Estrada [23] is denoted as the VE test. The test of normality developed by
Zhou and Shao [24] is denoted as the ZS test.

The HZ test, the R test and the DH test are implemented using the
corresponding functions in the R (R version 4.1.1 (2021-08-10), [25]) package
MVN [26]. The $\text{KTO}_{\text{S}}$ test and the $\text{KTO}_{\text{K}}$ are
implemented using the functions in the R package ICS [27]. The BS test and the
ZS test are implemented using their functions in the R package mvnormalTest
[28]. The VE test is implemented using its function in the R package
mvShapiroTest [29]. The Mardia’s skewness and kurtosis tests are implemented
based on the asymptotic distributions of the test statistics derived in [1]
using the unbiased estimate of the population covariance matrix.

Figure 2: Estimated powers of MaxS, MaxK (dashed line), MaxSK (dotted line),
MS, MK (dashed line), HZ, R, DH, $\text{KTO}_{\text{S}}$,
$\text{KTO}_{\text{K}}$ (dashed line), BS, VE and ZS tests for nominal level
5% (horizontal dashed line) in Model 1 (top), Model 2 (middle) and Model 3
(bottom) for $n=200$, $p=5$, $q=2$ based on 1000 replicates.

For the comparison of performances of the tests, we consider three simulation
models. In each of the models, the non-Gaussian feature is supported on a
small part of the data. Let $\boldsymbol{\Sigma}=(\sigma_{ij})$ with
$\sigma_{ij}=0.5+0.5\mathbb{I}(i=j)$ and $\boldsymbol{\alpha}=\alpha{\bf
1}_{q}$. The dimension of the matrix $\boldsymbol{\Sigma}$ is to be determined
based on the context. Let
$\mathbf{X}=(\mathbf{X}_{q}^{\top},\mathbf{X}_{p-q}^{\top})^{\top}$, where
$\mathbf{X}_{q}$ and $\mathbf{X}_{p-q}$ are independent and $\mathbf{X}_{p-q}$
follows ${\cal N}_{p-q}(\bf 0,\boldsymbol{\Sigma})$. Then:

  * •

Model 1: $\mathbf{X}_{q}\sim{\cal SN}_{q}(\bf
0,\boldsymbol{\Sigma},\boldsymbol{\alpha})$, where
$\boldsymbol{\alpha}=\alpha{\bf 1}_{q}$;

  * •

Model 2: $\mathbf{X}_{q}\sim t_{q}(\bf 0,\boldsymbol{\Sigma},\nu)$, where the
degree of freedom is $\nu$;

  * •

Model 3: $\mathbf{X}_{q}\sim{\cal ST}_{q}(\bf
0,\boldsymbol{\Sigma},\boldsymbol{\alpha},\nu)$, where
$\boldsymbol{\alpha}=\alpha{\bf 1}_{q}$ and $\alpha=1/\nu$.

It can be seen that the class of distributions in Model 1 is skewed Gaussian
distributions, while in Model 2, the class of distributions is symmetric heavy
tailed. In Model 3, the non-Gaussian distributions are both skewed and heavy
tailed.

We fix the sample size $n=200$, $p=5$ and $q=2$. Then we vary the values of
$\alpha$ or $\nu$ to investigate the changes in power of the tests in the
distributions. The plots of the estimated power curves of the tests are
presented in Figure 2. In the panel of Model 1, the kurtosis-based tests,
e.g., MaxK test, Madia’s kurtosis (MK) test and the $\text{KTO}_{\text{K}}$
test are not included, as Model 1 is concerned with skewness only. Similarly,
in the panel for Model 2, the skewness-based tests like the MaxS test,
Mardia’s skewness (MS) test and the $\text{KTO}_{\text{S}}$ test are not
included, as Model 2 is concerned with kurtosis only. However, in the plot for
Model 3, all the tests are included.

From Figure 2, it can be seen that in Model 1, the performances of MaxS and
MaxST tests are significantly better than all other tests. The power of the
MaxST test is slightly lower than the MaxS test. This is because MaxST
combines the MaxS and the MaxT tests using Bonferroni correction. If one of
the tests does not exhibit a high power, then the power of MaxST would be
lower than the best performing test. Similar observations can be made in the
panel for Model 2, where the estimated powers of the MaxK test and the MaxST
test are found to be better than all other tests, and the power curve of the
MaxSK test is slightly below the power curve of the MaxK test. In Model 3, the
estimated powers of the MaxSK test and the MaxK test is higher than other
testing procedures, while the power exhibited by the MaxS test is considerably
low in this model (this is because the skewness in Model 3 is rather weak). In
all the cases, it can be clearly seen that our tests significantly outperform
the other tests most of the time. The test for Gaussianity, i.e., the MaxSK
test, always performs quite well and exhibits better performance than all the
other tests.

###  7.3 Detection of Sub-Dimensions Supporting Skewness and Excess Kurtosis

The testing procedures described earlier can be used to detect the sub-
dimensions supporting non-Gaussian features in the data. Suppose the data are
skewed, but the skewness is supported on only a small sub-dimension of the
data. Then, to detect the sub-dimension supporting skewness, we can first
conduct the MaxS test. If the p-value of the MaxS test is small (say, lower
than 5%), then there is statistical evidence of presence of skewness in the
sample. Next, we find the sub-dimension corresponding to the maximum
$\tilde{b}_{1,q,i}$, which is the detected sub-dimension supporting skewness
in the data. Similarly, if a heavy-tailed component is present in a small sub-
dimension of the data, and we wish to detect that sub-dimension, we conduct
the MaxK test and if it rejects Gaussianity, we find the sub-dimension which
corresponds to maximum $|\tilde{b}_{2,q,i}|$. If we want to detect the sub-
dimension supporting a non-Gaussian distribution, we can use the MaxSK test in
the following way. We first conduct the MaxSK test. If it rejects Gaussianity,
then we find which p-value, whether for MaxS or MaxK, caused the rejection. If
it is only one of the tests, say MaxS, then we detect the sub-dimension
corresponding to maximum $\tilde{b}_{1,q,i}$. Otherwise, if the p-values of
both the MaxS test and the MaxK test are below 2.5% (due to Bonferroni
correction on the nominal level 5%), then we find the sub-dimensions
corresponding to maximum $\tilde{b}_{1,q,i}$ and maximum
$|\tilde{b}_{2,q,i}|$. The union of these two sub-dimensions is the detected
sub-dimension supporting the non-Gaussian distribution.

Table 3: Indices of all sub-dimensions for $p=5$. Sub-dimension | Index | Sub-dimension | Index | Sub-dimension | Index
---|---|---|---|---|---
(1) | 1 | (2) | 2 | (3) | 3
(4) | 4 | (5) | 5 | (1, 2) | 6
(1, 3) | 7 | (1, 4) | 8 | (1, 5) | 9
(2, 3) | 10 | (2, 4) | 11 | (2, 5) | 12
(3, 4) | 13 | (3, 5) | 14 | (4, 5) | 15
(1, 2, 3) | 16 | (1, 2, 4) | 17 | (1, 2, 5) | 18
(1, 3, 4) | 19 | (1, 3, 5) | 20 | (1, 4, 5) | 21
(2, 3, 4) | 22 | (2, 3, 5) | 23 | (2, 4, 5) | 24
(3, 4, 5) | 25 | (1, 2, 3, 4) | 26 | (1, 2, 3, 5) | 27
(1, 2, 4, 5) | 28 | (1, 3, 4, 5) | 29 | (2, 3, 4, 5) | 30
(1, 2, 3, 4, 5) | 31 |  |  |  |

To investigate the performance of the detection procedure described above, we
consider the three models described in subsection 7.2. In Model 1, we fix
$\alpha=5$ and conduct the detection procedure 1000 times on independent
replicates. The proportion of times each of the sub-dimensions is detected as
the one supporting the skewed distribution is computed. Also, the size of the
sub-dimensions (denoted as $q$) thus detected to support the skewed
distribution is also recorded, and the proportions for the $q$ values are
computed. All the possible sub-dimensions from $p=5$ variables are assigned
indices, which are presented in Table 3, and the proportions thus computed are
plotted against the indices in the first row of Figure 3. The estimated power
of the MaxS test there is 0.922. It can be clearly seen that the highest
proportions in the respective histograms are attained for the true sub-
dimension and the true $q$.

Figure 3: Histograms of the detection rate of sub-dimensions in Model 1 using
the centered and scaled skewness measure (first row), in Model 2 using the
centered and scaled kurtosis measure (second row), in Model 3 using the
procedure for detecting non-Gaussianity (third row). In all the cases,
$n=200$, $p=5$ and the true sub-dimension supporting skewness or excess
kurtosis is $q=2$. Results based on 1000 replicates.

Similarly, the procedure to detect the sub-dimension supporting the heavy
tailed distribution using the MaxK test is also carried out based on 1000
replicates in Model 2 with $\nu=5$. The estimated power of the MaxK test there
is 0.985. The histograms are presented in the second row of Figure 3, and we
again see that the highest proportions are attained for the true sub-dimension
and true $q$.

Finally, the procedure to detect the sub-dimension supporting a non-Gaussian
distribution using the MaxSK test is carried out in Model 3 taking
$\nu=1/\alpha=5$ and based on 1000 independent replicates. The estimated power
of the MaxSK test is 0.974 there. The histograms are presented in the third
row of Figure 3, where we again see that the highest proportions are attained
for the true sub-dimension and true value of $q$.

##  8 Sub-Dimensional Data Analysis

###  8.1 Fisher’s Iris Data

We revisit Fisher’s iris dataset discussed in section 1 where we considered
the part of these data related to the species ‘Iris setosa’ to demonstrate
that Mardia’s test of skewness fails to detect skewed features in sub-
dimensions. In Table 4, the p-values of the tests are presented. It can be
seen that our test can detect skewness in the data while Mardia’s test fails.
Next, we consider the whole Fisher’s iris dataset and compute the p-values of
all the tests. We find that our test of kurtosis can detect the deviation of
kurtosis from Gaussian kurtosis, while Mardia’s test of kurtosis fails.

Fisher’s iris data are generally modelled using a Gaussian distribution.
However, our findings point to the non-Gaussianity of the data, and thus it
may be judicious to use non-Gaussian and skewed distributions while analyzing
this dataset.

Table 4: Estimated p-values of the tests for the two data examples. | Iris setosa | Iris | Wind
---|---|---|---
Test name | $p=4$ | $p=4$ | $p=3$
MS | 0.236 | 0.000 | 0.004
MK | 0.448 | 0.611 | 0.338
MaxS | 0.004 | 0.000 | 0.026
MaxK | 0.349 | 0.006 | 0.007
MaxSK | 0.008 | 0.000 | 0.014
$\text{MaxS}_{1}$ | 0.001 | 0.327 | 0.022
$\text{MaxS}_{2}$ | 0.060 | 0.000 | 0.014
$\text{MaxS}_{3}$ | 0.102 | 0.000 | 0.002
$\text{MaxS}_{4}$ | 0.246 | 0.000 | –
$\text{MaxK}_{1}$ | 0.244 | 0.000 | 0.006
$\text{MaxK}_{2}$ | 0.229 | 0.076 | 0.126
$\text{MaxK}_{3}$ | 0.578 | 0.253 | 0.348
$\text{MaxK}_{4}$ | 0.438 | 0.621 | –
$\text{MaxSK}_{1}$ | 0.002 | 0.000 | 0.012
$\text{MaxSK}_{2}$ | 0.120 | 0.000 | 0.028
$\text{MaxSK}_{3}$ | 0.204 | 0.000 | 0.004
$\text{MaxSK}_{4}$ | 0.492 | 0.000 | –
HZ | 0.050 | 0.000 | 0.097
R | 0.000 | 0.000 | 0.011
DH | 0.000 | 0.000 | 0.022
$\text{KT}_{S}$ | 0.221 | 0.040 | 0.710
$\text{KT}_{K}$ | 0.875 | 0.046 | 0.349
BS | 0.060 | 0.000 | 0.020
VG | 0.012 | 0.000 | 0.057
ZS | 0.000 | 0.060 | 0.050

###  8.2 Wind Speed Data in Saudi Arabia

We consider a trivariate windspeed dataset produced by Yip [30] with the
Weather Research and Forecasting (WRF) model. The three components correspond
to bi-weekly mid-day windspeed during the period 2009-2014 at three locations
near Dumat Al Jandal, the first wind farm currently under construction in
Saudi Arabia. It is important to study the distributional properties of this
trivariate windspeed vector because they are crucial for understanding wind
patterns that will influence the production of electricity by the nearby wind
farm. In particular, it is of interest to assess whether a Gaussian
distribution is suitable, or a non-Gaussian model needs to be developed.

The dataset consists of $n=156$ trivariate windspeed vectors. A Ljung-Box test
reveals no indication of serial dependence, hence the data are treated as a
random sample from a three-dimensional distribution. The p-values of the
various tests are listed in Table 4. At the 5% level, Mardia’s tests support
skewness but do not reject a Gaussian kurtosis. Our global tests, however,
reject both symmetry and Gaussian kurtosis, suggesting that a non-Gaussian
distribution would be more suitable to model these data. Looking at our sub-
dimensional tests, we observe that skewness is rejected in all sub-dimensions,
whereas Gaussian kurtosis is rejected only for the $q=1$ dimensional
marginals. Among the other seven tests of normality, four of them reject
Gaussianity whereas the other three do not.

In summary, our new tests suggest to use a non-Gaussian distribution to model
these data and provide additional information about the non-Gaussian behavior
in sub-dimensional components of the trivariate distribution.

##  9 Discussion

We have developed some new tests of skewness and kurtosis which take into
account the skewness and excess kurtosis present in the sub-dimensions of the
data. It was demonstrated through analyses of simulated and real data that our
tests outperform the classical Mardia’s tests of skewness and kurtosis when
the skewness and the excess kurtosis are present in a small sub-dimension of
the whole data. Moreover, our tests can also be used as tests of Gaussianity,
and it is observed that as such, they outperform several popular tests of
Gaussianity. We have further developed a methodology to detect the true sub-
dimension when the skewness and the excess kurtosis are supported on a small
sub-dimension of the data.

One limitation of our methodology is that it considers all the possible sub-
dimensions, which is $2^{p}-1$, to detect skewness or excess kurtosis. The
number $2^{p}-1$ becomes large for even moderate values of $p$. So, the
methodology is computationally intensive. Future research needs to develop
suitable computational methods when the dimension $p$ of the multivariate data
is high to reduce the computation burden.

## Acknowledgements

This research was supported by the King Abdullah University of Science and
Technology (KAUST).

## References

  * [1] K. V. Mardia, Measures of multivariate skewness and kurtosis with applications, Biometrika 57 (3) (1970) 519–530.
  * [2] K. V. Mardia, Applications of some measures of multivariate skewness and kurtosis in testing normality and robustness studies, Sankhyā: The Indian Journal of Statistics, Series B 36 (2) (1974) 115–128.
  * [3] R. A. Fisher, The use of multiple measurements in taxonomic problems, Annals of Eugenics 7 (2) (1936) 179–188.
  * [4] A. Azzalini, A. Dalla Valle, The multivariate skew-normal distribution, Biometrika 83 (4) (1996) 715–726.
  * [5] A. Azzalini, A. Capitanio, Statistical applications of the multivariate skew normal distribution, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61 (3) (1999) 579–602.
  * [6] A. Azzalini, A. Capitanio, The Skew-Normal and Related Families, Vol. 3, Cambridge University Press, 2014.
  * [7] M. G. Genton, N. M. R. Loperfido, Generalized skew-elliptical distributions and their quadratic forms, Annals of the Institute of Statistical Mathematics 57 (2) (2005) 389–401.
  * [8] L. Baringhaus, N. Henze, Limit distributions for Mardia’s measure of multivariate skewness, The Annals of Statistics 20 (4) (1992) 1889–1902.
  * [9] J. L. Gastwirth, M. E. B. Owens, On classical tests of normality, Biometrika 64 (1) (1977) 135–139.
  * [10] G. G. Gregory, Large sample theory for $U$-statistics and tests of fit, The Annals of Statistics 5 (1) (1977) 110–123.
  * [11] A. Erdélyi, Higher Transcendental Functions, Vol. 2, McGraw-Hill Book Company, 1953.
  * [12] A. Higuchi, Symmetric tensor spherical harmonics on the $N$-sphere and their application to the de Sitter group $\text{SO}(N,1)$, Journal of Mathematical Physics 28 (7) (1987) 1553–1566.
  * [13] M. A. Blanco, M. Flórez, M. Bermejo, Evaluation of the rotation matrices in the basis of real spherical harmonics, Journal of Molecular Structure: THEOCHEM 419 (1–3) (1997) 19–27.
  * [14] N. Henze, On Mardia’s kurtosis test for multivariate normality, Communications in Statistics– Theory and Methods 23 (4) (1994) 1031–1045.
  * [15] N. Henze, B. Zirkler, A class of invariant consistent tests for multivariate normality, Communications in Statistics– Theory and Methods 19 (10) (1990) 3595–3617.
  * [16] J. P. Royston, An extension of Shapiro and Wilk’s $W$ test for normality to large samples, Journal of the Royal Statistical Society: Series C (Applied Statistics) 31 (2) (1982) 115–124.
  * [17] J. Royston, Some techniques for assessing multivarate normality based on the Shapiro-Wilk $W$, Journal of the Royal Statistical Society: Series C (Applied Statistics) 32 (2) (1983) 121–133.
  * [18] P. Royston, Approximating the Shapiro-Wilk $W$-test for non-normality, Statistics and Computing 2 (3) (1992) 117–119.
  * [19] P. Royston, Remark AS R94: A remark on algorithm AS 181: The $W$-test for normality, Journal of the Royal Statistical Society. Series C (Applied Statistics) 44 (4) (1995) 547–551.
  * [20] J. A. Doornik, H. Hansen, An omnibus test for univariate and multivariate normality, Oxford Bulletin of Economics and Statistics 70 (2008) 927–939.
  * [21] A. Kankainen, S. Taskinen, H. Oja, Tests of multinormality based on location vectors and scatter matrices, Statistical Methods and Applications 16 (3) (2007) 357–379.
  * [22] K. O. Bowman, L. R. Shenton, Omnibus test contours for departures from normality based on $\sqrt{b_{1}}$ and $b_{2}$, Biometrika 62 (2) (1975) 243–250.
  * [23] J. A. Villasenor Alva, E. González Estrada, A generalization of Shapiro–Wilk’s test for multivariate normality, Communications in Statistics– Theory and Methods 38 (11) (2009) 1870–1883.
  * [24] M. Zhou, Y. Shao, A powerful test for multivariate normality, Journal of Applied Statistics 41 (2) (2014) 351–363.
  * [25] R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria (2017).
URL https://www.R-project.org/

  * [26] S. Korkmaz, D. Goksuluk, G. Zararsiz, MVN: An R package for assessing multivariate normality., The R Journal 6 (2) (2014) 151–162.
URL https://journal.r-project.org/archive/2014-2/korkmaz-goksuluk-zararsiz.pdf

  * [27] K. Nordhausen, H. Oja, D. E. Tyler, Tools for exploring multivariate data: The package ICS, Journal of Statistical Software 28 (6) (2008) 1–31.
URL http://www.jstatsoft.org/v28/i06/

  * [28] Y. Zhang, M. Zhou, Y. Shao, mvnormalTest: Powerful Tests for Multivariate Normality, R package version 1.0.0 (2020).
URL https://CRAN.R-project.org/package=mvnormalTest

  * [29] E. Gonzalez-Estrada, J. A. Villasenor-Alva, mvShapiroTest: Generalized Shapiro-Wilk test for multivariate normality, R package version 1.0 (2013).
URL https://CRAN.R-project.org/package=mvShapiroTest

  * [30] C. M. A. Yip, Statistical characteristics and mapping of near-surface and elevated wind resources in the Middle East, Ph.D. thesis, KAUST, http://dx.doi.org/10.25781/KAUST-02YSW (2018).