PDF of the product of two independent Gamma random variables

Let $X \sim \text{Gamma}(a,b)$ with pdf $f(x)$:

and $Y \sim \text{Gamma}(\alpha,\beta)$ be independent with pdf $g(y)$:

Then, the pdf of the product $Z = X Y$ can be obtained as $h(z)$:

where I am using the TransformProduct function from mathStatica/Mathematica to do the nitty-gritties, and BesselK[n,z] denotes the modified Bessel function of the second kind. This is much simpler than requiring MeijerG functions. I should note that I am one of the authors of the software function used.

Quick Monte Carlo check

• against theoretical solution derived above when $a =2$, $b = 3$, $\alpha = 4$, and $\beta = 1.1$

Looks fine :)