Introduction

Background

The package MultiStatM provides general formulae for set partitions, multivariate moments and cumulants, vector Hermite polynomials. It provides theoretical formulae for some important symmetric and asymmetric multivariate distributions and well as estimation functions for multivariate moments and cumulants and connected measures of multivariate skewness and kurtosis.

The formulae implemented in the package can be found in the book “Multivariate Statistical Methods - Going Beyond the Linear”, Springer 2021 by Gy.Terdik and are fully general. For example, in the conversion formulae from multivariate moment to multivariate cumulants, given any list of (numerical) multivariate moments up to order \(k\), the conversion formula provides all multivariate cumulants up to order \(k\); this differs to a large degree from the formulae provided in the package kStatistics (Di Nardo and Guarino, 2022) which calculates one by one (individually) the cumulants of order \(r\) which are the entries of our cumulant vectors.

The packages MaxSkew and MultiSkew (Franceschini and Loperfido (2017a,b)) for detecting, measuring and removing multivariate skewness, computes the third multivariate cumulant of either the raw, centered or standardized data; s the main measures of multivariate skewness, together with their bootstrap distributions and provides orthogonal data projections with maximal skewness.

The package matrixcalc (Novomestky (2021)) provides the Commutation matrix, Elimination matrix, Duplication matrix for Cartesian tensor products of two vectors, which are particular cases of those provided in the package MultiStatM.

The package sn ( Azzalini (2022)) discusses for the skew-normal and the skew-t distributions, statistical methods are provided for data fitting and model diagnostics, in the univariate and the multivariate case. Random numbers generator for multivariate skew distributions are provided. In the package MultiStatM complete formulae for theoretical multivariate moments and cumulants of any order are implemented.

The package moments (Komsta and Novomestky (2022)) deals with functions to calculate moments, cumulants, Pearson’s kurtosis, Geary’s kurtosis and skewness; tests related to them from univariate data.

A careful study of the cumulants is a necessary and typical part of nonlinear statistics. Such a study of cumulants for multivariate distributions is made complicated by the index notations. One solution to this problem is the usage of tensor analysis. In this package (and the connected book) we offer an alternate method, which we believe is simpler to follow. The higher-order cumulants with the same degree for a multivariate vector can be collected together and kept as a vector. To be able to do so, we introduce a particular differential operator on a multivariate function, called the T -derivative, and use it to obtain cumulants and provide results which are somewhat analogous to well-known results in the univariate case.

More specifically, with the symbol \(\otimes\) denoting the Cartesian tensor product, consider the operator \(D_{\boldsymbol{\lambda}}^{\otimes}\), which we refer to as the \(\operatorname{T}\)-derivative; see Jammalamadaka et al. (2006) for details. For any function \(\boldsymbol{\phi}(\boldsymbol{\lambda})\), the~\(\operatorname{T}\)-derivative is defined as \[\begin{equation}\label{Tderiv} D_{\boldsymbol{\lambda}}^{\otimes}\boldsymbol{\boldsymbol{\phi}}% (\boldsymbol{\lambda})=\operatorname{vec}\left(\left( \frac{\partial\boldsymbol{\phi }(\boldsymbol{\lambda})}{\partial\boldsymbol{\lambda}^{\top}}\right) ^{\top }\right)=\boldsymbol{\phi}(\boldsymbol{\lambda})\otimes\frac{\partial}{\partial \boldsymbol{\lambda}}.% \end{equation}\] \({\boldsymbol{\phi}}\) is \(k\)-times differentiable, with~its \(k\)-th \(\operatorname{T}\)-derivative \(D_{\boldsymbol{\lambda}}^{\otimes k}\boldsymbol{\boldsymbol{\phi} }(\boldsymbol{\lambda})=D_{\boldsymbol{\lambda}}^{\otimes}\left( D_{\boldsymbol{\lambda}}^{\otimes k-1}\boldsymbol{\boldsymbol{\phi} }(\boldsymbol{\lambda})\right)\).

In the following we demonstrate the use of this technique through the characterization of several multivariate distributions via their cumulants and by extending the discussion to statistical inference for multivariate skewness and kurtosis.

We note that Kollo (2006) provides formulae for cumulants in terms of matrices; however, retaining a matrix structure for all higher-order cumulants leads to high-dimensional matrices with special symmetric structures which are quite hard to follow notionally and computationally. McCullagh (2018) provides quite an elegant approach using tensor methods; however, tensor methods are not very well known and computationally not so~simple.

The method discussed here is based on relatively simple calculus. Although the tensor product of Euclidean vectors is not commutative, it has the advantage of permutation equivalence and allows one to obtain general results for cumulants and moments of any order, as it will be demonstrated in this paper, where general formulae, suitable for algorithmic implementation through a computer software, will be provided.

Methods based on a matrix approach do not provide this type of result; see also (Ould-Baba (2015), which goes as far as the sixth-order moment matrices, whereas there is no such limitation in our derivations and our results. For further discussion, one can see also Kolda (2009) and Qi (2006).

Set Partitions

The package MultiStatM provides several functions dealing with set partitions. Such functions provide some basic tools used to built the multivariate formulae for moments and cumulants in the following sections.

Generally a set of \(N\) elements can be split into a set of disjoint subsets, i.e. it can be partitioned. The set of \(N\) elements will correspond to set \(1 : N = \{1, 2, \dots ,N\}\). If \(\cal{K} = \{b_1, b_2, \dots , b_r \}\) where each \(b_j \subset 1 : N\), then \(\cal{K}\) is a partition provided \(\cup b_j = 1 : n\), each \(b_j\) is non-empty and \(b_j \cap b_i = \emptyset\) (the empty set) is disjoint whenever \(j \neq i\). The subsets \(b_j\), \(j = 1, 2, \dots, r\) are called the blocks of \(\cal{K}\). We will call \(r\) (the number of the blocks in partition \(\cal{K}\)), the size of \(\cal{K}\), and denote it by \(|\cal{K}| = r\), and a partition with size \(r\) will be denoted by \(\cal{K}_{\{r\}}\). Let us denote the set of all partitions of the numbers \(1 : N\) by \(\cal{P}_N\).

Consider next a partition \(\cal{K}_{\{r\}} = \{b_1, b_2, \dots , b_r \} \in \cal{P}_N\), with size \(r\). Denote the cardinality \(k_j\) of a block in the partition \(\cal{K}_{\{r\}}\), i.e. \(k_j =|b_j|\). The type of a partition \(\cal{K}_{\{r\}}\) is \(l = [ l_1, \dots , l_N]\), if \(\cal{K}_{\{r\}}\) contains exactly \(l_j\) blocks with cardinality \(j\). A partition with size \(r\) and type \(l\) will be denoted by \(\cal{K}_{\{r|l\}}\). It is clear that \(l_j ≥ 0\), and \(\sum_j jl_j=N\), and \(\sum_j l_j=r\). Naturally, some \(l_j\) ’s are zero.

The basic function is Partition_Type_All which provides complete information on the partition of a set of N elements, namely

• Part.class: the list of all possible partitions given as partition matrices

• S_N_r: a vector with the number of partitions of size r=1, r=2, etc. (Stirling numbers of second kind )

• eL_r: a list of partition types with respect to partitions of size r=1, r=2, etc.

• S_r_j: vectors of number of partitions with given types grouped by partitions of size r=1, r=2, etc.

Example 1. Consider the case where N=4 and run the following

PTA<-Partition_Type_All(4)

All the partition matrices are listed in Part.class. See the first three below

PTA$Part.class[[1]]
#> [1] 1 1 1 1
PTA$Part.class[[2]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    1    1    0
#> [2,]    0    0    0    1
PTA$Part.class[[3]] 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    1    0    1
#> [2,]    0    0    1    0

The lists S_N_r and eL_r contain respectively the number of partitions with r=1 blocks, r=2 blocks, etc. and the types of partitions for each partition of size r

PTA$S_N_r 
#> [1] 1 7 6 1
PTA$eL_r
#> [[1]]
#> [1] 0 0 0 1
#> 
#> [[2]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    1    0
#> [2,]    0    2    0    0
#> 
#> [[3]]
#> [1] 2 1 0 0
#> 
#> [[4]]
#> [1] 4 0 0 0

From the results above we see that there are 1 partition of 1 block (r=1), 7 partitions of two blocks (r=2), 6 partitions of 3 blocks and 1 partition of 4 blocks. From PTA$eL_r[[2]] we see that there are two types of partition with r=2: the first is of type \((l_1=1,l_3=1)\) and the second is of type \((l_2=2)\).

If one needs to know how many partitions of type \((l_1,l_3)\) are there for r=2 then consider the list S_r_j, i.e.

PTA$S_r_j[[2]] ## Partitions with size r=2, includes two types (above) each with number
#> [1] 4 3

That is, for partitions with r=2 blocks, 4 are of type \((l_1=1,l_3=1)\) and 3 are of type \((l_2=2)\).

Other available partition functions are:

• Partition_2perm which provides the permutation of N elements according to a partition matrix \(pU\);

• Partition_Indecomposable, which provides the list of all indecomposable partitions with respect to a partition matrix L

• Partition_Pairs, which provides the list of partitions dividing into pairs the set of N elements.

• Partition_DiagramsClosedNoLoops Closed Diagrams without Loops

• Permutation_Inverse which provides the inverse of a permutation

For further details see Terdik (2021, chapter 1).

Commutators, symmetrizer and selection matrices

The matr group of functions produce commutators and selection matrices. The use of matrices allows represent as linear combinations problems of permutation and powers of T-products. On the other side, the size of these matrix can quickly become quite important. To deal with this issues and option for sparse matrices is always provided; also a corresponding indx group of functions is provided; these function provide selection vectors which give equivalent results ans the corresponding functions in the group matr.

The function matr_Commutator_Kmn produces a commutation matrix, \(\mathbf{K}_{m \cdot n}\) of dimension \(mn \times mn\) such that, given a matrix \(\mathbf{A}\) \(m\times n\), \(\mathbf{K}_{m \cdot n} \operatorname{vec}\mathbf{A}=\operatorname{vec}\mathbf{A}'\) (see Terdik (2021, p.8,(1.12)) ). The same result can be obtained using indx_Commutator_Kmn

matr_Commutator_Kperm and indx_Commutator_Kperm produce any permutation of kronecker products of vectors of any length.

Example 2. For the product of vectors \(\mathbf{a}_1 \otimes \mathbf{a}_2 \otimes\mathbf{a}_3\) of dimensions \(d_1\) to \(d_3\) respectively. matr_Commutator_Kperm(c(3,1,2),c(d1,d2,d3)) produces \(\mathbf{a}_3 \otimes \mathbf{a}_1 \otimes\mathbf{a}_2\).

a1<-c(1,2)
a2<-c(2,3,4)
a3<-c(1,3)
p1<-a1%x%a2%x%a3
as.vector(matr_Commutator_Kperm(c(3,1,2),c(2,3,2))%*%p1) ## same result as below
#>  [1]  2  3  4  4  6  8  6  9 12 12 18 24
a3%x%a1%x%a2   
#>  [1]  2  3  4  4  6  8  6  9 12 12 18 24

The same result can be obtained by using indx_Commutator_Kperm

p1[indx_Commutator_Kperm(c(3,1,2),c(2,3,2))]
#>  [1]  2  3  4  4  6  8  6  9 12 12 18 24

The matr_Commutator_Mixing is exploited for deriving the covariance matrix of Hermite polynomials; see Terdik (2021, 4.6).

The Elimination and Qplication matrices- related functions respectively eliminate and restore duplicated or q-plicated elements in powers of T-products.

a<-c(1,2)
a3<-a%x%a%x%a
a3
#> [1] 1 2 2 4 2 4 4 8
as.vector(matr_Elimination(2,3)%*%a3)
#> [1] 1 2 4 8
as.vector(matr_Qplication(2,3)%*%matr_Elimination(2,3)%*%a3)
#> [1] 1 2 2 4 2 4 4 8

Closely connected to the above matrices are the functions indx_UnivMomCum and indx_Elimination. The former provides a vector of indexes to select univariate moments or cumulants of the single elements of a d-vector X from available vector of T-moments and T-cumulants. The latter eliminates the duplicated/q-plicated elements in a T-vector of multivariate moments and cumulants. The function indx_Elimination produces the same results as matr_Elimination and it is less demanding in terms of memory. The use of matr_Elimination can be preferable is one wishes to deal with linear combination of matrices. See examples 4 and 6 below for the use of indx_UnivMomCum and indx_Elimination.

The symmetrizer matrix, a \(d^n \times d^n\) matrix for the symmetrization of a T-product of \(n\) vectors with the same dimension \(d\) which overcomes the difficulties arising from the non commutative property of the Kronecker product, and simplifies considerably the computation formulae for multivariate polynomials and their derivatives (see Holmquist (1996) for details). The symmetrizer for a T-product of \(q\) vectors of dimension \(d\) is defined as \[ \mathbf{S}_{d \mathbf{1}q}=\frac{1}{q} \sum_{p \in \cal{P}_q} \mathbf{K}_p \] where \(\cal{P}_q\) is the set of all permutations of the numbers \(1:q\) and \(\mathbf{K}_p\) is the commutator matrix of for the permutation \(p \in \cal{P}_q\), (i.e. the matr_Commutator_Kperm of the package). Note that, by definition, computing the symmetrizer requires \(q!\) operations; in the package, the computational complexity is overcome by exploiting the Chacon and Duong (2015) efficient recursive algorithms for functionals based on higher order derivatives.

Multivariate T-Hermite Polynomials

Consider a Gaussian vector \(\mathbf{X}\) of dimension \(d\) with \(\operatorname{E}\mathbf{X}\) and \(\mathbf{\Sigma}=\operatorname{Cov}(\mathbf{X})=\operatorname{E}\mathbf{X X}'\) and define the generator function \[\begin{split} \Psi(\mathbf{X}; \mathbf{a})&=\exp \left(\mathbf{a}'\mathbf{X} - \frac{1}{2} \mathbf{a}' \mathbf{\Sigma} \mathbf{a}\right) \\ &=\exp \left(\mathbf{a}'\mathbf{X} - \frac{1}{2} \boldsymbol{\kappa}_2^{\otimes\prime} \mathbf{a}^{\otimes 2} \right) \\ \end{split} \] where \(\mathbf{a}\) is a \(d\)-vector of constants and \(\boldsymbol{\kappa}_2^{\otimes}=\operatorname{vec}\mathbf{\Sigma}\). The vector Hermite polynomials is defined via the T-derivative of the generator function, viz. \[ \mathbf{H}_n(\mathbf{X}) = D_\mathbf{a}^{\otimes n}\Psi(\mathbf{X};\mathbf{a})\big|_{\mathbf{a}=0} \] For example one has \[ \mathbf{H}_1(\mathbf{X})=\mathbf{X}, \quad \mathbf{H}_2(\mathbf{X})=\mathbf{X}^{\otimes 2} - \boldsymbol{\kappa}_2^{\otimes} \] Note that the multivariate T-Hermite polynomial \(\mathbf{H}_n(\mathbf{X})\) is a vector of dimension \(d^n\) which contains the \(n\)-th order polynomials of the vector \(\mathbf{X}^{\otimes n}\). For example the entries of \(\mathbf{H}_2(\mathbf{X})\) are the second order Hermite polynomials \(H_2(X_i,X_j)\), \(i,j=1,2, \dots d\); for \(d=2\) \[ \mathbf{H}_2(\mathbf{X}) = \left( (X_1^2 - \sigma_{11}), (X_1 X_2 - \sigma_{12}), (X_2 X_1 - \sigma_{21}), (X_2^2 - \sigma_{22})\right)^\prime. \] Note that \(\mathbf{H}_n(\mathbf{X})\) is \(n\)-symmetric, i.e. \(\mathbf{H}_2(\mathbf{X}) = \mathbf{S}_{d \mathbf{1}_n} \mathbf{H}_2(\mathbf{X})\) where \(\mathbf{S}_{d \mathbf{1}_n}\) is the symmetrizer defined in […]. From this one can get useful recursion formulae \[ \mathbf{H}_n(\mathbf{X})=\mathbf{S}_{d \mathbf{1}_n}\left( \mathbf{H}_{n-1}(\mathbf{X}) \otimes \mathbf{X}- (n-1) \mathbf{H}_{n-2}(\mathbf{X}) \otimes \boldsymbol{\kappa}_2^{\otimes} \right). \]

For further details, consult Terdik (2021, 4.5).

The definition of the \(d\)-variate Hermite polynomial requires the covariance matrix \(\mathbf{\Sigma}\) of the vector \(\mathbf{X}\). The Hermite_Poly_HN and Hermite_Poly_HN_multi functions compute the univariate and d-variate Hermite polynomials and their inverses up to a given order N evaluated at \(x\) for a given covariance matrix Sig2. By default Sig2=\(I_\mathbf{d}\).

Example 3 The first and the second \(3\)-variate Hermite polynomials evaluated at x<-c(1,2,3) where \(x\) is the realization of \(\mathbf{X} \sim N(\mathbf{0}, I_{\mathbf{3}})\) is

x<-c(1,2,3)
H2<-Hermite_Poly_HN_Multi(x,N=2)
H2[[1]]
#> [1] 1 2 3
H2[[2]]
#> [1] 0 2 3 2 3 6 3 6 8

If x is the realization of \(\mathbf{X} \sim N(\mathbf{0}, 4I_\mathbf{2})\)

H2<-Hermite_Poly_HN_Multi(x,Sig2=4*diag(3),N=2)
H2[[1]]
#> [1] 1 2 3
H2[[2]]
#> [1] -3  2  3  2  0  6  3  6  5

One can recover the vector x from H2 with the inverse function:

Hermite_Poly_NH_Multi_Inv(H2,N=2,Sig2=4*diag(3))[[1]] 
#> [1] 1 2 3

The function Hermite_N_Cov_X1_X2 can be exploited to obtain the covariance matrix of \(H_N(\mathbf{X}_1)\) and \(H_N(\mathbf{X}_2)\) for vectors \(\mathbf{X}_1\) and \(\mathbf{X}_2\) having covariance matrix \(\mathbf{\Sigma_{12}}\).

Covmat<-matrix(c(1,0.8,0.3,0.8,2,1,0.3,1,2),3,3)
Cov_X1_X2 <- Hermite_N_Cov_X1_X2(SigX12=Covmat,N=3)

Multivariate T-moments and T-cumulants

Multivariate moments and cumulants of all orders of a random vector \(\mathbf{X}\) in \(d\)-dimensions, with mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\mathbf{\Sigma}\) can be obtained by applying the T-derivative respectively to the characteristic function and the log of the CF.

More formally, let \(\boldsymbol{\lambda}\) a \(d\)-vector of real constants; \(\phi_{\mathbf{X}}(\boldsymbol{\lambda})\) and \(\psi_{\mathbf{X}}(\boldsymbol{\lambda})=\log\phi_{\mathbf{X}}(\boldsymbol{\lambda})\) denote, respectively, the characteristic function and the cumulant- function of \(\mathbf{X}\).

Then the \(k\)-th order moments and cumulants of the vector \(\mathbf{X}\) are obtained as \[ \boldsymbol{\mu}^\otimes_{\mathbf{X},k} = (-\mathbf{i})^k D_{\boldsymbol{\lambda}}^{\otimes k}\boldsymbol{\boldsymbol{\phi}% }_{\mathbf{X}}(\boldsymbol{\lambda}) \big|_{\boldsymbol{\lambda}=0}. \] \[ \boldsymbol{\kappa}^\otimes_{\mathbf{X},k} = \underline{\operatorname{Cum}}_k(\mathbf{X})= (-\mathbf{i})^k D_{\boldsymbol{\lambda}}^{\otimes k}\boldsymbol{\boldsymbol{\psi}% }_{\mathbf{X}}(\boldsymbol{\lambda}) \big|_{\boldsymbol{\lambda}=0}. \] Note that \(\boldsymbol{\mu}_{\mathbf{X},k} = \operatorname{E}\mathbf{X}^{\otimes k}\) that is a vector of dimension \(d^k\) that contains all possible moments of order order \(k\) formed by \(X_1, \dots, X_d\). This approach has the advantage of being straightforwardly extendable to any \(k\)-th order moment. An analogous discussion can be done for cumulants.

Note that one has \(\boldsymbol{\kappa}_{\mathbf{X},2} =\operatorname{vec} \mathbf{\Sigma}\).

The package MultiStatM contains functions which obtains moments from cumulants and vice-versa as well as function which provide theoretical moments and cumulants for some important multivariate distributions.

The conv_Cum2Mom and conv_Mom2Cum either for the univariate and multivariate cases provide conversion formulae for cumlants from moments and viceversa given any list of (theoretical) moments (or cumulants).

The conversion formula from moments to cumulants (see Terdik (2001, 3.4)) is given by \[\begin{split} \boldsymbol{\mu}_{n}^\otimes &= \sum_{\cal{K} \in \cal{P}_n} \mathbf{K}^{-1}_{p(\cal{K})} \prod^\otimes_{b_j \in \cal{K}} \kappa^\otimes_{|b_j|}\\ &= \mathbf{S}_{d \mathbf{1}_n}\left( \sum_{r=1}^n \sum_{\sum l_j =r, \sum j l_j = n} \frac{n!}{\prod_{j=1}^n l_j! (j!)^{l_j}} \prod_{j=1:n-r+1}^\otimes \kappa^{\otimes l_j}_j\right) \end{split} \] where the summation is over all partitions \(\cal{K} = \{b_1, b_2,\dots, b_k\}\) of \(1 : n\); \(|b_j|\) denotes the cardinality of block \(b_j\). The simpler second formula, exploiting the symmetrizer matrix, derives from symmetry of \(\boldsymbol{\mu}_{n}^\otimes\).

As far as the formula from cumulants to moments (Terdik (2021, 3.4)) is concerned, \[ \boldsymbol{\kappa}_{n}^\otimes = \mathbf{S}_{d \mathbf{1}_n}\left( \sum_{r=1}^n (-1)^{r-1} (r-1)!\sum_{\sum l_j =r, \sum j l_j = n} \prod_{j=1:n-r+1}^\otimes \frac{1}{{l_j}!}\left( \frac{1}{j!}\boldsymbol{\mu}^{\otimes}_j\right)^{l_j}\right) \]

Example 4. Consider the case of the 2-variate standard normal distribution with null mean vector and covariance matrix with unit elements on the main diagonal and off-diagonal elements equal to 0.5; in this case the the first four moments are given in the vector mu below

mu<-list(c(1,1),c(2,1.5,1.5,2),c(4,3,3,3,3,3,3,4),c(10,7,7,6.5,7,6.5,6.5,7,7,6.5,6.5,7,6.5,7,7,10))
cum<-conv_Mom2CumMulti(mu)
cum
#> [[1]]
#> [1] 1 1
#> 
#> [[2]]
#> [1] 1.0 0.5 0.5 1.0
#> 
#> [[3]]
#> [1] 0 0 0 0 0 0 0 0
#> 
#> [[4]]
#>  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Getting back to moments

conv_Cum2MomMulti(cum)
#> [[1]]
#> [1] 1 1
#> 
#> [[2]]
#> [1] 2.0 1.5 1.5 2.0
#> 
#> [[3]]
#> [1] 4 3 3 3 3 3 3 4
#> 
#> [[4]]
#>  [1] 10.0  7.0  7.0  6.5  7.0  6.5  6.5  7.0  7.0  6.5  6.5  7.0  6.5  7.0  7.0
#> [16] 10.0

I one wishes to select only the distinct moments from the vector of third moments, then

mu[[3]][indx_Elimination(2,3)]
#> [1] 4 3 3 4

Alternatively one can also use the elimination matrix

r.mu<-matr_Elimination(2,3)%*% mu[[3]]
as.vector(r.mu)
#> [1] 4 3 3 4

Note that matr_Elimination does not correspond the the function unique, rather it individuates the duplicated elements from the symmetry of the Kronecker product. This allow to recover the whole vector when needed.

as.vector(matr_Qplication(2,3)%*%r.mu)
#> [1] 4 3 3 3 3 3 3 4

The same result by using indx_Qplication

indx_Qplication(r.mu,2,3)
#> [1] 4 3 3 3 3 3 3 4

The MomCum functions provide theoretical moments and cumulants for some common multivariate distributions: Skew-normal. Canonical Fundamental Skew-normal (CFUSN), Uniform distribution on the sphere, central folded Normal distribution (univariate and multivariate); for detail on the multivariate formulae used see (Jammalamadaka et al., 2021).

Expressions for theoretical moments and cumulants are provided by the distr_NAME_MomCum_Th functions. In particular

• A \(d\)-vector \(\mathbf{U}\) having uniform distribution on the sphere \(\mathbb{S}_{d-1}\). Moments and cumulants of all orders are provided for \(\mathbf{U}\) (NAME=Uni) and its modulus (NAME=UniAbs). Recall that any \(d\)-vector, say \(\mathbf{W}\) has a spherically symmetric distribution if that distribution is invariant under the group of rotations in \(\mathbb{R}^{d}\). This is equivalent to saying that \(\mathbf{W}\) has the stochastic representation \(\mathbf{W}=R\mathbf{U}\) where \(R\) is a non negative random variable. Moments and cumulants of \(\mathbf{W}\) can be obtained by its stochastic representation as discussed in Jammalamadaka et al.(2021a, Theorem 1) and Jammalamadaka et al. (2021c, Lemma 1). Furthermore a \(d\)-vector \(\mathbf{X}\) has an elliptically symmetric distribution if it has the representation \[ \mathbf{X}=\boldsymbol{\mu}+\boldsymbol{\Sigma}^{1/2}\mathbf{W}% \] where \(\boldsymbol{\mu}\in\mathbb{R}^{d}\), \(\boldsymbol{\Sigma}\) is a variance-covariance matrix and \(\mathbf{W}\) is spherically distributed. Hence the cumulants of \(\mathbf{X}\) are just constant times the cumulants of \(\mathbf{W}\) except for the mean i.e. \[ \underline{\operatorname*{Cum}}_{m}\left( \mathbf{X}\right) =\left( \boldsymbol{\Sigma}^{1/2}\right) ^{\otimes m}\underline{\operatorname*{Cum}% }_{m}\left( \mathbf{W}\right) , \]

• If \(\mathbf{Z}\) denotes a \(d\)-vector with \(d\)-variate Normal distribution, the functions with NAME=Zabs and NAME=ZabsM provide the moments and cumulants of \(|\mathbf{Z}|\) respectively in the univariate (\(d=1\)) and multivariate case.

• The Multivariate Skew Normal distribution introduced by Azzalini and Dalla Valle (1996), whose marginal densities are scalar skew-normals. A \(d\)-dimensional random vector \(\mathbf{X}\) is said to have a multivariate skew-normal distribution, \(\text{SN}_{d}\left(\boldsymbol{\mu},\boldsymbol{\Omega},\boldsymbol{\alpha}\right)\) with shape parameter \(\boldsymbol{\alpha}\) if it has the density function \[ 2\varphi\left( \mathbf{X};\boldsymbol{\mu},\boldsymbol{\Omega}\right) \Phi\left( \boldsymbol{\alpha}^{\top}\left( \mathbf{X}-\boldsymbol{\mu }\right) \right) , \quad\mathbf{X} \in\mathbb{R}^{d}. \] where \(\varphi\left(\mathbf{X};\boldsymbol{\mu},\boldsymbol{\Omega}\right)\) is the \(d\)-dimensional normal density with mean \(\boldsymbol{\mu}\) and correlation matrix \(\boldsymbol{\Omega}\); here \(\varphi\) and \(\Phi\) denote the univariate standard normal density and the cdf. For a general formula for cumulants, see Jammalamadaka et al. (2021a, Lemma 4).

• Arellano-Valle and Genton (2005) introduced the CFUSN distribution (cf. their Proposition 2.3), to include all existing definitions of SN distributions. The marginal stochastic representation of \(\mathbf{X}\) with distribution \(\text{CFUSN}_{d,m}\left(\boldsymbol{\Delta}\right)\) is given by \[ \mathbf{X}=\boldsymbol{\Delta}\left\vert \mathbf{Z}_{1}\right\vert +\left( \mathbf{I}_{d}-\boldsymbol{\Delta\Delta}^{\top}\right) ^{1/2}\mathbf{Z}_{2} \] where \(\boldsymbol{\Delta}\), is the \(d\times m\) skewness matrix such that \(\left\Vert \boldsymbol{\Delta}\underline{a}\right\Vert <1\), for all \(\left\Vert \underline{a}\right\Vert =1\), and \(\mathbf{Z}_{1}\in\mathcal{N}\left( 0,\mathbf{I}_{m}\right)\) and \(\mathbf{Z}_{2}\in\mathcal{N}\left( 0,\mathbf{I}_{d}\right)\) are independent (Proposition 2.2. Arellano-Valle and Genton (2005)). A simple construction of \(\boldsymbol{\Delta}\) is \(\boldsymbol{\Delta}=\boldsymbol{\Lambda}\left(\mathbf{I}_{m}\mathbf{+}\boldsymbol{\Lambda}^{\top}\boldsymbol{\Lambda}\right)^{-1/2}\) with some real matrix \(\boldsymbol{\Lambda}\) with \(d\times m\). The \(\text{CFUSN}_{d,m}\left(\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\Delta}\right)\) can be defined via the linear transformation \(\boldsymbol{\mu}+\boldsymbol{\Sigma}^{1/2}\mathbf{X}\). For a general formula for cumulants, see Jammalamadaka et al. (2021a, Lemma 5).

The Rand functions provide random number generators for multivariate distributions.

The EVSK functions compute the theoretical values of the mean vector, covariance, skewness vector, total skenwness, kurtosis vector and total kurtosis for given multivariate distributions: Uniform on the sphere, modulus of the Uniform distribution on the sphere, Skew-normal distribution.

Example 5. For a skew-normal distribution with \(\alpha=(10,5,0)\) and correlation function \(\Omega= \text{diag} (1,1,1)\) we have the third moments and cumulants are

alpha<-c(10,5,0)
omega<-diag(3)
MSN<-distr_SkewNorm_MomCum_Th(r=3,omega,alpha,nMu=TRUE)
round(MSN$Mu[[3]],3)
#>  [1] 1.568 0.073 0.000 0.073 0.570 0.000 0.000 0.000 0.711 0.073 0.570 0.000
#> [13] 0.570 0.996 0.000 0.000 0.000 0.355 0.000 0.000 0.711 0.000 0.000 0.355
#> [25] 0.711 0.355 0.000
round(MSN$CumX[[3]],3)
#>  [1] 0.154 0.077 0.000 0.077 0.039 0.000 0.000 0.000 0.000 0.077 0.039 0.000
#> [13] 0.039 0.019 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
#> [25] 0.000 0.000 0.000

As another example, for the modulus of the Uniform distribution on the sphere, the fourth cumulant is:

distr_Uni_EVSK_Th(3, nCum = TRUE)$Kurt.U
#>  [1] -1.2  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0 -0.4  0.0 -0.4  0.0  0.0
#> [16]  0.0  0.0  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0
#> [31] -0.4  0.0  0.0  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -1.2  0.0  0.0  0.0 -0.4
#> [46]  0.0  0.0  0.0  0.0  0.0 -0.4  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0  0.0  0.0
#> [61] -0.4  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -0.4  0.0 -0.4  0.0 -0.4  0.0  0.0
#> [76]  0.0 -0.4  0.0  0.0  0.0 -1.2

Estimation

Estimating functions starting from a vector of multivariate data are available: multivariate moments and cumulants, skewness and kurtosis vectors Mardia’s skewness and kurtosis indexes, Mori, Rohatgi, Szekely (MRSz’s) skewness vector and kurtosis matrices.

A complete picture of skewness is provided by the third-order T-cumulant (skewness vector) of a standardized \(\mathbf{X}\); set \(\mathbf{Y}=\mathbf{\Sigma}^{-1/2}(\mathbf{X}-\boldsymbol{\mu})\), then the skewness vector is \[ \boldsymbol{\kappa}_{\mathbf{Y},3}^\otimes =\underline{\operatorname{Cum}}_3 \left( \mathbf{Y}\right)=\left(\mathbf{\Sigma}^{-1/2}\right)^{\otimes 3} \boldsymbol{\kappa}_{\mathbf{X},3}^\otimes. \] The total skewness of \(\mathbf{X}\) is defined by the square norm of the skewness vector: \(\gamma_{1,d}=\|\boldsymbol{\kappa}_{\mathbf{Y},3}^\otimes\|^2\). This definition guarantees that skewness is invariant under the shifting and orthogonal transformations, in other words it is affine invariant.

We note that Mardia’s multivariate skewness index (Mardia (1970)), denote it by \(\beta_{1,d}\), coincides with the total skewness \(\gamma_{1,d}\) since the third-order central moments and third-order cumulants are equal.

Mori, Rohatgi, Szekely (MRSz’s) skewness vector (Mori et al. (1994)) can also be recovered from the skewness vector as \[ \mathbf{b}(\mathbf{Y})= \left( \operatorname{vec}' \mathbf{I}_d \otimes \mathbf{I}_d \right)\boldsymbol{\kappa}_{\mathbf{Y},3}^\otimes \] Note that \(\operatorname{vec}' \mathbf{I}_d \otimes \mathbf{I}_d\) is a matrix of dimension \(d \times d^3\), which contains \(d\) unit values per-row, whereas all the others are 0; as a consequence, this measure does not take into account the contribution of cumulants of the type \(\operatorname{Cum}_3 (X_j,X_k,X_l)\), where all the three indices \(j\), \(k\), \(l\) are different from each other. The corresponding scalar measure of multivariate skewness is \(b(\mathbf{Y}) = \| \mathbf{b}(\mathbf{Y}) \|^2\).

The fourth-order T-cumulant of the standardized \(\mathbf{X}\), i.e. \(\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes\), will be called kurtosis vector of \(\mathbf{X}\); its square norm will be called the total kurtosis of \(\mathbf{X}\) \[ \gamma_{2,d}=\| \boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes \|^2 \] Mardia’s kurtosis index \(\beta_{2,d}= \operatorname{E}\left( \mathbf{Y}'\mathbf{Y} \right)^2\) is related to the kurtosis vector by the formula \[ \beta_{2,d}= \left( \operatorname{vec}' \mathbf{I}_{d^2} \right)\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes +d(d+2) \] A consequence of this is that Mardia’s measure does not depend on all the entries of \(\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes\) which has \(d(d +1)(d +2)(d +3)/24\) distinct elements, while \(\beta_{2,d}\) includes only \(d^2\) elements among them. We note that if \(\mathbf{X}\) is Gaussian, then \(\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes=\mathbf{0}\).

Cardoso, Mori, Szekely, Rothagi define what we will call the CMRS kurtosis matrix \[ \mathbf{B}(Y) =\operatorname{E}\left( \mathbf{Y}\mathbf{Y}' \mathbf{Y}\mathbf{Y}' \right) -(d+2)\mathbf{I}_d \] which can be expressed in terms of the kurtosis vector as \[ \operatorname{vec}\mathbf{B}(Y)\left( \mathbf{I}_{d^2}\otimes \operatorname{vec}' \mathbf{I}_d \right)\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes \] Note also that \(\operatorname{tr} \mathbf{B}(Y) = \beta_{2,d}\).

For further discussion on the above indexes and further multivariate indexes of skewness and kurtosis, as well as their asymptotic theory one can consult Terdik (2021, section 6) and Jammalamadaka et al. (2021a,b).

The function Esti_Variance_Skew_Kurt provides estimates of the covariance matrix of the data-estimated skewness and kurtosis vectors (Terdik (2021), formulae 6.13 and 6.22).

Example 6. Consider a multivariate data vector of dimension \(d=3\) and \(n=250\) from the multivariate skew-normal distribution of Example 5. The estimated first four cumulants are listed in the object EsMSN obtained by the Esti_EVSK function; the corresponding theoretical values are in the object ThMSN obtained by the istr_SkewNorm_EVSK_Th function.

data<-distr_SkewNorm_Rand(1000,omega,alpha)
EsMSN<-Esti_EVSK(data)
ThMSN<-distr_SkewNorm_EVSK_Th(omega,alpha)

Compare the distinct elements of the estimated skewness vector and the theoretical ones using indx_Elimination.

EsMSN$estSkew[indx_Elimination(3,3)]   
#>  [1]  0.55146965  0.37564103  0.05982354  0.20275976  0.03764258  0.04127893
#>  [7]  0.09830014  0.03208461 -0.06483243  0.12021715
ThMSN$SkewX[indx_Elimination(3,3)]   
#>  [1] 0.68927167 0.34463583 0.00000000 0.17231792 0.00000000 0.00000000
#>  [7] 0.08615896 0.00000000 0.00000000 0.00000000

If one wishes to recover the estimated univariate skewness and kurtosis of the components \(X1\), \(X2\) and \(X3\) of \(X\), then, using indx_UnivMomCum,

EsMSN$estSkew[indx_UnivMomCum(3,3)]  ## Get univariate skewness for X1,X2,X3
#> [1] 0.55146965 0.09830014 0.12021715
EsMSN$estKurt[indx_UnivMomCum(3,4)]  ## Get univariate kurtosis for X1,X2,X3
#> [1]  0.23230281 -0.01456662  0.01875074

An estimate of Mardia’s skewness index is provided together with the p-value under the null hypothesis of normality. The theoretical value of Mardia’s skewness can be recovered from the element SkewX.tot in the object ThMSN.

Esti_Skew_Mardia(data)
#> $Mardia.Skewness
#> [1] 0.9149352
#> 
#> $p.value
#> [1] 1.145771e-27
ThMSN$SkewX.tot
#> [1] 0.9279208

The MRS skewness vector and index are provided together with the p-value for the skewness index under the null hypothesis of normality, The theoretical value, for the distribution at hand, can be computed using formula […]

Esti_Skew_MRSz(data)
#> $MRSz.Skewness.Vector
#> [1] 0.7955083 0.4091087 0.2121253
#> 
#> $MRSz.Skewness.Index
#> [1] 0.8452006
#> 
#> $p.value
#> [1] 3.28982e-18

as.vector(t(c(diag(3))%x%diag(3))%*%ThMSN$SkewX)  ## Theoretical MRS skewness vector
#> [1] 0.8615896 0.4307948 0.0000000

References

Arellano-Valle, R.B., Genton, M.G. (2005) On fundamental skew distributions. Journal of Multivariate Analysis 96(1), 93–116.

Azzalini, A. (2022). The R package ‘sn’: The Skew-Normal and Related Distributions such as the Skew-t and the SUN (version 2.0.2). URL http://azzalini.stat.unipd.it/SN/, https://cran.r-project.org/package=sn

Azzalini, A,, Dalla Valle, A. (1996) The multivariate skew-normal distribution. Biometrika 83(4), 715–726

Chacón, J. E., & Duong, T. (2015). Efficient recursive algorithms for functionals based on higher order derivatives of the multivariate Gaussian density. Statistics and Computing, 25(5), 959-974.

Di Nardo, E. and Guarino, G. (2022). kStatistics: Unbiased Estimators for Cumulant Products and Faa Di Bruno’s Formula. R package version 2.1.1. https://CRAN.R-project.org/package=kStatistics

Franceschini, C. and Loperfido, N. (2017a). MultiSkew: Measures, Tests and Removes Multivariate Skewness. R package version 1.1.1. https://CRAN.R-project.org/package=MultiSkew

Franceschini, C. and Loperfido, N. (2017b). MaxSkew: Orthogonal Data Projections with Maximal Skewness. R package version 1.1. https://CRAN.R-project.org/package=MaxSkew

Holmquist, B. (1996). The d-variate vector Hermite polynomial of order. Linear Algebra and Its Applications, 237/238, 155–190.

Jammalamadaka, S. R., Subba Rao, T. and Terdik, Gy. (2006). Higher order cumulants of random vectors and applications to statistical inference and time series. Sankhya (Ser A, Methodology) 68, 326–356.

Jammalamadaka, S. R., Taufer, E. & Terdik, Gy. H. (2021a). On multivariate skewness and kurtosis. Sankhya A, 83(2), 607-644.

Jammalamadaka, S. R., Taufer, E. & Terdik, Gy. H. (2021b). Asymptotic theory for statistics based on cumulant vectors with applications. Scandinavian Journal of Statistics, 48(2), 708-728.

Jammalamadaka,S.R. , Taufer,E. & Terdik, Gy. (2021c). Cumulants of Multivariate Symmetric and Skew Symmetric Distributions., Symmetry 13, 1383.

Kolda, T.G.; Bader, B.W. (2009). Tensor decompositions and applications. SIAM Rev. 51, 455–500.

Kollo, T. (2008). Multivariate skewness and kurtosis measures with an application in ICA. Journal of Multivariate Analysis 99, 10, 2328–2338.

Komsta, L. and Novomestky F. (2022). moments: Moments, Cumulants, Skewness, Kurtosis and Related Tests. R package version 0.14.1. https://CRAN.R-project.org/package=moments

Novomestky, F. (2021). matrixcalc: Collection of Functions for Matrix Calculations. R package version 1.0-5. https://CRAN.R-project.org/package=matrixcalc

Qi, L. (2006). Rank and eigenvalues of a supersymmetric tensor, the multivariate homogeneous polynomial and the algebraic hypersurface it defines. J. Symb. Comput. 41, 1309–1327.

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

McCullagh, P. (2018). Tensor methods in statistics. Chapman and Hall/CRC.

Móri, T. F., Rohatgi, V. K., & Székely, G. J. (1994). On multivariate skewness and kurtosis. Theory of Probability & Its Applications, 38(3), 547–551.

Ould-Baba, H.; Robin, V.; Antoni, J. (2015). Concise formulae for the cumulant matrices of a random vector. Linear Algebra Appl. 485, 392–416.

Terdik, Gy. (2021). Multivariate statistical methods - going beyond the linear. Springer.