diff --git a/hquantlib.cabal b/hquantlib.cabal --- a/hquantlib.cabal +++ b/hquantlib.cabal @@ -1,5 +1,5 @@ name: hquantlib -version: 0.0.3.3 +version: 0.0.4.0 license: LGPL license-file: LICENSE author: Pavel Ryzhov @@ -19,7 +19,7 @@ source-repository this type: git location: https://github.com/paulrzcz/hquantlib.git - tag: 0.0.2.4 + tag: 0.0.3.4 flag optimize description : Enable optimizations for library and benchmarks @@ -48,6 +48,7 @@ QuantLib.Position QuantLib.Options QuantLib.Methods.MonteCarlo + QuantLib.Methods.Pricer other-modules: QuantLib.Currencies.America @@ -61,19 +62,21 @@ QuantLib.Time.Date QuantLib.Time.DayCounter QuantLib.Math.InverseNormal + QuantLib.Stochastic.PureMT build-depends: - base >3 && <5, - time >= 1.4.0.0 && < 1.7.0.0, - containers >= 0.5.0.0 && < 0.6.0.0, - hmatrix >= 0.17.0.0 && < 0.19.0.0, - hmatrix-gsl >= 0.17.0.0 && < 0.19.0.0, - hmatrix-special >= 0.4.0 && < 0.5.0, - parallel >= 3.2.0.0 && < 3.3.0.0, - mersenne-random >= 1.0.0.1 && < 2.0.0.0, - statistics >= 0.13.0.0 && < 0.14.0.0, - vector >= 0.11.0.0 && < 0.12.0.0, - vector-algorithms >= 0.7.0.0 && < 0.8.0.0 + base >3 && <5, + random >= 1.0 && < 2.0, + time >= 1.4.0.0 && < 1.7.0.0, + containers >= 0.5.0.0 && < 0.6.0.0, + hmatrix >= 0.17.0.0 && < 0.19.0.0, + hmatrix-gsl >= 0.17.0.0 && < 0.19.0.0, + hmatrix-special >= 0.4.0 && < 0.5.0, + parallel >= 3.2.0.0 && < 3.3.0.0, + mersenne-random-pure64 >= 0.2.0.0 && < 0.3.0.0, + statistics >= 0.13.0.0 && < 0.15.0.0, + vector >= 0.11.0.0 && < 0.13.0.0, + vector-algorithms >= 0.7.0.0 && < 0.8.0.0 hs-source-dirs: src ghc-options: -Wall @@ -87,6 +90,7 @@ ghc-options : -threaded -rtsopts other-modules : QuantLib.Math.InverseNormal QuantLib.Methods.MonteCarlo + QuantLib.Methods.Pricer QuantLib.Stochastic QuantLib.Stochastic.Discretize QuantLib.Stochastic.Process @@ -94,8 +98,9 @@ build-depends : base, hquantlib, parallel, - mersenne-random, - containers + mersenne-random-pure64, + containers, + time Test-Suite main-test default-language: Haskell2010 diff --git a/src/QuantLib/Methods/MonteCarlo.hs b/src/QuantLib/Methods/MonteCarlo.hs --- a/src/QuantLib/Methods/MonteCarlo.hs +++ b/src/QuantLib/Methods/MonteCarlo.hs @@ -13,37 +13,33 @@ -- | Summary type class aggregates all priced values of paths class PathPricer p => Summary m p | m->p where -- | Updates summary with given priced pathes - sSummarize :: m->[p]->m + sSummarize :: m -> [p] -> m -- | Defines a metric, i.e. calculate distance between 2 summaries - sNorm :: m->m->Double + sNorm :: m -> m -> Double -- | Path generator is a stochastic path generator class PathGenerator m where pgMkNew :: m->IO m - pgGenerate :: m->IO Path + pgGenerate :: Integer -> m->Path -- | Path pricer provides a price for given path class PathPricer m where - ppPrice :: m->Path->m + ppPrice :: m -> Path -> m -- | Monte Carlo engine function -monteCarlo :: (Summary s p, PathGenerator g) => PathMonteCarlo s p g->Int->IO s -monteCarlo (PathMonteCarlo s p g) size = do - priced <- mapM (const pricing) [1..size] - return $ sSummarize s priced - where pricing = do - !path <- pgGenerate g - return $! ppPrice p path +monteCarlo :: (Summary s p, PathGenerator g) => PathMonteCarlo s p g->Int->s +monteCarlo (PathMonteCarlo s p g) size = sSummarize s priced + where + !priced = map pricing [1..size] + pricing seed = ppPrice p (pgGenerate (fromIntegral seed) g) -- | Monte Carlo engine function. Parallelized version -monteCarloParallel :: (Summary s p, PathGenerator g) => PathMonteCarlo s p g->Int->IO s -monteCarloParallel (PathMonteCarlo s p g) size = do - priced <- mapM (const pricing) [1..size] `using` rpar - return $ sSummarize s priced - where pricing = do - !path <- pgGenerate g - return $! ppPrice p path +monteCarloParallel :: (Summary s p, PathGenerator g) => PathMonteCarlo s p g->Int->s +monteCarloParallel (PathMonteCarlo s p g) size = sSummarize s priced + where + !priced = map pricing [1..size] `using` rpar + pricing seed = ppPrice p (pgGenerate (fromIntegral seed) g) -- | Path-dependant Monte Carlo engine data PathMonteCarlo s p g = @@ -53,12 +49,6 @@ pmcGenerator :: g } --- | This pricer gets the last point of path -newtype LastPointPricer = LastPointPricer Dot - -instance PathPricer LastPointPricer where - ppPrice _ path = LastPointPricer (last path) - -- | Stochastic process generator data ProcessGenerator sp b d = ProcessGenerator { @@ -73,4 +63,5 @@ pgMkNew (ProcessGenerator start len process rnd d) = do newRnd <- ngMkNew rnd return $! ProcessGenerator start len process newRnd d - pgGenerate (ProcessGenerator start len sp b d) = generatePath b d sp len start + pgGenerate seed (ProcessGenerator start len sp b d) = generatePath newB d sp len start + where (_, newB) = ngSplitWithSeed seed b diff --git a/src/QuantLib/Methods/Pricer.hs b/src/QuantLib/Methods/Pricer.hs new file mode 100644 --- /dev/null +++ b/src/QuantLib/Methods/Pricer.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE BangPatterns #-} +module QuantLib.Methods.Pricer + ( MaxMinClosePricer (..) + , LastPointPricer (..) + , LogLastPointPricer (..) + ) where + +import QuantLib.Methods.MonteCarlo (PathPricer (..)) +import QuantLib.Stochastic.Process (Dot (..)) + +data MaxMinClosePricer = MMCP { + mmcpHigh :: Double, + mmcpLow :: Double, + mmcpClose :: Double + } deriving (Show) + +instance PathPricer MaxMinClosePricer where + ppPrice _ path = MMCP high low close + where !close = last xs + !high = maximum xs + !low = minimum xs + xs = map getX path + +-- | This pricer gets the last point of path +newtype LastPointPricer = LastPointPricer Double + +instance PathPricer LastPointPricer where + ppPrice _ = LastPointPricer <$> getX . last + +-- | This pricer estimates the log of difference between start and end of process +newtype LogLastPointPricer = LogLastPointPricer Double + +instance PathPricer LogLastPointPricer where + ppPrice _ path = LogLastPointPricer (log (lastX / firstX)) + where + lastX = getX $ last path + firstX = getX $ head path diff --git a/src/QuantLib/Stochastic.hs b/src/QuantLib/Stochastic.hs --- a/src/QuantLib/Stochastic.hs +++ b/src/QuantLib/Stochastic.hs @@ -3,3 +3,4 @@ import QuantLib.Stochastic.Process as Q import QuantLib.Stochastic.Discretize as Q import QuantLib.Stochastic.Random as Q +import QuantLib.Stochastic.PureMT as Q diff --git a/src/QuantLib/Stochastic/Process.hs b/src/QuantLib/Stochastic/Process.hs --- a/src/QuantLib/Stochastic/Process.hs +++ b/src/QuantLib/Stochastic/Process.hs @@ -1,9 +1,9 @@ -{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} module QuantLib.Stochastic.Process ( module QuantLib.Stochastic.Process ) where -import Control.Monad (foldM) import Data.List (foldl') import QuantLib.Stochastic.Random (NormalGenerator (..)) @@ -15,9 +15,9 @@ -- | 1D Stochastic process class StochasticProcess a where - drift :: a->Dot->Double - diff :: a->Dot->Double - evolve :: Discretize b=> b->a->Dot->Double->Dot + drift :: a -> Dot -> Double + diff :: a -> Dot -> Double + evolve :: Discretize b=> b -> a -> Dot -> Double -> Dot evolve discr p dot dw = Dot newT newX where !newT = getT dot + dDt p discr dot !newX = getX dot + dDrift p discr dot + dDiff p discr dot * dw @@ -30,16 +30,15 @@ type Path = [Dot] -- | Generates sample path for given stochastic process under discretization and normal generator for given amount of steps, starting from x0 -generatePath :: (StochasticProcess a, NormalGenerator b, Discretize c) => b->c->a->Int->Dot->IO Path -generatePath rnd discr sp steps x0 = do - (!list, _) <- foldM generator ([], rnd) [1..steps] - let !path = foldl' evolver [x0] list - return $! reverse path - where evolver p dw = evolve discr sp (head p) dw : p - generator (list, r) _ = do - (!p, newRnd) <- ngGetNext r - return (p:list, newRnd) - +generatePath :: (StochasticProcess a, NormalGenerator b, Discretize c) => b->c->a->Int->Dot->Path +generatePath rnd discr sp steps x0 = reverse path + where + (!list, _) = foldl' generator ([], rnd) [1..steps] + !path = foldl' evolver [x0] list + evolver p dw = evolve discr sp (head p) dw : p + generator (l, r) _ = (p:l, newRnd) + where + (!p, newRnd) = ngGetNext r -- | Geometric Brownian motion data GeometricBrownian = GeometricBrownian { diff --git a/src/QuantLib/Stochastic/PureMT.hs b/src/QuantLib/Stochastic/PureMT.hs new file mode 100644 --- /dev/null +++ b/src/QuantLib/Stochastic/PureMT.hs @@ -0,0 +1,36 @@ +module QuantLib.Stochastic.PureMT + ( + PureMT + , newPureMT + , randomDouble + , splitMT + , splitMTwithSeed + ) where + +import Data.Time.Calendar +import Data.Time.Clock +import System.CPUTime +import qualified System.Random.Mersenne.Pure64 as P + +data PureMT = PureMT P.PureMT Integer + +newPureMT :: IO PureMT +newPureMT = do + ct <- getCPUTime + t <- getCurrentTime + let seed = toModifiedJulianDay (utctDay t) + diffTimeToPicoseconds (utctDayTime t) + ct + return $ PureMT (P.pureMT $ fromIntegral seed) seed + +randomDouble:: PureMT -> (Double, PureMT) +randomDouble (PureMT mt seed) = (r, PureMT newMt seed) + where + (r, newMt) = P.randomDouble mt + +splitMT :: PureMT -> (PureMT, PureMT) +splitMT = splitMTwithSeed 1 + +splitMTwithSeed :: Integer -> PureMT -> (PureMT, PureMT) +splitMTwithSeed addedSeed mt@(PureMT _ seed) = (mt, PureMT newMt newSeed) + where + newSeed = seed + addedSeed + newMt = P.pureMT $ fromIntegral newSeed diff --git a/src/QuantLib/Stochastic/Random.hs b/src/QuantLib/Stochastic/Random.hs --- a/src/QuantLib/Stochastic/Random.hs +++ b/src/QuantLib/Stochastic/Random.hs @@ -1,7 +1,7 @@ {-# LANGUAGE BangPatterns #-} module QuantLib.Stochastic.Random ( BoxMuller - , createNormalGen +-- , createNormalGen , mkNormalGen , NormalGenerator (..) , InverseNormal @@ -9,22 +9,32 @@ ) where import QuantLib.Math.InverseNormal -import System.Random.Mersenne +import QuantLib.Stochastic.PureMT +class RandomGenerator a where + create :: IO a + next :: a -> (Double, a) + splitWithSeed :: Integer -> a -> (a, a) + +instance RandomGenerator PureMT where + create = newPureMT + next = randomDouble + splitWithSeed = splitMTwithSeed + -- | Box-Muller method -data BoxMuller = BoxMuller { +data BoxMuller a = BoxMuller { bmFirst :: Bool, bmSecondValue :: Double, - bmRng :: MTGen + bmRng :: a } -mkNormalGen :: IO BoxMuller +mkNormalGen :: RandomGenerator a => IO (BoxMuller a) mkNormalGen = do - rng <- newMTGen Nothing + rng <- create return $! createNormalGen rng -- | Creates normally distributed generator -createNormalGen :: MTGen->BoxMuller +createNormalGen :: RandomGenerator a => a -> BoxMuller a createNormalGen r = BoxMuller { bmFirst = True, bmSecondValue = 0.0, @@ -33,46 +43,55 @@ -- | Normally distributed generator class NormalGenerator a where - ngGetNext :: a -> IO (Double, a) + ngGetNext :: a -> (Double, a) ngMkNew :: a -> IO a + ngSplit :: a -> (a, a) + ngSplit = ngSplitWithSeed 1 + ngSplitWithSeed :: Integer -> a -> (a, a) -instance NormalGenerator BoxMuller where +instance RandomGenerator a => NormalGenerator (BoxMuller a) where ngMkNew _ = mkNormalGen ngGetNext = boxMullerGetNext + ngSplitWithSeed seed x = (x { bmRng = rng1 }, x { bmRng = rng2 }) + where + (rng1, rng2) = splitWithSeed seed (bmRng x) -boxMullerGetNext :: BoxMuller -> IO (Double, BoxMuller) -boxMullerGetNext (BoxMuller True _ rng) = do - (!r, !s1, !s2) <- getRs - let !ratio = boxMullerRatio r - let !bm = BoxMuller { - bmFirst = False, - bmSecondValue = s2*ratio, - bmRng = rng - } - return (s1*ratio, bm) - where getRs = do - x1 <- random rng :: IO Double - x2 <- random rng :: IO Double - let !s1 = 2.0*x1-1.0 - let !s2 = 2.0*x2-1.0 - let !r = s1*s1 + s2*s2 - if r>=1.0 || r<=0.0 then getRs else return (r, s1, s2) -boxMullerGetNext (BoxMuller False !s !r) = return (s, BoxMuller True s r) +boxMullerGetNext :: RandomGenerator a => BoxMuller a -> (Double, BoxMuller a) +boxMullerGetNext (BoxMuller True _ rng) = (s1*ratio, BoxMuller { + bmFirst = False, + bmSecondValue = s2*ratio, + bmRng = g2 + }) + where + (x1, g1) = next rng + (x2, g2) = next g1 + (r, s1, s2) = getRs + ratio = boxMullerRatio r + getRs = + let + as1 = 2.0*x1-1.0 + as2 = 2.0*x2-1.0 + ar = s1*s1 + s2*s2 + in + if r>=1.0 || r<=0.0 then getRs else (ar, as1, as2) +boxMullerGetNext (BoxMuller False !s !r) = (s, BoxMuller True s r) {-# ANN boxMullerRatio "NoHerbie" #-} boxMullerRatio :: Double -> Double boxMullerRatio r = sqrt (-2.0 * log r / r) -- | Normal number generation using inverse cummulative normal distribution -data InverseNormal = InverseNormal MTGen +newtype InverseNormal a = InverseNormal a -mkInverseNormal :: IO InverseNormal +mkInverseNormal :: RandomGenerator a => IO (InverseNormal a) mkInverseNormal = do - rng <- newMTGen Nothing + rng <- create return $! InverseNormal rng -instance NormalGenerator InverseNormal where +instance RandomGenerator a => NormalGenerator (InverseNormal a) where ngMkNew _ = mkInverseNormal - ngGetNext gen@(InverseNormal rng) = do - x <- random rng :: IO Double - return (inverseNormal x, gen) + ngGetNext (InverseNormal rng) = (inverseNormal x, InverseNormal newRng) + where (x, newRng) = next rng + ngSplitWithSeed seed (InverseNormal x) = (InverseNormal x1, InverseNormal x2) + where + (x1, x2) = splitWithSeed seed x diff --git a/src/Tests/McTest.hs b/src/Tests/McTest.hs --- a/src/Tests/McTest.hs +++ b/src/Tests/McTest.hs @@ -7,22 +7,10 @@ import Data.List import qualified Data.Map as M import QuantLib.Methods.MonteCarlo +import QuantLib.Methods.Pricer (MaxMinClosePricer (..)) import QuantLib.Stochastic -data MaxMinClosePricer = MMCP { - mmcpHigh :: Double, - mmcpLow :: Double, - mmcpClose :: Double - } deriving (Show) - -instance PathPricer MaxMinClosePricer where - ppPrice _ path = MMCP high low close - where !close = last xs - !high = maximum xs - !low = minimum xs - xs = map getX path - -data HistoSummary = HS (M.Map Double Int) +newtype HistoSummary = HS (M.Map Double Int) deriving (Show) toDouble :: Int -> Double @@ -54,9 +42,9 @@ let start = Dot 0.0 1.0 let sp = GeometricBrownian 0.0 0.005 let discrete= Euler 0.01 - rng <- mkInverseNormal + rng <- mkInverseNormal :: IO (InverseNormal PureMT) let pg = ProcessGenerator start 1000 sp rng discrete let pmc = PathMonteCarlo summary mmcp pg - s <- monteCarlo pmc 50000 - -- printMap s + let s = monteCarlo pmc 100000 + printMap s print (getHsSize s)