Implementing recurrence relations on State monads (in Haskell or Scala)

Question

I am working on a new implementation of the operators in http://www.thalesians.com/archive/public/academic/finance/papers/Zumbach_2000.pdf EDIT: clearer explanation here: https://www.olseninvest.com/customer/pdf/paper/001207-emaOfEma.pdf

Briefly, it's a whole bunch of cool time series operators based on the recurrence relation of the exponential moving average, where each application of the ema() operator takes the new value and the previous result of the ema. I can't seem to do latex on this stack exchange, but anyway my problem now is a software problem.

I implemented this in Scala by hiding a var deep inside the thunks that create EMA functions. This all works, but it's super tricky, because calling ema(5) and then ema(5) again will naturally lead to a different result. I'd like to try redoing all of this using State Monads, but I'm quickly getting lost in the weeds.

For example, I have the following simplified EMA State monad in Haskell:

import Control.Monad.State

type EMAState = Double
type Tau = Double

ema :: Tau -> Double -> State EMAState Double
ema tau x = state $ \y ->
  let alpha = 1 / tau
      mu = exp(-alpha)
      mu' = 1 - mu
      y' = (mu * y) + (mu' * x)
  in (y', y')

which I can readily test in GHCI:

*Main Control.Monad.State> runState (ema 5 10) 0
(1.8126924692201818,1.8126924692201818)

applying the input 10 to a 5-period EMA initialized to 0. This is all well and good, using forM I can apply multiple input values etc. Now, the next step is to implement an "iterated EMA", which is an EMA applied to itself N times.

iEMA[n](x) = EMA(iEMA[n-1](x))

Each of these intermediate EMAs will need to have their own state, aka previous result, to correctly calculate the vector of iterated EMAs. So, what I am looking for, is something which like this, (I think):

iema :: Int -> Tau -> Double -> State [EMAState] [Double]

Which is essentially a daisy chain of EMAs:

iEMA[3](x) = EMA(EMA(EMA(x,s1),s2),s3) = (x, [s1,s2,s3]) -> ([y1,y2,y3], [s1',s2',s3'])

And if all I care about is the 3rd iterated EMA ...

... -> (y3, [s1', s2', s3'])

The paper moves on from there, creating ever more complex operators built on iterated EMAs and averages of them etc, so I want to be able to functionally and purely compose these stateful operators building ever more complex states, but still quite simple input and output.

I really feel like this is what functional programming is good at, but I don't yet have the expertise to see how to put together these State monads in the correct way. Could someone please point me in the right direction with these iterated recurrence operators?

EDIT:

A couple of helpful folks have suggested repeated application of the same ema operator to the input data, but this is not sufficient. Each ema operator needs to maintain it's own previous value. Here's an example:

tau 5               
mu  0.818730753             
muprime 0.181269247             
        ema1    ema2    ema3     
    x   0       0       0       <- States_0
    1   0.1812  0.03285 0.00595 <- States_1
    5   1.0547  0.21809 0.04441 <- States_2

The x column is the raw input, ema1 uses its left for input and it's up for recurrence/state. ema2 uses its left for input (not x!) and it's up for state. It's an ema (ema (x) ). Ditto ema3 = ema (ema (ema (x) ) ). What I would like to do, which I think must be possible, is given an ema state monad, compose the ema3 state monad, or even better, the [ema] state monad with each each subsequent ema operating on the output of the previous.

Answer 1

Updated answer...

Define:

combine :: [ a -> State s a ] -> a -> State [s] a
combine fs a = state $ \ys ->
  let zs = zipWith (\f y a -> runState (f a) y) fs ys
      pairs = chain a zs
      as' = map fst pairs
      a' = last as'         -- we are only returning one result in this case
      ys' = map snd pairs
  in (a', ys')

chain :: a -> [ a -> (a,s) ] -> [ (a,s) ]
chain a [] = []
chain a (f:fs) = let (a',s) = f a
                 in (a',s) : chain a' fs

ema3 t = combine $ replicate 3 (ema t)

ghci> runState (ema3 5 1) [0,0,0]
(5.956242778945897e-3,[0.18126924692201818,3.2858539879675595e-2,5.956242778945897e-3])

ghci> runState (do ema3 5 1; ema3 5 5) [0,0,0]
(4.441089130249448e-2,[1.0547569416524334,0.21809729359983737,4.441089130249448e-2])

The combine is easily modified to return all of the results - just return as' instead of a'.

Original answer:

 combine :: (a -> State s b) -> (b -> State t c) -> (a -> State (s,t) c)
 combine f g a = state $ \(s,t) ->
   let (b,s') = runState (f a) s
       (c,t') = runState (g b) t
   in (c,(s',t'))

Then:

ema3 tau = ema tau `combine` ema tau `combine` ema tau

and em3 has type:

ema3 :: Tau -> Double -> State ((EMAState, EMAState), EMAState) Double

For instance:

ghci> runState (ema3 5 1) ((0,0),0)
(5.956242778945897e-3,((0.18126924692201818,3.2858539879675595e-2),5.956242778945897e-3))

Note that the state type of ema3 is ((Double,Double),Double) and not a 3-tuple or list.

In your example you run (ema3 5) first with input x = 1 and then with input x = 5 with initial state ((0,0),0):

ghci> runState (do ema3 5 1; ema3 5 5) ((0,0),0)
(4.441089130249448e-2,((1.0547569416524334,0.21809729359983737),4.441089130249448e-2))

and that gives you the second row in the table.

Answer 2

I may not be fully understanding your use case, but possibly you are looking for something like this:

ema' _ [] = get >>= return
ema' tau (x:xs) = do
  y <- get
  let alpha = 1 / tau
      mu = exp $ negate alpha
      mu' = 1 - mu
      y' = (mu * y) + (mu' * x)
  put y'
  ema' tau xs

It is like your original function except it accepts a list of x values, and it recursively executes for each one, updating y each time. When none are left, it returns the value of y as the answer.

It can be run like so:

*Main> evalState (ema' 5 [10]) 0
1.8126924692201818
*Main> evalState (ema' 5 [10, 10]) 0
3.2967995396436076
*Main> evalState (ema' 5 [10, 10, 10]) 0
4.511883639059737

When using the State monad, you don't need to wrap your functions in the state $ \y -> ... business. You can simply enclose your monadic code in a do block and use put and get to access the state. In this case, for each recursive execution of the function, I grab the last y with get, and then use put after doing math to update the state.

I think that in your version, you are including the State monad without actually getting anything for it (since you don't use put or get).

Also, the State monad may be overkill for this; you could accomplish the same thing using a fold over a list of x values.

Answer 3

Update based on the comments...

Three iterations of ema can be written using the monadic bind operator >>= like this:

ema3 tau x = ema tau x >>= ema tau >>= ema tau

or using the Kleisli arrow:

ema3 tau = ema tau >=> ema tau >=> ema tau

As a diagram the computation flows like this:

          y1          /---------\
           |          |         |
           v          |         v
  x  -->  EMA   -->  EMA  -->  EMA  -->  x' = y3'
          tau        tau       tau 
           |          ^         |
           |          |         v
           \----------/         y3'

(Original answer)

This is not a complete answer, but perhaps the OP comment on whether this is going in the right direction.

Here is what I understand the computation looks like:

          y1         y2        y3
           |          |         |
           v          v         v
  x  -->  EMA   -->  EMA  -->  EMA  -->  x'
          tau1       tau2      tau3
           |          |         |
           v          v         v
          y1'        y2'       y3'

The question is whether there is an elegant way to express this as a composition of EMA blocks, e.g. something like:

ema tau1 >o> ema tau2 >o> ema tau3

for some operator >o>.

Answer 4

Let's build the handy old Mealy machine

data Mealy i o where
  Mealy :: (i -> s -> (i, s)) -> s -> Mealy i o

which has all kinds of instances

instance Arrow Mealy
instance ArrowChoice Mealy
instance ArrowApply Mealy
instance Strong Mealy
instance Choice Mealy
instance Profunctor Mealy
instance Category * Mealy
instance Monad (Mealy a)
instance Functor (Mealy a)
instance Applicative (Mealy a)
instance Pointed (Mealy a)

We can use it to build recurrence relations

recur :: (a -> a -> a) -> a -> Mealy a a
recur f a0 = Mealy (\inp prior -> let post = f inp prior in (post, post)) a0

we can iterate them with our Category instance

iter :: Int -> Mealy a a -> Mealy a a
iter 0 _ = id
iter 1 m = m
iter n m = m >>> iter (n-1) m

and then, with all this machinery, we can create an infinite stream of iterated Mealy machines

data Stream a = Stream a (Stream a) deriving Functor

instance Functor Stream
instance Applicative Stream
instance Foldable Stream
instance Traversable Stream

ints :: Stream Int
ints = go 0 where go n = Stream n (go $ n + 1)

jet :: Mealy a a -> Stream (Mealy a a)
jet m = fmap (`iter` m) ints

All of these together give us, essentially, your desired structure. But it's a little difficult to interact with directly. We'll give it its own instances to help

newtype MealyJet i o = MealyJet { runMealyJet :: Stream (Mealy i o) }

instance Profunctor MealyJet
instance Applicative (MealyJet i)

instance Category MealyJet where
  id = MealyJet (pure id) -- technically this should be `jet id`, but it's equal to pure
  MealyJet f . MealyJet g = MealyJet (liftA2 (.) f g)

viewMealyJet :: MealyJet i o -> Mealy i (Stream o)
viewMealyJet (MealyJet m) = sequenceA m

And now, we can write these EMAs as needed

type Tau = Double

ema :: Tau -> Mealy Double Double
ema tau = recur $ \fresh prior -> 
  let alpha = 1 / tau
      mu    = exp (negate alpha)
      mu'   = 1 - mu
   in (mu * y) + (mu' * x)

emaJet :: Tau -> MealyJet Double Double
emaJet = MealyJet . jet . ema

emaComp :: MealyJet Double Double
emaComp = emaJet 1 >>> emaJet 2 >>> emaJet 3 >>> emaJet 4 >>> emaJet 5

fiveStack :: Mealy Double (Stream Double)
fiveStack = viewMealyJet emaComp