# How can I generate data with a prespecified correlation matrix?

I’m trying to generate correlated random sequence with mean = $0$, variance = $1$, correlation coefficient = $0.8$. In the code below, I use s1 & s2 as the standard deviations, and m1 & m2 as the means.

p = 0.8
u = randn(1, n)
v = randn(1, n)
x = s1 * u + m1
y = s2 * (p * u + sqrt(1 - p^2) * v) + m2


This gives me the correct corrcoef() of 0.8 between x and y. My question is how can I generate a series means if I want z that is also correlated with y (with the same correlation $r=0.8$), but not with x. Is there a particular formula I need to know? I found one but couldn't understand it. See this link.

-
–  gung Sep 10 at 20:10

It appears that you're asking how to generate data with a particular correlation matrix.

A useful fact is that if you have a random vector ${\bf x}$ with covariance matrix $\Sigma$, then the random vector ${\bf Ax}$ has mean ${\bf A} E({\bf x})$ and covariance matrix $\Omega = {\bf A} \Sigma {\bf A}^{T}$. So, if you start with data that has mean zero, multiplying by ${\bf A}$ will not change that, so your first requirement is easily satisfied.

Let's say you start with (mean zero) uncorrelated data (i.e. the covariance matrix is diagonal) - since we're talking about the correlation matrix, let's just take $\Sigma = I$. You can transform this to data with a given covariance matrix by choosing ${\bf A}$ to be the cholesky square root of $\Omega$ - then ${\bf Ax}$ would have the desired covariance matrix $\Omega$.

In your example, you appear to want something like this:

$$\Omega = \left( \begin{array}{ccc} 1 & .8 & 0 \\ .8 & 1 & .8 \\ 0 & .8 & 1 \\ \end{array} \right)$$

Unfortunately that matrix is not positive definite, so it cannot be a covariance matrix - you can check this by seeing that the determinant is negative. Perhaps, instead

$$\Omega = \left( \begin{array}{ccc} 1 & .8 & .3 \\ .8 & 1 & .8 \\ .3 & .8 & 1 \\ \end{array} \right) \ \ \ \ {\rm or} \ \ \ \Omega = \left( \begin{array}{ccc} 1 & 2/3 & 0 \\ 2/3 & 1 & 2/3 \\ 0 & 2/3 & 1 \\ \end{array} \right)$$

would suffice. I'm not sure how to calculate the cholesky square root in matlab (which appears to be what you're using) but in R you can use the chol() function.

In this example, for the two $\Omega$s listed above the proper matrix multiples(respectively) would be

$${\bf A} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ .8 & .6 & 0 \\ .3 & .933 & .1972 \\ \end{array} \right) \ \ \ \ {\rm or} \ \ \ {\bf A} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 2/3 & .7453 & 0 \\ 0 & .8944 & .4472 \\ \end{array} \right)$$

The R code used to arrive at this was:

x = matrix(0,3,3)
x[1,]=c(1,.8,.3)
x[2,]=c(.8,1,.8)
x[3,]=c(.3,.8,1)
t(chol(x))

[,1]      [,2]      [,3]
[1,]  1.0 0.0000000 0.0000000
[2,]  0.8 0.6000000 0.0000000
[3,]  0.3 0.9333333 0.1972027

x[1,]=c(1,2/3,0)
x[2,]=c(2/3,1,2/3)
x[3,]=c(0,2/3,1)
t(chol(x))

[,1]      [,2]      [,3]
[1,] 1.0000000 0.0000000 0.0000000
[2,] 0.6666667 0.7453560 0.0000000
[3,] 0.0000000 0.8944272 0.4472136

-
The MATLAB function is also called chol. Note that this can be quite numerically unstable if $\Omega$ is almost singular. In that case, using the symmetric square-root obtained, e.g., via the SVD, is often a better choice in terms of numerical stability. :) –  cardinal Jul 13 '12 at 13:09
Of course that's right @cardinal - lots of theoretically justified things go bad when you're trying to do things numerically with nearly singular matrices. I was (conveniently) imagining the situation where the target correlation matrix was not in the realm where this was an issue. It's good that you pointed this out - thanks (and thanks for the edit to my other answer) –  Macro Jul 13 '12 at 13:34
The main reason I was thinking about this was due to your keen eye in recognizing that the OP's first suggestion was not even positive definite. And, hopefully the edit to the other question was not overzealous; I like both these answers. –  cardinal Jul 13 '12 at 13:46

If you're using R, you can also use the mvrnorm function from the MASS package, assuming you want normally distributed variables. The implementation is similar to Macro's description above, but uses the eigenvectors of the correlation matrix instead of the cholesky decomposition and scaling with a singular value decomposition (if the empirical option is set to true).

If $X$ is a matrix with entries drawn from a normal distribution, $\Sigma$ is a positive definite correlation matrix with eigenvectors $\gamma$, and $\lambda$ is a square matrix with the square root eigen values from $\Sigma$ along the diagonal then:

$X' = \gamma\lambda X^{T}$

Where X' is a normally distributed matrix with correlation matrix of $\Sigma$ and column means are the same as $X$.

Note that the correlation matrix have to be positive definite, but converting it with the nearPD function from the Matrix package in R will be useful.

-

An alternative solution without cholesky factorization is the following. Let $\Sigma_y$ the desired covariance matrix and suppose you have data $x$ with $\Sigma_x = I$. Suppose $\Sigma_y$ is positive definite with $\Lambda$ the diagonal matrix of the eigenvalues and $V$ the matrix of column eigenvectors .

You can write $\Sigma_y = V \Lambda V^T = ( V \sqrt{\Lambda} ) (\sqrt{\Lambda}^T V^T ) = A A^T$.

$y=Ax$ generate the desired data.

-