Experimenting with F# (II) - The Normal CDF

Continuing on from the previous post, my next step is to evaluate the cumulative density function. There is no closed form formula, but there are several widely used approximations with varying degrees of accuracy. For this post I'll use Hart’s 1968 approximation, described by Graeme West in Better approximations to cumulative normal functions.

I've used this algorithm before in a C# application, so for side by side comparison here's a C# implementation.

public static double CDF(double x)
{
    bool   n = x < 0.0;
    double y = n ? -x : x;

    double c = 0.0;

    if (y <= 37.0)
    {
        double e = Math.Exp(-(y * y) / 2);

        if (y < 7.07106781186547)
        {
            double a = 0.0352624965998911 * y + 0.700383064443688;
            a = a * y + 6.37396220353165;
            a = a * y + 33.912866078383;
            a = a * y + 112.079291497871;
            a = a * y + 221.213596169931;
            a = a * y + 220.206867912376;

            double b = 0.0883883476483184 * y + 1.75566716318264;
            b = b * y + 16.064177579207;
            b = b * y + 86.7807322029461;
            b = b * y + 296.564248779674;
            b = b * y + 637.333633378831;
            b = b * y + 793.826512519948;
            b = b * y + 440.413735824752;
            c = e * a / b;
        }
        else
        {
            double a = y + 0.65;
            a = y + 4 / a;
            a = y + 3 / a;
            a = y + 2 / a;
            a = y + 1 / a;
            c = e / (a * 2.506628274631);
        }
    }

    return n ? c : 1.0 - c;
}

Pretty straightforward. With a few syntaxtical changes I could drop exactly the same code in to the F# source, but I'll re-write it to see how it looks in a more functional style and without the use of mutable variables:

/// Approximation of the Standard Normal CDF using Hart's 1968 approach. 
let CDF x =
    match abs(x) with 
    | y when y < 7.07106781186547 ->
        let a = ((((((0.0352624965998911 * y + 0.700383064443688) 
                                         * y + 6.37396220353165)
                                         * y + 33.912866078383)
                                         * y + 112.079291497871)
                                         * y + 221.213596169931)
                                         * y + 220.206867912376)

        let b = (((((((0.0883883476483184 * y + 1.75566716318264) 
                                          * y + 16.064177579207) 
                                          * y + 86.7807322029461) 
                                          * y + 296.564248779674) 
                                          * y + 637.333633378831) 
                                          * y + 793.826512519948) 
                                          * y + 440.413735824752)

        let r = Math.Exp(-0.5*y*y) * a / b;
        if x <= 0. then r else 1.0-r
        
    | y when y < 37. ->
        let a = (y+1.)/((y+2.)/((y+3.)/((y+4.)/(y+0.65))))
        let r = Math.Exp(-0.5*y*y) / (a * 2.506628274631)
        if x <= 0. then r else 1.0-r
        
    | _ (* Anything >= 37 *) -> 0.

Personally, I find the imperative version more readable, but that could be a sign that I'm more used to reading imperative code or that I'm not yet sufficiently familiar with F# to write pretty code. To some extent the layout of the F# code was dictated by the desire not to use mutable and the compiler's indentation requirements.

Compared to the previous post, this method threw up fewer surprises: the only strange thing that I noticed looking at the IL was that when the conditions were reversed and reordered (y >= 37., y >= 7, _) then the F# compiler decided to split the CDF function in two, generating an extra method CDF$cont@22: Double -> Double -> Unit -> Double that performed two of the matches. I may have done something stupid, but I would guess this is just a compiler wrinkle that will be ironed out by the time VS 2010 ships.

In the next post I'll use the Normal PDF and CDF implemented in the last two posts to implement some option pricing routines in F#.