how to generate specific random covariance matrices? 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…
 A: 
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.
A: In this particular case there's a simple, easy, completely general method.
Break the problem down into two parts:

*

*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}.$


*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))

