cbnd
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

//@version=5
// @description a standalone Cumulative Bivariate Normal Distribution (CBND) functions that do not require any external libraries. 
//	 This includes 3 differnt CBND calculations: Drezner(1978), Drezner and Wesolowsky (1990), and Genz (2004)
// comments:
// 		The standardized cumulative normal distribution function returns the probability that one random 
// 		variable is less than a and that a second random variable is less than b when the correlation 
// 		between the two variables is p. Since no closed-form solution exists for the bivariate cumulative 
// 		normal distribution, we present three approximations. The first one is the well-known 
// 		Drezner (1978) algorithm. The second one is the more efficient Drezner and Wesolowsky (1990) 
// 		algorithm. The third is the Genz (2004) algorithm, which is the most accurate one and therefore
// 		our recommended algorithm. West (2005b) and Agca and Chance (2003) discuss the speed and 
// 		accuracy of bivariate normal distribution approximations for use in option pricing in 
// 		ore detail.
// reference:
//      The Complete Guide to Option Pricing Formulas, 2nd ed. (Espen Gaarder Haug)
library("cbnd")

//import Cumulative Normal Distribution
import loxx/cnd/1 

// @function Returns the Cumulative Bivariate Normal Distribution (CBND) using Drezner 1978 Algorithm
// @param A float, 
// @param b float, 
// @param rho float, 
// @returns float.
export CBND1(float A, float b, float rho)=>
    float g = 0
    float P = 0
    float sum = 0
    float CBND2 = 0

    array<float> x = array.new<float>(5, 0)
    array.set(x, 0, 0.018854042)
    array.set(x, 1, 0.038088059)
    array.set(x, 2, 0.0452707394)
    array.set(x, 3, 0.038088059)
    array.set(x, 4, 0.018854042)

    array<float> y = array.new<float>(5, 0)
    array.set(y, 0, 0.04691008)
    array.set(y, 1, 0.23076534) 
    array.set(y, 2, 0.5)
    array.set(y, 3, 0.76923466)
    array.set(y, 4, 0.95308992)
    
    sum := 0
    for i = 0 to 4
        P := array.get(y, i) * rho
        g := 1 - P * P
        sum += array.get(x, i) * math.exp((2 * A * b * P - A * A - b * b) / g / 2) / math.sqrt(g)
    CBND2 := rho * sum + cnd.CND1(A) * cnd.CND1(b)
    CBND2

// @function Returns the Cumulative Bivariate Normal Distribution (CBND) using Drezner and Wesolowsky (1990) function
// @param A float, 
// @param b float, 
// @param rho float, 
// @returns float.
export CBND2(float A, float b, float rho)=>
// modified/corrected from the second function in Drez & Wes paper pg. 105
// 0/0 case resolved by l'H rule
    float h1 = 0
    float h2 = 0
    float LH = 0
    float h12 = 0 
    float h3 = 0
    float h5 = 0
    float h6 = 0
    float h7 = 0
    float h8 = 0
    float r1 = 0
    float r2 = 0
    float r3 = 0
    float rr = 0
    float AA = 0
    float ab = 0
    float CBND4 = 0

    array<float> x = array.new<float>(5, 0)
    array.set(x, 0, 0.04691008)
    array.set(x, 1, 0.23076534)
    array.set(x, 2, 0.5)
    array.set(x, 3, 0.76923466)
    array.set(x, 4, 0.95308992)

    array<float> W = array.new<float>(5, 0)
    array.set(W, 0, 0.018854042)
    array.set(W, 1, 0.038088059)
    array.set(W, 2, 0.0452707394)
    array.set(W, 3, 0.038088059)
    array.set(W, 4, 0.018854042)

    h1 := A
    h2 := b
    h12 := (h1 * h1 + h2 * h2) / 2

    if math.abs(rho) >= 0.7 
        r2 := 1 - rho * rho
        r3 := math.sqrt(r2)
        if rho < 0  
            h2 := -h2
        h3 := h1 * h2
        h7 := math.exp(-h3 / 2)
        if math.abs(rho) < 1 
            h6 := math.abs(h1 - h2)
            h5 := h6 * h6 / 2
            h6 := h6 / r3
            AA := 0.5 - h3 / 8
            ab := 3 - 2 * AA * h5
            LH := 0.13298076 * h6 * ab * (1 - cnd.CND1(h6)) - math.exp(-h5 / r2) * (ab + AA * r2) * 0.053051647
            for i = 0 to 4
                r1 := r3 * array.get(x, i)
                rr := r1 * r1
                r2 := math.sqrt(1 - rr)
                if h7 == 0 
                    h8 := 0
                else
                    h8 := math.exp(-h3 / (1 + r2)) / r2 / h7
                LH := LH - array.get(W, i) * math.exp(-h5 / rr) * (h8 - 1 - AA * rr)

        CBND4 := LH * r3 * h7 + cnd.CND1(math.min(h1, h2))
        if rho < 0 
            CBND4 := cnd.CND1(h1) - CBND4
    else
        h3 := h1 * h2
        if rho != 0 
            for i = 0 to 4
                r1 := rho * array.get(x, i)
                r2 := 1 - r1 * r1
                LH := LH + array.get(W, i) * math.exp((r1 * h3 - h12) / r2) / math.sqrt(r2)
        CBND4 := cnd.CND1(h1) * cnd.CND1(h2) + rho * LH
    CBND4


// @function Returns the Cumulative Bivariate Normal Distribution (CBND) using Genz (2004) algorithm (this is the preferred method)
// @param x float, 
// @param y float, 
// @param rho float, 
// @returns float.
export CBND3(float x, float y, float rho)=>
// 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, (1990),
// On the computation of the bivariate normal integral,
// Journal of Statist. Comput. Simul. 35, pp. 101-107,
//  with major modifications for double precision, and for |R| close to 1.
// This code was originally transelated into VBA by Graeme West

    int LG = 0
    int NG = 0

    float h = 0
    float k = 0
    float hk = 0
    float hs = 0
    float BVN = 0
    float Ass = 0
    float asr = 0
    float sn = 0
    float A = 0
    float b = 0
    float bs = 0
    float c = 0
    float d = 0
    float xs = 0
    float rs = 0
    float CBND = 0

    matrix<float> XX = matrix.new<float>(11, 4, 0)
    matrix<float> W = matrix.new<float>(11, 4, 0)

    matrix.set(W,  1, 1, 0.17132449237917)
    matrix.set(XX, 1, 1, -0.932469514203152)
    matrix.set(W,  2, 1, 0.360761573048138)
    matrix.set(XX, 2, 1, -0.661209386466265)
    matrix.set(W,  3, 1, 0.46791393457269)
    matrix.set(XX, 3, 1, -0.238619186083197)

    matrix.set(W,  1, 2, 0.0471753363865118)
    matrix.set(XX, 1, 2, -0.981560634246719)
    matrix.set(W,  2, 2, 0.106939325995318)
    matrix.set(XX, 2, 2, -0.904117256370475)
    matrix.set(W,  3, 2, 0.160078328543346)
    matrix.set(XX, 3, 2, -0.769902674194305)
    matrix.set(W,  4, 2, 0.203167426723066)
    matrix.set(XX, 4, 2, -0.587317954286617)
    matrix.set(W,  5, 2, 0.233492536538355)
    matrix.set(XX, 5, 2, -0.36783149899818)
    matrix.set(W,  6, 2, 0.249147045813403)
    matrix.set(XX, 6, 2, -0.125233408511469)

    matrix.set(W,  1, 3, 0.0176140071391521)
    matrix.set(XX, 1, 3, -0.993128599185095)
    matrix.set(W,  2, 3, 0.0406014298003869)
    matrix.set(XX, 2, 3, -0.963971927277914)
    matrix.set(W,  3, 3, 0.0626720483341091)
    matrix.set(XX, 3, 3, -0.912234428251326)
    matrix.set(W,  4, 3, 0.0832767415767048)
    matrix.set(XX, 4, 3, -0.839116971822219)
    matrix.set(W,  5, 3, 0.10193011981724)
    matrix.set(XX, 5, 3, -0.746331906460151)
    matrix.set(W,  6, 3, 0.118194531961518)
    matrix.set(XX, 6, 3, -0.636053680726515)
    matrix.set(W,  7, 3, 0.131688638449177)
    matrix.set(XX, 7, 3, -0.510867001950827)
    matrix.set(W,  8, 3, 0.142096109318382)
    matrix.set(XX, 8, 3, -0.37370608871542)
    matrix.set(W,  9, 3, 0.149172986472604)
    matrix.set(XX, 9, 3, -0.227785851141645)
    matrix.set(W,  10, 3, 0.152753387130726)
    matrix.set(XX, 10, 3, -0.0765265211334973)
        
    if math.abs(rho) < 0.3 
        NG := 1
        LG := 3
    else if  math.abs(rho) < 0.75 
        NG := 2
        LG := 6
    else
        NG := 3
        LG := 10
        
    h := -x
    k := -y
    hk := h * k
    BVN := 0
        
    if math.abs(rho) < 0.925 
        if math.abs(rho) > 0 
            hs := (h * h + k * k) / 2
            asr := math.asin(rho)
            for i = 1 to LG
                for ISs = -1 to 1 by 2
                    sn := math.sin(asr * (ISs * matrix.get(XX, i, NG) + 1) / 2)
                    BVN := BVN + matrix.get(W, i, NG) * math.exp((sn * hk - hs) / (1 - sn * sn))

            BVN := BVN * asr / (4 * math.pi)
        BVN := BVN + cnd.CND1(-h) * cnd.CND1(-k)
    else
        if rho < 0 
            k := -k
            hk := -hk
        if math.abs(rho) < 1 
            Ass := (1 - rho) * (1 + rho)
            A := math.sqrt(Ass)
            bs := math.pow(h - k, 2)
            c := (4 - hk) / 8
            d := (12 - hk) / 16
            asr := -(bs / Ass + hk) / 2
            if asr > -100  
                BVN := A * math.exp(asr) * (1 - c * (bs - Ass) * (1 - d * bs / 5) / 3 + c * d * Ass * Ass / 5)
            if -hk < 100 
                b := math.sqrt(bs)
                BVN := BVN - math.exp(-hk / 2) * math.sqrt(2 * math.pi) * cnd.CND1(-b / A) * b * (1 - c * bs * (1 - d * bs / 5) / 3)
            A := A / 2
            for i = 1 to LG
                for ISs = -1 to 1 by 2
                    xs := math.pow(A * (ISs * matrix.get(XX, i, NG) + 1), 2)
                    rs := math.sqrt(1 - xs)
                    asr := -(bs / xs + hk) / 2
                    if asr > -100 
                        BVN := BVN + A * matrix.get(W, i, NG) * math.exp(asr) * (math.exp(-hk * (1 - rs) / (2 * (1 + rs))) / rs - (1 + c * xs * (1 + d * xs)))
            BVN := -BVN / (2 * math.pi)
        if rho > 0 
            BVN := BVN + cnd.CND1(-math.max(h, k))
        else
            BVN := -BVN
            if k > h  
                BVN := BVN + cnd.CND1(k) - cnd.CND1(h)
    CBND := BVN
    CBND

x = input.float(0.5)
y = input.float(0.5)
rho = input.float(0.95)

out1 = CBND1(x, y, rho)
out2 = CBND2(x, y, rho)
out3 = CBND3(x, y, rho)

plot(out1, "Drezner 1978 ")
plot(out2, "Drezner-Weso-90a")
plot(out3, "Genze")
