cumulative normal distribution in Python

Question

Hi,

I am looking for a function in Numpy or Scipy (or any rigorous Python library) that will give me the cumulative normal distribution function in Python. This is for Black Scholes option pricing. I can program it myself but it is a (decent) approximation and I'd like to test if there's something (even) better as as I still get a few decimals of error which I'd like to reduce?

Answer 1

Adapted from here http://mail.python.org/pipermail/python-list/2000-June/039873.html

from math import *
def erfcc(x):
    """Complementary error function."""
    z = abs(x)
    t = 1. / (1. + 0.5*z)
    r = t * exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
     t*(.09678418+t*(-.18628806+t*(.27886807+
     t*(-1.13520398+t*(1.48851587+t*(-.82215223+
     t*.17087277)))))))))
    if (x >= 0.):
     return r
    else:
     return 2. - r

def ncdf(x):
    return 1. - 0.5*erfcc(x/(2**0.5))

Answer 2

Here's an example:

>>> from scipy.stats import norm
>>> norm.cdf(1.96)
array(0.97500210485177952)

If you need the inverse CDF:

>>>> norm.ppf(norm.cdf(1.96))
array(1.9599999999999991)

Answer 3

To build upon Unknown's example, the Python equivalent of the function normdist() implemented in a lot of libraries would be:

def normcdf(x, mu, sigma):
    t = x-mu;
    y = 0.5*erfcc(-t/(sigma*sqrt(2.0)));
    if y>1.0:
        y = 1.0;
    return y

def normpdf(x, mu, sigma):
    u = (x-mu)/abs(sigma)
    y = (1/(sqrt(2*pi)*abs(sigma)))*exp(-u*u/2)
    return y

def normdist(x, mu, sigma, f):
    if f:
        y = normcdf(x,mu,sigma)
    else:
        y = normpdf(x,mu,sigma)
    return y