{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE QuantifiedConstraints #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TupleSections         #-}
module Q.Stochastic.Process
        where
import           Control.Monad
import           Control.Monad.State
import           Data.List             (foldl')
import           Data.RVar
import           Data.Random
import           Numeric.LinearAlgebra

rwalkState :: RVarT (State Double) Double
rwalkState :: RVarT (State Double) Double
rwalkState = do
    Double
prev <- State Double Double -> RVarT (State Double) Double
forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
(MonadTrans t, Monad m) =>
m a -> t m a
lift State Double Double
forall s (m :: * -> *). MonadState s m => m s
get
    Double
change  <- Normal Double -> RVarT (State Double) Double
forall (d :: * -> *) t (n :: * -> *).
Distribution d t =>
d t -> RVarT n t
rvarT Normal Double
forall a. Normal a
StdNormal

    let new :: Double
new = Double
prev Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
change
    State Double () -> RVarT (State Double) ()
forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
(MonadTrans t, Monad m) =>
m a -> t m a
lift (Double -> State Double ()
forall s (m :: * -> *). MonadState s m => s -> m ()
put Double
new)
    Double -> RVarT (State Double) Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
new

type Time = Double

-- Dont know why this wasn't done.
-- Is there an easier way to do this where we either lift or return?
instance (Num a) => Num (RVarT m a) where
  + :: RVarT m a -> RVarT m a -> RVarT m a
(+) = (a -> a -> a) -> RVarT m a -> RVarT m a -> RVarT m a
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> a
forall a. Num a => a -> a -> a
(+)
  (-) = (a -> a -> a) -> RVarT m a -> RVarT m a -> RVarT m a
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (-)
  * :: RVarT m a -> RVarT m a -> RVarT m a
(*) = (a -> a -> a) -> RVarT m a -> RVarT m a -> RVarT m a
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> a
forall a. Num a => a -> a -> a
(*)
  abs :: RVarT m a -> RVarT m a
abs = (a -> a) -> RVarT m a -> RVarT m a
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM a -> a
forall a. Num a => a -> a
abs
  signum :: RVarT m a -> RVarT m a
signum = (a -> a) -> RVarT m a -> RVarT m a
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM a -> a
forall a. Num a => a -> a
signum
  fromInteger :: Integer -> RVarT m a
fromInteger Integer
x = a -> RVarT m a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> RVarT m a) -> a -> RVarT m a
forall a b. (a -> b) -> a -> b
$ Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
x



-- |Discretization of stochastic process over given interval
class (Num b) => Discretize d b where
  -- |Discretization of the drift process.
  dDrift  :: (StochasticProcess a b) => a -> d -> (Time, b) -> RVar b
  -- |Discretization of the diffusion process.
  dDiff   :: (StochasticProcess a b) => a -> d -> (Time, b) -> RVar b
  -- |dt used.
  dDt     :: (StochasticProcess a b) => a -> d -> (Time, b) -> Time


-- |A stochastic process of the form \(dX_t = \mu(X_t, t)dt + \sigma(S_t, t)dB_t \)
class (Num b) => StochasticProcess a b where
  -- |The process drift.
  pDrift  :: a -> (Time, b) -> RVar b
  -- |The process diffusion.
  pDiff   :: a -> (Time, b) -> RVar b

  -- |Evolve a process from a given state to a given time.
  pEvolve :: (Discretize d b) => a         -- ^The process
                             -> d         -- ^Discretization scheme
                             -> (Time, b) -- ^Initial state
                             -> Time      -- ^Target time t.
                             -> RVar b    -- ^\(dB_i\).
                             -> RVar b    -- ^\(X(t)\).
  pEvolve a
p d
disc s0 :: (Double, b)
s0@(Double
t0, b
x0) Double
t RVar b
dw = do
    if Double
t0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
t then b -> RVar b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x0 else do
      s' :: (Double, b)
s'@(Double
t', b
b') <- a -> d -> (Double, b) -> RVar b -> RVar (Double, b)
forall a b d.
(StochasticProcess a b, Discretize d b, Num b) =>
a -> d -> (Double, b) -> RVar b -> RVar (Double, b)
pEvolve' a
p d
disc (Double, b)
s0 RVar b
dw
      if Double
t' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
t then b -> RVar b
forall (m :: * -> *) a. Monad m => a -> m a
return b
b' else a -> d -> (Double, b) -> Double -> RVar b -> RVar b
forall a b d.
(StochasticProcess a b, Discretize d b) =>
a -> d -> (Double, b) -> Double -> RVar b -> RVar b
pEvolve a
p d
disc (Double, b)
s' Double
t RVar b
dw

  -- |Similar to evolve, but evolves the process with the discretization scheme \(dt\).
  pEvolve' :: (Discretize d b, Num b) => a -> d -> (Time, b) -> RVar b -> RVar (Time, b)
  pEvolve' a
process d
discr s :: (Double, b)
s@(Double
t, b
b) RVar b
dw = do
    let !newT :: Double
newT = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ a -> d -> (Double, b) -> Double
forall d b a.
(Discretize d b, StochasticProcess a b) =>
a -> d -> (Double, b) -> Double
dDt a
process d
discr (Double, b)
s
        !newX :: RVar b
newX = do
               b
drift <- a -> d -> (Double, b) -> RVar b
forall d b a.
(Discretize d b, StochasticProcess a b) =>
a -> d -> (Double, b) -> RVar b
dDrift a
process d
discr (Double, b)
s
               b
diff  <- a -> d -> (Double, b) -> RVar b
forall d b a.
(Discretize d b, StochasticProcess a b) =>
a -> d -> (Double, b) -> RVar b
dDiff a
process d
discr (Double, b)
s
               b
dw' <- RVar b
dw
               b -> RVar b
forall (m :: * -> *) a. Monad m => a -> m a
return (b -> RVar b) -> b -> RVar b
forall a b. (a -> b) -> a -> b
$ b
b b -> b -> b
forall a. Num a => a -> a -> a
+ b
drift b -> b -> b
forall a. Num a => a -> a -> a
+ b
diff b -> b -> b
forall a. Num a => a -> a -> a
* b
dw'
        newX :: RVar b

    (Double
newT,) (b -> (Double, b)) -> RVar b -> RVar (Double, b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>  RVar b
newX

-- |Geometric Brownian motion
data GeometricBrownian = GeometricBrownian {
    GeometricBrownian -> Double
gbDrift :: Double -- ^Drift
  , GeometricBrownian -> Double
gbDiff  :: Double -- ^Vol
} deriving (Int -> GeometricBrownian -> ShowS
[GeometricBrownian] -> ShowS
GeometricBrownian -> String
(Int -> GeometricBrownian -> ShowS)
-> (GeometricBrownian -> String)
-> ([GeometricBrownian] -> ShowS)
-> Show GeometricBrownian
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [GeometricBrownian] -> ShowS
$cshowList :: [GeometricBrownian] -> ShowS
show :: GeometricBrownian -> String
$cshow :: GeometricBrownian -> String
showsPrec :: Int -> GeometricBrownian -> ShowS
$cshowsPrec :: Int -> GeometricBrownian -> ShowS
Show)


instance StochasticProcess GeometricBrownian Double where
--  pDrift :: GeometricBrownian -> (Time, Double) -> RVar Double
  pDrift :: GeometricBrownian -> (Double, Double) -> RVar Double
pDrift GeometricBrownian
p (Double
_, Double
x) = Double -> RVar Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> RVar Double) -> Double -> RVar Double
forall a b. (a -> b) -> a -> b
$ GeometricBrownian -> Double
gbDrift GeometricBrownian
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x -- drift is prpotional to the spot.
  pDiff :: GeometricBrownian -> (Double, Double) -> RVar Double
pDiff  GeometricBrownian
p (Double
_, Double
x) = Double -> RVar Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> RVar Double) -> Double -> RVar Double
forall a b. (a -> b) -> a -> b
$ GeometricBrownian -> Double
gbDiff GeometricBrownian
p  Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x -- diffisuion is also prportional to the spot.


-- | Ito process
data ItoProcess = ItoProcess {
        ItoProcess -> (Double, Double) -> Double
ipDrift :: (Time, Double) -> Double,
        ItoProcess -> (Double, Double) -> Double
ipDiff  :: (Time, Double) -> Double
}