Simulation of PDF of sum of correlated Gamma random variables (in R) 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")

 A: 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? 
