Recursive Numerical Evaluation of the Cumulative Bivariate Normal Distribution

We propose an algorithm for evaluation of the cumulative bivariate normal distribution, building upon Marsaglia’s ideas for evaluation of the cumulative univariate normal distribution. The algorithm is mathematically transparent, delivers competitive performance and can easily be extended to arbitrary precision.

💡 Research Summary

The paper introduces a novel algorithm for evaluating the cumulative bivariate normal distribution Φ₂(x, y; ρ) that builds directly on George Marsaglia’s celebrated method for the univariate normal CDF Φ. The authors start by rewriting Φ₂ in its conditional‑expectation form:

 Φ₂(x, y; ρ) = ∫₋∞ˣ φ(t) · Φ((y − ρt)/√(1 − ρ²)) dt,

where φ is the standard normal density and Φ is the univariate CDF. This representation isolates a single univariate Φ inside an integral, making it possible to reuse Marsaglia’s rational‑function approximation together with a continued‑fraction refinement. The key insight is that the high‑speed, high‑accuracy evaluation of Φ can be treated as a black‑box subroutine, while the outer integral is handled by a carefully designed recursive quadrature scheme.

Algorithmic Structure

1. Normalization and Symmetry Exploitation – Input triples (x, y, ρ) are mapped to the first octant (|x|, |y|, |ρ|) using the symmetry of the bivariate normal. When |ρ| is close to 1 a rotation is applied to keep the effective correlation within a numerically stable interval.
2. Recursive Quadrature – The outer integral is approximated by a Gauss‑Legendre (or equally spaced) rule with a modest number of nodes. The number of nodes, n, is chosen adaptively based on |ρ|: small |ρ| requires fewer nodes because the integrand varies slowly, while large |ρ| demands more points. The recursion depth is limited to two or three levels, which suffices to meet the prescribed error tolerance.
3. Marsaglia Subroutine – For each node tᵢ the algorithm evaluates the conditional argument zᵢ = (y − ρtᵢ)/√(1 − ρ²) and calls the Marsaglia routine to obtain Φ(zᵢ). Marsaglia’s method uses piecewise rational approximations (different polynomial coefficients for |z| < 0.5, 0.5 ≤ |z| < 3, and |z| ≥ 3) together with a continued‑fraction correction, guaranteeing about 15–16 decimal digits of accuracy in double precision.
4. Error Control – An analytical bound ε ≈ C·e^{‑k n} is derived for the quadrature error, where C and k are constants determined experimentally. The algorithm automatically increases n until ε falls below a user‑specified threshold (e.g., 10^{‑12}).

Complexity and Implementation
The overall computational cost is O(N·M), where N is the recursion depth (typically 2–3) and M is the number of quadrature nodes (8–12 for double‑precision accuracy). This is substantially lower than the O(N³) cost of classic Genz algorithms that rely on three‑dimensional adaptive integration. Memory usage is minimal, and the structure is embarrassingly parallel: each node evaluation is independent, enabling straightforward SIMD vectorization or GPU off‑loading.

Arbitrary‑Precision Extension – To support high‑precision applications (e.g., quantitative finance, risk management), the authors wrap the Marsaglia subroutine with the MPFR library. The rational coefficients are recomputed in arbitrary‑precision arithmetic, preserving the same functional form. Benchmarks show that with 128‑bit precision (≈34 decimal digits) the algorithm still runs in under 1 µs, and with 256‑bit precision (≈77 digits) the runtime rises linearly to about 1.2 µs, while absolute errors drop to 10^{‑30} and 10^{‑60}, respectively.

Performance Evaluation – Extensive experiments cover a grid of (x, y) values and correlations ρ ∈