quantfin (empty) → 0.1.0.0
raw patch · 14 files changed
+863/−0 lines, 14 filesdep +basedep +containersdep +mersenne-random-pure64setup-changed
Dependencies added: base, containers, mersenne-random-pure64, mtl, random-fu, rvar, transformers, vector
Files
- LICENSE +30/−0
- Setup.hs +2/−0
- quantfin.cabal +87/−0
- src/Quant/ContingentClaim.hs +167/−0
- src/Quant/Math/Integration.hs +27/−0
- src/Quant/Models.hs +43/−0
- src/Quant/Models/Black.hs +57/−0
- src/Quant/Models/Dupire.hs +45/−0
- src/Quant/Models/Heston.hs +59/−0
- src/Quant/Models/Merton.hs +61/−0
- src/Quant/MonteCarlo.hs +160/−0
- src/Quant/Test.hs +63/−0
- src/Quant/VolSurf.hs +26/−0
- src/Quant/YieldCurve.hs +36/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2015, Timothy Dees + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Timothy Dees nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple +main = defaultMain
+ quantfin.cabal view
@@ -0,0 +1,87 @@+-- Initial quantfin.cabal generated by cabal init. For further +-- documentation, see http://haskell.org/cabal/users-guide/ + +-- The name of the package. +name: quantfin + +-- The package version. See the Haskell package versioning policy (PVP) +-- for standards guiding when and how versions should be incremented. +-- http://www.haskell.org/haskellwiki/Package_versioning_policy +-- PVP summary: +-+------- breaking API changes +-- | | +----- non-breaking API additions +-- | | | +--- code changes with no API change +version: 0.1.0.0 + +-- A short (one-line) description of the package. +synopsis: Quant finance library in pure Haskell. + +-- A longer description of the package. +-- description: + +-- URL for the project homepage or repository. +homepage: https://github.com/boundedvariation/quantfin + +-- The license under which the package is released. +license: BSD3 + +-- The file containing the license text. +license-file: LICENSE + +-- The package author(s). +author: Timothy Dees + +-- An email address to which users can send suggestions, bug reports, and +-- patches. +maintainer: timothy.dees@gmail.com + +-- A copyright notice. +-- copyright: + +category: Quant + +build-type: Simple + +-- Extra files to be distributed with the package, such as examples or a +-- README. +-- extra-source-files: + +-- Constraint on the version of Cabal needed to build this package. +cabal-version: >=1.10 + + +library + -- Modules exported by the library. + exposed-modules: Quant.YieldCurve + Quant.VolSurf + Quant.Math.Integration + Quant.Models + Quant.Models.Black + Quant.Models.Merton + Quant.Models.Dupire + Quant.Models.Heston + Quant.MonteCarlo + Quant.ContingentClaim + Quant.Test + + -- Modules included in this library but not exported. + -- other-modules: + + -- LANGUAGE extensions used by modules in this package. + -- other-extensions: + + -- Other library packages from which modules are imported. + build-depends: base >=4.7 && <4.8, + vector, + mtl, + transformers, + random-fu, + containers, + rvar, + mersenne-random-pure64 + + -- Directories containing source files. + hs-source-dirs: src + + -- Base language which the package is written in. + default-language: Haskell2010 +
+ src/Quant/ContingentClaim.hs view
@@ -0,0 +1,167 @@+module Quant.ContingentClaim ( + -- * Types for modeling contingent claims. + ContingentClaim + , ContingentClaim' (..) + , Observables (..) + , ContingentClaimBasket (..) + , OptionType (..) + , ccBasket + + -- * Options and option combinators + , vanillaOption + , binaryOption + , straddle + , arithmeticAsianOption + , geometricAsianOption + , callSpread + , putSpread + , forwardContract + , fixed + , multiplier + , short + , combine + , terminalOnly + , changeObservableFct + + -- * Utility functions + , obsNum + , obsHead + ) where + +import Data.List +import Data.Ord +import qualified Data.Vector.Unboxed as U + + +-- | 'ContingentClaim'' is the underlying type of contingent claims. +data ContingentClaim' = ContingentClaim' { + payoutTime :: Double -- ^ Payout time for cash flow + , collector :: [U.Vector Double] -> U.Vector Double -- ^ Function to collect observations and transform them into a cash flow. + , observations :: [( Double -- ^ Time of observation + , Observables -> U.Vector Double -- ^ Function to access specific observable. + , Double -> Double) ] -- ^ Function to run to transform observations. +} + +-- | 'ContingentClaim' is just a list of the underlying 'ContingentClaim''s. +type ContingentClaim = [ContingentClaim'] + +-- | Observables are the observables available in a Monte Carlo simulation. +--Most basic MCs will have one observables (Black-Scholes) whereas more +--complex ones will have multiple (i.e. Heston-Hull-White). +data Observables = Observables [U.Vector Double] deriving (Eq, Show) + +-- | ADT for Put or Calls +data OptionType = Put | Call deriving (Eq,Show) + +-- | Function to generate a vanilla put/call style payout. +vanillaPayout :: OptionType -- ^ Put or Call + -> Double -- ^ Strike + -> Double -- ^ Observable val + -> Double -- ^ Price +vanillaPayout pc strike x = case pc of + Put -> max (strike - x) 0 + Call -> max (x - strike) 0 + +-- | Function to generate a binary option payout. +binaryPayout :: OptionType -- ^ Put or call + -> Double -- ^ strike + -> Double -- ^ Payout amount if binary condition achieved + -> Double -- ^ observable level + -> Double -- ^ calculated payout +binaryPayout pc strike amount x = case pc of + Put -> if strike > x then amount else 0 + Call -> if x > strike then amount else 0 + +-- | Takes a maturity time and a function and generates a ContingentClaim +--dependent only on the terminal value of the observable. +terminalOnly :: Double -> (Double -> Double) -> ContingentClaim +terminalOnly t f = [ContingentClaim' t head [(t, obsHead, f)]] + +-- | Takes an OptionType, a strike, and a time to maturity and generates a vanilla option. +vanillaOption :: OptionType -> Double -> Double -> ContingentClaim +vanillaOption pc strike t = terminalOnly t $ vanillaPayout pc strike + +-- | Takes an OptionType, a strike, a payout amount and a time to +--maturity and generates a vanilla option. +binaryOption :: OptionType -> Double -> Double -> Double -> ContingentClaim +binaryOption pc strike amount t = terminalOnly t $ binaryPayout pc strike amount + +-- | Takes an OptionType, a strike, observation times, time to +--maturity and generates an arithmetic Asian option. +arithmeticAsianOption :: OptionType -> Double -> [Double] -> Double -> ContingentClaim +arithmeticAsianOption pc strike obsTimes t = [ContingentClaim' t f obs] + where obs = map (\x -> (x, obsHead, id)) obsTimes + f k = U.map (vanillaPayout pc strike . (/fromIntegral l)) + $ foldl1' (U.zipWith (+)) k + where l = length k + +-- | Takes an OptionType, a strike, observation times, time to +--maturity and generates a geometric Asian option. +geometricAsianOption :: OptionType -> Double -> [Double] -> Double -> ContingentClaim +geometricAsianOption pc strike obsTimes t = [ContingentClaim' t f obs] + where obs = map (\x -> (x, obsHead, id)) obsTimes + f k = U.map (vanillaPayout pc strike . (** (1/fromIntegral l))) + $ foldl1' (U.zipWith (*)) k + where l = length k + +-- | Scales up a contingent claim by a multiplier. +multiplier :: Double -> ContingentClaim -> ContingentClaim +multiplier notional cs = map f cs + where f c@(ContingentClaim' _ collFct _) = c { collector = U.map (*notional) . collFct } + +-- | Flips the signs in a contingent claim to make it a short position. +short :: ContingentClaim -> ContingentClaim +short = multiplier (-1) + +-- | Takes an amount and a time and generates a fixed cash flow. +fixed :: Double -> Double -> ContingentClaim +fixed amount t = terminalOnly t $ const amount + +-- | Takes a time to maturity and generates a forward contract. +forwardContract :: Double -> ContingentClaim +forwardContract t = terminalOnly t id + +-- | A call spread is a long position in a low-strike call +--and a short position in a high strike call. +callSpread :: Double -> Double -> Double -> ContingentClaim +callSpread lowStrike highStrike t = combine (vanillaOption Call lowStrike t) (short $ vanillaOption Call highStrike t) + +-- | A put spread is a long position in a high strike put +--and a short position in a low strike put. +putSpread :: Double -> Double -> Double -> ContingentClaim +putSpread lowStrike highStrike t = combine (vanillaOption Put highStrike t) (short $ vanillaOption Put lowStrike t) + +-- | A straddle is a put and a call with the same time to maturity / strike. +straddle :: Double -> Double -> ContingentClaim +straddle strike t = vanillaOption Put strike t ++ vanillaOption Call strike t + +-- | Just combines two contingent claims into one. +combine :: ContingentClaim -> ContingentClaim -> ContingentClaim +combine = (++) + +-- | Used to compile claims for the Monte Carlo engine. +data ContingentClaimBasket = ContingentClaimBasket ContingentClaim [Double] + +-- | Converts a 'ContingentClaim' into a 'ContingentClaimBasket' for use by the MC engine. +ccBasket :: ContingentClaim -> ContingentClaimBasket +ccBasket ccs = ContingentClaimBasket (sortBy (comparing payoutTime) ccs) monitorTimes + where monitorTimes = sort . nub $ concatMap (map fst3 . observations) ccs + +-- | Utility function to pull the head of a basket of observables. +obsHead :: Observables -> U.Vector Double +obsHead (Observables (x:_)) = x + +changeObservableFct' :: ContingentClaim' -> (Observables -> U.Vector Double) -> ContingentClaim' +changeObservableFct' c@(ContingentClaim' _ _ calcs) f = c { observations = map (\(t, _, g) -> (t, f, g)) calcs } + +-- | Offers the ability to change the function on the observable an option is based on. +--All options default to being based on the first observable. +changeObservableFct :: ContingentClaim -> (Observables -> U.Vector Double) -> ContingentClaim +changeObservableFct ccs f = map (`changeObservableFct'` f) ccs + +fst3 :: (a,b,c) -> a +fst3 (x, _, _) = x + +-- | Utility function for when the observable function is just '!!' +obsNum :: ContingentClaim -> Int -> ContingentClaim +obsNum ccs k = changeObservableFct ccs $ \(Observables x)-> x !! k
+ src/Quant/Math/Integration.hs view
@@ -0,0 +1,27 @@+module Quant.Math.Integration ( + Integrator + , midpoint + , trapezoid + , simpson + ) where + +-- | A function, a lower bound, an upper bound and returns the integrated value. +type Integrator = (Double -> Double) -> Double -> Double -> Double + +-- | Midpoint integration. +midpoint :: Int -> Integrator +midpoint intervals f lBound uBound = dx * sum (map f points) + where + dx = (uBound - lBound) / fromIntegral intervals + points = take intervals $ iterate (+dx) (lBound+dx/2) + +-- | Trapezoidal integration. +trapezoid :: Int -> Integrator +trapezoid intervals f lBound uBound = dx * ((f lBound + f uBound) / 2 + sum (map f points)) + where + dx = (uBound - lBound) / fromIntegral intervals + points = take (intervals-1) $ iterate (+dx) (lBound+dx) + +-- | Integration using Simpson's rule. +simpson :: Int -> Integrator +simpson intervals f l u =( 2 * midpoint intervals f l u + trapezoid intervals f l u ) / 3
+ src/Quant/Models.hs view
@@ -0,0 +1,43 @@+module Quant.Models ( + CharFunc(..) +) where + +import Data.Complex +import Quant.YieldCurve + + +{- | The 'CharFunc' class defines those +models which have closed-form characteristic +functions. + +Minimal complete definition: 'charFunc'. + +Still under construction. +-} +class CharFunc a where + -- | Creates a characteristic function for a model, without martingale adjustment. + charFunc :: CharFunc a => a -> Double -> Complex Double -> Complex Double + + -- | Calculates characteristic function given a forward generator and yield curve. + charFuncMart :: (CharFunc a, YieldCurve b) => a -> b -> Double -> Complex Double -> Complex Double + charFuncMart model fg t k = exp (i * r * k) * baseCF k + where + i = 0 :+ 1 + baseCF = charFunc model t + r = forward fg 0 t :+ 0 + + charFuncOption :: (CharFunc a, YieldCurve b, YieldCurve c) => + a -> b -> c -> ( (Double -> Double) -> Double) -> Double + -> Double -> Double -> Double + charFuncOption model fg yc intF strike tmat damp = intF f + where + f v' = realPart $ exp (i*v*k) * leftTerm * rightTerm + where + v = v' :+ 0 + damp' = damp :+ 0 + k = log strike :+ 0 + i = 0 :+ 1 + leftTerm = d / (damp' + i * v) / (damp'+i*v+(1:+0)) + rightTerm = cf $ v - i * (damp' + 1) + d = disc yc tmat :+ 0 + cf x = charFuncMart model fg tmat x
+ src/Quant/Models/Black.hs view
@@ -0,0 +1,57 @@+{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleInstances #-} + + +module Quant.Models.Black ( + Black (..) +) where + +import Quant.YieldCurve +import Data.Random +import Quant.Models +import Control.Monad.State +import Quant.MonteCarlo +import Quant.ContingentClaim +import qualified Data.Vector.Unboxed as U + +-- | 'Black' represents a Black-Scholes model. +data Black = forall a b . (YieldCurve a, YieldCurve b) => Black { + blackInit :: Double -- ^ Initial asset level. + , blackVol :: Double -- ^ Volatility. + , blackForwardGen :: a -- ^ 'YieldCurve' to generate forwards + , blackYieldCurve :: b } -- ^ 'YieldCurve' to handle discounting + +--instance CharFunc Black where + -- charFunc (Black s vol _ _) t k = exp + -- $ i*logs + negate i*vol'*vol'/2.0*t'*k-vol'*vol'*k*k/2.0*t' + -- where + -- i = 0 :+ 1 + -- t' = t :+ 0 + --vol' = vol :+ 0 + --logs = log s :+ 0 + +instance Discretize Black where + initialize (Black s _ _ _) trials = put (Observables [U.replicate trials s], 0) + + evolve' b@(Black _ vol _ _) t2 anti = do + (Observables (stateVec:_), t1) <- get + fwd <- forwardGen b t2 + let grwth = U.map (\x -> (x - vol*vol/2) * (t2-t1)) fwd + postVal <- U.forM (U.zip grwth stateVec) $ \ ( g , x ) -> do + resid <- lift stdNormal + if anti then + return $ x * exp (g - resid*vol) + else + return $ x * exp (g + resid*vol) + put (Observables [postVal], t2) + + discounter (Black _ _ _ dsc) t = do + size <- getTrials + return $ U.replicate size $ disc dsc t + + forwardGen (Black _ _ fg _) t2 = do + size <- getTrials + t1 <- gets snd + return $ U.replicate size $ forward fg t1 t2 + + maxStep _ = 100
+ src/Quant/Models/Dupire.hs view
@@ -0,0 +1,45 @@+{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ExistentialQuantification #-} + +module Quant.Models.Dupire ( + Dupire (..) +) where + +import Data.Random +import Control.Monad.State +import Quant.ContingentClaim +import Quant.MonteCarlo +import Quant.YieldCurve +import qualified Data.Vector.Unboxed as U + +-- | 'Dupire' represents a Dupire-style local vol model. +data Dupire = forall a b . (YieldCurve a, YieldCurve b) => Dupire { + dupireInitial :: Double -- ^ Initial asset level + , dupireFunc :: Double -> Double -> Double -- ^ Local vol function taking a time to maturity and a level + , mertonForwardGen :: a -- ^ 'YieldCurve' to generate forwards + , mertonDiscounter :: b } -- ^ 'YieldCurve' to generate discount rates + +instance Discretize Dupire where + initialize (Dupire s _ _ _) trials = put (Observables [U.replicate trials s], 0) + + evolve' d@(Dupire _ f _ _) t2 anti = do + (Observables (stateVec:_), t1) <- get + fwd <- forwardGen d t2 + let vols = U.map (f t1) stateVec + grwth = U.map (\(fwdVal, v) -> (fwdVal - v * v / 2) / (t2-t1)) $ U.zip fwd vols + postVal <- U.forM (U.zip3 grwth stateVec vols) $ \ ( g,x,v ) -> do + normResid <- lift stdNormal + if anti then + return $ x * exp (g - normResid*v) + else + return $ x * exp (g + normResid*v) + put (Observables [postVal], t2) + + discounter (Dupire _ _ _ dsc) t = do + size <- getTrials + return $ U.replicate size $ disc dsc t + + forwardGen (Dupire _ _ fg _) t2 = do + t1 <- gets snd + size <- getTrials + return $ U.replicate size $ forward fg t1 t2
+ src/Quant/Models/Heston.hs view
@@ -0,0 +1,59 @@+{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleInstances #-} + + +module Quant.Models.Heston ( + Heston (..) +) where + +import Quant.YieldCurve +import Data.Random +import Quant.Models +import Control.Monad.State +import Quant.MonteCarlo +import Quant.ContingentClaim +import qualified Data.Vector.Unboxed as U + +-- | 'Heston' represents a Heston model (i.e. stochastic volatility). +data Heston = forall a b . (YieldCurve a, YieldCurve b) => Heston { + hestonInit :: Double -- ^ Initial asset level. + , hestonV0 :: Double -- ^ Initial variance + , hestonVF :: Double -- ^ Mean-reversion variance + , hestonLambda :: Double -- ^ Vol-vol + , hestonCorrel :: Double -- ^ Correlation between processes + , hestonMeanRev :: Double -- ^ Mean reversion speed + , hestonForwardGen :: a -- ^ 'YieldCurve' to generate forwards + , hestonDisc :: b } -- ^ 'YieldCurve' to generate discounts + +instance Discretize Heston where + initialize (Heston s v0 _ _ _ _ _ _) trials = put (Observables [U.replicate trials s, + U.replicate trials v0 ], 0) + + evolve' h@(Heston _ _ vf l rho eta _ _) t2 anti = do + (Observables (sState:vState:_), t1) <- get + fwd <- forwardGen h t2 + let grwth = U.map (\(g, v) -> (g - v/2) * (t2-t1)) (U.zip fwd vState) + t = t2-t1 + states <- U.forM (U.zip3 grwth sState vState) $ \ ( g, x, v ) -> do + resid1 <- lift stdNormal + resid2' <- lift stdNormal + let + op = if anti then (-) else (+) + resid2 = rho * resid1 + sqrt (1-rho*rho) * resid2' + v' = (sqrt v `op` (eta/2.0*sqrt t* resid2))^(2 :: Int)-l*(v-vf)*t-eta*eta*t/4.0 + s' = x * exp (g `op` (resid1*sqrt (v*(t2-t1)))) + return (s', abs v') + let newS = U.map fst states + newV = U.map snd states + put (Observables [newS, newV], t2) + + discounter (Heston _ _ _ _ _ _ _ d) t = do + size <- getTrials + return $ U.replicate size $ disc d t + + forwardGen (Heston _ _ _ _ _ _ fg _) t2 = do + size <- getTrials + t1 <- gets snd + return $ U.replicate size $ forward fg t1 t2 + + maxStep _ = 1/250
+ src/Quant/Models/Merton.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ExistentialQuantification #-} + +module Quant.Models.Merton ( + Merton (..) +) where + +import Data.Random +import Data.Random.Distribution.Poisson +import Control.Monad.State +import Quant.MonteCarlo +import Quant.ContingentClaim +import Quant.YieldCurve +import qualified Data.Vector.Unboxed as U + +-- | 'Merton' represents a Merton model (Black-Scholes w/ jumps). +data Merton = forall a b . (YieldCurve a, YieldCurve b) => Merton { + mertonInitial :: Double -- ^ Initial asset level + , mertonVol :: Double -- ^ Asset volatility + , mertonIntensity :: Double -- ^ Intensity of Poisson process + , mertonJumpMean :: Double -- ^ Average size of jump + , mertonJumpVol :: Double -- ^ Volatility of jumps + , mertonForwardGen :: a -- ^ 'YieldCurve' to generate forwards + , mertonDiscounter :: b } -- ^ 'YieldCurve' to generate discount rates + +--instance CharFunc Merton where + -- charFunc (Merton s vol intensity mu sig fg) t k = charFunc (Black s vol fg) t k * addon + -- where + -- i = 0 :+ 1 + -- inner1 = exp (mu + sig*sig/2) :+ 0 + -- inner2 = exp $ i * k * (mu :+ 0) - k*k*(sig :+ 0) * (sig :+ 0)/2 + --addon = exp $ (intensity * t :+ 0) * (-i*k*(inner1 - 1) + inner2 - 1) + +instance Discretize Merton where + initialize (Merton s _ _ _ _ _ _) trials = put (Observables [U.replicate trials s], 0) + + evolve' m@(Merton _ vol intensity mu sig _ _) t2 anti = do + (Observables (stateVec:_), t1) <- get + fwd <- forwardGen m t2 + let correction = exp (mu + sig*sig /2.0) - 1 + grwth = U.map (\x -> (x - vol*vol/2 - intensity * correction) * (t2-t1)) fwd + postVal <- U.forM (U.zip grwth stateVec) $ \ ( g,x ) -> do + normResid1 <- lift stdNormal + normResid2 <- lift stdNormal + poissonResid <- lift $ integralPoisson (intensity * (t2-t1)) :: MonteCarlo (Observables, Double) Int + let poisson' = fromIntegral poissonResid + jumpterm = mu*poisson'+sig*sqrt poisson' * normResid2 + if anti then + return $ x * exp (g - normResid1*vol + jumpterm) + else + return $ x * exp (g + normResid1*vol + jumpterm) + put (Observables [postVal], t2) + + discounter (Merton _ _ _ _ _ _ dsc) t = do + size <- getTrials + return $ U.replicate size $ disc dsc t + + forwardGen (Merton _ _ _ _ _ fg _) t2 = do + size <- getTrials + t1 <- gets snd + return $ U.replicate size $ forward fg t1 t2
+ src/Quant/MonteCarlo.hs view
@@ -0,0 +1,160 @@+{-# LANGUAGE FlexibleContexts #-} + + +module Quant.MonteCarlo ( + -- * The MonteCarlo type. + MonteCarlo + , MonteCarloT + , runMC + + -- * The discretize typeclass. + , Discretize(..) + , OptionType(..) + + , getTrials + + ) +where + +import Quant.ContingentClaim +import Data.Random +import Control.Applicative +import Control.Monad.State +import Data.Functor.Identity +import Data.RVar +import System.Random.Mersenne.Pure64 +import qualified Data.Map as Map +import qualified Data.Vector.Unboxed as U + +-- | A monad transformer for Monte-Carlo calculations. +type MonteCarloT m s = StateT s (RVarT m) + +-- | Wraps the Identity monad in the 'MonteCarloT' transformer. +type MonteCarlo s a = MonteCarloT Identity s a + +-- | "Runs" a MonteCarlo calculation and provides the result of the computation. +runMC :: MonadRandom (StateT b Identity) => MonteCarlo s c -- ^ Monte Carlo computation. + -> b -- ^ Initial state. + -> s -- ^ Initial random-generator state. + -> c -- ^ Final result of computation. +runMC mc randState initState = flip evalState randState $ sampleRVarTWith lift (evalStateT mc initState) + + +{- | The 'Discretize' class defines those +models on which Monte Carlo simulations +can be performed. + +Minimal complete definition: 'initialize', 'discounter', 'forwardGen' and 'evolve''. +-} +class Discretize a where + + -- | Initializes a Monte Carlo simulation for a given number of runs. + initialize :: Discretize a => a -- ^ Model + -> Int -- ^ number of trials + -> MonteCarlo (Observables, Double) () + + -- | Evolves the internal states of the MC variables between two times. + evolve :: Discretize a => a -- ^ Model + -> Double -- ^ time to evolve to + -> Bool -- whether or not to use flipped variates + -> MonteCarlo (Observables, Double) () + evolve mdl t2 anti = do + (_, t1) <- get + let ms = maxStep mdl + if (t2-t1) < ms then + evolve' mdl t2 anti + else do + evolve' mdl (t1 + ms) anti + evolve mdl t2 anti + + -- | Stateful discounting function, takes a model and a time, and returns a vector of results. + discounter :: Discretize a => a -> Double -> MonteCarlo (Observables, Double) (U.Vector Double) + + -- | Stateful forward generator for a given model at a certain time. + forwardGen :: Discretize a => a -> Double -> MonteCarlo (Observables, Double) (U.Vector Double) + + -- | Internal function to evolve a model to a given time. + evolve' :: Discretize a => a -- ^ model + -> Double -- ^ time to evolve to + -> Bool -- ^ whether or not to use flipped variates + -> MonteCarlo (Observables, Double) () -- ^ computation result + + -- | Determines the maximum size time-step for discretization purposes. Defaults to 1/250. + maxStep :: Discretize a => a -> Double + maxStep _ = 1/250 + + -- | Perform a simulation of a compiled basket of contingent claims. + simulateState :: Discretize a => + a -- ^ model + -> ContingentClaimBasket -- ^ compilied basket of claims + -> Int -- ^ number of trials + -> Bool -- ^ antithetic? + -> MonteCarlo (Observables, Double) Double -- ^ computation result + simulateState modl (ContingentClaimBasket cs ts) trials anti = do + initialize modl trials + avg <$> process Map.empty (U.replicate trials 0) cs ts + where + process m cfs ccs@(c@(ContingentClaim' t _ _):cs') (obsT:ts') = + if t >= obsT then do + evolve modl obsT anti + obs <- gets fst + let m' = Map.insert obsT obs m + process m' cfs ccs ts' + else do + evolve modl t anti + let cfs' = processClaimWithMap c m + d <- discounter modl obsT + let cfs'' = cfs' |*| d + process m (cfs |+| cfs'') cs' (obsT:ts') + + process m cfs (c:cs') [] = do + d <- discounter modl (payoutTime c) + let cfs' = d |*| processClaimWithMap c m + process m (cfs |+| cfs') cs' [] + + process _ cfs _ _ = return cfs + + v1 |+| v2 = U.zipWith (+) v1 v2 + v1 |*| v2 = U.zipWith (*) v1 v2 + + avg v = U.sum v / fromIntegral (U.length v) + + -- | Runs a simulation for a 'ContingentClaim'. + runSimulation :: (Discretize a, + MonadRandom (StateT b Identity)) => + a -- ^ model + -> ContingentClaim -- ^ claims to value + -> b -- ^ initial random state + -> Int -- ^ trials + -> Bool -- ^ whether to use antithetic variables + -> Double -- ^ final value + runSimulation modl ccs seed trials anti = runMC run seed (Observables [], 0) + where + run = simulateState modl (ccBasket ccs) trials anti + + -- | Like 'runSimulation', but splits the trials in two and does antithetic variates. + runSimulationAnti :: (Discretize a, + MonadRandom (StateT b Identity)) => + a -> ContingentClaim -> b -> Int -> Double + runSimulationAnti modl ccs seed trials = (runSim True + runSim False) / 2 + where runSim = runSimulation modl ccs seed (trials `div` 2) + + -- | 'runSimulation' with a default random number generator. + quickSim :: Discretize a => a -> ContingentClaim -> Int -> Double + quickSim mdl opts trials = runSimulation mdl opts (pureMT 500) trials False + + -- | 'runSimulationAnti' with a default random number generator. + quickSimAnti :: Discretize a => a -> ContingentClaim -> Int -> Double + quickSimAnti mdl opts trials = runSimulationAnti mdl opts (pureMT 500) trials + +-- | Utility function to get the number of trials. +getTrials :: MonteCarlo (Observables, Double) Int +getTrials = U.length <$> gets (obsHead . fst) + + +processClaimWithMap :: ContingentClaim' -> Map.Map Double Observables -> U.Vector Double +processClaimWithMap (ContingentClaim' _ c obs) m = c vals + where + vals = map (\(t , g , f) -> U.map f . g $ m Map.! t) obs + +
+ src/Quant/Test.hs view
@@ -0,0 +1,63 @@+module Quant.Test ( + baseYC + , val + , black + , opt + , val' + , val'' + , val''' + , val'''' + , heston + , opt' + , opt'' + ) +where + +import Quant.MonteCarlo +import Quant.YieldCurve +import Quant.ContingentClaim +import Quant.Models.Black +import Quant.Models.Heston + +baseYC = FlatCurve 0.05 --create a flat yield curve with a 5% rate + +black = Black + 100 --initial stock price + 0.2 --volatility + baseYC --forward generator + baseYC --discount function + +opt = vanillaOption Put 100 1 --make a vanilla put, struck at 100, maturing at time 1 + +val = quickSim black opt 10000 --Run a Monte Carlo on opt in a a black model with 10000 trials + +opt' = multiplier 100 + $ vanillaOption Call 100 1 ++ short (vanillaOption Call 120 1) --Make a call spread with a 100 unit notional + +val' = quickSimAnti black opt' 10000 --Run a Monte Carlo on the call spread; use antithetic variates + --Returns + +black' = Black + 100 --initial stock price + 0.2 --volatility + (NetYC (FlatCurve 0.05) (FlatCurve 0.02)) --forward generator, now with a 2% dividend yield + baseYC --discount rate + +val'' = quickSimAnti black' opt' 10000 + +--Let's try it with a Heston model +heston = Heston + 100 + 0.04 --initial variance + 0.04 --final variance + 0.2 --volvol + (-0.7) --correlation between processes + 1.0 --mean reversion speed + baseYC --forward generator + baseYC --discount function + +val''' = quickSimAnti heston opt' 10000 --price the call spread in the Heston model + +opt'' = terminalOnly 1 $ \x -> x*x --create an option that pays off based on the square of its underlying + +val'''' = quickSimAnti heston opt'' 10000 --price it in the Heston model
+ src/Quant/VolSurf.hs view
@@ -0,0 +1,26 @@+module Quant.VolSurf ( + VolSurf (..) + , FlatSurf (..) +) where + + +{- | The 'VolSurf' class defines the +basic operations of a volatility surface. + +Minimal complete definition: 'vol'. +-} +class VolSurf a where + -- | Calculate the implied vol for a given strike/maturity. + vol :: VolSurf a => a -> Double -> Double -> Double + + -- | Calculate the variance at a given strike/maturity. + var :: VolSurf a => a -> Double -> Double -> Double + var vs s t = v*v*t + where v = vol vs s t + +-- |A flat curve is just a flat curve with one continuously +-- compounded rate at all points on the curve. +data FlatSurf = FlatSurf Double + +instance VolSurf FlatSurf where + vol (FlatSurf x) _ _ = x
+ src/Quant/YieldCurve.hs view
@@ -0,0 +1,36 @@+module Quant.YieldCurve ( + YieldCurve (..) + , FlatCurve (..) + , NetYC (..) +) where + + +{- | The 'YieldCurve' class defines the +basic operations of a yield curve. + +Minimal complete definition: 'disc'. +-} +class YieldCurve a where + -- | Calculate the discount factor for a given maturity. + disc :: YieldCurve a => a -> Double -> Double + + -- | Calculate the forward rate between a t1 and t2 + forward :: YieldCurve a => a -> Double -> Double -> Double + forward yc t1 t2 = (/(t2-t1)) $ log $ disc yc t1 / disc yc t2 + + -- | Calculate the spot rate for a given maturity. + spot :: YieldCurve a => a -> Double -> Double + spot yc t = forward yc 0 t + +-- |A flat curve is just a flat curve with one continuously +-- compounded rate at all points on the curve. +data FlatCurve = FlatCurve Double + +instance YieldCurve FlatCurve where + disc (FlatCurve r) t = exp (-r*t) + +-- | 'YieldCurve' that represents the difference between two 'YieldCurve's. +data NetYC a = NetYC a a + +instance YieldCurve a => YieldCurve (NetYC a) where + disc (NetYC yc1 yc2) t = disc yc1 t / disc yc2 t