diff --git a/hquantlib.cabal b/hquantlib.cabal --- a/hquantlib.cabal +++ b/hquantlib.cabal @@ -1,5 +1,5 @@ name: hquantlib -version: 0.0.3.0 +version: 0.0.3.1 license: LGPL license-file: LICENSE author: Pavel Ryzhov @@ -73,13 +73,25 @@ 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 + vector-algorithms >= 0.7.0.0 && < 0.8.0.0, + HerbiePlugin == 0.2.0.0 hs-source-dirs: src ghc-options: -Wall if flag(optimize) - ghc-options: -funbox-strict-fields -O2 -fspec-constr -fdicts-cheap + ghc-options: -funbox-strict-fields -O2 -fspec-constr -fdicts-cheap -fplugin=Herbie +executable mctest + default-language: Haskell2010 + main-is : Tests/McTest.hs + hs-source-dirs : src + ghc-options : -threaded -rtsopts + build-depends : base, + hquantlib, + parallel, + mersenne-random, + containers + Test-Suite main-test default-language: Haskell2010 type : exitcode-stdio-1.0 @@ -89,5 +101,5 @@ test-framework >= 0.8 && < 0.9, test-framework-hunit >= 0.3.0 && < 0.4.0, test-framework-quickcheck2 >= 0.3.0.0 && < 0.4.0, - QuickCheck >= 2.8.0 && < 2.9.0, + QuickCheck >= 2.8.0 && < 3.0.0, HUnit >= 1.2.5.2 && < 2.0.0.0 diff --git a/src/QuantLib/Math/Copulas.hs b/src/QuantLib/Math/Copulas.hs --- a/src/QuantLib/Math/Copulas.hs +++ b/src/QuantLib/Math/Copulas.hs @@ -54,6 +54,7 @@ {- Private implementations -} +{-# ANN aliMikhailHaqCopula "NoHerbie" #-} aliMikhailHaqCopula :: (Fractional a, Ord a) => a -> a -> a -> Maybe a aliMikhailHaqCopula theta x y | theta >= -1.0 && theta <= 1.0 @@ -61,6 +62,7 @@ | otherwise = Nothing +{-# ANN farlieGumbelMorgenstern "NoHerbie" #-} farlieGumbelMorgenstern :: (Fractional a, Ord a) => a -> a -> a -> Maybe a farlieGumbelMorgenstern theta x y | theta >= -1.0 && theta <= 1.0 @@ -72,6 +74,7 @@ implemented in Haskell by Nicholas Pezolano npezolano "at" gmail.com -} +{-# ANN claytonCopula "NoHerbie" #-} claytonCopula :: Double -> Double -> Double -> Maybe Double claytonCopula theta x y | theta ==0 @@ -92,6 +95,7 @@ | theta == 0.0 = Nothing | otherwise = Just (-1.0/theta * log (1 + (exp (-theta*x) - 1.0) * (exp (-theta*y) -1.0) / (exp (-theta) - 1.0) )) +{-# ANN galambosCopula "NoHerbie" #-} galambosCopula :: (Floating a, Ord a) => a -> a -> a -> Maybe a galambosCopula theta x y | theta <= 0.0 = Nothing @@ -102,6 +106,7 @@ | rho >= -1.0 && rho <= 1.0 = undefined | otherwise = Nothing +{-# ANN gumbelCopula "NoHerbie" #-} gumbelCopula :: (Floating a, Ord a) => a -> a -> a -> Maybe a gumbelCopula theta x y | theta >= 1.0 = Just (exp ( - ( (-log x) ** theta + (-log y) ** theta) ** (1.0/theta))) @@ -120,6 +125,7 @@ | a >= 0.0 && b >= 0.0 = Just (min (y * (x ** (1-a))) (x * (y ** (1-b)))) | otherwise = Nothing +{-# ANN plackettCopula "NoHerbie" #-} plackettCopula :: (Floating a, Ord a) => a -> a -> a -> Maybe a plackettCopula theta x y | theta >= 0.0 && theta /= 1.0 = Just $ (sumXyTheta1 - sqrt (sumXyTheta1 * sumXyTheta1 - 4.0 * x * y * theta * theta1))/(2*theta1) diff --git a/src/QuantLib/Math/InverseNormal.hs b/src/QuantLib/Math/InverseNormal.hs --- a/src/QuantLib/Math/InverseNormal.hs +++ b/src/QuantLib/Math/InverseNormal.hs @@ -64,7 +64,7 @@ -- | Computes the inverse cumulative standard normal distribution N(0, 1) inverseNormal :: Double -> Double -inverseNormal x +inverseNormal x | x < xlow = inverseInLowerRegion z | x <= xhigh = inverseInCentralRegion z | otherwise = inverseInHigherRegion z @@ -83,19 +83,21 @@ tolerance = 42*ulp {-# INLINE inverseRecovery #-} +{-# ANN inverseInLowerRegion "NoHerbie" #-} inverseInLowerRegion :: Double -> Double inverseInLowerRegion x = (((((c1*z+c2)*z+c3)*z+c4)*z+c5)*z+c6) / ((((d1*z+d2)*z+d3)*z+d4)*z+1.0) where z = sqrt (-2.0*log x) {-# INLINE inverseInLowerRegion #-} +{-# ANN inverseInCentralRegion "NoHerbie" #-} inverseInCentralRegion :: Double -> Double inverseInCentralRegion x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*z / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1.0) where r = z*z z = x - 0.5 {-# INLINE inverseInCentralRegion #-} +{-# ANN inverseInHigherRegion "NoHerbie" #-} inverseInHigherRegion :: Double -> Double inverseInHigherRegion x = -(((((c1*z+c2)*z+c3)*z+c4)*z+c5)*z+c6) / ((((d1*z+d2)*z+d3)*z+d4)*z+1.0) where z = sqrt (-2.0 * log (1.0 - x)) {-# INLINE inverseInHigherRegion #-} - 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 @@ -54,7 +54,7 @@ } -- | This pricer gets the last point of path -data LastPointPricer = LastPointPricer Dot +newtype LastPointPricer = LastPointPricer Dot instance PathPricer LastPointPricer where ppPrice _ path = LastPointPricer (last path) diff --git a/src/QuantLib/Models/Volatility.hs b/src/QuantLib/Models/Volatility.hs --- a/src/QuantLib/Models/Volatility.hs +++ b/src/QuantLib/Models/Volatility.hs @@ -1,16 +1,16 @@ -module QuantLib.Models.Volatility +module QuantLib.Models.Volatility ( Volatility , Estimation (..) , VolatilityEstimator (..) , VolatilityEstimatorAlgorithm (..) ) where -import QuantLib.Prices (IntervalPrice(..)) -import QuantLib.TimeSeries (IntervalPriceSeries) +import QuantLib.Prices (IntervalPrice (..)) +import QuantLib.TimeSeries (IntervalPriceSeries) -import qualified Data.Map as M -import Statistics.Sample (stdDev, fastVarianceUnbiased) +import qualified Data.Map as M import qualified Data.Vector.Unboxed as U +import Statistics.Sample (fastVarianceUnbiased, stdDev) -- | Volatility type type Volatility = Double @@ -52,6 +52,7 @@ simple :: IntervalPriceSeries -> Estimation simple = Estimation . stdDev . toLogArray +{-# ANN simpleDriftLess "NoHerbie" #-} simpleDriftLess :: IntervalPriceSeries -> Estimation simpleDriftLess = Estimation . sqrt . divByN . U.foldl' accum (T 0.0 0) . toLogArray where @@ -78,12 +79,13 @@ varRS = combine . M.foldl' point (T 0.0 0) where combine (T a n) = a / fromIntegral n - point (T a n) (IntervalPrice o h l c) = + point (T a n) (IntervalPrice o h l c) = T (a + logBase c h * logBase o h + logBase c l * logBase o l) (n + 1) toSimpleLogWith :: (IntervalPrice -> Double) -> IntervalPriceSeries -> U.Vector Double toSimpleLogWith f = U.fromList . map (f . snd) . M.toAscList +{-# ANN yangZhang "NoHerbie" #-} yangZhang :: IntervalPriceSeries -> Estimation yangZhang prices = Estimation $ sqrt (varO + k * varC + (1.0 - k) * varRS prices) where diff --git a/src/QuantLib/PricingEngines/BlackFormula.hs b/src/QuantLib/PricingEngines/BlackFormula.hs --- a/src/QuantLib/PricingEngines/BlackFormula.hs +++ b/src/QuantLib/PricingEngines/BlackFormula.hs @@ -2,10 +2,10 @@ ( blackFormulaImpliedStdDev ) where -import Data.Maybe -import QuantLib.Options -import Numeric.GSL.Root -import Numeric.GSL.Special.Erf +import Data.Maybe +import Numeric.GSL.Root +import Numeric.GSL.Special.Erf +import QuantLib.Options blackFormulaImpliedStdDev :: OptionType->Double->Double->Double->Double->Double->Maybe Double->Double->Int->Maybe Double blackFormulaImpliedStdDev opType strike forward blackPrice discount displacement guess accuracy maxIter @@ -21,6 +21,7 @@ blackFunction = blackImpliedStdDevHelper opType strike forward blackPrice displacement ([stdDev], _) = root DNewton accuracy maxIter blackFunction [realGuess] +{-# ANN blackImpliedStdDevHelper "NoHerbie" #-} blackImpliedStdDevHelper :: OptionType-> Double-> Double-> Double-> Double-> [Double]-> [Double] blackImpliedStdDevHelper opType strike forward blackPrice displacement [x] = [max 0.0 result - blackPrice] 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 @@ -8,14 +8,14 @@ , mkInverseNormal ) where -import System.Random.Mersenne -import QuantLib.Math.InverseNormal +import QuantLib.Math.InverseNormal +import System.Random.Mersenne -- | Box-Muller method data BoxMuller = BoxMuller { - bmFirst :: Bool, - bmSecondValue :: Double, - bmRng :: MTGen + bmFirst :: Bool, + bmSecondValue :: Double, + bmRng :: MTGen } mkNormalGen :: IO BoxMuller @@ -38,25 +38,31 @@ instance NormalGenerator BoxMuller where ngMkNew _ = mkNormalGen - ngGetNext (BoxMuller True _ rng) = do - (!r, !s1, !s2) <- getRs - let !ratio = sqrt (-2.0 * log r / 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) - - ngGetNext (BoxMuller False !s !r) = return (s, BoxMuller True s r) + ngGetNext = boxMullerGetNext +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) + +{-# 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 @@ -64,7 +70,7 @@ mkInverseNormal = do rng <- newMTGen Nothing return $! InverseNormal rng - + instance NormalGenerator InverseNormal where ngMkNew _ = mkInverseNormal ngGetNext gen@(InverseNormal rng) = do diff --git a/src/Tests/McTest.hs b/src/Tests/McTest.hs new file mode 100644 --- /dev/null +++ b/src/Tests/McTest.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +module Main where + +import Control.Monad +import Data.List +import qualified Data.Map as M +import QuantLib.Methods.MonteCarlo +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) + deriving (Show) + +toDouble :: Int -> Double +toDouble = fromIntegral + +addOnePath :: HistoSummary->MaxMinClosePricer->HistoSummary +addOnePath (HS m) (MMCP _ _ close) = HS newM + where (_, !newM) = M.insertLookupWithKey inserter roundedClose 1 m + !roundedClose = toDouble (round (close*10000))/10000 + inserter _ new_value old_value = old_value+new_value + +instance Summary HistoSummary MaxMinClosePricer where + sNorm _ _ = 0.0 -- we don't care about convergence now + sSummarize = foldl' addOnePath + +printMap :: HistoSummary->IO () +printMap (HS m) = forM_ list printPlain + where + printPlain (a, b) = putStrLn $ show a ++ "," ++ show b + list = M.toList m + +getHsSize :: HistoSummary -> Int +getHsSize (HS m) = M.size m + +main :: IO () +main = do + let summary = HS M.empty + let mmcp = MMCP 0.0 0.0 0.0 + let start = Dot 0.0 1.0 + let sp = GeometricBrownian 0.0 0.005 + let discrete= Euler 0.01 + rng <- mkInverseNormal + let pg = ProcessGenerator start 1000 sp rng discrete + let pmc = PathMonteCarlo summary mmcp pg + s <- monteCarlo pmc 50000 + -- printMap s + print (getHsSize s)