Vectorize/speed-up exponentially weighted covariance calculation in R
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.