Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Monte Carlo options pricing in Haskell
{-# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment