Using Auto Differentiation to compute the Greeks.

AD_Greeks.hs

module Main where


import Control.Comonad.Cofree
import Control.Lens
import Numeric.AD
import Numeric.AD.Types
import Text.Printf (printf)


data OptionType = Call | Put
    deriving (Eq, Show)

data Greeks a = Greeks
            { delta :: a
            , rho   :: a
            , vega  :: a
            , theta :: a
            , gamma :: a
            , charm :: a
            }
    deriving (Show, Eq)

pdf :: Floating a => a -> a
pdf x = recip (sqrt (2*pi)) * exp (-(x^(2::Int)) / 2)

cdf :: Floating a => a -> a
cdf x = 1 - pdf x * sum (zipWith (*) [b1,b2,b3,b4,b5] ts)
  where
    t  = 1 / (1 + b0 * x)
    ts = map (t^) ([1..] :: [Int])
    b0 = 0.2316419
    b1 = 0.319381530
    b2 = -0.356563782
    b3 = 1.781477937
    b4 = -1.821255978
    b5 = 1.330274429


-- Black Scholes probabilities
d1 s k r v t = (log (s/k) + (r + v^(2::Int) / 2) * t) / (v * sqrt t)
d2 s k r v t = d1 s k r v t - v * sqrt t

blackScholesPrice :: (Floating a, Ord a) => a -> a -> a -> a -> a -> OptionType -> a
blackScholesPrice s k sigma r t optionType
    | s < 0 = error "Negative Price"
    | t <= 0 = max 0 (if optionType == Call then s-k else k-s)
    | otherwise = if optionType == Call
                  then cdf d1' * s - cdf d2' * k * exp (-r * t)
                  else cdf (-d2') * k * exp (-r * t) - cdf (-d1') * s
        where
          d1' = d1 s k r sigma t
          d2' = d2 s k r sigma t

bsp, bsd, bsg, bsv, bst, bsc :: (Floating a, Ord a) => a -> [a] -> a
bsp k [spot, r, vol, t] = blackScholesPrice spot k vol r t Call
bsd k [spot, r, vol, t] = blackScholesDelta spot k vol r t Call
bsg k [spot, r, vol, t] = pdf d1' / (spot * vol * sqrt t)
  where d1' = d1 spot k r vol t
bsv k [spot, r, vol, t] = blackScholesVega  spot k vol r t
bsr k [spot, r, vol, t] = k * t * exp(-r * t) * cdf (d2 spot k r vol t)
bst k [spot, r, vol, t] = -spot * pdf d1' * vol / (2*sqrt t) - r*k*exp (-r * t) * cdf d2'
  where
    d1' = d1 spot k r vol t
    d2' = d2 spot k r vol t
bsc k [spot, r, vol, t] = -pdf d1' * (2*r*t - d2'*vol*sqrt t) / (2*t*vol*sqrt t)
  where
    d1' = d1 spot k r vol t
    d2' = d2 spot k r vol t

blackScholesVega :: Floating a => a -> a -> a -> a -> a -> a
blackScholesVega s k sigma r t = s * pdf d1' * sqrt t
  where
    d1' = d1 s k r sigma t

blackScholesDelta :: Floating a => a -> a -> a -> a -> a -> OptionType -> a
blackScholesDelta s k sigma r t optionType =
    if optionType == Call
    then cdf d1'
    else -cdf (-d1')
  where
     d1' = d1 s k r sigma t


-- g is a Cofree [] a, which works a bit like a tree.
-- unwrapped cuts off the root, then you are left with a list of branches
-- to select from, ordered by the argument with respect to which the
-- derivative is taken.
getGreeks :: Num a => Cofree [] a -> Greeks a
getGreeks g = Greeks delta rho vega theta gamma charm
  where
    dS = unwrapped . element 0
    dR = unwrapped . element 1
    dV = unwrapped . element 2
    dT = unwrapped . element 3
    delta = g^.dS.extracted
    rho   = g^.dR.extracted
    vega  = g^.dV.extracted
    gamma = g^.dS.dS.extracted
    -- tau = T-t, so if t increases, then tau *decreases*. Since these are
    -- really derivatives w.r.t. tau, we flip the sign.
    theta = negate $ g^.dT.extracted
    charm = negate $ g^.dS.dT.extracted

tableHeader :: String -> String -> String -> String
tableHeader = printf "%-8s %11s %11s"
tableRow :: String -> (Double, Double) -> String
tableRow label (x, y) =
    printf "%-8s %11.5f  %10.4e" label x y

main :: IO ()
main = do
    let spot = 100.0
    let k = 100.0
    let r = 0.06
    let vol = 0.2
    let t = 1.0
    let args = [spot, r, vol, t]
    let g = getGreeks $ grads (bsp (AD $ lift k)) args
    let f a b = (a, (a-b) / b) -- Relative error
    let optionDescription = "%s\nSpot: %.0f   Strike: %.0f   Interest: %.2f   Vol: %.2f   Tenor(years): %.1f"
    putStrLn $ printf optionDescription (show Call) spot k r vol t
    putStrLn $ tableHeader "" "Value" "Rel Err"
    putStrLn $ tableRow "Price" (f (bsp k args) (bsp k args))
    putStrLn $ tableRow "Delta" (f (bsd k args) (delta g))
    putStrLn $ tableRow "Vega"  (f (bsv k args) (vega g))
    putStrLn $ tableRow "Theta" (f (bst k args) (theta g))
    putStrLn $ tableRow "Rho"   (f (bsr k args) (rho g))
    putStrLn $ tableRow "Gamma" (f (bsg k args) (gamma g))
    putStrLn $ tableRow "Charm" (f (bsc k args) (charm g))


-- Output
---------------------------------------------------------------------------
-- Call
-- Spot: 100   Strike: 100   Interest: 0.06   Vol: 0.20   Tenor(years): 1.0
--                Value     Rel Err
-- Price       10.98955    0.0000e0
-- Delta        0.65542  -8.5164e-6
-- Vega        36.82701   4.9243e-6
-- Theta       -6.95586  -2.2077e-6
-- Rho         54.55262  -1.0232e-5
-- Gamma        0.01841   2.1945e-5
-- Charm       -0.07365   2.8790e-5