Product of two normally distributed random variables: Asymptotic approximation of tail probability

Question

Suppose $X \sim N(1,1)$, $Y \sim N(1,1)$ are iid normal random variables. My research problem is finding out the asymptotics of the tail function of XY (since the explicit formula is too complicated). If I'm not mistaken, the expression should look like: $$\begin{align}&\mathbb{P}(XY > z) = \mathbb{P}(XY > z, X > 0) + \mathbb{P}(XY > z, X < 0) \\&= \mathbb{P}\left(Y > \frac{z}{X}, X > 0\right) + \mathbb{P}\left(Y < \frac{z}{X}, X < 0\right) \\&= \int_0^\infty f_X(x) \left( \int_{\frac{z}{x}}^\infty f_Y(y) \, dy \right) dx + \int_{-\infty}^0 f_X(x) \left( \int_{-\infty}^{\frac{z}{x}} f_Y(y) \, dy \right) dx \\&=\int_0^\infty f(x) \overline{F}\left(\frac{z}{x}\right) dx + \int_{-\infty}^0 f(x) {F}\left(\dfrac{z}{x}\right) dx :=I_1+I_2. \end{align}$$ Using the well-known asymptotics of $\overline{F}$ for integral $I_1$ we get $$ I_1=\frac{1}{2\pi}\int_0^\infty \left(\frac{z}{x}-1\right)^{-1} e^{-\frac{1}{2}((x-1)^2+(\frac{z}{x}-1)^2)} dx. $$ Basically, I want to analyze this integral. The main issue is that I'm unsure which method I should use. The saddle point method this time doesn't work to my concern. Maybe someone has some ideas?

Answer 1

It is better to use the standard transformation: the variables $X_1 = (X + Y)/2\sim N(1,1/2)$, $X_2 = (X - Y)/2\sim N(0,1/2)$ are also independent and $XY = X_1^2 - X_2^2$, which allows to consider a sum instead of product. This simple idea is the most important ingredient, so I'll omit the technical detail. Rewrite $$ \mathrm{P}(XY>z) = \mathrm{P}(X_1^2 - X_2^2>z) = \frac{2}{\sqrt{\pi}}\int_0^\infty \mathrm{P}(X_1^2 >z + t^2) e^{-t^2}dt \\ = \frac{2}{\sqrt{\pi}}\int_0^\infty \left(\mathrm{P}(X_1 -1 >\sqrt{z + t^2}-1) + \mathrm{P}(1-X_1 >\sqrt{z + t^2}+1) \right) e^{-t^2}dt $$ Note that $X_1-1$ is centered, so we can use the standard asymptotics. The second probability is negligible compared to the first one, which is $$ \mathrm{P}(1-X_1 >\sqrt{z + t^2}-1)\sim \frac{e^{-(\sqrt{z + t^2}-1)^2}}{2\sqrt{\pi}(\sqrt{z + t^2}-1)}\sim \frac{e^{-z - t^2 + 2\sqrt{z} - 1}}{2\sqrt{\pi z}}, z\to\infty. $$ As a result, $$ \mathrm{P}(XY>z) \sim \frac{1}{\pi \sqrt{z}} e^{-z + 2\sqrt{z} -1} \int_0^\infty e^{-2 t^2}dt = \frac{1}{2e \sqrt{2\pi }} z^{-1/2}e^{-z+2\sqrt{z}}, z\to \infty, $$ which agrees with @heather-milkim's answer.

Answer 2

Your approach does work.

1.

We are considering two independent random variables $X, Y \sim N(1,1)$. We want to find the asymptotic behavior of the tail probability $P(XY > z)$ as $z \to \infty$. The probability can be decomposed based on the sign of $X$: $$ P(XY > z) = P(XY > z, X > 0) + P(XY > z, X < 0) $$ $$ P(XY > z) = P(Y > z/X, X > 0) + P(Y < z/X, X < 0) $$ Let $f(x) = \frac{1}{\sqrt{2\pi}} e^{-(x-1)^2/2}$ be the probability density function of $X$ (and $Y$). Let $F(y)$ be the cumulative distribution function (CDF) of $Y$, and $\overline{F}(y) = 1 - F(y)$ be the tail function. In terms of the standard normal CDF $\Phi$ and tail function $\overline{\Phi}$, we have $F(y) = \Phi(y-1)$ and $\overline{F}(y) = \overline{\Phi}(y-1)$. The probability can be written as: $$ P(XY > z) = \int_0^\infty f(x) \overline{F}(z/x) dx + \int_{-\infty}^0 f(x) F(z/x) dx $$ $$ P(XY > z) = \underbrace{\int_0^\infty f(x) \overline{\Phi}(z/x - 1) dx}_{I_1(z)} + \underbrace{\int_{-\infty}^0 f(x) \Phi(z/x - 1) dx}_{I_2(z)} $$

2.

As $z \to \infty$, the tail probability has the following asymptotic behavior: $$ P(XY > z) \sim \frac{e^{-1}}{2 \sqrt{2\pi}} z^{-1/2} e^{-z + 2\sqrt{z}} $$

3. Analysis of $I_1(z)$ (Main Contribution)

This integral represents the contribution from $X > 0$. $$ I_1(z) = \int_0^\infty \frac{1}{\sqrt{2\pi}} e^{-\frac{(x-1)^2}{2}} \overline{\Phi}(z/x - 1) dx $$ We use the asymptotic expansion for the standard normal tail function: for $t \to \infty$, $$ \overline{\Phi}(t) = \frac{\phi(t)}{t} (1 + O(t^{-2})) = \frac{1}{\sqrt{2\pi} t} e^{-t^2/2} (1 + O(t^{-2})) $$ Let $t(x) = z/x - 1$. Substituting the expansion into $I_1(z)$: $$ I_1(z) = \int_0^z \frac{1}{\sqrt{2\pi}} e^{-\frac{(x-1)^2}{2}} \frac{1}{\sqrt{2\pi} t(x)} e^{-\frac{t(x)^2}{2}} (1 + R(x,z)) dx + \int_z^\infty f(x) \overline{\Phi}(t(x)) dx $$ where $R(x,z) = O(t(x)^{-2}) = O((z/x-1)^{-2})$ for $t(x) > 0$ (i.e., $x < z$). The integral from $z$ to $\infty$ involves $t(x) < 0$; $\overline{\Phi}$ is bounded, and $f(x)$ decays super-exponentially, making this second integral negligible (it decays faster than any power of $z$). We focus on the first integral, denoted $J(z)$. $$ J(z) = \int_0^z \frac{1}{2\pi (z/x - 1)} e^{-\frac{1}{2} [(x-1)^2 + (z/x - 1)^2]} (1 + R(x,z)) dx $$ Let $g(x,z) = (x-1)^2 + (z/x - 1)^2$ and $h(x,z) = \frac{1}{2\pi (z/x - 1)}$. $$ J(z) = \int_0^z h(x,z) e^{-\frac{1}{2} g(x,z)} (1 + R(x,z)) dx $$ We apply Laplace's method. First, find the minimum of $g(x,z)$ for $x > 0$. The derivative is $g'(x,z) = 2(x-1) - \frac{2z}{x^3}(z-x)$. Setting $g'(x,z)=0$ leads to $x^4 - x^3 + zx - z^2 = 0$. As shown previously, $x_0 = \sqrt{z}$ is the unique positive solution for large $z$. The relevant values at the minimum $x_0 = \sqrt{z}$ are:

• $g(x_0, z) = 2(\sqrt{z}-1)^2 = 2z - 4\sqrt{z} + 2$
• $g'(x_0, z) = 0$
• $g''(x_0, z) = 8 - 4/\sqrt{z}$
• $h(x_0, z) = \frac{1}{2\pi (\sqrt{z}-1)}$

Let $\delta = z^{-1/4}$ and define the neighborhood $N_\delta = (x_0(1-\delta), x_0(1+\delta)) = (\sqrt{z}-z^{1/4}, \sqrt{z}+z^{1/4})$. For large $z$, $N_\delta \subset (0, z)$. We split $J(z)$ into integrals over $N_\delta$ and $(0, z) \setminus N_\delta$.

• Contribution from $N_\delta$: Let $x = x_0 + u = \sqrt{z} + u$, where $|u| \le z^{1/4}$.

  • Taylor expansion of $g(x,z)$ around $x_0$: $g(x,z) = g(x_0, z) + \frac{1}{2} g''(x_0, z) u^2 + O(g'''(x_0) u^3)$ $g(x,z) = 2(\sqrt{z}-1)^2 + \frac{1}{2} (8 - 4/\sqrt{z}) u^2 + O(z^{-3/2} u^3)$
  • Taylor expansion of $h(x,z)$ around $x_0$: $h(x,z) = h(x_0, z) (1 + \frac{u}{\sqrt{z}-1} + O(u^2/z)) = h(x_0, z) (1 + O(u/\sqrt{z}))$
  • The integral over $N_\delta$ of the main term ($R(x,z)=0$): $$ J_{N_\delta} = \int_{-z^{1/4}}^{z^{1/4}} h(x_0+u, z) e^{-\frac{1}{2} g(x_0+u, z)} du $$ $$ J_{N_\delta} \approx h(x_0, z) e^{-\frac{1}{2} g(x_0, z)} \int_{-z^{1/4}}^{z^{1/4}} (1 + O(u/\sqrt{z})) e^{-\frac{1}{4} g''(x_0, z) u^2} e^{O(z^{-3/2} u^3)} du $$ The $O(u/\sqrt{z})$ term integrates to zero by symmetry. The $e^{O(z^{-3/2} u^3)}$ term is $1+O(z^{-3/2} u^3)$; the integral of the $O(z^{-3/2} u^3)$ part is negligible. We extend the integration limits to $\pm \infty$ with exponentially small error: $$ J_{N_\delta} \approx h(x_0, z) e^{-\frac{1}{2} g(x_0, z)} \int_{-\infty}^{\infty} e^{-\frac{1}{4} g''(x_0, z) u^2} du $$ $$ J_{N_\delta} \approx h(x_0, z) e^{-\frac{1}{2} g(x_0, z)} \sqrt{\frac{4\pi}{g''(x_0, z)}} $$ $$ J_{N_\delta} \approx \frac{1}{2\pi (\sqrt{z}-1)} e^{-(\sqrt{z}-1)^2} \sqrt{\frac{4\pi}{8 - 4/\sqrt{z}}} $$ $$ J_{N_\delta} \sim \frac{1}{2\pi \sqrt{z}} e^{-(z - 2\sqrt{z} + 1)} \sqrt{\frac{\pi}{2}} = \frac{e^{-1}}{2 \sqrt{2\pi}} z^{-1/2} e^{-z + 2\sqrt{z}} $$

• Contribution outside $N_\delta$: For $x \in (0, z) \setminus N_\delta$, $|x - x_0| > z^{1/4}$. $g(x,z) \ge g(x_0, z) + C (x-x_0)^2$ for some $C > 0$ (related to $g''$). $g(x,z) > g(x_0, z) + C (z^{1/4})^2 = g(x_0, z) + C \sqrt{z}$. The integral over $(0, z) \setminus N_\delta$ is bounded by: $$ \int_{(0, z) \setminus N_\delta} |h(x,z)| e^{-\frac{1}{2} g(x,z)} dx < \text{Poly}(z) \times e^{-\frac{1}{2} g(x_0, z)} e^{-\frac{C}{2} \sqrt{z}} $$ This is exponentially smaller than $J_{N_\delta}$ (which is $\propto z^{-1/2} e^{-g(x_0)/2}$) because of the extra $e^{-C\sqrt{z}/2}$ factor.

• Contribution from Remainder $R(x,z)$: We need to estimate $\int_0^z h(x,z) e^{-\frac{1}{2} g(x,z)} R(x,z) dx$. In the dominant region $N_\delta$, $x \approx \sqrt{z}$, so $t(x) = z/x - 1 \approx \sqrt{z} - 1$. $R(x,z) = O(t(x)^{-2}) = O((\sqrt{z}-1)^{-2}) = O(z^{-1})$. The integral involving $R(x,z)$ over $N_\delta$ is approximately $J_{N_\delta} \times O(z^{-1})$. The integral outside $N_\delta$ is even smaller. Thus, the contribution from the remainder term is asymptotically negligible compared to $J_{N_\delta}$.

Combining these parts, we conclude: $$ I_1(z) \sim J_{N_\delta} \sim \frac{e^{-1}}{2 \sqrt{2\pi}} z^{-1/2} e^{-z + 2\sqrt{z}} $$

4. Analysis of $I_2(z)$ (Negligible Contribution)

This integral represents the contribution from $X < 0$. $$ I_2(z) = \int_{-\infty}^0 f(x) \Phi(z/x - 1) dx $$ For $x < 0$, as $z \to \infty$, the argument $t = z/x - 1 \to -\infty$. We use the asymptotic expansion for $\Phi(t)$ for $t \to -\infty$: $$ \Phi(t) = \frac{\phi(t)}{|t|} (1 + O(t^{-2})) = \frac{1}{\sqrt{2\pi} |t|} e^{-t^2/2} (1 + O(t^{-2})) $$ Substituting this gives an integral involving $e^{-\frac{1}{2} g(x,z)}$ where $g(x,z) = (x-1)^2 + (z/x - 1)^2$ is the same function as before. $$ I_2(z) \approx \int_{-\infty}^0 \frac{1}{2\pi (1 - z/x)} e^{-\frac{1}{2} g(x,z)} dx $$ The minimum of $g(x,z)$ for $x < 0$ occurs at $x_1 = -\sqrt{z}$. The value of $g$ at this minimum is $g(x_1, z) = 2(\sqrt{z}+1)^2 = 2z + 4\sqrt{z} + 2$. A full Laplace analysis is not needed here. We only need the exponential order of $I_2(z)$, which is determined by the minimum value of the exponent: $$ I_2(z) \propto e^{-\frac{1}{2} g(x_1, z)} = e^{-(\sqrt{z}+1)^2} = e^{-(z + 2\sqrt{z} + 1)} $$ (The notation $\propto$ means "is roughly of the order of").

5.

We compare the asymptotic behavior of $I_1(z)$ and $I_2(z)$:

• $I_1(z) \sim \frac{e^{-1}}{2 \sqrt{2\pi}} z^{-1/2} e^{-z + 2\sqrt{z}}$
• $I_2(z) \propto e^{-z - 2\sqrt{z} - 1}$

The ratio $I_2(z) / I_1(z)$ behaves like: $$ \frac{e^{-z - 2\sqrt{z} - 1}}{z^{-1/2} e^{-z + 2\sqrt{z} - 1}} = z^{1/2} e^{-4\sqrt{z}} $$ This ratio tends to 0 exponentially fast as $z \to \infty$. Therefore, $I_2(z)$ is negligible compared to $I_1(z)$, and $P(XY > z) = I_1(z) + I_2(z) \sim I_1(z)$.

The tail probability $P(XY > z)$ for independent $X, Y \sim N(1,1)$ is asymptotically dominated by the contribution from the region where $X > 0$. The application of Laplace's method to the integral $I_1(z)$ yields the final result: $$ P(XY > z) \sim \frac{e^{-1}}{2 \sqrt{2\pi}} z^{-1/2} e^{-z + 2\sqrt{z}} \quad \text{as } z \to \infty. $$

Answer 3

From the general result given in Theorem 2.3 of the following 2025 paper for the tail probability of the product of two random variables with a bivariate normal distribution:

Asymptotic approximations for the distribution of the product of correlated normal random variables

one can get as $z\to \infty$:

$$\mathrm{P}(XY>z) = \frac{\frac{1}{e}}{ \sqrt{2\pi z}} e^{-z} \cosh \left (2\sqrt{z} \right )\left (1+O(z^{-1/2}) \right) \\ =\frac{1}{2e \sqrt{2\pi }} z^{-1/2} \left (e^{-z+2\sqrt{z}}+ e^{-z-2\sqrt{z}} \right )\left (1+O(z^{-1/2}) \right)$$

by setting $\mu_X=\mu_Y=\sigma_X=\sigma_Y=1$ and $\rho=0$ in (2.7).

Answer 4

From $$\Pr (X Y<z)=\Pr \left(X<\frac{z}{Y}\right)=\int_{-\infty }^{\infty } \text{dPr}(Y<y) \Pr \left(X<\frac{z}{y}\right) $$

and

$$\text{dPr}(X Y<z)=\text{dPr}\left(X<\frac{z}{Y}\right)=\int_{-\infty \ }^{\infty } \text{dPr}(Y<y) \quad \text{dPr}\left(X<\frac{z}{y}\right) \, \ =\frac{\text{dz}}{2 \pi \sigma \tau } \ \left(\int_0^{\infty } \ \frac{e^{-\frac{y^2}{2 \sigma ^2}} e^{-\frac{\left(\frac{z}{y}\right)^2}{2 \tau ^2}}}{y^2} \, dy+\ldots \right) $$

$$\left( \int_0^{\infty } e^{-\frac{y^2}{2 \sigma ^2}} \ e^{-\frac{\left(\frac{z}{y}\right)^2} {2 \tau ^2}} \frac{1}{y^2} \, dy\right)^2 \ = \ \frac{\pi \tau ^2 e^{-\frac{2 z}{\sigma \tau }}}{2 z^2} $$

and similar with inverted signs for the negative axis.