The Bivariate Normal Integral via Owen's T Function as a Modified Euler's Arctangent Series

The Owen's T function is presented in four new ways, one of them as a series similar to the Euler's arctangent series divided by $2\pi$, which is its majorant series. All possibilities enable numerically stable and fast convergent computation of the bivariate normal integral with simple recursion. When tested $\Phi_\varrho^2(x,y)$ computation on a random sample of one million parameter triplets with uniformly distributed components and using double precision arithmetic, the maximum absolute error was $3.45\cdot 10^{-16}$. In additional testing, focusing on cases with correlation coefficients close to one in absolute value, when the computation may be very sensitive to small rounding errors, the accuracy was retained. In rare potentially critical cases, a simple adjustment to the computation procedure was performed - one potentially critical computation was replaced with two equivalent non-critical ones. All new series are suitable for vector and high-precision computation, assuming they are supplemented with appropriate efficient and accurate computation of the arctangent and standard normal cumulative distribution functions. They are implemented by the R package Phi2rho, available on CRAN. Its functions allow vector arguments and are ready to work with the Rmpfr package, which enables the use of arbitrary precision instead of double precision numbers. A special test with up to 1024-bit precision computation is also presented.