## Exponential Moving Average (EMA) Rates, part 2

In the last post, we simulated some Poisson data and then verified it by looking at its histogram and some descriptive statistics. We also built a basic sliding-window implementation and graphed its output.

To continue on, we’ll need to build a more realistic implementation, along with a method to feed it the simulated events. With that in hand, we’ll build an EMA function specialized for Poisson events.

### Incremental SMA implementation

First, we’ll make a new SMA implementation. We’ll construct an SMA object with the size of window as a parameter. It will have methods to incrementally update its state, or read the current average.

make.sma1 <- function (width) { force(width) data <- NULL clipped <- function (time) { data[data > time - width] } list( update=function (time) { data <<- c(clipped(time), time) }, read=function (time) { data <<- clipped(time) length(data) / width } ) }

Together with a function to feed in the observations, we can then graph the results:

run.it <- function (ma) { all.t <- seq(0, 480, 1) i <- 1 samples <- NULL for (t in all.t) { while (arrival.times[i] <= t) { ma$update(arrival.times[i]) i <- i + 1 } samples <- c(samples, ma$read(t)) } data.frame(t=all.t, r=samples) }

### Incremental EMA implementation

Let’s start building an online (EMA) approach. The simplest arrangement is to keep track of how many observations have been recorded, along with the time of the last observation. When a new observation arrives, the time interval between it and the last observation is calculated (`this.interval`

), along with the interval from the last observation to the new beginning of the sliding window (`last.interval`

). Then, the previous count is rescaled from the full width of the window to only `last.interval`

, and the new observation accumulated.

If we’re only taking a reading, and not recording an observation, we still rescale the count in that manner, but skip the final step of accumulating.

make.ema0 <- function (width) { force(width) count <- 0 last <- 0 advance <- function (time) { rate <- count / width this.interval <- time - last last.interval <- width - this.interval count <<- rate * last.interval last <<- time } list( update=function (time) { advance(time) count <<- count + 1 }, read=function (time) { advance(time) count / width } ) }

### Inconsistency

At this point, the code has a significant flaw (other than the possibility of `last.interval`

going negative): updating `count`

is highly dependent on the frequency with which the window is advanced. Here’s an example:

```
ema <- make.ema0(2)
ema$update(1)
ema$read(2)
```

0.25

If we add an intervening read, we get a different result:

```
ema <- make.ema0(2)
ema$update(1)
ema$read(1.5)
ema$read(2)
```

0.28125

The problem is that by splitting the window into two intervals (`this.interval`

and `last.interval`

) according to time-of-observation (or time-of-read), we’re effectively compounding the discount rate at intervals that are arbitrary. It’s similar to accumulating interest in a bank account weekly for one period, then monthly for another.

The solution is to compound continuously to find the true discount for the intervals.

We can get a feel for what a corrected answer would look like by inserting several small advances into the example:

ema <- make.ema0(2) ema$update(1) for (i in seq(1.01, 1.99, by=0.01)) { ema$read(i) } ema$read(2)

0.302885218245364

(Incidentally, either inserting several small advances as above, or the limit that doing so approaches, would eliminate the problem with `last.interval`

going negative.)

### Fixing the EMA

We start be rewriting `advance`

. If every temporary value in the function is substitued, it can be written as:

advance <- function (time) { count <<- count * (1 + (last - time) / width) last <<- time }

The expression for `count`

now follows the form of compounding for one period. If we take the limit as we call `advance`

over more and smaller intervals, we get continuous compounding. The complete code looks like:

make.ema1 <- function (width) { force(width) count <- 0 last <- 0 advance <- function (time) { count <<- count * exp((last - time) / width) last <<- time } list( update=function (time) { advance(time) count <<- count + 1 }, read=function (time) { advance(time) count / width } ) }

Trying again the examples shows consistent results, and close to our final approximation:

```
ema <- make.ema1(2)
ema$update(1)
ema$read(2)
```

0.303265329856317

```
ema <- make.ema1(2)
ema$update(1)
ema$read(1.5)
ema$read(2)
```

0.303265329856317

Finally, graphing the EMA against the SMA shows that the EMA is tuned to roughly the same period as the SMA, and shows the characteristic “warm-up” of most EMA implementations:

The next post explains how to start this EMA to fix the warm-up period.

## Comments

Display comments as Linear | Threaded