Given two independent random variables X∼Gamma(αX,βX) and Y∼Gamma(αY,βY), what is the distribution of the difference, i.e. D=X−Y?
If the result is not wellknown, how would I go about deriving the result?
Answer
I will outline how the problem can be approached and state
what I think the end result will be for the special case
when the shape parameters are integers, but not fill in the
details.

First, note that X−Y takes on values in (−∞,∞)
and so fX−Y(z) has support (−∞,∞). 
Second, from the standard results that the
density of the sum of two independent continuous random variables is the
convolution of their densities, that is,
fX+Y(z)=∫∞−∞fX(x)fY(z−x)dx
and that the density of the random variable −Y is
f−Y(α)=fY(−α), deduce that
fX−Y(z)=fX+(−Y)(z)=∫∞−∞fX(x)f−Y(z−x)dx=∫∞−∞fX(x)fY(x−z)dx. 
Third, for nonnegative random variables X and Y, note that the
above expression simplifies to
fX−Y(z)={∫∞0fX(x)fY(x−z)dx,z<0,∫∞0fX(y+z)fY(y)dy,z>0. 
Finally, using parametrization Γ(s,λ) to mean a
random variable with density
λ(λx)s−1Γ(s)exp(−λx)1x>0(x),
and with
X∼Γ(s,λ) and Y∼Γ(t,μ) random variables,
we have for z>0 that
fX−Y(z)=∫∞0λ(λ(y+z))s−1Γ(s)exp(−λ(y+z))μ(μy)t−1Γ(t)exp(−μy)dy=exp(−λz)∫∞0p(y,z)exp(−(λ+μ)y)dy.
Similarly, for z<0,
fX−Y(z)=∫∞0λ(λx)s−1Γ(s)exp(−λx)μ(μ(x−z))t−1Γ(t)exp(−μ(x−z))dx=exp(μz)∫∞0q(x,z)exp(−(λ+μ)x)dx.
These integrals are not easy to evaluate but for the special case
s=t, Gradshteyn and Ryzhik, Tables of Integrals, Series, and Products,
Section 3.383, lists the value of
∫∞0xs−1(x+β)s−1exp(−νx)dx
in terms of polynomial, exponential and Bessel functions of β
and this can be used to write down explicit expressions for fX−Y(z).
From here on, we assume that s and t are integers so
that p(y,z) is a polynomial in y and z of degree (s+t−2,s−1)
and q(x,z) is a polynomial in x and z of degree (s+t−2,t−1).

For z>0, the integral (1)
is the sum of s Gamma integrals with respect to y with coefficients
1,z,z2,…zs−1. It follows that the density of
X−Y is proportional to a mixture density of
Γ(1,λ),Γ(2,λ),⋯,Γ(s,λ)
random variables for z>0. Note that this result
will hold even if t is not an integer. 
Similarly, for z<0,
the density of
X−Y is proportional to a mixture density of
Γ(1,μ),Γ(2,μ),⋯,Γ(t,μ)
random variables flipped over, that is,
it will have terms such as (μz)k−1exp(μz)
instead of the usual (μz)k−1exp(−μz).
Also, this result will hold even if s is not an integer.
Attribution
Source : Link , Question Author : FBC , Answer Author : Dilip Sarwate