Friday 17 October 2008

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

No comments: