# 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. – Glen_b Jun 19 '13 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. – Glen_b Jun 19 '13 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 – S. Punky Jun 20 '13 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! – S. Punky Jun 20 '13 at 1:23
• "not limiting myself to finding a method that only satisfies Mardia's definitions" -- method of doing what? – Glen_b Jun 20 '13 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? – Sextus Empiricus Jun 29 '18 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 – S. Punky Jun 22 '13 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. – S. Punky Jun 23 '13 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? – S. Punky Jun 26 '13 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 – Quartz Jun 27 '13 at 10:10

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 Feb 27 '14 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. – S. Punky Feb 28 '14 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. – Francis Smart Feb 28 '14 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 :) – S. Punky Mar 1 '14 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.