I have shape parameters $k_X, k_Y$ and scale parameters $\theta_X, \theta_Y$, as well as a covariance $\sigma_{XY}$.
How do I generate random variables $(X,Y)$ such that the marginals are gamma distributed, $X\sim\Gamma(k_X,\theta_X)$ and $Y\sim\Gamma(k_Y,\theta_Y)$, and the covariance matrix satisfies $$ \Sigma=\begin{pmatrix} k_X\theta_X^2 & \sigma_{XY} \\ \sigma_{XY} & k_Y\theta_Y^2 \end{pmatrix}? $$
Any help involving R would be most appreciated.
Judging from earlier answers here and here, one of which links to this R-bloggers post, the canonical idea seems to be to use copulas, which allows a certain control of the correlation (and therefore the covariance, since we know the marginal variances).
However, the correlations we "get out" of the copula random generation are different than the ones we "put in". For instance, "feeding in" $\rho_{\text{Normal}}=0.80$ yields $\rho_{\text{Gamma}}\approx 0.78$, and "feeding in" $\rho_{\text{Normal}}=-0.80$ yields $\rho_{\text{Gamma}}\approx -0.66$:
library(copula)
nn_sim <- 1e6
set.seed(1)
my_copula <- normalCopula(0.8)
my_population <- mvdc(my_copula, margins=rep("gamma",2),
paramMargins=list(list(shape=2,scale=1),list(shape=2,scale=1)))
sims <- rMvdc(nn_sim,my_population)
cor(sims)
# [,1] [,2]
# [1,] 1.0000000 0.7840255
# [2,] 0.7840255 1.0000000
set.seed(1)
my_copula <- normalCopula(-0.8)
my_population <- mvdc(my_copula, margins=rep("gamma",2),
paramMargins=list(list(shape=2,scale=1),list(shape=2,scale=1)))
sims <- rMvdc(nn_sim,my_population)
cor(sims)
# [,1] [,2]
# [1,] 1.0000000 -0.6560454
# [2,] -0.6560454 1.0000000
This seems to be a common effect in copula RNGs.
Do I need to play around with the correlation I "put into" the copula RNG until I get the covariance I actually want, or is there a way to set this "input" correlation beforehand? Or is there some completely different way?
