Simulating covariance matrices with nonzero correlation

Question

How would you simulate a covariance matrix of 1,000 stocks where each pair has nonzero correlation?

I have literally no idea how to start with this.

Any suggestions?

Answer 1

What does 'simulate a covariance matrix' mean?

• If the question means, generate an arbitrary correlation matrix for 1000 stocks, then we can choose any symmetric matrix with all 1s down the diagonal, so long as every element is between -1 and 1 and the matrix is positive semi-definite. The large size of the matrix means that putting random values in every cell will almost certainly fail the positive semi-definiteness test, so I would start with a 1000*1000 identity matrix and add a small random positive or negative amount to a random cell (and its reflection) and check if the new matrix passes the test, then repeat this process to gradually build up a valid matrix. Then we map it from a correlation matrix to a covariance matrix by multiplying each entry by the product of the square roots of the variances of the two corresponding price series

• If the question is then asking us to simulate stock prices that obey the given correlation matrix, we need to generate uncorrelated price series, then do a Cholesky Decomposition on the correlation matrix and apply that to transform the uncorrelated prices (actually probably on the covariance matrix but it should be straight-forward to translate between the two as descibed above...)

• If the question is simply asking us to calculate the covariance matrix for 1000 stocks given their price histories, it's just a case of calculating every single pairwise covariance and putting them into a matrix (up to issues like whether the stock prices are sampled at the same times...)

Answer 2

I will just clarify Point 2 in StackG excellent answer. (It's really a comment, but it's too long and has too much math symbols to fit in the comment field.)

Suppose you're given a covariance matrix $C$ for the returns of $n$ assets. (1000 $\times$ 1000 is 1 million entries - should not be too large for modern computers to work with, but do be mindful of your memory requirements.)

You want to simulate the returns of the assets that would be consistent with the volatilities and correlations in $C$. You assume normal distribution for the return of each asset $N(\mu,\sigma^2)$.

$f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{\frac{(x-\mu)^2}{2\sigma}}$, where $\mu$ is the mean, which you will assume to be 0, and $\sigma$ is the standard deviation (on $C$'s diagonal).

You compute $n\times n$ matrix $H$ such that $H \times H^T = C$. A convenient numerical method for doing that is Choleski decomposition. However Choleski requires $C$ to be positive definite. In practice, you might have some asset that's a linear combination of other assets; or you might have numeric noise; or some historical time series shorter than others. If $C$ is not, but is not very far from being positive definite, then there are methods of tweaking it a little to get a positive definite matrix, so that you can use Choleski.

You generate $Z$ of $n$ normally distributed random numbers. You can do it by first generating uniformly distributed random numbers using your favorite pseudo-random number generator and then using Box-Muller transformation. Observe that $Z$ has a normal distribution with mean 0 and no correlation.

Then you just multiply $Y=HZ$. Each scenario $Y$ is normally distributed with mean 0 and covariance $C$, because:

$\overline{Y} = H\overline{Z} = 0.$

$C_Y = \overline{(Y-\overline{Y})(Y-\overline{Y})^T} = \overline{(HZ)(HZ)^T} = \overline{HZZ^TH^T} = H \overline{ZZ^T} H^T = HIH^T = HH^T = C.$

Answer 3

I would use Numpy (a library of Python) to do it. There's a function called numpy.random.multivariate_normal. It takes in 2 main arguments, an array of means (expected returns of the stocks) and an array (matrix) of covariances of the stocks.