# Method for generating correlated non-normal data

I'm interested in finding out a method for generating correlated, non-normal data. So ideally some sort of distribution that takes in a covariance (or correlation) matrix as a parameter and generates data that approximates it. But here's the catch: the method I'm trying to find should have the flexibility to also control its multivariate skewness and / or kurtosis.

I'm familiar Fleishman's method and the use of the power method of normal variates, but I believe most of those extensions only allow the user for certain combinations of marginal skewness and kurtosis, leaving the multivariate skewness / kurtosis just out there. What I was wondering is if there is a method that helps specify the multivariate skewness and / or kurtosis, alongside with some correlation / covariance structure.

About a year ago I took a seminar on copula distributions and I remember the professor casually mentioning that through the use of vine copulas, one could generate data which is, say, symmetric in each one of its 1-D marginals but jointly skewed and vice-versa. Or, even further, that any lower-dimensional margins could have some skewness or kurtosis while keeping the highest dimensions symmetric (or not). I was marveled by the idea that such flexibility could exist I've been trying to find some sort of article or conference paper that describes said method but I have been unsuccessful :(. It doesn't have to be through the use of copulas, I'm open to anything that works.

Edit: I have added some R code to try to show what I mean. So far I am only well-acquainted with Mardia's definition of multivariate skewness and kurtosis. When I first approached my problem I naively thought that if I used a symmetric copula (Gaussian in this case) with skewed marginals (beta, in this example), univariate tests on the marginals would yield significance but Mardia's test for multivarite skewness/kurtosis would be non-significant. I tried that and it didn't come out as I had expected:

library(copula)
library(psych)
set.seed(101)

cop1 <- {mvdc(normalCopula(c(0.5), dim=2, dispstr="un"),
c("beta", "beta"),list(list(shape1=0.5, shape2=5),
list(shape1=0.5, shape2=5)))}

Q1 <- rmvdc(cop1, 1000)
x1 <- Q1[,1]
y1 <- Q1[,2]

cop2 <- {mvdc(normalCopula(c(0.5), dim=2, dispstr="un"),
c("norm", "norm"),list(list(mean=0, sd=1),
list(mean = 0, sd=1)))}

Q2 <- rmvdc(cop2, 1000)
x2 <- Q2[,1]
y2 <- Q2[,2]

mardia(Q1)

Call: mardia(x = Q1)

Mardia tests of multivariate skew and kurtosis
Use describe(x) the to get univariate tests
n.obs = 1000   num.vars =  2
b1p =  10.33   skew =  1720.98  with probability =  0
small sample skew =  1729.6  with probability =  0
b2p =  22.59   kurtosis =  57.68  with probability =  0

mardia(Q2)
Call: mardia(x = Q2)

Mardia tests of multivariate skew and kurtosis
Use describe(x) the to get univariate tests
n.obs = 1000   num.vars =  2
b1p =  0.01   skew =  0.92  with probability =  0.92
small sample skew =  0.92  with probability =  0.92
b2p =  7.8   kurtosis =  -0.79  with probability =  0.43


Upon inspecting the contours for 'cop1' VS 'cop2' as well as the empirical bivariate density plots, I can also see that none of them look symmetric at all. That's when I realized this is probably a little more complicated than I thought.

I know that Mardia's is not the only definition of multivariate skewness/kurtosis, so I'm not limiting myself to finding a method that only satisfies Mardia's definitions.

thank you!

• +1 A most interesting question. Could you be more specific about what 'jointly skewed' means in this context (particularly a bivariate one)? While I can picture forms of joint distribution that are in some sense "different" in the four quadrants (about axes placed at the means, say), I'm not familiar with what "jointly skewed" might specifically refer to. Commented Jun 19, 2013 at 22:45
• As for using copulas; plainly a copula can be symmetric (in various senses) while the marginals are skewed, since the copula is transformed to marginal uniformity. So even something as simple as a multivariate lognormal ($\exp X$ where $X$ is multivariate normal) has skewed margins and a 'symmetric' copula (in the senses that the copula of a multivariate normal is symmetric, at least). But that's not (I assume) what you're asking about. Commented Jun 19, 2013 at 22:49
• hello. thank you very much for taking an interest in my question. this is the first time i post here so i hope i am doing things correctly. i will elaborate more on the comment section because the character limit prevents me from using R code to try and convey what i'm doing Commented Jun 20, 2013 at 1:14
• yes, i just realized that and added more detail. i appreciate you taking the time to guide me as far as how to use this board. thanks! Commented Jun 20, 2013 at 1:23
• "not limiting myself to finding a method that only satisfies Mardia's definitions" -- method of doing what? Commented Jun 20, 2013 at 1:33

After much searching, jumping around online forums, consulting with professors and doing A LOT of literature review, I have come to the conclusion that probably THE only way to address this problem is through the use of vine copulas indeed. It gives you some control over the pairwise skewness and kurtosis (or any higher moments) - for a p-variate random vector and the freedom to specify p-1 pair of copulas and the remaining p*(p-1)/2 - (p-1) dimensions can be specified in some kind of conditional copula.

I welcome other methods people might've come across but at least I'm going to leave this pointer towards an answer because i cannot, for the life of me, find any other ways to address this.

• What is a vine copula? Commented Jun 29, 2018 at 14:17

You might be able to solve this by modifying Ruscio and Kaczetow's (2008) algorithm. Their paper provides an iterative algorithm (with R code) that minimizes the difference between the actual and intended marginal shapes. You might be able to modify it so that its targeting the multivariate (rather than marginal) moments.

Ruscio, J., & Kaczetow, W. (2008). Simulating multivariate nonnormal data using an iterative algorithm. Multivariate Behavioral Research, 43(3), 355‐381. doi:10.1080/00273170802285693

• OMG! THANK YOU! i thought for a moment that this question would just be swallowed into oblivion Commented Jun 22, 2013 at 0:07
• well... i have reviewed the Ruscio & Kaczetow (2008) article. sadly it's just another (yet more flexible) implementation of the NORTA (NORmal To Anything) family of algorithms which is known to not work well with multivariate 3rd & 4th moments. i guess i'm back to square one on this one. Commented Jun 23, 2013 at 5:30

You might want to check the Generalized Elliptical Distribution, which allows for a "classical" shape matrix with flexibility for other features.

• Thank you! I will make sure to check this link out. Now, aren't elliptical distributions symmetric? So one can control the kurtosis but the skewness must remain at 0? Commented Jun 26, 2013 at 21:59
• Sure, but GE does not imply elliptical. For some skew elliptical variations check also here: stat.tamu.edu/~genton/STAT689/TAMU2009SE.pdf Commented Jun 27, 2013 at 10:10

I know this question was asked several years ago so the response is too late, but in case somebody else has a similar issue the following reference may be useful:

Qu et al. (2020). A method of generating multivariate non-normal random numbers with desired multivariate skewness and kurtosis. Behavior Research Methods volume 52, 939-946

– whuber
Commented Jan 23, 2021 at 22:05

I have come up with a simple method for doing this that does not involve coplas and other complex designs. I am afraid I do not have any formal reference though the method appears to be highly effective.

The idea is simple. 1. Draw any number of variables from a joint normal distribution. 2. Apply the univariate normal CDF of variables to derive probabilities for each variable. 3. Finally apply the inverse CDF of any distribution to simulate draws from that distribution.

I came up with this method in 2012 and demonstrated using Stata. I have also written a recent post showing the same method using R.

• (1) What is a "Spearman normal distribution"? (2) What distinction are you making, if any, between a CDF and a "normal CDF"? (3) Could you explain how this method introduces any correlation at all? I'm afraid your general uses of "variable" and "distribution" make your description rather vague, so it is hard to tell what it's really doing. Could you reword your answer to be more precise?
– whuber
Commented Feb 27, 2014 at 19:08
• thank you for your post! by following the links one can see more information on the method. it doesn't quite do what i was hoping to attain (i.e. control over the higher-order, higher-dimensional moments of the distribution) but still a very valuable approach. Commented Feb 28, 2014 at 9:04
• Unsurprisingly, I did not come up with a new method see: Cario, Marne C., and Barry L. Nelson. Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix. Technical Report, Department of Industrial Engineering and Management Sciences, Northwestern University, Evanston, Illinois, 1997. Yahav, Inbal, and Galit Shmueli. "On generating multivariate poisson data in management science applications." Robert H. Smith School Research Paper No. RHS (2009): 06-085. Commented Feb 28, 2014 at 9:12
• even if it's not a 'new method', i would still like to thank you for taking the time to look over my question and add something insightful :) Commented Mar 1, 2014 at 22:05

I believe the method presented in the following papers permits generating random multivariates with any (feasible) combination of mean, variance, skewness, and kurtosis.

1. Stanfield, P.M., Wilson, J.R., and Mirka, G.A. 1996. Multivariate Input Modeling with Johnson Distributions, Proceedings of the 1996 Winter Simulation Conference, eds. Charnes, J.M, Morrice, D.J., Brunner, D.T., and Swain, J.J., 1457-1464.
2. Stanfield, P.M., Wilson, J.R., and King, R.E. 2004. Flexible modelling of correlated operation times with application in product re-use facilities, International Journal of Production Research, Vol 42, No 11, 2179–2196.

Disclaimer: I am not one of the authors.