How can I generate data with a prespecified correlation matrix?

Question

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.

Answer 1

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

Answer 2

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.

Answer 3

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.