Heston.hs

{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances #-}


module Quant.Models.Heston (
    Heston (..)
) where

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 {
    hestonInit       :: Double  -- ^ Initial asset level.
  , hestonV0         :: Double  -- ^ Initial variance
  , hestonVF         :: Double  -- ^ Mean-reversion variance
  , hestonLambda     :: Double  -- ^ Vol-vol
  , hestonCorrel     :: Double  -- ^ Correlation between processes
  , hestonMeanRev    :: Double  -- ^ Mean reversion speed
  , hestonForwardGen :: a       -- ^ 'YieldCurve' to generate forwards
  , 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)

    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)

    discounter (Heston _ _ _ _ _ _ _ d) t = do
        size <- getTrials
        return $ U.replicate size $ disc d t

    forwardGen (Heston _ _ _ _ _ _ fg _) t2 = do
        size <- getTrials
        t1 <- gets snd
        return $ U.replicate size $ forward fg t1 t2

    maxStep _ = 1/250