Systematic negative bias in auto-correlations of random time-series

Question

I generated completely random time-series (using MATLAB "rand" function or Python "numpy.random.rand" function) and computed the auto-correlation of these time-series using either MATLAB xcorr, or Python numpy.correlate or just the definition of correlation coefficient. Before, computing the correlation I subtracted the signal mean from it.

For a random time-series we expect to have zero correlations; however, when I generate 1000 random arrays each with 10000 samples and compute the average auto-correlation for each of them, the distribution of average correlation of trials is always biased toward negative values. I tried this with both MATLAB and Python and even for more or longer trails but the negative bias is always there (although reduces a bit by increasing the system size).

I was wondering if anyone knows where this bias is coming from and how I can remove it?

Answer 1

Thanks for all the answers, after searching through literature I found out, indeed there is a systematic bias in estimating the autocorrelation of time-series with a finite-size. There is a series of old statistics papers discussing this bias and even they derived its analytical form for some cases like: Marriott, F. H. C., and J. A. Pope. "Bias in the estimation of autocorrelations." Biometrika 41.3/4 (1954): 390-402 (https://www.jstor.org/stable/2332719)

It seems, there are two major sources for the bias: - mean estimation of the finite-size time-series - correlation between covariance and variance (normalization part)

If we have a long enough time-series to better estimate the mean or use the true mean that we generated the data with, the negative bias that I observed in figures above will disappear.

Answer 2

The effect of subtracting the mean...

... the same type of correlation occurs in residuals of a linear regression, which are correlated.

An image that illustrates this is in the question Why are the residuals in $\mathbb{R}^{n-p}$?

In the example below, you see an illustration for the fitting of $\mathbf{y} = a + b\mathbf{x}$ with only three points.

The error is a vector perpendicular to the surface spanned by $x_1$ and $x_2$. For any observation, the error will point in the same direction and can be seen to be a multiple of a line (a 1D space).

In the above image, no matter what the data-point is, the residual will be along a 1D vector.

In a 2-sample case the effect is the strongest and the residuals are fully correlated $$\begin{array}{rrr}x_1 - \bar{x}& = &\frac{1}{2} (x_1 - x_2) \\ x_2 - \bar{x} &= & -\frac{1}{2} (x_1 - x_2)\end{array}$$.

For multiple samples the residuals can be computed with the residual maker matrix

$$A = I - X(X^TX)^{-1}X^T = \begin{bmatrix} 1-\frac{1}{n} & -\frac{1}{n} & -\frac {1}{n}& \dots & -\frac{1}{n} \\ -\frac{1}{n} & 1-\frac{1}{n} & -\frac {1}{n}& \dots & -\frac{1}{n} \\ -\frac{1}{n} & -\frac{1}{n} & 1-\frac{1}{n} & \dots & -\frac{1}{n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -\frac{1}{n} & -\frac{1}{n} & 1-\frac{1}{n} & \dots & 1-\frac{1}{n} \\ \end{bmatrix}$$

For the offdiagonal elements $A_{ij} = -1/n$ and for the diagonal elements $A_{ii} = 1-1/n$.

If you start with Gaussian white noise, then the residuals will be multivariate normal distributed with covariance matrix $$\Sigma_{ij} = AA^T = \begin{cases} 1-\frac{1}{n} & \quad {\text{if $i=j$}}\\ -\frac{1}{n} & \quad {\text{if $i\neq j$}}\\ \end{cases}$$

The correlation between different terms in the time-series will be $$\rho = \frac{cov(x_i,x_j)}{\sqrt{var(x_i)var(x_j)}} = - \frac{1}{n-1}$$

For your case with $n=10000$ we would expect a correlation $\approx 10^{-4}$. I am not sure what causes the discrepancy with your results.