Computing EWMA 2

Months ago I released a post entitled Computing EWMA exploring the functional approach to implement EWMA calculations in R and also compared that with what I called loop-oriented implementation. That was quite ammazing to see that the functional approach was almost 2 times faster.

After a while I found the function filter of the package stats.

## function (x, filter, method = c("convolution", "recursive"), 
##     sides = 2L, circular = FALSE, init = NULL) 
## NULL

This function applies linear filtering to univariate or multivariate time series. You can get a moving average by setting method="convolution" or an autoregression for method="recursive". Since EWMA is an autoregressive model, the "recursive" method is appropriate. The argument x receives the time series and filter a vector of coefficients. The "recursive" method uses the init argument to specify the initial values of the time series. Providing a series x and a filter $\phi$ the response y is much like the sequence below:

$$ \begin{align} y_1 & = 0 \hspace{0.5cm} \\ y_2 & = \phi y_1 + x_1 \hspace{0.5cm} \\ y_2 & = \phi y_2 + x_2 \hspace{0.5cm} \\ \vdots & \hspace{0.5cm} \\ y_{n+1} & = \phi y_n + x_n \end{align} $$

Taking a look at EWMA's dynamics

$$ \begin{align} \hat\sigma^2_0 & = 0 \hspace{0.5cm} \\ \hat\sigma^2_1 & = \lambda\hat\sigma^2_0 + (1 - \lambda)r^2_0 \hspace{0.5cm} \\ \hat\sigma^2_2 & = \lambda\hat\sigma^2_1 + (1 - \lambda)r^2_1 \hspace{0.5cm} \\ \vdots & \hspace{0.5cm} \\ \hat\sigma^2_{t+1} & = \lambda\hat\sigma^2_t + (1 - \lambda)r^2_t \end{align} $$

we observe that the input $x_i = (1 - \lambda)r^2_i$, the filter $\phi = \lambda$ and $y_i = \sigma^2_i$ is the response. The initial value can be set to $0$. The code implementation follows below:

ewma.filter <- function(rets, lambda) {
    r2 <- (1 - lambda) * rets^2
    sqrt(filter(r2, lambda, "recursive", init = 0))
}

See that this implementation is even simpler than the others, although it is not so easy to comprehend. Running the same test done before we see that this is not only simpler but also faster. It runs almost 2 times faster than the functional approach.

system.time(replicate(10000, ewma.loop(rets, lambda)))

##    user  system elapsed 
##   4.101   0.026   4.127

system.time(replicate(10000, ewma.func(rets, lambda)))

##    user  system elapsed 
##   2.293   0.009   2.302

system.time(replicate(10000, ewma.filter(rets, lambda)))

##    user  system elapsed 
##   1.062   0.009   1.072

Here it has the full code.