Black and Scholes formula

Here is an implementation of the Black and Scholes formula for simple and digital European options in F#:

//pricing.fs file

#light

open System
open Options

// Abramowitz and Stegun formula for the
// standard normal cumulative distribution
let normalCDF x =
 let b1 =  0.319381530
 let b2 = -0.356563782
 let b3 =  1.781477937
 let b4 = -1.821255978
 let b5 =  1.330274429
 let p  =  0.2316419
 let c  =  0.39894228
 let RationalApproximation t =
  c * exp( -x * x / 2.0 ) *
  t *( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 )
 if (x >= 0.0)
 then 1.0 - RationalApproximation(1.0 / (1.0 + p * x))
 else RationalApproximation(1.0 / (1.0 - p * x))

//Black and Scholes model parameters
type BlackScholesParameters =
{ S0: float;     //current underlying price
  sigma: float;  //volatility
  r: float;}     //risk free rate


let blackSholesPricing myOption parameters =
  match myOption with
 | SimpleOption (payoffType, K, T) ->
     let d1 = (Math.Log(parameters.S0/K)
       + (parameters.r + 0.5*parameters.sigma**2.0)*T)
       / (parameters.sigma * Math.Sqrt(T))
     let d2 = d1 - parameters.sigma * Math.Sqrt(T)
     match payoffType with
    | Call -> parameters.S0 * normalCDF(d1)
           - K * Math.Exp(-parameters.r * T) * normalCDF(d2)
    | Put  -> K * Math.Exp(-parameters.r * T) * normalCDF(-d2)
           - parameters.S0 * normalCDF(-d1)

 | DigitalOption (payoffType, K, T) ->
     let d2 = (Math.Log(parameters.S0/K)
          + (parameters.r + 0.5*parameters.sigma**2.0)*T)
          / (parameters.sigma * Math.Sqrt(T))
          - parameters.sigma * Math.Sqrt(T)
     match payoffType with
    | Call -> Math.Exp(-parameters.r * T) * normalCDF(d2)
    | Put  -> Math.Exp(-parameters.r * T) * normalCDF(-d2)

We can now price an option using Black and Scholes formula, and test the convergence of our Monte Carlo pricer:

//exe.fs file

#light

open System
open Options
open Pricing
open MonteCarlo

let pricingParameters = {S0=100.0; sigma=0.19; r=0.06}

let K = 102.5 // Strike
let T = 0.3   // Maturity in years
let option = SimpleOption(Call, K, T)

//BS price
let priceBS = blackSholesPricing option pricingParameters

let monteCarlo = simulPayoffs option pricingParameters.S0
             pricingParameters.r T pricingParameters.sigma

//MC price (10,000 simulations)
let priceMC1 = mean (monteCarlo 10000)
//MC price (100,000 simulations)
let priceMC2 = mean (monteCarlo 100000)
//MC price (1,000,000 simulations)
let priceMC3 = mean (monteCarlo 1000000)

//MC relative errors
let error1 = (priceMC1 - priceBS) / priceBS
let error2 = (priceMC2 - priceBS) / priceBS
let error3 = (priceMC3 - priceBS) / priceBS

printf "Black and Scholes price: %f\n" priceBS
printf "Monte Carlo (10,000 simulations) price: %f,
error = %f\n" priceMC1 error1
printf "Monte Carlo (100,000 simulations) price: %f,
error = %f\n" priceMC2 error2
printf "Monte Carlo (1,000,000 simulations) price: %f,
error = %f\n" priceMC3 error3

Output:

Black and Scholes price: 3.836593
MC (10,000 simulations) price: 3.784829, error = -0.013492
MC (100,000 simulations) price: 3.883840, error = 0.012315
MC (1,000,000 simulations) price: 3.833737, error = -0.000744

Computing confidence intervals

The function simulPayoffs of our first program generates a list of n simulated payoffs for a given option. Then, to estimate the price of the option, we compute its mean using the following function:

//Compute the mean of a list
let mean list = Seq.fold (+) 0.0 list / float (Seq.length list)

In order to estimate the statistical error of this result, we define a function which computes the standard deviation of a list:

//Compute the estimated standard deviation of a list
let standardDeviation list =
  let m = mean list
  Seq.sumByFloat(fun x -> (x - m)**2.0) list
  / float (Seq.length list - 1)

And using the central limit theorem we can also define function which computes a confidence interval for our price estimation:

//Compute a 1-p confidence interval for p in (0,1)
let confidenceInterval list p =
  let m = mean list
  let sigma = standardDeviation list
  let quantile = 1.0-(p/2.0) |> NormalCDFInverse
  let n = Seq.length list
  let size = quantile*sigma/Math.Sqrt(float n)
 (m - size, m + size)

Let's test these functions on a simple European option:

let T = 0.3         //maturity
let strike = 102.5  //strike
let S0 = 100.0      //current underlying price
let r = 0.06        //risk free rate
let sigma = 0.19    //volatility
let n = 100000      //number of simulations

let call = SimpleOption(Call, strike, T)

let simul = simulPayoffs call S0 r T sigma n //simulated payoffs

let price = mean simul                   // estimated price
let sd = standardDeviation simul         // standard deviation
let ic95 = confidenceInterval simul 0.05 // 95% CI
let ic99 = confidenceInterval simul 0.01 // 99% CI


Console.WriteLine ("price: " + price.ToString())
Console.WriteLine ("standard deviation: " + sd.ToString())
Console.Write "95% confidence interval: "
print_any ic95
Console.Write "\n99% confidence interval: "
print_any ic99


Output:
price: 3.86677460654713
standard deviation: 39.7853613645219
95% confidence interval: (3.620132697, 4.113416516)
99% confidence interval: (3.542652267, 4.190896946)

Adding a new option type

Adding a new option type to this Monte Carlo pricer is very easy: For example, let's suppose that we want to price an option with two barriers which pays one if the price of the underlying is between these barriers, zero otherwise. We just need to add one line to option.fs to define a new type for two barriers options, and add a new case to our payoff function in payoff.fs:

//option.fs file
#light

type optionType = Call | Put

//oneBarrierOption: (call or put, strike, expiry)
type oneBarrierOption = optionType * float * float

//twoBarriersOption:
//   (call or put, lowerBarrier, upperBarrier, expiry)
type twoBarriersOption = optionType * float * float * float

type EuropeanOption = SimpleOption of oneBarrierOption
            | DigitalOption of oneBarrierOption
            | TwoBarriersDigitalOption of twoBarriersOption


//payoffs.fs file
#light
open System
open Options

let payoff myOption spot =
match myOption with

| SimpleOption (payoffType, strike, _) -> match payoffType with
| Call -> max(spot - strike) 0.0
| Put -> max(strike - spot) 0.0
                   
| DigitalOption (payoffType, strike, _) -> match payoffType with
| Call -> if(spot > strike) then 1.0 else 0.0
| Put -> if(spot < strike) then 1.0 else 0.0
               
| TwoBarriersDigitalOption (_, lower, upper, _) ->
   if(spot > lower && spot < upper)
   then 1.0
   else 0.0

We can now price this kind of option:

//main.fs file
#light

open System
open Options
open Payoffs
open Montecarlo

let T = 0.3             //maturity
let strikeInf = 99.5    //lower barrier
let strikeSup = 102.5   //upper barrier
let S0 = 100.0          //current underlying price
let r = 0.06            //risk free rate
let sigma = 0.19        //volatility
let n = 100000          //number of simulations

let myOption
= TwoBarriersDigitalOption(Call, strikeInf, strikeSup, T)

let price = mean (simulPayoffs myOption S0 r T sigma n)

print_any price

output:
0.1111708073

Monte Carlo simulation

The aim of this blog is to show how one could implement financial mathematics models in a functional programming language. The examples are given in F#, a new .NET language based on Ocaml.
In this first post, I will give an example of a simple Monte Carlo model for European options as we can find in the first chapters of Mark Joshi’s book, C++ design patterns and derivatives.

We first define the option types and the payoff functions for simple and digital options:

//option.fs file
#light
type optionType = Call | Put

//oneBarrierOption: (call or put, strike, expiry)
type oneBarrierOption = optionType * float * float

type EuropeanOption = SimpleOption of oneBarrierOption
          | DigitalOption of oneBarrierOption


//payoffs.fs file
#light
open System
open Options

let payoff myOption spot =
match myOption with

| SimpleOption (payoffType, strike, _) -> match payoffType with
 | Call -> max(spot - strike) 0.0
 | Put -> max(strike - spot) 0.0

| DigitalOption (payoffType, strike, _) -> match payoffType with
 | Call -> if(spot > strike) then 1.0 else 0.0
 | Put -> if(spot < strike) then 1.0 else 0.0

We define the type EuropeanOption as a discriminated union: it can be either a simple option or a digital option.
The oneBarrierOption type is simply defined as a tuple.

With this EuropeanOption type we define a payoff function which given a European option and a double, returns the payoff value for simple and digital options, both for calls and puts. As the EuropeanOption type is defined as a discriminated union, we can use pattern matching to return the good value for each kind of options. This pattern matching is type-safe: everything is statically typed.
The equivalent code in C++ would be much longer and more complicated to code: We would typically first define a base class option and derive from this class the different options classes. Then we would define a base class payoff and derive the payoffs from this class. And finally we would need to find a way to communicate between the different options and payoffs classes.

We can then write the functions for the pricing of these options using Monte Carlo technique:

//montecarlo.fs file
#light

open System
open Options
open Payoffs

//generate n random number U(0,1)
let gen = System.Random() in
let standardUniformSimulation n =
 [for i in 1..n -> gen.NextDouble()]

//Abramowitz and Stegun formula for the inverse
//standard normal cumulative distribution
let NormalCDFInverse p =
  let c0 = 2.515517
  let c1 = 0.802853
  let c2 = 0.010328
  let d1 = 1.432788
  let d2 = 0.189269
  let d3 = 0.001308
  let RationalApproximation t
    = t - ((c2*t + c1)*t + c0) / (((d3*t + d2)*t + d1)*t + 1.0)
  if (p >= 0.5)
  then RationalApproximation(sqrt(-2.0*log(1.0-p)))
  else -RationalApproximation(sqrt(-2.0*log(p)))

//generate n random number N(0,1)
let standardNormalSimulation n =
 standardUniformSimulation n |> List.map(NormalCDFInverse)

//simulate n geometric Brownian motions
let simulGeometricBrownian S0 r sigma T n =
standardNormalSimulation n
 |> List.map(fun x ->
  S0 * Math.Exp(r*T - 0.5*sigma*sigma*T
             + sigma*Math.Sqrt(T)*x))

//simulate n payoffs
let simulPayoffs myOption S0 r T sigma n =
 simulGeometricBrownian S0 r sigma T n
 |> List.map(payoff myOption)
 |> List.map(( * ) (Math.Exp((-r) * T)))

//Compute the mean of a list
let mean list = Seq.fold (+) 0.0 list / float (Seq.length list)

The first function uses the .NET pseudo-random number generator to simulate the uniform distribution.
We define a function standardNormalSimulation which simulates the normal distribution by applying the normal inverse cumulative distribution to a list of uniform random numbers.
We can then easily write two functions which simulate a list of geometric Brownian motions and payoffs.

Finally we can use these functions to price an option:

//main.fs file
#light

open System
open Options
open Payoffs
open Montecarlo

let expiry = 0.5
let myOption = SimpleOption(Call, 102.5, expiry)

let S0 = 100.0      //current underlying price
let r = 0.4         //risk free rate
let sigma = 0.11    //volatility
let n = 100000     //number of simulations

let simulation = simulPayoffs myOption S0 r expiry sigma n
let price = mean simulation

print_any price

output:
>>2.880850219