Product distribution

We first write the cumulative distribution function of ${\displaystyle Z}$ starting with its definition

${\displaystyle {\begin{aligned}F_{Z}(z)&\,{\stackrel {\text{def}}{=}}\ \mathbb {P} (Z\leq z)\\&=\mathbb {P} (XY\leq z)\\&=\mathbb {P} (XY\leq z,X\geq 0)+\mathbb {P} (XY\leq z,X\leq 0)\\&=\mathbb {P} (Y\leq z/X,X\geq 0)+\mathbb {P} (Y\geq z/X,X\leq 0)\\&=\int _{0}^{\infty }f_{X}\left(x\right)\int _{-\infty }^{z/x}f_{Y}\left(y\right)\,dy\,dx+\int _{-\infty }^{0}f_{X}\left(x\right)\int _{z/x}^{\infty }f_{Y}\left(y\right)\,dy\,dx\end{aligned}}}$

We find the desired probability density function by taking the derivative of both sides with respect to ${\displaystyle z}$. Since on the right hand side, ${\displaystyle z}$ appears only in the integration limits, the derivative is easily performed using the fundamental theorem of calculus and the chain rule. (Note the negative sign that is needed when the variable occurs in the lower limit of the integration.)

${\displaystyle {\begin{aligned}f_{Z}(z)&=\int _{0}^{\infty }f_{X}(x)f_{Y}(z/x){\frac {1}{x}}\,dx-\int _{-\infty }^{0}f_{X}(x)f_{Y}(z/x){\frac {1}{x}}\,dx\\&=\int _{0}^{\infty }f_{X}(x)f_{Y}(z/x){\frac {1}{|x|}}\,dx+\int _{-\infty }^{0}f_{X}(x)f_{Y}(z/x){\frac {1}{|x|}}\,dx\\&=\int _{-\infty }^{\infty }f_{X}(x)f_{Y}(z/x){\frac {1}{|x|}}\,dx.\end{aligned}}}$

where the absolute value is used to conveniently combine the two terms.

A faster more compact proof begins with the same step of writing the cumulative distribution of ${\displaystyle Z}$ starting with its definition:

${\displaystyle {\begin{aligned}F_{Z}(z)&{\overset {\underset {\mathrm {def} }{}}{=}}\ \ \mathbb {P} (Z\leq z)\\&=\mathbb {P} (XY\leq z)\\&=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f_{X}\left(x\right)f_{Y}\left(y\right)u\left(z-xy\right)\,dy\,dx\end{aligned}}}$

where ${\displaystyle u\left(.\right)}$ is the Heaviside step function and serves to limit the region of integration to values of ${\displaystyle x}$ and ${\displaystyle y}$ satisfying ${\displaystyle xy\leq z}$.

We find the desired probability density function by taking the derivative of both sides with respect to ${\displaystyle z}$.

${\displaystyle {\begin{aligned}f_{Z}(z)&=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f_{X}\left(x\right)f_{Y}\left(y\right)\delta (z-xy)\,dy\,dx\\&=\int _{-\infty }^{\infty }f_{X}\left(x\right)f_{Y}\left(z/x\right)\left[\int _{-\infty }^{\infty }\delta (z-xy)\,dy\right]\,dx\\&=\int _{-\infty }^{\infty }f_{X}\left(x\right)f_{Y}\left(z/x\right){\frac {1}{|x|}}\,dx.\end{aligned}}}$

where we utilize the translation and scaling properties of the Dirac delta function ${\displaystyle \delta }$.

A more intuitive description of the procedure is illustrated in the figure below. The joint pdf ${\displaystyle f_{X}(x)f_{Y}(y)}$ exists in the ${\displaystyle x}$-${\displaystyle y}$ plane and an arc of constant ${\displaystyle z}$ value is shown as the shaded line. To find the marginal probability ${\displaystyle f_{Z}(z)}$ on this arc, integrate over increments of area ${\displaystyle dx\,dy\;f(x,y)}$ on this contour.

Diagram to illustrate the product distribution of two variables.

Starting with ${\displaystyle y={\frac {z}{x}}}$, we have ${\displaystyle dy=-{\frac {z}{x^{2}}}\,dx=-{\frac {y}{x}}\,dx}$. So the probability increment is ${\displaystyle \delta p=f(x,y)\,dx\,|dy|=f_{X}(x)f_{Y}(z/x){\frac {y}{|x|}}\,dx\,dx}$. Since ${\displaystyle z=yx}$ implies ${\displaystyle dz=y\,dx}$, we can relate the probability increment to the ${\displaystyle z}$-increment, namely ${\displaystyle \delta p=f_{X}(x)f_{Y}(z/x){\frac {1}{|x|}}\,dx\,dz}$. Then integration over ${\displaystyle x}$, yields ${\displaystyle f_{Z}(z)=\int f_{X}(x)f_{Y}(z/x){\frac {1}{|x|}}\,dx}$.

Let ${\displaystyle X\sim f(x)}$ be a random sample drawn from probability distribution ${\displaystyle f_{x}(x)}$. Scaling ${\displaystyle X}$ by ${\displaystyle \theta }$ generates a sample from scaled distribution ${\displaystyle \theta X\sim {\frac {1}{|\theta |}}f_{x}\left({\frac {x}{\theta }}\right)}$ which can be written as a conditional distribution ${\displaystyle g_{x}(x|\theta )={\frac {1}{|\theta |}}f_{x}\left({\frac {x}{\theta }}\right)}$.

Letting ${\displaystyle \theta }$ be a random variable with pdf ${\displaystyle f_{\theta }(\theta )}$, the distribution of the scaled sample becomes ${\displaystyle f_{X}(\theta x)=g_{X}(x\mid \theta )f_{\theta }(\theta )}$ and integrating out ${\displaystyle \theta }$ we get ${\displaystyle h_{x}(x)=\int _{-\infty }^{\infty }g_{X}(x|\theta )f_{\theta }(\theta )d\theta }$ so ${\displaystyle \theta X}$ is drawn from this distribution ${\displaystyle \theta X\sim h_{X}(x)}$. However, substituting the definition of ${\displaystyle g}$ we also have ${\displaystyle h_{X}(x)=\int _{-\infty }^{\infty }{\frac {1}{|\theta |}}f_{x}\left({\frac {x}{\theta }}\right)f_{\theta }(\theta )\,d\theta }$ which has the same form as the product distribution above. Thus the Bayesian posterior distribution ${\displaystyle h_{X}(x)}$ is the distribution of the product of the two independent random samples ${\displaystyle \theta }$ and ${\displaystyle X}$.

For the case of one variable being discrete, let ${\displaystyle \theta }$ have probability ${\displaystyle P_{i}}$ at levels ${\displaystyle \theta _{i}}$ with ${\displaystyle \sum _{i}P_{i}=1}$. The conditional density is ${\displaystyle f_{X}(x\mid \theta _{i})={\frac {1}{|\theta _{i}|}}f_{x}\left({\frac {x}{\theta _{i}}}\right)}$. Therefore ${\displaystyle f_{X}(\theta x)=\sum {\frac {P_{i}}{|\theta _{i}|}}f_{X}\left({\frac {x}{\theta _{i}}}\right)}$.

The distribution of the product of two random variables which have lognormal distributions is again lognormal. This is itself a special case of a more general set of results where the logarithm of the product can be written as the sum of the logarithms. Thus, in cases where a simple result can be found in the list of convolutions of probability distributions, where the distributions to be convolved are those of the logarithms of the components of the product, the result might be transformed to provide the distribution of the product. However this approach is only useful where the logarithms of the components of the product are in some standard families of distributions.

Let ${\displaystyle Z}$ the product of two independent variables ${\displaystyle Z=X_{1}X_{2}}$ each uniformly distributed on the interval [0,1], possibly the outcome of a copula transformation. As noted in "Lognormal Distributions" above, PDF convolution operations in the Log domain correspond to the product of sample values in the original domain. Thus, making the transformation ${\displaystyle u=\ln(x)}$, such that ${\displaystyle p_{U}(u)\,|du|=p_{X}(x)\,|dx|}$, each variate is distributed independently on u as

${\displaystyle p_{U}(u)=p_{X}(x)/|du/dx|={\frac {1}{x^{-1}}}=e^{u},\;\;-\infty <u\leq 0}$.

and the convolution of the two distributions is the autoconvolution

${\displaystyle c(y)=\int _{u=0}^{y}e^{u}e^{y-u}du=-\int _{u=y}^{0}e^{y}du=-ye^{y},\;\;-\infty <y\leq 0}$

Next retransform the variable to ${\displaystyle z=e^{y}}$ yielding the distribution

${\displaystyle c_{2}(z)=c_{Y}(y)/|dz/dy|={\frac {-ye^{y}}{e^{y}}}=-y=\ln(1/z)}$ on the interval [0,1]

For the product of multiple ( >2 ) independent samples the characteristic function route is favorable. If we define ${\displaystyle {\tilde {y}}=-y}$ then ${\displaystyle c({\tilde {y}})}$ above is a Gamma distribution of shape 1 and scale factor 1, ${\displaystyle c({\tilde {y}})={\tilde {y}}e^{-{\tilde {y}}}}$ , and its known CF is ${\displaystyle (1-it)^{-1}}$. Note that ${\displaystyle |d{\tilde {y}}|=|dy|}$ so the Jacobian of the transformation is unity.

The convolution of ${\displaystyle n}$ independent samples from ${\displaystyle {\tilde {Y}}}$ therefore has CF ${\displaystyle (1-it)^{-n}}$ which is known to be the CF of a Gamma distribution of shape ${\displaystyle n}$:

${\displaystyle c_{n}({\tilde {y}})=\Gamma (n)^{-1}{\tilde {y}}^{(n-1)}e^{-{\tilde {y}}}=\Gamma (n)^{-1}(-y)^{(n-1)}e^{y}}$.

Making the inverse transformation ${\displaystyle z=e^{y}}$ we get the PDF of the product of the n samples:

${\displaystyle f_{n}(z)={\frac {c_{n}(y)}{|dz/dy|}}=\Gamma (n)^{-1}{\Big (}-\log z{\Big )}^{n-1}e^{y}/e^{y}={\frac {{\Big (}-\log z{\Big )}^{n-1}}{(n-1)!\;\;\;}}}$

The following, more conventional, derivation from Stackexchange[6] is consistent with this result. First of all, letting ${\displaystyle Z_{2}=X_{1}X_{2}}$ its CDF is

${\displaystyle {\begin{aligned}F_{Z_{2}}(z)=\Pr {\Big [}Z_{2}\leq z{\Big ]}&=\int _{x=0}^{1}\Pr {\Big [}X_{2}\leq {\frac {z}{x}}{\Big ]}f_{X_{1}}(x)\,dx\\&=\int _{x=0}^{z}dx+\int _{x=z}^{1}{\frac {z}{x}}\,dx\\&=z-z\log z,\;\;0<z\leq 1\end{aligned}}}$

The density of ${\displaystyle z_{2}{\text{ is then }}f(z_{2})=-\log(z_{2})}$

Multiplying by a third independent sample gives distribution function

${\displaystyle {\begin{aligned}F_{Z_{3}}(z)=\Pr {\Big [}Z_{3}\leq z{\Big ]}&=\int _{x=0}^{1}\Pr {\Big [}X_{3}\leq {\frac {z}{x}}{\Big ]}f_{Z_{2}}(x)\,dx\\&=-\int _{x=0}^{z}\log(x)\,dx-\int _{x=z}^{1}{\frac {z}{x}}\log(x)\,dx\\&=-z{\Big (}\log(z)-1{\Big )}+{\frac {1}{2}}z\log ^{2}(z)\end{aligned}}}$

Taking the derivative yields ${\displaystyle f_{Z_{3}}(z)={\frac {1}{2}}\log ^{2}(z),\;\;0<z\leq 1.}$

The author of the note conjectures that, in general, ${\displaystyle f_{Z_{n}}(z)={\frac {(-\log z)^{n-1}}{(n-1)!\;\;\;}},\;\;0<z\leq 1}$

The geometry of the product distribution of two random variables in the unit square.

The figure illustrates the nature of the integrals above. The shaded area within the unit square and below the line z = xy, represents the CDF of z. This divides into two parts. The first is for 0 < x < z where the increment of area in the vertical slot is just equal to dx. The second part lies below the xy line, has y-height z/x, and incremental area dx z/x.

The product of two independent Normal samples follows a modified Bessel function. Let ${\displaystyle x,y}$ be samples from a Normal(0,1) distribution and ${\displaystyle z=xy}$. Then

${\displaystyle p_{Z}(z)={\frac {K_{0}(|z|)}{\pi }},\;\;\;-\infty <z<+\infty }$

The variance of this distribution could be determined, in principle, by a definite integral from Gradsheyn and Ryzhik,[7]

${\displaystyle \int _{0}^{\infty }x^{\mu }K_{\nu }(ax)\,dx=2^{\mu -1}a^{-\mu -1}\Gamma {\Big (}{\frac {1+\mu +\nu }{2}}{\Big )}\Gamma {\Big (}{\frac {1+\mu -\nu }{2}}{\Big )},\;\;a>0,\;\nu +1\pm \mu >0}$

thus ${\displaystyle \int _{-\infty }^{\infty }{\frac {z^{2}K_{0}(|z|)}{\pi }}\,dz={\frac {4}{\pi }}\;\Gamma ^{2}{\Big (}{\frac {3}{2}}{\Big )}=1}$

A much simpler result, stated in a section above, is that the variance of the product of zero-mean independent samples is equal to the product of their variances. Since the variance of each Normal sample is one, the variance of the product is also one.

The product of correlated Normal samples case was recently addressed by Nadarajaha and Pogány.[8] Let ${\displaystyle X{\text{, }}Y}$ be zero mean, unit variance, normally distributed variates with correlation coefficient ${\displaystyle \rho {\text{ and let }}Z=XY}$

Then

${\displaystyle f_{Z}(z)={\frac {1}{\pi {\sqrt {1-\rho ^{2}}}}}\exp \left({\frac {\rho z}{1-\rho ^{2}}}\right)K_{0}\left({\frac {|z|}{1-\rho ^{2}}}\right)}$

Mean and variance: For the mean we have ${\displaystyle \operatorname {E} [Z]=\rho }$ from the definition of correlation coefficient. The variance can be found by transforming from two unit variance zero mean uncorrelated variables U, V. Let

${\displaystyle X=U,\;\;Y=\rho U+{\sqrt {(1-\rho ^{2})}}V}$

Then X, Y are unit variance variables with correlation coefficient ${\displaystyle \rho }$ and

${\displaystyle (XY)^{2}=U^{2}{\bigg (}\rho U+{\sqrt {(1-\rho ^{2})}}V{\bigg )}^{2}=U^{2}{\bigg (}\rho ^{2}U^{2}+2\rho {\sqrt {1-\rho ^{2}}}UV+(1-\rho ^{2})V^{2}{\bigg )}}$

Removing odd-power terms, whose expectations are obviously zero, we get

${\displaystyle \operatorname {E} [(XY)^{2}]=\rho ^{2}\operatorname {E} [U^{4}]+(1-\rho ^{2})\operatorname {E} [U^{2}]\operatorname {E} [V^{2}]=3\rho ^{2}+(1-\rho ^{2})=1+2\rho ^{2}}$

Since ${\displaystyle (\operatorname {E} [Z])^{2}=\rho ^{2}}$ we have

${\displaystyle \operatorname {Var} (Z)=\operatorname {E} [Z^{2}]-(\operatorname {E} [Z])^{2}=1+2\rho ^{2}-\rho ^{2}=1+\rho ^{2}}$

High correlation asymptote In the highly correlated case, ${\displaystyle \rho \rightarrow 1}$ the product converges on the square of one sample. In this case the ${\displaystyle K_{0}}$ asymptote is ${\displaystyle K_{0}(x)\rightarrow {\sqrt {\tfrac {\pi }{2x}}}e^{-x}{\text{ in the limit as }}x={\frac {|z|}{1-\rho ^{2}}}\rightarrow \infty }$ and

${\displaystyle {\begin{aligned}p(z)&\rightarrow {\frac {1}{\pi {\sqrt {1-\rho ^{2}}}}}\exp \left({\frac {\rho z}{1-\rho ^{2}}}\right){\sqrt {\frac {\pi (1-\rho ^{2})}{2z}}}\exp \left(-{\frac {|z|}{1-\rho ^{2}}}\right)\\&={\frac {1}{\sqrt {2\pi z}}}\exp {\Bigg (}{\frac {-|z|+\rho z}{(1-\rho )(1+\rho )}}{\Bigg )}\\&={\frac {1}{\sqrt {2\pi z}}}\exp {\Bigg (}{\frac {-z}{1+\rho }}{\Bigg )},\;\;z>0\\&\rightarrow {\frac {1}{\Gamma ({\tfrac {1}{2}}){\sqrt {2z}}}}e^{-{\tfrac {z}{2}}},\;\;{\text{ as }}\rho \rightarrow 1\\\end{aligned}}}$

which is a Chi-squared distribution with one degree of freedom.

Multiple correlated samples. Nadarajaha et. al. further show that if ${\displaystyle Z_{1},Z_{2},..Z_{n}{\text{ are }}n}$ iid random variables sampled from ${\displaystyle f_{Z}(z)}$ and ${\displaystyle {\bar {Z}}={\tfrac {1}{n}}\sum Z_{i}}$ is their mean then

${\displaystyle f_{\bar {Z}}(z)={\frac {n^{n/2}2^{-n/2}}{\Gamma ({\frac {n}{2}})}}|z|^{n/2-1}\exp \left({\frac {\beta -\gamma }{2}}z\right){W}_{0,{\frac {1-n}{2}}}(|z|),\;\;-\infty <z<\infty .}$

where W is the Whittaker function while ${\displaystyle \beta ={\frac {n}{1-\rho }},\;\;\gamma ={\frac {n}{1+\rho }}}$.

Using the identity ${\displaystyle W_{0,\nu }(x)={\sqrt {\frac {x}{\pi }}}K_{\nu }(x/2),\;\;x\geq 0}$, see for example the DLMF compilation. eqn(13.13.9),[9] this expression can be somewhat simplified to

${\displaystyle f_{\bar {z}}(z)={\frac {n^{n/2}2^{-n/2}}{\Gamma ({\frac {n}{2}})}}|z|^{n/2-1}\exp \left({\frac {\beta -\gamma }{2}}z\right){\sqrt {{\frac {\beta +\gamma }{\pi }}|z|}}\;K_{\frac {1-n}{2}}\left({\frac {\beta +\gamma }{2}}|z|\right),\;\;-\infty <z<\infty .}$

The pdf gives the distribution of a sample covariance.

Multiple non-central correlated samples. The distribution of the product of correlated non-central normal samples was derived by Cui et.al.[10] and takes the form of an infinite series of modified Bessel functions of the first kind.

Moments of product of correlated central normal samples

For a central normal distribution N(0,1) the moments are

${\displaystyle \operatorname {E} \left[X^{p}\right]={\frac {1}{\sigma {\sqrt {2\pi }}}}\int _{-\infty }^{\infty }x^{p}\exp(-{\tfrac {x^{2}}{2\sigma ^{2}}})\,dx={\begin{cases}0&{\text{if }}p{\text{ is odd,}}\\\sigma ^{p}(p-1)!!&{\text{if }}p{\text{ is even.}}\end{cases}}}$

where ${\displaystyle n!!}$ denotes the double factorial.

If ${\displaystyle X,Y\sim {\text{Norm}}(0,1)}$ are central correlated variables, the simplest bivariate case of the multivariate normal moment problem described by Kan,[11] then

${\displaystyle \operatorname {E} [X^{p}Y^{q}]={\begin{cases}0&{\text{if }}p+q{\text{ is odd,}}\\{\frac {p!q!}{2^{\tfrac {p+q}{2}}}}\sum _{k=0}^{t}{\frac {(2\rho )^{2k}}{{\Big (}{\frac {p}{2}}-k{\Big )}!\;{\Big (}{\frac {q}{2}}-k{\Big )}!\;(2k)!}}&{\text{if }}p{\text{ and }}q{\text{ are even}}\\{\frac {p!q!}{2^{\tfrac {p+q}{2}}}}\sum _{k=0}^{t}{\frac {(2\rho )^{2k+1}}{{\Big (}{\frac {p-1}{2}}-k{\Big )}!\;{\Big (}{\frac {q-1}{2}}-k{\Big )}!\;(2k+1)!}}&{\text{if }}p{\text{ and }}q{\text{ are odd}}\end{cases}}}$

where

${\displaystyle \rho }$ is the correlation coefficient and ${\displaystyle t=\min([p,q]/2)}$

[needs checking]

The distribution of the product of non-central correlated normal samples was derived by Cui et al.[10] and takes the form of an infinite series.

These product distributions are somewhat comparable to the Wishart distribution. The latter is the joint distribution of the four elements (actually only three independent elements) of a sample covariance matrix. If ${\displaystyle x_{t},y_{t}}$ are samples from a bivariate time series then the ${\displaystyle W=\sum _{t=1}^{K}{\dbinom {x_{t}}{y_{t}}}{\dbinom {x_{t}}{y_{t}}}^{T}}$ is a Wishart matrix with K degrees of freedom. The product distributions above are the unconditional distribution of the aggregate of K > 1 samples of ${\displaystyle W_{2,1}}$.

Let ${\displaystyle u_{1},v_{1},u_{2},v_{2}}$ be independent samples from a normal(0,1) distribution.
Setting ${\displaystyle z_{1}=u_{1}+iv_{1}{\text{ and }}z_{2}=u_{2}+iv_{2}{\text{ then }}z_{1},z_{2}}$ are independent zero-mean complex normal samples with circular symmetry. Their complex variances are ${\displaystyle \operatorname {Var} |z_{i}|=2.}$

The density functions of

${\displaystyle r_{i}\equiv |z_{i}|=(u_{i}^{2}+v_{i}^{2})^{\frac {1}{2}},\;\;i=1,2}$ are Rayleigh distributions defined as:

${\displaystyle f_{r}(r_{i})=r_{i}e^{-r_{i}^{2}/2}{\text{ of mean }}{\sqrt {\tfrac {\pi }{2}}}{\text{ and variance}}{\frac {4-\pi }{2}}}$

The variable ${\displaystyle y_{i}\equiv r_{i}^{2}}$ is clearly Chi-squared with two degrees of freedom and has PDF

${\displaystyle f_{y_{i}}(y_{i})={\tfrac {1}{2}}e^{-y_{i}/2}{\text{ of mean value }}2}$

Wells et. al.[12] show that the density function of ${\displaystyle s\equiv |z_{1}z_{2}|}$ is

${\displaystyle f_{s}(s)=sK_{0}(s),\;\;s\geq 0}$

and the cumulative distribution function of ${\displaystyle s}$ is

${\displaystyle P(a)=\Pr[s\leq a]=\int _{s=0}^{a}sK_{0}(s)ds=1-aK_{1}(a)}$

Thus the polar representation of the product of two uncorrelated complex Gaussian samples is

${\displaystyle f_{s,\theta }(s,\theta )=f_{s}(s)p_{\theta }(\theta ){\text{ where }}p(\theta ){\text{ is uniform on }}[0,2\pi ]}$.

The first and second moments of this distribution can be found from the integral in Normal Distributions above

${\displaystyle m_{1}=\int _{0}^{\infty }s^{2}K_{0}(s)\,dx=2\Gamma ^{2}({\tfrac {3}{2}})=2({\tfrac {\sqrt {\pi }}{2}})^{2}={\frac {\pi }{2}}}$

${\displaystyle m_{2}=\int _{0}^{\infty }s^{3}K_{0}(s)\,dx=2^{2}\Gamma ^{2}({\tfrac {4}{2}})=4}$

Thus its variance is ${\displaystyle \operatorname {Var} (s)=m_{2}-m_{1}^{2}=4-{\frac {\pi ^{2}}{4}}}$.

Further, the density of ${\displaystyle z\equiv s^{2}={|r_{1}r_{2}|}^{2}={|r_{1}|}^{2}{|r_{2}|}^{2}=y_{1}y_{2}}$ corresponds to the product of two independent Chi-square samples ${\displaystyle y_{i}}$ each with two DoF. Writing these as scaled Gamma distributions ${\displaystyle f_{y}(y_{i})={\tfrac {1}{\theta \Gamma (1)}}e^{-y_{i}/\theta }{\text{ with }}\theta =2}$ then, from the Gamma products below, the density of the product is

${\displaystyle f_{Z}(z)={\tfrac {1}{2}}K_{0}({\sqrt {z}}){\text{ with expectation }}\operatorname {E} (z)=4}$

The product of non-central independent complex Gaussians is described by O’Donoughue and Moura[13] and forms a double infinite series of modified Bessel functions of the first and second types.

The product of two independent Gamma samples, ${\displaystyle z=x_{1}x_{2}}$, defining ${\displaystyle \Gamma (x;k_{i},\theta _{i})={\frac {x^{k_{i}-1}e^{-x/\theta _{i}}}{\Gamma (k_{i})\theta _{i}^{k_{i}}}}}$, follows[14]

${\displaystyle {\begin{aligned}p_{Z}(z)&={\frac {2}{\Gamma (k_{1})\Gamma (k_{2})}}{\frac {z^{{\frac {k_{1}+k_{2}}{2}}-1}}{\left(\theta _{1}\theta _{2}\right)^{\frac {k_{1}+k_{2}}{2}}}}K_{k_{1}-k_{2}}\left(2{\sqrt {\frac {z}{\theta _{1}\theta _{2}}}}\right)\\\\&={\frac {2}{\Gamma (k_{1})\Gamma (k_{2})}}{\frac {y^{{\frac {k_{1}+k_{2}}{2}}-1}}{\theta _{1}\theta _{2}}}K_{k_{1}-k_{2}}\left(2{\sqrt {y}}\right){\text{ where }}y={\frac {z}{\theta _{1}\theta _{2}}}\\\end{aligned}}}$

Nagar et. al.[15] define a correlated bivariate beta distribution

${\displaystyle f(x,y)={\frac {x^{a-1}y^{b-1}(1-x)^{b+c-1}(1-y)^{a+c-1}}{B(a,b,c)(1-xy)^{a+b+c}}},\;\;\;0<x,y<1}$

where

${\displaystyle B(a,b,c)={\frac {\Gamma (a)\Gamma (b)\Gamma (c)}{\Gamma (a+b+c)}}}$

Then the pdf of Z = XY is given by

${\displaystyle f_{Z}(z)={\frac {B(a+c,b+c)z^{a-1}(1-z)^{c-1}}{B(a,b,c)}}{_{2}F_{1}}(a+c,a+c;a+b+2c;1-z),\;\;\;0<z<1}$

where ${\displaystyle {_{2}F_{1}}}$ is the Gauss hypergeometric function defined by the Euler integral

${\displaystyle {_{2}F_{1}}(a,b,c,z)={\frac {\Gamma (c)}{\Gamma (a)\Gamma (c-a)}}\int _{0}^{1}v^{a-1}(1-v)^{c-a-1}(1-vz)^{-b}\,dv}$