Solving a double integral over transformations of joint bivariate standard normal values

My problem is about calculating the covariance between transformations of two test statistics based on the correlation between these test statistics. Let $$X$$ and $$Y$$ be two test statistics whose joint distribution is a standard bivariate normal distribution with a known correlation. Then, the joint distribution of these two variables can be formulated as follows $$f_{X, Y}\left(x, y; \mu = \left[\begin{array}{c} 0 \\ 0 \end{array}\right], \Sigma = \left[\begin{array}{cc} 1 & \rho\\ \rho & 1\end{array}\right]\right) = \frac{1}{2\pi\sqrt{1 - \rho^2}} \exp\left(-\frac{1}{2(1 - \rho^2)}(x^2+y^2-2\rho xy)\right)$$ The transformation that I want to apply actually transforms $$X$$ and $$Y$$ back to a standard normal distribution in two steps. The first step of the transformation is to convert the test statistics into $$p$$-values, whereas the second step transforms them back to standard normal values. However, there are two types of this transformation. The first type converts the test statistics into one-sided $$p$$-values which are converted back to standard normal values in the second step. This transformation can be shown as follows $$X^* = \Phi^{-1}(1 - \Phi(X)) \\ Y^* = \Phi^{-1}(1 - \Phi(Y)),$$ where $$X^*$$ and $$Y^*$$ are transformations of $$X$$ and $$Y$$, respectively, and $$\Phi$$ is the distribution function of standard normal distribution and $$\Phi^{-1}$$ is its inverse. This transformation seems trivial for one-sided $$p$$-values; however, it gets a little more tricky with two-sided $$p$$-values. For two-sided $$p$$-values the transformation can be done as $$X^* = \Phi^{-1}(1 - 2\Phi(1 - |X|)) \\ Y^* = \Phi^{-1}(1 - 2\Phi(1 - |Y|)).$$
Although the transformation is different in two-sided $$p$$-values case, $$X^*$$ and $$Y^*$$ are still distributed with a standard bivariate normal distribution, but the covariance between $$X^*$$ and $$Y^*$$ is unknown. In a similar problem, Brown (1975) showed the calculation of covariance between transformations of two joint standard normal values. The same approach shown in Brown (1975) can be applied to my question by $$\mbox{Cov}(X^*, Y^*) = \mbox{E}[X^*Y^*] - \mbox{E}[X^*]\mbox{E}[Y^*].$$ Assuming that $$X^*$$ and $$Y^*$$ jointly follow bivariate normal distribution $$\mbox{E}[X^*]$$ and $$\mbox{E}[Y^*]$$ are equal to $$0$$. $$\mbox{E}[X^*Y^*]$$, however, can be calculated by solving the following double integral (for two-sided $$p$$-values) $$\int_{-l}^{l}\int_{-l}^{l}X^*Y^*f_{X, Y}\left(x, y; \mu = \left[\begin{array}{c} 0 \\ 0 \end{array}\right], \Sigma = \left[\begin{array}{cc} 1 & \rho\\ \rho & 1\end{array}\right]\right)dxdy \\ = \int_{-l}^{l}\int_{-l}^{l}\Phi^{-1}(1 - 2\Phi(1 - |x|)) \Phi^{-1}(1 - 2\Phi(1 - |y|)) \frac{1}{2\pi\sqrt{1 - \rho^2}} \exp\left(-\frac{1}{2(1 - \rho^2)}(x^2+y^2-2\rho xy)\right) dxdy,$$ where $$\Phi(x) = \frac{1}{2}[1 + erf(x/\sqrt{2})]$$ ($$erf(.)$$ is the error function) and $$l$$ is the scalar limit for the test statistics, $$X$$ and $$Y$$.

I am trying to solve this double integral with numerical integration method implemented in R with adaptIntegrate function in package cubature where I have a problem. The problem arises from the calculation of the transformed values for a specific value of $$X$$ (or $$Y$$). When $$X$$ (or $$Y$$) is 0, the two-sided $$p$$-value becomes 1 (you can try 2 * pnorm(abs(0), lower.tail = FALSE)). Then, since the $$p$$-value is 1, when I convert it back to a standard normal value (qnorm(2 * pnorm(abs(0), lower.tail = FALSE), lower.tail = FALSE)) results in -Inf and it does not allow me to calculate the integral above, because the result converges to infinity.

Do you have any ideas about how I can overcome this problem? Or, do you think the multiplication in the double integral has a closed form so that I can solve it exactly without numerical methods? I'd be glad to hear any help or ideas. Thanks in advance.

• In your double integrals are you using "$X$" to mean "$x$" and "$Y$" to mean "$y$"? If not, then you will need to indicate what they do mean so that we can make sense of the question. – whuber Dec 10 '18 at 16:00
• Yes, you are right. My mistake, I'll edit the post. Thanks. – Ozan C. Dec 12 '18 at 13:39

Given $$\left(\begin{matrix} X\\Y \end{matrix}\right) \sim N\left[(\left(\begin{matrix} 0\\0 \end{matrix}\right),\Sigma_{2\times 2}\right]$$ , then $$\Sigma_{2\times 2}^{-\frac 12}\left(\begin{matrix} X\\Y \end{matrix}\right) \sim N\left[(\left(\begin{matrix} 0\\0 \end{matrix}\right),\mathrm{I}_{2\times 2}\right]$$
So you just need to find $$\Sigma_{2\times 2}^{-\frac 12}$$.