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.Comment: 12 pages, 1 figur

Meyer, Christian

English

arXiv

Christian Meyer

Crossref

Journal of Statistical Software

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 delivers competitive performance and can easily be extended to arbitrary precision

Journal of Statistical Software 5and%x = −αx√1 + α2x, αx =1√1− %2(yx− %),%y = −αy√1 + α2y, αy =1√1− %2(xy− %).From the axis y = 0 to the diagonal x = y we get by applying the formula Meyer (2013,Equation 3.18)Φ2(x, 0; %) ={12 · Φ2(x, x; 1− 2%2), % < 0,Φ(x)− 12 · Φ2(x, x; 1− 2%2), % ≥ 0.(12)Specifically, we obtain1− 2%2x = 1−2 · (%x− y)2x2 + y2 − 2%xy= 1− 2 · ax1 + ax(13)withax =(%x− yx ·√1− %2)2(14)=(x− yx ·√1− %2−√1− %1 + %)2(15)=(x+ yx ·√1− %2−√1 + %1− %)2(16)where, in an implementation, Equation 15 should be used for % −→ 1, and Equation 16 for% −→ −1, in order to avoid catastrophic cancellation. Note also that%x < 0 ⇐⇒ αx > 0 ⇐⇒yx> %.In a last step, if necessary to ensure x ≤ 0 and % ≥ 0, we apply the formulas Meyer (2013,Equations 2.15 and 3.27)Φ2(x, x; %) = 2 · Φ(x)− 1 + Φ2(−x,−x; %), (17)Φ2(x, x; %) = 2 · Φ(x) · Φ (λ(%) · x)− Φ2 (λ(%) · x, λ(%) · x,−%) . (18)Specifically, we obtain|λ(1− 2%2x) · x| = |x| ·√%2x1− %2x=|%x− y|√1− %2(19)=∣∣∣∣∣ x− y√1− %2 − x ·√1− %1 + %∣∣∣∣∣ (20)=∣∣∣∣∣ x+ y√1− %2 − x ·√1 + %1− %∣∣∣∣∣ , (21)6 Recursive Numerical Evaluation of the Cumulative Bivariate Normal Distributionwhere, in an implementation, Equation 20 should be used for % −→ 1, and Equation 21 for% −→ −1, in order to avoid catastrophic cancellation.It will be favorable to work with2%2x =2 · ax1 + axinstead of 1− 2%2x. If Equation 18 has to be applied (i.e., if 1− 2%2x < 0, which is equivalentwith ax > 1), correspondingly we will work with1− (−(1− 2%2x)) = 2− 2%2x =21 + ax. (22)3. ImplementationIn the following we will discuss an implementation of the algorithm derived in Section 2. TheC++ language has been chosen because it is the market standard in quantitative finance, oneof the fields frequently requiring evaluation of normal distributions.3.1. Evaluation on the diagonalSource code (in C++) for evaluation of Φ2(x; %) as in Section 2.1, for x ≤ 0 and % ≥ 0, willbe provided at the end of this section. In the following we will comment on some details ofthe implementation.Equations 2–10 show that it is reasonable to provide a = 1− %, instead of %, as input for theevaluation of Φ2(x; %). Moreover, cf. Equations 13 and 22, double inversion (i.e., computationof 1− (1− z) instead of z) is to be avoided in the reduction algorithm.Values for px = Φ(x) and for pxs = Φ(λ(%) · x) are also expected as input parameters. Thismakes sense because the values are needed by the reduction algorithm as well (and henceshould not be computed twice).Evaluation of arcsin(%) is to be avoided for % −→ 1 and has been replaced (without optimiza-tion of the cutoff point) by arcsin(%) = arccos(√1− %2).Constants (c1 = 2/π, c2 =√π/2, c3 = π/2, c4 = 1/(2π)) have been pre-computed indouble precision. The recursion stops if a new term does not change the computed sum. Ifthe a priori bound for the absolute error, given by Equation 1, is less than 5 · 10−17, theupper bound is returned (relative accuracy on the diagonal may be increased by droppingthis condition but overall relative accuracy will still be determined by the reduction to thediagonal, cf. Section 3.2), and by the accuracy of the implementation of Φ. The final resultis always checked against the upper and lower bound.In order to avoid cancellation, d2k and d2k+1 will be bracketed before summation (cf. Sec-tion 4.2).double Phi2diag( double x, double a, double px, double pxs ){if( a <= 0.0 ) return px;if( a >= 1.0 ) return px * px;Journal of Statistical Software 7double b = 2.0 - a, sqrt_ab = sqrt( a * b );double c1 = 6.36619772367581343e-001;double c2 = 1.25331413731550025;double c3 = 1.57079632679489662;double c4 = 1.591549430918953358e-001;double asr = ( a > 0.1 ? asin( 1.0 - a ) : acos( sqrt_ab ) );double comp = px * pxs;if( comp * ( 1.0 - a - c1 * asr ) < 5e-17 )return b * comp;double tmp = c2 * x;double alpha = a * x * x / b;double a_even = -tmp * a;double a_odd = -sqrt_ab * alpha;double beta = x * x;double b_even = tmp * sqrt_ab;double b_odd = sqrt_ab * beta;double delta = 2.0 * x * x / b;double d_even = ( 1.0 - a ) * c3 - asr;double d_odd = tmp * ( sqrt_ab - a );double res = 0.0, res_new = d_even + d_odd;int k = 2;while( res != res_new ){d_even = ( a_odd + b_odd + delta * d_even ) / k;a_even *= alpha / k;b_even *= beta / k;k++;a_odd *= alpha / k;b_odd *= beta / k;d_odd = ( a_even + b_even + delta * d_odd ) / k;k++;res = res_new;res_new += d_even + d_odd;}res *= exp( -x * x / b ) * c4;return max( ( 1.0 + c1 * asr ) * comp, b * comp - max( 0.0, res ) );}3.2. Reduction to the diagonalThe following C++ function Phi2() evaluates Φ2(x, y; %) using Equation 11. It is assumedthat Phi(x) evaluates Φ(x). The special cases |%| = 1 and x = y = 0 are dealt with separately.8 Recursive Numerical Evaluation of the Cumulative Bivariate Normal Distributiondouble Phi2( double x, double y, double rho ){if( ( 1.0 - rho ) * ( 1.0 + rho ) <= 0.0 )if( rho > 0.0 )return Phi( min( x, y ) );elsereturn max( 0.0, min( 1.0, Phi( x ) + Phi( y ) - 1.0 ) );if( x == 0.0 && y == 0.0 )if( rho > 0.0 )return Phi2diag( 0.0, 1.0 - rho, 0.5, 0.5 );elsereturn 0.5 - Phi2diag( 0.0, 1.0 + rho, 0.5, 0.5 );return max( 0.0,min( 1.0,Phi2help( x, y, rho ) + Phi2help( y, x, rho ) ) );}The auxiliary function Phi2help() evaluates Φ2(x, 0; %x)− δx. Source code is provided in thefollowing. The special cases |%| = 1 and x = y = 0 are dealt with in Phi2(). Therefore, inPhi2help() there is no check against 1.0 - rho == 0.0, 1.0 + rho == 0.0 or s == 0.0.It is assumed that Phi(x) evaluates Φ(x), and that sqr(x) evaluates x*x. The cutoff points|%| = 0.99 in Phi2help() have been set by visual inspection and might be optimized.double Phi2help( double x, double y, double rho ){double s = sqrt( ( 1.0 - rho ) * ( 1.0 + rho ) );double a = 0.0, b1 = -fabs( x ), b2 = 0.0;if( rho > 0.99 ){double tmp = sqrt( ( 1.0 - rho ) / ( 1.0 + rho ) );b2 = -fabs( ( x - y ) / s - x * tmp );a = sqr( ( x - y ) / x / s - tmp );}else if( rho < -0.99 ){double tmp = sqrt( ( 1.0 + rho ) / ( 1.0 - rho ) );b2 = -fabs( ( x + y ) / s - x * tmp );a = sqr( ( x + y ) / x / s - tmp );}else{b2 = -fabs( rho * x - y ) / s;a = sqr( b2 / x );}Journal of Statistical Software 9double p1 = Phi( b1 ), p2 = Phi( b2 ), q = 0.0;if( a <= 1.0 )q = 0.5 * Phi2diag( b1, 2.0 * a / ( 1.0 + a ), p1, p2 );elseq = p1 * p2 - 0.5 * Phi2diag( b2, 2.0 / ( 1.0 + a ), p2, p1 );int c1 = ( y / x >= rho ), c2 = ( x < 0.0 ), c3 = c2 && ( y >= 0.0 );return ( c1 && c3 ? q - 0.5: c1 && c2 ? q: c1 ? 0.5 - p1 + q: c3 ? p1 - q - 0.5: c2 ? p1 - q: 0.5 - q );}4. DiscussionWe will now discuss convergence of the algorithm in theory (Section 4.1), the possibility ofcancellation (Section 4.2), and performance of the algorithm in practice (Section 4.3).4.1. ConvergenceIn this section we will discuss the convergence of the recursion for evaluation of the functionD (cf. Section 2.1). We will show that convergence is at least linear, and we will derive botha priori and a posteriori error bounds.Throughout the section we will assume that x ≤ 0, % ≥ 0. We start by noting thatfor k ≥ 0 even : ak+1 ≤ 0, bk+1 ≥ 0, ak+1 + bk+1 ≥ 0, dk ≥ 0;for k ≥ −1 odd : ak+1 ≥ 0, bk+1 ≤ 0, ak+1 + bk+1 ≤ 0, dk ≤ 0.Proofs are easy by induction. In particular, observe thata1 + b1 = x2 ·√1− %2 · 2%1 + %≥ 0and that d0 ≥ 0 due to the convexity of arcsin on [0, 1]. Moreover, for k ≥ 0 even, we haveak+3 + bk+3 =α · ak+1 + β · bk+1k + 3≥ α · (ak+1 + bk+1)k + 3≥ 0.Analogously, we finda0 + b0 =(√1− % ·(√1 + %−√1− %))·√π2· x ≤ 0and, for k ≥ −1 odd,ak+3 + bk+3 =α · ak+1 + β · bk+1k + 3≤ β · (ak+1 + bk+1)k + 3≤ 0.10 Recursive Numerical Evaluation of the Cumulative Bivariate Normal DistributionNow, for k ≥ 0 even, we conclude that1 + δk + 2· (dk + bk+1) ≥bk+1 + δ · dkk + 2+δ · bk+1k + 2≥ ak+1 + bk+1 + δ · dkk + 2+δβ· k + 3k + 2· bk+3 ≥ dk+2 + bk+3and, by induction,∑i≥2m,i evendi ≤ (d2m + b2m+1) ·(1 +1 + δ2m+ 2+(1 + δ)2(2m+ 2)(2m+ 4)+ . . .)= (d2m + b2m+1) ·∞∑k=0m! · (1 + δ)k2k · (m+ k)!≤ (d2m + b2m+1) ·∞∑k=0(1+δ2)kk!= (d2m + b2m+1) · exp(1 + δ2).Analogously, for k ≥ −1 odd, we find that1 + δk + 2· (dk + bk+1) ≤ dk+2 + bk+3and, by induction,∑i≥2m,i odddi ≥ (d2m+1 + b2m+2) ·(1 +1 + δ2m+ 3+(1 + δ)2(2m+ 3)(2m+ 5)+ . . .)≥ (d2m+1 + b2m+2) ·(1 +1 + δ2m+ 2+(1 + δ)2(2m+ 2)(2m+ 4)+ . . .)≥ (d2m+1 + b2m+2) · exp(1 + δ2).Combining the above inequalities, we find the following a posteriori error bound:∣∣∣∣∣∣∑i≥2mdi∣∣∣∣∣∣ ≤ max ∑i≥2m,i even|di|,∑i≥2m,i odd|di|≤ max (|d2m + b2m+1|, |d2m+1 + b2m+2|) · exp(1 + δ2)Applying the inequalities again, an a priori error bound can be derived as well:∣∣∣∣∣∣∑i≥2mdi∣∣∣∣∣∣ ≤ max (|d0 + b1|, |d1 + b2|) ·(1+δ2)mm!· exp(1 + δ2)Consequently, the order of convergence of the recursion is at least that of an exponential series.That is, convergence is at least linear. The above error bounds, of course, are rather weak.In particular, they are ignoring the alternating signs of the dk. We will therefore complementthe analysis by investigating the possibility of cancellation (Section 4.2), and the performanceof the algorithm in practice (Section 4.3).Journal of Statistical Software 114.2. CancellationSince the signs of the dk are alternating, cancellation may become an issue. However, it canbe avoided by choosing an appropriate order of summation. Numerical evidence leads to thefollowing conjecture:Conjecture: The sequence (dk)k≥−1 is log-concave, i.e., for all k ≥ −1 the following Turán-type inequality holds:dk+2 · dk ≤ d2k+1 (23)A proof of this conjecture is outside the scope of this paper. Note that related statementshave been proved for Mills’ ratio, cf. Baricz (2008).If the conjecture holds and if there is m ≥ 0 with sm := d2m + d2m+1 ≥ 0 then−1 ≥ d2md2m+1≥ d2m+2d2m+3and therefore sm+1 = d2m+2 + d2m+3 ≥ 0 as well. That is, the sequence (sm)m≥0 changessigns not more than once. More precisely, we haves0 ≤ 0 ⇐⇒ x ≥ x(%) :=√2π· arcsin(%)− % · π/2%− 1 +√1− %2.The function x(%) is decreasing in % withlim%−→0x(%) =√2π−√π2≈ −0.455, lim%−→1x(%) = −√2π≈ −0.798.In order to avoid cancellation, it is therefore reasonable to group the dk into pairs.4.3. PerformanceEvaluation of Φ2(x, y; %) as in Section 3 will require (at most) four calls to an implementationof the cumulative standard normal distribution Φ (Phi() in the code). The actual choicemay well determine both precision and running time of the algorithm. For testing purposesI have been using a hybrid method, calling the algorithm from West (2005, Figur 2) (a C++implementation is available from West (2006)) for absolute value of the argument larger than0.5, and Phi() from Marsaglia (2004) else. Besides Phi(), exp() will be called two times,asin() or acos() two times, and sqrt() six times. Everything else is elementary arithmetic.Due to the reduction algorithm, the final result will be a sum. Therefore, very high precisionin terms of relative error can not be expected. Consequently, evaluation of the diagonal aimsat absolute error as well.The Phi2diag() function is behaving as it may be expected from an approximation by aTaylor series around zero: (absolute) error increases with decreasing x. For x < −7 (or% −→ 0 or % −→ 1) the error bounds from Equation 1 are taking over, and absolute errordecreases again. The maximum absolute error is obtained for x ≈ −7, % ≈ 0.8 (maximum errorof the upper bound is obtained for % =√1− 4/π2 ≈ 0.7712, cf. Meyer (2013, Theorem 5.2)).In general, assuming that all numerical fallacies in the reduction algorithm have been takencare of, the diagonal is expected to provide a worst case because the errors of the two calls12 Recursive Numerical Evaluation of the Cumulative Bivariate Normal Distribution0.E+001.E-162.E-163.E-164.E-165.E-166.E-167.E-16-10 -8 -6 -4 -2 0 2 4 6 8 10Phi2 (max)Genz (max)Phi2 (99%)Genz (99%)Figure 1: Absolute error of implementations of Φ2 in a simulation study.to Phi2diag() will not cancel. With respect to the reduction algorithm, the case %x ≈ y,%y ≈ x, implying |%| ≈ 1, is most critical.In order to give an impression of the algorithm’s behaviour, we will discuss the results of asimulation study. For each n ∈ {0, . . . , 200}, m ∈ {1, . . . , 106}, the value of Φ2(xnm, ynm; %nm)has been computed via the Phi2() function from Section 3.2 where xnm has been drawn froma uniform distribution on [xn − 0.05, xn + 0.05] with xn := n/10 − 10, ynm has been drawnfrom a uniform distribution on [−10, 10], and %nm := 2Φ(rnm)− 0.5 where rnm has been drawnfrom a uniform distribution on [−10, 10] as well.The C++ implementation available from West (2006) of the Algorithm of Genz (2004) hasbeen serving as a competitor. Both functions have been evaluated against a quad-doubleprecision version of Phi2(), implemented using the QD library of Bailey, Hida, and Li (2010)and quad-double precision constants.The diagram in Figure 1 is displaying, for n ∈ {0, . . . , 200}, the 99% quantile and the maxi-mum of the absolute difference between the double precision algorithms (Phi2 and Genz) andthe quad-double precision algorithm.Apart from a shift due to subtractions for positive x, errors of Phi2 are rather symmetricaround zero. The peaks at |xn| ≈ 7 are due to the Taylor expansion around zero; the peaks at|xn| ≈ 2 are due to Taylor expansion after transformation of the argument. The characteristicsof the 99% quantile, in particular the little peaks at |xn| ≈ 0.7, are already visible in the errorof the Φ function used. The maximum error of Genz almost always stays below the one ofJournal of Statistical Software 13Phi2. Note that the maximum error of Genz is determined by the case % −→ −1 and mightbe reduced by careful consideration of that case.In the simulation study, Phi2 was a little slower than Genz: it took approximately five min-utes and four minutes to perform the 201 · 106 evaluations on a fairly standard office PC(and it took two days to perform the corresponding quad-double precision evaluations). Thenumber of recursion steps used by Phi2diag is increasing with |x|. It should be easy to findan appropriate trade-off between speed and precision (maybe even equal precision with lessrunning time) by replacing the condition terminating the recursion.AcknowledgmentsThe author wishes to thank Axel Vogt for helpful discussion.ReferencesAǧca Ş, Chance DM (2003). “Speed and Accuracy Comparison of Bivariate Normal Dis-tribution Approximations for Option Pricing.” Journal of Computational Finance, 6(4),61–96.Bailey DH, Hida Y, Li XS (2010). “QD (C++/Fortran-90 Double-Double and Quad-DoublePackage).” URL http://crd.lbl.gov/~dhbailey/mpdist/.Balakrishnan N, Lai CD (2009). Continuous Bivariate Distributions. 2nd edition. Springer-Verlag, New York.Baricz Á (2008). “Mills’ Ratio: Monotonicity Patterns and Functional Inequalities.” Journalof Mathematical Analysis and Applications, 340, 1362–1370.Divgi DR (1979). “Calculation of Univariate and Bivariate Normal Probability Functions.”The Annals of Statistics, 7(4), 903–910.Drezner Z (1978). “Computation of the Bivariate Normal Integral.” Mathematics of Compu-tation, 32(141), 277–279.Drezner Z, Wesolowsky GO (1990). “On the Computation of the Bivariate Normal Integral.”Journal of Statistical Computation and Simulation, 35, 101–107.Genz A (2004). “Numerical Computation of Rectangular Bivariate and Trivariate Normaland t Probabilities.” Statistics and Computing, 14(3), 151–160.Genz A, Bretz F (2009). Computation of Multivariate Normal and t Probabilities. Springer-Verlag, New York.Hull J (2008). Options, Futures, and Other Derivatives. 6th edition. Prentice Hall.Kotz S, Balakrishnan N, Johnson, Lloyd N (2000). Continuous Multivariate Distributions,volume 1: Models and Applications. 2nd edition. John Wiley & Sons.14 Recursive Numerical Evaluation of the Cumulative Bivariate Normal DistributionMarsaglia G (2004). “Evaluating the Normal Distribution.” Journal of Statistical Software,11(4), 1–11. URL http://www.jstatsoft.org/v11/i04/.Meyer C (2013). “The Bivariate Normal Copula.” Communications in Statistics – Theory andMethods. Forthcoming, preprint available at http://arxiv.org/abs/0912.2816.Patel JK, Read CB (1996). Handbook of the Normal Distribution. 2nd edition. Dekker.Terza JV, Welland U (1991). “A Comparison of Bivariate Normal Algorithms.” Journal ofStatistical Computation and Simulation, 39(1-2), 115–127.Vogt A (2008). “Computing the Cumulative Bivariate Normal Distribution.” URL http://www.axelvogt.de/axalom/.Wang MC, Kennedy WJ (1990). “Comparison of Algorithms for Bivariate Normal Probabil-ity over a Rectangle Based on Self-Validated Results from Interval Analysis.” Journal ofStatistical Computation and Simulation, 37(1-2), 13–25.West G (2005). “Better Approximations to Cumulative Normal Functions.” Wilmott Maga-zine, pp. 70–76.West G (2006). “Code for Univariate, Bivariate and Trivariate Cumulative Normal Functions.”URL http://www.finmod.co.za/research.htm.Affiliation:Christian MeyerDZ BANK AGPlatz der Republik60265 Frankfurt, GermanyE-mail: christian.meyer@dzbank.deJournal of Statistical Software http://www.jstatsoft.org/published by the American Statistical Association http://www.amstat.org/Volume 52, Issue 10 Submitted: 2010-04-28March 2013 Accepted: 2012-12-18

Directory of Open Access Journals

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 mathemat-ically transparent, delivers competitive performance and can easily be extended to arbitrary precision. 

Recursive Numerical Evaluation of the Cumulative Bivariate Normal Distribution

Abstract

Similar works

Full text

Available Versions

Crossref

Journal of Statistical Software

Directory of Open Access Journals

CiteSeerX