How to calculate the PDF of product of two random variables from generalized gamma distributions?

Question

Namely, we want to find the explicit formula of PDF of double Generalized Gamma distribution.

MellinConvolve[
 PDF[GammaDistribution[φ1,(Ω1/φ1)^(1/k1),k1,0],x],
 PDF[GammaDistribution[φ2,(Ω2/φ2)^(1/k2),k2,0],x],
x,y]

TransformedDistribution[
  x*y
 , {
     x \[Distributed]  GammaDistribution[φ1,(Ω1/φ1)^(1/k1), k1, 0], 
     y \[Distributed]  GammaDistribution[φ2,(Ω2/φ2)^(1/k2), k2, 0]
   }
] // PDF[#, t] &

pdfx = (E^(-((x^k1 φ1)/Ω1))  k1 x^(-1 +  k1 φ1) (φ1/Ω1)^φ1) /Gamma[φ1]

pdfy = (E^(-((y^k2  φ2)/Ω2)) k2  y^(-1 + k2  φ2)  (φ2/Ω2)^φ2)/Gamma[φ2] /. {y -> z/x};

Assuming[
   {k1 > 0, φ1 > 0, Ω1 > 0}, 
   Integrate[
      pdfx*pdfy*1/x
      , {x, 0, Infinity}
      , PrincipalValue -> True
   ]
]

None of them can produce explicit result.

---

The expected result is (I have VERIFIED some test cases under some parameter settings, it's true):

$$ f_T(t)=\frac{k_2 \alpha \beta^{\varphi_1-1 / 2} \alpha^{\varphi_2-1 / 2}(2 \pi)^{1-(\alpha+\beta) / 2}}{t \Gamma\left(\varphi_1\right) \Gamma\left(\varphi_2\right)} \times G_{q, 0}^{0, q}\left[\left.\frac{\alpha^\alpha \beta^\beta \Omega_1^\beta \Omega_2^\alpha}{t^{k_2 \alpha} \varphi_1^\beta \varphi_2^\alpha} \right\rvert\, \begin{array}{c} \Delta(\beta,1-\varphi_1),\Delta(\alpha,1-\varphi_2) & - \\ - & - \end{array}\right] $$

Where $q=\alpha+\beta$, $\alpha,\beta \in \mathbb{N}^+$ and $\alpha/\beta=k_1/k_2$; $\Delta(m, n) \triangleq \frac{n}{m}, \frac{n+1}{m}, \cdots, \frac{n+m-1}{m}$

On[Assert];
Assert[ResourceFunction["RationalQ"][Rationalize[k1/k2]]]
MyDeltaSign[m_Integer /; (m > 0), n_] := 
 Table[(n + j)/m, {j, 0, m - 1}]
{\[Alpha], \[Beta]} = Rationalize[k1/k2] // NumeratorDenominator
k2*\[Alpha]*\[Beta]^(\[CurlyPhi]1 - 1/2)*\[Alpha]^(\[CurlyPhi]2 - 
      1/2)*(2*Pi)^(1 - (\[Alpha] + \[Beta])/2)/(t*Gamma[\[CurlyPhi]1]*
      Gamma[\[CurlyPhi]2])*
   MeijerG[{MyDeltaSign[\[Beta], 1 - \[CurlyPhi]1]~Join~
      MyDeltaSign[\[Alpha], 
       1 - \[CurlyPhi]2], {}}, {{}, {}}, (\[Alpha]^\[Alpha]*\[Beta]^\
\[Beta]*\[CapitalOmega]1^\[Beta]*\[CapitalOmega]2^\[Alpha])/(t^(k2*\
\[Alpha])*\[CurlyPhi]1^\[Beta]*\[CurlyPhi]2^\[Alpha])] 

Rewrite the above integral expression in FoxH: (firsly replace /. {x :> \[Tau]^(1/k1)}, then using formula (07.34.21.0012.01) from MeijerG.pdf)

k2 (\[CurlyPhi]1/\[CapitalOmega]1)^(1/k1) (\[CurlyPhi]2/\[CapitalOmega]2)^(1/k2) Gamma[\[CurlyPhi]1] Gamma[\[CurlyPhi]2]*FoxH[
(*Upper List*){
    (*Upper Front list*){{1-(-(1/k1)+\[CurlyPhi]1), -k2/k1}},
    (*Upper Rear List*){}
  },
(*Lower List*){
    (*Lower Front List*){{-(1/k2)+\[CurlyPhi]2, 1}},
    (*Lower Rear List*){}
  }, 
(t^k2 \[CurlyPhi]2 (\[CurlyPhi]1/\[CapitalOmega]1)^(k2/k1))/\[CapitalOmega]2]

Answer 1

Even for concrete values of parameters the result is in terms of MeijerG:

\[CurlyPhi]1 = 1; \[CurlyPhi]2 = 2; k1 = 1; k2 = 2; \[CapitalOmega]1 
= 3; \[CapitalOmega]2 = 4; pr = TransformedDistribution[x*y, 
{x \[Distributed]GammaDistribution[\[CurlyPhi]1,(\[CapitalOmega]1/\[CurlyPhi]1)^(1/k1), k1, 0], 
y \[Distributed]GammaDistribution[\[CurlyPhi]2,(\[CapitalOmega]2/\[CurlyPhi]2)^(1/k2), k2, 0]}] ;
PDF[pr, t]

Piecewise[ {{(t^3*MeijerG[{{}, {}}, {{-3/2, -1, 0}, {}}, t^2/72])/ (2592*Sqrt[Pi]), t > 0}}, 0]

Answer 2

If we help Mathematica a little, we will achieve the goal and obtain the solution to this integral.

I use Mathematica version 13.3, because version 14.0 have some Bugs in FunctionExpand to calculate a numeric values FoxH function see here and here

You want calculate this integral:

$$\int_0^{\infty } \frac{e^{-\frac{x^{\text{k1}} \text{$\varphi $1}}{\text{$\Omega $1}}-\frac{\left(\frac{t}{x}\right)^{\text{k2}} \text{$\varphi $2}}{\text{$\Omega $2}}} \text{k1} \text{k2} x^{-2+\text{k1} \text{$\varphi $1}} \left(\frac{t}{x}\right)^{-1+\text{k2} \text{$\varphi $2}} \left(\frac{\text{$\varphi $1}}{\text{$\Omega $1}}\right)^{\text{$\varphi $1}} \left(\frac{\text{$\varphi $2}}{\text{$\Omega $2}}\right)^{\text{$\varphi $2}} \Gamma (\text{$\varphi $1})}{\Gamma (\text{$\varphi $2})} \, dx$$

Using Mellin Transform and Inverse Mellin Transform then:

$Version
(*"13.3.0 for Microsoft Windows (64-bit) (June 3, 2023)"*)

F1 = MellinTransform[(E^(-A*(x^k1 \[CurlyPhi]1)/\[CapitalOmega]1 - (
 t^k2 x^-k2 \[CurlyPhi]2)/\[CapitalOmega]2)
 k1 k2 t^(-1 + k2 \[CurlyPhi]2)
 x^(-1 + k1 \[CurlyPhi]1 - 
 k2 \[CurlyPhi]2) \[CurlyPhi]1^\[CurlyPhi]1 \[CurlyPhi]2^\
 \[CurlyPhi]2 \[CapitalOmega]1^-\[CurlyPhi]1 \[CapitalOmega]2^-\
 \[CurlyPhi]2 Gamma[\[CurlyPhi]1])/Gamma[\[CurlyPhi]2], A, s] // 
 PowerExpand(*Where: A=1*)

 (*(E^(-((t^k2 x^-k2 \[CurlyPhi]2)/\[CapitalOmega]2)) k1 k2 t^(-1 + 
 k2 \[CurlyPhi]2) x^(-1 - k1 s + k1 \[CurlyPhi]1 - 
 k2 \[CurlyPhi]2) \[CurlyPhi]1^(-s + \[CurlyPhi]1) \[CurlyPhi]2^\
 \[CurlyPhi]2 \[CapitalOmega]1^(
 s - \[CurlyPhi]1) \[CapitalOmega]2^-\[CurlyPhi]2 Gamma[s] Gamma[\[CurlyPhi]1])/Gamma[\[CurlyPhi]2]*)

 
 F2 = Integrate[F1, {x, 0, \[Infinity]}, Assumptions -> {\[CapitalOmega]2 > 0, \[CapitalOmega]1 > 0, k1 > 0, 
 k2 > 0, \[CurlyPhi]1 > 0, \[CurlyPhi]2 > 0, t > 0}]

 (*ConditionalExpression[(k1 t^(-1 + 
 k2 \[CurlyPhi]2) (\[CurlyPhi]1/\[CapitalOmega]1)^(-s + \
 \[CurlyPhi]1) (\[CurlyPhi]2/\[CapitalOmega]2)^\[CurlyPhi]2 ((
 t^k2 \[CurlyPhi]2)/\[CapitalOmega]2)^(-((
 k1 s - k1 \[CurlyPhi]1 + k2 \[CurlyPhi]2)/k2))
 Gamma[s] Gamma[\[CurlyPhi]1] Gamma[(k1 (s - \[CurlyPhi]1))/
 k2 + \[CurlyPhi]2])/Gamma[\[CurlyPhi]2], 
 k1 \[CurlyPhi]1 < k2 \[CurlyPhi]2 + k1 Re[s]]*)

 InverseMellinTransform[F2[[1]], s, A] /. A -> 1(*Can't compute.Probably in newer version like 14.1 can computed*)

Using another method to compute Inverse Mellin Transform:

 Sum[Assuming[m \[Element] Integers && m >= 0, 
 Residue[F2[[1]]*A^(-s) /. A -> 1, {s, -m}]], {m, 0, 
 Infinity}](*Can't compute.Probably in newer version like 14.1 I hope can computed*)

$\sum _{m=0}^{\infty } \frac{(-1)^m \text{k1} t^{-1+\text{k2} \text{$\varphi $2}} \left(\frac{\text{$\varphi $1}}{\text{$\Omega $1}}\right)^{m+\text{$\varphi $1}} \left(\frac{\text{$\varphi $2}}{\text{$\Omega $2}}\right)^{\text{$\varphi $2}} \left(\frac{t^{\text{k2}} \text{$\varphi $2}}{\text{$\Omega $2}}\right)^{\frac{\text{k1} m}{\text{k2}}+\frac{\text{k1} \text{$\varphi $1}}{\text{k2}}-\text{$\varphi $2}} \Gamma (\text{$\varphi $1}) \Gamma \left(-\frac{\text{k1} (m+\text{$\varphi $1})}{\text{k2}}+\text{$\varphi $2}\right)}{m! \Gamma (\text{$\varphi $2})}=\frac{\text{k1} t^{-1+\text{k2} \text{$\varphi $2}} \left(\frac{\text{$\varphi $1}}{\text{$\Omega $1}}\right)^{\text{$\varphi $1}} \left(\frac{\text{$\varphi $2}}{\text{$\Omega $2}}\right)^{\text{$\varphi $2}} \left(\frac{t^{\text{k2}} \text{$\varphi $2}}{\text{$\Omega $2}}\right)^{\frac{\text{k1} \text{$\varphi $1}}{\text{k2}}-\text{$\varphi $2}} \sum _{m=0}^{\infty } \frac{\left(-\frac{\text{$\varphi $1} \left(\frac{t^{\text{k2}} \text{$\varphi $2}}{\text{$\Omega $2}}\right)^{\text{k1}/\text{k2}}}{\text{$\Omega $1}}\right)^m \Gamma \left(-\frac{\text{k1} (m+\text{$\varphi $1})}{\text{k2}}+\text{$\varphi $2}\right)}{m!}}{\Gamma (\text{$\varphi $2})}$

Using:

 Sum[ Gamma[a1 + A1 m]*z^m/m!, {k, 0, Infinity}] == FoxH[{{{1 - a1, A1}}, {}}, {{{0, 1}}, {}}, -z]
 (*And again, it's impossible to compute, it's a shame.
 I have the impression that Mathematica can't do anything.*)

Then we have the Solution expressed by FoxH function:

 F3=(k1*t^(-1 + k2*\[CurlyPhi]2)*(\[CurlyPhi]1/\[CapitalOmega]1)^\[CurlyPhi]1*(\[CurlyPhi]2/\[CapitalOmega]2)^\[CurlyPhi]2*((t^k2*\[CurlyPhi]2)/\[CapitalOmega]2)^((k1*\[CurlyPhi]1)/k2 - \[CurlyPhi]2)*FoxH[{{{1 + (k1*\[CurlyPhi]1)/k2 - \[CurlyPhi]2, -(k1/k2)}}, {}}, {{{0, 1}}, {}}, (\[CurlyPhi]1*((t^k2*\[CurlyPhi]2)/\[CapitalOmega]2)^(k1/k2))/\[CapitalOmega]1]*Gamma[\[CurlyPhi]1])/Gamma[\[CurlyPhi]2]

 \[CurlyPhi]1 = 1; \[CurlyPhi]2 = 2; k1 = 1; k2 = 2; 
 \[CapitalOmega]1 = 3; \[CapitalOmega]2 = 4; t = 1;

  F3 // FunctionExpand // N
  (*0.162355*)

  pr = TransformedDistribution[x*y, {x \[Distributed] 
  GammaDistribution[\[CurlyPhi]1, \
  (\[CapitalOmega]1/\[CurlyPhi]1)^(1/k1), k1, 0], 
  y \[Distributed] 
  GammaDistribution[\[CurlyPhi]2, \
  (\[CapitalOmega]2/\[CurlyPhi]2)^(1/k2), k2, 0]}];
  PDF[pr, 1.0]
   (*0.162355*)

Mathematica 13.3 or 14.0 can't convert FoxH to MeijerG function!