We can solve our differential equation using the Implicit Euler method which is unconditionally stable. We can also take this opportunity to use the Vector Package rather than Arrays as it has a richer set of combinators and to tidy up the code to make the payoff explicit (thanks to suggestions by Ben Moseley).

First let’s get our imports arranged so the code isn’t too cluttered with qualified names.

```
{-# LANGUAGE RecordWildCards #-}
import Prelude hiding
( length
, scanl
, scanr
, zip
, tail
, init
, last
)
import qualified Data.Vector as V
import Data.Vector
( (!)
, generate
, length
, Vector
, prescanl
, scanl
, scanr
, zip
, tail
, init
, last
, unfoldr
)
```

Define a type synomym for the payoff function and put all the parameters for the trade in one place.

```
type Payoff = Double -> Double
data Config = Config {
r :: Double,
sigma :: Double,
t :: Double,
m :: Int, -- Space steps
n :: Int, -- Time steps
xMax :: Double,
deltaX :: Double,
deltaT :: Double
}
defaultConf :: Config
defaultConf = Config {
r = 0.05,
sigma = 0.20,
t = 3.0,
m = 80,
n = 800,
xMax = 150,
deltaX = xMax defaultConf / fromIntegral (m defaultConf),
deltaT = t defaultConf / fromIntegral (n defaultConf)
}
```

The usual(!) *Comonad* class:

```
class Comonad c where
coreturn :: c a -> a
(=>>) :: c a -> (c a -> b) -> c b
```

Instead of using an Array to represent a "time slice" we use a Vector. Notice we can use the Vector library functions to implement the cobind.

```
data PointedArray a = PointedArray Int (Vector a)
deriving Show
instance Functor PointedArray where
fmap f (PointedArray i x) = PointedArray i (fmap f x)
instance Comonad PointedArray where
coreturn (PointedArray i x) = x!i
(PointedArray i x) =>> f =
PointedArray i (V.map (f . flip PointedArray x) (generate (length x) id))
```

We start with our partial differential equation:

Again we can approximate this by a difference equation but in this case we approximate the time derivative by .

Re-arranging (and not forgetting we are going backwards in time so that *n* is earlier than *n* − 1), we obtain an implicit equation for *z*_{n} in terms *z*_{n − 1}:

where

Thus we have a matrix equation which we need to invert. Fortunately the matrix is tridiagonal so we can use the Thomas Algorithm to invert it.

Let’s represent our tridiagonal matrix as a Vector of triples

```
data TriDiMatrix a = TriDiMatrix (Vector (a, a, a))
deriving Show
```

Here’s our tridiagonal matrix. Note that we have encoded two of the boundary conditions in this (the case for *S*(*t*, 0) and for *S*(*t*, *x*) as *x* → ∞).

```
triDi =
TriDiMatrix $ unfoldr h 0
where
h :: Int -> Maybe ((Double, Double, Double), Int)
h k | k == m' + 1 = Nothing
h k | k == m' = Just (( 0.0, 1.0, 0.0), k + 1)
h k | k == 0 = Just (( 0.0, 1.0, 0.0), k + 1)
h k = Just ((afn k, bfn k, cfn k), k + 1)
afn i = deltaT' * (r' * (fromIntegral i) - sigma'^2 * (fromIntegral i)^2) / 2
bfn i = 1 + deltaT' * (r' + sigma'^2 * (fromIntegral i)^2)
cfn i = negate $ deltaT' * (r' * (fromIntegral i) + sigma'^2 * (fromIntegral i)^2) / 2
m' = m defaultConf
deltaT' = deltaT defaultConf
r' = r defaultConf
sigma' = sigma defaultConf
```

Next we implement the Thomas Algorithm. *solveTriDi* takes a (tridiagonal) matrix and a vector and returns the application of the inverse of the matrix to the vector.

```
solveTriDi :: TriDiMatrix Double -> Vector Double -> Vector Double
solveTriDi (TriDiMatrix m) d = z
where
-- Forward sweep
c' = prescanl f (let (a1, b1, c1) = m!0 in c1 / b1) (tail m)
where f ci' (ai, bi, ci) = ci / (bi - ci'*ai)
d' = scanl g d1' $ zip (tail m) (zip c' (tail d))
where g di' ((ai, bi, ci), (ci', di)) = (di - di'*ai) / (bi - ci'*ai)
d1 = d!0
(_, b1, _) = m!0
d1' = d1 / b1
-- Backward substitution
z = scanr h xn $ zip c' (init d')
where
h (ci', di') xi1 = di' - ci'*xi1
xn = last d'
```

Now set up the boundary condition for our equation which is now a function of the payoff.

```
pricePArrAtT :: Config -> Payoff -> Int -> PointedArray Double
pricePArrAtT cfg@(Config {..}) fn i =
PointedArray i (V.fromList [ fn (deltaX * (fromIntegral j)) | j <- [0..m] ])
```

Finally we can price our call option without having to worry about the stability of the solution.

```
callPrices :: Int -> [Double]
callPrices i = map coreturn $ iterate (=>> oneStep) $ pricePArrAtT defaultConf fn i
where
oneStep (PointedArray j x) = (solveTriDi triDi x)!j
fn x = max 0 (x - k) -- A call
k = 50.0 -- The strike
callPrice i = (callPrices i)!!(n defaultConf - 1)
```