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