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
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