PDF of the product of two independent Gamma random variables

I am trying to find out the density of the product $XY$ of two independent Gamma random variables $X \sim \mathrm{Gamma}(k_1, \theta_1)$ and $Y \sim \mathrm{Gamma}(k_2, \theta_2)$, where $k_i$'s are the shape parameters and $\theta_i$'s are the scale parameters.

I know the formula for the density of the product of two independent RV https://en.wikipedia.org/wiki/Product_distribution, but I couldn't evaluate the integral.

I know the resulting density will involve Meijer G-functions, and I found this paper: http://epubs.siam.org/doi/abs/10.1137/0118065, but the Gamma distribution density in this paper is given as $$f_i(x_i)=\frac{1}{\Gamma{\left(b_i\right)}} x_i^{b_i-1} e^{-x_i}$$ with unity scale parameter.

Is there a formula to find the density of the product of two independent Gamma random variables with arbitrary shape and scale parameters?

Any help is highly appreciated.

• What is amma RV? Commented Aug 17, 2017 at 2:24
• Just in case someone stumbles on this in the future and needs the product of a Gamma rv and a power of a second Gamma, i.e., $X^n Y$, that's where the Meijer G functions are involved and apparently can't be simplified the way they can for $XY$. See math.stackexchange.com/questions/4273438/… for the result of the integral for specific $n$. I haven't yet filled in the details so if you need the density of $X^n Y$ there's some more work to do but that should get you close! Commented Oct 12, 2021 at 3:42

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 :)

• This is perfect! Thank you so much :) Can you point to any documentation about how this function TransformProduct works? Commented Aug 17, 2017 at 23:18
• I should ask this in a new question, but in case you see this comment—can mathStatica handle the generalization where $Z = X^n Y$? That is, instead of the product of two Gamma random variables, you raise one of them to an integer power $n$ and then multiply by the other? $X^n$ will be Generalized-Gamma-distributed. My efforts to simplify the result by hand/Wolfram Alpha are failing and I wonder how mathStatica would fare. (I don't even need the exact density; up to a normalizing constant is enough.) Commented Oct 11, 2021 at 6:41

In order to provide an analytical derivation and confirm the result, it is possible to proceed as follows: Let's start by defining independent gamma product $$Z_{2}=XY$$, $$X=Z_{1}$$ and $$Y=\frac{Z_{2}}{Z_{1}}$$. Then we can write joint density using jacobian, $$$$f(z_{1},z_{2})=f_{x}(z_{1})f_{y}(z_{2})|J|$$$$ where jacobian $$|J|$$ is, $$\begin{bmatrix}\frac{\partial x }{\partial z_{1}}&\frac{\partial x }{\partial z_{2}}\\\frac{\partial y }{\partial z_{1}}&\frac{\partial y }{\partial z_{2}}\end{bmatrix}$$ $$\begin{bmatrix}1&0\\\frac{-z_{2}}{z_{1}^{2}}&\frac{1}{z_{1}}\end{bmatrix}$$ Then this equals, $$$$f(z_{1},z_{2})=f_{x}(z_{1})f_{y}(z_{2})\frac{1}{z_{1}}$$$$ $$$$f(z_{1},z_{2})=\frac{\beta_{x}^{-\alpha_{x}}z_{1}^{\alpha_{x}-1}e^{\frac{-z_{1}}{\beta_{x}}}}{\Gamma(\beta_{x})}\frac{\beta_{y}^{-\alpha_{y}}(\frac{z_{2}}{z_{1}})^{\alpha_{y}-1}e^{\frac{-z_{2}/z_{1}}{\beta_{y}}}}{\Gamma(\beta_{y})}$$$$ In order to get marginal density of $$Z_{2}=XY$$ the product, we integrate out the joint density we wrote over $$Z_{1}$$

$$$$f(z_{2})=\beta_{x}^{-\alpha_{x}}\beta_{y}^{-\alpha_{y}}z_{2}^{\alpha_{y}-1}\int_{0}^{\infty}\left[\frac{z_{1}^{\alpha_{x}-\alpha_{y}-1}e^{\frac{-z_{1}}{\beta_{x}}}e^{\frac{-z_{2}/z_{1}}{\beta_{y}}}}{\Gamma(\beta_{x})\Gamma(\beta_{y})}dz_{1}\right]$$$$ Using the fact that, $$$$\int_{0}^{\infty}x^{\alpha-1}\exp(-ax-bx^{-1})dx=2\left(\frac{b}{a}\right)^{\alpha/2}K_{\alpha}\left(2\sqrt{ab} \right)$$$$ We will have, $$$$f\left(z_{2}=xy\right)=\frac{z_{2}^{\alpha_{y}-1}\beta_{x}^{-\alpha_{x}}\beta_{y}^{-\alpha_{y}}2\left(\frac{z_{2}\beta_{x}}{\beta_{y}}\right)^{\frac{\left(\alpha_{x}-\alpha_{y}\right)}{2}}K_{\alpha_{x}-\alpha_{y}}\left(2\sqrt{\frac{z_{2}}{\beta_{x}\beta_{y}}} \right)}{\Gamma(\alpha_{x})\Gamma(\alpha_{y})}$$$$

Another Monte-Carlo check confirms @wolfies nice result above,

Plus an integral on the pdf confirms the accuracy of it to be the correct density. After integrating in R we obtain,

integrate(gamprod,lower = 0.,upper = Inf,ax = 2,bx = 3,ay = 4,by = 1.1) 1 with absolute error < 6.8e-05