# Sample autocorrelation function from multiple realisations of a stochastic processs using the convolution theorem

Consider a wide-sense stationary stochastic process $X(t)$ with constant mean $\mu$ and standard deviation $\sigma$. The (auto)correlation function $\rho(\tau)$ of $X(t)$ is defined as:

$$\rho(\tau) = \frac{\mathbb{E}[(X(t)-\mu)(X(t+\tau)-\mu)]}{\sigma^2}$$

My question relates to the estimation of the correlation function from a discretised realisation $\boldsymbol{X} = \{X_1,X_2,\dots,X_n\}$ of $X(t)$. A simple estimator $\hat{\rho}(\tau)$ is to compute Pearson's correlation coefficient for all pairs of lags separated by $k$:

$$\hat{\rho}(k)=\frac{1}{(n-k) \sigma^2} \sum_{t=1}^{n-k} (X_t-\mu)(X_{t+k}-\mu)$$

Now suppose that instead of a single realisation of $X(t)$ being available there are $m$ realisations of $X(t)$ available. The estimator is easily extended to this case by considering all discretised points seperated by lag $k$ for all realisations:

$$\hat{\rho}(k)=\frac{1}{j(n-k) \sigma^2} \sum_{j=1}^m \sum_{t=1}^{n-k} (X_{t}^j-\mu)(X_{t+k}^j-\mu)$$

However these estimators are computationally inefficient in comparison to an approach depending on the signal processing definition of autocorrelation as the convolution of the signal with itself:

$$\rho(\tau) = X(t) * \bar{X}(t) = \int_{-\infty}^\infty X(t)X(t-\tau)\, {\rm d}t$$

It is then possible to use the convolution theorem to write:

$$X(t) * \overline{{X}(t)} = \mathcal{F}^{-1}\big\{\mathcal{F}\{X(t)\}\cdot\overline{\mathcal{F}\{X(t)\}}\big\}$$

Where $\mathcal{F}\{X(t)\}$ is the Fourier transform of $X(t)$ and $\overline{\mathcal{F}\{X(t)\}}$ is its complex conjugate. This leads to a computationally efficient estimator for $\rho(\tau)$ that makes use of optimised fast Fourier transform algorithms:

$$\hat{\rho}(k) = \mathcal{F}^{-1}\big\{\mathcal{F}\{\boldsymbol{X}\}\cdot\overline{\mathcal{F}\{\boldsymbol{X}\}}\big\}$$

Where $\boldsymbol{X}$ is a discretised single realisation of $X(t)$ (as above). My question is whether it possible to extend this approach to the case where there are $m$ realisations of $X(t)$? Say that $\boldsymbol{X}$ is an $n$ by $m$ matrix instead of a vector.

• If $\mathbb{E}[X_t]$ and $R_k = \mathbb{E}[(X_t-\mu)( X_{t+k}-\mu)]$ is constant then $R_k =\lim_{A \to \infty} \frac{1}{2A} \mathbb{E}[\int_{-A}^A (X_t-\mu)( X_{t+k}-\mu)dt] = \lim_{N \to \infty} \frac{1}{2N} \mathbb{E}[\sum_{n=-N}^N (X_n-\mu)(X_{n+k}- \mu)]=\frac{1}{2N} \sum_{n=-N}^N \mathbb{E}[(X_n-\mu)(X_{n+k}- \mu)]$. Aug 29 '17 at 16:12
• For this last one, a good estimator is $\widehat{R}_k=\frac{1}{2N m}\sum_{j=1}^m \sum_{n=-N}^N (x_n^{(j)}-\mu)(x_{n+k}^{(j)}- \mu)$ where $x^{(j)},j=1 \ldots m$ are realizations of your process. You can compute $\widehat{R}_k$ (for $k = 1 \ldots K$) efficiently from the FFT of $\frac{1}{m} \sum_{j=1}^m x^{(j)}_n-\mu, n = -N \ldots N-K$ Aug 29 '17 at 16:13
• Did you get what I wrote ? Aug 31 '17 at 22:33