Simulation of PDF of sum of correlated Gamma random variables (in R)

Question

My question is very related to the general sum of Gamma RVs question found in the following link:

[The sum of two independent gamma random variables

There is some helpful R code there for generating a distribution function for a sum of Gamma random variables. However, I am wondering if there is a more general form that allows for arbitrary correlations between the variables.

Specifically, suppose I have:

$X_i$ ~ $Gamma(\alpha_i,\beta_i), \;\;\;\;\;\; i=1,...,n$

but $cov(X_i,X_j)=\rho_{ij}\neq 0$

Can the R code below be augmented to account for this in simulating the empirical distribution of $\Sigma_{i=1}^n X_i$? If these were normally distributed it would be easy, but I cannot find any generalization for Gamma random variables.

shape <- 1:3 #ki
scale <- 1:3 # thetai
make_cumgenfun  <-  function(shape, scale) {
  # we return list(shape, scale, K, K', K'')
  n  <-  length(shape)
  m <-   length(scale)
  stopifnot( n == m, shape > 0, scale > 0 )
  return( list( shape=shape,  scale=scale, 
                Vectorize(function(s) {-sum(shape * log(1-scale * s) ) }),
                Vectorize(function(s) {sum((shape*scale)/(1-s*scale))}) ,
                Vectorize(function(s) { sum(shape*scale*scale/     (1-s*scale)) }))    )
}

solve_speq  <-  function(x, cumgenfun) {
      # Returns saddle point!
      shape <- cumgenfun[[1]]
      scale <- cumgenfun[[2]]
      Kd  <-   cumgenfun[[4]]
      uniroot(function(s) Kd(s)-x,lower=-100,
              upper = 0.3333, 
              extendInt = "upX")$root
}

make_fhat <-  function(shape,  scale) {
cgf1  <-  make_cumgenfun(shape, scale)
K  <-  cgf1[[3]]
Kd <-  cgf1[[4]]
Kdd <- cgf1[[5]]
# Function finding fhat for one specific x:
fhat0  <- function(x) {
    # Solve saddlepoint equation:
    s  <-  solve_speq(x, cgf1)
    # Calculating saddlepoint density value:
    (1/sqrt(2*pi*Kdd(s)))*exp(K(s)-s*x)
}
# Returning a vectorized version:
return(Vectorize(fhat0))
} #end make_fhat

 fhat  <-  make_fhat(shape, scale)
plot(fhat, from=0.01,  to=40, col="red", main="unnormalized saddlepoint approximation\nto sum of three gamma variables")

"If these were normally distributed it would be easy," This is because if $X_i $ are (marginally) normal and the correlation coefficients are known, then there is a _unique_ multivariate normal density that has the specified marginal densities and the specified correlations. But, given Gamma random variables (or most other kinds of variables too), there is no unique joint density that has the desired marginals and the desired correlations, and we need the joint density to get a reasonably accurate simulation. — Dilip Sarwate, Nov 14 '15 at 03:13
@DilipSarwate, so there is nothing that can be done? What if I have an estimate of the correlation between the $X_i$? More generally, how could I model the sum of correlated Gamma RVs...surely it cant be "impossible"...even if it is imperfect. — Justin, Nov 14 '15 at 03:17
It is not whether the correlation is estimated or known. If $X_i\sim\Gamma(t_i,\lambda_i),i=1,2$ and $\rho$ is given, we cannot say _which_ joint density $f_{X_1,X_2}(x_1,x_2)$ we are working with; there are _numerous_ possible joint densities that could have produced these correlated Gamma random variables. _Which_ one do you want to simulate? and what makes this choice better or more suitable than all other possibilities? In contrast, for normal random variables, there is a unique bivariate normal density for given means, variances and correlation coefficient, and we can simulate from that. — Dilip Sarwate, Nov 14 '15 at 04:25
There is a (not "the") practical solution using Gaussian copulas: If $G(x;t_i,\lambda_i)$ are your Ga$(t_i,\lambda_i)$ cdfs, consider the Gaussian vector of the $\Phi^{-1}(G(X_i;t_i,\lambda_i))$ and determine by trial-and-error the correlation matrix of this vector for the $X_i$ to have the appropriate correlations. There is no close form formula for this. (A keyword is _copula_.) — Xi'an, Nov 14 '15 at 09:16
@Xi'an, I just started looking toward copulas, but I'm admittedly not fully there with understanding. I apologize if this is overly tedious, but could you lay out a more detailed way to approach this? What if I use the basic Gaussian copula? — Justin, Nov 15 '15 at 17:37
Maybe something could be done if you could assume some multivariate gamma distribution. — kjetil b halvorsen, Apr 13 '18 at 21:39

score 1 · Answer 1 · answered Dec 18 '18 at 21:24

Partially answered in comments:

If these were normally distributed it would be easy, This is because if $X_i$ are (marginally) normal and the correlation coefficients are known, then there is a unique multivariate normal density that has the specified marginal densities and the specified correlations. But, given Gamma random variables (or most other kinds of variables too), there is no unique joint density that has the desired marginals and the desired correlations, and we need the joint density to get a reasonably accurate simulation. – Dilip Sarwate

( @DilipSarwate, so there is nothing that can be done? What if I have an estimate of the correlation between the $X_i$? More generally, how could I model the sum of correlated Gamma RVs ... surely it cant be impossible ... even if it is imperfect. – Justin )

It is not whether the correlation is estimated or known. If $Xi∼Γ(t_i,λ_i),i=1,2 $and $ρ$ is given, we cannot say which joint density $f_{X1,X2}(x1,x2)$ we are working with; there are numerous possible joint densities that could have produced these correlated Gamma random variables. Which one do you want to simulate? and what makes this choice better or more suitable than all other possibilities? In contrast, for normal random variables, there is a unique bivariate normal density for given means, variances and correlation coefficient, and we can simulate from that. – Dilip Sarwate

Another point is that some multivariate gamma distribution might not even be determined by all bivariate correlations, more information might be needed.

There is a (not "the") practical solution using Gaussian copulas: If $G(x;t_i,λ_i)$ are your $Ga(t_i,λ_i)$ cdfs, consider the Gaussian vector of the $Φ^{−1}(G(X_i;t_i,λ_i))$ and determine by trial-and-error the correlation matrix of this vector for the $X_i$ to have the appropriate correlations. There is no close form formula for this. (A keyword is copula.) – Xi'an

Some other copula than the gaussian could also be used ---

A simple simulation example with a bivariate gaussian copula in R: How to simulate from a Gaussian copula?

Simulation of PDF of sum of correlated Gamma random variables (in R)

1 Answers1