CDF of the normal distribution in .net or c++ math libraries?
- Thread starter thenmmm
- Start date
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;
}
}
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;
}
}
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.
.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.
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.
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?
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?
How did a great & informative discussion about probabilities & related implementation details turn into this? 