MODULE Bivariate_Trivariate_Normal

! Appendix: Fortran Source Code
 
! Code converted using TO_F90 by Alan Miller
! Date: 2001-02-25  Time: 11:57:13
! Latest revision - 25 February 2001
! amiller @ bigpond.net.au
 
!     Code from Alan Genz's web site for bivariate & trivariate
!     normal integrals, adapted by Alan Miller.

!     http://www.sci.wsu.edu/math/faculty/genz/papers/

!     Alan C Genz
!     2001-01-22


IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

! COMMON /trvbkd/b2p, b3p, rho21, rho31, rho
REAL (dp), SAVE  :: b2p, b3p, rho21, rho31, rho


CONTAINS


FUNCTION bvn( sh, sk, r ) RESULT(fn_val)

!     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, INTENT(IN)  :: sh, sk, r
REAL              :: fn_val

REAL, PARAMETER     :: zero = 0, twopi = 6.283185
INTEGER, PARAMETER  :: lg = 5
REAL, PARAMETER     :: x(lg) =  &
                       (/ -0.9061798, -0.5384693, 0.0, 0.5384693, 0.9061798 /)
!          X(I) = +-SQRT((5+-SQRT(10/7))/9), 0
REAL, PARAMETER     :: w(lg) =  &
                       (/ 0.2369269,  0.4786287, 0.5688889, 0.4786287, 0.2369269 /)

INTEGER  :: i
REAL     :: as, a, b, c, rs, xs
REAL     :: sn, asr, h, k, bs, hs, hk

h = sh
k = sk
hk = h*k
fn_val = 0

!     Absolute value of the correlation coefficient is less than 0.8

IF ( ABS(r) < 0.8 ) THEN
  hs = ( h*h + k*k )/2
  asr = ASIN( r )/2
  DO i = 1, lg
    sn = SIN( asr*(x(i)+1) )
    fn_val = fn_val + w(i)*EXP( ( sn*hk - hs )/( 1 - sn*sn ) )
  END DO
  fn_val = asr*fn_val/twopi + phi(-h)*phi(-k)
ELSE
  
!     For larger correlation coefficient
  
  IF ( r < 0 ) THEN
    k = -k
    hk = -hk
  END IF
  IF ( ABS(r) < 1 ) THEN
    as = ( 1 - r )*( 1 + r )
    a = SQRT( as )
    bs = ( h - k )**2
    c = ( 4 - hk )/8
    fn_val = a*EXP( -(bs/as + hk)/2 )*( 1 - c*(bs - as)/3 )
    IF ( -hk < 100 ) THEN
      b = SQRT(bs)
      fn_val = fn_val - EXP(-hk/2)*SQRT(twopi)*phi(-b/a)*b*(1-c*bs/3)
    END IF
    a = a/2
    DO i = 1, lg
      xs = ( a*( x(i) + 1 ) )**2
      rs = SQRT( 1 - xs )
      asr = -( bs/xs + hk )/2
      IF ( asr > -100 ) fn_val = fn_val + a*w(i)*   &
                        ( EXP(-bs/(2*xs) - hk/(1+rs))/rs - EXP(asr)*(1+c*xs) )
    END DO
    fn_val = -fn_val/twopi
  END IF
  IF ( r > 0 ) fn_val =  fn_val + phi( -MAX( h, k ) )
  IF ( r < 0 ) fn_val = -fn_val + MAX( zero, phi(-h) - phi(-k) )
END IF

RETURN
END FUNCTION bvn



FUNCTION phi(x) RESULT(fn_val)

REAL, INTENT(IN)  :: x
REAL              :: fn_val

fn_val = phid( REAL(x, KIND=dp) )
RETURN
END FUNCTION phi



FUNCTION bvnd( dh, dk, r ) RESULT(fn_val)

!     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 REAL (dp), and for |R| close to 1.

! BVND - calculate the probability that X is larger than DH, and Y is
!       larger than DK.

! Parameters

!   DH  REAL (dp), integration limit
!   DK  REAL (dp), integration limit
!   R   REAL (dp), correlation coefficient


REAL (dp), INTENT(IN)  :: dh, dk, r
REAL (dp)              :: fn_val

REAL (dp), PARAMETER  :: zero = 0.0_dp, twopi = 6.283185307179586_dp

REAL (dp), PARAMETER  :: x(10,3) = RESHAPE(  &
                             !     Gauss Legendre Points, N =  6
     (/ -0.9324695142031522_dp, -0.6612093864662647_dp,  &
        -0.2386191860831970_dp,  0.0_dp, 0.0_dp, 0.0_dp, &
         0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp,  &
                             !     Gauss Legendre Points, N = 12
        -0.9815606342467191_dp, -0.9041172563704750_dp,  &
        -0.7699026741943050_dp, -0.5873179542866171_dp,  &
        -0.3678314989981802_dp, -0.1252334085114692_dp,  &
         0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp,  &
                             !     Gauss Legendre Points, N = 20
        -0.9931285991850949_dp, -0.9639719272779138_dp,  &
        -0.9122344282513259_dp, -0.8391169718222188_dp,  &
        -0.7463319064601508_dp, -0.6360536807265150_dp,  &
        -0.5108670019508271_dp, -0.3737060887154196_dp,  &
        -0.2277858511416451_dp, -0.07652652113349733_dp /), (/ 10, 3 /) )
REAL (dp), PARAMETER  :: w(10,3) = RESHAPE(  &
                             !     Gauss Legendre Weights, N =  6
  (/ 0.1713244923791705_dp, 0.3607615730481384_dp, 0.4679139345726904_dp, &
     0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp,  &
                          !     Gauss Legendre Weights, N = 12
     0.04717533638651177_dp, 0.1069393259953183_dp, 0.1600783285433464_dp, &
     0.2031674267230659_dp,  0.2334925365383547_dp, 0.2491470458134029_dp, &
     0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp,  &
                          !     Gauss Legendre Weights, N = 20
     0.01761400713915212_dp, 0.04060142980038694_dp, 0.06267204833410906_dp, &
     0.08327674157670475_dp, 0.1019301198172404_dp, 0.1181945319615184_dp, &
     0.1316886384491766_dp,  0.1420961093183821_dp, 0.1491729864726037_dp, &
     0.1527533871307259_dp /), (/ 10, 3 /) )

INTEGER    :: i, is, lg, ng
REAL (dp)  :: as, a, b, c, d, rs, xs, bvn
REAL (dp)  :: sn, asr, h, k, bs, hs, hk

IF ( ABS(r) < 0.3_dp ) THEN
  ng = 1
  lg = 3
ELSE IF ( ABS(r) < 0.75_dp ) THEN
  ng = 2
  lg = 6
ELSE
  ng = 3
  lg = 10
END IF
h = dh
k = dk
hk = h*k
bvn = 0
IF ( ABS(r) < 0.925_dp ) 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 < 0 ) THEN
    k = -k
    hk = -hk
  END IF
  IF ( ABS(r) < 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 > -100 ) bvn = a*EXP(asr)  &
                       *( 1 - c*( bs - as )*( 1 - d*bs/5 )/3 + c*d*as*as/5 )
    IF ( -hk < 100 ) THEN
      b = SQRT(bs)
      bvn = bvn - EXP( -hk/2 )*SQRT(twopi)*phid(-b/a)*b  &
            *( 1 - c*bs*( 1 - d*bs/5 )/3 )
    END IF
    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 > -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
  END IF
  IF ( r > 0 ) bvn =  bvn + phid( -MAX( h, k ) )
  IF ( r < 0 ) bvn = -bvn + MAX( zero, phid(-h) - phid(-k) )
END IF

fn_val = bvn
RETURN
END FUNCTION bvnd



FUNCTION tvnd( limit, sigma ) RESULT(fn_val)

!     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.

!     Bivariate normal distribution function BVND is required.

! Parameters

!   LIMIT  REAL (dp) array of three upper integration limits.
!   SIGMA  REAL (dp) 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.


REAL (dp), INTENT(IN)  :: limit(3), sigma(3)
REAL (dp)              :: fn_val

! EXTERNAL trvfnd
! COMMON /trvbkd/b2p, b3p, rho21, rho31, rho

LOGICAL   :: tail
REAL (dp) :: sq21, sq31
REAL (dp) :: b1, b2, b3, rho32
REAL (dp), PARAMETER :: sqtwpi = 2.506628274631000_dp, xcut = -8.5_dp,  &
                        eps = 5D-16

b1 = limit(1)
b2 = limit(2)
b3 = limit(3)
rho21 = sigma(1)
rho31 = sigma(2)
rho32 = sigma(3)
IF ( ABS(b2) >= MAX( ABS(b1), ABS(b3) ) ) THEN
  b1 = b2
  b2 = limit(1)
  rho31 = rho32
  rho32 = sigma(2)
ELSE IF ( ABS(b3) >= MAX( ABS(b1), ABS(b2) ) ) THEN
  b1 = b3
  b3 = limit(1)
  rho21 = rho32
  rho32 = sigma(1)
END IF
tail = .false.
IF ( b1 > 0 ) THEN
  tail = .true.
  b1 = -b1
  rho21 = -rho21
  rho31 = -rho31
END IF
IF ( b1 > xcut ) THEN
  IF ( 2*ABS(rho21) < 1 ) THEN
    sq21 = SQRT( 1 - rho21**2 )
  ELSE
    sq21 = SQRT( ( 1 - rho21 )*( 1 + rho21 ) )
  END IF
  IF ( 2*ABS(rho31) < 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
  fn_val = adoned( trvfnd, xcut, b1, eps )/sqtwpi
ELSE
  fn_val = 0
END IF
IF ( tail ) fn_val = bvnd( -b2, -b3, rho32 ) - fn_val

RETURN
END FUNCTION tvnd



FUNCTION trvfnd(t) RESULT(fn_val)

REAL (dp), INTENT(IN)  :: t
REAL (dp)              :: fn_val

! COMMON /trvbkd/b2, b3, rho21, rho31, rho

fn_val = EXP( -t*t/2 )*bvnd(t*rho21 - b2p, t*rho31 - b3p, rho)
RETURN
END FUNCTION trvfnd



FUNCTION adoned( f, a, b, tol ) RESULT(fin)

!     One Dimensional Globally Adaptive Integration Function

REAL (dp), INTENT(IN)  :: a, b, tol
REAL (dp)              :: fin

! EXTERNAL f

INTERFACE
  FUNCTION f(x) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)  :: x
    REAL (dp)              :: fn_val
  END FUNCTION f
END INTERFACE

INTEGER, PARAMETER  :: nl = 100
REAL (dp)  :: ei(nl), ai(nl), bi(nl), fi(nl), ERR
INTEGER    :: i, im, ip

ip = 1
ai(ip) = a
bi(ip) = b
CALL krnrdd( ai(ip), bi(ip), f, ei(ip), fi(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)
CALL krnrdd( ai(ip), bi(ip), f, ei(ip), fi(ip) )
CALL krnrdd( ai(im), bi(im), f, ei(im), fi(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) > ei(ip) ) ip = i
  fin = fin + fi(i)
  ERR = ERR + ei(i)
END DO
IF ( ERR > tol .AND. im < nl ) GO TO 10

RETURN
END FUNCTION adoned



SUBROUTINE krnrdd( a, b, f, abserr, fn_val)

!     Kronrod Rule

REAL (dp), INTENT(IN)   :: a, b
REAL (dp), INTENT(OUT)  :: abserr, fn_val

! EXTERNAL f

INTERFACE
  FUNCTION f(x) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)  :: x
    REAL (dp)              :: fn_val
  END FUNCTION f
END INTERFACE

REAL (dp)  :: abscis, center, 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
INTEGER, PARAMETER :: n = 11

REAL (dp), PARAMETER  :: wg(0:n/2) = (/  &
     0.2729250867779007_dp, 0.05566856711617449_dp, 0.1255803694649048_dp,  &
     0.1862902109277352_dp, 0.2331937645919914_dp, 0.2628045445102478_dp  /)
REAL (dp), PARAMETER  :: xgk(0:n) = (/  &
     0.0000000000000000_dp, 0.9963696138895427_dp, 0.9782286581460570_dp,   &
     0.9416771085780681_dp, 0.8870625997680953_dp, 0.8160574566562211_dp,   &
     0.7301520055740492_dp, 0.6305995201619651_dp, 0.5190961292068118_dp,   &
     0.3979441409523776_dp, 0.2695431559523450_dp, 0.1361130007993617_dp /)
REAL (dp), PARAMETER  :: wgk(0:n) = (/  &
     0.1365777947111183_dp, 0.9765441045961290D-02, 0.2715655468210443D-01, &
     0.4582937856442671D-01, 0.6309742475037484D-01, 0.7866457193222764D-01, &
     0.9295309859690074D-01, 0.1058720744813894_dp, 0.1167395024610472_dp,  &
     0.1251587991003195_dp, 0.1312806842298057_dp, 0.1351935727998845_dp /)

!     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 ) == 0 ) resltg = resltg + wg(j/2)*funsum
END DO
fn_val = resltk*hflgth
abserr = 3*ABS( ( resltk - resltg )*hflgth )

RETURN
END SUBROUTINE krnrdd



FUNCTION phid(z) RESULT(fn_val)

! Normal distribution probabilities accurate to 1.e-15.
! Z = no. of standard deviations from the mean.
! P, Q = probabilities to the left & right of Z.   P + Q = 1.
!       PDF = the probability density.

!       Based upon algorithm 5666 for the error function, from:
!       Hart, J.F. et al, 'Computer Approximations', Wiley 1968

!       Programmer: Alan Miller, modified by Alan Genz

! Latest revision of Fortran 77 version - 30 March 1986
! Latest revision of Fortran 90 version - 23 June 1997

IMPLICIT NONE
INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(12, 60)
REAL (dp), INTENT(IN)  :: z
REAL (dp)              :: fn_val

! Local variables
REAL (dp), PARAMETER :: p0 = 220.2068679123761_dp, p1 = 221.2135961699311_dp,  &
                        p2 = 112.0792914978709_dp, p3 = 33.91286607838300_dp,  &
                        p4 = 6.373962203531650_dp, p5 = .7003830644436881_dp,  &
                        p6 = 0.03526249659989109_dp,  &
                        q0 = 440.4137358247522_dp, q1 = 793.8265125199484_dp,  &
                        q2 = 637.3336333788311_dp, q3 = 296.5642487796737_dp,  &
                        q4 = 86.78073220294608_dp, q5 = 16.06417757920695_dp,  &
                        q6 = 1.755667163182642_dp, q7 = .8838834764831844D-1, &
                        cutoff = 7.071067811865475_dp

REAL (dp)            :: zabs, expntl

zabs = ABS(z)

! |Z| > 37.

IF (zabs > 37._dp) THEN
    fn_val = 0._dp
  ELSE

! |Z| <= 37.

    expntl = EXP(-0.5_dp*zabs**2)

! |Z| < CUTOFF = 10/sqrt(2).

  IF (zabs < cutoff) THEN
      fn_val = 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
    fn_val = expntl/(zabs + 1._dp/(zabs + 2._dp/(zabs + 3._dp/(zabs +  &
                     4._dp/(zabs + 0.65_dp)))))
  END IF
END IF

IF (z > 0.0_dp) fn_val = 1.0_dp - fn_val

RETURN
END FUNCTION phid

END MODULE Bivariate_Trivariate_Normal