# Vectorize/speed-up exponentially weighted covariance calculation in R

I have the following 2-column data.frame:

Data <- data.frame(matrix(rnorm(100),ncol = 2))


I would like to calculate the variance-covariance matrix as the exponentially weighted average of the past squared/cross observations (on an expanding window):

Where delta is set to be 60/61. Thus, at t=4, the sum runs from 1 to 4, and the averages are computed over the first 4 observations. I have managed to implement a function calculating the upper triangular part of the covariance matrix, but it is quite slow. Does anyone have ideas to improve the speed?

rollEWCov <- function(Data){
res <- c()
for(i in 1:nrow(Data)){ # Expanding window of data used to calc covariances
means <-  colMeans(Data[1:i,])
Cov <- matrix(0,nrow = ncol(Data),ncol=ncol(Data))
for(k in 1:ncol(Data)){ # first data-column
l <- k
while(l <= ncol(Data)){ # second data-column
Sum <- 0
for(j in 1:i){ # calc sum of exponentially weighted average of past returns
Sum <- Sum + (1-60/61)*(60/61)^(j-1)*(Data[i-(j-1),k] - means[k])*(Data[i-(j-1),l] - means[l])
}
Cov[k,l] <- Sum
l <- l + 1
}
}
res[[i]] <- Cov
}
return(res)
}


A key performance bottleneck in your code is your calculation of the Cov_t^{k,l} values:

Sum <- 0
for(j in 1:i){ # calc sum of exponentially weighted average of past returns
Sum <- Sum + (1-60/61)*(60/61)^(j-1)*(Data[i-(j-1),k] - means[k])*(Data[i-(j-1),l] - means[l])
}
Cov[k,l] <- Sum


You are summing from j=1 through j=i by looping through the values individually and adding them up. In R you can get significant speedups by using the vectorized sum function:

Cov[k,l] <- sum((1-60/61)*(60/61)^(0:(i-1))*(Data[i:1,k] - means[k])*(Data[i:1,l] - means[l]))


Additionally, this part of the code hard-codes the choice of delta=60/61, which may make it harder to see what the code is doing or update to another delta value. It would be best to pass delta to the function and then use delta instead of 60/61 in this expression. Here is the updated function:

rollEWCovMod <- function(Data, delta){
res <- c()
for(i in 1:nrow(Data)){ # Expanding window of data used to calc covariances
means <-  colMeans(Data[1:i,])
Cov <- matrix(0,nrow = ncol(Data),ncol=ncol(Data))
for(k in 1:ncol(Data)){ # first data-column
l <- k
while(l <= ncol(Data)){ # second data-column
# calc sum of exponentially weighted average of past returns
Cov[k,l] <- sum((1-delta)*delta^(0:(i-1))*(Data[i:1,k] - means[k])*(Data[i:1,l] - means[l]))
l <- l + 1
}
}
res[[i]] <- Cov
}
return(res)
}


To see the performance improvements of the vectorized sum, we can benchmark with a 500 x 2 matrix:

Data <- data.frame(matrix(rnorm(1000),ncol = 2))
system.time(rollEWCov(Data))
#    user  system elapsed
#  11.300   0.103  12.031
system.time(rollEWCovMod(Data, 60/61))
#    user  system elapsed
#   0.142   0.006   0.149
all.equal(rollEWCov(Data), rollEWCovMod(Data, 60/61))
# [1] TRUE


We obtained a 80x speedup and got identical results to the original function.

While using sum() helps a lot (as josliber did), you can get another 2 orders of magnitude of performance improvement by removing two of the loops. You can replace them with a call to crossprod(). Using a matrix instead of a data.frame is also helpful for performance. Functions and timings below (the other two functions are also byte-compiled).

rollEWCovMod2 <- compiler::cmpfun(function(Data, delta) {
mat <- as.matrix(Data)  # matrix is often faster than data.frame
nc <- ncol(mat)         # only calculate number of cols once
nr <- nrow(mat)         # only calculate number of rows once
out <- vector("list", nr)  # pre-allocate output

# Expanding window of data used to calc covariances
for (i in seq_len(nr)) {
# sequence of observation locations
one..i <- seq_len(i)
# subset matrix one time
mat_subset <- mat[one..i, , drop = FALSE]
# column means
means <- matrix(.colMeans(mat_subset, i, nc), i, nc, byrow = TRUE)
# subtract the column mean from each column
scaled <- mat_subset - means
#scaled <- scale(mat_subset, center = TRUE, scale = FALSE)
# calculate exponential weights
wts <- (1-delta) * delta^(-one..i + i)
# caclulate crossproduct
res <- crossprod(scaled*wts, scaled)
dimnames(res) <- NULL     # remove row and column names
res[lower.tri(res)] <- 0  # set lower values to zero
# store result
out[[i]] <- res
}
return(out)
})


Reproducible example and timings.

set.seed(21)
Data <- data.frame(matrix(rnorm(4000), ncol = 20))
delta <- 60/61
system.time(r1 <- rollEWCov(Data))
#    user  system elapsed
# 120.788   0.008 120.879
system.time(r2 <- rollEWCovMod(Data, delta))
#    user  system elapsed
#     2.1     0.0     2.1
system.time(r3 <- rollEWCovMod2(Data, delta))
#    user  system elapsed
#   0.024   0.000   0.025
all.equal(r1, r2)
# [1] TRUE
all.equal(r1, r3)
# [1] TRUE
all.equal(r2, r3)
# [1] TRUE