quantfin 0.1.0.1 → 0.1.0.2

raw patch · 18 files changed

+587/−374 lines, 18 filesdep +quantfindep +randomdep +random-sourcedep ~basenew-component:exe:examplePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: quantfin, random, random-source

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• Test.hs +96/−0
• quantfin.cabal +22/−6
• src/Quant/ContingentClaim.hs +106/−89
• src/Quant/Math/Integration.hs +4/−1
• src/Quant/Math/Interpolation.hs +59/−0
• src/Quant/Math/Utilities.hs +41/−0
• src/Quant/Models.hs +0/−43
• src/Quant/Models/Black.hs +14/−14
• src/Quant/Models/Dupire.hs +16/−19
• src/Quant/Models/Heston.hs +19/−25
• src/Quant/Models/Merton.hs +19/−21
• src/Quant/Models/Processes.hs +20/−0
• src/Quant/MonteCarlo.hs +77/−68
• src/Quant/Test.hs +0/−63
• src/Quant/Time.hs +21/−0
• src/Quant/Types.hs +24/−0
• src/Quant/VolSurf.hs +40/−19
• src/Quant/YieldCurve.hs +9/−6

@@ -0,0 +1,96 @@+
+import Quant.Time
+import Data.Monoid
+import Quant.MonteCarlo
+import Quant.YieldCurve
+import Quant.ContingentClaim
+import Quant.Models.Black
+import Quant.Models.Heston
+
+--create a flat yield curve with a 5% rate
+baseYC :: FlatCurve
+baseYC = FlatCurve 0.05 
+
+black :: Black
+black = Black 
+			100     --initial stock price
+			0.2     --volatility
+			baseYC  --forward generator
+			baseYC  --discount function
+
+--make a vanilla put, struck at 100, maturing at time 1
+vanopt :: ContingentClaim
+vanopt = vanillaOption Call 100 (Time 1) --built in function
+
+vanopt' :: ContingentClaim
+vanopt' = specify $ do
+	x <- monitor (Time 1)
+	return $ CashFlow (Time 1) (max (x - 100) 0) --roll your own
+
+--Run a Monte Carlo on opt in a a black model with 10000 trials
+vanoptPrice :: Double
+vanoptPrice = quickSim black vanopt 100000 
+
+--Make a call spread with a 100 unit notional, using some handy combinators.
+cs :: ContingentClaim
+cs =  multiplier 100 
+   $  vanillaOption Call 100 (Time 1) 
+   <> short (vanillaOption Call 120 (Time 1)) 
+
+--Run a Monte Carlo on the call spread; use antithetic variates
+csPrice :: Double
+csPrice = quickSim black' cs 100000 
+
+black' :: Black
+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
+
+callSpreadAnti :: Double
+callSpreadAnti = quickSimAnti black' cs 100000
+
+--Let's try it with a Heston model
+heston :: Heston
+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
+
+--price the call spread in the Heston model
+csHeston :: Double
+csHeston = quickSimAnti heston cs 100000
+
+--create an option that pays off based on the square of its underlying
+squareOpt :: ContingentClaim
+squareOpt = terminalOnly (Time 1) $ \x -> x*x  --using the built in function
+squareOpt' = specify $ do --roll your own
+	x <- monitor (Time 1)
+	return $ CashFlow (Time 1) $ x*x
+squareOptPrice :: Double
+squareOptPrice = quickSimAnti black squareOpt 100000
+
+--create an option with a bizarre payoff
+bizarre :: ContingentClaim
+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
+bizarrePrice :: Double
+bizarrePrice = quickSimAnti black bizarre 100000
+
+main :: IO ()
+main = do
+	print vanoptPrice
+	print csPrice
+	print callSpreadAnti
+	print csHeston
+	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.1
+version:             0.1.0.2
 
 -- A short (one-line) description of the package.
 synopsis:            Quant finance library in pure Haskell.
@@ -37,7 +37,7 @@ -- A copyright notice.
 -- copyright:           
 
-category:            Quant
+category:            Finance
 
 build-type:          Simple
 
@@ -54,14 +54,19 @@   exposed-modules:     Quant.YieldCurve
                        Quant.VolSurf
                        Quant.Math.Integration
-                       Quant.Models
+                       Quant.Math.Interpolation
+                       Quant.Math.Utilities
+                       --Quant.Models
                        Quant.Models.Black
                        Quant.Models.Merton
                        Quant.Models.Dupire
                        Quant.Models.Heston
+                       Quant.Models.Processes
+                       --Quant.RNG.MWC64X
                        Quant.MonteCarlo
                        Quant.ContingentClaim
-                       Quant.Test
+                       Quant.Types
+                       Quant.Time
   
   -- Modules included in this library but not exported.
   -- other-modules:       
@@ -77,11 +82,22 @@                        random-fu,
                        containers,
                        rvar,
-                       mersenne-random-pure64
+                       mersenne-random-pure64,
+                       random-source,
+                       random
 
   -- Directories containing source files.
   hs-source-dirs:      src
   
   -- Base language which the package is written in.
   default-language:    Haskell2010
-  +  
+  ghc-options:         -Wall -fexcess-precision
+
+executable "example"
+  build-depends:
+    base, 
+    quantfin
+  main-is: Test.hs
+  ghc-options: -Wall -fexcess-precision -rtsopts
+  default-language: Haskell2010

@@ -1,13 +1,17 @@ module Quant.ContingentClaim (
     -- * Types for modeling contingent claims.
-    ContingentClaim
-  , ContingentClaim' (..)
+    ContingentClaim (..)
+  , CCProcessor (..)
   , Observables (..)
-  , ContingentClaimBasket (..)
+  , MCObservables
   , OptionType (..)
-  , ccBasket
+  , CashFlow (..)
+  , CCBuilder
 
   -- * Options and option combinators
+  , specify
+  , monitor
+  , monitorByNum
   , vanillaOption
   , binaryOption
   , straddle
@@ -16,46 +20,59 @@   , callSpread
   , putSpread
   , forwardContract
-  , fixed
+  , zcb
+  , fixedBond
   , multiplier
   , short
   , combine
   , terminalOnly
-  , changeObservableFct
 
-  -- * Utility functions
-  , obsNum
-  , obsHead
-        )  where
+)  where
 
-import Data.List
-import Data.Ord
-import qualified Data.Vector.Unboxed as U
+import Control.Monad.Reader
+import Control.Monad.Writer.Strict
+import Quant.Types
+import Quant.Time
+import qualified Data.Map as M
 
+type MCMap = M.Map Time MCObservables
+type PayoffFunc a = MCMap -> a
 
--- | '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 
-  , observations :: [( Double                            
-                     , Observables  -> U.Vector Double   
-                     , Double       -> Double) ]         {- ^ List containing:
-                                                         -- Time of observation, 
-                                                         -- Function to access specific observable, 
-                                                         -- Function to collect observations and transform them into a cash flow. -}
+
+data CCProcessor = CCProcessor  {
+                      monitorTime      :: Time
+                    , payoutFunc       :: Maybe [PayoffFunc CashFlow]
 }
 
--- | '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)
+type CCBuilder w r a = WriterT w (Reader r) a
 
--- | ADT for Put or Calls
-data OptionType = Put | Call deriving (Eq,Show)
+monitor :: Time -> CCBuilder ContingentClaim MCMap Double
+monitor = monitorByNum 0
 
+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.
+
+specify :: CCBuilder ContingentClaim MCMap CashFlow -> ContingentClaim
+specify x = w `mappend` ContingentClaim [CCProcessor (last0 w') (Just [f])]
+  where
+    w  = runReader (execWriterT x) M.empty
+    f  = runReader . liftM fst $ runWriterT x
+    w' = 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
@@ -77,94 +94,94 @@ 
 -- | 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)]]
+terminalOnly :: Time -> (Double -> Double) -> ContingentClaim
+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 -> Double -> ContingentClaim
+vanillaOption :: OptionType -> Double -> Time -> 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 :: OptionType -> Double -> Double -> Time -> 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
+arithmeticAsianOption :: OptionType -> Double -> [Time] -> Time -> ContingentClaim
+arithmeticAsianOption pc strike obsTimes t = specify $ do
+  x <- mapM monitor obsTimes
+  let avg = sum x / fromIntegral (length obsTimes)
+  return $ CashFlow t $ vanillaPayout pc strike avg
 
 -- | 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
+--maturity and generates an arithmetic Asian option.
+geometricAsianOption :: OptionType -> Double -> [Time] -> Time -> ContingentClaim
+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 notional cs = map f cs
-    where f c@(ContingentClaim' _ collFct _) = c { collector = U.map (*notional) . collFct }
+multiplier notional cs = ContingentClaim $ map f (unCC cs)
+    where f (CCProcessor t g) = CCProcessor t $ fmap (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 = 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
+zcb :: Time -> Double -> ContingentClaim
+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 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 :: Double -> ContingentClaim
-forwardContract t = terminalOnly t id
+forwardContract :: Time -> ContingentClaim
+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 -> Double -> ContingentClaim
-callSpread lowStrike highStrike t = combine (vanillaOption Call lowStrike t) (short $ vanillaOption Call highStrike t)
+callSpread :: Double -> Double -> Time -> ContingentClaim
+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 -> Double -> ContingentClaim
-putSpread lowStrike highStrike t = combine (vanillaOption Put highStrike t) (short $ vanillaOption Put lowStrike t)
+putSpread :: Double -> Double -> Time -> ContingentClaim
+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 -> Double -> ContingentClaim
-straddle strike t = vanillaOption Put strike t ++ vanillaOption Call strike t
+straddle :: Double -> Time -> ContingentClaim
+straddle strike t = vanillaOption Put strike t <> vanillaOption Call strike t
 
--- | Just combines two contingent claims into one. 
+-- | 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+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
+      | monitorTime cc1 > monitorTime cc2 = cc2 : combine' (cc1:ccs1) ccs2
+      | otherwise = cc1 : combine' ccs1 (cc2:ccs2)
+    combine' [] [] = []
+    combine' cs [] = cs
+    combine' [] cs = cs

@@ -14,6 +14,7 @@     where
         dx = (uBound - lBound) / fromIntegral intervals
         points = take intervals $ iterate (+dx) (lBound+dx/2)
+{-# INLINE midpoint #-}
 
 -- | Trapezoidal integration.
 trapezoid :: Int -> Integrator
@@ -21,7 +22,9 @@     where
         dx = (uBound - lBound) / fromIntegral intervals
         points = take (intervals-1) $ iterate (+dx) (lBound+dx)
+{-# INLINE trapezoid #-}
 
 -- | 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+simpson intervals f l u =( 2 * midpoint intervals f l u  + trapezoid intervals f l u ) / 3
+{-# INLINE simpson #-}

@@ -0,0 +1,59 @@+module Quant.Math.Interpolation (
+    linearInterpolator
+  , logLinearInterpolator
+  , linearVarianceInterpolator
+  , cSplineInterpolator
+  , Interpolator1d
+  --, cubicSpline
+) where
+
+import Quant.Math.Utilities (tdmaSolver)
+import Data.List (zipWith5)
+
+type Interpolator1d = [Double] -> [Double] -> Double -> Double
+
+linearInterpolator :: Interpolator1d
+linearInterpolator (x1:x2:xs) (y1:y2:ys) x
+    | x >= x2   = linearInterpolator (x2:xs) (y2:ys) x
+    | x <= x1   = y1
+    | otherwise = wt1 * y1 + wt2 * y2
+        where
+          wt1 = (x2-x) / (x2-x1)
+          wt2 = (x-x1) / (x2-x1)
+linearInterpolator _ [y] _ = y
+{-# INLINE linearInterpolator #-}
+
+logLinearInterpolator :: Interpolator1d
+logLinearInterpolator x1 x2 x = exp $ linearInterpolator x1 (map log x2) x
+{-# INLINE logLinearInterpolator #-}
+
+linearVarianceInterpolator :: Interpolator1d
+linearVarianceInterpolator xs ys = linearInterpolator xs 
+                                 . map (\(x, y) -> y*y*x) 
+                                 $ zip xs ys
+{-# INLINE linearVarianceInterpolator #-}
+
+cSplineInterpolator :: Interpolator1d
+cSplineInterpolator xs ys x = evalSpline xs ys moments
+  where
+    h = zipWith (-) (tail xs) (init xs)
+    lambda = 0.0 : zipWith (\a b->a/(a+b)) (tail h) (init h)
+    mu = map (\a->1.0-a) lambda++[0.0]
+    dj hj hj1 yj1 yj yjm1= 6.0/(hj+hj1)*((yj1-yj)/hj1-(yj-yjm1)/hj)
+    d = 0.0 : zipWith5 dj (init h) (tail h) 
+        (tail $ tail ys) (tail $ init ys) 
+        (init $ init ys) ++ [0.0] 
+    moments = tdmaSolver mu (replicate (length d) 2) lambda d
+    evalSpline (x1:x2:xs') (y1:y2:ys') (m1:m2:ms)
+      | x <= x1   = y1
+      | x >= x2   = evalSpline (x2:xs') (y2:ys') (m2:ms)
+      | otherwise = y1+beta*term+gamma*term*term+
+                      delta*term*term*term
+          where
+            gamma = m1/2.0
+            beta = (y2-y1)/h'-(2*m1+m2)*h'/6.0
+            delta = (m2-m1)/6.0/h'
+            term = x-x1
+            h' = x2-x1
+    evalSpline _ (y':_) _ = y'
+{-# INLINE cSplineInterpolator #-}

@@ -0,0 +1,41 @@+module Quant.Math.Utilities (
+    tdmaSolver
+) where
+
+import Control.Monad
+import Control.Monad.ST
+import qualified Data.Vector.Mutable as M
+import qualified Data.Vector as V
+
+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 
+        runST $ do
+            c' <- V.thaw c
+            M.write c' 0 (V.head c / V.head b)
+            forM_ [1..V.length c-1] $ \x -> do
+                let ai = a V.! x
+                    bi = b V.! x
+                    ci = c V.! x
+                ci1' <- M.read c' (x-1)
+                M.write c' x $ ci / (bi-ai*ci1')
+            cf <- V.unsafeFreeze c'
+            d' <- V.thaw d
+            M.write d' 0 (V.head d / V.head b)
+            forM_ [1..V.length d-1] $ \x -> do
+                let ai  = a  V.! x
+                    bi  = b  V.! x
+                    di  = d  V.! x
+                    ci1 = cf V.! (x-1)
+                di1' <- M.read d' (x-1)
+                M.write d' x $ (di-ai*di1') / (bi-ai*ci1)
+            df <- V.unsafeFreeze d'
+            xn <- M.new $ V.length d
+            M.write xn (V.length d-1) $ V.last df
+            forM_ (reverse [0..V.length df-2]) $ \ x-> do
+                let ci = cf V.! x
+                    di = df V.! x
+                xi1 <- M.read xn $ x+1
+                M.write xn x $ di - ci*xi1
+            V.unsafeFreeze xn
+{-# INLINE tdmaSolver #-}

@@ -1,43 +0,0 @@-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

@@ -1,17 +1,16 @@ {-# LANGUAGE ExistentialQuantification #-}
-{-# LANGUAGE FlexibleInstances #-}
 
 
 module Quant.Models.Black (
     Black (..)
 ) where
 
+import Quant.Time
 import Quant.YieldCurve
 import Data.Random
 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 {
@@ -30,27 +29,28 @@             --logs = log s :+ 0
 
 instance Discretize Black where
-    initialize (Black s _ _ _) trials = put (Observables [U.replicate trials s], 0)
+    initialize (Black s _ _ _)  = put (Observables [s], Time 0)
+    {-# INLINE initialize #-}
 
     evolve' b@(Black _ vol _ _) t2 anti = do
-        (Observables (stateVec:_), t1) <- get
+        (Observables (stateVal:_), 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
+        let  grwth = (fwd - vol*vol/2) * timeDiff t1 t2
+        postVal <- do
              resid <- lift stdNormal
              if anti then
-                return $ x * exp (g - resid*vol)
+                return $ stateVal * exp (grwth - resid*vol)
              else
-                return $ x * exp (g + resid*vol)
+                return $ stateVal * exp (grwth + resid*vol)
         put (Observables [postVal], t2)
+    {-# INLINE evolve' #-}
 
-    discounter (Black _ _ _ dsc) t = do
-        size <- getTrials
-        return $ U.replicate size $ disc dsc t
+    discount (Black _ _ _ dsc) t = disc dsc t
+    {-# INLINE discount #-}
 
     forwardGen (Black _ _ fg _) t2 = do
-        size <- getTrials
-        t1 <- gets snd
-        return $ U.replicate size $ forward fg t1 t2
+      (_, t1) <- get
+      return $ forward fg t1 t2
+    {-# INLINE forwardGen #-}
 
     maxStep _ = 100

@@ -1,45 +1,42 @@-{-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE ExistentialQuantification #-}
 
 module Quant.Models.Dupire (
     Dupire (..)
 ) where
 
+import Quant.Time
 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
+ , dupireFunc        ::  Time -> 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)
+    initialize (Dupire s _ _ _) = put (Observables [s], Time 0)
+    {-# INLINE initialize #-}
 
     evolve' d@(Dupire _ f _ _) t2 anti = do
-        (Observables (stateVec:_), t1) <- get
+        (Observables (stateVal:_), 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)
+        let vol   = f t1 stateVal
+            grwth = (fwd - vol * vol / 2) * timeDiff t1 t2
+        normResid <- lift stdNormal
+        let s' | anti      = stateVal * exp (grwth - normResid*vol)
+               | otherwise = stateVal * exp (grwth - normResid*vol)
+        put (Observables [s'], t2)
+    {-# INLINE evolve' #-}
 
-    discounter (Dupire _ _ _ dsc) t = do
-        size <- getTrials
-        return $ U.replicate size $ disc dsc t
+    discount (Dupire _ _ _ dsc) t = disc dsc t
+    {-# INLINE discount #-}
 
     forwardGen (Dupire _ _ fg _) t2 = do
         t1 <- gets snd
-        size <- getTrials
-        return $ U.replicate size $ forward fg t1 t2+        return $ forward fg t1 t2
+    {-# INLINE forwardGen #-}

@@ -1,18 +1,16 @@ {-# LANGUAGE ExistentialQuantification #-}
-{-# LANGUAGE FlexibleInstances #-}
 
 
 module Quant.Models.Heston (
     Heston (..)
 ) where
 
+import Quant.Time
 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 {
@@ -26,34 +24,30 @@   , 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)
+    initialize (Heston s v0 _ _ _ _ _ _) = put (Observables [s, v0], Time 0)
+    {-# INLINE initialize #-}
 
     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)
+        let grwth = (fwd - vState/2) * t
+            t = timeDiff t1 t2
+        resid1  <- lift stdNormal
+        resid2' <- lift stdNormal
+        let 
+          op = if anti then (-) else (+)
+          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)
+    {-# INLINE evolve' #-}
 
-    discounter (Heston _ _ _ _ _ _ _ d) t = do
-        size <- getTrials
-        return $ U.replicate size $ disc d t
+    discount (Heston _ _ _ _ _ _ _ d) t = disc d t
+    {-# INLINE discount #-}
 
     forwardGen (Heston _ _ _ _ _ _ fg _) t2 = do
-        size <- getTrials
         t1 <- gets snd
-        return $ U.replicate size $ forward fg t1 t2
+        return $ forward fg t1 t2
+    {-# INLINE forwardGen #-}
 
-    maxStep _ = 1/250+    maxStep _ = 1/12

@@ -1,17 +1,16 @@-{-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE ExistentialQuantification #-}
 
 module Quant.Models.Merton (
     Merton (..)
 ) where
 
+import Quant.Time
 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 {
@@ -32,30 +31,29 @@             --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)
+    initialize (Merton s _ _ _ _ _ _) = put (Observables [s], Time 0)
+    {-# INLINE initialize #-}
 
     evolve' m@(Merton _ vol intensity mu sig _ _) t2 anti = do
-        (Observables (stateVec:_), t1) <- get
+        (Observables (stateVal:_), 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)
+            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
+        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)
+    {-# INLINE evolve' #-}
 
-    discounter (Merton _ _ _ _ _ _ dsc) t = do
-        size <- getTrials
-        return $ U.replicate size $ disc dsc t
+    discount (Merton _ _ _ _ _ _ dsc) t = disc dsc t
+    {-# INLINE discount #-}
 
     forwardGen (Merton _ _ _ _ _ fg _) t2 = do
-        size <- getTrials
         t1 <- gets snd
-        return $ U.replicate size $ forward fg t1 t2+        return $ forward fg t1 t2
+    {-# INLINE forwardGen #-}

@@ -0,0 +1,20 @@+
+module Quant.Models.Processes (
+    ProcessSpec (..)
+  , lognormal
+) where
+
+data ProcessSpec = ProcessSpec {
+    procInit    :: Double
+  , procGrowth  :: Double
+  , procElapsed :: Double
+}
+
+lognormal :: ProcessSpec -> Double -> Double -> Double 
+lognormal (ProcessSpec initVal r t) vol normRand = initVal * exp ( g + sig * normRand )
+  where
+    g   = (r - vol*vol/2) * t
+    sig = vol * sqrt t
+
+--cir :: ProcessSpec -> Double 
+--cir (ProcessSpec initVal r t) vol normRand 

@@ -11,8 +11,6 @@   , Discretize(..)
   , OptionType(..)
 
-  , getTrials
-
   )
 where
 
@@ -21,10 +19,11 @@ import Control.Applicative
 import Control.Monad.State
 import Data.Functor.Identity
+import Quant.Time
 import Data.RVar
+import Data.Foldable (foldl')
 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)
@@ -50,87 +49,109 @@ 
     -- | Initializes a Monte Carlo simulation for a given number of runs.
     initialize :: Discretize a => a   -- ^ Model
-                               -> Int -- ^ number of trials
-                               -> MonteCarlo (Observables, Double) ()
+                               -> MonteCarlo (MCObservables, Time) ()
 
     -- | 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 :: Discretize a => a            -- ^ Model
+                           -> Time         -- ^ time to evolve to
+                           -> Bool         -- whether or not to use flipped variates
+                           -> MonteCarlo (MCObservables, Time) ()
     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
+        unless (t2==t1) $
+          if timeDiff t1 t2 < ms then 
+              evolve' mdl t2 anti
+          else do
+              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.
-    discounter :: Discretize a => a -> Double -> MonteCarlo (Observables, Double) (U.Vector Double)
+    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
+
     -- | Stateful forward generator for a given model at a certain time.
-    forwardGen :: Discretize a => a -> Double -> MonteCarlo (Observables, Double) (U.Vector Double)
+    forwardGen :: Discretize a => a -> Time -> MonteCarlo (MCObservables, Time) 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
+    evolve' :: Discretize a => a          -- ^ model
+                            -> Time       -- ^ time to evolve to
+                            -> Bool       -- ^ whether or not to use flipped variates
+                            -> MonteCarlo (MCObservables, Time) () -- ^ computation result
 
     -- | Determines the maximum size time-step for discretization purposes. Defaults to 1/250.
     maxStep :: Discretize a => a -> Double
     maxStep _ = 1/250
+    {-# INLINE maxStep #-}
 
     -- | Perform a simulation of a compiled basket of contingent claims.
     simulateState :: Discretize a => 
-            a                      -- ^ model
-        ->  ContingentClaimBasket  -- ^ compilied basket of claims
+           a                       -- ^ model
+        -> ContingentClaim         -- ^ 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')
+        -> MonteCarlo (MCObservables, Time) Double -- ^ computation result
+    simulateState modl (ContingentClaim ccb) trials anti = avg <$> replicateM trials singleTrial
+          where 
+            singleTrial = initialize modl >> 
+                            process (0 :: Double) Map.empty ccb []
 
-                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
+            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
+                  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)
 
-                v1 |+| v2 = U.zipWith (+) v1 v2
-                v1 |*| v2 = U.zipWith (*) v1 v2
+            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 [])                       
 
-                avg v = U.sum v / fromIntegral (U.length v)
+            process discCFs obsMap [] (cf:cfs) = do
+              evolve modl (cfTime cf) anti
+              d <- discountState modl $ cfTime cf
+              process (discCFs+d*cfAmount cf) obsMap [] cfs
 
+            process discCFs _ _ _ = return $! discCFs
+
+            insertCF (CashFlow t amt) (CashFlow t' amt':cfs)
+              | t > t' = CashFlow t' amt' : insertCF (CashFlow t amt) cfs
+              | otherwise = CashFlow t amt : CashFlow t' amt' : cfs
+            insertCF cf [] = [cf]
+
+            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 [], 0)
+                                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 (ccBasket ccs) trials anti
+            run = simulateState modl ccs trials anti
 
     -- | Like 'runSimulation', but splits the trials in two and does antithetic variates.
     runSimulationAnti :: (Discretize a,
@@ -145,16 +166,4 @@ 
     -- | '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
-
-
+    quickSimAnti mdl opts trials = runSimulationAnti mdl opts (pureMT 500) trials

@@ -1,63 +0,0 @@-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

@@ -0,0 +1,21 @@+
+module Quant.Time (
+	Time(..)
+  , timeDiff
+  , timeOffset
+  , timeFromZero
+  ) where
+
+data Time = Time Double deriving (Eq,Show,Ord)
+
+timeDiff :: Time -> Time -> Double
+timeDiff (Time x) (Time y) = y - x
+{-# INLINE timeDiff #-}
+
+timeOffset :: Time -> Double -> Time
+timeOffset (Time x) y = Time (x+y)
+{-# INLINE timeOffset #-}
+
+timeFromZero :: Time -> Double
+timeFromZero (Time x) = x
+{-# INLINE timeFromZero #-}

@@ -0,0 +1,24 @@+
+module Quant.Types (
+    CashFlow(..)
+  , Observables(..)
+  , MCObservables
+  , OptionType(..)
+  ) where
+
+import Quant.Time
+
+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)

@@ -1,10 +1,13 @@ module Quant.VolSurf (
     VolSurf (..)
  ,  FlatSurf (..)
+ ,  GridSurf (..)
 ) where
 
-
+import Quant.Math.Interpolation
 import Quant.YieldCurve
+import Quant.Time
+import qualified Data.Map as M
 
 {- | The 'VolSurf' class defines the
 basic operations of a volatility surface.
@@ -13,20 +16,22 @@ -}
 class VolSurf a where
     -- | Calculate the implied vol for a given strike/maturity.
-    vol :: VolSurf a => a -> Double -> Double -> Double
+    vol :: VolSurf a => a -> Double -> Time -> 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
+    var :: VolSurf a => a -> Double -> Time -> Double
+    var vs s t = v*v*t'
+        where 
+            v  = vol vs s t
+            t' = timeFromZero t
 
     -- | Calculates Dupire local vol for a given strike/maturity/forward generating yield curve.
-    localVol :: (VolSurf a, YieldCurve b) => a       -- ^ Volatility surface
-                                          -> Double  -- ^ Initial stock price
-                                          -> b       -- ^ 'YieldCurve' to generate forwards
-                                          -> Double  -- ^ Current stock level
-                                          -> Double  -- ^ Time
-                                          -> Double  -- ^ Local volatility
+    localVol :: (VolSurf a, YieldCurve b) => a           -- ^ Volatility surface
+                                          -> Double      -- ^ Initial stock price
+                                          -> b           -- ^ 'YieldCurve' to generate forwards
+                                          -> Double      -- ^ Current stock level
+                                          -> Time        -- ^ Time
+                                          -> Double      -- ^ Local volatility
     localVol v s0 rcurve k t 
         | w==0.0 || solution<0.0 = sqrt dwdt
         | otherwise = sqrt solution
@@ -35,23 +40,39 @@                 f = s0/dr
                 y = log $ k/f
                 dy = 1.0E-6
-                [kp, km] = [k*exp dy, k/exp dy]
+                kp = k*exp dy
+                km = k/exp dy
                 [w, wp, wm] = map (\x->var v (x/s0) t) [k, kp, km]
                 dwdy = (wp-wm)/2.0/dy
                 d2wdy2 = (wp-2.0*w+wm)/dy/dy
-                dt = min 0.0001 (t/2.0)
+                dt = min 0.0001 (timeFromZero t/2.0)
                 dwdt = let 
                     strikept = k*dr/drpt
                     strikemt = k*dr/drmt
-                    drpt = disc rcurve $ t+dt
-                    drmt = disc rcurve $ t-dt
-                        in case t of
-                            0 -> (var v (strikept/s0) (t+dt) -w)/dt
-                            _ -> (var v (strikept/s0) (t+dt)-var v (strikemt/s0) (t-dt))/2.0/dt                       
+                    drpt = disc rcurve $ timeOffset t dt
+                    drmt = disc rcurve $ timeOffset t (-dt)
+                        in case timeFromZero t of
+                            0 -> (var v (strikept/s0) (timeOffset t dt) -w)/dt
+                            _ -> (var v (strikept/s0) (timeOffset t dt)-var v (strikemt/s0) (timeOffset t (-dt)))/2.0/dt                       
                 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
 
 instance VolSurf FlatSurf where
-    vol (FlatSurf x) _ _ = x+    vol (FlatSurf x) _ _ = x
+
+data GridSurf = GridSurf {
+    gridStrikes             :: [Double]
+  , gridMaturities          :: [Time]
+  , gridQuotes              :: M.Map (Double, Time) Double
+  , gridStrikeInterpolator  :: Interpolator1d
+  , gridTimeInterpolator    :: Interpolator1d
+}
+
+instance VolSurf GridSurf where
+    vol (GridSurf sts mats quotes vInterp tInterp) strike t = tInterp mats' interpolatedVs $ timeFromZero t
+        where
+            mats' = map timeFromZero mats
+            interpolatedVs = map ((\k -> vInterp sts k strike) . allForT) mats
+            allForT t' = map (\ x -> quotes M.! (x, t')) sts

@@ -4,6 +4,7 @@  ,  NetYC (..)
 ) where
 
+import Quant.Time
 
 {- | The 'YieldCurve' class defines the
 basic operations of a yield curve.
@@ -12,22 +13,24 @@ -}
 class YieldCurve a where
     -- | Calculate the discount factor for a given maturity.
-    disc :: YieldCurve a => a -> Double -> Double
+    disc :: YieldCurve a => a -> Time -> 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
+    forward :: YieldCurve a => a -> Time -> Time -> Double
+    forward yc t1 t2 = (/(timeFromZero t2-timeFromZero t1)) $ log $ disc yc t1 / disc yc t2
+    {-# INLINE forward #-}
 
     -- | Calculate the spot rate for a given maturity.
-    spot :: YieldCurve a => a -> Double -> Double
-    spot yc t = forward yc 0 t
+    spot :: YieldCurve a => a -> Time -> Double
+    spot yc t = forward yc (Time 0) t
+    {-# INLINE spot #-}
 
 -- |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)
+    disc (FlatCurve r) t = exp (-r*timeFromZero t)
 
 -- | 'YieldCurve' that represents the difference between two 'YieldCurve's.
 data NetYC a = NetYC a a