Vectorized implementation of exponentially weighted moving standard deviation using R?

Question

I am trying to implement a vectorized exponentially weighted moving standard deviation using R. Is this the correct approach?

ewma <- function (x, alpha) {
  c(stats::filter(x * alpha, 1 - alpha, "recursive", init = x[1]))
}
ewmsd <- function(x, alpha) {
  sqerror <- na.omit((x - lag(ewma(x, alpha)))^2)
  ewmvar <- c(stats::filter(sqerror * alpha, 1 - alpha, "recursive", init = 0))
  c(NA, sqrt(ewmvar))
}

I'm guessing it's not, since its output is different from Python's pandas.Series.ewm.std() function.

When I run

ewmsd(x = 0:9, alpha = 0.96)

the output is

 [1]        NA 0.2236068 0.4874679 0.7953500 1.1353903 1.4993855 1.8812961 2.2764708 2.6812160 3.0925367

However, with

pd.Series(range(10)).ewm(alpha = 0.96).std()

the output is

0         NaN
1    0.707107
2    0.746729
3    0.750825
4    0.751135
5    0.751155
6    0.751156
7    0.751157
8    0.751157
9    0.751157

Answer 1

According to the documentation for Pandas, the pandas.Series.ewm() function receives an adjust parameter, which defaults to TRUE. When adjust == TRUE, the exponentially weighted moving average from pandas.Series.ewm.mean() is calculated through weights, not recursively. Naturally, this affects the standard deviation output as well. See this Github issue and this question for more info.

Here's a vectorized solution in R:

   ewmsd <- function(x, alpha) {
      n <- length(x)
      sapply(
        1:n,
        function(i, x, alpha) {
          y <- x[1:i]
          m <- length(y)
          weights <- (1 - alpha)^((m - 1):0)
          ewma <- sum(weights * y) / sum(weights)
          bias <- sum(weights)^2 / (sum(weights)^2 - sum(weights^2))
          ewmsd <- sqrt(bias * sum(weights * (y - ewma)^2) / sum(weights))
        },
        x = x,
        alpha = alpha
      )
    }