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 )