Correlate normal shocks

Question

I am trying to generate some random standard normal variables and correlate them In particular I want:

$$ \bf Y \sim \mathcal N(0, \Sigma) $$

where $\textbf{Y} = (Y_1,\dots,Y_n)$ is the vector I want to simulate, and $\Sigma$ the given covariance matrix.

I am doing the following in matlab:

1. Simulate a vector of uncorrelated Gaussian random variables, $\bf Z $
2. then find a square root of $\Sigma$, i.e. a matrix $\bf C$ such that $\bf C \bf C^\intercal = \Sigma$.

Then the target vector is given by $$ \bf Y = \bf C \bf Z. $$

Here is my matlab code:

N = 500000
u_1 = normrnd(zeros(N,1),1);
u_2 = normrnd(zeros(N,1),1);
u_3 = normrnd(zeros(N,1),1);
u_4 = normrnd(zeros(N,1),1);

rv = [u_1 '; u_2'; u_3'; u_4'];


VarCov = [1 -0.87 0.0 -0.6;
          -0.87 1 0.0 0.7;
          0.0 0.0 1 0.0
          -0.6 0.0 0.0 1];

ch = chol(VarCov);
result = ch * allshocks;

However, when I then compute the means and covariances of the resulting vector, the means are not identical to zero (not even close) and the covariances are are also off the ballpark. Am I doing something wrong?

Answer 1

The theory is correct; the code isn't. $\Sigma$ has to be a symmetric positive semi-definite while in your code, VarCov isn't even symmetric. You should get the expected results once this is fixed - assuming there aren't any other bugs in your code (its been a while since I used Matlab).