mjhjmtu

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

$$int_0^infty int_0^infty (1 + alpha u)^-1(1 + v)^-1 Phileft(fracbetasqrtgamma + uvright) textdu textdv,qquad gamma > 0, quad beta in mathbbR, quad alpha in mathbbN^*,$$
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:

$$textEleft Phileft(fracbetasqrtgamma + X^2Y^2right)right,qquad gamma > 0, quad beta in mathbbR,$$
where
$$X sim textHalf-Cauchy(0,1), qquad Y sim textHalf-Cauchy(0,alpha^-1/2), quad alpha in mathbbN^* , (textlarge).$$

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 Phileft(fracbetasqrtgamma + x^2y^2right) textdx textdy,quad gamma > 0, beta in mathbbR, alpha in mathbbN^*, 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?)
$$textEleft Phileft(fracbetasqrtgamma + X^2Y^2right)right - Phileft(fracbetasqrtgammaright),qquad gamma > 0, quad beta in mathbbR,$$
would it be sensible to approximate the difference
$$Phileft(fracbetasqrtgamma + X^2Y^2right) - Phileft(fracbetasqrtgammaright),$$
with a Taylor series in zero or is this unreasonable given that the half-Cauchy distributions are heavy-tailed?

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

$$beta/sqrtgamma + uv ge min(beta/sqrtgamma, 0) = delta gt -infty.$$

Therefore, because $Phi$ is a CDF for a distribution supported on $(-infty,infty),$ $$Phileft(fracbetasqrtgamma + uvright)ge Phi(delta) = epsilon gt 0.$$

Consequently

$$eqalign
&int_0^infty int_0^infty (1+alpha u)^-1(1+v)^-1 Phileft(fracbetasqrtgamma + uvright) mathrmdumathrmdv \
&ge
epsilonint_0^infty int_0^infty (1+alpha u)^-1(1+v)^-1 mathrmdumathrmdv \
&= lim_Mtoinftylim_Ntoinftyepsilonint_0^M (1+alpha u)^-1mathrmduint_0^N (1+v)^-1 mathrmdv \
&=fracepsilonalpha lim_Mtoinftylim_Ntoinfty log(1 + Malpha)log(1 + N),
$$

which diverges to $+infty.$

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)))