(BM) estimators [@chen:seila:1987], spectral variance estimators [@hannan:1970], regenerative estimators [@seil:1982], and initial sequence estimators [@dai:jon:2017; @koso:2000]. [@vats:rob:fle:jon:2020] recommend the BM estimator due to its ease of implementation, universality of application, computational efficiency, and well established asymptotic properties.

Let $\hat{\Sigma}_k$ denote the BM estimator of $\Sigma$ from the $k$th chain. A default combined estimator of $\Sigma$, which we call the average batch means (ABM) estimator is, $\hat{\Sigma} = m^{-1} \sum_{k=1}^{m} \hat{\Sigma}_{k}\,.$ ABM naturally retains asymptotic properties of $\hat{\Sigma}_k$, since the $m$ chains are independent. This averaging of the variance is reasonable for independent and identically distributed samples. However, before each chain has adequately explored the state space, chains may dominate different parts of ${\mathcal{X}}$, so that each $\hat{\Sigma}_{k}$ significantly underestimates the truth. As a result, the ABM estimator also underestimates the truth and many of the benefits of running parallel chains from dispersed starting values are lost.

In the steady-state literature for univariate $\phi$-mixing processes ($p = 1$), [@arg:and:2006] introduced a *replicated batch means* (RBM) estimator that estimates $\Sigma$ by essentially centering the chain around the overall mean $\hat{\mu}$, instead of the $k$th chain mean, $\hat{\mu}_k$. We develop the RBM estimator in this multivariate setting and obtain (strong) consistency, large-sample bias 