Vectorize/speed-up exponentially weighted covariance calculation in R

Question

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)
}

Answer 1

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.

Answer 2

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