F#: a Monte Carlo Pricer for European Calls and Puts

Here we implement Monte Carlo simulation for pricing European Calls and Puts in a Black-Scholes world. Although we already have closed-form solutions, I think it's a good, simple but non-trivial example to show how Monte Carlo simulation can be nicely and concisely done in F#, thanks to the following functional programming techniques:

• List processing
• Higher-order function
• Function currying
• Function pipelining

There are a few things to note about the code example:

• Normal random numbers are generated using Box-Muller transform.
• The number of simulation runs is 10,000,000.
• The time function is used to measure CPU time and absolute time.
• The seed to initiate System.Random is deliberately fixed as 1 so that I can get the same series of normal random numbers every time I run the code. If you don't want this, please remember to remove the seed.
• The underlying dynamics and European call option details are as follows:
• S0 (spot price) = 50.0
• vol (annualized volatility) = 20%
• risk free rate = 1%
• strike = 40.0
• expiration (in years) = 0.25
• As a sanity check, the closed-form  theoretical price based on the parameters above is 10.11842.

After running the code on my 4-year-old laptop, I got the following output :

option_price = 10.118465
CPU time = 5709ms
Absolute time = 7115ms

It took about 7 seconds to run 10 million simulations. And the simulated price agrees with the theoretical price up to 4 decimal places.

Finally, here is the code:
///////////////////////////////////////////////////////////////////////////////

open System
open System.Diagnostics

let N = 10000000
let rnd = new Random(1)
let randomPairs = List.init N (fun i -> (rnd.NextDouble(), rnd.NextDouble()))
let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)
let randomNormals = randomPairs |> List.map box_muller_transform

let call_payout K S = max (S-K) 0.0

let European_Option_Price (payout:float->float) S0 r vol T =
    let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x)
    randomNormals |> List.map S_T |> List.map payout |> List.average 
                  |> (*) (exp(-r*T))

let S0 = 50.0
let K = 40.0
let r = 0.01
let vol = 0.2
let T = 0.25
let price() = 
    let price = European_Option_Price (call_payout K) S0 r vol T
    printfn "option_price = %f" price

let time f = 
    let proc = Process.GetCurrentProcess()
    let cpu_time_stamp = proc.TotalProcessorTime
    let timer = new Stopwatch()
    timer.Start()
    try
        f()
    finally
        let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds
        printfn "CPU time = %dms" (int64 cpu_time)
        printfn "Absolute time = %dms" timer.ElapsedMilliseconds

time price