Exponential moving average of a vector in R

Question

I have a simple vector as follows:

x = c(14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)

I am trying to find the rolling EMA of this vector using the following function -

library(TTR)
y = EMA(x, 5)

I am getting the result as follows -

 [1]     NA     NA     NA     NA 13.33400 13.22267 13.52844 14.44563 16.51042 16.88695

However, I want the result as below -

 [1]     14.24 14.03 13.06 13.43 13.33400 13.22267 13.52844 14.44563 16.51042 16.88695

1. The first value should be the same as in the original vector
2. The second value should be the EMA of the first and second value
3. The third value should be the EMA of the initial three values in the vector
4. The fourth value should be the EMA of the initial four values in the vector

The rest of the calculation is correctly handled by the function EMA

Solutions I tried -

1. Running the following command - zoo::rollapplyr(x, width = 5, FUN = EMA, partial = TRUE) will give error as EMA has its own rolling window.

2. Using function stats::filter works but the answer is not correct as I am not sure about the right value of the ratio parameter. Fast R implementation of an Exponentially Weighted Moving Average? Here is a custom function.

ema_2 <- function (k, width) {
  ratio <- 2/(width + 1)
  c(stats::filter(k * ratio, 1 - ratio, "convolution", init = k[1]))
}

The ideal solution should take at most twice the computation time as taken by EMA function of TTR library.

Here are the performance results of the 2 solutions shared by Waldi and Andre.

              expr     min       lq     mean   median       uq      max neval cld
    TTR::EMA(x, 5) 433.593 457.5815 500.9478 477.0535 530.7105  1128.49  1000   a
        EMA3(x, 5) 445.388 468.9585 515.2009 490.4345 546.5025  1843.46  1000   a
 rollmeanEMA(x, 5) 436.689 461.0885 535.7035 481.8815 538.3150 33122.75  1000   a

Thanks!

Answer 1

Looking at C source code of EMA shows that the first value is the mean of the averaging window:

    /* Raw mean to start EMA */
    double seed = 0.0;
    for(i = first; i < first + i_n; i++) {
      d_result[i] = NA_REAL;
      seed += d_x[i] / i_n;
    }
    d_result[first + i_n - 1] = seed;

This can easily be calculated to replace the NAs :

x = c(14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)

EMA2 <- function(x,n) {
  y = TTR::EMA(x, n)
  noNA <- which.min(is.na(x))
  y[noNA:(noNA+n-2)] <- cumsum(x[noNA:(noNA+n-2)])/1:(n-1)
  y
}

EMA2(x,5)
#>  [1] 14.24000 14.03000 13.60333 13.43250 13.33400 13.22267 13.52844 14.44563
#>  [9] 16.51042 16.88695

this also works with leading NAs:

x = c(NA, NA, 14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)
EMA2(x,5)
#> [1]       NA       NA 14.24000 14.03000 13.60333 13.43250 13.33400 13.22267 13.52844 14.44563
#> [11] 16.51042 16.88695

The calculation overhead on this short vector is minimal, and this should be even better on a longer vector :

microbenchmark::microbenchmark(TTR::EMA(x,5),EMA2(x,5),times=1000)

#> Unit: microseconds
#>           expr   min    lq     mean median     uq   max neval cld
#> TTR::EMA(x, 5) 157.7 161.8 181.6156  164.0 180.55 593.5  1000  a 
#>     EMA2(x, 5) 164.2 167.5 193.0643  170.6 193.20 857.1  1000   b

Answer 2

This gives the desired result:

require(TTR)

x <- c(14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)

rollmeanEMA <- function(vec, len) {
  c(cumsum(vec[1:(len-1)]) / seq_along(vec[1:(len-1)]),
    EMA(vec, len)[len:length(vec)])
}

rollmeanEMA(x,5)
#[1] 14.24000 14.03000 13.60333 13.43250 13.33400 13.22267 13.52844 14.44563
#[9] 16.51042 16.88695

Edit: As I suggested in the comments, replacing the NA part with mean(). This gives a tremendous speedup. Plus, removed the surrounding condition.

y <- rnorm(1000000)

system.time( rollmeanEMA(y,10000) )
#   user  system elapsed
#  0.031   0.003   0.034

system.time( EMA(y,10000) )
#   user  system elapsed
#  0.018   0.002   0.019

Added NA "handling":

rollmeanEMA <- function(vec, len) {
  v_n <- !is.na(vec)
  c( vec[is.na(vec)],
     cumsum(vec[v_n][1:(len-1)]) / seq_along(vec[v_n][1:(len-1)]),
     EMA(vec[v_n], len)[len:length(vec[v_n])])
}