The distribution of the product of independent variance-gamma random variables

Let $X$ and $Y$ be independent variance-gamma random variables with zero location parameter; then the exact probability density function of the product $XY$ is derived. Some basic distributional properties are also derived, including formulas for the cumulative distribution function and the characteristic function, as well as asymptotic approximations for the density, tail probabilities and the quantile function. As special cases, we deduce some key distributional properties for the product of two independent asymmetric Laplace random variables as well as the product of four jointly correlated zero mean normal random variables with a particular block diagonal covariance matrix. As a by-product of our analysis, we deduce some new reduction formulas for the Meijer $G$-function.


Introduction
The variance-gamma (VG) distribution with parameters m > −1/2, 0 ≤ |β| < α, µ ∈ R, denoted by VG(m, α, β, µ), has probability density function (PDF) where the normalising constant is given by Here K m (x) is a modified Bessel function of the second kind, which is defined in Appendix A. Alternative parametrisations are given in [6,17,21].Other names include the generalized Laplace distribution [17, Section 4.1], Bessel function distribution [23] and the McKay Type II distribution [16].The VG distribution has an attractive distributional theory containing as special and limiting cases the normal, gamma and asymmetric Laplace distributions, as well as the product of correlated zero mean normal random variables amongst others [5].The VG distribution has been widely used in financial modelling [21,22] and further application areas and distributional properties are given in the review [5] and Chapter 4 of the book [17].In this paper, we set µ = 0. Let X ∼ VG(m, α 1 , β 1 , 0) and Y ∼ VG(n, α 2 , β 2 , 0), m, n > −1/2, 0 ≤ |β i | < α i , i = 1, 2, be independent VG random variables.In this paper, we derive an exact formula, expressed in terms of the Meijer G-function (defined in Appendix A), for the PDF of the product XY .Our formula generalises one of [13] for the PDF of the product XY in the case β 1 = β 2 = 0 (a product of symmetric VG random variables).We derive our formula for the PDF using the theory of Mellin transforms.A systematic account of the application of Mellin transforms to derive formulas for the PDF of products of independent random variables was first given by [4], and this approach has since been used to find exact formulas for the PDFs of products of a number of important continuous random variables; see, for example, [18,19,29,30,31,34].
With our formula for the PDF of the product XY at hand, we are able to derive a number of other key distributional properties, including formulas for the cumulative distribution function (CDF) and the characteristic function, as well as asymptotic approximations for the PDF, tail probabilities and the quantile function, none of which were given in [13].Our asymptotic analysis of the PDF of the product XY at the origin implies that the distribution is unimodal with mode at the origin.All these results are stated and proved in Section 2.
In the case m = 1/2, µ = 0 the VG distribution corresponds to the asymmetric Laplace distribution (with zero location parameter), and further setting β = 0 yields the classical symmetric Laplace distribution.A comprehensive account of the distributional theory of the Laplace and asymmetric Laplace distributions is given in the excellent book [17]; however, there appears to be a surprising gap in the literature regarding the distributional theory of the product of independent asymmetric Laplace distributions, with, to our best knowledge, only formulas for the PDF and CDF of the symmetric Laplace distribution (with zero location parameter) available in [25].We fill in this gap in Section 3.1 by providing formulas for the PDF, CDF and characteristic function of the product of two independent asymmetric Laplace random variables (with zero location parameter).
Also, in the case m = µ = 0, the VG distribution corresponds to the distribution of the product of two correlated zero mean normal random variables, which itself has numerous applications that date back to the 1930's with the works of [2,35]; see [10] for an overview.Our formula for the PDF corresponding to this case is given in Corollary 3.6, which yields an exact formula for the PDF of the product N 1 N 2 N 3 N 4 of four jointly correlated zero mean normal random variables (N 1 , N 2 , N 3 , N 4 ) ⊺ for which (N 1 , N 2 ) ⊺ and (N 3 , N 4 ) ⊺ are independent.To the best of our knowledge, this is one of the first explicit formulas in the literature for the PDF of the product of three or more zero mean correlated normal random variables; see Remark 3.8 for a further discussion.
The formulas obtained in this paper are expressed in terms of a number of standard special functions, all of which are defined in Appendix A. As a by-product of our analysis, we obtain some new reduction formulas for the Meijer G-function, which are collected in Section 4.

The variance-gamma product distribution 2.1 Probability density function
In the following theorem, we provide an explicit formula for the PDF of the product Z = XY as an infinite series involving the Meijer G-function.The PDF of the product Z is plotted for a range of parameter values in Figures 1 and 2. Figure 1 shows the effect of varying the shape parameters m and n, whilst Figure 2 demonstrates the effect of varying the skewness parameters 0) are independent, and let Z = XY denote their product.Then, for z ∈ R, which is in agreement with equation (37) of [13].If one of the skewness parameters is equal to zero, then the PDF (2.2) reduces to a single infinite series.Without loss of generality, let We observe that the PDF f Z (z) is symmetric about the origin if β 1 = 0 or β 2 = 0. We remark that the exact PDF of the ratio X/Y , where X ∼ VG(m, α 1 , β 1 , 0) and Y ∼ VG(n, α 2 , β 2 , 0) are independent, is expressed as a double sum involving the Gaussian hypergeometric function when both β 1 and β 2 are non-zero [12], and the PDF simplifies to a single infinite series of hypergeometric functions if one of the skewness parameters is zero [12], and to just a single hypergeometric function if β 1 = β 2 = 0 [24].The increase in complexity of the PDF of the product XY for non-zero skewness parameters thus mirrors that of the ratio X/Y .A similar behaviour also occurs for the product of two correlated normal random variables as one transitions from zero to non-zero means [3].
Proof of Theorem 2.1.We first consider the case z > 0. For a random variable V , we define its positive part V + by V + = max{V, 0} and its negative part V − by V − = max{−V, 0}.Utilising the Taylor expansion of the exponential function and evaluating the resulting integrals with formula (A.58), we derive the Mellin transform for X + : Similarly, we obtain the Mellin transform for X − : The Mellin transform of Z + is then given by Taking the inverse Mellin transform yields the PDF of Z + , and we need to find: where c > 0 and we evaluated the contour integral using a change of variables and the contour integral definition (A.69) of the Meijer G-function.On combining (2.4), (2.5), (2.6) and (2.7) we deduce a formula for the PDF of Z for z > 0. Arguing similarly we can obtain a formula for the PDF of Z for z < 0, and overall we get that, for z ∈ R, where a ij = (1 + (−1) i+j )/2, i, j ≥ 0, so that a ij = 1 if i + j is even, and a ij = 0 if i + j is odd.The formula (2.8) can be easily rewritten as (2.2), and this completes the proof. 2 Remark 2.3.As noted by one of the reviewers, the series representations (2.4) and (2.5) of M f X + (s) and M f X − (s), respectively, can also be represented in terms of the Fox-Wright function (defined in Appendix A): Using special properties of the modified Bessel function of the second kind, simpler formulas for the PDF of Z = XY can be obtained.We illustrate this in the following Proposition 2.4.Using reduction formulas for the Meijer G-function it is also possible to obtain simpler formulas for the PDF; see Section 3 in which we provide simplified formulas for certain special cases of the VG product distribution.
Proposition 2.4.Suppose that m − 1/2 ≥ 0 and n − 1/2 ≥ 0 are non-negative integers.Then, (2.2) can be simplified to a finite double sum expressed in terms of the modified Bessel function of the second kind: Proof.The modified Bessel function of the second kind takes on an elementary form (A.57) in this case.Consequently, the PDF of Z = XY can be expressed as Evaluating the integrals using (A.56) yields (2.9).

Cumulative distribution function
Let F Z (z) = P(Z ≤ z) denote the CDF of the product Z = XY .In the following theorem, we provide an exact formula for the CDF of Z for the case of independent symmetric VG random variables ( Proof.For ease of exposition, we set α 1 = α 2 = 1, with general case following from a simple rescaling.As the PDF (as given by (2.3) in this case) is symmetric about the origin when β 1 = β 2 = 0, we also only consider the case z > 0. Letting 1/C m,n = 4πΓ(m + 1/2)Γ(n + 1/2), we have, for z > 0, where in the third step we used the indefinite integral formula (A.77) and in the final step we used that the Meijer G-function in (2.12) evaluated at u = 0 is equal to zero, which is easily seen from the contour integral representation (A.69) of the G-function.Finally, we deduce (2.11) from (2.10) using the relation (A.71).
When m − 1/2 ≥ 0 and n − 1/2 ≥ 0 are non-negative integers, a finite double sum representation of the CDF can be given, even for non-zero skewness parameters.For µ where tµ,ν (x) is a normalisation of the modified Lommel function of the first kind t µ,ν (x), which is defined in Appendix A.
A simple formula can be given for the probability P(Z ≤ 0).We will make use of a recent result of [11], which states that, for X ∼ VG(m, α, β, 0), where and the Gaussian hypergeometric function is defined in Appendix A.
Proof.A simple calculation that exploits the fact that X and Y are independent shows that , and the result now follows from the formula (2.16).

Characteristic function
Let φ Z (t) = E[e itZ ] denote the characteristic function of the product Z = XY .In the following theorem, we provide an exact formula for the characteristic function for the case of independent symmetric VG random variables (β 1 = β 2 = 0).We note that the moment generating function M Z (t) = E[e tZ ] is not defined for t ̸ = 0, which is easily seen from the limiting forms (2.22) and (2.23).
Theorem 2.9.Let (2.18) Proof.Since the PDF (2.3) of Z is symmetric about the origin we have that where we evaluated the integral using formula 5.6.3(18) of [20] and obtained the final equality using (A.70).Finally, we obtain (2.18) from (2.17) using the relation (A.71).
When m − 1/2 ≥ 0 and n − 1/2 ≥ 0 are non-negative integers, a finite double sum representation of the characteristic function can be given, even for non-zero skewness parameters.The formula is expressed in terms of the Whittaker function, which is defined in Appendix A.
Proposition 2.10.Suppose m − 1/2 ≥ 0 and n − 1/2 ≥ 0 are non-negative integers and that z) dz using the formula (2.9) for the PDF and then compute the resulting integral using the integral formula (A.80).

Asymptotic behaviour of the density, tail probabilities and quantile function
The representations of the PDF and CDF from Sections 2.1 and 2.2 are rather complicated and difficult to parse on first inspection.Some insight can be gained from the following asymptotic approximations.
Remark 2.14.The limiting forms for the CDF and PDF of Z given in Proposition 2.11 and Corollary 2.12 simplify for certain parameter values.For example, if as z → ∞, and therefore the limiting form (2.20) simplifies to Remark 2.15.The VG(m, α, β, 0) distribution is also unimodal; however, the mode is at the origin if and only if m ≤ 1/2 (see [5, Section 2.8]).We see that the growth of the singularity of the PDF f Z (z) in the limit z → 0 increases as m and n decrease below 0, which is to be expected given that VG(m, α, β, 0) PDF is bounded for m > 0, has a logarithmic singularity at the origin for m = 0, and a power law singularity at the origin for −1/2 < m < 0 [5, Section 2.1].
For 0 < p < 1, let Q(p) = F −1 (p) denote the quantile function of the product Z.In the following proposition, we give asymptotic approximations for the quantile function.
We now prove Proposition 2.11.We will make use of Lemma 2.1 of [1]: Lemma 2.17 (Arendarczyk and Dȩbicki [1]).Let X 1 and X 2 be independent, non-negative random variables such that, as x → ∞, Proof of Proposition 2.11.We will denote the positive and negative parts of a random variable X by X + and X − , respectively.Let X ∼ VG(m, α 1 , β 1 , 0) and Y ∼ VG(n, α 2 , β 2 , 0) be independent.The asymptotic behaviour of the CDFs of X + , X − and Y + , Y − are given in [5] as follows: as x, y → ∞, Applying Lemma 2.17 now gives that, as z → ∞, The limiting form (2.20) now follows from the formula FXY (z) = FX + Y + (z) + FX − Y − (z).The limiting form (2.21) is obtained similarly, since Proof of Corollary 2.12.Since f Z (z) = − F ′ Z (z) and f Z (z) = F ′ Z (z), we deduce the limiting forms (2.22) and (2.23) from differentiating the limiting forms (2.20) and (2.21), respectively, and retaining the highest order terms. 2 Proof of Proposition 2.13.Recall that the PDF (2.2) is a double infinite series of Meijer Gfunctions.To ease notation, we will let x = α 1 α 2 z/4.To determine the behaviour of f Z (z) as z → 0 we study the asymptotic behaviour as x ↓ 0 of the functions where i, j ≥ 0 are non-negative integers such that j ≤ 2i, and it suffices to study the limit x ↓ 0 since g i,j,m,n (x) is an even function in x.Using the contour integral definition of the Meijer G-function (A.69) followed by an application of the residue theorem we can write where the integration path L is a loop that encircles the poles of the gamma functions and A is the collection of these poles.We now compute the residues in (2.30) for the five different cases of m and n as given in the statement of the proposition.
1. Let m, n > 0. Suppose first that i = j = 0. We need to compute for all poles s * of the gamma functions Γ(s), Γ(m + s) and Γ(n + s).First, we consider the pole at s = 0, for which is suffices to find the coefficient of s −1 in the Laurent expansion of (Γ(s)) 2 Γ(m + s)Γ(n + s)x −2s . (2.32) We will make use of the Taylor expansion of x −2s about s = 0, as given by as well as the standard asymptotic expansion where γ is the Euler-Mascheroni constant.Using (2.33) and (2.34) we find after a simple calculation that the coefficient of s −1 in the Laurent expansion of (2.32) is given by Γ(m)Γ(n)[−2 ln(x)+ 2γ], and therefore It is easily seen that for all other poles besides the one at s = 0 the residue (2.31) is of a smaller asymptotic order in the limit x ↓ 0 than the residue (2.35).Similarly, it is readily seen that for i ≥ 1 or j ≥ 1 the residues in (2.30) are of a smaller asymptotic order in the limit x ↓ 0 than the residue (2.35).Finally, we conclude that, as x → 0, (where we can state the limiting form for x → ∞ as it sufficed to consider the limit x ↓ 0) and that g i,j,m,n (x) is of smaller asymptotic order in the limit x → 0 (equivalently in the limit z → 0) for i ≥ 1 or j ≥ 1.Thus, as z → 0, 2. Let m = 0, n > 0. Arguing similarly to in part 1, we find that the dominant residue in (2.30) in the limit x ↓ 0 is given by with all other residues being of lower asymptotic order in the limit x ↓ 0. Arguing as in part 1, we deduce the limiting form (2.24).
4. Let m < 0 and m < n.Since m < 0, the dominant residue now is located at the pole s = −m, and is given by and we thus obtain the limiting form (2.26).
5. Let m = n < 0. As m = n < 0, the dominant pole is again located at the pole s = −m, and is this time given by We thus obtain the limiting form (2.27).
The final assertion that the distribution of Z is unimodal with mode 0 for all parameter constellations follows because parts 1-5 imply that the PDF f Z (z) has a singularity at the origin for all parameter values, and the PDF is bounded everywhere except for the singularity at z = 0. 2 To prove Proposition 2.16 we will make use of the following lemma.
2 ) and Y ∼ Laplace(α 2 ) are mutually independent random variables.Denote their product by where the Struve function K 0 is defined as in (A.61).
We now consider the more general setting that X ∼ VG(m, α 1 , 0, 0) and Y ∼ Laplace(α 2 ) are independent random variables.In this case, a simple formula can still be given for the characteristic function of the product Z = XY .The formula (3.48) is valid provided m + 1/2 ∈ R + \ Z; if this condition is not met, then one can apply Proposition 2.10 instead.

Product of correlated zero mean normal random variables
Let (U i , V i ) ⊺ , i = 1, 2, be independent bivariate normal random vectors with zero mean vector, variances (σ 2 U i , σ 2 V i ) ⊺ and correlation coefficient ρ i .It was noted by [6,8] that the product which yielded a new derivation of the PDF of the product W i which was obtained independently by [15,26,35].In the following corollary, we provide formulas for the PDF, CDF and characteristic function of the product Z = W 1 W 2 .
Corollary 3.6.Let the previous notations prevail.Also, let s (iii) Suppose that ρ 1 = ρ 2 = 0.Then, for t ∈ R, ) ⊺ is a multivariate normal random vector with zero mean vector and covariance matrix for which s = σ 1 σ 2 σ 3 σ 4 .In the case ρ 1 = ρ 2 = 0, we have that Z is thus distributed as the product of four independent zero mean normal random variables, and setting ρ 1 = ρ 2 = 0 in the PDF (3.50) gives the formula which is a special case of a more general formula of [31] for the product of k ≥ 2 independent zero mean normal random variables.Whilst formulas for the PDF for the product of two correlated zero mean normal random variables and the product of k independent zero mean normal random variables are known in the literature, the case of the product of three or more correlated zero mean normal random variables is not well-understood.To the best of our knowledge, the formula (3.50) is one of the first such explicit formulas, and we hope that the result inspires other researchers to find further explicit formulas for the PDF of the product of three of more correlated zero mean normal random variables.
2. We note that formulas, expressed in terms of the Meijer G-function, for the CDF of the product of k independent standard normal random variables have previously been given by [32].The results can be summarised as follows.Let X 1 , . . ., X k be independent normal random variables with X i ∼ N (0, σ 2 i ), i = 1, . . ., k. Denote the characteristic function of the product In the cases k = 2, 3, 4, this formula simplifies to We could not find reduction formulas in the literature to obtain further simplifications for k ≥ 5.

Reduction formulas for the Meijer G-function
As a by-product of our analysis, we obtain the following reduction formulas for the Meijer G-function, which we believe to be new.where the function G µ,ν (x) is defined as in (2.13).In the case a = b = 0, the formula (4.54) simplifies to

A Special functions
In this appendix, we define several special functions, and state some of their basic properties that are needed in this paper.Unless otherwise stated, these and further properties can be found in the standard references [14,20,27].

A.1 Bessel, Struve and Lommel functions
The modified Bessel function of the second kind can be defined, for ν ∈ R and x > 0, by

.45) Remark 3 . 3 .
Formulas (3.42) and(3.43)for the PDF and CDF of the product of two independent symmetric Laplace random variables are in agreement with the formulas of[25] for the PDF and CDF of the product |XY |.However, to the best of our knowledge, the other formulas of Corollaries 3.1 and 3.2 are new.Proof of Corollary 3.1.(i): Set m = n = 1/2 in the formula of Proposition 2.4.

Remark 3 . 7 .
.51) Proof.(i) Combine (2.2) with (3.49).(ii) Combine (2.10) with (3.49) (with ρ 1 = ρ 2 = 0).(iii) We obtain (3.51) by applying the reduction formula (A.72) (with a = b = 0) to (2.17).The formula (3.51) can be connected to the modified Bessel function K 0 by Nicholson's formula (A.59): 3. A formula, in terms of the Meijer G-function, for the characteristic function of Z = N 1 N 2 N 3 N 4 (with ρ 1 = ρ 2 = 0) has previously been given by [7, Propostion 2.2], as a special case of a general formula for the characteristic function of the product of k independent zero mean normal random variables.Formula (3.51) provides a simpler expression for the characteristic function for the case of the product of k = 4 independent zero mean normal random variables and complements other simpler formulas for the cases k = 2[33, Example 11.22] and k = 3 [7, Propostion 2.2].