## The Implicit Euler Method

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

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)

coreturn (PointedArray i x) = x!i
(PointedArray i x) =>> f =
PointedArray i (V.map (f . flip PointedArray x) (generate (length x) id))

$\frac{\partial V}{\partial t} + \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + r S\frac{\partial V}{\partial S} - rV = 0$

Again we can approximate this by a difference equation but in this case we approximate the time derivative by $\frac{z^n_i - z^{n-1}_i}{\Delta t}$.

\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 (and not forgetting we are going backwards in time so that n is earlier than n − 1), we obtain an implicit equation for zn in terms zn − 1:

\begin{aligned} z^{n-1}_i = A_i z^n_{i-1} + B_i z^n_i + C_i z^n_{i+1}\end{aligned}

where

\begin{aligned} A_i &= \frac{\Delta t}{2}(ri - \sigma^2 i^2) \\ B_i &= 1 + \Delta t(r + \sigma^2 i^2) \\ C_i &= -\frac{\Delta t}{2}(ri + \sigma^2 i^2)\end{aligned}

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)