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

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:


\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)
About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s