Difference of Gamma random variables

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 $(-\infty,\infty)$
  and so $f_{X-Y}(z)$ has support $(-\infty,\infty)$.

• 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,
  $$f_{X+Y}(z) = \int_{-\infty}^\infty f_X(x)f_Y(z-x)\,\mathrm dx$$
  and that the density of the random variable $-Y$ is
  $f_{-Y}(\alpha) = f_Y(-\alpha)$, deduce that
  $$f_{X-Y}(z) = f_{X+(-Y)}(z) = \int_{-\infty}^\infty f_X(x)f_{-Y}(z-x)\,\mathrm dx = \int_{-\infty}^\infty f_X(x)f_Y(x-z)\,\mathrm dx.$$

• Third, for non-negative random variables $X$ and $Y$, note that the
  above expression simplifies to
  $$f_{X-Y}(z) = \begin{cases} \int_0^\infty f_X(x)f_Y(x-z)\,\mathrm dx, & z < 0,\\ \int_{0}^\infty f_X(y+z)f_Y(y)\,\mathrm dy, &#038; z > 0. \end{cases}$$

• Finally, using parametrization $\Gamma(s,\lambda)$ to mean a
  random variable with density
  $\lambda\frac{(\lambda x)^{s-1}}{\Gamma(s)}\exp(-\lambda x)\mathbf 1_{x>0}(x)$,
  and with
  $X \sim \Gamma(s,\lambda)$ and $Y \sim \Gamma(t,\mu)$ random variables,
  we have for $z > 0$ that
  $$\begin{align*}f_{X-Y}(z) &= \int_{0}^\infty \lambda\frac{(\lambda (y+z))^{s-1}}{\Gamma(s)}\exp(-\lambda (y+z)) \mu\frac{(\mu y)^{t-1}}{\Gamma(t)}\exp(-\mu y)\,\mathrm dy\\ &= \exp(-\lambda z) \int_0^\infty p(y,z)\exp(-(\lambda+\mu)y)\,\mathrm dy.\tag{1} \end{align*}$$
  Similarly, for $z < 0$,
  $$\begin{align*}f_{X-Y}(z) &= \int_{0}^\infty \lambda\frac{(\lambda x)^{s-1}}{\Gamma(s)}\exp(-\lambda x) \mu\frac{(\mu (x-z))^{t-1}}{\Gamma(t)}\exp(-\mu (x-z))\,\mathrm dx\\ &= \exp(\mu z) \int_0^\infty q(x,z)\exp(-(\lambda+\mu)x)\,\mathrm dx.\tag{2} \end{align*}$$

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
$$\int_0^\infty x^{s-1}(x+\beta)^{s-1}\exp(-\nu x)\,\mathrm dx$$
in terms of polynomial, exponential and Bessel functions of $\beta$
and this can be used to write down explicit expressions for $f_{X-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, z^2, \ldots z^{s-1}$. It follows that the density of
  $X-Y$ is proportional to a mixture density of
  $\Gamma(1,\lambda), \Gamma(2,\lambda), \cdots, \Gamma(s,\lambda)$
  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
  $\Gamma(1,\mu), \Gamma(2,\mu), \cdots, \Gamma(t,\mu)$
  random variables flipped over, that is,
  it will have terms such as $(\mu|z|)^{k-1}\exp(\mu z)$
  instead of the usual $(\mu z)^{k-1}\exp(-\mu z)$.
  Also, this result will hold even if $s$ is not an integer.