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)))
[[ 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?
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.