how to generate specific random covariance matrices?

Question

I am trying to create random covariance matrices for three joint gaussian variables. My goal is to sample random covariances matrices that always have correlation between 0.7 and 0.9 (or 0 if there isn’t).

So far I am doing it manually with a repeat until is.positive.definite is true… But I am unable to achieve it, my repeat takes a lot of time because most of my matrices samples return false for the positive.definite.

Is there a library to do this or an simpler approach for this?

On the math side I know I can have correlation between: $X_1$ and $X_2$. $X_2$ and $X_3$. $X_1$ and $X_3$ If I am not mistaken, I can have correlation between the three pair or just one pair, there shouldn’t be any issue. But if there is correlation between two of them, the remaining correlation couldn’t be 0, otherwise the matrix would never be positive definite…

Answer 1

In this particular case there's a simple, easy, completely general method.

Break the problem down into two parts:

1. Generate random variances $\sigma_i^2,$ $i=1,2,3.$ These define a diagonal matrix $\Sigma = \pmatrix{\sigma_1&0&0\\0&\sigma_2&0\\0&0&\sigma_3}.$

2. Generate a random correlation matrix $R = \pmatrix{1&\rho_3&\rho_2\\\rho_3&1&\rho_1\\\rho_2&\rho_1&1}.$

The resulting random covariance is $\Sigma R \Sigma.$ It is symmetric by construction. It will be positive-definite if and only if $R$ is, which is equivalent to $|\rho_3|\le 1,$ $|\rho_2|\le 1,$ and $R$ has positive determinant.

What happens if you generate $(\rho_1,\rho_2,\rho_3)$ using any distribution you like supported on the cube $[0.7,0.9]^3$? The only condition you need to check concerns the determinant. But since

$$\det R = 1 - (\rho_1^2+\rho_2^2+\rho_3^2) + 2\rho_1\rho_2\rho_3,$$

we may do a little bit of Calculus and establish that the minimum value of the determinant is attained when one of the $\rho_i$ equals $0.7$ and the other two equal $0.9,$ with a value of $24/1000\gt 0.$ Consequently

no matter how $\rho_1, \rho_2, \rho_3$ are generated, $\det R$ is always positive. Therefore, provided the $\sigma_i$ are positive, $\Sigma R \Sigma$ is a positive-definite covariance matrix.

---

As an example, you could generate the $\sigma_i^2$ independently with (say) some Gamma distribution and generate the $\rho_i$ uniformly. I created $100,000$ such covariance matrices this way; it took less than two seconds. Here's a summary of the results on which are superimposed the intended distribution density functions, showing the method works as intended.

It is clear that

When $\sigma_1, \ldots, \rho_3$ are drawn from any six-dimensional distribution supported on $(0,\infty)^3\times (0.7,0.9)^3,$ $\Sigma R \Sigma$ is a valid covariance matrix with all correlations between $0.7$ and $0.9.$ Conversely, any distribution of covariance matrices with these properties determines such a distribution of $\sigma_1, \ldots, \rho_3.$

You can even introduce dependencies between the $\sigma_i$ and the $\rho_j$ if you like.

---

This is the R code to reproduce the figure. rcov generates an array of n such covariance matrices (referenced by a third index).

rcov <- function(n=1, shape=1, rate=1) {
  sigma <- matrix(rgamma(3*n, shape, rate), 3)
  rho <- matrix(runif(3*n, 0.7, 0.9), 3)
  array(sapply(1:n, function(i) {
    diag(sigma[,i]) %*% matrix(c(1, rho[3,i], rho[2,i],
                                rho[3,i], 1, rho[1,i],
                                rho[2,i], rho[1,i], 1), 3, 3) %*% diag(sigma[,i])
  }), c(3,3,n))
}

shape <- c(2, 5, 10)
rate <- shape
set.seed(17)
system.time(rho <- apply(Sigma <- rcov(1e5, shape, rate), 3, cov2cor)[c(2, 3, 6), ])

gray <- "#f0f0f0"
par(mfrow=c(1,4))
hist(rho, freq=FALSE, col=gray,
     main=expression(paste("Histogram of all ", rho[i])), xlab="Value")
abline(h=1 / (0.9 - 0.7), lwd=2)
for (i in 1:3) {
  hist(sqrt(Sigma[i,i,]), freq=FALSE, breaks=30, col=gray,
       main=bquote(sigma[.(i)]), xlab="Value")
  curve(dgamma(x, shape[i], rate[i]), lwd=2, add=TRUE)
}
par(mfrow=c(1,1))

Answer 2

The G-Wishart distribution (Letac & Massam, 2007) is a distribution on positive definite matrices with fixed zeros corresponding to the missing edges of a graph $\mathcal G$ with nodes the indices $(i,j)$ of the associated variates. It has a density of the same form as the Wishart distribution: $$p(\Sigma|\delta,\Xi)\propto|\Sigma|^{(\delta-2)/2}\exp\left\{-\frac{1}{2}\text{tr}(\Sigma^\text{T}\Xi)\right\}$$and enjoys the most useful property that the conditional distributions of the submatrices of $\Sigma$ associated with the cliques of the graph all are standard Wishart, which allows for a Gibbs sampling approach to its simulation.

This distribution is implemented in R via the function rgwish. The graph $\mathcal G$ is described by an adjacency upper-triangular matrix adj that is made of 0's and 1's, with 0's indicating the fixed zeroes of the matrix.

In the current question, this R function can be called until all constraints are satisfied. The matrix $\Xi$ (denoted D in rgwish) can be chosen towards favouring the constraints to be met.

Answer 3

Another possible solution is using eigendecomposition. You can use the following code to randomly generate an orthonormal matrix

In R:

Q_matrix = function(d=2){
  require(pracma)
  A = matrix(runif(d^2), d, d)
  Q = gramSchmidt(A)$Q # get orthonormal matrix by Gram Schmidt procedure (s.t. Q%*%t(Q) = I)
  return(Q)
}

Then specify the diagonal matrix containing all the eigenvalues, L, then the covariance matrix will be

L = diag(1:2)
Q = Q_matrix(2)
S = Q%*%L%*%t(Q)

Answer 4

Here is how I would proceed to develop random normal deviates (Y, X, Z) which I would use to compute corresponding covariance matrices.

In the special case of a bivariate multinormal distribution, there is a formula that states for variable Y, which assumes a value of y (so Y=y), there is a conditional univariate normal distribution for X with mean as a function of y and the correlation between X and Y, and a formula for the variance. Here is a reference: https://www2.stat.duke.edu/courses/Spring12/sta104.1/Lectures/Lec22.pdf .

$E[X|Y = y] = µ_X + σ_X ρ\frac{(y − µ_Y)}{σ_Y} $

$Var[X|Y = y] = σ_X (1 − ρ^2) $

Assuming that the largest correlation exist between Y and X and also between Y and Z, then similarly:

$E[Z|Y = y] = µ_Z + σ_Z ρ\frac{(y − µ_Y)}{σ_Y} $

$Var[Z|Y = y] = σ_Z (1 − ρ^2) $

where the supplied correlations are applicable to Y and X, and also Y and Z, respectively.

So, using a package normal random deviate function, you simulate Y and use the value Y=y, to derive conditional univariate normal random deviate function values for X and Z, respectively.

[EDIT] To address a comment, my suggested approach generates data from which to construct a sample random covariance matrix. I would note, to quote Wikipedia:

Another issue is the robustness to outliers, to which sample covariance matrices are highly sensitive.[2][3][4]

However, if you have simulated point data, one may be able to address outliers (another topic) and perhaps relatedly mitigate this issue with your sample population of generated random covariance matrices.

[EDIT EDIT] Using a freely available google spreadsheet, I have a simulated 100 normal, and two sets (of 100) conditional normal random deviates (namely, X, Y and Z) with targeted parameters. Upon valuing the covariance matrix (and the associated correlation matrix), one can build a database of random simulated covariance matrices. Further, I can control the generating process to create a mixture of outliers, to investigate robustness estimates for covariance and correlations.

Here is a rough extract from the spreadsheet, please excuse the formating which is lost on pasting from the spreadsheet into a document text file.

THEORETICAL X Y Z

Mean 0.25 0.1 0.2

SE 0.3 0.2 0.1

CORREL 1.0 0.7 0.5

OBS Mean 0.2177 0.0849 0.2028

OBS SE 0.3377 0.1576 0.0899

OBS VAR 0.1141 0.0081

OBS Correl 1 0.800 0.432

Example of underlying cell formulas employed:

${\text=NORMINV(RAND(),B$3,B$4) , =NORMINV(RAND(),0,C$4*SQRT(1-$C$5*$C$5)+C13 , =C$3+C$4*C$5*($B14-$B$3)/$B$4 , =E$3+E$4*E$5*($B14-$B$3)/$B$4 }$

LIVE FORMULA VARIANCE_COVARIANCE MATRIX
0.1129 0.0538 0.0174
0.0538 0.0400 0.0077
0.0174 0.0077 0.0080

LIVE FORMULA CORRELATION MATRIX
1.0000 0.8003 0.5777
0.8003 1.0000 0.4356
0.5777 0.4356 1.0000

Here is an output example of ranged valued Random Data Matrices produced:

/////////////////////////// DATA RUN DATABASE ////////////

VARIANCE_COVARIANCE MATRIX
0.0973 0.0529 0.0190
0.0529 0.0485 0.0102
0.0190 0.0102 0.0091

CORRELATION MATRIX
1.0000 0.7700 0.6373
0.7700 1.0000 0.5000
0.6373 0.500 1.0000

So far, the random matrice set I have created look pretty rationale without major swings from the targetted parameters. Has anyone seen the output from the Wishart approach, and would it be at all useful in testing robustness of say correlation matrices with a select noise component (as I can readily accomplish)?