Determining standard error of the mean from a correlated, stationary time series using known autocorrelation without block averaging

Question

I'd like to determine the SEM of measurements taken from a stationary time series. SEM calculation using all measurements isn't accurate because adjacent measurements may be highly correlated, so the number of independent samples is lower than the number of actual measurements. The use-case is energy measurements from molecular dynamics simulation, but I'm using simulated data to control the autocorrelation first. I'd like to achieve this using the estimated autocorrelation and not by block averaging.

---

Here's a simulated dataset with autocorrelation=0.95, using python:

import numpy as np

T = 1000
y = np.zeros((T,))
autocorr = 0.95
for i in range(1,T):
    y[i] = autocorr * y[i-1] + np.random.normal()

One could then perform block averaging to visualize how the SEM approaches the true value while using larger and larger block sizes. Using skimage.measure.block_reduce to do the blocking shows the following, and you might guesstimate the SEM to be around 0.4:

from scipy.stats import sem
from skimage.measure import block_reduce
import matplotlib.pyplot as plt

plt.plot([sem(block_reduce(y, (i,), func=np.mean)) for i in range(1,200)], '-o')
plt.xlabel('Block size (no. data points)')
plt.ylabel('SEM')

However I'd rather use the autocorrelation since block averaging has a bunch of arbitrary choices - when do you decide the SEM has converged? From which measurement do you start the blocking? Also block averaging requires removing the data that doesn't fit into a block, which seems wasteful.

Alternatively, autocorrelation can be estimated by a bunch of methods, for example in numpy:

ac = np.correlate(y-y.mean(), y-y.mean(), mode='full') #correlation of y against itself
ac /= ac.max() #normalize
ac = ac[int(ac.shape[0]/2):] #only return the autocorrelation starting from timelag 0.
print(ac[1])

out: 0.9304976

Pretty close to the true value. Or perhaps autocorrelation is found by fitting an AR(1) autoregressive model. I used pymc3 to do that but the code isn't essential for the question, and it generates an estimate for the autocorrelation consistent with the above.

---

So at this point, how does one determine the SEM given the estimated autocorrelation? Thanks!

Answer 1

Seems as if "effective sample size" is a suitable concept for you to look at, in particular the formulae for $Var(\hat{\mu})$ and $n_{eff}$.

https://en.wikipedia.org/wiki/Effective_sample_size

A remaining problem however is the fact that your (auto)correlations are not uniform but lag dependent. To solve it you may employ the magnitude of the sample correlation matrix (i.e. the summed entries of the inverse correlation matrix), which is obtained immediately from the results of the autocorrelation analysis.

More concretely, if $\{y_q,...,y_n\}$ denotes the time series in question and $\it{R}$ the corresponding sample correlation matrix, then the optimal weight vector $\mathbf{w} = (w_1,...,w_n)$ yielding the largest effective sample size (and thus the lowest-variance (unbiased) linear estimator of the mean) is given by the solution of: \begin{align} \it{R}\mathbf{w}= \left( \begin{array}{c} 1\\ \vdots\\ 1\\ \end{array} \right) \end{align} giving $ \hat{\mu} := \frac{1}{|\it{R}|}\sum_{i=1}^n w_iy_i$ as the lowest variance estimator of the mean. $|\it{R}| = \sum_{i=1}^n w_i$ denotes the magnitude of $\it{R}$.

The overall approach is nicely explained in great detail here

https://golem.ph.utexas.edu/category/2014/12/effective_sample_size.html

which is also a reference in the above wiki link.

Answer 2

After posting this I went on a long tour of all the techniques to calculate SEM for autocorrelated data like timeseries. I wrote up the full version including a comparison at: https://ljmartin.github.io/technical-notes/stats/estimators-autocorrelated/

The simplest way that is also quite robust is probably to first calculate the autocorrelation (but not by using the autocorrelation function). Use the statsmodels AutoReg function to estimate that parameter, called rho. Then, using the estimated value for rho (autocorrelation), apply a correction factor to the SEM.

The full code for this is:

def correction_factor(rho, n):
    d = ((n-1)*rho - n*rho**2 + rho**(n+1)) / (1-rho)**2
    k = np.sqrt( (1 + (2*d)/n) / ( 1 - (2*d)/(n*(n-1))  ) )
    return k

timeseries = **YOUR TIMESERIES**
fit_result = AutoReg(timeseries-timeseries.mean(), lags = [1]).fit()
estimated_rho = fit_result.params[1]
n=len(timeseries) #actual number of samples
correction = correction_factor(estimated_rho, n) #correction factor
naive_sem = stats.sem(timeseries) #first calculate SEM naively
sem = naive_sem*correction #now correct the naive SEM by multiplying by correction

The correction factor is known analytically from the equation for an autoregressive AR(1) process (see ref below) and works well when you have sampled for at least about 10-20 multiples of the autocorrelation time. Knowing the autocorrelation time is a different question but in practice try halving your data and see if the SEM is that different - if it is then you don't have enough.

Ref:

Bence, James R. "Analysis of short time series: correcting for autocorrelation." Ecology 76.2 (1995): 628-639.