How to use the Cholesky decomposition, or an alternative, for correlated data simulation

Question

I use Cholesky decomposition to simulate correlated random variables given a correlation matrix. The thing is, the result never reproduces the correlation structure as it is given. Here is a small example in Python to illustrate the situation.

import numpy as np    

n_obs = 10000
means = [1, 2, 3]
sds = [1, 2, 3] # standard deviations 

# generating random independent variables 
observations = np.vstack([np.random.normal(loc=mean, scale=sd, size=n_obs)
                   for mean, sd in zip(means, sds)])  # observations, a row per variable

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])

L = np.linalg.cholesky(cor_matrix)

print(np.corrcoef(L.dot(observations))) 

This prints:

[[ 1.          0.34450587  0.57515737]
 [ 0.34450587  1.          0.1488504 ]
 [ 0.57515737  0.1488504   1.        ]]

As you can see, the post-hoc estimated correlation matrix drastically differs from the prior one. Is there a bug in my code, or is there some alternative to using the Cholesky decomposition?

Edit

I beg your pardon for this mess. I didn't think there was an error in the code and/or in the way Cholesky decomposition was applied due to some misunderstanding of the material I had studied before. In fact I was sure that the method itself was not meant to be precise and I had been okay with that up until the situation that made me post this question. Thank you for pointing at the misconception I had. I've edited the title to better reflect the real situation as proposed by @Silverfish.

Answer 1

The approach based on the Cholesky decomposition should work, it is described here and is shown in the answer by Mark L. Stone posted almost at the same time that this answer.

Nevertheless, I have sometimes generated draws from the multivariate Normal distribution $N(\vec\mu, \Sigma)$ as follows:

$$ Y = Q X + \vec\mu \,, \quad \hbox{with}\quad Q=\Lambda^{1/2}\Phi \,, $$

where $Y$ are the final draws, $X$ are draws from the univariate standard Normal distribution, $\Phi$ is a matrix containing the normalized eigenvectors of the target matrix $\Sigma$ and $\Lambda$ is a diagonal matrix containing the eigenvalues of $\Sigma$ arranged in the same order as the eigenvectors in the columns of $\Phi$.

Example in R (sorry I'm not using the same software you used in the question):

n <- 10000
corM <- rbind(c(1.0, 0.6, 0.9), c(0.6, 1.0, 0.5), c(0.9, 0.5, 1.0))
set.seed(123)
SigmaEV <- eigen(corM)
eps <- rnorm(n * ncol(SigmaEV$vectors))
Meps <- matrix(eps, ncol = n, byrow = TRUE)    
Meps <- SigmaEV$vectors %*% diag(sqrt(SigmaEV$values)) %*% Meps
Meps <- t(Meps)
# target correlation matrix
corM
#      [,1] [,2] [,3]
# [1,]  1.0  0.6  0.9
# [2,]  0.6  1.0  0.5
# [3,]  0.9  0.5  1.0
# correlation matrix for simulated data
cor(Meps)
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.6002078 0.8994329
# [2,] 0.6002078 1.0000000 0.5006346
# [3,] 0.8994329 0.5006346 1.0000000

---

You may be also interested in this post and this post.

Answer 2

People would be likely find your error much faster if you explained what you did with words and algebra rather than code (or at least writing it using pseudocode).

You appear to be doing the equivalent of this (though possibly transposed):

1. Generate an $n\times k$ matrix of standard normals, $Z$

2. multiply the columns by $\sigma_i$ and add $\mu_i$ to get nonstandard normals

3. calculate $Y=LX$ to get correlated normals.

where $L$ is the left Cholesky factor of your correlation matrix.

What you should do is this:

1. Generate an $n\times k$ matrix of standard normals, $Z$

2. calculate $X=LZ$ to get correlated normals.

3. multiply the columns by $\sigma_i$ and add $\mu_i$ to get nonstandard normals

There's many explanations of this algorithm on site. e.g.

How to generate correlated random numbers (given means, variances and degree of correlation)?

Can I use the Cholesky-method for generating correlated random variables with given mean?

This one discusses it directly in terms of the desired covariance matrix, and also gives an algorithm for getting a desired sample covariance:

Generating data with a given sample covariance matrix

Answer 3

There's nothing wrong with the Cholesky factorization. There is an error in your code. See edit below.

Here is MATLAB code and results, first for n_obs = 10000 as you have, then for n_obs = 1e8. For simplicity, since it doesn't affect the results, I don't bother with means, i.e., I make them zeros. Note that MATLAB's chol produces an upper triangular Cholesky factor R of the matrix M such that R' * R = M. numpy.linalg.cholesky produces a lower triangular Cholesky factor, so an adjustment vs. my code is needed; but I believe your code is fine in that respect.

   >> correlation_matrix = [1.0, 0.6, 0.9; 0.6, 1.0, 0.5;0.9, 0.5, 1.0];
   >> SD = diag([1 2 3]);
   >> covariance_matrix = SD*correlation_matrix*SD
   covariance_matrix =
      1.000000000000000   1.200000000000000   2.700000000000000
      1.200000000000000   4.000000000000000   3.000000000000000
      2.700000000000000   3.000000000000000   9.000000000000000
   >> n_obs = 10000;
   >> Random_sample = randn(n_obs,3)*chol(covariance_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.599105015695768   0.898395949647890
      0.599105015695768   1.000000000000000   0.495147514173305
      0.898395949647890   0.495147514173305   1.000000000000000
   >> n_obs = 1e8;
   >> Random_sample = randn(n_obs,3)*chol(covariance_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.600101477583914   0.899986072541418
      0.600101477583914   1.000000000000000   0.500112824962378
      0.899986072541418   0.500112824962378   1.000000000000000

Edit: I found your mistake. You incorrectly applied the standard deviation. This is the equivalent of what you did, which is wrong.

   >> n_obs = 10000;
   >> Random_sample = randn(n_obs,3)*SD*chol(correlation_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.336292731308138   0.562331469857830
      0.336292731308138   1.000000000000000   0.131270077244625
      0.562331469857830   0.131270077244625   1.000000000000000
   >> n_obs=1e8;
   >> Random_sample = randn(n_obs,3)*SD*chol(correlation_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.351254525742470   0.568291702131030
      0.351254525742470   1.000000000000000   0.140443281045496
      0.568291702131030   0.140443281045496   1.000000000000000

Answer 4

CV is not about code, but I was intrigued to see how this would look after all the good answers, and specifically @Mark L. Stone contribution. The actual answer to the question is provided on his post (please credit his post in case of doubt). I'm moving this appended info here to facilitate retrieval of this post in the future. Without playing down any of the other excellent answers, after Mark's answer, this wraps up the issue by correcting the post in the OP.

Source

IN PYTHON:

import numpy as np

no_obs = 1000             # Number of observations per column
means = [1, 2, 3]         # Mean values of each column
no_cols = 3               # Number of columns

sds = [1, 2, 3]           # SD of each column
sd = np.diag(sds)         # SD in a diagonal matrix for later operations

observations = np.random.normal(0, 1, (no_cols, no_obs)) # Rd draws N(0,1) in [3 x 1,000]

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])          # The correlation matrix [3 x 3]

cov_matrix = np.dot(sd, np.dot(cor_matrix, sd))   # The covariance matrix

Chol = np.linalg.cholesky(cov_matrix)             # Cholesky decomposition

array([[ 1.        ,  0.        ,  0.        ],
       [ 1.2       ,  1.6       ,  0.        ],
       [ 2.7       , -0.15      ,  1.29903811]])

sam_eq_mean = Chol .dot(observations)             # Generating random MVN (0, cov_matrix)

s = sam_eq_mean.transpose() + means               # Adding the means column wise
samples = s.transpose()                           # Transposing back

print(np.corrcoef(samples))                       # Checking correlation consistency.

[[ 1.          0.59167434  0.90182308]
 [ 0.59167434  1.          0.49279316]
 [ 0.90182308  0.49279316  1.        ]]

---

IN [R]:

no_obs = 1000             # Number of observations per column
means = 1:3               # Mean values of each column
no_cols = 3               # Number of columns

sds = 1:3                 # SD of each column
sd = diag(sds)         # SD in a diagonal matrix for later operations

observations = matrix(rnorm(no_cols * no_obs), nrow = no_cols) # Rd draws N(0,1)

cor_matrix = matrix(c(1.0, 0.6, 0.9,
                      0.6, 1.0, 0.5,
                      0.9, 0.5, 1.0), byrow = T, nrow = 3)     # cor matrix [3 x 3]

cov_matrix = sd %*% cor_matrix %*% sd                          # The covariance matrix

Chol = chol(cov_matrix)                                        # Cholesky decomposition

     [,1] [,2]      [,3]
[1,]    1  1.2  2.700000
[2,]    0  1.6 -0.150000
[3,]    0  0.0  1.299038

sam_eq_mean = t(observations) %*% Chol          # Generating random MVN (0, cov_matrix)

samples = t(sam_eq_mean) + means

cor(t(samples))

          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6071067 0.8857339
[2,] 0.6071067 1.0000000 0.4655579
[3,] 0.8857339 0.4655579 1.0000000

colMeans(t(samples))
[1] 1.035056 2.099352 3.065797
apply(t(samples), 2, sd)
[1] 0.9543873 1.9788250 2.8903964

Answer 5

As others have already shown: cholesky works. Here a piece of code which is very short and very near to pseudocode: a codepiece in MatMate:

Co = {{1.0, 0.6, 0.9},  _
      {0.6, 1.0, 0.5},  _
      {0.9, 0.5, 1.0}}           // make correlation matrix


chol = cholesky(co)              // do cholesky-decomposition           
data = chol * unkorrzl(randomn(3,100,0,1))  
                                 // dot-multiply cholesky with random-
                                 // vectors with mean=0, sdev=1  
                                 //(refined by a "decorrelation" 
                                 //to remove spurious/random correlations)   


chk = data *' /100               // check the correlation of the data
list chk

1.0000  0.6000  0.9000
0.6000  1.0000  0.5000
0.9000  0.5000  1.0000

Answer 6

Python code:

import numpy as np

# desired correlation matrix
cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])

L = np.linalg.cholesky(cor_matrix)

# build some signals that will result in the desired correlation matrix
X = L.dot(np.random.normal(0,1,(3,1000))) # the more the sample (1000) the better

# estimate their correlation matrix
np.corrcoef(X)
array([[1.        , 0.58773667, 0.8978625 ],
       [0.58773667, 1.        , 0.47424997],
       [0.8978625 , 0.47424997, 1.        ]])

# Very good approxiamation :)