David Owen’s blog

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.