# How to simulate data with values of Mardia's Kurtosis?

I have two (related) questions regarding the simulation of correlated, non-normal data with a user-specified value of Mardia’s 1970 multivariate kurtosis.

(1) It is very common in my field to use the method developed by Vale and Maurelli (1983) of doing a polynomial transformation of standard, normal random variables ($Y=a + bX + cX^{2} + dX^{3}$ where $X \sim N(0,1)$) where the coefficients c and d control the (univariate) skewness and kurtosis estimates. The non-normality is then achieved through the 1-dimensional marginal distributions but there is no control over the multivariate skewness/kurtosis estimates (as defined in Mardia (1970). So here comes my first question:

How can I obtain the population values of multivariate kurtosis if I have only the population values of all the univariate marginals?

Here's my attempt at a solution. It's using the semTools package in R (calculates Mardia's measure) and the lavaan package (implements the Vale & Maurelli method)

    library(semTools)
library(lavaan)

### three variables all correlated at 0.5
model <- 'x1 ~~ 0.5*x2
x2 ~~ 0.5*x3
x3 ~~ 0.5*x1'

### vector to store Mardia's kurtosis values
mardia_values <- double(100)

for (i in 1:100){

### generate simulated data with N=100,000 population skewness of 2 and kurtosis of 7
dataz <- simulateData(model, sample.nobs=100000, skewness=c(2,2,2), kurtosis=c(7,7,7))

mardia_values[i] <- as.numeric(mardiaKurtosis(datz))
}

mean(mardia_values)


So essentially I'm approximating it via a simulation (at a large sample size)

Would this method be acceptable? I'm very open to suggestions (particularly if it can be derived analytically

(2) Now, on a related note... does anyone know of any method to simulate data where the value of Mardia Kurtosis can be specified by the user? So far, the only thing I have been able to come up with is, once again, using the Vale and Maurelli (1983) method and sort of just go by trail and error until I get the value that I want. Since Mardia defined these measures in the population, I assume it should be possible to have some method to generate data where the value can be specified in advance, but I do not know how.

Thank you!

References

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

Vale, C. D., & Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.

• Have you considered using Yuan, Chan, and Bentler's (2000) multivariate transformation approach? Their goal was to transform toward multivariate normality, but you might be able to reverse their approach, that is, simulate multivariate normals, and then transform to your desired values of multivariate kurtosis.www3.nd.edu/~kyuan/papers/Yuan-Chan-Bentler-BJMSP00.pdf – Anthony Jan 8 '15 at 15:14