CDF of the normal distribution in .net or c++ math libraries?

Discussion in 'Options' started by thenmmm, Nov 23, 2010.

  1. thenmmm

    thenmmm

    Is anyone aware if this function is available as a standard namespace or a header file? Currently i am approximating it with a stirling approximation which gives very small error...but still.
    I might post the c# code later....thanks.
     
  2. sjfan

    sjfan

    Here's mine - works to practically useful precision

    public class NormalDistribution
    {
    private static double A1 = -3.969683028665376e+01;
    private static double A2 = 2.209460984245205e+02;
    private static double A3 = -2.759285104469687e+02;
    private static double A4 = 1.383577518672690e+02;
    private static double A5 = -3.066479806614716e+01;
    private static double A6 = 2.506628277459239e+00;

    private static double B1 = -5.447609879822406e+01;
    private static double B2 = 1.615858368580409e+02;
    private static double B3 = -1.556989798598866e+02;
    private static double B4 = 6.680131188771972e+01;
    private static double B5 = -1.328068155288572e+01;

    private static double C1 = -7.784894002430293e-03;
    private static double C2 = -3.223964580411365e-01;
    private static double C3 = -2.400758277161838e+00;
    private static double C4 = -2.549732539343734e+00;
    private static double C5 = 4.374664141464968e+00;
    private static double C6 = 2.938163982698783e+00;

    private static double D1 = 7.784695709041462e-03;
    private static double D2 = 3.224671290700398e-01;
    private static double D3 = 2.445134137142996e+00;
    private static double D4 = 3.754408661907416e+00;

    private static double P_LOW = 0.02425;
    // P_high = 1 - p_low
    private static double P_HIGH = 0.97575;

    /// <summary>
    /// Find the Z' such that P = Pr[Z less than Z'];
    /// </summary>
    /// <param name="P"></param>
    public static double InverseNormal(double p)
    {
    double x = 0;
    double q, r;

    if ((0 < p) && (p < P_LOW))
    {
    q = System.Math.Sqrt(-2 * System.Math.Log(p));
    x = (((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);
    }
    else
    {
    if ((P_LOW <= p) && (p <= P_HIGH))
    {
    q = p - 0.5;
    r = q * q;
    x = (((((A1 * r + A2) * r + A3) * r + A4) * r + A5) * r + A6) * q / (((((B1 * r + B2) * r + B3) * r + B4) * r + B5) * r + 1);
    }
    else
    {
    if ((P_HIGH < p) && (p < 1))
    {
    q = System.Math.Sqrt(-2 * System.Math.Log(1 - p));
    x = -(((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);
    }
    }
    }


    return x;
    }

    /// <summary>
    /// Support function for the bivariate normal
    /// </summary>
    /// <param name="x"></param>
    /// <param name="y"></param>
    /// <param name="aprime"></param>
    /// <param name="bprime"></param>
    /// <param name="rho"></param>
    /// <returns></returns>
    private static double f(double x, double y, double aprime, double bprime, double rho)
    {
    double r = aprime * (2 * x - aprime) + bprime * (2 * y - bprime) + 2 * rho * (x - aprime) * (y - bprime);
    return Math.Exp(r);
    }

    public static double NormPDF(double x)
    {
    return NormalDistribution.NormPDF(x, 0, 1);
    }
    public static double NormPDF(double x, double mu, double sigma)
    {
    double xn = (x - mu) / sigma;
    return System.Math.Exp(-0.5 * xn * xn) / (System.Math.Sqrt(2 * System.Math.PI) * sigma);
    }

    /// <summary>
    /// This is the cumulative bivariate normal distribution. Based on
    /// http://finance-old.bi.no/~bernt/gcc_prog/recipes/recipes/node23.html
    /// </summary>
    /// <param name="a"></param>
    /// <param name="b"></param>
    /// <param name="rho"></param>
    /// <returns></returns>
    public static double BivariateNormalCDF(double a, double b, double rho)
    {
    if ((a <= 0.0) && (b <= 0.0) && (rho <= 0.0))
    {
    double aprime = a / Math.Sqrt(2.0 * (1.0 - rho * rho));
    double bprime = b / Math.Sqrt(2.0 * (1 - rho * rho));
    double[] A = new double[4] { 0.3253030, 0.4211071, 0.1334425, 0.006374323 };
    double[] B = new double[4] { 0.1337764, 0.6243247, 1.3425378, 2.2626645 };

    double sum = 0.0;

    for (int i = 0; i < 4; i++)
    for (int j = 0; j < 4; j++)
    sum = sum + A * A[j] * f(B, B[j], aprime, bprime, rho);

    sum = sum * (Math.Sqrt(1.0 - rho * rho) / Math.PI);
    return sum;
    }
    else if (a * b * rho <= 0.0)
    {
    if ((a <= 0.0) && (b >= 0.0) && (rho >= 0.0))
    return NormalCDF(a) - BivariateNormalCDF(a, -1.0 * b, -1.0 * rho);
    else if ((a >= 0.0) && (b <= 0.0) && (rho >= 0.0))
    return NormalCDF(b) - BivariateNormalCDF(-1.0 * a, b, -1.0 * rho);
    else if ((a >= 0.0) && (b >= 0.0) && (rho <= 0.0))
    return NormalCDF(a) + NormalCDF(b) - 1.0 + BivariateNormalCDF(-1 * a, -1 * b, rho);
    else
    throw new Exception("This part of the code should never be reached.");

    }
    else if (a * b * rho >= 0.0)
    {
    double denum = Math.Sqrt(a * a - 2 * rho * a * b + b * b);
    double rho1 = ((rho * a - b) * Math.Sign(a)) / denum;
    double rho2 = ((rho * b - a) * Math.Sign(b)) / denum;
    double delta = (1.0 - Math.Sign(a) * Math.Sign(b)) / 4.0;
    return BivariateNormalCDF(a, 0.0, rho1) +
    BivariateNormalCDF(b, 0.0, rho2) - delta;
    }
    else
    throw new Exception("This part of the code should never be reached.");

    }

    public static double NormalCDF(double z)
    {
    if (z > 6.0) return 1;
    if (z < -6.0) return 0;

    double b1 = 0.31938153;
    double b2 = -0.356563782;
    double b3 = 1.781477937;
    double b4 = -1.821255978;
    double b5 = 1.330274429;
    double p = 0.2316419;
    double c2 = 0.3989423;

    double a = System.Math.Abs(z);
    double t = 1.0 / (1.0 + a * p);
    double b = c2 * System.Math.Exp((-1 * z) * (z / 2.0));
    double n = ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;
    n = 1.0 - b * n;
    if (z < 0.0) n = 1.0 - n;
    return n;
    }
    }
     
  3. thenmmm

    thenmmm

    10x, i'll take a look...basically that's what i have so far (the stirling's approximation gives a very small error and it's used only in the case when i have to calculate the CDF in the d1, d2 calculations in black-scholes - in other words...the error is even smaller (the funny variable names are meant...for a reason :).

    private void button1_Click(object sender, EventArgs e)
    {
    double ju = Convert.ToDouble(textBox1.Text); // price
    double pathcitup = Convert.ToDouble(textBox3.Text); //strike
    double jd = Convert.ToDouble(textBox4.Text); //riskfree
    double kooo = Convert.ToDouble(textBox5.Text); //volat.


    try
    {
    dateTimePicker1.MinDate = DateTime.Today;
    DateTime today = DateTime.Today.AddDays(0);
    DateTime varDate = dateTimePicker1.Value;
    System.TimeSpan diffResult = varDate - today;
    Double nigas = Convert.ToDouble(diffResult.TotalDays);
    Double tubi = nigas / 365;


    double d1 = (Math.Log(ju/pathcitup)+ (jd + Math.Pow(kooo, 2) / 2) * tubi) / (kooo * (Math.Sqrt(tubi)));


    double d2 = d1 - kooo * Math.Sqrt(tubi);
    double nanananananananayee = 3.141592;
    double a = (8 * (nanananananananayee - 3)) / (3 * nanananananananayee * (4 - nanananananananayee)); // (8*($pi-3))/(3*$pi*(4-$pi));
    double x = d1 / Math.Sqrt(2); // $x = $d1/Sqrt(2);
    double x2 = Math.Pow(x, 2); // $x2 = Pow($x,2);
    double term = (-x2 * ((4 / 3.141592 + a * x2) / (1 + a * x2))); // $term = (-$x2*((4/3.141592+$a*$x2)/(1+$a*$x2)));
    double calc = (x / Math.Abs(x)) * Math.Sqrt(1 - Math.Pow(2.71828, term)); // ($x/abs($x))*Sqrt(1-Pow( 2.71828,$term));
    double nd1 = 0.5 * (1 + calc);

    double xd1 = d2 / Math.Sqrt(2); // $xd1 = $d2/Sqrt(2);
    double x22 = Math.Pow(xd1, 2); //$x22 = Pow($xd1,2);
    double term1 = (-x22 * ((4 / 3.141592 + a * x22) / (1 + a * x22))); // $term1 = (-$x22*((4/3.141592+$a*$x22)/(1+$a*$x22)));
    double calc1 = (xd1 / Math.Abs(xd1)) * Math.Sqrt(1 - Math.Pow(2.71828, term1)); // ($xd1/abs($xd1))*Sqrt(1-Pow( 2.71828,$term1));
    double nd2 = 0.5 * (1 + calc1);// $nd2 = 0.5*(1+$calc1);

    double exp = -jd * tubi; // $exp = -$r*$t;
    double ko = (ju * nd1) - pathcitup * Math.Exp(exp) * nd2; // $pachanga = ($s*$nd1)-$k*exp($exp)*$nd2;
    string johny = ko.ToString();
    textBox2.Text = johny;
    double ko3 = pathcitup * Math.Pow(2.71828, -jd * tubi) - ju + ko; //$donadoni = $k*Pow(2.71828,-$r*$t) - $s + $pachanga;
    textBox6.Text = ko3.ToString();


    }
    catch
    {
    }


    ---------------------------------------------------------------

    anyways...here's a screenshot so far. Basically the most usefull thing in the code above is the fact that it calculates times automatically and tells you what day of the week will be.

    http://oi52.tinypic.com/15pk27q.jpg

    According to the Espen Haug's excell black-scholes approximation, the real price has a difference of .00...something compared to my calculation, now that i've tested it among 20-30 theoretical prices. Not a big difference at all, especially since real market prices barely correspond to BS.
     
  4. sjfan

    sjfan

    Huh? market price always equates to BS at the implied IV of the market price by definition. Also - CBOE has a javascript bs calculator (or they used to) that's useful to check answers against.

     
  5. thenmmm

    thenmmm

    An interesting read by the way:

    http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1012075&rec=1&srcabs=283308

    In addition, I recommend you watching the movie "trillion dollar bet", an entire doumentary dedicated to BS...with some obvious negativism expressed in it.

    Finally, it will be very rational to claim that option prices correspond only to the so called mysterious forces that no one has seen and no one can price: The forces of supply and demand.

    How do you think option prices were priced during the crash in 1987? Like...option traders were like: "Hm...the market just felt 20%...I've lost everything...now before I commit a suicide i will check in my calculator what Gamma is". Bollox....

    So...prices have some correlation with black-scholes, no doubt about it . But are just not made of black-scholes entirely.

    CBOE by the way have a usefull margin calculator...the most usefull tool on their website perhaps. But i am yet to find an online tool which prices also vanna, color etc. - which is what i am going to program later today.
     
  6. Craig66

    Craig66

  7. thenmmm

    thenmmm

  8. sjfan

    sjfan

    Not sure you understood my point = in the options market, you can either quote the option as a price or as an implied volatility (much like the bond market can quote either price of yield). The standard pricing equation is agreed upon as being black scholes in equity (other models are standard in bond options stuff) so that you are indifferent to either looking at BS implied vol or to market price. So the two matches by definition.

    So, quoted black scholes implied vol == quoted price. It's a quoting convention. Note: no one calibrates their BS with historical vol; BS is used to imply a market traded vol.

    Your supply and demand exerts on either the price or BS implied vol - same business.

    Maybe your shallow understanding of how all this is actually applied in the market explains why someone with two degrees as you do can't find a job in this market. (I now see that my attempt to help by providing code only rewarded me with your arrogance - in recommending that I watch a documentary). Care to return to the other thread on vol calculation where you completely failed at understanding how to scale vol?

     
  9. thenmmm

    thenmmm



    Actually no, because again i am not someone like you without job, being in the forum 24/7 - a wild guess is that you're the md here - otherwise it will be...somewhat weird. AGain...i gotta go now...
     
  10. How did a great & informative discussion about probabilities & related implementation details turn into this? :(
     
    #10     Nov 23, 2010