Suppose we want to find the price of a European call option. Then we need to solve the BlackScholes 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):
Rearranging 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!(j1) + 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
Pingback: The Implicit Euler Method « Idontgetoutmuch’s Weblog
Pingback: Option Pricing Using Haskell Parallel Arrays « Idontgetoutmuch’s Weblog