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

import Data.Array

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