cnd
// 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 normal distrubtion routine that does not require any external libraries. 
//	 This includes 3 differnt CND calculations: Hart (1968), Abromowitz and Stegun (1974) Polynomial Approximation, and
//	 Newton-Cotes using Boole's rule.
// comments:
//		The cumulative normal distribution function  integral has no closed-form solution, so a numerical approximation must be used. 
//	 	Recent research has shown that a high precision approximation can be of great importance in options valuation (West, 2005b). 
//		For this reason, we recommend the use of the Hart algorithm, which offers double precision. I have also included third- and 
//		fifth-degree polynomial approximations that are well known from the options literature. The polynomial approximations are known to 
//	 	suffer from inaccuracies far out in the tails and can in some special situations result in completely wrong options prices.
// reference:
//      The Complete Guide to Option Pricing Formulas, 2nd ed. (Espen Gaarder Haug)

library("cnd")
// @function Returns the Cumulative Normal Distribution (CND) using the Hart (1968) method.
// @param x float, 
// @returns float.
export CND1(float x)=>
    float Exponential = 0
    float SumA = 0
    float SumB = 0
    float CND = 0
    float y = math.abs(x)
    if y > 37 
        CND := 0
    else
        Exponential := math.exp(-math.pow(y, 2) / 2)
        if y < 7.07106781186547 
            SumA := 0.0352624965998911 * y + 0.700383064443688
            SumA := SumA * y + 6.37396220353165
            SumA := SumA * y + 33.912866078383
            SumA := SumA * y + 112.079291497871
            SumA := SumA * y + 221.213596169931
            SumA := SumA * y + 220.206867912376
            SumB := 0.0883883476483184 * y + 1.75566716318264
            SumB := SumB * y + 16.064177579207
            SumB := SumB * y + 86.7807322029461
            SumB := SumB * y + 296.564248779674
            SumB := SumB * y + 637.333633378831
            SumB := SumB * y + 793.826512519948
            SumB := SumB * y + 440.413735824752
            CND := Exponential * SumA / SumB
        else
            SumA := y + 0.65
            SumA := y + 4 / SumA
            SumA := y + 3 / SumA
            SumA := y + 2 / SumA
            SumA := y + 1 / SumA
            CND := Exponential / (SumA * 2.506628274631)
    CND := x > 0 ? 1 - CND : CND
    CND

// @function Returns the Cumulative Normal Distribution (CND) using the Abromowitz and Stegun (1974) Polynomial Approximation.
// @param x float, 
// @returns float.
export CND2(float x)=>
    float CND2 = 0
    if x == 0 
        CND2 := 0.5
    else
        float L = 0
        float k  = 0
        var float a1 = 0.31938153
        var float a2 = -0.356563782
        var float a3 = 1.781477937
        var float a4 = -1.821255978
        var float a5 = 1.330274429
        L := math.abs(x)
        k := 1 / (1 + 0.2316419 * L)
        CND2 := 1 - 1 / math.sqrt(2 * math.pi) 
			 * math.exp(-math.pow(L, 2) / 2) 
			 * (a1 * k 
			 + a2 * math.pow(k, 2) 
			 + a3 * math.pow(k, 3) 
			 + a4 * math.pow(k, 4) 
			 + a5 * math.pow(k, 5))
        if x < 0
            CND2 :=  1 - CND2 
    CND2

// internal normal distribution function
ND(float x)=> 
    float out = math.exp(-x * x * 0.5) / math.sqrt(2 * math.pi)
    out 
	
// internal Newton-Cotes method, Boole’s rule
Boole(float StartPoint, float EndPoint, int n)=>
    float[] X = array.new<float>(n + 1 , 0)
    float[] Y = array.new<float>(n + 1 , 0)
    float delta_x = (EndPoint - StartPoint) / n
    for i = 0 to n
        array.set(X, i, StartPoint + i * delta_x)
        array.set(Y, i, ND(array.get(X, i)))
    float sum = 0
    for t = 0 to (n - 1) / 4 
        int ind = 4 * t
        sum += (1 / 45.0) * 
             (14 * array.get(Y, ind) 
             + 64 * array.get(Y, ind + 1) 
             + 24 * array.get(Y, ind + 2) 
             + 64 * array.get(Y, ind + 3) 
             + 14 * array.get(Y, ind + 4)) 
             * delta_x
    sum

// @function Returns the Cumulative Normal Distribution (CND) using Newton-Cotes method, Boole’s rule
// @param x float, 
// @returns float.
export CND3(float x)=> 
    float out = Boole(-10.0, x, 240)
	out

// error function
erf(float xin)=>
    float[] cof = array.new<float>(28, 0)
    array.set(cof, 0, -1.3026537197817094)
    array.set(cof, 1, 6.4196979235649026e-1)
    array.set(cof, 2, 1.9476473204185836e-2)
    array.set(cof, 3, -9.561514786808631e-3)
    array.set(cof, 4, -9.46595344482036e-4)
    array.set(cof, 5, 3.66839497852761e-4)
    array.set(cof, 6, 4.2523324806907e-5)
    array.set(cof, 7, -2.0278578112534e-5)
    array.set(cof, 8, -1.624290004647e-6)
    array.set(cof, 9, 1.303655835580e-6)
    array.set(cof, 10, 1.5626441722e-8)
    array.set(cof, 11, -8.5238095915e-8)
    array.set(cof, 12, 6.529054439e-9)
    array.set(cof, 13, 5.059343495e-9)
    array.set(cof, 14, -9.91364156e-10)
    array.set(cof, 15, -2.27365122e-10)
    array.set(cof, 16, 9.6467911e-11)
    array.set(cof, 17, 2.394038e-12)
    array.set(cof, 18, -6.886027e-12)
    array.set(cof, 19,  8.94487e-13)
    array.set(cof, 20, 3.13092e-13)
    array.set(cof, 21, -1.12708e-13)
    array.set(cof, 22, 3.81e-16)
    array.set(cof, 23, 7.106e-15)
    array.set(cof, 24, -1.523e-15)
    array.set(cof, 25, -9.4e-17)
    array.set(cof, 26, 1.21e-16)
    array.set(cof, 27, -2.8e-17)

    bool isneg = false
    float d = 0
    float dd = 0
    float x = xin

    if xin < 0
        x := -xin
        isneg := true

    float t = 2 / (2 + x)
    float ty = 4 * t - 2

    for j = array.size(cof) - 1 to 1  
        float tmp = d
        d := ty * d - dd + array.get(cof, j)
        dd := tmp

    float res = t * math.exp(-x * x + 0.5 * (array.get(cof, 0) + ty * d) - dd)
    out = isneg ? res - 1 : 1 - res
    out

// Complementary Error Function
erfc(float x)=>
    out = 1 - erf(x)
    out

// @function Returns the standard normal cumulative distribution function. The distribution has a mean of 0 (zero) and a standard deviation of one. Use this function in place of a table of standard normal curve areas.
// @param Z float, The value for which you want the distribution.
// @returns float, standard normal cumulative distribution function
export CND4(float a)=>
    float NPY_SQRT1_2 = 0.707106781186547524400844362104849039 
    float x = 0
    float y = 0
    float z = 0
    if na(a) 
        runtime.error("ndtr") 
    x := a * NPY_SQRT1_2
    z := math.abs(x)
    if z < NPY_SQRT1_2
        y := 0.5 + 0.5 * erf(x)
    else 
        y := 0.5 * erfc(z)
        if x > 0
            y := 1.0 - y
    y

// example usage:
inx = input.float(0.5)
out1 = CND1(inx)
out2 = CND2(inx)
out3 = CND3(inx)
out4 = CND4(inx)

plot(out1, "out1")
plot(out2, "out2")
plot(out3, "out3")
plot(out4, "out4")
