Wednesday, 23 December 2009

Surface plotting redux

Time has been short for the last fortnight, so haven't made as much progress as I'd like on this project. This short post covers one topic: adding labels to the plot from the last post. Here's an example of the result, using a 'pretty' function from F# for Visualization:

An Example Surface

To render text in a Viewport3D, this code follows roughly the same approach described by Eric Sink:

  1. Calculate the extents of the label in 3D space
  2. Create a rectangular "billboard" at the desired location
  3. "Texture Map" the desired text onto the billboard using a VisualBrush and a TextBlock

Positioning the text

To align the center of the label with the grid lines, the width of the text needs to be estimated. This can be done using GlyphTypeface, which provides the advance widths for individual characters:

module Wpf = 
  /// Attempts to retrieve a GlyphTypeface for the specified family,
  /// style and weight. If no matching typeface can be obtained then
  /// this raises an invalid argument exception.
  let GetTypeface (family:FontFamily)(style:FontStyle)(weight:FontWeight)=
    let t = Typeface(family, style, weight, FontStretches.Normal)
        
    let success, typeface = t.TryGetGlyphTypeface()
    if not success then
      invalidArg "family" "Cannot find a GlyphTypeface for the specified font family, style and weight"
            
    typeface

  /// Measures the length of a line of text in the specified typeface, 
  /// assuming a size of 1.0. Callers can scale the resulting width 
  /// according to the size they are using.
  let MeasureLine (text:string) (f:GlyphTypeface) = 
    Seq.fold (fun w (c:char) -> w + f.AdvanceWidths.[f.CharacterToGlyphMap.[int(c)]]) 0.0 text

Creating the billboard

The position of the label in three dimensional space is done using its center point and two vectors - one that points along the line of text and one that points from the bottom to the top of the line. Together the two vectors define the plane on which the text is rendered.

Texture coordinates are assigned to vertices, so in order to allow the text to be viewed the "right way around" from above and underneath we need to create two rectangles on which to display the text, with distinct vertices.

The extents of the rectangle are calculated by moving half of the text width away from the specified center point along the 'along' axis and half of the height on the 'up' axis. Two rectangular billboards are then created and painted using a visual brush containing a text block.

/// Creates a label oriented on a plane defined by the center point of the
/// text and two unit vectors defining its orientation and plane.
let CreateLabel (text:string) (center:Point3D) (up:Vector3D) (along:Vector3D) =           
    let height      = 0.035 // hard-coded font height, 
    let width       = height * Wpf.MeasureLine text this.typeface
    let bottomLeft  = center - width / 2.0 * along - height / 2.0 * up;
    let bottomRight = bottomLeft + (along * width);
    let topLeft     = bottomLeft + (up * height);
    let topRight    = topLeft + (along * width);

    let points = Point3DCollection()
    let rectangle = MeshGeometry3D()
    rectangle.Positions <- points

    // Shorthand for adding a point and associated UV coordinate
    let point p u v =
        points.Add(p)
        rectangle.TextureCoordinates.Add(Point(u, v))

    // Shorthand for creating a triangle from three point indices
    let triangle p0 p1 p2 = 
        rectangle.TriangleIndices.Add(p0)
        rectangle.TriangleIndices.Add(p1)
        rectangle.TriangleIndices.Add(p2)

    // Add the points and texture coordinates for the underside 
    point bottomLeft  0.0 1.0
    point topLeft     0.0 0.0
    point bottomRight 1.0 1.0
    point topRight    1.0 0.0
        
    // Add the points and texture coordinates for the topside
    point bottomLeft  1.0 1.0
    point topLeft     1.0 0.0
    point bottomRight 0.0 1.0
    point topRight    0.0 0.0

    // Define the triangles that make up the two sides of the rectangle
    triangle 0 3 1
    triangle 0 2 3
    triangle 4 5 7
    triangle 4 7 6        

    let block = TextBlock(Run(text))
    block.Foreground <- this.brush
    block.FontFamily <- this.family

    // Hint that the visual brush is not going to change, so that it is
    // not re-rendered on every frame.
    let brush = VisualBrush(block)
    RenderOptions.SetCachingHint(brush, CachingHint.Cache)

    let visual = ModelVisual3D()
    visual.Content <- GeometryModel3D(rectangle, DiffuseMaterial(brush))
    visual        

Performance

Having written the small amount of code described above, I fired up the application to be met with abysmally jittery rendering. Three small changes restored bearable performance:

Enabling caching for the VisualBrush

By default, the contents of the VisualBrush are rendered for every frame. WPF allows the application to hint that the Visual underlying the VisualBrush is not frequently changing and can be cached. If the hint is obeyed, then the Visual is rendered off-screen, reducing the cost to regular texture mapping.

This has an advantage over rendering the text to a bitmap yourself, since WPF will re-render the Visual at different resolutions as its projection onto the screen changes in size to ensure that the appearance doesn't significantly degrade.

To use this hint requires a single line, although the corresponding MSDN page describes how the results can be tuned:

let brush = VisualBrush(textblock)
RenderOptions.SetCachingHint(brush, CachingHint.Cache)

Disabling hit-testing

Leaving hit-testing enabled has a noticable impact on performance. Since it's not used in this application, it can be disabled on the viewport with:

viewport.IsHitTestVisible <- false

Disabling clipping

There isn't a significant amount of content outside the window and the Viewport3D fills the window, so clipping can be disabled for an additional speed up. The Maximize WPF 3D Performance page on MSDN indicates that the Viewport3D's clipping is slow. Disabling it is a single-line change:

viewport.ClipToBounds <- false

Updated code

The updated code can be downloaded from here.

Saturday, 5 December 2009

Plotting a surface using WPF 3D from F#

One of the features that drew me to F# was the ease of interoperability with code written in other CLR based languages, allowing me to easily re-use existing code with F# as a concise and expressive glue language to write the core of the application – something I’ve used IronPython for in the past.

Continuing with the option pricing theme, I’m going to put together an application to visualize the sensitivity of the option price to the various inputs in the pricing. Windows Presentation Foundation seemed to be an easy way of achieving this, so this post walks through the creation of a small F# application to plot a function of two variables as a 3D surface.

Getting started

Being familiar with 3D graphics, but never having used WPF from F# before (or indeed any other language) it took me a little while to get my bearings. After creating a new F# project the following steps were necessary before I could use get going:

  • Add references to WindowsBase, PresentationCore, PresentationFramework and System.Xml
  • Downloaded 3D Tools from and added a reference to that library

3DTools provides two pieces of functionality that will be used in this application – a trackball that allows the camera to be rotated by clicking and dragging in the viewport, and a mechanism by which to draw lines of constant screen-space width (ScreenSpaceLines3D).

Constructing an approximation of the surface

The API for displaying a surface will accept a set of points on the X-axis, a second set of points on the Z-axis and a function y = f(x,z) to calculate the Y-axis displacement for each (X,Z) point. For each (Xi, Zi) two triangles will be created in the square made up of (Xi, Zi), (Xi+1,Zi), (Xi+1,Zi+1) and (Xi,Zi+1), e.g.

The following code constructs a mesh as described above:

/// Evaluates the surface defined by f at each point in the cartesian 
/// product of xs and zs. xs and zs are required to be ordered in ascending
/// numerical order.
let EvaluateSurface (xs:float array) (zs:float array) (f: float -> float -> float) =
    if zs.Length < 2 then invalidArg "zs" "Cannot construct surface with fewer than two points on the z-axis"
    if xs.Length < 2 then invalidArg "xs" "Cannot construct surface with fewer than two points on the x-axis"

    // Find the smallest and largest value on each axis; the members of 
    // each axis array must be passed in ascending order.
    let xMin   = xs.[0]
    let xMax   = xs.[xs.Length-1]
    let xRange = xMax - xMin
    let zMin   = zs.[0]
    let zMax   = zs.[zs.Length-1]
    let zRange = zMax - zMin

    // Evaluate the function f for each pairing of points in xs and ys. Can
    // be done with List.iter but the imperative form generates simpler, 
    // more efficient IL
    let points  = Point3DCollection()
    let mutable yMin = Double.MaxValue
    let mutable yMax = Double.MinValue
    for x in xs do
        for z in zs do
            let y = f x z
            if y < yMin then yMin <- y elif y > yMax then yMax <- y
            points.Add(Point3D(x, y, z))
    
    // Reposition the points in object space so (yMax + yMin)/2 is 
    // positioned at y=0, and similarly for the X and Y axes.
    let yRange      = yMax - yMin
    let translation = Vector3D( -(xMax+xMin)/2.0, -(yMax+yMin)/2.0, -(zMax+zMin)/2.0 )
    let scale       = Vector3D( 1.0 / xRange, 1.0 / yRange, 1.0 / zRange )
    let texture     = PointCollection()
    for i in 0..points.Count-1 do
        let p = points.[i]
        points.[i] <- Point3D( (p.X + translation.X) * scale.X,
                               (p.Y + translation.Y) * scale.Y,
                               (p.Z + translation.Z) * scale.Z ) 
        texture.Add(Point((p.Y-yMin)/yRange, 0.0))
    
    // We now have |xs| * |zs| points. We need to form a triangle mesh 
    // over the top. For each index i, we construct two triangles from 
    // points (i, i+1, i+|xs|) and (i+1, i+|xs|, i+|xs|+1)
    let numSquares = xs.Length * (zs.Length-1)
    let triangles  = Int32Collection(numSquares*6)
    for i in 0..zs.Length-2 do
        let row = i * xs.Length
        for j in 0..xs.Length-2 do
            // For convenience of notation, give the indices of the four
            // corners of the square labels.
            let bl = row+j
            let br = bl+1     
            let tl = bl+xs.Length  
            let tr = tl+1      
            // Construct two triangles (bl, tl, br) and (tl, tr, br)
            triangles.Add(bl)
            triangles.Add(tl)
            triangles.Add(br)
            triangles.Add(tl)
            triangles.Add(tr)
            triangles.Add(br)
            
    let mesh = MeshGeometry3D()
    mesh.Positions          <- points
    mesh.TextureCoordinates <- texture
    mesh.TriangleIndices    <- triangles
    mesh

Adding contour lines

To make it easy to see the relative height of different points on the surface and to highlight the shape of the surface, it will be coloured according to its displacement on the Y-axis. An easy way of doing this is by assigning UV texture coordinates to each point based on the ratio of the Y coordinate to the Y-axis range - (Y-YMIN)/(YMAX-YMIN) – as seen in the code above. A GradientBrush can then be used to specify the colour for each Y value. Contour lines can be added by adding thin bands of a single colour to the gradient brush. The following code will construct such a brush:

/// Creates a gradient brush containing equally-sized stripes of each 
/// colour in the list, separated by a thin gray contour line. If the 
/// blend parameter is true then each band will blend smoothly from one
/// colour to the next. If the blend parameter is false then each band
/// will be a solid colour.
let CreateGradientBrush (stripes:int list) blend = 
    /// Utility method to convert from an integer xxRRGGBB representation 
    /// to a WPF Color
    let inline ColorFromRRGGBB hex =
        let r = byte((hex >>> 16) &&& 0xff)
        let g = byte((hex >>> 8) &&& 0xff)
        let b = byte(hex &&& 0xff)
        Color.FromRgb(r, g, b)

    let gradient = LinearGradientBrush()
    gradient.StartPoint <- Point(0.0, 0.5)
    gradient.EndPoint   <- Point(1.0, 0.5)
        
    // A small change in texture space used to create a visually sharp
    // edge between two adjacent colours that shouldn't appear blended
    // into each other.
    let epsilon    = 0.0000001 
    let lineWidth  = 0.003
    let lineColour = Color.FromRgb(0x66uy, 0x66uy, 0x66uy)
        
    let limit = float(stripes.Length)    
    for i in 0..stripes.Length-2 do
        let c   = ColorFromRRGGBB stripes.[i]
        let c'  = if blend then ColorFromRRGGBB stripes.[i+1] else c
        let bottom   = float(i) / limit
        let top      = (float(i+1) / limit) - epsilon
        let lineMin  = top - lineWidth
        let colorMax = lineMin - epsilon
        gradient.GradientStops.Add(GradientStop(c,  bottom))
        gradient.GradientStops.Add(GradientStop(c', colorMax))
        gradient.GradientStops.Add(GradientStop(lineColour, lineMin))
        gradient.GradientStops.Add(GradientStop(lineColour, top))

    // Handle the final colour separately so that there isn't a gray
    // contour line at the top - it just finishes in the final colour.
    let numStops   = gradient.GradientStops.Count
    let lastOffset = gradient.GradientStops.[numStops-1].Offset
    let topColour  = ColorFromRRGGBB stripes.[stripes.Length-1]
    gradient.GradientStops.Add(GradientStop(topColour, lastOffset+epsilon))
    gradient.GradientStops.Add(GradientStop(topColour, 1.0))
    gradient

Rendering the mesh

Having constructed the model and the texture that we are going to apply, then there are not many more steps required to render it on the screen:

1. Apply the gradient to the mesh

To apply the gradient to the mesh, we need to create a GeometryModel3D that binds together the mesh with the material. In order for both sides of the mesh to be rendered it is necessary to set both the front and the back material:

let mesh     = EvaluateSurface xs ys f
let gradient = CreateGradientBrush [0x62a4ce;0x79aed0;0x90b9d3;0xa7c3d6;0xbeced9] false
let model    = GeometryModel3D()
model.Geometry     <- mesh
model.Material     <- DiffuseMaterial(gradient)
model.BackMaterial <- model.Material

2. Define the scene

In order to see the surface, it is necessary to add some lighting to the scene. To apply a uniform light to all faces an AmbientLight can be used. The model and the light are bound together in a single Model3DGroup (because we don't need to move or transform them separately) and then wrapped in a ModelVisual3D for display.

let light = AmbientLight()
light.Color <- Colors.White

let group = Model3DGroup()
group.Children.Add( model )
group.Children.Add( light )

let visual = ModelVisual3D()
visual.Content <- group

3. Create the user interface

The final steps are to create a Viewport3D, add the ModelVisual3D to it, add a camera and place it in a Window to be displayed. Wrapping the Viewport3D in an instance of the 3DTools TrackballDecorator is all that is required to allow the user to rotate the camera around the origin:
// Create a camera through which to view the model
let camera = PerspectiveCamera()
camera.FieldOfView   <- 60.0
camera.Position      <- Point3D(-1., 1., 1.5)
camera.LookDirection <- Vector3D(1., -1. -1.5)

// Create a viewport to display the visual
let viewport = Viewport3D()
viewport.Camera <- camera
viewport.Children.Add( visual )

let trackball = TrackballDecorator()
trackball.Content <- viewport

let window = Window()
window.Width   <- 800.0
window.Height  <- 800.0
window.Title   <- "Surface plot"
window.Content <- trackball
window.Show()

Application().Run() |> ignore

The results

You can download the code here and do with it as you wish. If you make any enhancements it would be good if you could let me know. Here are a couple of examples of the output:
open System
open Psi.Charting
open System.Windows

open Utils

[<EntryPoint>]
[<STAThread>]
let main args = 
    // Display two windows - both plotting the same function but over
    // different ranges and using different colour schemes.
    let f     = (fun x y -> Math.Cos(x) + Math.Cos(y))
    let p1 = plot_surface (linspace -1.0 1.0 40) (linspace -1.0 1.0 40) f  
    let p2 = plot_surface (linspace -4.5 4.5 100) (linspace -4.5 4.5 100) f
    p2.SetColourScheme(AlternativeColourScheme, true)

    Application().Run() |> ignore
    0

Plot 1: Showing cos(x)+cos(y) over the range -1..1 on both axes

Plot 2: Showing cos(x)+cos(y) over the range -4.5 to 4.5 on both axes

Saturday, 28 November 2009

The Black-Scholes formula in F#

A 30-second introduction to options

I’m not going to go into much depth here, since this blog has a technical focus, but it doesn’t hurt to provide a bit of background and define a few terms.

There are many different types of options, but in the vanilla case, an option is a contract wherein one party sells (for an amount referred to as the Premium) another party the right to purchase (a Call option) or sell (a Put option) an asset on pre-agreed terms (for example, at a fixed price referred to as the Strike), at certain points in time. The option buyer can then choose to exercise this right if they wish, but they are not obliged to do so.

If the option can be exercised only on a single date it is referred to as a European option, if it can be exercised on any date it can be referred to as an American option, and if it can be exercised on a fixed list of dates then it is referred to as Bermudan. These names have nothing to do with where the options are traded.

The asset underlying an option can take many different forms: it may be a stock (an Equity Option), a currency (an FX option), another financial contract (e.g. an Interest Rate Future Option – an option on an Interest Rate Futures contract) or a more complex deal (IR Swaptions exercise into Interest Rate Swap, Credit Default Swaptions into CDS, etc). Consequently there are many different pricing methods in existence, catering for different features and peculiarities of specific option types. For simple options there are closed-form approximations, for more exotic options simulation techniques are necessary.

The payoff of a vanilla call option on the exercise date is Max(0, S – K), sometimes written as (S-K)+ – the excess of the asset price over the strike price, if positive, or zero. Why? Because the call option allows you to pay K to receive the asset: if the price of the asset is greater than K then you can profit to the tune of S-K by exercising the option and selling the asset. If the price of the asset is less than K then you won’t exercise the option and your payoff is 0.

For a vanilla put option the option purchaser has a right to sell at K. If the price S is less than K on the exercise date then the option buyer can profit to the tune of K-S by buying the asset for S and exercising the option to sell it at the higher price K.

The Black-Scholes model

In 1973, Fischer Black and Myron Scholes published a paper deriving a formula, which has since come to bear their names, that calculates a theoretical price for a European option in a certain model of the world. That model of the world makes a number of assumptions, including:
  • It is possible to trade continuously, in any quantity, in order to hedge and without incurring transaction costs,
  • there is no arbitrage opportunity,
  • cash can be borrowed and lent at a constant risk-free rate of interest,
  • the underlying stock may be sold short,
  • the underlying stock price follows a Geometric Brownian Motion with a known constant volatility and drift,
  • the underlying stock does not pay a dividend.

Over the years various people have developed extensions to the original model that relax one or more of the assumptions, giving rise to a family of related models. The book The Complete Guide to Option Pricing Formulas by Espen Gaarder Haug catalogues many of these variations.

The Black-Scholes formula

To start with, we'll look at the original Black-Scholes formula. In a Black-Scholes world, the formula for a Call option, with the payoff (S-K)+ described above, on a non-dividend paying stock is:
The Black-Scholes formula for a call option
Where:
  • N is the Standard Normal CDF, as described in the previous post
  • S is the spot (= current = time t) price of the underlying asset
  • K is the strike price of the option – the price payable to purchase the asset at maturity
  • r is the continuously compounded annual "risk free" interest rate
  • T is the time to maturity, expressed as a fraction of a year
And:
Parts of the Black-Scholes formula: d1 and d2
Recall that the payoff of the put option, (K-S)+, is the reverse of the call payoff (S-K)+; this reversal is mirrored in the Black-Scholes formula for a put option:
The Black-Scholes formula for a put option

An F# implementation

After all of that, the F# implementation is trivial and not that interesting:
module Options.BlackScholesMerton

open Distributions

let BlackScholes call S K T r v = 
    let sqrtT  = Math.Sqrt(T)
    let vSqrtT = v * sqrtT
    let d1     = ( Math.Log( S / K ) + ( r + v * v / 2.0 ) * T ) / vSqrtT
    let d2     = d1 - vSqrtT
    let KeRT   = K * Math.Exp(-r * T)
        
    match call with 
    | true  -> S * Normal.CDF(d1) - KeRT * Normal.CDF(d2)
    | false -> KeRT * Normal.CDF(-d2) - S * Normal.CDF(-d1)

As an example of how this can be used, let's look at the RBS (Royal Bank of Scotland) December Call with a strike of 32 on Euronext. This is an American style option, but since RBS isn’t in a position to pay a dividend right now it can be priced as a European option. As of the 27th of November the mid price of the stock was 34.7225. That gives us the first three parameters of the function.

The next parameter, T, is the time to the expiry of the contract expressed as a fraction of a year. The contract specification specifies that the contract expires at 18:30 on the Third Friday of the Expiry month – 18th of December in our case. For simplicity here we can calculate the time to maturity on the basis of a 365 day year and the actual number of days to maturity (we could alternatively count the minutes to today's close, the minutes on the final trading day and the number of whole days in between).

In the contexts I've worked in, the rate r is usually interpolated from a LIBOR/Swap curve. In this case for simplicity we could just interpolate between short term deposit rates such as those on the Financial Times' Markets pages which as of today would give a mid rate of around 0.59%.

The final parameter – the volatility of the underlying asset over the remaining lifetime of the option – has a large impact on the price of the option. Depending on what we were using the option pricing formula for, we could take a couple of different approaches.

We might state our expectation of volatility over the lifetime of the option – perhaps based on historical data, although there is no guarantee the same volatility will be realised again. In this case we would just plug in our expectation as the volatility parameter.

Alternatively, we could note that we have fixed all of the other parameters to the model and can observe the traded price of the option: the only unknown is the volatility. Consequently we can run the Black-Scholes formula in reverse to obtain what is referred to as an implied volatility – the volatility that equates the Black-Scholes model theoretical price with the observed option prices on the market. This will be the topic of another post, as it will make for some interesting code and give me a chance to experiment with some visualizations. For the time being let's assume that we’ve somehow arrived at the conclusion that the volatility will be 20% over the remaining lifetime of the option.

Putting all of that together, a simple example of the usage of the code in FSI is:

> let expiry = DateTime( 2009, 12, 18 );;
val expiry : DateTime = 18/12/2009 00:00:00

> let T = float((expiry - DateTime.Today).Days) / 365.0;;
val T : float = 0.05479452055

> let r = 0.0059;;
val r : float = 0.0059

> let v = 0.20;;
val v : float = 0.2

> BlackScholesMerton.BlackScholes true 34.7225 32.0 T r v;;
val it : float = 2.758012344

That price is lower than the observed price of around 3.75. To some extent this is explained by the 20% volatility estimate being too low, but there is likely to be difference between the theoretical price and the market price due to other factors such as (perhaps) low trading volume which will make it more difficult to hedge.

Unfortunately there's not much technical content in this post, but it does provide some building blocks that will be used in future posts. Finally, I'm sure I shouldn't need to say anything like this, but any use of the above information or code that you make you do at your own risk: I am not qualified to give advice regarding financial matters and you should do your own research if considering putting any money on the line.

Sunday, 22 November 2009

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#.

Experimenting with F#

The F# language has been around for a couple of years now but has been evolving gradually from a research project to a production-quality language. It already has some industry adoption, with projects in a few large banks being built using the language. With its inclusion in Visual Studio 2010 I expect that we’re going to see it being more widely used in the next few years. At first sight it seems a pleasant language to work with, but it’s been years since I’ve written even a moderately sized application in a functional language so it will be interesting to try.

I’m going to start out fairly slowly, to get a feel for how the code I write gets transformed into CIL, to understand how it’s implemented and so reassure myself that I’m not writing F# in a style that is going to compile down to badly performing IL.

In the first couple of posts I’m going to put together a simple option pricer using the closed-form Black-Scholes approach. This formula requires the standard Normal distribution’s PDF and CDF which, as a learning exercise, I’ll implement below.

For reference, I will be using the October 2009 CTP (1.9.7.8).

The Normal Probability Density function

The Standard Normal probability density function is:
So a straightfoward implementation, factoring out the constant portion of the formula, is:
namespace Distributions

open System

module Normal = 
    /// Internal constant used in the evaluation of the PDF. 
    let private OneOverSqrt2PI = 1.0/Math.Sqrt(2.0 * Math.PI)

    /// Evaluates the standard normal PDF at point x. 
    let n x = Math.Exp( -0.5*x*x ) * OneOverSqrt2PI
Firing up Reflector we see the following IL for the Normal.n method:
.method public static float64 n(float64 x) cil managed
{
    .maxstack 4
    L_0000: nop 
    L_0001: ldc.r8 -0.5
    L_000a: ldarg.0 
    L_000b: mul 
    L_000c: ldarg.0 
    L_000d: mul 
    L_000e: call float64 [mscorlib]System.Math::Exp(float64)
    L_0013: call float64 Distributions.Normal::get_OneOverSqrt2PI()
    L_0018: mul 
    L_0019: ret 
}

That was a little surprising: I had assumed that the private OneOverSqrt2PI would have been translated into a private static double in the class implementing the module. Instead it had been implemented as a private property, with the constant 1/sqrt(2*PI) being retrieved by a call to the property’s getter.

A further surprise was that the double field backing the OneOverSqrt2PI property and the .cctor initialising it were not contained within the class implementing the module: they were placed in a private class named StartupCode$LibraryName.$Distributions (after the namespace). The property's getter retrieved the calculated value from this class. It seems that the F# compiler systematically generates a "startup class" for each namespace within the application.

Going one step further, it’s interesting to see what the JIT does with this: will the method retrieving the constant be inlined, or will we incur the cost of an extra method call each time we evaluate the PDF?

By default Visual Studio will suppress the JIT optimisations when the code is being debugged, which would disable the inlining even if the JIT would normally inline the method. To run with optimisation enabled, it is first necessary to go to Tools -> Options -> debugging and untick Suppress JIT optimization on module load and Enable just my code.

Having done that and hit a breakpoint at Normal.n, it's possibly to get a nice disassembly of the optimized JIT'd code using SOS.dll. To do this type .load sos in the immediate debug window, then use !u to disassemble the code at the instruction pointer EIP (copied from the registers debug window):

push  ebp
mov   ebp,esp
sub   esp,8
fld   qword ptr [ebp+8]
fld   st(0)
fmul  dword ptr ds:[01519424h]
fmulp st(1),st
sub   esp,8
fstp  qword ptr [esp]
call  675D746F (System.Math.Exp(Double))
fstp  qword ptr [ebp-8]
call  dword ptr ds:[00143828h] (Distributions.Normal.get_OneOverSqrt2PI())
fld   qword ptr [ebp-8]
fmulp st(1),st
mov   esp,ebp
pop   ebp
ret   8

In the disassembly above, it’s can be seen that the retrieval of the constant OneOverSqrt2PI is not inlined. This will still run pretty quickly, but is there any way we can improve on this?

The Literal attribute?

Being an F# novice, the first possible solution that jumped out at me when skimming the specification was the [<Literal>] attribute. The documentation states that a let binding annotated with [<Literal>] would be treated as if it were a literal at the point of use. Trying it out, when the let is annotated with [<Literal>] its value is inlined in the generated IL. Unfortunately to use literal the value must be a compile-time constant. Changing the F# source to:

namespace Distributions

open System

module Normal = 
    /// Internal constant used in the evaluation of the PDF: equal to 
    /// 1.0/sqrt(2 * PI)
    [<Literal>]
    let private OneOverSqrt2PI = 0.398942280401433

    /// Evaluates the standard normal PDF at point x. 
    let n x = Math.Exp( -0.5*x*x ) * OneOverSqrt2PI
results in the following IL, where it can be seen that the constant has been inlined at compile time:
.method public static float64 n(float64 x) cil managed
{
    .maxstack 4
    L_0000: nop 
    L_0001: ldc.r8 -0.5
    L_000a: ldarg.0 
    L_000b: mul 
    L_000c: ldarg.0 
    L_000d: mul 
    L_000e: call float64 [mscorlib]System.Math::Exp(float64)
    L_0013: ldc.r8 0.398942280401433
    L_001c: mul 
    L_001d: ret 
}

However, when the value of the let binding is a compile-time constant, then it it inlined even without the [<Literal>] attribute. Reading Chris Smith's article, F# Zen - The Literal Attribute, it appears that the purpose of the literal attribute is to cause the binding to be interpreted as its value in constructs where it would otherwise be considered to be a new local binding of the same name - for example, to use it as a pattern to match against, not a name to bind the value to.

Wrapping up...

I had expected to cover both the PDF and the CDF in a single post, but having spent this many column-inches on just the simple PDF I'll save looking at the CDF for another post.

This post leaves me with a few questions still open:

  • Why are the Module-level constant values are held in a per-namespace StartupCode class? (scope?)
  • Is there an alternative to using a compile-time constant that avoids a non-inlined method call?
  • Why get_OneOverSqrt2PI wasn't inlined - it would seem a good candidate given it simply returns a readonly field
I'm sure some of the answers will become apparent when I have the time to finish digesting the F# specification, but if anyone reading this post knows the answer to any of the above I'd be interested to know.