Multivariate sample kurtosis due to Koziol (1989) is defined by
$$\widetilde{b}_{n,d}^{(2)}=\frac{1}{n^2}\sum_{j,k=1}^n(Y_{n,j}^\top Y_{n,k})^4,$$
where \(Y_{n,j}=S_n^{-1/2}(X_j-\overline{X}_n)\), \(j=1,\ldots,n\), are the scaled residuals, \(\overline{X}_n\) is the sample mean and \(S_n\) is the sample covariance matrix of the random vectors \(X_1,\ldots,X_n\). To ensure that the computation works properly
\(n \ge d+1\) is needed. If that is not the case the function returns an error. Note that for \(d=1\), we have a measure proportional to the squared sample kurtosis.