Efficient calculation of var-covar matrix in R

I'm looking for efficiency gains in calculating the (auto)covariance matrix from individual measurements over time t with t, t-1, etc..

In the data matrix, each row represents an individual and each column represents monthly measurements (the columns are in time order). Similar to the following data (although with some more co-variance).

# simulate data
set.seed(1)
periods <- 70L
ind <- 90000L
mat <- sapply(rep(ind, periods), rnorm)

Below is the (ugly) code I came up with to get the covariance matrix for measurements/ lagged measurements. It takes almost 4 seconds to run. I'm sure that by moving to data.table, thinking more and not relying on loops I could cut the time by a big amount. But since covariance matrices are ubiquitous I suspect there already exists a standard (and efficient) way to do this in R that I should know about first.

# Get variance covariance matrix for 0-5 lags    
n_lags <- 5L # Number of lags
vcov <- matrix(0, nrow = n_lags + 1L, ncol = n_lags + 1)
for (i in 0L:n_lags) {
  for (j in i:n_lags) {
    vcov[j + 1L, i + 1L] <- 
      sum(mat[, (1L + (j - i)):(periods - i)] *
          mat[, 1L:(periods - j)]) /
      (ind * (periods - j) - 1)
  }
}
round(vcov, 3)

       [,1]   [,2]  [,3]  [,4]  [,5]  [,6]
[1,]  1.001  0.000 0.000 0.000 0.000 0.000
[2,]  0.000  1.001 0.000 0.000 0.000 0.000
[3,]  0.000  0.000 1.001 0.000 0.000 0.000
[4,]  0.000  0.000 0.000 1.001 0.000 0.000
[5,] -0.001  0.000 0.000 0.000 1.001 0.000
[6,]  0.000 -0.001 0.000 0.000 0.000 1.001