The Transformed BVN Problem

If the substitution $r=\sin{\theta}$ is used, then

$$L(h,k,\rho) = \Phi(-h)\Phi(-k) + \frac{1}{2\pi}\int_0^{\sin^... ... e^{-\frac{h^2+k^2-2hk\sin(\theta)}{2\cos(\theta)^2}}d\theta ,$$ (3)

but there is a singularity in the integrand when $r=1$ with both formulas, and the Drezner and Wesolowsky fixed-rule numerical integration method loses accuracy when $\vert\rho\vert \simeq 1$. In order to avoid this singularity, Drezner and Wesolowsky integrated between $\rho$ and $s = sign(\rho)$ rather than between $0$ and $\rho$ to obtain

$$L(h,k,\rho) = L(h,k,s) - \frac{s}{2\pi}\int_\rho^s \frac{1} {\sqrt{1-r^2}}e^{-\frac{h^2+k^2-2hkr}{2(1-r^2)}}dr,$$

where

$$L(h,k,s) = \left\{ \begin{array}{cl} \Phi(-max(h,k)) & \mb... ...(0,\Phi(-h)-\Phi(k)) & \mbox{if s = -1} \end{array} \right. .$$

If the substitution $x = \sqrt{1-r^2}$ is used, then

$$L(h,k,\rho) = L(h,k,s) - \frac{s}{2\pi} \int_0^{\sqrt{1-\rho... ...1}{\sqrt{1-x^2}} e^{-\frac{h^2+k^2-2shk\sqrt{1-x^2}}{2x^2}}dx.$$ (4)