simulateVolatility :: Float -> Int -> Float -> Float -> Float -> IO [Float]
simulateVolatility sigma kbar gamma_kbar b m0 = do
  gen <- newStdGen
  let ns = take 1 $ randoms gen :: [Int]
  markovStateVector <- generateMarkovStateVector m0 kbar

  -- forM monad will keep the state of current markov state vector.
  -- Each iteration will return one point from time series simulation
  -- and update markov state vector.
  let sim = flip evalState markovStateVector $ do
              forM [0..] $ \_ -> do
                -- Get Markov state vector from state
                stateVector <- get

                -- updateMarkovStateVector already uses unsafe IO, so we might as
                -- well abuse it a bit more.
                let nextStateVector = unsafePerformIO $
                                      updateMarkovStateVector stateVector gamma_v m0 kbar
                -- Save updated state vector for next step
                put nextStateVector

                -- Calculate next step in the simulation
                let result = sigma * (sqrt  (product stateVector))
                return result

  return sim
    where gamma_v = gammaVector gamma_kbar kbar b


-- Computes gamma vector from gamma_kbar by inverting the ralationship
-- defined on page 54 of Calvet-Fisher.
simulateVolatility :: Float -> Int -> Float -> Float -> Float -> IO [Float]
simulateVolatility sigma kbar gamma_kbar b m0 = do
  gen <- newStdGen
  let ns = take 1 $ randoms gen :: [Int]
  markovStateVector <- generateMarkovStateVector m0 kbar

  -- forM monad will keep the state of current markov state vector.
  -- Each iteration will return one point from time series simulation
  -- and update markov state vector.
  let sim = flip evalState markovStateVector $ do
              forM [0..] $ \_ -> do
                -- Get Markov state vector from state
                stateVector <- get

                -- updateMarkovStateVector already uses unsafe IO, so we might as
                -- well abuse it a bit more.
                let nextStateVector = unsafePerformIO $
                                      updateMarkovStateVector stateVector gamma_v m0 kbar
                -- Save updated state vector for next step
                put nextStateVector

                -- Calculate next step in the simulation
                let result = sigma * (sqrt  (product stateVector))
                return result

  return sim
    where gamma_v = gammaVector gamma_kbar kbar b


-- Computes gamma vector from gamma_kbar by inverting the ralationship
-- defined on page 54 of Calvet-Fisher.
simulateVolatility sigma kbar gamma_kbar b m0 = do
  gen <- newStdGen
  let ns = take 1 $ randoms gen :: [Int]
  markovStateVector <- generateMarkovStateVector m0 kbar

  -- forM monad will keep the state of current markov state vector.
  -- Each iteration will return one point from time series simulation
  -- and update markov state vector.
  let sim = flip evalState markovStateVector $ do
              forM [0..] $ \_ -> do
                -- Get Markov state vector from state
                stateVector <- get

                -- updateMarkovStateVector already uses unsafe IO, so we might as
                -- well abuse it a bit more.
                let nextStateVector = unsafePerformIO $
                                      updateMarkovStateVector stateVector gamma_v m0 kbar
                -- Save updated state vector for next step
                put nextStateVector

                -- Calculate next step in the simulation
                let result = sigma * (sqrt  (product stateVector))
                return result

  return sim
    where gamma_v = gammaVector gamma_kbar kbar b


-- Computes gamma vector from gamma_kbar by inverting the ralationship
-- defined on page 54 of Calvet-Fisher.