Monte Carlo options pricing in Haskell

optionModel.hs

{-# OPTIONS_GHC -funbox-strict-fields #-}
{-# OPTIONS_GHC -fexcess-precision #-}
{-# OPTIONS_GHC -Odph #-}
{-# OPTIONS_GHC -O2 #-}
-- {-# LANGUAGE BangPatterns    #-}
-- {-# LANGUAGE DoAndIfThenElse #-}
{-# LANGUAGE RankNTypes      #-}

-- GistID: 8970316

{- This has some theoretical issues, the largest of which is the random
 - number generation scheme. We should really generate a huge vector of
 - standard normals up front and then partition it and feed the batches to
 - our simulator. Reseeding the generator inside the simulator seems likely
 - to ruin it.
 -}

module Main where

import           Control.Monad.Identity
import qualified Control.Monad.ST                as ST
import qualified Data.Array.Repa                 as R
import qualified Data.Array.Repa.Eval            as Re
import           Data.Maybe                      (fromJust)
import qualified Data.Vector.Unboxed             as U
import qualified Data.Vector.Generic             as V
import           Statistics.Distribution         (cumulative)
import           Statistics.Distribution.Normal  (normalDistr)
import           System.Clock
import           System.Console.GetOpt
import           System.Environment              (getArgs)
import           System.Exit
import           System.IO
import           System.Random.MWC
import           System.Random.MWC.Distributions
import           Text.Printf                     (PrintfArg, printf)

-- import Control.Monad.Mersenne.Random as MR
-- import System.Random.Mersenne.Pure64


exampleStock, sx :: Stock
exampleStock = Stock 100 0.2 0.06 (1/365)
sx = exampleStock

exampleVanilla, vx :: VanillaOption
exampleVanilla = Vanilla 0 1 99 (exampleStock {tick = 1})
vx = exampleVanilla

exampleBarrier, bx :: BarrierOption
exampleBarrier = Barrier exampleVanilla 101 0
bx = exampleBarrier

exampleAsian, ax :: AsianOption
exampleAsian = Asian (Vanilla 0 1 99 exampleStock) Arithmetic 365
ax = exampleAsian


programDescription :: String
programDescription = "Price options with Monte Carlo methods."

data Options   = O {optSamples :: Int
                   ,optSeed    :: Int
                   ,optVerbose :: Bool
                   ,optHelp    :: Bool
                   } deriving Show

defaultOptions :: Options
defaultOptions = O { optSamples = 30000
                   , optSeed    = 0xDEADBEEF
                   , optVerbose = False
                   , optHelp    = False
                   }


main :: IO ()
main = do
    options <- getArgs >>= parseArgs
    when (optVerbose options) $ hPrint stderr options

    -- let theOption = vx {vanillaUnderlying = sx}
    -- let theControl = vx {vanillaUnderlying = sx}
    -- let theControl = asVanilla bx
    let theOption = ax
    let theControl = ax {asianAveraging = Geometric}
    let hexSeed = optSeed options
    let n = optSamples options


    printf "%.4f Analytical (if relevant)\n" (solution theOption)
    (wmc, cpumc, cimc)    <- timeit $ return (confidenceInterval . normalMC theOption n $ hexSeed)
    (wac, cpuac, ciac)    <- timeit $ return (confidenceInterval . antitheticMC theOption n $ hexSeed)
    (wcc, cpucc, cicc)    <- timeit $ return (confidenceInterval . controlMC normalMC theOption theControl n $ hexSeed)
    (wcac, cpucac, cicac) <- timeit $ return (confidenceInterval . controlMC antitheticMC theOption theControl n $ hexSeed)
    printf "%fs %fs %s\n" wmc cpumc cimc
    printf "%fs %fs %s\n" wac cpuac ciac
    printf "%fs %fs %s\n" wcc cpucc cicc
    printf "%fs %fs %s\n" wcac cpucac cicac


normalMC :: forall a. Option a => a -> Int -> Int -> MonteResult Double
normalMC theOption = monte (payoff theOption)
    where normals = genNormals . round . recip . tick . underlying
          payoff o = realize o . normals o -- Seed -> Option Payoff


antitheticMC :: forall a. Option a => a -> Int -> Int -> MonteResult Double
antitheticMC theOption = monte (payoff theOption) . (`div` 2)
    where normals = genNormals . round . recip . tick . underlying
          payoff o seed = let positives = normals o seed
                              negatives = V.map negate positives
                          in (realize o positives + realize o negatives) / 2


controlMC :: (Option a, Option b, Solvable b) =>
    (forall c. Option c => c -> Int -> Int -> MonteResult Double) ->
        a -> b -> Int -> Int -> MonteResult Double
controlMC doMC theOption theControl n hexSeed = result
  where
    result = MonteResult controlledMean controlledVar controlledStderr Nothing
    n' = fromIntegral n
    MonteResult mx varx _ x = doMC theOption n hexSeed
    MonteResult my vary _ y = doMC theControl n hexSeed
    controlledMean = mx + beta * (my - ey)
      where
        ey = solution theControl
        beta  = -covarxy / vary
    controlledStderr = 1.96 * sqrt (controlledVar / n')
    controlledVar = varx - covarxy*covarxy / vary
    covarxy = V.sum (V.zipWith (*) xdiff ydiff) / (n' - 1)
    xdiff = V.map (subtract mx) (fromJust x)
    ydiff = V.map (subtract my) (fromJust y)



data MonteResult t = MonteResult { mean :: t
                                 , var :: t
                                 , serr :: t
                                 , sample :: Maybe (U.Vector t)
                                 }

confidenceInterval :: (PrintfArg t, Floating t, Ord t) => MonteResult t -> String
confidenceInterval (MonteResult m _ s _)
    | m < 1e-1  = printf "%.4e +- %.4e (%.4e, %.4e)" m s (m-s) (m+s)
    | otherwise = printf "%.4f +- %.4e (%.4f, %.4f)" m s (m-s) (m+s)


monte :: (Floating t, U.Unbox t, Re.Elt t) =>
             (Int -> t) -> Int -> Int -> MonteResult t
monte f n seed = runIdentity $ do
    let seeds = R.fromUnboxed (R.Z R.:. n) $ genInts n seed
    samples <- R.computeUnboxedP $ R.map f seeds
    ms <- R.sumAllP samples
    let m = ms / fromIntegral n -- Mean
    vs <- R.sumAllP . R.map (\x -> (x - m)^(2::Int)) $ samples
    let v = recip (fromIntegral (n-1)) * vs -- Variance
    let s = 1.96 * sqrt (v / fromIntegral n) -- Stderr @ 95%
    return $! MonteResult m v s (Just $ R.toUnboxed samples)


data Average = Arithmetic | Geometric
    deriving (Show, Eq)

data Stock = Stock { spot         :: {-# UNPACK #-} !Double
                   , volatility   :: {-# UNPACK #-} !Double
                   , riskFreeRate :: {-# UNPACK #-} !Double
                   , tick         :: {-# UNPACK #-} !Double
                   } deriving (Show)

data VanillaOption = Vanilla { vanillaInception  :: {-# UNPACK #-} !Double
                             , vanillaMaturity   :: {-# UNPACK #-} !Double
                             , vanillaStrike     :: {-# UNPACK #-} !Double
                             , vanillaUnderlying :: {-# UNPACK #-} !Stock
                             } deriving Show

data AsianOption = Asian { asianAsVanilla :: {-# UNPACK #-} !VanillaOption
                         , asianAveraging :: !Average
                         , asianPeriods   :: {-# UNPACK #-} !Int
                         } deriving (Show)

data BarrierOption = Barrier { barrierAsVanilla :: {-# UNPACK #-} !VanillaOption
                             , barrierCeiling   :: {-# UNPACK #-} !Double
                             , barrierFloor     :: {-# UNPACK #-} !Double
                             } deriving (Show)


class Brownian bm where
    drift      :: bm -> Double
    variance   :: bm -> Double
    resolution :: bm -> Double
    path       :: bm -> U.Vector Double -> U.Vector Double

instance Brownian Stock where
    drift            = riskFreeRate
    variance stock   = volatility stock^(2::Int)
    resolution       = tick
    path s normals   = prices
      where
        prices    = V.scanl' (*) (spot s) brownians
        brownians = V.map factor normals
        factor z  = exp ((r - 0.5 * sigma^(2::Int)) * dt + sigma * sqrt dt * z)
          where
            dt = tick s
            r = riskFreeRate s
            sigma = volatility s

class Option opt where
    asVanilla  :: opt -> VanillaOption
    realize    :: opt -> U.Vector Double -> Double
    inception  :: opt -> Double
    maturity   :: opt -> Double
    strike     :: opt -> Double
    underlying :: opt -> Stock
    periods    :: opt -> Int
    periods = const 1


instance Option VanillaOption where
    asVanilla    = id
    inception    = vanillaInception
    strike       = vanillaStrike
    maturity     = vanillaMaturity
    underlying   = vanillaUnderlying
    realize o ns = discount * max 0 (V.last prices - strike o)
      where
        prices = path (underlying o) ns
        discount = exp $ -rfr * maturity o
        rfr = riskFreeRate (underlying o)

instance Option BarrierOption where
    asVanilla    = barrierAsVanilla
    inception    = inception . asVanilla
    strike       = strike . asVanilla
    maturity     = maturity . asVanilla
    underlying   = underlying . asVanilla
    realize o ns = discount * payoff
      where
        payoff | V.any hit prices = 0
               | otherwise        = max 0 (V.last prices - strike o)
        prices = path (underlying o) ns
        hit x = x > barrierCeiling o || x < barrierFloor o
        discount = exp $ -rfr * maturity o
        rfr = riskFreeRate (underlying o)

instance Option AsianOption where
    asVanilla    = asianAsVanilla
    inception    = inception . asVanilla
    strike       = strike . asVanilla
    maturity     = maturity . asVanilla
    underlying   = underlying . asVanilla
    periods      = asianPeriods
    realize o ns = discount * max 0 (avg prices - strike o)
      where
        discount = exp $ -rfr * maturity o
          where rfr = riskFreeRate (underlying o)
        avg x = if asianAveraging o == Arithmetic
                then mean' x
                else exp . mean' $ V.map log x
                  where mean' = (/ fromIntegral (V.length x)) . V.sum
        prices = V.drop upToWindow $ path (underlying o) ns
        total = round $ (maturity o - inception o) / tick (underlying o)
        upToWindow = total - periods o



class Solvable a where
    solution :: a -> Double


instance Solvable VanillaOption where
    solution o = s0 * normalCDF d1 - exp (-r * t) * k * normalCDF d2
      where
        d1 = (log (s0/k) + (r + sigma^(2 :: Int) / 2) * t) / (sigma * sqrt t)
        sigma = volatility stock
        stock = underlying o
        s0 = spot stock
        k = strike o
        r = riskFreeRate stock
        t = maturity o
        d2 = d1 - sigma * sqrt t
        normalCDF = cumulative (normalDistr 0 1)


instance Solvable AsianOption where
    solution o = exp ((rho-r)*t) * euroPrice
      where
        euroPrice = solution Vanilla { vanillaStrike = strike o
                                     , vanillaMaturity = maturity o
                                     , vanillaInception = inception o
                                     , vanillaUnderlying = stock' }
        s0 = spot stock
        stock = underlying o
        stock' = Stock { spot = s0
                   , volatility = sigmaZ
                   , riskFreeRate = rho
                   , tick = tick stock }
        rho = ((2*r - sigma^(2::Int)) * (m'+1)/(2*m') + sigmaZ^(2::Int)) / 2
        sigmaZ = sigma * sqrt (((m' + 1) * (2 * m' + 1)) / (6 * m' * m'))
        m' = fromIntegral $ periods o
        sigma = volatility stock
        r = riskFreeRate stock
        t = maturity o



---------------- BEGIN BORING STUFF ---------------------


opts :: [OptDescr (Options -> Options)]
opts = [ Option "n" ["samples"]
            (ReqArg (\s opts' -> opts' {optSamples = read s}) "10000")
            "Number of Monte Carlo Samples"
       , Option "s" ["seed"]
            (ReqArg (\s opts' -> opts' {optSeed = read s}) "0xDEADBEEF")
            "PRNG Seed"
       , Option "v" ["Verbose"]
            (NoArg (\opts' -> opts' {optVerbose = True}))
            "Verbose output"
       , Option "h" ["help"]
            (NoArg (\opts' -> opts' {optHelp = True}))
            "Print this help message"
       ]

parseArgs :: [String] -> IO Options
parseArgs argv =
    case getOpt Permute opts argv of
    (args, _nonargs, []) -> do
        let opts' = foldl (flip id) defaultOptions args
        if optHelp opts'
        then hPutStrLn stderr helpMessage  >> exitSuccess
        else return opts'
    (_, _, errs) -> do
        mapM_ (hPutStrLn stderr) (errs ++ [helpMessage])
        exitWith $ ExitFailure 1
    where
      helpMessage = usageInfo programDescription opts


genNormals :: Int -> Int -> U.Vector Double
genNormals = genNormals'

genInts :: Int -> Int -> U.Vector Int
genInts = genInts'


-- For using Mersenne Twister
-- genNormals'' :: Int -> Int -> U.Vector Double
-- genNormals'' n seed = MR.evalRandom (V.replicateM n computation) gen
    -- where
        -- computation = getDouble
        -- gen = pureMT (fromIntegral seed)

-- genInts'' :: Int -> Int -> U.Vector Int
-- genInts'' n seed = MR.evalRandom (V.replicateM n computation) gen
    -- where
        -- computation = getInt
        -- gen = pureMT (fromIntegral seed)


genNormals' :: Int -> Int -> U.Vector Double
genNormals' n seed = ST.runST $ do
    gen <- initialize (U.singleton $ fromIntegral seed)
    let computation = standard gen
    V.replicateM n computation

genInts' :: Int -> Int -> U.Vector Int
genInts' n seed = ST.runST $ do
    gen <- initialize (U.singleton $ fromIntegral seed)
    uniformVector gen n

diffSpecToSec :: TimeSpec -> TimeSpec -> Double
diffSpecToSec t u = fromIntegral (sec t - sec u) + fromIntegral (nsec t - nsec u) / 10^(9::Int)

timeit :: IO a -> IO (Double, Double, a)
timeit action = do
    wallTic <- getTime Monotonic
    cpuTic  <- getTime ThreadCPUTime
    result <- action
    wallToc <- getTime Monotonic
    cpuToc  <- getTime ThreadCPUTime
    let wallTime = diffSpecToSec wallToc wallTic
    let cpuTime = diffSpecToSec cpuToc cpuTic
    return (wallTime, cpuTime, result)