EWMA Volatility in R using data frames

Question

I am trying to get EWMA volatility from a series of stock daily returns from a data frame called base_retorno_diario

  Data        IBOV         ABEV3       AEDU3 ALLL3   BBAS3        BBDC3        BBDC4
1 2000-01-04 -0.063756245  0.00000000     0     0 -0.029935852 -0.080866107 -0.071453347
2 2000-01-05  0.024865308 -0.03762663     0     0 -0.008082292  0.043269231  0.060889055
3 2000-01-06 -0.008510238 -0.03157895     0     0  0.014074074  0.014285714  0.008098592
4 2000-01-07  0.012557359 -0.02484472     0     0 -0.022644266  0.017719219  0.000000000
5 2000-01-10  0.043716564  0.00000000     0     0  0.050074738  0.005357143  0.006985679
6 2000-01-11 -0.026401514 -0.02388535     0     0 -0.008540925 -0.059058615 -0.046479362

First row of the new data frame (n_row and n_col is the number of rows and columns on the returns data frame base_retorno_diario)

EWMA_VARIANCE = as.data.frame( base_retorno_diario[1,2:n_col]^2 )

then I created the following loop

i = 2
while(i<=n_row){
    EWMA_VARIANCE = rbind(EWMA_VARIANCE,
        EWMA_VARIANCE[(i-1), 1:(n_col-1)] * DECAY_FACTOR +
        (1-DECAY_FACTOR) * base_retorno_diario[i,2:n_col]^2
    )
    i=i+1
}

It works fine but it is taking too long (the original data frame has 3560 obs of 101 variables), is there anyway to avoid the loop in this case ? DECAY_FACTOR = 0.97

Answer 1

You can avoid this loop with some matrix algebra. Let's assume the raw data is a vector (a_1, a_2, a_3, ..., a_n) and we want to create the EWMA variance (x_1, x_2, x_3, ..., x_n) according to your definition. Let d be the decay factor. If i understood your code correctly, you currently have a recursive definition

which makes things difficult. I believe this non-recursive definition is identical

This allows us to take advantage of some linear algebra to get the job done with matrix multiplication. For brevity, I will assign shorter variable names to your data.frame and decay factor

dd <- base_retorno_diario
d <- DECAY_FACTOR

Now we begin by calculating all the squared values first, and then take the pairwise difference that we can see are part of the non-recursive definition.

asquare <- as.matrix(dd[,2:7])^2
asqdiffs <-sapply(data.frame(asquare), diff)

Now we create an appropriate matrix with the values of d to take are of the summing part of the non-recursive definition and then perform the subtraction (with a little offset for the initial term)

dx <- outer(1:nrow(asqdiffs), 1:nrow(asqdiffs), FUN=function(x,y) 
    ifelse(x>=y, d^(x-y+1),0 )
)
EWMA_VARIANCE <- asquare - rbind(0, dx %*% asqdiffs)

This method seem to produce the same results are yours, but it is about 20x faster in my tests.