Initialised repo

.gitignore

@@ -0,0 +1,6 @@
+*.hi
+*.o
+*.swp
+cabal-dev/*
+dist/*
+

LICENSE

@@ -0,0 +1,30 @@
+Copyright (c) 2012, Ben Edwards
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Ben Edwards nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Setup.hs

@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain

hawkes.cabal

@@ -0,0 +1,22 @@
+-- Initial hawkes.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name:                hawkes
+version:             0.1.0.0
+-- synopsis:            
+-- description:         
+license:             BSD3
+license-file:        LICENSE
+author:              Ben Edwards
+maintainer:          [email protected]
+-- copyright:           
+category:            Finance
+build-type:          Simple
+cabal-version:       >=1.8
+
+executable hawkes-generate
+  main-is:             HawkesGen.hs
+  hs-source-dirs:      src
+  build-depends:       base ==4.5.*,
+                       vector,
+                       mersenne-random

src/HawkesGen.hs

@@ -0,0 +1,85 @@
+{-# LANGUAGE BangPatterns #-}
+
+module Main
+  where
+
+import qualified Data.Vector.Unboxed as U
+
+import Control.Arrow ((&&&))
+import Control.Monad.Loops (unfoldrM)
+import Data.List (foldl')
+import System.Random.Mersenne 
+  ( MTGen
+  , newMTGen
+  , random
+  )
+
+main :: IO () 
+main = do gen <- newMTGen Nothing
+          evs <- generateHawkes gen hc 1000
+          print evs
+  where hc = HC 1.2 0.6 0.8 
+
+-- | This is the context that a user can pass around in order to
+--   generate hawkes processes. It is for simple univariate processes
+data HawkesContext = HC { hcLambda :: {-# UNPACK #-} !Double 
+                        , hcAlpha  :: {-# UNPACK #-} !Double 
+                        , hcBeta   :: {-# UNPACK #-} !Double 
+                        } deriving (Show, Eq)
+
+-- | This is an internal structure for the unfold
+data HawkesSeed = HS { hsGen     :: MTGen
+                     , hsHc      :: HawkesContext
+                     , hsHorizon :: Double 
+                     , hsEvents  :: [Double]
+                     , hsLStar   :: Double
+                     } 
+
+untilM :: Monad m => (a ->  m (Maybe a)) -> a -> m a
+untilM f x = f x >>= maybe (return x) (\x -> untilM f x)
+
+generateHawkes :: MTGen         -- ^ Generator for even timings 
+               -> HawkesContext -- ^ Context for the process
+               -> Double        -- ^ The event time horizon
+               -> IO [Double]   -- ^ GeneratedProcess
+generateHawkes g hc t = 
+  do let lstar = hcLambda hc
+     ev <- (\x -> -(1 / lstar * log x)) `fmap` random g
+     let hs = HS g hc t [ev] lstar
+     if ev <= t then hsEvents `fmap` untilM generateHawkes' hs
+                else return []
+
+
+generateHawkes' :: HawkesSeed -> IO (Maybe HawkesSeed)
+generateHawkes' hs@(HS g hc t evs lstar) =
+  do ev' <- eventloop g hc t ev evs lstar
+     return $! update hs `fmap` ev'
+  where ev = head evs
+        lstar = intensity hc ev (tail evs) + (hcAlpha hc)
+        update hs x = hs { hsEvents = x:evs }
+
+
+eventloop :: MTGen         -- ^ PRNG for event times
+          -> HawkesContext -- ^ Context
+          -> Double        -- ^ Event horizon
+          -> Double        -- ^ Last generated event
+          -> [Double]      -- ^ Event times
+          -> Double        -- ^ Max intensity
+          -> IO (Maybe Double)  -- ^ Newly generated event
+eventloop g hc t ev evs lstar = do
+  ev' <- (\x -> ev - (1 / lstar) * log x) `fmap` random g 
+  if ev' >= t 
+    then return Nothing
+    else do d <- random g
+            let i = intensity hc ev' evs
+            if d <= i / lstar
+              then return $ Just ev'
+              else eventloop g hc t ev' evs i
+
+intensity :: HawkesContext -- ^ Process Context
+          -> Double        -- ^ Event time        
+          -> [Double]      -- ^ Events
+          -> Double        -- ^ Intensity
+intensity (HC l a b) ev evs = foldl' binop l evs
+  where binop z x = z + a * (exp $ -b * (ev - x))
+