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

About these ads

This entry was posted on April 1, 2012 at 2:57 pm and is filed under category theory, Haskell, Probability. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

2 Responses to “Solving a Partial Differential Equation Comonadically”

The Implicit Euler Method « Idontgetoutmuch’s Weblog Says:
April 22, 2012 at 2:11 pm | Reply
[...] Idontgetoutmuch’s Weblog Mathematical Notes « Solving a Partial Differential Equation Comonadically [...]
Option Pricing Using Haskell Parallel Arrays « Idontgetoutmuch’s Weblog Says:
January 5, 2013 at 2:45 pm | Reply
[...] let us translate our pricer using the Explicit Euler Method to use [...]

Idontgetoutmuch’s Weblog

Solving a Partial Differential Equation Comonadically

Like this:

2 Responses to “Solving a Partial Differential Equation Comonadically”

Leave a Reply Cancel reply

Follow “Idontgetoutmuch's Weblog”