Standard error estimation by an automated blocking method

The sample mean $\overline{X}$ is probably the most popular estimator of the expected value in all sciences and $\mathrm{var}\left(\overline{X}\right)$ measures the error (standard- and mean-square-errors). Here, an alternative approach to estimation of $\mathrm{var}\left(\overline{X}\right)$ for time series data is presented. The method has an accuracy similar to dependent bootstrapping, but scales in $O\left(n\right)$ time, and applies to stationary time series, including stationary Markov chains. The computational complexity is bounded by $12n$ floating point operations, but this can be reduced to $n+O\left(1\right)$ in large computations. Convergence in relative error squared is faster than $n^{-1/2}$ and the method is insensitive to the probability distribution of the observations. It is proven that a small part of the correlation structure is relevant to the convergence rate of the method. From this, proof of the Blocking method [Flyvbjerg and Petersen, J. Chem. Phys. 91, 461 (1989)] follows as a corollary. The result is also used to propose a hypothesis test surveying the relevant part of the correlation structure. It yields a fully automatic method which is sufficiently robust to operate without supervision. An algorithm and sample code showing the implementation is available for Python, C++, and R [www.github.com/computative/block]. Method validation using autoregressive AR(1) and AR(2) processes and physics applications is included. Method self-evaluation is provided by bias and mean-square-error statistics. The method is easily adapted to multithread applications and data larger than computing cluster memory, such as ultralong time series or data streams. This way, the paper provides a stringent and modern treatment of the Blocking method using rigorous linear algebra, multivariate probability theory, real analysis, and Fisherian statistical inference.

Figure 1

Whenever $H_0$ is true, the pdf of $M_j$ is known and plotted above. The test concludes that $H_0$ is false if the observed value of $M_j$ is sufficiently unlikely. That is, if the observed value of $M_j$ is larger than $100(1-\alpha )$ percentile for a suitable $\alpha$; the shaded area represents a probability of $\alpha$. The value $d-j=〈M_j〉$ is the expected value of $M_j$ whenever $H_0$ is true since then $M_j\sim {\chi }_{d-j}^2$ is ${\chi }^2$ distributed.

Figure 2

Flow chart of algorithm. The idea is to return the estimate of $\mathrm{var}\left(\overline{X}\right)$, as ${\widehat{\sigma }}_k^2/n_k$ for the smallest value of $k$ such that there is no evidence that ${\gamma }_j\left(1\right)\ne 0$ for $j>k$. This is sensible because then, according to Corollary t7, there is no reason to believe that error $e_k$ is reduced by further iterations of the method.

Figure 3

Relative error squared of two autoregressive models versus observations per time autocorrelation time constant. There is exponential convergence rate for two common correlation structures in natural sciences. In the tail of the distribution, there are outlayers that failed the probabilistic test. The outlayers are not of major concern because (1) they are rare and (2) induced mistakes are small in magnitude. (a), (c) Represent test results for an AR(1) time series with an autocorrelation that is positive with exponential decay, typical of Metropolis-type Markov chains, where it is expected that the observations correlate positively. The process was distributed Gamma(1,1). (b), (d) Represent test results for an AR(2) process. The autocorrelation has exponential decay, but oscillates. The process was distributed multivariate standard normal. It was found that the method was insensitive to the distribution of the observations ($p\ge 0.34$). Consequently, the difference in behavior of the method is attributed to $\gamma \left(h\right)$, as explained by Corollary t7. The expected relative error squared ${\epsilon }^2$ was modeled by gamma regression, $log\left({\epsilon }^2\right)={\beta }_0+{\beta }_{1}log(n/\tau )$. Deviance explained was $50.65\%$ and $65.39\%$ for AR(1) and AR(2) on 6078 degrees of freedom, respectively. Dashed lines give 95% confidence intervals of the expected relative error squared. The plots indicate that it is reasonable to expect the first digit of the method was correct for some $n\gtrsim 20\tau$, and two digits correct for some $n\gtrsim 25000\tau$.

Figure 4

Case study of two textbook physics applications for mean variance estimation. The variance of the mean energy was estimated using manual and automatic Blocking methods and compared with dependent bootstrapping using geometric simulation. The estimates are plotted in (a) and (b), while the autocorrelation is depicted in (c) and (d). (a), (c) Represent test results for a two electron quantum dot with trail energy of a Slater-Jastrow–type state function for a theory of a harmonic oscillator potential with Coulomb repulsion ($\omega =1$). The importance sampling/Hastings-Metropolis theorem implementation implies acceptance rate $>99.999$. The autocorrelation time constant was $\tau =360$. It is clear that $\mathrm{var}\left(\overline{X}\right)=1.46\times {10}^{-8}$ at ${10}^8$ samples according to all three models, and by the above discussion, the error of the variance estimate is correct to the second digit. (b), (d) Represent test results for the Ising model for a $20\times 20$ grid of spins with periodic boundary conditions at temperature $T=2.4$. The energy was sampled from a Boltzman distribution at significance level 95% according to a ${\chi }^2$ test using a Metropolis-type Markov chain. The autocorrelation time constant was measured at $\tau =16$, so at ${10}^6$ samples $\mathrm{var}\left(\overline{X}\right)=1.49\times {10}^{-1}$ according to all three methods, and there is reason to believe that the first two digits are correct.

Figure 5

In the case that the time series is too large for memory, or is so large that it is not saved in a single file, it is possible to reduce the size of each chunk of the time series by repeated use of the matrices $T_j$. This is convenient in the case that the time series is generating by a multithreaded program. The procedure is as follows: Choose one of the chunks of the time series, number $i$. (Step 1) Transform chunk number $i$ by applying the matrix $T_{D-d}{T}_{D-d-1}\cdots T_2{T}_1$. Define $Y_i$ to be the result of this transformation. (Step 2) Repeat for all $j\ne i$. (Step 3) Concatenate all the chunks after they have been formed into one long vector ${({\mathrm{Y}}_1,{\mathrm{Y}}_2,\cdots ,{\mathrm{Y}}_{2^k})}^{\mathsf{T}}$. This vector will now have size $2^d$, thus is small enough to (Step 4) be handled on a single node by using the algorithm of Fig. 2. See Sec. 3d for more details including the definitions of the numbers $D,d$, and $k$.

Figure 6

Python implementation of the algorithm. The code is purposefully verbose to aid implementation in languages of lower level of abstraction, such as C. In practice, the implementation can be optimized and shrunk to about 10–15 lines of code. The most recent implementations for Python, C++, and R are available (see [8]).