Option Base 1

'Bivariate Cumulative Normal Distribution Approximation (As given by Hull, etc)
Public Function BivN(a As Double, b As Double, rho As Double) As Double

    Dim aarray, barray
    Dim rho1 As Double, rho2 As Double, delta As Double
    Dim ap As Double, bp As Double
    Dim Sum As Double
    Dim i As Integer, j As Integer, Pi As Double
    Dim denum As Double
    
    aarray = Array(0.24840615, 0.39233107, 0.21141819, 0.03324666, 0.00082485334)
    barray = Array(0.10024215, 0.48281397, 1.0609498, 1.7797294, 2.6697604)
    
    Pi = 3.14159265358979
    
    ap = a / Sqr(2 * (1 - rho ^ 2))
    bp = b / Sqr(2 * (1 - rho ^ 2))
    
    If a <= 0 And b <= 0 And rho <= 0 Then
        Sum = 0
        For i = 1 To 5
            For j = 1 To 5
                Sum = Sum + aarray(i) * aarray(j) * f(barray(i), barray(j), ap, bp, rho)
            Next
        Next
        BivN = Sum * Sqr(1 - rho ^ 2) / Pi
        
    ElseIf a <= 0 And b >= 0 And rho >= 0 Then
        BivN = Norm(a) - BivN(a, -b, -rho)
        
    ElseIf a >= 0 And b <= 0 And rho >= 0 Then
        BivN = Norm(b) - BivN(-a, b, -rho)
    ElseIf a >= 0 And b >= 0 And rho <= 0 Then
        BivN = Norm(a) + Norm(b) - 1 + BivN(-a, -b, rho)
    ElseIf a * b * rho >= 0 Then
        denum = Sqr(a ^ 2 - 2 * rho * a * b + b ^ 2)
        rho1 = (rho * a - b) * Sgn(a) / denum
        rho2 = (rho * b - a) * Sgn(b) / denum
        delta = (1 - Sgn(a) * Sgn(b)) / 4
        BivN = BivN(a, 0, rho1) + BivN(b, 0, rho2) - delta
    End If
    
End Function


Public Function f(X, y, ap, bp, rho)

Dim r
r = ap * (2 * X - ap) + bp * (2 * y - bp) + 2 * rho * (X - ap) * (y - bp)
f = Exp(r)
End Function

Public Function Norm(x As Double) As Double

Norm = Application.NormSDist(x)

End Function