2

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
chol(sigma)

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:

eigen(sigma)$values
# [1]  0.0051090288  0.0025000000 -0.0001090288

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

MichaelChirico
  • 1,270
  • 1
  • 9
  • 20
  • 1
    I think you mean $F=LZ$ – Jonathan Lisic Nov 05 '16 at 01:14
  • Possible duplicate? http://stats.stackexchange.com/questions/124538/how-to-generate-a-large-full-rank-random-correlation-matrix-with-some-strong-cor – Matthew Gunn Nov 05 '16 at 11:47
  • Please consult our [many threads on this topic](http://stats.stackexchange.com/search?q=correlation+positive+definite). – whuber Nov 05 '16 at 20:29
  • 1
    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. – whuber Nov 05 '16 at 20:44
  • @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. – MichaelChirico Nov 05 '16 at 20:49

1 Answers1

4

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$

Jonathan Lisic
  • 1,342
  • 7
  • 16
  • 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? – MichaelChirico Nov 05 '16 at 01:22
  • 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. – Jonathan Lisic Nov 05 '16 at 03:18
  • 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... – MichaelChirico Nov 05 '16 at 15:14
  • 1
    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. – whuber Nov 05 '16 at 20:50
  • Good point, I'll update the answer when I get a chance. – Jonathan Lisic Nov 05 '16 at 21:08