hquantlib 0.0.3.3 → 0.0.4.0

raw patch · 8 files changed

+182/−106 lines, 8 filesdep +mersenne-random-pure64dep +randomdep −mersenne-randomdep ~hmatrixdep ~hmatrix-gsldep ~hmatrix-specialPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: mersenne-random-pure64, random

Dependencies removed: mersenne-random

Dependency ranges changed: hmatrix, hmatrix-gsl, hmatrix-special, time

API changes (from Hackage documentation)

Files

• hquantlib.cabal +20/−15
• src/QuantLib/Methods/MonteCarlo.hs +16/−25
• src/QuantLib/Methods/Pricer.hs +37/−0
• src/QuantLib/Stochastic.hs +1/−0
• src/QuantLib/Stochastic/Process.hs +14/−15
• src/QuantLib/Stochastic/PureMT.hs +36/−0
• src/QuantLib/Stochastic/Random.hs +53/−34
• src/Tests/McTest.hs +5/−17

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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 {

@@ -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

@@ -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

@@ -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)