Double integral involving the normal CDF

Question

I need to compute (or best approximate?) the following integral

$$\int_0^\infty \int_0^\infty (1 + \alpha u)^{-1}(1 + v)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right) \text{d}u \text{d}v,\qquad \gamma > 0, \quad \beta \in \mathbb{R}, \quad \alpha \in \mathbb{N}^*,$$ where $\Phi(\cdot)$ is the standard normal cumulative distribution function and $\alpha$ is typically a large natural number. What strategy would you advise for this?

I don't know if this is helpful but I can also reformulate my problem (up to multiplicative constants) as computing the following expectation:

$$\text{E}\left\{ \Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right)\right\},\qquad \gamma > 0, \quad \beta \in \mathbb{R},$$ where $$X \sim \text{Half-Cauchy}(0,1), \qquad Y \sim \text{Half-Cauchy}(0,\alpha^{-1/2}), \quad \alpha \in \mathbb{N}^* \, (\text{large}).$$

EDIT: I made a mistake in the first integral, it should be:

$$\int_0^\infty \int_0^\infty (1 + \alpha x^2)^{-1}(1 + y^2)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + x^2y^2}}\right) \text{d}x \text{d}y,\quad \gamma > 0, \beta \in \mathbb{R}, \alpha \in \mathbb{N}^*\, \text{(large)},$$ which now corresponds to the above expectation.

EDIT 2: Another way to solve my problem would be to characterise the difference (hopefully small?) $$\text{E}\left\{ \Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right)\right\} - \Phi\left(\frac{\beta}{\sqrt{\gamma}}\right),\qquad \gamma > 0, \quad \beta \in \mathbb{R},$$ would it be sensible to approximate the difference $$\Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right) - \Phi\left(\frac{\beta}{\sqrt{\gamma}}\right),$$ with a Taylor series in zero or is this unreasonable given that the half-Cauchy distributions are heavy-tailed?

Answer 1

Under the conditions that $\alpha, \gamma, u,$ and $v$ are all positive,

$$\beta/\sqrt{\gamma + uv} \ge \min(\beta/\sqrt{\gamma}, 0) = \delta \gt -\infty.$$

Therefore, because $\Phi$ is a CDF for a distribution supported on $(-\infty,\infty),$ $$\Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right)\ge \Phi(\delta) = \epsilon \gt 0.$$

Consequently

$$\eqalign{ &\int_0^\infty \int_0^\infty (1+\alpha u)^{-1}(1+v)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right) \mathrm{d}u\mathrm{d}v \\ &\ge \epsilon\int_0^\infty \int_0^\infty (1+\alpha u)^{-1}(1+v)^{-1} \mathrm{d}u\mathrm{d}v \\ &= \lim_{M\to\infty}\lim_{N\to\infty}\epsilon\int_0^M (1+\alpha u)^{-1}\mathrm{d}u\int_0^N (1+v)^{-1} \mathrm{d}v \\ &=\frac{\epsilon}{\alpha} \lim_{M\to\infty}\lim_{N\to\infty} \log(1 + M\alpha)\log(1 + N), }$$

which diverges to $+\infty.$

Answer 2

If approximation is OK, you could simulate the expected value through an average. A law of large numbers implies that the average converges to the expected value - provided the latter exists, which is however not obvious to me. I did simulate a few runs, though, and got very similar results each time, which is not indicative of an issue with heavy tails.

library(LaplacesDemon)

draws <- 1e6
alpha <- 3
gamma <- 1
beta <- 1
x <- rhalfcauchy(draws, scale=1)
y <- rhalfcauchy(draws, scale=1/sqrt(alpha))

mean(pnorm(beta/sqrt(gamma+x^2*y^2)))