# Calculating error of mean of time series

This is probably a stupid question. I'm mathematically literate, but my knowledge of statistical methods is embarrassingly hopeless.

I have multiple time series representing various quantities for different systems, obtained through computational (molecular dynamics) simulation. I would like to state the mean and standard deviation of each of these quantities for each system, along with an estimate of the error of the mean of the time series.

I've been averaging the values over all time steps to get the mean, and calculating standard deviation using the basic sample formula. I've then used this standard deviation to obtain the standard error of the mean as $$\mathrm{SE}_{\bar{x}} = \frac{\sigma_x}{\sqrt{n}}$$ (where $\sigma$ = standard deviation and $n$ = number of time steps, which is very large).

But it occurs to me that the data in a time series is likely to be correlated, so this isn't a valid way to obtain the standard error of the mean. Is this correct?

If so, that's where I get stuck. What statistical tests can I perform to obtain an estimate of the error of the mean?

• Look at ARIMA
– Carl
Apr 19 '17 at 18:33

Clearly you have good statistical intuition, because you are exactly right! Because of correlations between the individual terms, the standard error of the mean of the observations is not an accurate estimate of the error bar on the population mean from time series data.

The actual variance of the sample mean $m$ is $$(\delta m)^2 = \frac{1}{n} \left[ g_0 + 2 \sum_{k=1}^{n-1} \frac{n-k}{n} g_k \right]$$ where $g_{k}$ is the co-variance between $x_i$ and $x_{i-k}$.

It turns out to be a bit of a pain to apply this result. If you just plug in the estimated co-variances for the $g_k$, you get quite wrong results, essentially because of correlations between the $g_k$ estimators. There are a number of different ways you can proceed, with various pros and cons. One relatively simple approach without too much downside is to just drop the higher co-variances, for which you don't have good estimates anyway; it turns out that using a cutoff $k_{\rm max} \approx \sqrt{n}$ works out alright. For more discussion of these issues, see papers like Ryo Okui, “Asymptotically Unbiased Estimation of Autocovariances and Autocorrelations with Long Panel Data”, Econometric Theory (2010) 26: 1263.

An earlier comment-er suggested just doing an ARIMA fit and taking the error bar on the mean from the error bar on the mean parameter of the ARIMA model. That's fine if the data are actually well-fit by an ARIMA model. But the approach I am suggesting here is model-independent.

• You are giving too much credit for a three-word comment. But what you suggest in your interpretation of the comment is valuable. Apr 20 '17 at 11:38
• Are both approaches actually compatible or do they reflect a different view of the data? (I saw this question and answer while formulating my own, related question: stats.stackexchange.com/questions/322289.) Jan 9 '18 at 17:30

I'm answering my own question for the reference of future people who might find it helpful. I've accepted David Wright's answer, though, because it contains the actual solution and he did all the work.

In this case, I used ARIMA, specifically the auto.arima function in the R package forecast. I wanted to be able to implement my solution computationally due to the large size of my data sets, so R was very useful. Using this function, the mean and standard error thereof were simply printed to the screen, saving me a lot of work.

Thank you, Carl, for mentioning ARIMA in the first place, and David Wright for expanding upon the process. You were both very helpful. I also enjoyed David's explanation of the model-independent method, and I'll definitely look into that more.

I was only able to obtain an estimate for the mean when $d=0$ (that is, the ARIMA model featured no differencing). That makes sense, since the mean would be reduced to zero (or something close?) by differencing, but it's worth mentioning.

As it transpired, the optimal ARIMA model (restricted to $(p, 0, q)$, as above) for most — but not all — of my data sets turned out to be $(0, 0, 0)$. If I'm not mistaken (which I may well be!), this means the variation in the data can be represented by white noise and wasn't correlated after all, so I was being over-cautious. Regardless, it's good to have confirmed that, and this question might prove useful to someone with a more strongly correlated time series at some point.

Maybe "the blocking method" proposed by H. Flyvbjerg and H. G. Petersen in "Error Estimates on Averages of Correlated Data", J. Chem. Phys. 91, 461 (1989) would be useful.