Skip to content

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)
}

sma1.svg

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
        }
    )
}

ema0.svg

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:

both.svg

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

Trackbacks

No Trackbacks

Comments

Display comments as Linear | Threaded

No comments

The author does not allow comments to this entry

Add Comment

E-Mail addresses will not be displayed and will only be used for E-Mail notifications.
To leave a comment you must approve it via e-mail, which will be sent to your address after submission.
Enclosing asterisks marks text as bold (*word*), underscore are made via _word_.
Form options