# Generating pairs of random variables with given covariance and gamma marginals

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?

• The short answer is yes, and the reason is that the relationship between the correlation coefficient of the copula and the correlation of the Gamma variates depends on the copula.
– whuber
Dec 18, 2019 at 14:33
• Thanks! Can one say something for specific copulas? I'm quite open to using other ones. Dec 18, 2019 at 15:03
• Tautologically, you could create a copula by starting with correlated Gamma marginals and transforming them (with the Probability Integral Transformation)! For any other copula I would expect the analysis to be difficult unless the shape parameters are the same or the scale parameters are the same. Thus, for general purpose use, a numerical procedure to find a solution (given a parameterized family of copulas) would be preferable.
– whuber
Dec 18, 2019 at 17:03

Using Nadarajah & Gupta (2006, Applied Mathematics Letters), I found the following solution, but it works only if $$c > \max\{k_X,k_Y\}$$:

nsims <- 100000
kX <- 2
kY <- 3
sigma_XY <- 7
thetaX <- 2
thetaY <- 3
r <- thetaY/thetaX
c <- kX*kY*thetaX*thetaY/sigma_XY
lX <- c-kX
lY <- c-kY
U <- rbeta(nsims, kX, lX)
V <- rbeta(nsims, kY, lY)
W <- rgamma(nsims, shape = c, scale = thetaX)
X <- U*W
Y <- V*W*r
cov(X,Y) # ~ 7

layout(t(c(2,1)))
plot(ecdf(X), col = "red")
curve(pgamma(x, shape = kX, scale = thetaX),
add = TRUE, lty = "dashed", lwd = 3)
plot(ecdf(Y), col = "red")
curve(pgamma(x, shape = kY, scale = thetaY),
add = TRUE, lty = "dashed", lwd = 3)


• Thanks! This is a nice non-approximate solution for special cases (the most important limitation being that it only works for positively correlated data). It doesn't seem like there will be much more forthcoming, so I'll accept this. Thanks again! Dec 28, 2019 at 7:20