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.

  • 1
    Cholesky works just fine, and this is really a "can you find the bug in my code" type question. The title and content of the question, as it is originally written, are basically "Cholesky doesn't work, what's an alternative"? That will be very confusing to users searching this site. Should this question be edited to reflect this? (The downside is that javlacalle's answer would be less relevant. The upside is the question text would then reflect what searchers would actually find on the page.) – Silverfish Jul 6 '15 at 10:06
  • @Antoni Parellada Yes, i think you have translated my MATLAB code for the (a) correct way of doing it into Python numpy, complete with adjustment for np.linalg.cholesky being lower triangular vs. MATLAB's chol being upper triangular. I had already translated the OP's incorrect code to its MATLAB equivalent and duplicated his incorrect results. – Mark L. Stone Sep 1 '16 at 2:13
up vote 9 down vote accepted

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.

  • To make the reproduced correlation-matrix precise one should remove the spurious correlations in the random-data from the random-generator before applying it to the data-generation-procedure. For instance, check the correlation of your random-data in eps to see that spurious correlations first. – Gottfried Helms Jul 19 '15 at 11:44

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

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

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.


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

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

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.