Fast R implementation of an Exponentially Weighted Moving Average?

Question

I'd like to perform an exponentially weighted moving average (with parameterization defined here) on a vector in R. Is there a better implementation than my first attempt below?

My first attempt was:

ewma <- function(x, a) {
  n <- length(x)
  s <- rep(NA,n)
  s[1] <- x[1]
  if (n > 1) {
    for (i in 2:n) {
      s[i] <- a * x[i] + (1 - a) * s[i-1]
    }
  }
  return(s)
}

y <- 1:1e7
system.time(s <- ewma(y,0.5))
#user  system elapsed 
#   2.48    0.00    2.50

In my second attempt, I thought I could do better by vectorizing:

ewma_vectorized <- function(x,a) {
  a <- 0.1
  n <- length(x)
  w <- cumprod(c(1, rep(1-a, n-1)))
  x1_contribution <- w * x[1]
  w <- a * w
  x <- x[-1]
  s <- apply(as.array(1:(n-1)), 1, function(i,x,w){sum(w[i:1] * x[1:i])}, x=x, w=w)
  s <- x1_contribution + c(0,s)
  return(s)
}

system.time(s <- ewma_vectorized(y,0.5))
# I stopped the program after it continued to run for 4min

I guess I shouldn't have been too surprised by the results in my second attempt. It was a pretty ugly attempt at vectorization. But there has to be something like this that improves on my first attempt...right?

UPDATE:

I did find a better implementation here and adapted it as follows:

ewma_vectorized_v2 <- function(x, a) {
  s1 <- x[1]
  sk <- s1
  s <- vapply(x[-1], function(x) sk <<- (1 - a) * x + a * sk, 0)
  s <- c(s1, s)
  return(s)
}

system.time(s <- ewma_vectorized_v2(y,0.5))
# user  system elapsed 
#   1.74    0.01    1.76

Are you just looking to use a fast implementation? [The TTR package implements it in C](https://github.com/joshuaulrich/TTR/blob/master/src/moving_averages.c). — Gregor Thomas, Mar 13 '17 at 21:52
I was trying to do this in R and avoid using any packages. But that is a helpful a link. — lowndrul, Mar 13 '17 at 22:00
@Gregor: that's a *terrible* package. :P The OP could use `stats::filter` instead. — Joshua Ulrich, Mar 13 '17 at 22:19
@JoshuaUlrich Ha! I'll certainly defer to you on that one. I'll admit I haven't used it, just followed [Gabor's link here](http://stackoverflow.com/q/39260128/903061). — Gregor Thomas, Mar 13 '17 at 22:29

score 10 · Accepted Answer · answered Mar 13 '17 at 22:28

You can do this with stats::filter:

ewma.filter <- function (x, ratio) {
  c(filter(x * ratio, 1 - ratio, "recursive", init = x[1]))
}
set.seed(21)
x <- rnorm(1e4)
all.equal(ewma.filter(x, 0.9), ewma(x, 0.9))
# [1] TRUE

This is a bit faster than the compiled version of your first attempt:

ewma <- compiler::cmpfun(function(x, a) {
  n <- length(x)
  s <- rep(NA,n)
  s[1] <- x[1]
  if (n > 1) {
    for (i in 2:n) {
      s[i] <- a * x[i] + (1 - a) * s[i-1]
    }
  }
  return(s)
})
microbenchmark(ewma.filter(x,0.9), ewma(x, 0.9))
Unit: microseconds
                expr      min        lq   median       uq      max neval
 ewma.filter(x, 0.9)  318.508  341.7395  368.737  473.254 1477.000   100
        ewma(x, 0.9) 1364.997 1403.4015 1458.961 1503.876 2221.252   100

digEmAll · Answer 2 · 2017-03-13T22:24:23.173

On my machine (R 3.3.2 windows) you first loop takes ~16 seconds. Enabling the jit compilation, by adding the line compiler::enableJIT(2) before the function definition, the code runs in ~1 second.

Anyway, if you really want to be fast, I think you should use C/C++, as you can see in the following example using Rcpp :

library(Rcpp)

sourceCpp(
  code = 
    "
     #include <Rcpp.h>
     // [[Rcpp::export]]
     Rcpp::NumericVector ewmaRcpp(Rcpp::NumericVector x, double a){
       int n = x.length();
       Rcpp::NumericVector s(n);
       s[0] = x[0];
       if (n > 1) {
         for (int i = 1; i < n; i++) {
           s[i] = a * x[i] + (1 - a) * s[i-1];
         }
       }
       return s;
     }

    ")

y <- 1:1e7
system.time(s2 <- ewmaRcpp(y,0.5))
# user  system elapsed 
# 0.06    0.01    0.07

Dirk Eddelbuettel · Answer 3 · 2017-03-13T22:38:06.817

@digEmAll was very kind with a Rcpp version, but also note that you could just use the TTR package, or, as its author notes, the stats::filter() approach I used in a post on the (now defunct) R Graph Gallery a decade ago.

Anyway, a quick shootout shootout shows the Rcpp version as much faster ... which probably means we got the parameterisation wrong:

R> sourceCpp("/tmp/ema.cpp")

R> library(TTR)

R> library(microbenchmark)

R> y <- as.numeric(1:1e6)   # else the sequence creates ints

R> microbenchmark(ewmaRcpp(y,0.5), EMA(y, n=10))
Unit: milliseconds
             expr      min       lq     mean   median       uq      max neval cld
 ewmaRcpp(y, 0.5)  2.43666  2.45705  3.06699  2.47283  2.51439  9.76883   100  a 
   EMA(y, n = 10) 15.13208 15.37910 21.36930 15.59278 17.26318 76.45934   100   b
R>

Actually, lambda=0.5 is an exceptionally strong decay which would correspond to a half-life on one day, or N=1. If I use that, the gap is even wider.

For completeness, the whole file which can just be Rcpp::sourceCpp()-ed:

#include <Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector ewmaRcpp(Rcpp::NumericVector x, double a){
  int n = x.length();
  Rcpp::NumericVector s(n);
  s[0] = x[0];
  if (n > 1) {
    for (int i = 1; i < n; i++) {
      s[i] = a * x[i] + (1 - a) * s[i-1];
    }
  }
  return s;
}

/*** R
library(TTR)
library(microbenchmark)

y <- as.numeric(1:1e6)   # else the sequence creates ints
microbenchmark(ewmaRcpp(y,0.5), EMA(y, n=1))
*/

I would set `TTR::EMA(y, ratio=0.5)` for the benchmark (not that I think it will make much difference). Also note that `TTR::EMA` is doing a bit more than `stats::filter` and this Rcpp version: it has some error-checking, handles leading `NA`, and uses the `try.xts` and `reclass` paradigm to internally handle many different types of objects. — Joshua Ulrich, Mar 13 '17 at 22:58

Fast R implementation of an Exponentially Weighted Moving Average?

3 Answers3

Linked