Variance of the Sum of Correlated Variables in R I'm looking to compute $n\text{Var}\left(\frac{1}{n}\sum_{i=1}^nX^{(i)}\right) = \frac{1}{n}\sum_{i=1}^n\sum_{j=1}^n\text{Cov}\left(X^{(i)},X^{(j)}\right)$ in R. Assuming the X's not to be iid, we get
$$ n\text{Var}\left(\frac{1}{n}\sum_{i=1}^nX^{(i)}\right) = \text{Var}\left(X^{(1)}\right) + 2\sum_{k=1}^{n-1}\left(1-\frac{k}{n}\right)\text{Cov}\left(X^{(i)},X^{(i+k)}\right).$$
But as far as I can see, $\text{Cov}\left(X^{(i)},X^{(i+k)}\right)$ is the auto-covariance $C_{XX}(k)$, i.e.
$$ n\text{Var}\left(\frac{1}{n}\sum_{i=1}^nX^{(i)}\right) = \text{Var}\left(X^{(i)}\right) + 2\sum_{k=1}^{n-1}\left(1-\frac{k}{n}\right)C_{XX}(k).$$
However, when I try to compute this value for a simple Markov chain I get very inconsistent results. So I'm just wondering if anyone could help me figure out where this has gone wrong? Or is this method jut infeasible and I should rather compute several Markov chains and compute the empirical variance?
I must add that I'd like to use this on a MCMC method on a multivariate distribution where I'd like to compute $n\text{Cov}\left(\frac{1}{n}\sum_{i=1}^nX^{(i)}\right)$.
# Generates MC with transition probability matri K
MC <- function(x0, K, nsim) {
  X <- array(dim=nsim)
  X[1] <- x0
  for (i in 2:nsim) {
    X[i] <- sample(c(1,2), 1, prob = K[X[i-1],])
  }
  I <- as.numeric(X == y)
}

# Computes nCov(1/n sum X)
ncov <- function(x) {
  n <- dim(x)[1]
  a <- acf(x, plot=FALSE, type='covariance', lag.max=Inf)$acf
  c <- sum(c(1,2*(1-seq(1,n-1)/n))*a)
  return(c)
}

K = matrix(c(0.9,0.2,0.1,0.8),nrow=2)
x <- matrix(MC(x0, K, 100000))
ncov(x)

 A: What you're doing, although theoretically valid, is not computationally stable.  This is because although your estimator is a consistent estimator, for large $k$ there are very few samples available to estimate the lag covariance, leading to bad estimates of this value.
There has been a rich study of estimating this variance term in the MCMC literature. There are estimators that weight and truncate the lag $k$, and these estimators are called spectral variance estimators. They have been shown to be strongly consistent estimators in same Flegal et.al (2010) paper, however, they are a tad bit slow to compute.
There is also the faster batch means estimator, mentioned in the same paper as above. They have also been shown to be consistent (in the same paper).
However, the above paper is only meant for univariate quantities, so these do not assess the covariance for a Markov chain from a multivariate distribution. You can find details of this covariance matrix in the post here. There has recently been work done in estimating the variance-covariance matrix for Markov chains from multivariate distributions. Extensions of the spectral variance estimators are in Vats et.al (2018) and extensions of the batch means estimator are in Vats et. al. (2017). Note that for this multivariate setting the spectral variance estimators are often too computationally expensive.
You can implement all these estimators using the R package mcmcse. This post has some more information about the implementation.
