One way of including autocorrelation in error estimates is to use blocking in combination with bootstrap of jack-knife. An alternative is given by computing the so-called autocorrelation function and the autocorrelation time following the ideas by Sokal.
The autocorrelation (also autocovariance) function \(C_X(t)\) of a stochastic process \(X_t\) is defined as \[\begin{equation} C_{X}(t)\ = \ \braket{(X_i-\mu)(X_{i+t}-\mu)} \end{equation}\] and its normalised version \(\Gamma_X\) \[\begin{equation} \Gamma_{X}(t)\ =\ \frac{C_{X}(t)}{C_{X}(0)}\,, \end{equation}\] which is the one we will work with. As it is apparent from this definition, \(C_X\) does not depend on \(i\), but only on (time) difference \(t\). Typically, the autocorrelation function decays exponentially like \[\begin{equation} C_{X}(t)\ \propto\ \sum_{n=0}^\infty A_n\ e^{-|t|/\tau_n}\,,\quad \tau_0>\tau_1>...\in\mathbb{R}^+\,. \end{equation}\] This does not necessarily need to be the case, but in the following we will always assume this to be the case. Then, for large enough \(|t|\) only the largest correlation time \(\tau_0\) survives. Hence, one defines the so-called exponential autocorrelation time as follows: \[\begin{equation} \tau_{\text{exp}}\ =\ \limsup_{t \rightarrow \infty}\frac{t}{-\ln|\Gamma_{X}(t)|}\,. \end{equation}\] In terms of \(C_X\) we can compute the variance of \(T_{\bar x}\), the standard estimator of the mean of \(X_t\): \[\begin{equation} \begin{split} \mathrm{Var}\left(T_{\bar x}\right)\ &=\ \frac{1}{N}\mathrm{Var}\left(X_1 + X_2 + ... + X_N\right)\\ &=\ \frac{1}{N^{2}}\sum_{r,s=1}^{N}C_X(r-s)\,.\\ \end{split} \end{equation}\] At this point one usually argues that terms \(\braket{X_i,X_{i+t}}\) are identically zero for \(t\neq 0\). This is no longer the case now, since the data is autocorrelated. Therefore, we have to proceed by including all terms in the sum and counting multiplicities for identical \(t=r-s\): \(t\) can take integer values in the range from \([-N+1,N-1]\). \(t=-N+1\) appears once, like \(t=N-1\). \(-N+2\) appears twice, much like \(N-2\). It is easy to see that this generalises to: \(t\) appears with a multiplicity of \(N-|t|\). Thus, we can write \[ \begin{split} \mathrm{Var}(T_{\bar x})\ &=\ \frac{1}{N^2} \sum_{t=-(N-1)}^{N-1} (N-|t|)C_{X}(t)\\ &=\ \frac{1}{N} \sum_{t=-(N-1)}^{N-1} \left(1-\frac{|t|}{N}\right)C_{X}(t)\\ &=\ \frac{1}{N}C_{X}(0)\sum_{t=-(N-1)}^{N-1}\left(1-\frac{|t|}{N}\right)\Gamma_{X}(t)\,.\\ \end{split} \] If we assume now that \(N\) is much larger than the autocorrelation time \(\tau_0\) in \(X_t\) and that \(\Gamma(t)\propto\exp(-|t|/\tau)\), we can make the following approximation \[\begin{equation} \mathrm{Var}(T_{\bar x})\ \stackrel{N \gg \tau}{\approx}\ \frac{2}{N}C_{X}(0)\ \tau_{\mathrm{int},X}\,, \end{equation}\] with the following definition: Integrated Autocorrelation Time \[\begin{equation} \label{eq:tauint} \tau_{\mathrm{int},X}\ =\ \frac{1}{2}\sum_{t=-\infty}^{\infty}\Gamma_{x}(t)\,. \end{equation}\] The factor \(1/2\) is conventionally inserted to obtain \(\tau_\mathrm{int}\approx\tau_\mathrm{exp}\) for \(\Gamma(t)\propto\exp(-|t|/\tau)\) with \(\tau\gg 1\). The error we make in this approximation is of the order \(\tau/N\). At this point one notes that \(C_X(0)=\mathrm{Var}(X)\), and, therefore, in case of autocorrelation the variance of \(\bar x\) is roughly a factor \(2\tau_{\text{int},X}\) larger than it would be if the \(X_t\) were independent. Stated differently, in case of autocorrelation the effective number of measurements is \(N/(2\tau_{\text{int},X})\).
An estimator for the autocorrelation function \(C(t)\) is given as follows \[\begin{equation} \bar{C}_{X}(t)\ =\ \frac{1}{N-|t|} \sum_{i=1}^{N-|t|} (x_{i}-\bar{x}_{N}) (x_{i+|t|}-\bar{x}_{N})\\ \end{equation}\] which has a bias of the order \(1/N\) and for \(\Gamma(t)\) \[\begin{equation} \bar\Gamma(t)\ =\ \frac{\bar C(t)}{\bar C(0)}\,. \end{equation}\] The variance of this estimator (see e.g. Refs. Priestley (1981); Anderson (2011)) can be computed to \[\begin{equation} \label{eq:gammavar} \mathrm{Var}(\bar C_X(t))\ =\ \frac{1}{N}\sum_{i=-\infty}^\infty [C^2(i) + C(i+t)C(i-t) +\kappa(t, i, i +t)] + o\left(\frac{1}{N}\right) \end{equation}\] with \(\kappa\) the connected only 4–point autocorrelation function. The integrated autocorrelation time should then be estimated using the following (biased) estimator \[\begin{equation} \overline{\tau}_{\text{int}}\ =\ \frac{1}{2}\sum_{t=-W}^{W} \frac{\overline{C}(t)}{\overline{C}(0)}\,, \qquad \tau \ll W \ll N \end{equation}\] because the variance of this estimator grows like \[ \mathrm{Var}\left(\overline{\tau}_{\text{int}}\right)\ \approx\ \frac{2(2W+1)}{N} \overline{\tau}_{\text{int}}^{2}\,. \] This windowing procedure introduces a bias \[ \operatorname{bias}\left(\overline{\tau}_{\text{int}}\right)\ =\ -\frac{1}{2}\sum_{|t|>W} \Gamma(t) + o\left(\frac{1}{N}\right)\,. \] However, we have to live with this bias because otherwise the relative variance of the estimator becomes order 1. A practical way to chose the cutoff \(W\) is to chose it as the first \(t\)–value where \(\Gamma(t)\) becomes zero within errors. The variance of \(\bar\Gamma\) can be estimated using and one finds \[\begin{equation} \label{eq:dGamma} \mathrm{Var}(\bar\Gamma(t))\ \approx\ \frac{1}{N}\sum_{i=1}^{t+\Lambda} \left[\bar\Gamma(i+t) + \bar\Gamma(i-t) - 2\bar\Gamma(i)\bar\Gamma(t)\right]^2\,, \end{equation}\] where it was used that the dominant contribution to the variance comes from the disconnected contributions Madras and Sokal (1988). So, only those were kept and again a suitable windowing procedure was adopted Lüscher (2005). We refer to appendix C of Ref. Madras and Sokal (1988) and to Ref. Sokal (1996) for more details.
Write a function which gets an (autocorrelated) scalar time series as input and computes an estimate for the autocorrelation function \(\Gamma(t)\) as defined above. Note that many scripting languages (like e.g. ‘R’) include such a function. In ‘R’ the corresponding function is called ‘acf’. The implementation should include a parameter ‘W.max’ for the maximal lag to compute the autocorrelation function for.
Estimate the integrated autocorrelation time from this autocorrelation function estimate. In order to prevent the summing up of only noise, you have to stop integrating at some point. At this point implement this cut-off as an input parameter ‘W’.
Implement also the estimator of the error for the estimated autocorrelation function, as discussed above. You will need an additional parameter \(\Lambda\) as a cut-off for the error estimate. How can you use the such estimated error to implement an automatic cut-off procedure for \(\tau_\mathrm{int}\)?
Here a possible solution in ‘R’ using the function ‘acf’ provided by ‘R’
myacf <- function(tseries, W.max, Lambda=100) {
N <- length(tseries)
if(missing(W.max)) W.max <- floor(sqrt(N))
Gamma.tmp <- rep(0, times=2*W.max+W.max+1)
dGamma <- rep(0., times=W.max+1)
## for t > W.max we set Gamma to 0 for simplicity
Gamma.tmp[1:(W.max+1)] <- acf(tseries, lag.max = W.max, plot=FALSE)$acf
## now we determine the error using (E.11) from Luescher, hep-lat/0409106
for(t in c(0:W.max)) {
k <- c(max(1,(t-Lambda)):(t+Lambda))
dGamma[t+1] <- sum((Gamma.tmp[(k+t+1)]+Gamma.tmp[(abs(k-t)+1)]-2*Gamma.tmp[t+1]*Gamma.tmp[(k+1)])^2);
dGamma[t+1] <- sqrt(dGamma[t+1]/N)
}
Gamma <- Gamma.tmp[1:(W.max+1)]
## and we determine where to stop summing Gamma -> W
## for a more sophisticated approach see Wolff, hep-lat/0306017
W <- 0
for(t in c(0:W.max)) {
W <- t
if(Gamma[t+1] < dGamma[t+1]) break
}
## compute integrated autocorrelation time
tau <- 0.5 + sum(Gamma[2:(W+1)])
## Madras, Sokal approximation for the error, (E.14) in Luescher, hep-lat/0409106
dtau <- sqrt((4*W + 2) * tau^2 / N)
res <- list(lags = c(0:(W.max)), Gamma=Gamma, dGamma=dGamma,
W.max=W.max, W=W, tseries=tseries, tau=tau, dtau=dtau)
return(invisible(res))
}Anderson, T.W. 2011. The Statistical Analysis of Time Series. Wiley. https://doi.org/10.1002/9781118186428.
Lüscher, Martin. 2005. “Schwarz-preconditioned HMC algorithm for two-flavour lattice QCD.” Comput. Phys. Commun. 165: 199–220. https://doi.org/10.1016/j.cpc.2004.10.004.
Madras, Neal, and Alan Sokal. 1988. “The Pivot Algorithm: A Highly Efficient Monte Carlo Method for the Self-Avoiding Walk.” Journal of Stat. Phys. 50: 109–86.
Priestley, M.B. 1981. Spectral Analysis and Time Series. Vol. I and II. Academic Press. https://doi.org/10.1002/for.3980010411.
Sokal, Alan. 1996. “Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms.” www.stat.unc.edu/faculty/cji/Sokal.pdf.