# Distribution of product of two independent inverse gaussian random variables

Let us consider the case of $$x,y \sim\text{i.i.d. IG}(\mu, \lambda)$$ - Inverse Gaussian Distribution (aka Wald distribution) with the pdf:

$$$$f(x)=\frac{\sqrt{\lambda}}{\sqrt{2 \pi x^3}} \exp\left(- \frac{\lambda(x-\mu)^2}{2 \mu^2 x} \right)$$$$

My question: What will be the distribution of $$z=x \times y$$?

I've tried taking the integral in order to get to the product distribution but I couldn't move further than this:

$$$$f_z(z) = \int_{0}^\infty \frac{\sqrt{\lambda}}{\sqrt{2 \pi x^3}} \exp\left(- \frac{\lambda(x-\mu)^2}{2 \mu^2 x} \right) \times \frac{\sqrt{\lambda}}{\sqrt{2 \pi \frac{z^3}{x^3}}} \exp\left(- \frac{\lambda \left(\frac{z}{x}-\mu \right)^2}{2 \mu^2 \frac{z}{x}} \right)\frac{1}{x} dx = \\ \int_{0}^\infty \frac{\lambda}{2 \pi \sqrt{z^3}} \exp\left(- \frac{\lambda z (x-\mu)^2 + \lambda x^2 \left(\frac{z}{x}-\mu\right)^2}{2 \mu^2 z x} \right) \frac{1}{x} dx = \\ \frac{\lambda}{2 \pi z^\frac{3}{2}} \exp \left( \frac{2\lambda}{\mu} \right) \int_{0}^\infty \frac{1}{x} \exp \left(- \lambda \frac{(x^2 + z)(z+\mu^2)}{2 \mu^2 x z} \right) dx$$$$

I have a feeling that this might have another IG distribution or something close to it, maybe something like: IG$$\left( \mu^2, \frac{\lambda^4}{\mu^2(\mu^4 + 2 \lambda^2)} \right)$$ - these were derived based on the mean and variance of products of independent random variables.

Is there an analytical solution for this or do I need to use numerical ones?

$$$$f_z(z) = \frac{\lambda}{\pi z^{\frac{3}{2}}} \exp\left(\frac{2\lambda}{\mu}\right) K_0 \left(\frac{\lambda(z+\mu^2)}{\sqrt{z}\mu^2}\right),$$$$
where $$K_0(x)=\int_0^\infty\frac{cos(xt)}{\sqrt{t^2+1}}dt$$ is the modified Bessel function of the third kind.
This is obtained after regrouping the elements in the exponent of the original question and noticing that the PDF looks like a product distribution of two variables following i.i.d. exponential distribution with the rate $$\frac{\lambda(z+\mu^2)}{2 z \mu^2}$$.