# Comptes Rendus
Mathématique

Robert E. Gaunt and Siqi Li

The variance-gamma ratio distribution

Volume 361 (2023), p. 1151-1161

Published online: 24 October 2023

https://doi.org/10.5802/crmath.495

<figure><img src="image_1.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
<figure><img src="image_2.png" /></figure>

# The variance-gamma ratio distribution

Robert E. Gaunt*, a and Siqi Li a

**Abstract.** Let X and Y be independent variance-gamma random variables with zero location parameter; then the exact probability density function of the ratio X/Y is derived. Some basic distributional properties are also derived, including identification of parameter regimes under which the density is bounded, asymptotic approximations of tail probabilities, and fractional moments; in particular, we see that the mean is undefined. In the case that X and Y are independent symmetric variance-gamma random variables, an exact formula is also given for the cumulative distribution function of the ratio X/Y.

**Keywords.** Variance-gamma distribution, ratio distribution, product of correlated normal random variables, hypergeometric function, Meijer G-function.

**2020 Mathematics Subject Classification.** 60E05, 62E15.

**Funding.** SL is supported by a University of Manchester Research Scholar Award.

Manuscript received 24 February 2023, accepted 10 March 2023.

## 1. Introduction

The variance-gamma (VG) distribution with parameters $m > -1/2$, $0 \le |\beta| < \alpha$, $\mu \in \mathbb{R}$ has probability density function (PDF)

$$f_X(x) = Me^{\beta(x-\mu)} |x-\mu|^m K_m(\alpha|x-\mu|), \quad x \in \mathbb{R}, \qquad (1)$$

where the normalising constant is given by

$$M = M_{m,\alpha,\beta} = \frac{\gamma^{2m+1}}{\sqrt{\pi}(2\alpha)^m \Gamma(m+1/2)},$$

with $\gamma^2 = \alpha^2 - \beta^2$. Here $K_m(x)$ is a modified Bessel function of the second kind, which is defined in Appendix A. If the random variable X has PDF (1), then we write $X \sim VG(m, \alpha, \beta, \mu)$. The VG distribution is also referred to as the Bessel function distribution [17], the generalized Laplace distribution [13] and the McKay Type II distribution [11]. Alternative parametrisations can be found in [4, 13, 15]. The VG distribution is widely used in financial modelling [15, 16], and further application areas and distributional properties are given in Chapter 4 of the book [13]. In this paper, we set $\mu = 0$.

Let X and Y be independent random variables. The distribution of the ratio X/Y for X ~ VG(m, 0, 1, 0) and Y ~ VG(n, 0, 1, 0), m, n > 1, was studied by [19] (they referred to X and Y

<footnote>a Department of Mathematics, The University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote>E-mails: robert.gaunt@manchester.ac.uk, siqi.li-8@postgrad.manchester.ac.uk</footnote>

<footnote>* Corresponding author.</footnote>
as Bessel function random variables). In particular, an exact formula, expressed in terms of the Gaussian hypergeometric function (defined in Appendix A), was obtained for the PDF. The study of [19] was motivated by the fact that the ratio of independent random variables arises throughout the sciences and the already mentioned fact that VG random variables appear in many applications, so that the ratio X/Y could, for example, represent the ratio of distributions of log-returns of two different financial assets. However, it has been observed, for example, that when fitting the VG distribution to log-returns of financial assets the skewness parameter $\beta$ is non-zero (see, for example, [22]). This provides motivation for this paper, in which we derive the exact distribution of the VG ratio X/Y for independent $X \sim \text{VG}(m, \alpha_1, \beta_1, 0)$ and $Y \sim \text{VG}(n, \alpha_2, \beta_2, 0)$, with $m, n > -1/2$, $0 \le |\beta_i| < \alpha_i$, $i = 1, 2$, thereby generalising results of [19]. In the case of a ratio of independent symmetric VG random variables ($\beta_1 = \beta_2 = 0$), we also derive an exact formula for the cumulative distribution function (CDF) of the ratio X/Y, a key distributional property that is not given in [19]. This formula is expressed in terms of the Meijer G-function, which is defined in Appendix A.

We stress that our results hold for $m, n > -1/2$, a wider range of validity than claimed by [19]. This is significant because the regime $m < 1$ is often encountered in applications; for example, when fitting the VG distribution to log returns of financial assets [22]. Another source of interest is that 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 dating back to 1936 with the work of [1]; see [6] for an overview. Results corresponding to this case are given in Corollary 13. Moreover, the case $m, n < 1$ warrants attention because the VG ratio distribution undergoes a significant change in behaviour in that the density is bounded for $m > 0$, but has a singularity at the origin for $-1/2 < m \le 0$; see Proposition 7. The tails of the distribution also become heavier if $-1/2 < n \le 0$; see Proposition 8 and Corollary 9, where in the corollary we give asymptotic approximations for tail probabilities. It follows from Proposition 8 that the mean of the VG ratio distribution is undefined.

## 2. The variance-gamma ratio distribution

### 2.1. Probability density function

In the following theorem, we provide an explicit expression for the PDF $f_Z(z)$ of $Z = X/Y$ in terms of an infinite series involving the Gaussian hypergeometric function. Throughout this paper, we let $\gamma_i^2 = \alpha_i^2 - \beta_i^2$, $i = 1, 2$. We also denote $a_{ij} = (1 + (-1)^{i+j})/2$, $i, j \ge 0$, so that $a_{ij} = 1$ if $i + j$ is even, and $a_{ij} = 0$ if $i + j$ is odd.

**Theorem 1.** Let $m, n > -1/2$, $0 \le |\beta_i| < \alpha_i$, $i = 1, 2$. Suppose $X \sim \text{VG}(m, \alpha_1, \beta_1, 0)$ and $Y \sim \text{VG}(n, \alpha_2, \beta_2, 0)$ are independent. Let $Z = X/Y$. Then, for $z \in \mathbb{R}$,

$$ f_Z(z) = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1} |z|^{-2n-2}}{\pi \alpha_1^{2m+2n+2} \Gamma(m+1/2) \Gamma(n+1/2)} \sum_{i,j=0}^{\infty} \frac{\beta_1^i \beta_2^j}{i! j!} \frac{2^{i+j} a_{ij} ((i+j)/2)! z^{-j}}{\alpha_1^{i+j} \Gamma(m+n+2+i+j)} \\ \times \Gamma\left(m+n+1+\frac{i+j}{2}\right) \Gamma\left(m+1+\frac{i+j}{2}\right) \Gamma\left(n+1+\frac{i+j}{2}\right) \\ \times {}_2F_1\left(m+n+1+\frac{i+j}{2}, n+1+\frac{i+j}{2}; m+n+2+i+j; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right). \quad (2) $$

**Remark 2.** When $\beta_1 = \beta_2 = 0$, the PDF $f_Z(z)$ can be expressed as a single hypergeometric function:

$$ f_Z(z) = \frac{\alpha_2^{2n+1} |z|^{-2n-2} \Gamma(m+1) \Gamma(n+1)}{\pi \alpha_1^{2n+1} (m+n+1) \Gamma(m+1/2) \Gamma(n+1/2)} {}_2F_1\left(m+n+1, n+1; m+n+2; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right), \quad (3) $$
in agreement with the formula given in Theorem 2.1 of [19]. We note that an application of the formula (26) yields the following further simplification when $m = n = 0$:

$$f_Z(z) = \frac{2\alpha_1}{\pi^2\alpha_2} \frac{\log|\alpha_1 z/\alpha_2|}{(\alpha_1 z/\alpha_2)^2 - 1}. \quad (4)$$

When $\beta_1 = 0$ the PDF $f_Z(z)$ simplifies to

$$f_Z(z) = \frac{\gamma_2^{2n+1}|z|^{-2n-2}}{\pi \alpha_1^{2n+1}} \sum_{k=0}^{\infty} \frac{k!}{(2k)!} \frac{\Gamma(m+n+1+k)\Gamma(m+1+k)\Gamma(n+1+k)}{\Gamma(m+1/2)\Gamma(n+1/2)\Gamma(m+n+2+2k)} \\ \times \left(\frac{2\beta_2}{\alpha_1 z}\right)^{2k} {}_2F_1\left(m+n+1+k, n+1+k; m+n+2+2k; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right),$$

and when $\beta_2 = 0$ the PDF simplifies to

$$f_Z(z) = \frac{\gamma_1^{2m+1} \alpha_2^{2n+1} |z|^{-2n-2}}{\pi \alpha_1^{2m+2n+1}} \sum_{k=0}^{\infty} \frac{k!}{(2k)!} \frac{\Gamma(m+n+1+k)\Gamma(m+1+k)\Gamma(n+1+k)}{\Gamma(m+1/2)\Gamma(n+1/2)\Gamma(m+n+2+2k)} \\ \times \left(\frac{2\beta_1}{\alpha_1}\right)^{2k} {}_2F_1\left(m+n+1+k, n+1+k; m+n+2+2k; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right).$$

Observe that the PDF $f_Z(z)$ is symmetric about the origin if $\beta_1 = 0$ or $\beta_2 = 0$.

We remark that the exact PDF of the product of correlated normal random variables is expressed as a double sum involving the modified Bessel function of the second kind when both normal random variables have non-zero mean [2], and the PDF simplifies to a single infinite series of modified Bessel functions of the second kind when one of the means is zero [2], and simply to the VG PDF when both means are zero [5]. Thus, the increase in complexity of the PDF of the VG ratio distribution for non-zero skewness parameters $\beta_1, \beta_2$ mirrors the increase in complexity of the product of two correlated normal random variables as one moves from zero means to non-zero means.

**Proof.** We consider the case $z > 0$; the case $z < 0$ is similar and thus omitted from this proof. The PDF of $Z = X/Y$ for $z > 0$ can be expressed as

$$f_Z(z) = \int_{-\infty}^{\infty} |y| f_X(yz) f_Y(y) dy = I_1 + I_2,$$

where

$$I_1 = M_1 M_2 \int_0^\infty y^{m+n+1} e^{\beta_1 yz} e^{\beta_2 y} z^m K_m(\alpha_1 yz) K_n(\alpha_2 y) dy, \quad (5)$$

$$I_2 = M_1 M_2 \int_{-\infty}^{0} (-y)^{m+n+1} e^{\beta_1 yz} e^{\beta_2 y} z^m K_m(-\alpha_1 yz) K_n(-\alpha_2 y) dy, \quad (6)$$

and $M_1 := M_{m,\alpha_1,\beta_1}$ and $M_2 := M_{n,\alpha_2,\beta_2}$ are the normalising constants of the VG($m, \alpha_1, \beta_1, 0$) and VG($n, \alpha_2, \beta_2, 0$) distributions. On using Taylor's expansion of the exponential function and
interchanging integration and summation in the first step, and evaluating the integral in the
second step using equation (6.576.4) of [9], we obtain that

$$
\begin{align}
I_1 &= M_1 M_2 z^m \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} \frac{\beta_1^i z^i \beta_2^j}{i! j!} \int_0^{\infty} y^{m+n+1+i+j} K_m(\alpha_1 yz) K_n(\alpha_2 y) dy \nonumber \\
&= \frac{M_1 M_2 2^{m+n-1} \alpha_2^n}{\alpha_1^{m+2n+2} z^{2n+2}} \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} \frac{\beta_1^i \beta_2^j}{i! j!} \frac{2^{i+j}}{\alpha_1^{i+j} z^j \Gamma(m+n+2+i+j)} \Gamma\left(m+n+1+\frac{i+j}{2}\right) \nonumber \\
&\qquad \times \Gamma\left(n+1+\frac{i+j}{2}\right) \Gamma\left(m+1+\frac{i+j}{2}\right) \Gamma\left(1+\frac{i+j}{2}\right) \nonumber \\
&\qquad \times {}_2F_1\left(m+n+1+\frac{i+j}{2}, 1+n+\frac{i+j}{2}; m+n+2+i+j; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right). \tag{7}
\end{align}
$$

By a similar calculation, we get that

$$
\begin{equation}
\begin{split}
I_2 = \frac{M_1 M_2 2^{m+n-1} \alpha_2^n}{\alpha_1^{m+2n+2} z^{2n+2}} & \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} \frac{(-1)^{i+j} \beta_1^i \beta_2^j}{i! j!} \frac{2^{i+j}}{\alpha_1^{i+j} z^j \Gamma(m+n+2+i+j)} \\
& \times \Gamma\left(m+n+1+\frac{i+j}{2}\right) \Gamma\left(n+1+\frac{i+j}{2}\right) \Gamma\left(m+1+\frac{i+j}{2}\right) \Gamma\left(1+\frac{i+j}{2}\right) \\
& \times {}_2F_1\left(m+n+1+\frac{i+j}{2}, n+1+\frac{i+j}{2}; m+n+2+i+j; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right).
\end{split}
\tag{8}
\end{equation}
$$

Summing up (7) and (8) now yields (2) for $z > 0$. $\square$

<figure><img src="image_3.png" /></figure>

Using special properties of the modified Bessel function of the second kind and the hyper-
geometric function, equivalent forms and elementary formulas for the PDF of Z = X/Y can be
derived; we illustrate this in the following Corollary 3 and Proposition 4.

**Corollary 3.** The PDF (2) can be expressed in the equivalent forms:

$$
\begin{equation}
\begin{aligned}
f_Z(z) = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1} |z|^{2m}}{\pi \alpha_2^{m+2n+2} \Gamma(m+1/2) \Gamma(n+1/2)} & \sum_{i,j=0}^{\infty} \frac{\beta_1^i \beta_2^j}{i! j!} \frac{2^{i+j} a_{ij} ((i+j)/2)! z^i}{\alpha_2^{i+j} \Gamma(m+n+2+i+j)} \\
& \times \Gamma\left(m+n+1+\frac{i+j}{2}\right) \Gamma\left(m+1+\frac{i+j}{2}\right) \Gamma\left(n+1+\frac{i+j}{2}\right) \\
& \times {}_2F_1\left(m+n+1+\frac{i+j}{2}, m+1+\frac{i+j}{2}; m+n+2+i+j; 1-\frac{\alpha_1^2 z^2}{\alpha_2^2}\right),
\end{aligned}
\tag{9}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
f_Z(z) = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1}}{\pi \alpha_1^{2m} \alpha_2^{2n+2} \Gamma(m+1/2) \Gamma(n+1/2)} & \sum_{i,j=0}^{\infty} \frac{\beta_1^i \beta_2^j}{i! j!} \frac{2^{i+j} a_{ij}((i+j)/2)! z^i}{\alpha_2^{i+j} \Gamma(m+n+2+i+j)} \\
& \times \Gamma\left(m+n+1+\frac{i+j}{2}\right) \Gamma\left(m+1+\frac{i+j}{2}\right) \Gamma\left(n+1+\frac{i+j}{2}\right) \\
& \times {}_2F_1\left(n+1+\frac{i+j}{2}, 1+\frac{i+j}{2}; m+n+2+i+j; 1-\frac{\alpha_1^2 z^2}{\alpha_2^2}\right),
\end{split}
\tag{10}
\end{equation}
$$
and

$$
\begin{equation}
\begin{aligned}
f_Z(z) = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1} z^{-2}}{\pi \alpha_1^{2m+2} \alpha_2^{2n} \Gamma(m+1/2) \Gamma(n+1/2)} & \sum_{i,j=0}^{\infty} \frac{\beta_1^i \beta_2^j}{i! j!} \frac{2^{i+j} a_{ij} ((i+j)/2)! z^{-j}}{\alpha_1^{i+j} \Gamma(m+n+2+i+j)} \\
& \times \Gamma\left(m+n+1+\frac{i+j}{2}\right) \Gamma\left(m+1+\frac{i+j}{2}\right) \Gamma\left(n+1+\frac{i+j}{2}\right) \\
& \times {}_2F_1\left(1+\frac{i+j}{2}, m+1+\frac{i+j}{2}; m+n+2+i+j; 1-\frac{\alpha_2^2}{\alpha_1^2 z^2}\right).
\end{aligned}
\tag{11}
\end{equation}
$$

Note that (9) and (10) hold for $-\sqrt{2\alpha_2}/\alpha_1 < z < \sqrt{2\alpha_2}/\alpha_1$, whilst (11) holds for $z < -\alpha_2/(\sqrt{2}\alpha_1)$ or $z > \alpha_2/(\sqrt{2}\alpha_1)$.

**Proof.** Apply the three transformation formulas given in equation (9.131.1) of [9] to the hyperge-
ometric functions in the formula (2). $\square$

**Proposition 4.** *Suppose that $m - 1/2 \ge 0$ and $n - 1/2 \ge 0$ are integers. Then, (2) can be simplified to the elementary form:*

$$
f_Z(z) = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1} |z|^{m-1/2}}{(2\alpha_1)^{m+1/2} (2\alpha_2)^{n+1/2}} \sum_{i=0}^{m-1/2} \sum_{j=0}^{n-1/2} \frac{(m+n-i-j)!}{(m-1/2-i)!(n-1/2-j)!} \binom{m-1/2+i}{i} \\
\times \binom{n-1/2+j}{j} (2\alpha_1|z|)^{-i} (2\alpha_2)^{-j} \left( \frac{1}{u_1^{m+n+1-i-j}} + \frac{1}{u_2^{m+n+1-i-j}} \right), \quad (12)
$$

where $u_1 = \alpha_1|z| + \alpha_2 + \beta_1z + \beta_2$ and $u_2 = \alpha_1|z| + \alpha_2 - \beta_1z - \beta_2$.

**Remark 5.** As noted by [13], the VG(1/2, α, β, 0) distribution corresponds to the asymmetric Laplace distribution. Hence, setting *m* = *n* = 1/2 in (12) yields the PDF of the ratio of independent asymmetric Laplace random variables. In the case *β*₁ = *β*₂ = 0, the formula further reduces to the PDF of the ratio of independent Laplace random variables; this distribution is known as the Lomax distribution (see, for example, [12]).

**Proof.** Apply the formula (24) to the modified Bessel functions of the second kind that appear in the integrands of $I_1$ and $I_2$ given in (5) and (6). Evaluating the resulting integrals using the standard formula $\int_0^\infty x^a e^{-bx} dx = a!b^{-a-1}$ for $a = 0, 1, 2, \dots, b > 0$, summing up $I_1 + I_2$, and simplifying yields (12). $\square$

**2.2. Cumulative distribution function**

Let $F_Z(z) = \mathbb{P}(Z \le z)$ denote the CDF of the VG ratio distribution. In the following theorem, we provide an exact formula for the CDF in the case of a ratio of independent symmetric VG random variables $(\beta_1 = \beta_2 = 0)$, the same setting as in the work of [19].

**Theorem 6.** Let $\beta_1 = \beta_2 = 0$, and $m, n > -1/2$, $\alpha_1, \alpha_2 > 0$. Then, for $z \in \mathbb{R}$,

$$
F_Z(z) = \frac{1}{2} + \frac{\alpha_1 z / \alpha_2}{2\pi \Gamma(m+1/2)\Gamma(n+1/2)} G_{3,3}^{2,3} \left( \frac{\alpha_1^2 z^2}{\alpha_2^2} \middle| \begin{array}{c} -n, \frac{1}{2}, 0 \\ m, 0, -\frac{1}{2} \end{array} \right). \quad (13)
$$

**Proof.** For ease of exposition, we set $\alpha_1 = \alpha_2 = 1$; the general case follows from a simple rescaling.
As the PDF is symmetric when $\beta_1 = \beta_2 = 0$, we also just consider the case $z > 0$. Using the
representation (33) of the hypergeometric function in terms of the Meijer G-function allows us to
write the PDF (3) in terms of the Meijer G-function:

$$
f_Z(z) = C_{m,n} z^{-2n-2} G_{2,2}^{2,2} \left( z^{-2} \left| \begin{array}{c} -m-n, -n \\ 0, -n \end{array} \right. \right) = C_{m,n} z^{-2n-2} G_{2,2}^{2,2} \left( z^2 \left| \begin{array}{c} 1, n+1 \\ m+n+1, n+1 \end{array} \right. \right),
$$
where $1/C_{m,n} = \pi\Gamma(m+1/2)\Gamma(n+1/2)$, and we used the relation (31) to obtain the second equality. For $z > 0$, the CDF of Z is thus given by

$$
\begin{align*}
F_Z(z) &= \frac{1}{2} + C_{m,n} \int_0^z y^{-2n-2} G_{2,2}^{2,2} \left( y^2 \middle| \begin{array}{c} 1, n+1 \\ m+n+1, n+1 \end{array} \right) dy \\
&= \frac{1}{2} + \frac{C_{m,n}}{2} \int_0^{z^2} u^{-n-3/2} G_{2,2}^{2,2} \left( u \middle| \begin{array}{c} 1, n+1 \\ m+n+1, n+1 \end{array} \right) du \\
&= \frac{1}{2} + \frac{C_{m,n}}{2} \left[ u^{-n-1/2} G_{3,3}^{2,3} \left( u \middle| \begin{array}{c} n+\frac{3}{2}, 1, n+1 \\ m+n+1, n+1, n+\frac{1}{2} \end{array} \right) \right]_0^{z^2} \\
&= \frac{1}{2} + \frac{C_{m,n}}{2} \left[ G_{3,3}^{2,3} \left( u \middle| \begin{array}{c} 1, \frac{1}{2}-n, \frac{1}{2} \\ m+\frac{1}{2}, \frac{1}{2}, 0 \end{array} \right) \right]_0^{z^2} = \frac{1}{2} + \frac{C_{m,n}}{2} G_{3,3}^{2,3} \left( z^2 \middle| \begin{array}{c} 1, \frac{1}{2}-n, \frac{1}{2} \\ m+\frac{1}{2}, \frac{1}{2}, 0 \end{array} \right),
\end{align*}
$$

where we used the integral formula (34) in the third step; the relation (32) in the fourth step; and in the final step we used that the Meijer G-function in the penultimate equality evaluated at $u = 0$ is equal to zero, which is readily deduced from the contour integral representation (30) of the G-function. Finally, we using (32) yields

$$
F_Z(z) = \frac{1}{2} + \frac{C_{m,n}}{2} \cdot z G_{3,3}^{2,3} \left( z^2 \middle| \begin{array}{c} \frac{1}{2}, -n, 0 \\ m, 0, -\frac{1}{2} \end{array} \right) = \frac{1}{2} + \frac{C_{m,n}}{2} \cdot z G_{3,3}^{2,3} \left( z^2 \middle| \begin{array}{c} -n, \frac{1}{2}, 0 \\ m, 0, -\frac{1}{2} \end{array} \right). \quad (14)
$$

The formula (14) is also readily seen to hold for $z < 0$. We therefore arrive at the formula (13), which conveniently holds for all $z \in \mathbb{R}$ □

2.3. Asymptotic behaviour of the density and tail probabilities, and fractional moments

The representations of the PDF $f_Z(z)$ given in Theorem 1, Corollary 3 and Proposition 4 are rather complicated and difficult to parse at first inspection. Some insight can be gained from the following propositions.

**Proposition 7.** Let $m, n > -1/2, 0 \le |\beta_i| < \alpha_i, i = 1, 2$.

(1) Suppose $m > 0$. Then

$$
f_Z(0) = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1} \Gamma(m) \Gamma(n+1)}{\pi \alpha_1^{2m} \alpha_2^{2n+2} \Gamma(m+1/2) \Gamma(n+1/2)} {}_2F_1 \left( 1, n+1; \frac{1}{2}; \frac{\beta_2^2}{\alpha_2^2} \right).
$$

(2) Suppose $m = 0$. Then, as $z \to 0$,

$$
f_Z(z) \sim - \frac{2\gamma_1\gamma_2^{2n+1}}{\pi^{3/2}\alpha_2^{2n+2}} \frac{\Gamma(n+1)\log|z|}{\Gamma(n+1/2)} {}_2F_1 \left( 1, n+1; \frac{1}{2}; \frac{\beta_2^2}{\alpha_2^2} \right).
$$

(3) Suppose $-1/2 < m < 0$. Then, as $z \to 0$,

$$
f_Z(z) \sim \frac{\gamma_1^{2m+1} \gamma_2^{2n+1}}{\alpha_2^{2m+2n+2} \sin(-m\pi)} \frac{\Gamma(m+n+1)|z|^{2m}}{\Gamma(m+1/2)\Gamma(n+1/2)} {}_2F_1 \left( m+1, m+n+1; \frac{1}{2}; \frac{\beta_2^2}{\alpha_2^2} \right).
$$

**Proposition 8.** Let $m, n > -1/2, 0 \le |\beta_i| < \alpha_i, i = 1, 2$.

(1) Suppose $n > 0$. Then, as $|z| \to \infty$,

$$
f_Z(z) \sim \frac{\gamma_1^{2m+1} \gamma_2^{2n+1} \Gamma(m+1) \Gamma(n) z^{-2}}{\pi \alpha_1^{2m+2} \alpha_2^{2n} \Gamma(m+1/2) \Gamma(n+1/2)} {}_2F_1 \left( 1, m+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2} \right). \quad (15)
$$

(2) Suppose $n = 0$. Then, as $|z| \to \infty$,

$$
f_Z(z) \sim \frac{2\gamma_1^{2m+1}\gamma_2}{\pi^{3/2}\alpha_1^{2m+2}} \frac{\Gamma(m+1)z^{-2}\log|z|}{\Gamma(m+1/2)} {}_2F_1\left(1, m+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right). \quad (16)
$$
(3) Suppose $-1/2 < n < 0$. Then, as $|z| \to \infty$,

$$f_Z(z) \sim \frac{\gamma_1^{2m+1} \gamma_2^{2n+1}}{\alpha_1^{2m+2n+2} \sin(-n\pi)} \frac{\Gamma(m+n+1)|z|^{-2-2n}}{\Gamma(m+1/2)\Gamma(n+1/2)} {}_2F_1\left(n+1, m+n+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right). \quad (17)$$

**Corollary 9.** Let $m, n > -1/2$, $0 \le |\beta_i| < \alpha_i$, $i = 1, 2$. Let $\bar{F}_Z(z) = \mathbb{P}(Z > z)$.

(1) Suppose $n > 0$. Then, as $z \to \infty$,

$$\bar{F}_Z(z) \sim \frac{\gamma_1^{2m+1} \gamma_2^{2n+1}}{\pi \alpha_1^{2m+2} \alpha_2^{2n}} \frac{\Gamma(m+1)\Gamma(n)z^{-1}}{\Gamma(m+1/2)\Gamma(n+1/2)} {}_2F_1\left(1, m+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right). \quad (18)$$

(2) Suppose $n = 0$. Then, as $z \to \infty$,

$$\bar{F}_Z(z) \sim \frac{2\gamma_1^{2m+1}\gamma_2}{\pi^{3/2}\alpha_1^{2m+2}} \frac{\Gamma(m+1)z^{-1}\log(z)}{\Gamma(m+1/2)} {}_2F_1\left(1, m+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right).$$

(3) Suppose $-1/2 < n < 0$. Then, as $z \to \infty$,

$$\bar{F}_Z(z) \sim \frac{\gamma_1^{2m+1} \gamma_2^{2n+1}}{\alpha_1^{2m+2n+2} \sin(-n\pi)} \frac{\Gamma(m+n+1)z^{-1-2n}}{(1+2n)\Gamma(m+1/2)\Gamma(n+1/2)} {}_2F_1\left(n+1, m+n+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right).$$

**Remark 10.** Since ${}_2F_1(a, b; c; 0) = 1$, the formulas of Proposition 7 simplify when $\beta_2 = 0$, whilst the formulas of Proposition 8 and Corollary 9 simplify when $\beta_1 = 0$.

**Remark 11.** Proposition 7 tells us that, as $z \to 0$, the PDF $f_Z(z)$ is bounded for $m > 0$, but has a singularity for $-1/2 < m \le 0$. Note that whether the density is bounded or unbounded is determined solely by the parameter $m$. We remark that this behaviour is consistent with the behaviour of the VG density $f_X(x)$, which is bounded for all $x \in \mathbb{R}$ if $m > 0$, but has a singularity as $x \to 0$ if $-1/2 < m \le 0$ (see [7, equation (2.4)]).

The asymptotic behaviour of the PDF as $z \to 0$ shows that the VG ratio distribution is unimodal if $-1/2 < m \le 0$, as in this case the density is bounded everywhere except for the singularity at $z = 0$. Also, the VG ratio distribution was shown to be unimodal in the $\beta_1 = \beta_2 = 0$ case by [19] (they proved it for $m, n \ge 1$, although their argument also applies for $m, n > -1/2$). Our numerical tests strongly suggest that the VG ratio distribution is unimodal for all parameter values, although we have been unable to prove this.

We see from Proposition 8 that the tails of the VG ratio distribution become heavier as $n \le 0$ decreases; note that now the rate of decay of the tails with respect to $z$ is solely determined by the parameter $n$. We remark that it is curious that, whilst the VG ratio distribution is only symmetric around the origin if $\beta_1 = 0$ or if $\beta_2 = 0$, for all parameter values we have $\lim_{z \to \infty} f_Z(z)/f_Z(-z) = 1$. This is in contrast to the VG distribution itself, for which $\lim_{x \to \infty} f_X(x)/f_X(-x) = 1$ if and only if $\beta = 0$ (see [7, p. 4]).

It follows from Proposition 8 that, for all possible parameter values, the mean of the VG ratio distribution is undefined. This is a common feature of ratio distributions; for example, it is well-known that the ratio of independent standard normal random variables follows the Cauchy distribution, for which the mean does not exist. Whilst the mean of the VG ratio distribution is not defined, the fractional moments $\mathbb{E}[|Z|^k]$ do exist for suitable $k < 1$, as seen in Proposition 12 below.

The PDF of the VG ratio distribution is plotted for a range of parameter values in Figures 1 and 2. Figure 1 demonstrates the effect of varying the shape parameters $m$ and $n$, whilst Figure 2 shows the effect of varying the skewness parameters $\beta_1$ and $\beta_2$. In both figures, $\alpha_1 = \alpha_2 = 1$. The figures confirm the assertions made in Remark 11.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** $\beta_1 = \beta_2 = 0$, (a): $m = -0.25$; (b): $m = 0$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** $m = n = 1.5$, $\beta_1 = 0.5$.</figcaption></figure>

**Proof of Proposition 7.** Suppose first that $m > 0$. We shall make use of the representation (10) of the PDF $f_Z(z)$. Let

$$C = \frac{\gamma_1^{2m+1} \gamma_2^{2n+1}}{\pi \alpha_1^{2m} \alpha_2^{2n+2} \Gamma(m+1/2) \Gamma(n+1/2)}$$

Then

$$\lim_{z \to 0} f_Z(z) = C \sum_{k=0}^{\infty} \left(\frac{2\beta_2}{\alpha_2}\right)^{2k} \frac{k! \Gamma(m+n+1+k) \Gamma(m+1+k) \Gamma(n+1+k)}{(2k)! \Gamma(m+n+2+2k)} \times {}_2F_1(n+1+k, 1+k; m+n+2+2k; 1). \quad (19)$$

Using the formula (27) to express the hypergeometric functions in (19) in terms of the gamma function, and cancelling the gamma functions gives that

$$\begin{aligned} \lim_{z \to 0} f_Z(z) &= C\Gamma(m) \sum_{k=0}^{\infty} \frac{k!}{(2k)!} \Gamma(n+1+k) \left(\frac{2\beta_2}{\alpha_2}\right)^{2k} \\ &= C\Gamma(m)\Gamma(n+1){}_2F_1\left(1, n+1; \frac{1}{2}; \frac{\beta_2^2}{\alpha_2^2}\right), \end{aligned}$$

where we obtained the second equality by using the basic formula $(u)_k = \Gamma(u+k)/\Gamma(u)$ to put the infinite series into the hypergeometric form (25).

The cases $m = 0$ and $-1/2 < m < 0$ are treated similarly, but instead of using the equality (27) we apply the limiting forms (28) and (29) for the hypergeometric function, respectively. In the case $-1/2 < m < 0$, we simplify our limiting form for $f_Z(z)$ by applying the standard formula $\Gamma(-m)\Gamma(m+1) = \pi/\sin(-m\pi)$. $\square$

**Proof of Proposition 8.** The proof is similar to that of Proposition 7, except we make use of the representation (11) of the PDF $f_Z(z)$ instead of the representation (10). $\square$
**Proof of Corollary 9.** Suppose first that $n > 0$. Then, by (15), as $z \to \infty$,

$$ \bar{F}_Z(z) \sim \frac{\gamma_1^{2m+1}\gamma_2^{2n+1}}{\pi\alpha_1^{2m+2}\alpha_2^{2n}} \frac{\Gamma(m+1)\Gamma(n)}{\Gamma(m+1/2)\Gamma(n+1/2)} {}_2F_1\left(1, m+1; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right) \int_z^\infty x^{-2} \mathrm{d}x, $$

and evaluating $\int_z^\infty x^{-2} \mathrm{d}x = z^{-1}$ yields (18). The cases $n=0$ and $-1/2 < n < 0$ are dealt with similarly making use of the limiting forms (16) and (17) and the formulas $\int_z^\infty x^{-2} \log(x) \mathrm{d}x = z^{-1}(\log(z)+1) \sim z^{-1}\log(z)$, as $z \to \infty$, and $\int_z^\infty x^{-2-2n} \mathrm{d}x = z^{-1-2n}/(1+2n)$. $\square$

**Proposition 12.** Let $m, n > -1/2$, $0 \le |\beta_i| < \alpha_i$, $i = 1, 2$. Then, for $\max\{-1, -2m-1\} < k < \min\{1, 2n+1\}$,

$$ \mathbb{E}[|Z|^k] = \frac{\alpha_2^k (1 - \beta_1^2 / \alpha_1^2)^{m+1/2} (1 - \beta_2^2 / \alpha_2^2)^{n+1/2}}{\alpha_1^k \cos(k\pi/2) \Gamma(m+1/2) \Gamma(n+1/2)} \Gamma\left(m + \frac{k+1}{2}\right) \Gamma\left(n + \frac{1-k}{2}\right) \\ \times {}_2F_1\left(\frac{k+1}{2}, m + \frac{k+1}{2}; \frac{1}{2}; \frac{\beta_1^2}{\alpha_1^2}\right) {}_2F_1\left(\frac{1-k}{2}, n + \frac{1-k}{2}; \frac{1}{2}; \frac{\beta_2^2}{\alpha_2^2}\right). $$

We will need the following formula of [8] for the absolute moments of the VG($m, \alpha, \beta, 0$) distribution. Then, for $k > \max\{-1, -2m-1\}$,

$$ \mathbb{E}[|X|^k] = \frac{2^k (1 - \beta^2 / \alpha^2)^{m+1/2}}{\sqrt{\pi} \alpha^k \Gamma(m+1/2)} \Gamma\left(m + \frac{k+1}{2}\right) \Gamma\left(\frac{k+1}{2}\right) {}_2F_1\left(\frac{k+1}{2}, m + \frac{k+1}{2}; \frac{1}{2}; \frac{\beta^2}{\alpha^2}\right). \quad (20) $$

**Proof.** By independence, $\mathbb{E}[|Z|^k] = \mathbb{E}[|X|^k]\mathbb{E}[|Y|^{-k}]$, and calculating the absolute moments $\mathbb{E}[|X|^k]$ and $\mathbb{E}[|Y|^{-k}]$ using equation (20) yields the desired formula for $\mathbb{E}[|Z|^k]$ following a simplification using the identity $\Gamma((k+1)/2)\Gamma((1-k)/2) = \pi \sec(k\pi/2)$. $\square$

## 2.4. 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 $(\sigma_{U_i}^2, \sigma_{V_i}^2)$ and correlation coefficient $\rho_i$. Set $s_i = \sigma_{U_i} \sigma_{V_i}$, and let $W_i = U_i V_i$. It was noted by [3, 5] that the product $W_i$ has a VG distribution,

$$ W_i \sim \text{VG}\left(0, \frac{1}{s_i(1 - \rho_i^2)}, \frac{\rho_i}{s_i(1 - \rho_i^2)}, 0\right), \quad (21) $$

which yielded a new derivation of the PDF of the product $W_i$ which was earlier obtained independently by [10] and [20]. Let $T = W_1/W_2$. On combining (2) and (4) with (21) we immediately deduce the following formulas for the PDF of the ratio $T$.

**Corollary 13.** Let the previous notations prevail. Then, for $t \in \mathbb{R}$,

$$ f_T(t) = \frac{(s_1/s_2)(1-\rho_1^2)^{3/2} t^{-2}}{\pi^2 \sqrt{1-\rho_2^2}} \sum_{i,j=0}^{\infty} \frac{a_{ij}((i+j)/2)!}{i!j!(i+j+1)!} (2\rho_1)^i \left(\frac{1-\rho_1^2}{1-\rho_2^2}\right)^j \left(\frac{2s_1\rho_2}{s_2t}\right)^j \\ \times {}_2F_1\left(1+\frac{i+j}{2}, 1+\frac{i+j}{2}; 2+i+j; 1-\frac{s_1(1-\rho_1^2)}{s_2(1-\rho_2^2)t}\right). \quad (22) $$

The PDF (22) simplifies as follows when $\rho_1 = \rho_2 = 0$. For $t \in \mathbb{R}$,

$$ f_T(t) = \frac{2s_2 \log|s_2 t/s_1|}{\pi^2 s_1 (s_2 t/s_1)^2 - 1}. \quad (23) $$

**Remark 14.** When $s_1 = s_2 = 1$, (23) is also the density of the product of two independent standard Cauchy random variables with PDF $f(x) = \pi^{-1}(1+x^2)^{-1}$, $x \in \mathbb{R}$. This follows from the well-known result that the quotient of independent standard normal random variables follows the standard Cauchy distribution. More generally, formulas for the PDF of the product of independent Student's $t$ random variables are given in [18].
## Appendix A. Special functions

In this appendix, we define the modified Bessel function of the second kind, the Gaussian hypergeometric function and the Meijer G-function, and state some basic properties that are needed in this paper. Unless otherwise stated, these and further properties can be found in the standard references [9, 14, 21].

The modified Bessel function of the second kind is defined, for $m \in \mathbb{R}$ and $x > 0$, by

$$K_m(x) = \int_0^\infty e^{-x \cosh(t)} \cosh(mt) \, dt.$$

For $m = 0, 1, 2, \dots$, we have the elementary representation

$$K_{m+1/2}(x) = \sqrt{\frac{\pi}{2x}} \sum_{j=0}^{m} \frac{(m+j)!}{(m-j)!j!} (2x)^{-j} e^{-x} \quad (24)$$

The Gaussian hypergeometric function is defined by the power series

$${}_2F_1(a, b; c; x) = \sum_{j=0}^{\infty} \frac{(a)_j (b)_j}{(c)_j j!} x^j, \quad (25)$$

for $|x| < 1$, and by analytic continuation elsewhere. Here $(u)_j = u(u+1)\cdots(u+k-1)$ is the ascending factorial.

A special case is the following:

$${}_2F_1(1, 1; 2; x) = -x^{-1} \log(1-x). \quad (26)$$

The hypergeometric function has the following behaviour as $x \to 1^{-}$.

$$(1) \quad c-a-b>0:$$

$${}_2F_1(a, b; c; 1) = \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}. \quad (27)$$

$$(2) \quad c-a-b=0:$$

$${}_2F_1(a, b; a+b; x) \sim -\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \log(1-x), \quad x \to 1^{-}. \quad (28)$$

$$(3) \quad c-a-b<0:$$

$${}_2F_1(a, b; c; x) \sim \frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)} (1-x)^{c-a-b}, \quad x \to 1^{-}. \quad (29)$$

The Meijer G-function is defined, for $x \in \mathbb{R}$, by the contour integral

$$G_{p,q}^{m,n} \left( x \middle| \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \right) = \frac{1}{2\pi i} \int_L \frac{\prod_{j=1}^m \Gamma(b_j-s) \prod_{j=1}^n \Gamma(1-a_j+s)}{\prod_{j=n+1}^p \Gamma(a_j-s) \prod_{j=m+1}^q \Gamma(1-b_j+s)} x^s \, ds, \quad (30)$$

where the integration path $L$ separates the poles of the factors $\Gamma(b_j - s)$ from those of the factors $\Gamma(1 - a_j + s)$, and we use the convention that the empty product is 1.

The G-function satisfies the relations

$$G_{p,q}^{m,n} \left( x \middle| \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \right) = G_{q,p}^{n,m} \left( x^{-1} \middle| \begin{matrix} 1-b_1, \dots, 1-b_q \\ 1-a_1, \dots, 1-a_p \end{matrix} \right), \quad (31)$$

$$x^{\alpha} G_{p,q}^{m,n} \left( x \middle| \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \right) = G_{p,q}^{m,n} \left( x \middle| \begin{matrix} a_1 + \alpha, \dots, a_p + \alpha \\ b_1 + \alpha, \dots, b_q + \alpha \end{matrix} \right). \quad (32)$$

The hypergeometric function can be expressed in terms of the Meijer G-function:

$${}_2F_1(a, b; c; 1-x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)\Gamma(c-a)\Gamma(c-b)} G_{2,2}^{2,2} \left( x \middle| \begin{matrix} 1-a, 1-b \\ 0, c-a-b \end{matrix} \right). \quad (33)$$
(see http://functions.wolfram.com/07.23.26.0007.01). Also, on combining equations 5.4 (1) and 5.4 (13) of [14], we obtain the indefinite integral formula

$$ \int x^{\alpha-1} G_{p,q}^{m,n} \left( x \mid \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \right) dx = x^{\alpha} G_{p+1,q+1}^{m,n+1} \left( x \mid \begin{matrix} 1-\alpha, a_1, \dots, a_p \\ b_1, \dots, b_q, -\alpha \end{matrix} \right). \quad (34) $$

## References

[1] C. C. Craig, "On the frequency function of xy", Ann. Math. Stat. 7 (1936), p. 1-15.

[2] G. Cui, X. Yu, S. Iommelli, L. Kong, "Exact distribution for the product of two correlated gaussian random variables", IEEE Sel. Top. Signal Process. 23 (2016), p. 1662-1666.

[3] R. E. Gaunt, "Rates of Convergence of Variance-Gamma Approximations via Stein's Method", PhD Thesis, University of Oxford, 2013.

[4] ——, "Variance-gamma approximation via Stein's method", Electron. J. Probab. 19 (2014), article no. 38 (33 pages).

[5] ——, "A note on the distribution of the product of zero mean correlated normal random variables", Stat. Neerl. 73 (2019), no. 2, p. 176-179.

[6] ——, "The basic distributional theory for the product of zero mean correlated normal random variables", Stat. Neerl. 76 (2022), p. 450-470.

[7] ——, "Stein factors for variance-gamma approximation in the Wasserstein and Kolmogorov distances", J. Math. Anal. Appl. 514 (2022), no. 1, article no. 126274 (32 pages).

[8] ——, "On the moments of the variance-gamma distribution", Stat. Probab. Lett. (2023), article no. 109884 (4 pages).

[9] I. S. Gradshteyn, I. M. Ryzhik, Table of Integrals, Series and Products, 7th ed., Academic Press Inc., 2007.

[10] L. P. Grishchuk, "Statistics of the microwave background anisotropies caused by the squeezed cosmological perturbations", Phys. Rev. D 53 (1996), p. 6784-6795.

[11] H. Holm, M.-S. Alouini, "Sum and difference of two squared correlated Nakagami variates with the McKay distribution", IEEE Trans. Commun. 52 (2004), no. 8, p. 1367-1376.

[12] N. L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distribution. Vol. 1, John Wiley & Sons, 1994.

[13] S. Kotz, T. J. Kozubowski, K. Podgórski, The Laplace distribution and generalizations. A revisit with applications to communications, economics, engineering, and finance, Birkhäuser, 2001.

[14] Y. L. Luke, The special functions and their approximations. Vol. I, Mathematics in Science and Engineering, vol. 53, Academic Press Inc., 1969.

[15] D. B. Madan, P. P. Carr, E. C. Chang, "The variance gamma process and option pricing", Eur. Finance Rev. 2 (1998), no. 1, p. 74-105.

[16] D. B. Madan, E. Seneta, "The Variance Gamma (V.G.) Model for Share Market Returns", J. Bus. 63 (1990), no. 4, p. 511-524.

[17] A. T. McKay, "A Bessel function distribution", Biometrika 24 (1932), p. 39-44.

[18] S. Nadarajah, "Exact distribution of the product of N Student's t RVs", Methodol. Comput. Appl. Probab. 14 (2012), no. 4, p. 997-1009.

[19] S. Nadarajah, S. Kotz, "The Bessel ratio distribution", C. R. Acad. Sci. Paris 343 (2006), no. 8, p. 531-534.

[20] S. Nadarajah, T. K. Pogány, "On the distribution of the product of correlated normal random variables", C. R. Acad. Sci. Paris 354 (2016), no. 2, p. 201-204.

[21] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark (eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, 2010.

[22] E. Seneta, "Fitting the variance-gamma model to financial data", J. Appl. Probab. 41A (2004), p. 177-187.