I thought it'd be possible to use any combination of correlations and standard deviations to generate correlated standard normal draws.

The intent is to get around checking whether the variance-covariance matrix $\Sigma$ is positive-definite.

In the three-variable case, my idea is, given a vector of standard deviations $(\sigma_X, \sigma_Y, \sigma_Z)$ and a triple of correlations, $(\rho_{XY}, \rho_{XZ}, \rho_{YZ})$, we construct $\Sigma$ as

$$\Sigma = \begin{bmatrix} \sigma_X & 0 & 0 \\ 0 & \sigma_Y & 0 \\ 0 & 0 & \sigma_Z \end{bmatrix} \begin{bmatrix} 1 & \rho_{XY} & \rho_{XZ} \\ \rho_{XY} & 1 & \rho_{YZ} \\ \rho_{XZ} & \rho_{YZ} & 1 \end{bmatrix} \begin{bmatrix} \sigma_X & 0 & 0 \\ 0 & \sigma_Y & 0 \\ 0 & 0 & \sigma_Z \end{bmatrix} = \begin{bmatrix} \sigma_X^2 & \sigma_{XY} & \sigma_{XZ} \\ \sigma_{XY} & \sigma_Y^2 & \sigma_{YZ} \\ \sigma_{XZ} & \sigma_{YZ} & \sigma_Z^2 \end{bmatrix} $$

With $\Sigma$ in hand, we draw standard normal variates $Z$ and correlate them with the transformation

$$ F = L\Sigma $$

Where $L$ is such that $LL^{T} = \Sigma$ (e.g., Cholesky decomposition). And voila, $F \sim N(0, \Sigma)$!

I thought we'd be able to specify whatever (valid) values of $\sigma_k$ and $\rho_{ij}$ we'd like, but in coding this up I seem to have created an invalid $\Sigma$:

sds = c(.05, .05, .05)
corr = c(0, .385, .97)

#generate center matrix
corr_mat = diag(3)
corr_mat[lower.tri(corr_mat)] = corr_mat[upper.tri(corr_mat)] = corr

#generate sigma
sigma = diag(sds) %*% corr_mat %*% diag(sds)

#get L

Error in chol.default(sigma) : the leading minor of order 3 is not positive definite

How can $\Sigma$ not be positive definite? Is there some restriction on the relative magnitude of standard deviations and correlations? I thought the answer is no...

$\Sigma$ is definitely not positive definite:

# [1]  0.0051090288  0.0025000000 -0.0001090288

I see the negative eigenvalue is very small... is this possibly just a numerical problem?

  • 1
    $\begingroup$ I think you mean $F=LZ$ $\endgroup$ Nov 5 '16 at 1:14
  • $\begingroup$ Possible duplicate? stats.stackexchange.com/questions/124538/… $\endgroup$ Nov 5 '16 at 11:47
  • $\begingroup$ Please consult our many threads on this topic. $\endgroup$
    – whuber
    Nov 5 '16 at 20:29
  • 1
    $\begingroup$ Letting $$v=\left(-\frac{77}{\sqrt{87130}}, -\frac{97\sqrt{2}}{\sqrt{43585}}, \frac{1}{\sqrt{2}}\right) \approx (-0.26086, -0.657231, 0.707107),$$ you may compute that $$v\,\Sigma\,v^\prime = 2 - \frac{\sqrt{8713}}{20\sqrt{5}}\approx -0.087223.$$Since variances cannot be negative, obviously $\Sigma$ cannot be a covariance matrix. The exact arithmetic used in this demonstration shows the issue is not floating point error. $\endgroup$
    – whuber
    Nov 5 '16 at 20:44
  • $\begingroup$ @whuber thanks. The point is I was thinking there's no restrictions on the correlation matrix, which is obviously not true. Things clicked as soon as Jonathan nudged me to question that assumption. $\endgroup$ Nov 5 '16 at 20:49

The issue is that your correlation matrix is not positive definite. You can get this in lot's of ways. One possible way occurs when $cor(X,Y) = .9$, $cor(X,Z) = .9$, $cor(Y,Z) = 0$ this isn't possible for normal variables, and will give you a non-positive definite matrix with eigenvalues (2.2727922, 1.0000000, -0.2727922).

We can show this with a bit more rigor in the case of the random variables $X$, $Y$, and $Z$ distributed $\mathcal{N}(0,1)$. For all three variables, we can write each as the sum of two independent random variables distributed $\mathcal{N}(0,1)$ and multiplied by scalars such that $X=\alpha_1X_1 + \alpha_2X_2$ and $var(X) = \alpha_1^2 + \alpha_2^2 =1$ :

In the case of $X$ and $Y$, where $\mathbb{E}(X_1Y_1)=1$, $\mathbb{E}(X_1Y_2)=0$, and $\mathbb{E}(X_2Y) = 0$ then $\mathbb{E}(XY)=\mathbb{E}((\alpha_1X_1 + \alpha_2X_2)(\beta_1Y_1 + \beta_2Y_2))=\alpha_1\beta_1$.

In the case of $X$ and $Z$, where $\mathbb{E}(X_1Z_1)=1$, $\mathbb{E}(X_1Z_2)=0$, and $\mathbb{E}(X_2Z) = 0$ then $\mathbb{E}(XZ)=\mathbb{E}((\alpha_3X_3 + \alpha_4X_4)(\gamma_1Z_1 + \gamma_2Z_2))=\alpha_3\gamma_1$.

The correlation between $X_1$ and $X_3$ can be zero, if $\alpha_1^2 + \alpha_3^2 < 1$, but if $\mathbb{E}(X,Y)$ and $\mathbb{E}(X,Z)$ are large this may not be possible.

To find values that work, you need to ensure the determinants of the principal minors are positive. To do this use the $2 \times 2$ and $3 \times 3$ formulas for the determinant, plug in the values you know (the diagonal), and check that the resulting value:

The $2\times 2$ Determinant

$\left| \begin{array}{cc} 1 & a \\ a & 1 \end{array} \right| = 1 - a^2 > 0$

The $3\times 3$ Determinant

$\left| \begin{array}{ccc} 1 & a & b \\ a & 1 & c \\ b & c & 1 \end{array} \right| = 1 + 2abc - b^2 - a^2 - c^2 > 0$

  • $\begingroup$ could you elaborate on why it's not possible? I guess if x,y and y,Z are highly correlated, so too must x,z. is there a way to make this relationship exact/bounded? $\endgroup$ Nov 5 '16 at 1:22
  • $\begingroup$ I've updated my answer, let me know if it is unclear. There are probably lots of ways to make the covariance matrix non-positive definite, I'm just thinking of the simplest one I came up with. $\endgroup$ Nov 5 '16 at 3:18
  • $\begingroup$ Your exposition is a bit unclear, though your point is well taken, and the determinant inequality is most useful (is there a geometric interpretation?). How are $X_1$ and $X_2$ constructed -- mainly, how do we know that $\mathbb{E}[X_1 Y_1] = 1$? In my head I was regressing, e.g, $Y$ on $X$ and separately on $Z$. There should be a relationship of the resulting coefficients... $\endgroup$ Nov 5 '16 at 15:14
  • 1
    $\begingroup$ Generally, you need to do more than check the determinant: the determinants of every minor of the matrix must be non-negative, too. For instance, the matrix $$\pmatrix{-2&1\\1&-2}$$ has positive determinant but is not positive-definite. (The $1\times 1$ minors--that is, the diagonal elements--are both negative.) For your $3\times 3$ matrix, the checks of the $2\times 2$ minors amount to confirming that $a,b,c$ all lie in the interval $[-1,1]$. For larger matrices, though, it's important to check minors of all dimensions. $\endgroup$
    – whuber
    Nov 5 '16 at 20:50
  • $\begingroup$ Good point, I'll update the answer when I get a chance. $\endgroup$ Nov 5 '16 at 21:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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