quantfin 0.1.0.2 → 0.2.0.0

raw patch · 17 files changed

+499/−173 lines, 17 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Test.hs +18/−6
• quantfin.cabal +3/−2
• src/Quant/ContingentClaim.hs +85/−49
• src/Quant/Math/Interpolation.hs +2/−0
• src/Quant/Math/Utilities.hs +11/−1
• src/Quant/Models/Black.hs +12/−13
• src/Quant/Models/Dupire.hs +11/−9
• src/Quant/Models/Heston.hs +7/−6
• src/Quant/Models/Merton.hs +10/−9
• src/Quant/Models/Processes.hs +5/−5
• src/Quant/MonteCarlo.hs +63/−61
• src/Quant/RNG/MWC64X.hs +88/−0
• src/Quant/Time.hs +1/−1
• src/Quant/Types.hs +101/−8
• src/Quant/VectorOps.hs +77/−0
• src/Quant/VolSurf.hs +1/−1
• src/Quant/YieldCurve.hs +4/−2

@@ -19,10 +19,10 @@ 			baseYC  --discount function
 
 --make a vanilla put, struck at 100, maturing at time 1
-vanopt :: ContingentClaim
+vanopt :: ContingentClaim1
 vanopt = vanillaOption Call 100 (Time 1) --built in function
 
-vanopt' :: ContingentClaim
+vanopt' :: ContingentClaim1
 vanopt' = specify $ do
 	x <- monitor (Time 1)
 	return $ CashFlow (Time 1) (max (x - 100) 0) --roll your own
@@ -31,8 +31,12 @@ vanoptPrice :: Double
 vanoptPrice = quickSim black vanopt 100000 
 
+vanoptPrice' :: Double
+vanoptPrice' = quickSim black vanopt' 100000 
+
+
 --Make a call spread with a 100 unit notional, using some handy combinators.
-cs :: ContingentClaim
+cs :: ContingentClaim1
 cs =  multiplier 100 
    $  vanillaOption Call 100 (Time 1) 
    <> short (vanillaOption Call 120 (Time 1)) 
@@ -68,29 +72,37 @@ csHeston = quickSimAnti heston cs 100000
 
 --create an option that pays off based on the square of its underlying
-squareOpt :: ContingentClaim
+squareOpt :: ContingentClaim1
 squareOpt = terminalOnly (Time 1) $ \x -> x*x  --using the built in function
+
+squareOpt' :: ContingentClaim1
 squareOpt' = specify $ do --roll your own
 	x <- monitor (Time 1)
 	return $ CashFlow (Time 1) $ x*x
 squareOptPrice :: Double
 squareOptPrice = quickSimAnti black squareOpt 100000
 
+squareOptPrice' :: Double
+squareOptPrice' = quickSimAnti black squareOpt' 100000
+
+
 --create an option with a bizarre payoff
-bizarre :: ContingentClaim
+bizarre :: ContingentClaim1
 bizarre = specify $ do
   x <- monitor (Time 1)   --check the price of asset 0 @ time 1
   y <- monitor (Time 2)   --check the price of asset 0 @ time 2
   z <- monitor (Time 3)   --check the price of asset 0 @ time 3
-  return $ CashFlow (Time 4) $ x ^ 3 / y ^ 2 - 3 * z --payoff @ time 4
+  return $ CashFlow (Time 4) $ sin x * cos y / (z ** sin x) --payoff @ time 4
 bizarrePrice :: Double
 bizarrePrice = quickSimAnti black bizarre 100000
 
 main :: IO ()
 main = do
 	print vanoptPrice
+	print vanoptPrice'
 	print csPrice
 	print callSpreadAnti
 	print csHeston
 	print squareOptPrice
+	print squareOptPrice'
 	print bizarrePrice

@@ -10,7 +10,7 @@ -- PVP summary:      +-+------- breaking API changes
 --                   | | +----- non-breaking API additions
 --                   | | | +--- code changes with no API change
-version:             0.1.0.2
+version:             0.2.0.0
 
 -- A short (one-line) description of the package.
 synopsis:            Quant finance library in pure Haskell.
@@ -62,11 +62,12 @@                        Quant.Models.Dupire
                        Quant.Models.Heston
                        Quant.Models.Processes
-                       --Quant.RNG.MWC64X
+                       Quant.RNG.MWC64X
                        Quant.MonteCarlo
                        Quant.ContingentClaim
                        Quant.Types
                        Quant.Time
+                       Quant.VectorOps
   
   -- Modules included in this library but not exported.
   -- other-modules:       

@@ -1,9 +1,13 @@+{-# LANGUAGE RankNTypes #-}
+
 module Quant.ContingentClaim (
     -- * Types for modeling contingent claims.
     ContingentClaim (..)
+  , ContingentClaim1
+  , ContingentClaim2
+  , ContingentClaim3
+  , ContingentClaim4
   , CCProcessor (..)
-  , Observables (..)
-  , MCObservables
   , OptionType (..)
   , CashFlow (..)
   , CCBuilder
@@ -11,7 +15,11 @@   -- * Options and option combinators
   , specify
   , monitor
-  , monitorByNum
+  , monitor1
+  , monitor2
+  , monitor3
+  , monitor4
+  , monitor5
   , vanillaOption
   , binaryOption
   , straddle
@@ -35,44 +43,78 @@ import Quant.Time
 import qualified Data.Map as M
 
-type MCMap = M.Map Time MCObservables
-type PayoffFunc a = MCMap -> a
+-- | Key type for building contingent claims.
+--Monoid instance allows for trivial combinations of
+--contingent claims.
+newtype ContingentClaim a = ContingentClaim { unCC :: [CCProcessor a] }
 
+-- | Contingent claims with one observable.
+type ContingentClaim1 = forall a . Obs1 a => ContingentClaim a
+-- | Contingent claims with two observables.
+type ContingentClaim2 = forall a . Obs2 a => ContingentClaim a
+-- | Contingent claims with three observables.
+type ContingentClaim3 = forall a . Obs3 a => ContingentClaim a
+-- | Contingent claims with four observables.
+type ContingentClaim4 = forall a . Obs4 a => ContingentClaim a
 
-data CCProcessor = CCProcessor  {
-                      monitorTime      :: Time
-                    , payoutFunc       :: Maybe [PayoffFunc CashFlow]
-}
+instance Monoid (ContingentClaim a) where
+  mempty  = ContingentClaim []
+  mappend = combine
 
+-- |Basic element of a `ContingentClaim`.  Each element contains
+--a Time.  Each Time, the observables are stored in the map.
+--Also, optionally a payout function may be applied at any time step.
+data CCProcessor a = CCProcessor  {
+    monitorTime      :: Time
+  , payoutFunc       :: [M.Map Time a -> CashFlow]
+}
 
 type CCBuilder w r a = WriterT w (Reader r) a
 
-monitor :: Time -> CCBuilder ContingentClaim MCMap Double
-monitor = monitorByNum 0
+-- | 'monitor' gets the value of the first observable at a given time.
+monitor :: Obs1 a => Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitor = monitor1
 
-monitorByNum :: Int -> Time -> CCBuilder ContingentClaim MCMap Double
-monitorByNum idx t = do
-  tell $ ContingentClaim [CCProcessor t Nothing]
-  m <- lift ask
-  return $ obsGet (m M.! t) !! idx  --I know, I know.
+-- | 'monitor1' gets the value of the first observable at a given time.
+monitor1 :: Obs1 a => Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitor1 = monitorGeneric get1
 
-specify :: CCBuilder ContingentClaim MCMap CashFlow -> ContingentClaim
-specify x = w `mappend` ContingentClaim [CCProcessor (last0 w') (Just [f])]
+-- | 'monitor2' gets the value of the second observable at a given time.
+monitor2 :: Obs2 a => Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitor2 = monitorGeneric get2
+
+-- | 'monitor3' gets the value of the third observable at a given time.
+monitor3 :: Obs3 a => Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitor3 = monitorGeneric get3
+
+-- | 'monitor4' gets the value of the fourth observable at a given time.
+monitor4 :: Obs4 a => Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitor4 = monitorGeneric get4
+
+-- | 'monitor5' gets the value of the fifth observable at a given time.
+monitor5 :: Obs5 a => Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitor5 = monitorGeneric get5
+
+-- | 'monitorGeneric' applies any getter at a given time.
+monitorGeneric :: (a -> Double) -> Time -> CCBuilder (ContingentClaim a) (M.Map Time a) Double
+monitorGeneric f t = do
+  tell $ ContingentClaim [CCProcessor t []]
+  m <- ask
+  return $ f (m M.! t)
+
+-- | Pulls a ContingentClaim out of the CCBuilder monad.
+specify :: CCBuilder (ContingentClaim a) (M.Map Time a) CashFlow -> ContingentClaim a
+specify x = w `mappend` ContingentClaim [CCProcessor lastW [f]]
   where
     w  = runReader (execWriterT x) M.empty
     f  = runReader . liftM fst $ runWriterT x
-    w' = map monitorTime $ unCC w
+    lastW = last0 . map monitorTime $ unCC w
     -- Equivalent to Prelude's last, but with a default of zero
     last0 [] = Time 0
     last0 [y] = y
     last0 (_:ys) = last0 ys
 
-newtype ContingentClaim = ContingentClaim { unCC :: [CCProcessor] }
 
-instance Monoid ContingentClaim where
-  mempty  = ContingentClaim []
-  mappend = combine
-
 -- | Function to generate a vanilla put/call style payout.
 vanillaPayout :: OptionType  -- ^ Put or Call
               -> Double      -- ^ Strike
@@ -94,23 +136,23 @@ 
 -- | Takes a maturity time and a function and generates a ContingentClaim 
 --dependent only on the terminal value of the observable.
-terminalOnly :: Time -> (Double -> Double) -> ContingentClaim
+terminalOnly :: Obs1 a => Time -> (Double -> Double) -> ContingentClaim a
 terminalOnly t g = specify $ do
   x <- monitor t
   return $ CashFlow t $ g x
 
 -- | Takes an OptionType, a strike, and a time to maturity and generates a vanilla option.
-vanillaOption :: OptionType -> Double -> Time -> ContingentClaim
+vanillaOption :: Obs1 a => OptionType -> Double -> Time -> ContingentClaim a
 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 -> Time -> ContingentClaim
+binaryOption :: Obs1 a => OptionType -> Double -> Double -> Time -> ContingentClaim a
 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 -> [Time] -> Time -> ContingentClaim
+arithmeticAsianOption :: Obs1 a => OptionType -> Double -> [Time] -> Time -> ContingentClaim a
 arithmeticAsianOption pc strike obsTimes t = specify $ do
   x <- mapM monitor obsTimes
   let avg = sum x / fromIntegral (length obsTimes)
@@ -118,70 +160,64 @@ 
 -- | Takes an OptionType, a strike, observation times, time to
 --maturity and generates an arithmetic Asian option.
-geometricAsianOption :: OptionType -> Double -> [Time] -> Time -> ContingentClaim
+geometricAsianOption :: Obs1 a => OptionType -> Double -> [Time] -> Time -> ContingentClaim a
 geometricAsianOption pc strike obsTimes t = specify $ do
   x <- mapM monitor obsTimes
   let avg = product x ** (1 / fromIntegral (length obsTimes))
   return $ CashFlow t $ vanillaPayout pc strike avg
 
 -- | Scales up a contingent claim by a multiplier.
-multiplier :: Double -> ContingentClaim -> ContingentClaim
+multiplier :: Double -> ContingentClaim a -> ContingentClaim a
 multiplier notional cs = ContingentClaim $ map f (unCC cs)
-    where f (CCProcessor t g) = CCProcessor t $ fmap (fmap (scale.)) g
+    where f (CCProcessor t g) = CCProcessor t $ (fmap (scale.)) g
           scale (CashFlow dt amt) = CashFlow dt (amt*notional)
 
 -- | Flips the signs in a contingent claim to make it a short position.
-short :: ContingentClaim -> ContingentClaim
+short :: ContingentClaim a -> ContingentClaim a
 short = multiplier (-1)
 
 -- | Takes an amount and a time and generates a fixed cash flow.
-zcb :: Time -> Double -> ContingentClaim
+zcb :: Time -> Double -> ContingentClaim a
 zcb t amt = specify $ return $ CashFlow t amt
 
 -- | Takes a face value, an interest rate, a payment frequency and makes a fixed bond
-fixedBond :: Double -> Double -> Double -> Int -> ContingentClaim
+fixedBond :: Double -> Double -> Double -> Int -> ContingentClaim a
 fixedBond faceVal intRate freq pmts = zcb (Time $ fromIntegral pmts * freq) faceVal 
                                    <> mconcat (map f [1..pmts])
   where
     f x = zcb (Time $ fromIntegral x * freq) (faceVal * intRate * freq) 
 
 -- | Takes a time to maturity and generates a forward contract.
-forwardContract :: Time -> ContingentClaim
+forwardContract :: Obs1 a => Time -> ContingentClaim a
 forwardContract t = specify $ do
   x <- monitor t
   return $ CashFlow t x
 
 -- | A call spread is a long position in a low-strike call
 --and a short position in a high strike call.
-callSpread :: Double -> Double -> Time -> ContingentClaim
+callSpread :: Obs1 a => Double -> Double -> Time -> ContingentClaim a
 callSpread lowStrike highStrike t = mappend (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 -> Time -> ContingentClaim
+putSpread :: Obs1 a => Double -> Double -> Time -> ContingentClaim a
 putSpread lowStrike highStrike t = mappend (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 -> Time -> ContingentClaim
+straddle :: Obs1 a => Double -> Time -> ContingentClaim a
 straddle strike t = vanillaOption Put strike t <> vanillaOption Call strike t
 
 -- | Combines two contingent claims into one. 
-combine :: ContingentClaim -> ContingentClaim -> ContingentClaim
+combine :: ContingentClaim a -> ContingentClaim a -> ContingentClaim a
 combine (ContingentClaim x) (ContingentClaim y) = ContingentClaim $ combine' x y
   where
     combine' (cc1:ccs1) (cc2:ccs2)
       | monitorTime cc1 == monitorTime cc2 = let
-          (CCProcessor t mf)  = cc1
-          (CCProcessor _ mf') = cc2 in
-            case mf of
-              Nothing -> cc2 : combine' ccs1 ccs2
-              Just a  -> case mf' of
-                Nothing -> cc1 : combine' ccs1 ccs2
-                Just b  -> CCProcessor t (Just (a ++ b)) : combine' ccs1 ccs2
+            CCProcessor t mf  = cc1
+            CCProcessor _ mf' = cc2 in
+          CCProcessor t (mf++mf') : combine' ccs1 ccs2
       | monitorTime cc1 > monitorTime cc2 = cc2 : combine' (cc1:ccs1) ccs2
       | otherwise = cc1 : combine' ccs1 (cc2:ccs2)
-    combine' [] [] = []
-    combine' cs [] = cs
-    combine' [] cs = cs+    combine' cs1 cs2 = cs1 ++ cs2

@@ -21,6 +21,7 @@           wt1 = (x2-x) / (x2-x1)
           wt2 = (x-x1) / (x2-x1)
 linearInterpolator _ [y] _ = y
+linearInterpolator _ _ _ = error "Empty interpolation list."
 {-# INLINE linearInterpolator #-}
 
 logLinearInterpolator :: Interpolator1d
@@ -56,4 +57,5 @@             term = x-x1
             h' = x2-x1
     evalSpline _ (y':_) _ = y'
+    evalSpline _ _ _ = error "Empty spline list."
 {-# INLINE cSplineInterpolator #-}

@@ -1,12 +1,15 @@ module Quant.Math.Utilities (
     tdmaSolver
+  , cotraverseVec
 ) where
 
 import Control.Monad
 import Control.Monad.ST
 import qualified Data.Vector.Mutable as M
 import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
 
+-- |Tridiagonal matrix solver.  Pretty simple.
 tdmaSolver :: (Fractional a, Ord a) => [a] -> [a] -> [a] -> [a] -> [a]
 tdmaSolver aL bL cL dL = V.toList $ 
     let [a,b,c,d] = map V.fromList [aL,bL,cL,dL] in 
@@ -38,4 +41,11 @@                 xi1 <- M.read xn $ x+1
                 M.write xn x $ di - ci*xi1
             V.unsafeFreeze xn
-{-# INLINE tdmaSolver #-}+{-# INLINE tdmaSolver #-}
+
+-- | Something similar to cotraverse in the Distributive package,
+--but specialized to unboxed vectors, which are not functors.
+cotraverseVec :: (U.Unbox b1, U.Unbox b, Functor f) =>
+                       (f b1 -> b) -> Int -> f (U.Vector b1) -> U.Vector b
+cotraverseVec f l m = U.map (\i -> f (fmap (U.!i) m)) (U.enumFromN 0 l)
+{-# INLINE cotraverseVec #-}

@@ -1,16 +1,17 @@ {-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 
 module Quant.Models.Black (
     Black (..)
 ) where
 
+import Quant.Types
 import Quant.Time
 import Quant.YieldCurve
 import Data.Random
 import Control.Monad.State
 import Quant.MonteCarlo
-import Quant.ContingentClaim
 
 -- | 'Black' represents a Black-Scholes model.
 data Black = forall a b  . (YieldCurve a, YieldCurve b) => Black {
@@ -28,24 +29,22 @@             --vol' = vol :+ 0
             --logs = log s :+ 0
 
-instance Discretize Black where
-    initialize (Black s _ _ _)  = put (Observables [s], Time 0)
+instance Discretize Black Observables1 where
+    initialize (Black s _ _ _)  = put (Observables1 s, Time 0)
     {-# INLINE initialize #-}
 
     evolve' b@(Black _ vol _ _) t2 anti = do
-        (Observables (stateVal:_), t1) <- get
+        (Observables1 stateVal, t1) <- get
         fwd <- forwardGen b t2
-        let  grwth = (fwd - vol*vol/2) * timeDiff t1 t2
-        postVal <- do
-             resid <- lift stdNormal
-             if anti then
-                return $ stateVal * exp (grwth - resid*vol)
-             else
-                return $ stateVal * exp (grwth + resid*vol)
-        put (Observables [postVal], t2)
+        let t = timeDiff t1 t2 
+            grwth = (fwd - vol*vol/2) * t
+        resid <- lift stdNormal
+        let resid' = if anti then -resid else resid
+            postVal = stateVal * exp (grwth + resid'*vol*sqrt t)
+        put (Observables1 postVal, t2)
     {-# INLINE evolve' #-}
 
-    discount (Black _ _ _ dsc) t = disc dsc t
+    discount (Black _ _ _ dsc) t = return $ disc dsc t
     {-# INLINE discount #-}
 
     forwardGen (Black _ _ fg _) t2 = do

@@ -1,13 +1,14 @@ {-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 module Quant.Models.Dupire (
     Dupire (..)
 ) where
 
 import Quant.Time
+import Quant.Types
 import Data.Random
 import Control.Monad.State
-import Quant.ContingentClaim
 import Quant.MonteCarlo
 import Quant.YieldCurve
 
@@ -18,22 +19,23 @@  , mertonForwardGen  ::  a  -- ^ 'YieldCurve' to generate forwards
  , mertonDiscounter  ::  b } -- ^ 'YieldCurve' to generate discount rates
 
-instance Discretize Dupire where
-    initialize (Dupire s _ _ _) = put (Observables [s], Time 0)
+instance Discretize Dupire Observables1 where
+    initialize (Dupire s _ _ _) = put (Observables1 s, Time 0)
     {-# INLINE initialize #-}
 
     evolve' d@(Dupire _ f _ _) t2 anti = do
-        (Observables (stateVal:_), t1) <- get
+        (Observables1 stateVal, t1) <- get
         fwd <- forwardGen d t2
         let vol   = f t1 stateVal
-            grwth = (fwd - vol * vol / 2) * timeDiff t1 t2
+            t     = timeDiff t1 t2
+            grwth = (fwd - vol * vol / 2) * t
         normResid <- lift stdNormal
-        let s' | anti      = stateVal * exp (grwth - normResid*vol)
-               | otherwise = stateVal * exp (grwth - normResid*vol)
-        put (Observables [s'], t2)
+        let s' | anti      = stateVal * exp (grwth - normResid*vol*sqrt t)
+               | otherwise = stateVal * exp (grwth - normResid*vol*sqrt t)
+        put (Observables1 s', t2)
     {-# INLINE evolve' #-}
 
-    discount (Dupire _ _ _ dsc) t = disc dsc t
+    discount (Dupire _ _ _ dsc) t = return $ disc dsc t
     {-# INLINE discount #-}
 
     forwardGen (Dupire _ _ fg _) t2 = do

@@ -1,4 +1,5 @@ {-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 
 module Quant.Models.Heston (
@@ -6,11 +7,11 @@ ) where
 
 import Quant.Time
+import Quant.Types
 import Quant.YieldCurve
 import Data.Random
 import Control.Monad.State
 import Quant.MonteCarlo
-import Quant.ContingentClaim
 
 -- | 'Heston' represents a Heston model (i.e. stochastic volatility).
 data Heston = forall a b  . (YieldCurve a, YieldCurve b) => Heston {
@@ -23,12 +24,12 @@   , hestonForwardGen :: a       -- ^ 'YieldCurve' to generate forwards
   , hestonDisc       :: b }     -- ^ 'YieldCurve' to generate discounts
 
-instance Discretize Heston where
-    initialize (Heston s v0 _ _ _ _ _ _) = put (Observables [s, v0], Time 0)
+instance Discretize Heston Observables2 where
+    initialize (Heston s v0 _ _ _ _ _ _) = put (Observables2 s v0, Time 0)
     {-# INLINE initialize #-}
 
     evolve' h@(Heston _ _ vf l rho eta _ _) t2 anti = do
-        (Observables (sState:vState:_), t1) <- get
+        (Observables2 sState vState, t1) <- get
         fwd <- forwardGen h t2
         let grwth = (fwd - vState/2) * t
             t = timeDiff t1 t2
@@ -39,10 +40,10 @@           resid2 = rho * resid1 + sqrt (1-rho*rho) * resid2'
           v' = (sqrt vState `op` (eta/2.0*sqrt t* resid2))^(2 :: Int)-l*(vState-vf)*t-eta*eta*t/4.0
           s' = sState * exp (grwth `op` (resid1*sqrt (vState*t)))
-        put (Observables [s', abs v'], t2)
+        put (Observables2 s' (abs v'), t2)
     {-# INLINE evolve' #-}
 
-    discount (Heston _ _ _ _ _ _ _ d) t = disc d t
+    discount (Heston _ _ _ _ _ _ _ d) t = return $ disc d t
     {-# INLINE discount #-}
 
     forwardGen (Heston _ _ _ _ _ _ fg _) t2 = do

@@ -1,15 +1,16 @@ {-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 module Quant.Models.Merton (
     Merton (..)
 ) where
 
 import Quant.Time
+import Quant.Types
 import Data.Random
 import Data.Random.Distribution.Poisson
 import Control.Monad.State
 import Quant.MonteCarlo
-import Quant.ContingentClaim
 import Quant.YieldCurve
 
 -- | 'Merton' represents a Merton model (Black-Scholes w/ jumps).
@@ -30,27 +31,27 @@           --  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 _ _ _ _ _ _) = put (Observables [s], Time 0)
+instance Discretize Merton Observables1 where
+    initialize (Merton s _ _ _ _ _ _) = put (Observables1 s, Time 0)
     {-# INLINE initialize #-}
 
     evolve' m@(Merton _ vol intensity mu sig _ _) t2 anti = do
-        (Observables (stateVal:_), t1) <- get
+        (Observables1 stateVal, t1) <- get
         fwd <- forwardGen m t2
         let correction = exp (mu + sig*sig /2.0) - 1
             grwth = (fwd - vol*vol/2 - intensity * correction) * t
             t = timeDiff t1 t2
         normResid1 <- lift stdNormal
         normResid2 <- lift stdNormal
-        poissonResid <- lift $ integralPoisson (intensity * t) :: MonteCarlo (MCObservables, Time) Int
+        poissonResid <- lift $ integralPoisson (intensity * t) :: MonteCarlo (Observables1, Time) Int
         let  poisson' = fromIntegral poissonResid
              jumpterm = mu*poisson'+sig*sqrt poisson' * normResid2
-             s' | anti      = stateVal * exp (grwth - normResid1*vol + jumpterm)
-                | otherwise = stateVal * exp (grwth + normResid1*vol + jumpterm)
-        put (Observables [s'], t2)
+             s' | anti      = stateVal * exp (grwth - normResid1*vol*sqrt t + jumpterm)
+                | otherwise = stateVal * exp (grwth + normResid1*vol*sqrt t + jumpterm)
+        put (Observables1 s', t2)
     {-# INLINE evolve' #-}
 
-    discount (Merton _ _ _ _ _ _ dsc) t = disc dsc t
+    discount (Merton _ _ _ _ _ _ dsc) t = return $ disc dsc t
     {-# INLINE discount #-}
 
     forwardGen (Merton _ _ _ _ _ fg _) t2 = do

@@ -1,14 +1,14 @@ 
 module Quant.Models.Processes (
     ProcessSpec (..)
+  , normal
   , lognormal
 ) where
 
-data ProcessSpec = ProcessSpec {
-    procInit    :: Double
-  , procGrowth  :: Double
-  , procElapsed :: Double
-}
+data ProcessSpec = ProcessSpec {-# UNPACK #-} !Double !Double !Double
+
+normal :: ProcessSpec -> Double -> Double -> Double
+normal (ProcessSpec initVal r t) vol normRand = initVal + vol * normRand + r * t
 
 lognormal :: ProcessSpec -> Double -> Double -> Double 
 lognormal (ProcessSpec initVal r t) vol normRand = initVal * exp ( g + sig * normRand )

@@ -1,4 +1,6 @@ {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
 
 
 module Quant.MonteCarlo (
@@ -11,6 +13,11 @@   , Discretize(..)
   , OptionType(..)
 
+  , runSimulation
+  , runSimulationAnti
+  , quickSim
+  , quickSimAnti
+
   )
 where
 
@@ -45,17 +52,17 @@ 
 Minimal complete definition: 'initialize', 'discounter', 'forwardGen' and 'evolve''.
 -}
-class Discretize a where
+class Discretize a b | a -> b where
 
     -- | Initializes a Monte Carlo simulation for a given number of runs.
-    initialize :: Discretize a => a   -- ^ Model
-                               -> MonteCarlo (MCObservables, Time) ()
+    initialize :: Discretize a b => a   -- ^ Model
+                                 -> MonteCarlo (b, Time) ()
 
     -- | Evolves the internal states of the MC variables between two times.
-    evolve :: Discretize a => a            -- ^ Model
-                           -> Time         -- ^ time to evolve to
-                           -> Bool         -- whether or not to use flipped variates
-                           -> MonteCarlo (MCObservables, Time) ()
+    evolve :: Discretize a b => a            -- ^ Model
+                             -> Time         -- ^ time to evolve to
+                             -> Bool         -- whether or not to use flipped variates
+                             -> MonteCarlo (b, Time) ()
     evolve mdl t2 anti = do
         (_, t1) <- get
         let ms = maxStep mdl
@@ -66,35 +73,30 @@               evolve' mdl (timeOffset t1 ms) anti
               evolve mdl t2 anti
 
-    -- | Stateful discounting function, takes a model and a time, and returns a vector of results.
-    discountState :: Discretize a => a -> Time -> MonteCarlo (MCObservables, Time) Double
-    discountState m t = return $ discount m t
-    {-# INLINE discountState #-}
-
     -- | Non-stateful discounting function...might need to find a better place to put this.
-    discount :: Discretize a => a -> Time -> Double
+    discount :: Discretize a b => a -> Time -> MonteCarlo (b, Time) Double
 
     -- | Stateful forward generator for a given model at a certain time.
-    forwardGen :: Discretize a => a -> Time -> MonteCarlo (MCObservables, Time) Double
+    forwardGen :: Discretize a b => a -> Time -> MonteCarlo (b, Time) Double
 
     -- | Internal function to evolve a model to a given time.
-    evolve' :: Discretize a => a          -- ^ model
-                            -> Time       -- ^ time to evolve to
-                            -> Bool       -- ^ whether or not to use flipped variates
-                            -> MonteCarlo (MCObservables, Time) () -- ^ computation result
+    evolve' :: Discretize a b => a          -- ^ model
+                              -> Time       -- ^ time to evolve to
+                              -> Bool       -- ^ whether or not to use flipped variates
+                              -> MonteCarlo (b, Time) () -- ^ computation result
 
     -- | Determines the maximum size time-step for discretization purposes. Defaults to 1/250.
-    maxStep :: Discretize a => a -> Double
+    maxStep :: Discretize a b => a -> Double
     maxStep _ = 1/250
     {-# INLINE maxStep #-}
 
     -- | Perform a simulation of a compiled basket of contingent claims.
-    simulateState :: Discretize a => 
-           a                       -- ^ model
-        -> ContingentClaim         -- ^ compilied basket of claims
-        -> Int                     -- ^ number of trials
-        -> Bool                    -- ^ antithetic?
-        -> MonteCarlo (MCObservables, Time) Double -- ^ computation result
+    simulateState :: Discretize a b => 
+           a                              -- ^ model
+        -> ContingentClaim b              -- ^ compilied basket of claims
+        -> Int                            -- ^ number of trials
+        -> Bool                           -- ^ antithetic?
+        -> MonteCarlo (b, Time) Double -- ^ computation result
     simulateState modl (ContingentClaim ccb) trials anti = avg <$> replicateM trials singleTrial
           where 
             singleTrial = initialize modl >> 
@@ -104,31 +106,27 @@             process discCFs obsMap c@(CCProcessor t mf:ccs) allcfs@(CashFlow cft amt:cfs) = 
               if t > cft then do
                   evolve modl cft anti
-                  d <- discountState modl cft
+                  d <- discount modl cft
                   process (discCFs+d*amt) obsMap c cfs
               else do
                   evolve modl t anti
                   obs <- gets fst
                   let obsMap' = Map.insert t obs obsMap
-                  case mf of
-                    Nothing -> process discCFs obsMap' ccs allcfs
-                    Just f -> let newCFs = map ($obsMap') f
-                                  insertCFList xs cfList = foldl' (flip insertCF) cfList xs in
-                        process discCFs obsMap' ccs (insertCFList newCFs allcfs)
+                      newCFs = map ($obsMap') mf
+                      insertCFList xs cfList = foldl' (flip insertCF) cfList xs
+                  process discCFs obsMap' ccs (insertCFList newCFs allcfs)
 
             process discCFs obsMap (CCProcessor t mf:ccs) [] = do
               evolve modl t anti
               obs <- gets fst
               let obsMap' = Map.insert t obs obsMap
-              case mf of
-                Nothing -> process discCFs obsMap' ccs []
-                Just f -> let newCFs = map ($obsMap') f
-                              insertCFList xs cfList = foldl' (flip insertCF) cfList xs in
-                        process discCFs obsMap' ccs (insertCFList newCFs [])                       
+                  newCFs = map ($obsMap') mf
+                  insertCFList xs cfList = foldl' (flip insertCF) cfList xs
+              process discCFs obsMap' ccs (insertCFList newCFs [])
 
             process discCFs obsMap [] (cf:cfs) = do
               evolve modl (cfTime cf) anti
-              d <- discountState modl $ cfTime cf
+              d <- discount modl $ cfTime cf
               process (discCFs+d*cfAmount cf) obsMap [] cfs
 
             process discCFs _ _ _ = return $! discCFs
@@ -140,30 +138,34 @@ 
             avg v = sum v / fromIntegral trials
 
-    -- | 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 [], Time 0)
-       where
-            run = simulateState modl ccs trials anti
+-- | Runs a simulation for a 'ContingentClaim'.
+runSimulation :: (Discretize a b,
+                         MonadRandom (StateT c Identity)) =>
+                            a                             -- ^ model
+                         -> ContingentClaim b             -- ^ claims to value
+                         -> c                             -- ^ initial random state
+                         -> Int                           -- ^ trials
+                         -> Bool                          -- ^ whether to use antithetic variables
+                         -> Double                        -- ^ final value
+runSimulation modl ccs seed trials anti = runMC run seed (undefined, Time 0)
+   where
+        run = simulateState modl ccs trials anti
+{-# INLINE runSimulation #-}
 
-    -- | 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)
+-- | Like 'runSimulation', but splits the trials in two and does antithetic variates.
+runSimulationAnti :: (Discretize a b,
+                         MonadRandom (StateT c Identity)) =>
+                        a -> ContingentClaim b -> c -> Int -> Double
+runSimulationAnti modl ccs seed trials = (runSim True + runSim False) / 2
+    where runSim x = runSimulation modl ccs seed (trials `div` 2) x
+{-# INLINE runSimulationAnti #-}
 
-    -- | '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
+-- | 'runSimulation' with a default random number generator.
+quickSim :: Discretize a b => a -> ContingentClaim b -> Int -> Double
+quickSim mdl opts trials = runSimulation mdl opts (pureMT 500) trials False
+{-# INLINE quickSim #-}
 
-    -- | 'runSimulationAnti' with a default random number generator.
-    quickSimAnti :: Discretize a => a -> ContingentClaim -> Int -> Double
-    quickSimAnti mdl opts trials = runSimulationAnti mdl opts (pureMT 500) trials+-- | 'runSimulationAnti' with a default random number generator.
+quickSimAnti :: Discretize a b => a -> ContingentClaim b -> Int -> Double
+quickSimAnti mdl opts trials = runSimulationAnti mdl opts (pureMT 500) trials
+{-# INLINE quickSimAnti #-}

@@ -0,0 +1,88 @@+module Quant.RNG.MWC64X (
+    MWC64X(..)
+  , randomWord32
+  , randomWord64
+  , randomInt
+  , randomDouble
+  , randomInt64
+  , skip
+) where
+
+import Data.Int
+import Data.Bits 
+import Data.Word 
+import Data.Random.Internal.Words
+import System.Random
+
+data MWC64X = MWC64X {-# UNPACK #-} !Word64 deriving (Eq,Show)
+
+randomWord32 :: MWC64X -> (Word32, MWC64X)
+randomWord32 (MWC64X state) = (x `xor` c, MWC64X state')
+    where c = fromIntegral $ state `shiftR` 32 :: Word32 
+          x = fromIntegral $ state .&. 0xFFFFFFFF :: Word32 
+          state' = fromIntegral x * aConst + fromIntegral c 
+
+randomInt :: MWC64X -> (Int,MWC64X)
+randomInt g = (fromIntegral i, g')
+        where (i, g') = randomWord64 g
+
+randomWord64 :: MWC64X -> (Word64, MWC64X) 
+randomWord64 x = (buildWord64'' y1 y2, x'') 
+    where (y1, x' )  = randomWord32 x 
+          (y2, x'') = randomWord32 x'
+
+randomDouble :: MWC64X -> (Double, MWC64X)
+randomDouble x = (fromIntegral (val `div` 2048) / 9007199254740992, x')
+    where (val, x') = randomWord64 x
+
+randomInt64 :: MWC64X -> (Int64,MWC64X)
+randomInt64 g = (fromIntegral i, g')
+        where (i, g') = randomWord64 g
+
+instance RandomGen MWC64X where
+    next g = (fromIntegral w, g')
+        where (w, g') = randomWord64 g
+    split g = (skip g skipConst, g)
+
+addMod64 :: Word64 -> Word64 -> Word64 -> Word64
+addMod64 a b m = (a+b) `mod` m
+
+--mulMod64 :: Word64 -> Word64 -> Word64 -> Word64
+--mulMod64 a b m = (a * b) `mod` m
+mulMod64 :: Word64 -> Word64 -> Word64 -> Word64
+mulMod64 a b m = f 0 a b
+    where f r a1 b1
+            | a1 == 0 = r
+            | otherwise = f r' a' b'
+                where r' = if a1 .&. 1 == 1 then addMod64 r b1 m else r
+                      b' = addMod64 b1 b1 m
+                      a' = a `shiftR` 1
+
+powMod64 :: Word64 -> Word64 -> Word64 -> Word64
+powMod64 a e m = f a 1 e
+    where f sqr acc e1
+            | e1 == 0 = acc
+            | otherwise = f sqr' acc' e' 
+                where acc' = if e1 .&. 1 == 1 then mulMod64 acc sqr m else acc
+                      sqr' = mulMod64 sqr sqr m
+                      e'   = e `shiftR` 1
+
+skip :: MWC64X -> Word64 -> MWC64X
+skip (MWC64X st) d = MWC64X st'
+    where
+        m   = powMod64 aConst d mConst
+        c   = st `shiftR` 32 
+        x   = st .&. 0xFFFFFFFF 
+        x'  = mulMod64 (x * aConst + c) m mConst
+        x'' = fromIntegral $ x' `div` aConst :: Word32
+        c'  = fromIntegral $ x' `mod` aConst :: Word32
+        st' = buildWord64'' c' x''
+
+mkWord64 :: Word32 -> Word32 -> Word64
+mkWord64 a b = (fromIntegral $ a `shiftL` 32) .&. fromIntegral b
+
+aConst, mConst, skipConst :: Word64
+aConst = 4294883355
+--bConst = 4077358422479273989
+mConst = 18446383549859758079
+skipConst = 1073741824

@@ -6,7 +6,7 @@   , timeFromZero
   ) where
 
-data Time = Time Double deriving (Eq,Show,Ord)
+data Time = Time {-# UNPACK #-} !Double deriving (Eq,Show,Ord)
 
 timeDiff :: Time -> Time -> Double
 timeDiff (Time x) (Time y) = y - x

@@ -1,24 +1,117 @@ 
 module Quant.Types (
     CashFlow(..)
-  , Observables(..)
-  , MCObservables
+  , Observables1(..)
+  , Observables2(..)
+  , Observables3(..)
+  , Observables4(..)
+  , Observables5(..)
   , OptionType(..)
+  , Obs1(..)
+  , Obs2(..)
+  , Obs3(..)
+  , Obs4(..)
+  , Obs5(..)
   ) where
 
 import Quant.Time
 
+-- | A CashFlow is just a time and an amount.
 data CashFlow = CashFlow {
     cfTime   :: Time
   , cfAmount :: Double
 }
 
--- | 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 a = Observables { obsGet :: [a] } deriving (Show)
 
-type MCObservables = Observables Double
-
 -- | Type for Put or Calls
 data OptionType = Put | Call deriving (Eq,Show)
+
+-- | Single-observable container.
+data Observables1 = Observables1 {-# UNPACK #-} !Double
+-- | Two observable container.
+data Observables2 = Observables2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+-- | Three observable container.
+data Observables3 = Observables3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double 
+                                 {-# UNPACK #-} !Double
+-- | Four observable container.
+data Observables4 = Observables4 {-# UNPACK #-} !Double {-# UNPACK #-} !Double 
+                                 {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+-- | Five observable container.
+data Observables5 = Observables5 {-# UNPACK #-} !Double {-# UNPACK #-} !Double 
+                                 {-# UNPACK #-} !Double {-# UNPACK #-} !Double 
+                                 {-# UNPACK #-} !Double
+
+class Obs1 a where
+    get1 :: a -> Double
+
+class (Obs1 a) => Obs2 a where
+	get2 :: a -> Double
+
+class (Obs2 a) => Obs3 a where
+	get3 :: a -> Double
+
+class (Obs3 a) => Obs4 a where
+	get4 :: a -> Double
+
+class (Obs4 a) => Obs5 a where
+	get5 :: a -> Double
+
+instance Obs1 Observables1 where
+	get1 (Observables1 x) = x
+	{-# INLINE get1 #-}
+
+instance Obs1 Observables2 where
+	get1 (Observables2 x _) = x
+	{-# INLINE get1 #-}
+
+instance Obs1 Observables3 where
+	get1 (Observables3 x _ _) = x
+	{-# INLINE get1 #-}
+
+instance Obs1 Observables4 where
+	get1 (Observables4 x _ _ _) = x
+	{-# INLINE get1 #-}
+
+instance Obs1 Observables5 where
+	get1 (Observables5 x _ _ _ _) = x
+	{-# INLINE get1 #-}
+
+instance Obs2 Observables2 where
+	get2 (Observables2 _ x) = x
+	{-# INLINE get2 #-}
+
+instance Obs2 Observables3 where
+	get2 (Observables3 _ x _) = x
+	{-# INLINE get2 #-}
+
+instance Obs2 Observables4 where
+	get2 (Observables4 _ x _ _) = x
+	{-# INLINE get2 #-}
+
+instance Obs2 Observables5 where
+	get2 (Observables5 _ x _ _ _) = x
+	{-# INLINE get2 #-}
+
+instance Obs3 Observables3 where
+	get3 (Observables3 _ _ x) = x
+	{-# INLINE get3 #-}
+
+instance Obs3 Observables4 where
+	get3 (Observables4 _ _ x _) = x
+	{-# INLINE get3 #-}
+
+instance Obs3 Observables5 where
+	get3 (Observables5 _ _ x _ _) = x
+	{-# INLINE get3 #-}
+
+instance Obs4 Observables4 where
+	get4 (Observables4 _ _ _ x) = x
+	{-# INLINE get4 #-}
+
+instance Obs4 Observables5 where
+	get4 (Observables5 _ _ _ x _) = x
+	{-# INLINE get4 #-}
+
+instance Obs5 Observables5 where
+	get5 (Observables5 _ _ _ _ x) = x
+	{-# INLINE get5 #-}

@@ -0,0 +1,77 @@+
+module Quant.VectorOps (
+    (.*.)
+  , (.*)
+  , (*.)
+  , (./.)
+  , (/.)
+  , (./)
+  , (.+.)
+  , (.+)
+  , (+.)
+  , (.-.)
+  , (.-)
+  , (-.)
+  ) where
+
+import qualified Data.Vector.Unboxed as U
+
+infixl 7 .*.
+(.*.) :: (U.Unbox a, Num a) => U.Vector a -> U.Vector a -> U.Vector a
+x .*. y = U.zipWith (*) x y
+{-# INLINE (.*.) #-}
+
+infixl 7 *.
+(*.) :: (U.Unbox a, Num a) => a -> U.Vector a -> U.Vector a
+x *. y = U.map (x*) y
+{-# INLINE (*.) #-}
+
+infixl 7 .*
+(.*) :: (U.Unbox a, Num a) => U.Vector a -> a -> U.Vector a
+x .* y = U.map (*y) x
+{-# INLINE (.*) #-}
+
+infixl 7 ./.
+(./.) :: (U.Unbox a, Fractional a) => U.Vector a -> U.Vector a -> U.Vector a
+x ./. y = U.zipWith (/) x y
+{-# INLINE (./.) #-}
+
+infixl 7 /.
+(/.) :: (U.Unbox a, Fractional a) => a -> U.Vector a -> U.Vector a
+x /. y = U.map (x/) y
+{-# INLINE (/.) #-}
+
+infixl 7 ./
+(./) :: (U.Unbox a, Fractional a) => U.Vector a -> a -> U.Vector a
+x ./ y = U.map (/y) x
+{-# INLINE (./) #-}
+
+infixl 6 .+.
+(.+.) :: (U.Unbox a, Num a) => U.Vector a -> U.Vector a -> U.Vector a
+x .+. y = U.zipWith (+) x y
+{-# INLINE (.+.) #-}
+
+infixl 6 +.
+(+.) :: (U.Unbox a, Num a) => a -> U.Vector a -> U.Vector a
+x +. y = U.map (x+) y
+{-# INLINE (+.) #-}
+
+infixl 6 .+
+(.+) :: (U.Unbox a, Num a) => U.Vector a -> a -> U.Vector a
+x .+ y = U.map (+y) x
+{-# INLINE (.+) #-}
+
+infixl 6 .-.
+(.-.) :: (U.Unbox a, Num a) => U.Vector a -> U.Vector a -> U.Vector a
+x .-. y = U.zipWith (-) x y
+{-# INLINE (.-.) #-}
+
+infixl 6 -.
+(-.) :: (U.Unbox a, Num a) => a -> U.Vector a -> U.Vector a
+x -. y = U.map (x-) y
+{-# INLINE (-.) #-}
+
+infixl 6 .-
+(.-) :: (U.Unbox a, Num a) => U.Vector a -> a -> U.Vector a
+x .- y = U.map (+negate y) x
+{-# INLINE (.-) #-}

@@ -57,7 +57,7 @@                 solution = dwdt/(1.0-y/w*dwdy+0.25*(-0.25-1.0/w+y*y/w/w)*dwdy*dwdy+0.5*d2wdy2)
 
 -- |A flat surface has one volatility at all times/maturities.
-data FlatSurf = FlatSurf Double
+data FlatSurf = FlatSurf {-# UNPACK #-} !Double
 
 instance VolSurf FlatSurf where
     vol (FlatSurf x) _ _ = x

@@ -27,13 +27,15 @@ 
 -- |A flat curve is just a flat curve with one continuously 
 -- compounded rate at all points on the curve.
-data FlatCurve = FlatCurve Double
+data FlatCurve = FlatCurve {-# UNPACK #-} !Double
 
 instance YieldCurve FlatCurve where
     disc (FlatCurve r) t = exp (-r*timeFromZero t)
+    {-# INLINE disc #-}
 
 -- | '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+    disc (NetYC yc1 yc2) t = disc yc1 t / disc yc2 t
+    {-# INLINE disc #-}