Vectorize/speed-up exponentially weighted covariance calculation in R I have the following 2-column data.frame:

Data &lt;- 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 &lt;- function(Data){      
  res &lt;- c()
  for(i in 1:nrow(Data)){ # Expanding window of data used to calc covariances
    means &lt;-  colMeans(Data[1:i,])
    Cov &lt;- matrix(0,nrow = ncol(Data),ncol=ncol(Data))
    for(k in 1:ncol(Data)){ # first data-column
      l &lt;- k
      while(l &lt;= ncol(Data)){ # second data-column
        Sum &lt;- 0
        for(j in 1:i){ # calc sum of exponentially weighted average of past returns
          Sum &lt;- 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] &lt;- Sum
        l &lt;- l + 1
      }
    }
    res[[i]] &lt;- Cov
  }
  return(res)
}
A key performance bottleneck in your code is your calculation of the Cov_t^{k,l} values:

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


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

Data &lt;- 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.