Solving a Partial Differential Equation Comonadically

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

\begin{aligned} \label {BlackScholeshDiffEqn0} \frac{\partial z}{\partial t} + \frac{1}{2}\sigma^2 x^2\frac{\partial^2 z}{\partial x^2} + rx\frac{\partial z}{\partial x} - rz = 0\end{aligned}

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

\begin{aligned} \frac{z^{n+1}_i - z^n_i}{-\Delta t} + \frac{1}{2} \sigma^2 (i\Delta x)^2\frac{z^n_{i+1} - 2z^n_i + z^n_{i-1}}{\Delta x^2} + r(i\Delta x)\frac{z^n_{i+1} - z^n_{i-1}}{2\Delta x} - rz^n_i = 0\end{aligned}

Re-arranging we obtain:

\begin{aligned} z^{n+1} = z^n_{i-1}\delta t(\frac{1}{2} \sigma^2i^2 - ri) + z^n_i(1 - \delta t(r + \delta^2i^2)) + z^n_{n+1}\delta t(\frac{1}{2} \sigma^2i^2 + ri)\end{aligned}

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
    \Delta T \le \frac{1}{\sigma^2n^2}

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

2 thoughts on “Solving a Partial Differential Equation Comonadically

  1. Pingback: The Implicit Euler Method « Idontgetoutmuch’s Weblog

  2. Pingback: Option Pricing Using Haskell Parallel Arrays « Idontgetoutmuch’s Weblog

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s