### Normal Distribution Best Approach

1

I'm trying to build a simple program to price call options using the black scholes formula http://en.wikipedia.org/wiki/Black%E2%80%93Scholes. I'm trying to figure our the best way to get probabilities from a normal distribution. For example if I were to do this by hand and I got the value of as d1=0.43 than I'd look up 0.43 in this table http://www.math.unb.ca/~knight/utility/NormTble.htm and get the value 0.6664.

I believe that there are no functions in c or objective-c to find the normal distribution. I'm also thinking about creating a 2 dimensional array and looping through until I find the desired value. Or maybe I can define 300 doubles with the corresponding value and loop through those until I get the appropriate result. Any thoughts on the best approach?

2012-04-04 02:07
by SNV7
Just for the record, a similar question about calculating the cumulative normal distribution in Objective-C has been asked before. @Jason S. has provided one way to do it, but you can also import a C++ library that has the normal CDF already implemented - Li-aung Yip 2012-04-04 04:04

5

You need to define what it is you are looking for more clearly. Based on what you posted, it appears you are looking for the cumulative distribution function or P(d < d1) where d1 is measured in standard deviations and d is a normal distribution: by your example, if d1 = 0.43 then P(d < d1) = 0.6664.

The function you want is called the error function `erf(x)` and there are some good approximations for it.

Apparently `erf(x)` is part of the standard `math.h` in C. (not sure about objective-c but I assume it probably contains it as well).

But `erf(x)` is not exactly the function you need. The general form P(d < d1) can be calculated from `erf(x)` in the following formula:

``````P(d<d1) = f(d1,sigma) = (erf(x/sigma/sqrt(2))+1)/2
``````

where sigma is the standard deviation. (in your case you can use sigma = 1.)

You can test this on Wolfram Alpha for example: f(0.43,1) = (erf(0.43/sqrt(2))+1)/2 = 0.666402 which matches your table.

There are two other things that are important:

1. If you are looking for P(d < d1) where d1 is large (greater in absolute value than about 3.0 * sigma) then you should really be using the complementary error function `erfc(x) = 1-erf(x)` which tells you how close P(d < d1) is to 0 or 1 without running into numerical errors. For d1 < -3*sigma, P(d < d1) = (erf(d1/sigma/sqrt(2))+1)/2 = erfc(-d1/sigma/sqrt(2))/2, and for d1 > 3*sigma, P(d < d1) = (erf(d1/sigma/sqrt(2))+1)/2 = 1 - erfc(d1/sigma/sqrt(2))/2 -- but don't actually compute that; instead leave it as 1 - K where K = erfc(d1/sigma/sqrt(2))/2. For example, if d1 = 5*sigma, then P(d < d1) = 1 - 2.866516*10-7

2. If for example your programming environment doesn't have `erf(x)` built into the available libraries, you need a good approximation. (I thought I had an easy one to use but I can't find it and I think it was actually for the inverse error function). I found this 1969 article by W. J. Cody which gives a good approximation for erf(x) if |x| < 0.5, and it's better to use erf(x) = 1 - erfc(x) for |x| > 0.5. For example, let's say you want erf(0.2) ≈ 0.2227025892105 from Wolfram Alpha; Cody says evaluate with x * R(x2) where R is a rational function you can get from his table.

If I try this in Javascript (coefficients from Table II of the Cody paper):

`````` // use only for |x| <= 0.5
function erf1(x)
{
var y = x*x;
return x*(3.6767877 - 9.7970565e-2*y)/(3.2584593 + y);
}
``````

then I get `erf1(0.2) = 0.22270208866303123` which is pretty close, for a 1st-order rational function. Cody gives tables of coefficients for rational functions up to degree 4; here's degree 2:

`````` // use only for |x| <= 0.5
function erf2(x)
{
var y = x*x;
return x*(21.3853322378 + 1.72227577039*y + 0.316652890658*y*y)
/ (18.9522572415 + 7.8437457083*y + y*y);
}
``````

which gives you `erf2(0.2) = 0.22270258922638206` which is correct out to 10 decimal places. The Cody paper also gives you similar formulas for erfc(x) where |x| is between 0.5 and 4.0, and a third formula for erfc(x) where |x| > 4.0, and you can check your results with Wolfram Alpha or known erfc(x) tables for accuracy if you like.

Hope this helps!

2012-04-04 02:12
by Jason S
Thanks, yes I should have been a little more clear in the question - SNV7 2012-04-04 02:18
Excellent it works! - SNV7 2012-04-04 02:30