You can create whole families of joint distributions on $(X,Y)$ such that $X\sim \Gamma(a_1,\alpha)$ and $Y\sim \Gamma(a_2,\beta)$ by using copulas like
$$
F_{(X,Y)}(x,y) = \mathbb{P}(X\le x,Y\le y) =
\dfrac{F_X(x)F_Y(y)}{1+\varrho (1-F_X(x))(1-F_Y(y))}
$$
for $-1\le \varrho \le 1$. The joint distribution is continuous, which means the event $X=Y$ has probability zero.
Now, if you have a specific reason for using McKay's bivariate distribution, with pdf
$$
f_{(X,Y)}(x,y) = \alpha^{p+q} x^{p-1} (y-x)^{q-1} \exp\{-\alpha y\} / [\Gamma(p) \Gamma(q)]\,\mathbb{I}_{0\le x\le y} \,,
$$
which gives
$$
X\sim \Gamma(p,\alpha)\,,\quad Y\sim \Gamma(p+q,\alpha)
$$
as marginals,
you must compute $\mathbb{E}[G(X,Y)]$ as
$$
\int_0^\infty \int_0^y G(x,y)\,\alpha^{p+q} x^{p-1} (y-x)^{q-1} \exp\{-\alpha y\} / [\Gamma(p) \Gamma(q)]\,\text{d}x\,\text{d}y\,.
$$
Since you know the density, you can just use fitdistr.
# Sample data
library(LaplacesDemon)
x <- rinvgamma(1000, 1,2)
library(MASS)
f <- function(x, rho, a, s)
1/(a*gamma(rho)) * (a / (x-s))^(rho+1) * exp( - a/(x-s) )
fitdistr( x, f, list(rho=1, a=1, s=0) )
Best Answer
You can find other versions of the bivariate gamma distribution in the following paper by Nadarajah:
http://www.emis.de/journals/HOA/MPE/2005/2151.pdf