Autocorrelation of a noisy linear map

Question

I am interested in calculating the autocorrelation function of a linear map with some noise (model given below) but am slightly confused in doing so.

At first, I did not realize there were two definitions of the autocorrelation function: one where you subtract the mean of the time series data (often used in statistics), and one where you do not.

Let $x(n) = \frac{1}{2}\left(a x(n) + b + \xi\right)$ where $\xi$ is some Gaussian noise with mean zero and standard deviation $\sigma$. One can analytically calculate the autocorrelation function as follows:

$$ \langle x(n+1)x(n)\rangle = \langle \frac{1}{2} (a x(n) + b + \xi_n)x(n)\rangle\\ = \frac{a}{2}\sigma_x^2 + \frac{b}{2}\mu_x $$

Where we have used that the noise is uncorrelated with the time series data by assumption.

Similarly,

$$ \langle x(n+2)x(n) \rangle = \frac{a}{2}\langle x(n+1)x(n) \rangle + \frac{b}{2}\langle x(n) \rangle \\ = \left(\frac{a}{2}\right)^2 \sigma_x^2 + \frac{a}{2}\frac{b}{2}\mu_x + \frac{b}{2}\mu_x $$

And then generalizing,

$$ \langle x(n+\tau)x(n) \rangle = \left(\frac{a}{2}\right)^\tau\sigma_x^2 + \frac{b}{2}\mu_x\sum_{l=0}^{\tau-1}\left(\frac{a}{2}\right)^l\\ = \left(\frac{a}{2}\right)^\tau\sigma_x^2 + \frac{b}{2}\mu_x\frac{1 - \left(\frac{a}{2}\right)^\tau}{1-\frac{a}{2}} $$

Where I've summed the geometric series to arrive at the final result.

Autocorrelations are computed using NumPy's fftconvolve function as acf = fftconvolve(x, x[::-1]) (i.e., convolve the time series data with itself reversed, for anyone not familiar with Python).

When I subtract the mean, I get an exponentially decaying autocorrelation function, as depicted in the plot below. This makes sense, given the first term in the analytically derived expression above. Clearly $\mu_x = 0$ after we have subtracted off the mean.

However, originally, I did not realize I was supposed to subtract the mean, and I found a linearly decaying plot as depicted below. By increasing the duration of the simulation, the line extends seemingly to infinity. Is there a simple explanation for this? It does not fall out of the analytical calculation, so I assume there must be something going on because it is a discrete time series? I am new to this, but I assume there is a simple explanation for this phenomenon.

Answer 1

I believe the issue is just that when using NumPy you're using finite duration data.

This example shows what I mean. It starts with a 100-vector of 1s and takes the correlation as you suggest. The resulting output is definitely not a constant, but a triangle.

This is because, in python, you're using finite duration data. So any constant (DC) value will show up like this.

 #
 # Example for correlation of a "constant" value.
 #
 #
 import numpy
 import scipy.signal

 x = numpy.empty(100); x.fill(1)

 acf = scipy.signal.fftconvolve(x, x[::-1])

 print acf

 # OUTPUT:
 #
 #[  1.   2.   3.   4.   5.   6.   7.   8.   9.  10.  11.  12.  13.  14.  15.
 #  16.  17.  18.  19.  20.  21.  22.  23.  24.  25.  26.  27.  28.  29.  30.
 #  31.  32.  33.  34.  35.  36.  37.  38.  39.  40.  41.  42.  43.  44.  45.
 #  46.  47.  48.  49.  50.  51.  52.  53.  54.  55.  56.  57.  58.  59.  60.
 #  61.  62.  63.  64.  65.  66.  67.  68.  69.  70.  71.  72.  73.  74.  75.
 #  76.  77.  78.  79.  80.  81.  82.  83.  84.  85.  86.  87.  88.  89.  90.
 #  91.  92.  93.  94.  95.  96.  97.  98.  99. 100.  99.  98.  97.  96.  95.
 #  94.  93.  92.  91.  90.  89.  88.  87.  86.  85.  84.  83.  82.  81.  80.
 #  79.  78.  77.  76.  75.  74.  73.  72.  71.  70.  69.  68.  67.  66.  65.
 #  64.  63.  62.  61.  60.  59.  58.  57.  56.  55.  54.  53.  52.  51.  50.
 #  49.  48.  47.  46.  45.  44.  43.  42.  41.  40.  39.  38.  37.  36.  35.
 #  34.  33.  32.  31.  30.  29.  28.  27.  26.  25.  24.  23.  22.  21.  20.
 #  19.  18.  17.  16.  15.  14.  13.  12.  11.  10.   9.   8.   7.   6.   5.
 #   4.   3.   2.   1.]

This is why sometimes, rather than using the biased estimator of autocorrelation: $$ \hat\rho[m] = \frac{1}{N} \sum_{n=0}^{N-1-|m|} x[n+m]x[n] $$ the unbiased estimator is used: $$ \hat\rho_{\tt u}[m] = \frac{1}{N-|m|} \sum_{n=0}^{N-1-m} x[n+|m|]x[n] $$ which effectively applies an inverse triangular window to the biased estimator.