Appendix: Fortran Source Code REAL FUNCTION BVN( SH, SK, R ) * * A function for computing bivariate normal probabilities. * This function uses an algorithm given in the paper * "Numerical Computation of Bivariate and * Trivariate Normal Probabilities", by * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenz@wsu.edu * * BVN - calculate the probability that X is larger than SH and Y is * larger than SK. * * Parameters * * SH REAL, integration limit * SK REAL, integration limit * R REAL, correlation coefficient * LG INTEGER, number of Gauss Rule Points and Weights * REAL SH, SK, R, ZERO, TWOPI INTEGER I, LG PARAMETER ( ZERO = 0, TWOPI = 6.283185, LG = 5 ) REAL X(LG), W(LG) , AS, A, B, C, RS, XS REAL PHI, SN, ASR, H, K, BS, HS, HK SAVE X, W * X(I) = +-SQRT((5+-SQRT(10/7))/9), 0 DATA X(1), X(2), X(3) / -0.9061798, -0.5384693, 0.0 / DATA W(1), W(2), W(3) / 0.2369269, 0.4786287, 0.5688889 / X(4) = -X(2) X(5) = -X(1) W(4) = W(2) W(5) = W(1) H = SH K = SK HK = H*K BVN = 0 * * Absolute value of the correlation coefficient is less than 0.8 * IF ( ABS(R) .LT. 0.8 ) THEN HS = ( H*H + K*K )/2 ASR = ASIN( R )/2 DO I = 1, LG SN = SIN( ASR*(X(I)+1) ) BVN = BVN + W(I)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) ) END DO BVN = ASR*BVN/TWOPI + PHI(-H)*PHI(-K) ELSE * * For larger correlation coefficient * IF ( R .LT. 0 ) THEN K = -K HK = -HK ENDIF IF ( ABS(R) .LT. 1 ) THEN AS = ( 1 - R )*( 1 + R ) A = SQRT( AS ) BS = ( H - K )**2 C = ( 4 - HK )/8 BVN = A*EXP( -(BS/AS + HK)/2 )*( 1 - C*(BS - AS)/3 ) IF ( -HK .LT. 100 ) THEN B = SQRT(BS) BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*PHI(-B/A)*B*(1-C*BS/3) ENDIF A = A/2 DO I = 1, LG XS = ( A*( X(I) + 1 ) )**2 RS = SQRT( 1 - XS ) ASR = -( BS/XS + HK )/2 IF ( ASR .GT. -100 ) BVN = BVN + A*W(I)* & ( EXP(-BS/(2*XS)-HK/(1+RS))/RS - EXP(ASR)*(1+C*XS) ) END DO BVN = -BVN/TWOPI ENDIF IF ( R .GT. 0 ) BVN = BVN + PHI( -MAX( H, K ) ) IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, PHI(-H) - PHI(-K) ) ENDIF END * REAL FUNCTION PHI(X) REAL X DOUBLE PRECISION PHID PHI = PHID( DBLE(X) ) END * DOUBLE PRECISION FUNCTION BVND( DH, DK, R ) * * A function for computing bivariate normal probabilities. * * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenz@wsu.edu * * This function is based on the method described by * Drezner, Z and G.O. Wesolowsky, (1989), * On the computation of the bivariate normal inegral, * Journal of Statist. Comput. Simul. 35, pp. 101-107, * with major modifications for double precision, and for |R| close to 1. * * BVND - calculate the probability that X is larger than DH and Y is * larger than DK. * * Parameters * * DH DOUBLE PRECISION, integration limit * DK DOUBLE PRECISION, integration limit * R DOUBLE PRECISION, correlation coefficient * DOUBLE PRECISION DH, DK, R, ZERO, TWOPI INTEGER I, IS, LG, NG PARAMETER ( ZERO = 0, TWOPI = 6.283185307179586D0 ) DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS, BVN DOUBLE PRECISION PHID, SN, ASR, H, K, BS, HS, HK * Gauss Legendre Points and Weights, N = 6 DATA ( W(I,1), X(I,1), I = 1,3) / & 0.1713244923791705D+00,-0.9324695142031522D+00, & 0.3607615730481384D+00,-0.6612093864662647D+00, & 0.4679139345726904D+00,-0.2386191860831970D+00/ * Gauss Legendre Points and Weights, N = 12 DATA ( W(I,2), X(I,2), I = 1,6) / & 0.4717533638651177D-01,-0.9815606342467191D+00, & 0.1069393259953183D+00,-0.9041172563704750D+00, & 0.1600783285433464D+00,-0.7699026741943050D+00, & 0.2031674267230659D+00,-0.5873179542866171D+00, & 0.2334925365383547D+00,-0.3678314989981802D+00, & 0.2491470458134029D+00,-0.1252334085114692D+00/ * Gauss Legendre Points and Weights, N = 20 DATA ( W(I,3), X(I,3), I = 1, 10 ) / & 0.1761400713915212D-01,-0.9931285991850949D+00, & 0.4060142980038694D-01,-0.9639719272779138D+00, & 0.6267204833410906D-01,-0.9122344282513259D+00, & 0.8327674157670475D-01,-0.8391169718222188D+00, & 0.1019301198172404D+00,-0.7463319064601508D+00, & 0.1181945319615184D+00,-0.6360536807265150D+00, & 0.1316886384491766D+00,-0.5108670019508271D+00, & 0.1420961093183821D+00,-0.3737060887154196D+00, & 0.1491729864726037D+00,-0.2277858511416451D+00, & 0.1527533871307259D+00,-0.7652652113349733D-01/ SAVE X, W IF ( ABS(R) .LT. 0.3D0 ) THEN NG = 1 LG = 3 ELSE IF ( ABS(R) .LT. 0.75D0 ) THEN NG = 2 LG = 6 ELSE NG = 3 LG = 10 ENDIF H = DH K = DK HK = H*K BVN = 0 IF ( ABS(R) .LT. 0.925D0 ) THEN HS = ( H*H + K*K )/2 ASR = ASIN(R) DO I = 1, LG DO IS = -1, 1, 2 SN = SIN( ASR*( IS*X(I,NG) + 1 )/2 ) BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) ) END DO END DO BVN = BVN*ASR/( 2*TWOPI ) + PHID(-H)*PHID(-K) ELSE IF ( R .LT. 0 ) THEN K = -K HK = -HK ENDIF IF ( ABS(R) .LT. 1 ) THEN AS = ( 1 - R )*( 1 + R ) A = SQRT(AS) BS = ( H - K )**2 C = ( 4 - HK )/8 D = ( 12 - HK )/16 ASR = -( BS/AS + HK )/2 IF ( ASR .GT. -100 ) BVN = A*EXP(ASR) & *( 1 - C*( BS - AS )*( 1 - D*BS/5 )/3 + C*D*AS*AS/5 ) IF ( -HK .LT. 100 ) THEN B = SQRT(BS) BVN = BVN - EXP( -HK/2 )*SQRT(TWOPI)*PHID(-B/A)*B & *( 1 - C*BS*( 1 - D*BS/5 )/3 ) ENDIF A = A/2 DO I = 1, LG DO IS = -1, 1, 2 XS = ( A*( IS*X(I,NG) + 1 ) )**2 RS = SQRT( 1 - XS ) ASR = -( BS/XS + HK )/2 IF ( ASR .GT. -100 ) THEN BVN = BVN + A*W(I,NG)*EXP( ASR ) & *( EXP( -HK*( 1 - RS )/( 2*( 1 + RS ) ) )/RS & - ( 1 + C*XS*( 1 + D*XS ) ) ) END IF END DO END DO BVN = -BVN/TWOPI ENDIF IF ( R .GT. 0 ) BVN = BVN + PHID( -MAX( H, K ) ) IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, PHID(-H) - PHID(-K) ) ENDIF BVND = BVN END * DOUBLE PRECISION FUNCTION TVND( LIMIT, SIGMA ) * * A function for computing trivariate normal probabilities. * This function uses an algorithm given in the paper * "Numerical Computation of Bivariate and * Trivariate Normal Probabilities", * by * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenz@wsu.edu * * TVND calculates the probability that X(I) < LIMIT(I), I = 1, 2, 3. * * Parameters * * LIMIT DOUBLE PRECISION array of three upper integration limits. * SIGMA DOUBLE PRECISION array of three correlation coefficients, * SIGMA should contain the lower left portion of the * correlation matrix R. * SIGMA(1) = R(2,1), SIGMA(2) = R(3,1), SIGMA(3) = R(3,2). * * TVND cuts the outer integral over -infinity to B1 to * an integral from -8.5 to B1 and then uses an adaptive * integration method to compute the integral of a bivariate * normal distribution function. * EXTERNAL TRVFND DOUBLE PRECISION LIMIT(3), SIGMA(3) LOGICAL TAIL * * Bivariate normal distribution function BVND is required. * DOUBLE PRECISION SQ21, SQ31, RHO, SQTWPI, XCUT, BVND, ADONED DOUBLE PRECISION B1, B2, B3, B2P, B3P, RHO21, RHO31, RHO32, EPS PARAMETER ( SQTWPI = 2.506628274631000D0, XCUT = -8.5D0 ) PARAMETER ( EPS = 5D-16 ) COMMON /TRVBKD/B2P, B3P, RHO21, RHO31, RHO B1 = LIMIT(1) B2 = LIMIT(2) B3 = LIMIT(3) RHO21 = SIGMA(1) RHO31 = SIGMA(2) RHO32 = SIGMA(3) IF ( ABS(B2) .GE. MAX( ABS(B1), ABS(B3) ) ) THEN B1 = B2 B2 = LIMIT(1) RHO31 = RHO32 RHO32 = SIGMA(2) ELSE IF ( ABS(B3) .GE. MAX( ABS(B1), ABS(B2) ) ) THEN B1 = B3 B3 = LIMIT(1) RHO21 = RHO32 RHO32 = SIGMA(1) END IF TAIL = .FALSE. IF ( B1 .GT. 0 ) THEN TAIL = .TRUE. B1 = -B1 RHO21 = -RHO21 RHO31 = -RHO31 END IF IF ( B1 .GT. XCUT ) THEN IF ( 2*ABS(RHO21) .LT. 1 ) THEN SQ21 = SQRT( 1 - RHO21**2 ) ELSE SQ21 = SQRT( ( 1 - RHO21 )*( 1 + RHO21 ) ) END IF IF ( 2*ABS(RHO31) .LT. 1 ) THEN SQ31 = SQRT( 1 - RHO31**2 ) ELSE SQ31 = SQRT( ( 1 - RHO31 )*( 1 + RHO31 ) ) END IF RHO = ( RHO32 - RHO21*RHO31 )/(SQ21*SQ31) B2P = B2/SQ21 RHO21 = RHO21/SQ21 B3P = B3/SQ31 RHO31 = RHO31/SQ31 TVND = ADONED( TRVFND, XCUT, B1, EPS )/SQTWPI ELSE TVND = 0 END IF IF ( TAIL ) TVND = BVND( -B2, -B3, RHO32 ) - TVND END * DOUBLE PRECISION FUNCTION TRVFND(T) DOUBLE PRECISION T, B2, B3, RHO21, RHO31, RHO, BVND COMMON /TRVBKD/B2, B3, RHO21, RHO31, RHO TRVFND = EXP( -T*T/2 )*BVND( T*RHO21 - B2, T*RHO31 - B3, RHO ) END * DOUBLE PRECISION FUNCTION ADONED( F, A, B, TOL ) * * One Dimensional Globally Adaptive Integration Function * EXTERNAL F DOUBLE PRECISION F, A, B, TOL INTEGER NL, I, IM, IP PARAMETER ( NL = 100 ) DOUBLE PRECISION EI(NL), AI(NL), BI(NL), FI(NL), FIN, ERR, KRNRDD IP = 1 AI(IP) = A BI(IP) = B FI(IP) = KRNRDD( AI(IP), BI(IP), F, EI(IP) ) IM = 1 10 IM = IM + 1 BI(IM) = BI(IP) AI(IM) = (AI(IP)+BI(IP))/2 BI(IP) = AI(IM) FIN = FI(IP) FI(IP) = KRNRDD( AI(IP), BI(IP), F, EI(IP) ) FI(IM) = KRNRDD( AI(IM), BI(IM), F, EI(IM) ) ERR = ABS( FIN - FI(IP) - FI(IM) )/2 EI(IP) = EI(IP) + ERR EI(IM) = EI(IM) + ERR IP = 1 ERR = 0 FIN = 0 DO I = 1, IM IF ( EI(I) .GT. EI(IP) ) IP = I FIN = FIN + FI(I) ERR = ERR + EI(I) END DO IF ( ERR .GT. TOL .AND. IM .LT. NL ) GO TO 10 ADONED = FIN END * * DOUBLE PRECISION FUNCTION KRNRDD( A, B, F, ABSERR ) * * Kronrod Rule * EXTERNAL F DOUBLE PRECISION & A, ABSCIS, ABSERR, B, CENTER, F, FC, FUNSUM, HFLGTH, & RESLTG, RESLTK * * The abscissae and weights are given for the interval (-1,1) * because of symmetry only the positive abscisae and their * corresponding weights are given. * * XGK - abscissae of the 2N+1-point Kronrod rule: * XGK(2), XGK(4), ... N-point Gauss rule abscissae; * XGK(1), XGK(3), ... abscissae optimally added * to the N-point Gauss rule. * * WGK - weights of the 2N+1-point Kronrod rule. * * WG - weights of the N-point Gauss rule. * INTEGER J, N PARAMETER ( N = 11 ) * DOUBLE PRECISION WG(0:(N+1)/2), WGK(0:N), XGK(0:N) SAVE WG, WGK, XGK DATA WG( 1)/ 0.5566856711617449D-01/ DATA WG( 2)/ 0.1255803694649048D+00/ DATA WG( 3)/ 0.1862902109277352D+00/ DATA WG( 4)/ 0.2331937645919914D+00/ DATA WG( 5)/ 0.2628045445102478D+00/ DATA WG( 0)/ 0.2729250867779007D+00/ * DATA XGK( 1)/ 0.9963696138895427D+00/ DATA XGK( 2)/ 0.9782286581460570D+00/ DATA XGK( 3)/ 0.9416771085780681D+00/ DATA XGK( 4)/ 0.8870625997680953D+00/ DATA XGK( 5)/ 0.8160574566562211D+00/ DATA XGK( 6)/ 0.7301520055740492D+00/ DATA XGK( 7)/ 0.6305995201619651D+00/ DATA XGK( 8)/ 0.5190961292068118D+00/ DATA XGK( 9)/ 0.3979441409523776D+00/ DATA XGK(10)/ 0.2695431559523450D+00/ DATA XGK(11)/ 0.1361130007993617D+00/ DATA XGK( 0)/ 0.0000000000000000D+00/ * DATA WGK( 1)/ 0.9765441045961290D-02/ DATA WGK( 2)/ 0.2715655468210443D-01/ DATA WGK( 3)/ 0.4582937856442671D-01/ DATA WGK( 4)/ 0.6309742475037484D-01/ DATA WGK( 5)/ 0.7866457193222764D-01/ DATA WGK( 6)/ 0.9295309859690074D-01/ DATA WGK( 7)/ 0.1058720744813894D+00/ DATA WGK( 8)/ 0.1167395024610472D+00/ DATA WGK( 9)/ 0.1251587991003195D+00/ DATA WGK(10)/ 0.1312806842298057D+00/ DATA WGK(11)/ 0.1351935727998845D+00/ DATA WGK( 0)/ 0.1365777947111183D+00/ * * * List of major variables * * CENTER - mid point of the interval * HFLGTH - half-length of the interval * ABSCIS - abscissae * RESLTG - result of the N-point Gauss formula * RESLTK - result of the 2N+1-point Kronrod formula * * HFLGTH = ( B - A )/2 CENTER = ( B + A )/2 * * compute the 2N+1-point Kronrod approximation to * the integral, and estimate the absolute error. * FC = F(CENTER) RESLTG = FC*WG(0) RESLTK = FC*WGK(0) DO J = 1, N ABSCIS = HFLGTH*XGK(J) FUNSUM = F( CENTER - ABSCIS ) + F( CENTER + ABSCIS ) RESLTK = RESLTK + WGK(J)*FUNSUM IF( MOD( J, 2 ) .EQ. 0 ) RESLTG = RESLTG + WG(J/2)*FUNSUM END DO KRNRDD = RESLTK*HFLGTH ABSERR = 3*ABS( ( RESLTK - RESLTG )*HFLGTH ) END * DOUBLE PRECISION FUNCTION PHID(Z) * * Normal distribution probabilities accurate to 1.e-15. * Z = no. of standard deviations from the mean. * * Based upon algorithm 5666 for the error function, from: * Hart, J.F. et al, 'Computer Approximations', Wiley 1968 * * Programmer: Alan Miller * * Latest revision - 30 March 1986 * DOUBLE PRECISION P0, P1, P2, P3, P4, P5, P6, & Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7, & Z, P, EXPNTL, CUTOFF, ROOTPI, ZABS PARAMETER( & P0 = 220.20 68679 12376 1D0, & P1 = 221.21 35961 69931 1D0, & P2 = 112.07 92914 97870 9D0, & P3 = 33.912 86607 83830 0D0, & P4 = 6.3739 62203 53165 0D0, & P5 = .70038 30644 43688 1D0, & P6 = .035262 49659 98910 9D0 ) PARAMETER( & Q0 = 440.41 37358 24752 2D0, & Q1 = 793.82 65125 19948 4D0, & Q2 = 637.33 36333 78831 1D0, & Q3 = 296.56 42487 79673 7D0, & Q4 = 86.780 73220 29460 8D0, & Q5 = 16.064 17757 92069 5D0, & Q6 = 1.7556 67163 18264 2D0, & Q7 = .088388 34764 83184 4D0 ) PARAMETER( ROOTPI = 2.5066 28274 63100 1D0 ) PARAMETER( CUTOFF = 7.0710 67811 86547 5D0 ) * ZABS = ABS(Z) * * |Z| > 37 * IF ( ZABS .GT. 37 ) THEN P = 0 ELSE * * |Z| <= 37 * EXPNTL = EXP( -ZABS**2/2 ) * * |Z| < CUTOFF = 10/SQRT(2) * IF ( ZABS .LT. CUTOFF ) THEN P = EXPNTL*( ( ( ( ( ( P6*ZABS + P5 )*ZABS + P4 )*ZABS & + P3 )*ZABS + P2 )*ZABS + P1 )*ZABS + P0 ) & /( ( ( ( ( ( ( Q7*ZABS + Q6 )*ZABS + Q5 )*ZABS & + Q4 )*ZABS + Q3 )*ZABS + Q2 )*ZABS + Q1 )*ZABS + Q0 ) * * |Z| >= CUTOFF. * ELSE P = EXPNTL/( ZABS + 1/( ZABS + 2/( ZABS + 3/( ZABS & + 4/( ZABS + 0.65D0 ) ) ) ) )/ROOTPI END IF END IF IF ( Z .GT. 0 ) P = 1 - P PHID = P END -------------------------------------------------------------------------------- Alan C Genz 2001-02-25