Suppose we want to find the price of a European call option. Then we need to solve the Black-Scholes equation:





Although this particular equation can be solved explicitly, under more realistic assumptions we have to rely on numerical methods. We can approximate the partial differential equation by a difference equation (the minus sign on the left hand side is because we are stepping backwards in time):





Re-arranging we obtain:





We can implement this in Haskell using a comonad.

import Data.Array class Comonad c where coreturn :: c a -> a (=>>) :: c a -> (c a -> b) -> c b

We make each "time slice" of the differential equation an array with a distinguished element.

data PointedArray i a = PointedArray i ( Array i a) deriving Show

We make this into a functor by using the functor instance of the underlying array.

instance Ix i => Functor ( PointedArray i) where fmap f ( PointedArray i a) = PointedArray i ( fmap f a)

An array with a distinguished element is a comonad in which the cobind updates each element of the array simultaneously.

instance Ix i => Comonad ( PointedArray i) where coreturn ( PointedArray i a) = a ! i ( PointedArray i a) => > f = PointedArray i (listArray (bounds a) ( map (f . flip PointedArray a) ( range $ bounds a)))

Now let’s set up the parameters for our equation.

The interest rate — 5% seems a bit high these days r = 0.05

The volatility — 20% seems a bit low these days sigma = 0.2

The strike k = 50.0

Tthe time horizon in years t = 3.0

The granularity of the space component of the grid — 3 times the strike should be good enough m = 80 xMax = 150 deltaX = xMax / ( fromIntegral m)

The granularity of time component of the grid. A necessary condition for the explicit Euler method (which we are using) to be stable is

n = 800 deltaT = t / ( fromIntegral n)

Now we can encode the backwards step of the difference equation and two of the boundary conditions:

f ( PointedArray j x) | j == 0 = 0.0 f ( PointedArray j x) | j == m = xMax - k f ( PointedArray j x) = a * x ! (j -1 ) + b * x ! j + c * x ! (j +1 ) where a = deltaT * (sigma ^ 2 * ( fromIntegral j) ^ 2 - r * ( fromIntegral j)) / 2 b = 1 - deltaT * (r + sigma ^ 2 * ( fromIntegral j) ^ 2 ) c = deltaT * (sigma ^ 2 * ( fromIntegral j) ^ 2 + r * ( fromIntegral j)) / 2

Finally we can set the boundary condition at maturity time to be the payoff of a call.

priceAtT :: PointedArray Int Double priceAtT = PointedArray 0 (listArray ( 0 , m) [ max 0 (deltaX * ( fromIntegral j) - k) | j <- [ 0 .. m] ])

And now we can iterate backwards in time to find the price today (but only at points on the grid).

prices = iterate ( => > f) priceAtT

If we want the price of the option for a given stock price today the we can just pull out the required value.

pricesAtM m ( PointedArray _ x) = x ! m price m = last $ map (pricesAtM m) $ take n $ iterate ( => > f) priceAtT