Everything you wish to know about time correlations but are afraid to ask

Abstract

We discuss the various definitions of time correlation functions and how to estimate them from experimental or simulation data. We start with the various definitions, both in real and in Fourier space, and explain how to extract from them a characteristic time scale. We then study how to estimate the correlation functions, i.e. how to obtain a good approximation to them from a sample of data obtained experimentally. Finally we discuss some practical details that arise in the actual computation of these estimates, and we describe some relevant algorithms.

1 Introduction

This chapter is about the definitions and practical computation of time correlation functions, i.e. the mathematical tools that enable us to find and study dynamical correlations in physical systems. Why do we care about correlations? Science is about understanding how things work (or, more ambitiously, how nature works [1]). The question of what “understanding” something really means is not one we plan to answer or even discuss here, but most of us would probably agree that it involves (possibly among other things) knowledge of a causation mechanism: to know how the state of the system at some time influences the behavior of the same system at a later time. Now of course correlation does not imply causality: if $A$ and $B$ are correlated, it might be that $A$ causes $B$, or that $B$ causes $A$, or that something else causes both $A$ and $B$. So correlations do not (directly) tell us about cause and effect. However, causality is not directly measurable, while correlations are. Correlations do not provide us with a causality mechanism, but do constrain the cause-effect relationships we might care to imagine: to explain how a system works, you are free to come up with whatever mechanism (theory) you wish, but if that mechanism does not produce the kind of correlations actually observed, then it cannot be right.

Correlations play a major role in the study of systems of many particles, where the prediction (and maybe even the observation) of detailed particle trajectories is out of the question, due to the large number of variables involved. Instead, one takes a statistical approach, predicting (macroscopic) temporal correlations rather than microscopic trajectories. Conversely, in a macroscopic experiment one does not measure particle trajectories, but records quantities that are the result of the collective instantaneous state of many microscopic components. The exact variations of the observed quantities are due to factors beyond our control (we say they are noisy), but the overall trends and the correlations between the quantities measured at successive times furnish meaningful information about the system’s dynamics. Moreover, time correlations provide information on the dynamics even in the absence of overall trends (e.g. when a physical system is in thermodynamic equilibrium). But time correlations are also greatly useful for the study of systems quite more complex than physical systems in thermodynamic equilbrium: in particular, of special interest for this volume, biological systems.

We will thus proceed to present the mathematical definitions of time correlation functions and their Fourier transforms (§2), to discuss the definition and meaning of the correlation time (§3) and finally to turn to the question of computing time correlations from experimental data (§4).

But before moving on, I would like to make two clarifications. First, this chapter is about questions you are afraid to ask, not because they are so advanced that they touch very dark and well-kept secrets, but because they are so basic that you are embarrassed to ask. No honest scientific question should be embarrassing, but you know. Second, this is perhaps not everything you want to know about time correlations. In particular, I have not included any material regarding their theoretical calculation. But I do discuss all (well, most) of what you need to know to compute them from experimental or simulation data, starting from their various possible definitions and touching many practical and sometimes nasty details. I apologize if these clarifications are disappointing, but you must agree that “almost all you need to know to be able to compute time correlations from experimental data, some of which is so basic that you are embarrassed to ask” would make for considerably less catchy title.

2 Definition of time correlation functions

2.1 Correlation and covariance

Let $x$ and $y$ be two random variables and $p(x,y)$ their joint probability density. We use $⟨\mathrm{\dots }⟩$ to represent the appropriate averages, e.g. the mean of $x$ is $⟨x⟩=\int xp(x)𝑑x$ (recall that the probability distribution of $x$ can be obtained from the joint probability, $p(x)=\int 𝑑yp(x,y)$, and in this context is called marginal probability) and its variance is ${\text{Var}}_x=⟨{(x-⟨x⟩)}^2⟩$. Let us write

$$C_{xy}=⟨xy⟩=\int xyp(x,y)𝑑x𝑑y.$$ (1)

Eq. (1) defines the correlation of $x$ and $y$. Their covariance is defined as

$${\text{Cov}}_{x,y}=⟨\left(x-⟨x⟩\right)\left(y-⟨y⟩\right)⟩=⟨xy⟩-⟨x⟩⟨y⟩.$$ (2)

A property of the covariance is that it is bounded by the product of the standard deviations (provided of course that they exist, which cannot always be taken for granted):

$${\text{Cov}}_{x,y}^2\le {\text{Var}}_x{\text{Var}}_y.$$ (3)

The Pearson correlation coefficient, or Pearson correlation for short, is

$$r_{x,y}=\frac{{\text{Cov}}_{x,y}}{\sqrt{{\text{Var}}_x{\text{Var}}_y}}.$$ (4)

From the inequality (3) it follows that the Pearson coefficient is bounded, $-1\le r_{x,y}\le 1$. It can be shown that the equality holds only when the relationship between $x$ and $y$ is linear [2, §2.12].

The variables are said to be uncorrelated if their covariance is null:

$${\text{Cov}}_{x,y}=0⟺⟨xy⟩=⟨x⟩⟨y⟩\text{(uncorrelated)}.$$ (5)

Absence of correlation implies that the variance of the sum is the sum of the variances, because

$${\text{Var}}_{x+y}={\text{Var}}_x+{\text{Var}}_y+2{\text{Cov}}_{x,y},$$ (6)

but is weaker than independence: independence means that $p(x,y)=p(x)p(y)$. In the particular case when the joint probability $p(x,y)$ is Gaussian, ${\text{Cov}}_{x,y}=0$ is equivalent to $x$ and $y$ being independent, but this not true in general. On the other hand it is clear that independence does imply absence of correlation. The covariance, or the correlation coefficient, can be said to measure the degree of linear association between $x$ and $y$, since it is possible to build a nonlinear dependence of $x$ on $y$ that yields zero covariance (see Ch. 2 of [2]).

2.2 Fluctuating quantities as stochastic processes

Consider now an experimentally observable quantity, $A$, that provides some useful information on a property of a system of interest. We will usually assume here for simplicity that $A$ is a scalar, but the present considerations can be rather easily generalized to vector quantities. $A$ can represent the magnetization of a material, the number of bacteria in a certain culture, the rainfall in a certain area, the prize of a share in the stock market, etc. We assume that $A$ can be measured repeatedly and locally in time, so that we actually deal with a function $A(t)$.

We wish to compare the values of $A$ measured at different times. However, we are interested in cases where $A$ is noisy, i.e. subject to random fluctuations that arise because of our incomplete knowledge of the variables affecting the system’s evolution, or because of our inability control them with enough precision (e.g. we do not know all variables that can affect the variation of the price of a stock market asset, we do not know all the interactions in a magnetic system, we cannot control all the degrees of freedom of a thermal bath that fixes the temperature of a sample). The quantity $A(t)$ is then a random variable, and to compare it at different values of its argument we will resort to the correlation and covariance defined in the previous section. This will lead us to define time correlation functions, and this section is devoted to stating their precise definitions.

Let us stress that a statement like “$A(t_1)$ is larger than $A(t_2)$” is useless in practice, because even if it is meaningful for a particular realization of the measurement process (experiment), the noisy character of the observable means that a different repetition of the experiment, under the same conditions, will yield a different function $A(t)$. Actually repeating the experiment may or may not be feasible depending on the case, but we assume that we know enough about the system to be able to assert that a hypothetical repetition of the experiment would not exactly reproduce the original $A(t)$. For clarity, it may be easier to imagine that several copies of the system are made and let evolve in parallel under the same conditions, each copy then producing a signal slightly different from that of the other copies.

We note that the expression “under the same conditions” is implicitly qualified in some system-dependent way. Clearly we expect that two strictly identical copies of a system evolving under exactly identical conditions will produce the same function $A(t)$. Same conditions here must be understood in a statistical way: the system is prepared by choosing a random initial configuration extracted from a well-defined probability distribution, or two identical copies evolve with a random dynamics with known statistical properties (e.g. coupled to a thermal bath at given temperature and pressure). We are excluding from consideration cases where the fluctuations are mainly due to experimental errors. If that where the only source of noise, one could in principle repeat the measurement enough times so that the average $⟨A(t)⟩$ is known with enough precision. $⟨A(t)⟩$ would then be an accurate description of the system’s actual evolution, and the correlations we are about to study would be dominated by properties of the measurement process rather than by the dynamics of the system itself. Instead we are interested in the opposite case: experimental error is negligible, and the fluctuations of the observable are due to some process intrinsic to the system. Indeed in many cases (such as macroscopic observables of systems in thermodynamic equilibrium) the average of the signal is uninteresting (it’s a constant), but the time correlation function unveils interesting details of the system’s dynamics.

2.2.1 Stochastic processes

From our discussion of $A$ as a fluctuating quantity, it is clear that $A(t)$ is not an ordinary function. Rather, at each particular value of time, $A(t)$ is a random variable, and $A(t)$ as a whole is a random function, or stochastic process. A stochastic process $A(t)$ is then a family of random variables indexed by $t$. Thus there must exist a family of probability densities $p(A,t)$ that allows to compute all moments $⟨A^n(t)⟩$ and in general any average $⟨f\left(A(t)\right)⟩$. However, $P(A,t)$ is not enough to characterize the stochastic process, because in general the variables $A(t)$ at different times are not independent. Thus $p(A,t)$ is actually a marginal distribution of some more complicated multivariate probability density.

It is natural to imagine a functional $P\left[A(t)\right]$ that gives the joint probability for all random variables $A(t)$. However, an infinite-dimensional probability distribution is a wild beast to ride, and we shall content ourselves with the (infinite) set of $n$-variable joint distributions

$$P_n(A_1,t_1,A_2,t_2,\mathrm{\dots },A_n,t_n).$$ (7)

Most stochastic processes can be completely specified by the set of all joint probabilities of the form (7) for all $n$ and all possible choices of $t_1,\mathrm{\dots },t_n$. The sense of “most” is highly technical [2], but should include all processes of interest to physicists. Here we shall need only the set corresponding to $n=2$ (which trivially gives also the set for $n=1$ marginalizing on the second variable).

The notion of stationary processes is an important one, connected with the thermodynamic idea of equilibrium. Completely stationary processes are those for which all the joint probability distributions that define the process are translation invariant, that is

$$P_n(A_1,t_1,\mathrm{\dots },A_n,t_n)=P(A_1,t_1+s,\mathrm{\dots },A_n,t_n+s),\forall s,n,t_i.$$ (8)

A less restrictive notion is that of stationary processes up order $M$. This is defined by the requirement that all the joint moments up to order $M$ exist and are time-translation invariant:

$$⟨A^{m_1}(t_1)A^{m_2}(t_2)\mathrm{\dots }{A}^{m_n}(t_n)⟩=⟨A^{m_1}(t_1+s)A^{m_2}(t_2+s)\mathrm{\dots }{A}^{m_n}(t_n+s)⟩$$ (9)

for all $s$, $n$, $\{t_1,\mathrm{\dots },t_n\}$ and $\{m_1,\mathrm{\dots },m_n\}$ such that $m_1+m_2+\mathrm{\dots }+m_n\le M$. This is less restrictive not only because of the bound on the number of random variables considered, but also because invariance is imposed only on the moments, and not on the joint distributions themselves.

In completely stationary processes the time origin is irrelevant: no matter when one performs an experiment, the statistical properties of the observed signal are always the same. In particular, this implies that the average $⟨A(t)⟩$ is constant in time (but does not mean that time correlations are trivial). This is the situation one expects when observing a system in thermodynamic equilibrium. For processes stationary up to order $M$, time origin is irrelevant for moments up to order $M$. For a stationary process up to first order, one can assert that $⟨A(t)⟩$ is a constant independent of $t$, and for a process stationary up to second order one can in addition assert that $⟨A^2(t)⟩$ (and hence the variance) is time-independent, and that $⟨A(t_1)A(t_2)⟩=⟨A(0)A(t_2-t_1)⟩$ (using $s=-t_1$), i.e. that all second-order moments are function only of the time difference.

2.3 Time correlation functions

The time correlation function is the correlation of the random variables $A(t_0)$ and $A(t_0+t)$, i.e. the second order moment

$$C(t_0,t)=⟨A^{\ast }(t_0)A(t_0+t)⟩=\int 𝑑A_1𝑑A_{2}P(A_1,t_0,A_2,t_0+t)A_1^{\ast }{A}_2,$$ (10)

where the star stands for complex conjugate. The time difference $t$ is sometimes called lag, and the function is called self correlation or autocorrelation to emphasize the fact that it is the correlation of the same observable quantity measured at different times. Note however that from the point of view of probability theory $A(t_0)$ and $A(t_0+t)$ are two different (though usually not independent) random variables. One can also define cross-correlations of two observables at different times:

$$C_{AB}(t_0,t)=⟨A^{\ast }(t_0)B(t_0+t)⟩.$$ (11)

The star again indicates complex conjugate. From now on however we shall restrict ourselves to real quantities and omit it in the following equations.

One can also define a correlation function of the fluctations $\delta A(t)=A(t)-⟨A(t)⟩$. This is called the connected time correlation function,

$$C_c(t_0,t)=⟨\delta A(t_0)\delta A(t_0+t)⟩=C(t_0,t)-⟨A(t_0)⟩⟨A(t_0+t)⟩.$$ (12)

Remembering (5) it is clear that this function is zero when the variables are uncorrelated (which we expect to happen for $t\to \mathrm{\infty }$). For this reason this function is often more useful than the correlation (10), which for uncorrelated variables takes the value $⟨A(t_0)⟩⟨A(t_0+t)⟩$.

The names employed here are usual in the physics literature (e.g. [3, 4]). In the mathematical statistics literature, the connected correlation function is called autocovariance (in fact it is just the covariance of $A(t_0)$ and $A(t_0+t)$), while the name autocorrelation is reserved for the normalized connected correlation defined below (15).

In what follows, unless otherwise stated, we will assume that we are dealing with processes stationary at least up to second order, so that the time correlation is a function only of the lag $t$. In this case the difference between the connected and nonconnected correlation functions is a constant,

$$C_c(t)=⟨A(0)A(t)⟩-{⟨A⟩}^2=C(t)-{⟨A⟩}^2,\text{(stationary process).}$$ (13)

Finally, let us define a normalized connected correlation so that its absolute value is bounded by 1, in analogy with the Pearson coefficient:

$\rho (t_0,t)$ $={\displaystyle \frac{C_c(t_0,t)}{\sqrt{C_c(t_0,0)C_c(t_0+t,0)}}},$ (14)

$\rho (t)$ $={\displaystyle \frac{C_c(t)}{C_c(0)}},$ (stationary). (15)

2.3.1 Properties of the time correlation function

From the definition of the correlation function, it is easy to see that in the stationary case the following hold:

1. $C_c(0)={\text{Var}}_A$.

2. $|C_c(t)|\le C_c(0)$ for all $t$.

3. If $A(t)$ is real-valued, the correlation is even: $C(t)=C(-t)$.

Also, for sufficiently long lags one expects that the variables become independent and correlation is lost, so that $C(t\to \mathrm{\infty })\to {⟨A⟩}^2$. Thus the connected correlation should tend to zero, and it will typically have a Fourier transform (see § 2.4)

2.3.2 Example

Let us conclude this section of definitions with an example. Fig. 1 shows on the left two synthetic stochastic signals, generated with the random process described in §4.2.2. On the right there are the corresponding connected correlations. The two signals look different, and this difference is captured by $C_c(t)$. We can say that fluctuations for the lower signal are more persistent: when some value of the signal is reached, it tends to stay at similar values for longer, compared to the other signal. This is reflected in the slower decay of $C_c(t)$.

Figure 1: Two stochastic signals (left) and their respective connected time correlations (right). Correlation times are $\tau \approx 4.5$ (top), $\tau \approx 100$ (bottom). The signals were generated through Eq. (44) with $\mu =0$, ${\sigma }^2=2$, $N={10}^5$, and $w=0.8$ (top), $w=0.99$ (bottom).

2.4 Fourier transforms of time correlations

Time correlation functions are often studied in frequency space, either because they are obtained experimentally in the frequency domain (e.g. in techniques like dielectric spectroscopy), because the Fourier transforms are easier to compute or handle analytically, because they provide an easier or alternative interpretation of the fluctuating process, or for other practical or theoretical reasons. Although the substance is the same, the precise definitions used can vary. One must pay attention to i) the convention used for the Fourier transform pairs and ii) the exact variables that are transformed.

Here we define the Fourier transform pairs as

$$\tilde{f}(\omega )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}𝑑tf(t)e^{i\omega t},f(t)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{d\omega }{2\pi }\tilde{f}(\omega )e^{-i\omega t},$$ (16)

but some authors choose the opposite sign for the forward transform and/or a different placement for the $1/2\pi$ factor (sometimes splitting it between the forward and inverse transforms). Depending on the convention a factor $2\pi$ can appear or disappear in some relations, like (20) below.

As for the second point above, time correlations are defined in terms of two times, $t_1$ and $t_2$ (which we chose to write as $t_0$ and $t_0+t$). One can transform one or both of $t_1$ and $t_2$ or $t_0$ and $t$. Let us consider two different choices, useful in different circumstances. First take the connected time correlation (12) and do a Fourier transform on $t$:

$$C_c(t_0,\omega )=\int 𝑑t{e}^{i\omega t}{C}_c(t_0,t_0+t)=e^{-i\omega t_0}⟨\delta A(t_0)\delta \tilde{A}(\omega )⟩,$$ (17)

where $\delta \tilde{A}(\omega )$ stands for the Fourier transform of $\delta A(t)$. This definition is convenient when there is an explicit dependence on $t_0$ but the evolution with $t_0$ is slow, as in the physical aging of glassy systems (see e.g. [5]); one then studies a time-dependent spectrum. If on the other hand $C_c$ is stationary, it is more convenient to write $t_1=t_0$, $t_2=t_0+t$ and do a double transform in $t_1$, $t_2$:

$$\begin{array}{c}{\tilde{C}}_c({\omega }_1,{\omega }_2)=\int 𝑑t_1𝑑t_2{e}^{i{\omega }_1{t}_1}{e}^{i{\omega }_2{t}_2}C(t_1,t_2-t_1)=\hfill \\ \hfill \int 𝑑t_1𝑑t_2{e}^{i{\omega }_1{t}_1}{e}^{i{\omega }_2{t}_2}⟨\delta A(t_1)\delta A(t_2)⟩=⟨\delta \tilde{A}({\omega }_1)\delta \tilde{A}({\omega }_2)⟩=\\ \hfill \int 𝑑t_1𝑑t_2{e}^{i({\omega }_1+{\omega }_2)t_1}{e}^{i{\omega }_2(t_2-t_1)}{C}_c(t_2-t_1)=(2\pi )\delta ({\omega }_1+{\omega }_2){\tilde{C}}_c({\omega }_2),\end{array}$$ (18)

where ${\tilde{C}}_c(\omega )$ is the Fourier transform of the stationary connected correlation with respect to $t$ and we have used the integral representation of Dirac’s delta, $\delta (\omega -{\omega }^{\prime })={(2\pi )}^{-1}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{it(\omega -{\omega }^{\prime })}𝑑t$. As the above shows, the transform is zero unless ${\omega }_1=-{\omega }_2$. For this reason it is useful to define the reduced correlation,

$$C_c^R(\omega )=⟨\delta \tilde{A}(-\omega )\delta \tilde{A}(\omega )⟩=⟨\delta {\tilde{A}}^{\ast }(\omega )\delta \tilde{A}(\omega )⟩,$$ (19)

where the rightmost equality holds when $A$ is real.

The transform ${\tilde{C}}_c(\omega )$ is a well-defined function, because the connected correlation decays to zero (and usually fast enough). We can then say

$${\tilde{C}}_c(\omega )=\frac{1}{2\pi \delta (0)}{C}_c^R(\omega ).$$ (20)

This relation furnishes a fast way to compute $C_c(t)$ in practice (§4.3). But let us comment on the at first sight baffling infinite factor relating the reduced correlation to ${\tilde{C}}_c(\omega )$. This originates in the fact that $\delta \tilde{A}(\omega )$ cannot exist as an ordinary function, since we have assumed that $A(t)$ is stationary. This implies that the signal has an infinite duration, and that, its statistical properties being always the same, it cannot decay to zero for long times as required for its Fourier transform to exist. The Dirac delta can be understood as originating from a limiting process where one considers signals of duration $T$ (with suitable, usually periodic, boundary conditions) and then takes $T\to \mathrm{\infty }$. Then $2\pi \delta (\omega =0)={\int 𝑑t{e}^{it\omega }|}_{\omega =0}=\int 𝑑t=T$. These considerations can be made rigorous by defining the signal’s power spectrum,

$$h(\omega )=\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}⟨A_T(\omega )A_T^{\ast }(\omega )⟩,A_T(\omega )={\int }_{-T/2}^{T/2}A(t)e^{i\omega t}𝑑t,$$ (21)

and then proving that $h(\omega )={\tilde{C}}_c(\omega )$ [2, §4.7, 4.8].

3 Correlation time

The connected correlation function measures how correlation is gradually lost as time elapses. One often seeks for a summary of this detailed information, in the form of a time scale that measures the interval it takes for significant decorrelation to happen: this is the correlation time $\tau$. While some definitions imply that after a correlation time the connected correlation has descended to a prescribed level, this quantity is usually better understood as a scale, which precise meaning depends on the details of the shape of the correlation function. It is most useful to compare correlation functions of similar shape, which change as some environmental conditions varies (e.g. when one studies a given system at a set of different temperatures). You should also keep in mind that correlation functions can be quite complicated, and it may be appropriate to describe them using more than one scale (e.g. they could be a superposition of exponential decays). Thus when one speaks of the correlation time one means (or should mean) the largest of them. For the purpose of comparing different decays to see which one is slower, it is less relevant which of the times one chooses as description of the decay as long as a) one chooses the same scale in all cases and b) one knows that all of them scale in the same way when the control variable is changed. We now discuss several possible definitions of $\tau$.

3.1 Correlation threshold

A simple, “quick and dirty”, way to define a correlation time is to choose it as the time it takes for the (normalized, connected) correlation to drop below a prescribed threshold $ϵ$:

$$\rho (t=\tau )=ϵ.$$ (22)

This definition is that is easy to apply to experimental data, and can be useful to compare correlation functions as long as they have the same shape. But it can be ambiguous if the decay displays oscillations, or if many time scales are present, and if it is inadequate when applied to functions of different shapes (see Fig. 2) or to power laws.

Figure 2: Different possible shapes for the decay of the time correlation. Top left: simple exponential. Top right: simple exponential (black curve) and double exponential (red curve) decay. In this case, the threshold criterion (here $ϵ=0.1$) labels the red curve as the fastest, but it clearly has a longer tail. Bottom left: stretched exponential. The average correlation time (26) is $⟨\tau ⟩\approx 22.7$ (black curve), $⟨\tau ⟩\approx 92.6$ (red curve). Bottom right: exponentially damped harmonic oscillations.

3.2 Fit parameter

Sometimes an analytical expression can be fit to the correlation function, and a time scale extracted from the fit parameters. For example, if one can fit an exponential decay,

$$C_c(t)=A{e}^{-t/\tau },$$ (23)

then the fitted value of the time scale $\tau$ can be used as correlation time (which in this particular case coincides with the threshold definition using $ϵ=1/e$). However, real-life time correlations are usually more complicated than a simple exponential. One can perhaps find more complicated functions that fit the correlation, but be wary of the proliferation of parameters (as when fitting the sum of two, three, $n$ exponentials).

You may may get away with fitting the only the last part of the decay with (23), if you can fit a sizeable part of the “tail”. The rationale is that the first part of the decay is dominated by fast microscopic processes one is (often) not interested in, so that by fitting the tail one obtains a good estimate of the correlation time of the slowest process. The problem with this strategy is how to decide when the tail starts; if the $\tau$ obtained this way is too sensitive to how this choice is made then it is probably no good.

A function widely used to describe non-exponential decays, employing only three parameters, is the Kolrausch-William-Watts or stretched exponential function,

$$C_c(t)=A{e}^{-{(t/\tau )}^{\beta }}.$$ (24)

Here $\tau$ is a time scale and the stretching exponent $\beta$ controls the shape ( gives a stretched exponential, i.e. a function with a longer tail compared to a simple exponential of the same $\tau$, while $\beta >1$ produces instead a compressed exponential). However comparing the $\tau$ of two functions with different $\beta$ can be misleading (see Fig. 2). A better description of the decay using a single number can be achieved by considering the correlation as a superposition of exponential processes,

$$\rho (t)={\int }_0^{\mathrm{\infty }}w(\tau )e^{-t/\tau }𝑑\tau ,$$ (25)

which defines the correlation time distribution function $w(\tau )$ as essentially the inverse Laplace transform of $\rho (t)$ [8]. Then the average correlation time is $⟨\tau ⟩=\int \tau w(\tau )𝑑\tau$. For the stretched exponential one finds [8]

$$⟨\tau ⟩=\frac{\tau }{\beta }\mathrm{\Gamma }\left(\frac{1}{\beta }\right),$$ (26)

where $\mathrm{\Gamma }(x)$ is Euler’s gamma function.

3.3 Integral time

A quite general way to define a correlation time is from the integral of the normalized connected correlation,

$${\tau }_{\text{int}}={\int }_0^{\mathrm{\infty }}\rho (t)𝑑t.$$ (27)

Clearly for a pure exponential decay ${\tau }_{\text{int}}=\tau$. In general, if $\rho (t)=f(t/\tau )$ then ${\tau }_{\text{int}}=\text{const}\tau$. The integral time is related to the variance of the estimate of the mean of the signal (see §4.2 below and [9, §2]). With some care, it can be computed from experimental or simulation data, avoiding the difficulties encountered when using thresholds or fitting functions. The procedure explained next is expected to work well if long enough sequences are at disposal and the decay does not display strong oscillations or anticorrelation [9].

If ${\widehat{C}}_k^{(c)}\approx C_c(k\mathrm{\Delta }t)$ is the estimate of the stationary connected correlation (obtained as explained in §4.2), then the integral can be approximated by a discrete sum, but the sum cannot run over all available values $k=0,\mathrm{\dots },N-1$, because as discussed in §4.2 the variance of ${\widehat{C}}_k^{(c)}$ for $k$ near $N-1$ is large, so that the sum ${\sum }_{k=0}^{N-1}{\widehat{C}}_k^{(c)}$ is dominated by statistical noise (more precisely, its variance does not go to zero for $N\to \mathrm{\infty }$). A way around this difficulty is [9, §3] to cut-off the integral at a finite time $t_c=c\mathrm{\Delta }t$ such that $c\ll N$ but the correlation is already small at $t=t_c$ (implying that $t_c$ is larger than a few times $\tau$). Thus ${\tau }_{\text{int}}$ is defined self-consistently as

$${\tau }_{\text{int}}={\int }_0^{\alpha {\tau }_{\text{int}}}\rho (t)𝑑t,$$ (28)

where $\alpha$ should be chosen larger than about $5$, and within a region of values such that ${\tau }_{\text{int}}(\alpha )$ is approximately independent of $\alpha$. Longer tails will require larger values of $\alpha$; we have found adequate values to be as large as 20 in some cases. To solve (28) you can compute $\tau (M)={\sum }_k^M{\widehat{C}}_k^{(c)}/{\widehat{C}}_0^{(c)}$ starting with $M=1$ and increasing $M$ until $\alpha \tau (M)>M$.

3.4 Correlation time from spectral content

Another useful definition of correlation time can be obtained from the Fourier representation of the normalized correlation,

$$\tilde{\rho }(\omega )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}𝑑t\rho (t)e^{i\omega t}.$$ (29)

Because $\rho (t)$ is normalized, ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{d\omega }{2\pi }\tilde{\rho }(\omega )=1$, a characteristic frequency ${\omega }_0$ (and a characteristic time ${\tau }_0=1/{\omega }_0$) can be defined such that $\tilde{\rho }(\omega )$ for $\omega \in [-{\omega }_0,{\omega }_0]$ holds half of the spectrum [10], i.e.

$${\int }_{-{\omega }_0}^{{\omega }_0}\frac{d\omega }{2\pi }\tilde{\rho }(\omega )=\frac{1}{2}.$$ (30)

This this definition of can be expressed directly in the time domain writing

$${\int }_{-{\omega }_0}^{{\omega }_0}\frac{d\omega }{2\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}𝑑t\rho (t)e^{i\omega t}=2{\int }_0^{\mathrm{\infty }}𝑑t\rho (t){\int }_{-{\omega }_0}^{{\omega }_0}\frac{d\omega }{2\pi }{e}^{i\omega t}=\frac{2}{\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}𝑑t\rho (t)\frac{\sin {\omega }_{0}t}{t},$$ (31)

where we have used the fact that $\rho (t)$ is even. Then the correlation time is defined by

$${\int }_0^{\mathrm{\infty }}\frac{dt}{t}\rho (t)\sin \left(\frac{t}{{\tau }_0}\right)=\frac{\pi }{4}.$$ (32)

It can be seen that if $\rho (t)=f(t/\tau )$, then ${\tau }_0$ is proportional to $\tau$ (it suffices to change the integration variable to $u=t/\tau$ in the integral above).

An advantage of this definition is that it copes well with the case when inertial effects are important and manifest in (damped) oscillations of the correlation function (see Fig. 2). In particular, for harmonic oscillations of frequency $\nu$, ${\tau }_0=1/(2\pi \nu )$ while ${\tau }_{\text{int}}$ is undefined.

4 Computing time correlations from experimental data

In this section we examine in detail how to compute in practice the time correlation function of a signal recorded in an experiment or produced in numerical simulation. Up till know we have discussed the theoretical definitions of correlation functions, which are given in terms of averages over some probability distribution, or ensemble. However, when dealing with experimental or simulation data we do not have direct access to the probability distribution, but only to a set of samples, i.e. results from experiments, distributed randomly according to an unknown distribution. We must try to compute the averages we want (the correlation functions), as accurately as possible, using these samples. This is what the field of statistics is about: building estimators that allow to compute the quantities of interest as best as possible from the available samples.

We assume that the experiment records a scalar signal with a uniform sampling interval $\mathrm{\Delta }t$, so that we are handled $N$ real-valued and time-ordered values forming a sequence $a_i$, with $i=1,\mathrm{\dots },N$. It is understood that if the data are digitally sampled from a continuous time process, proper filtering has been applied. In what follows we shall measure time in units of the sampling interval, so that in the formulas and algorithms below we shall make no use of $\mathrm{\Delta }t$. To recover the original time units one must simply remember that $a_i=A(t_i)$ with $t_i=t_0+(i-1)\mathrm{\Delta }t$. For the stationary case we shall write $C_k=C(t_k)$ where $t_k$ is the time difference, $t_k=k\mathrm{\Delta }t$ and $k=0,\mathrm{\dots },N-1$, and in the non-stationary case $C_{i,k}=C(t_i,t_k)$.

4.1 Stationary vs. non-stationary signals

It is clearly hopeless to attempt to estimate an ensemble average unless it is possible to obtain many samples under the same conditions (i.e., if one is throwing dice, one should throw many times the same dice). This implies there is a huge difference in how one can treat a sample from a stationary process vs. a sample from a non-stationary one. If the process is stationary, we can essentially consider the samples $a_i$ of one sequence as different repetitions of the same experiment. Estimation then basically consists in replacing ensemble averages with time averages, e.g. one estimates the mean $⟨A⟩$ as $\overline{a}$, with

$$\overline{a}=\frac{1}{N}\sum _{i=0}^{N-1}{a}_i,$$ (33)

and the stationary correlation with formula (37) below. There is more to be said (see §4.2 below), but in essence the situation is this.

On the other hand, if the process is not stationary it means that the samples of a single sequence are different experiments, in the sense that the conditions have changed from one sample to another (i.e. the system has evolved in some nontrivial way). In this case the correlation depends on two times (not just on their difference) and the mean itself can depend on time. The only way to estimate the mean or the correlation function is to obtain many sequences $a_i^{(n)}$, $n=1,\mathrm{\dots },M$ reinitializing the system to the same macroscopic conditions each time (in a simulation, one can for example restart the simulation with the same parameters but changing the random number seed). Time averaging is no good in this case, instead the ensemble average is replaced by average across the different sequences, i.e. the (time-dependent) mean is estimated as

$$⟨A(t_i)⟩\approx \overline{a_i}=\frac{1}{M}\sum _{n=1}^M{a}_i^{(n)},$$ (34)

and the time correlation as

$$C_c(t_i,t_k)\approx {\widehat{C}}_{i,k}^{(c)}=\frac{1}{M}\sum _n^M\delta a_i^{(n)}\delta a_{i+k}^{(n)},\delta a_i^{(n)}=a_i^{(n)}-\overline{a_i},$$ (35)

where the hat distinguishes an estimate from the actual quantity. These estimators have the desirable property that they “get better and better” when the number of samples $M$ grows, if the samples are independent and the ensemble has finite variance. More precisely, ${\widehat{C}}_{i,k}^{(c)}\to C(t_i,t_i)$ as $M\to \mathrm{\infty }$ (and similarly for $\overline{a}$). This property is called consistency, (see [2, §5.2]) and is guaranteed because $⟨{\widehat{C}}_{i,k}^{(c)}⟩=C(t_i,t_i)$ and $⟨\overline{a_i}⟩=⟨A(t_i)⟩$ (i.e. the estimators are unbiased) and their variance vanishes for $M\to \mathrm{\infty }$ (because the samples are independent).

If one has several sequences sampled from a stationary system, it is possible to combine the two averages: one first computes the stationary correlation estimate (37) for each sequence, and then averages the different correlation estimates (over $n$ at fixed lag $k$). It is clearly wrong to average the sequences themselves before computing the correlation, as this will tend to approach the (stationary) ensemble mean $⟨A⟩$ for all times, destroying the dynamical correlations.

Before turning to estimation in the stationary case, let us stress that it is always sensible to check whether the sequence is actually stationary. A first check on the mean can be done dividing the sequence in several sections and computing the average of each section, then looking for a possible systematic trend. If this check passes, then one should compute the time correlation of each section and then checking that all of them show the same decay (using fewer sections than in the first check, as longer sections will be needed to obtain meaningful estimates of the correlation function). It is important to note that this second check is necessary even if the first one passes; as we have noted above it is possible for the mean to be time-independent while the time correlation depends on two times (i.e. the stochastic process is only stationary to first order).

4.2 Estimation of the time correlation of a stationary process

If our samples $a_i$ are form a stationary signal, we build our estimators using time averages in lieu of ensemble averages, in particular the estimator for $⟨A⟩$ is (33). As we said above, the idea is that we regard the samples as different realizations of the same experiment. However, they are not independent realizations, but correlated realizations. The estimate (33) is still good in the sense that it is consistent, but it has higher error than the equivalent estimate built with independent samples, because it has a higher variance. The variance of the estimate with correlated samples is [9]

$${\text{Var}}_{\overline{a}}\approx \frac{2{\tau }_{\text{int}}}{N}\left[⟨A^2⟩-{⟨A⟩}^2\right],$$ (36)

i.e. $2{\tau }_{\text{int}}$ times larger than the variance in the independent sample case, where ${\tau }_{\text{int}}$ is the integral correlation time (27). In this sense $N/2{\tau }_{\text{int}}$ can be thought of as the number of “effectively independent” samples.

To estimate the connected correlation function we use

$${\widehat{C}}_k^{(c)}=\frac{1}{N-k}\sum _{j=1}^{N-k}\delta a_j\delta a_{j+k},\delta a_j=a_j-\overline{a},$$ (37)

where $\overline{a}$ is usually the estimate (33), although the true sample mean $⟨A⟩$ can be used if known. If the true mean is used, the estimator (37) is unbiased, otherwise it is asymptotically unbiased, i.e. the bias tends to zero for $N\to \mathrm{\infty }$, provided the Fourier transform of $C_c(t)$ exists. More precisely, $⟨{\widehat{C}}_{c,k}⟩=C_c(t_k)-\alpha /N$, where $\alpha =2\pi {\text{Var}}_a{\tilde{C}}_c(\omega =0)$. The variance of the estimate ${\widehat{C}}_{c,k}$ is $O(1/(N-k))$ [2, §5.3]. This is sufficient to show that, at fixed $k$, the estimator is consistent, i.e. ${\widehat{C}}_k^{(c)}\to C_c(t_k)$ for $N\to \mathrm{\infty }$. However the variance increases for increasing $k$, and thus the tail of ${\widehat{C}}_{c,k}$ is considerably noisy. In practice, for $k$ near $N$ the estimate is mostly useless, and the usual rule of thumb is to use ${\widehat{C}}_{c,k}$ only for $k\le N/2$.

It is natural to propose an estimate for the nonconnected correlation function doing a time average, analagous to (37),

$${\widehat{C}}_k=\frac{1}{N-k}\sum _{j=1}^{N-k}{a}_j{a}_{j+k}.$$ (38)

However, although ${\widehat{C}}_k$ is unbiased, it can have a large variance when the signal has a mean value larger than the typical fluctuation (see §4.2.1). In general, the connected estimator should be used (an exception may be when an experiment can only furnish many short independent samples of the signal, see the example below in §4.2.2).

Another asymptotically unbiased estimator of the connected correlation can be obtained by using $1/N$ instead of $1/(N-i)$ as a the prefactor of the sum in (37). Calling this estimator $C_{c,k}^{\ast }$, it holds that $⟨C_{c,k}^{\ast }⟩=C_c(t_k)-\alpha /N-kC(t_k)/N-\alpha k/N^2$, where $\alpha$ is defined as before, and again $\alpha =0$ if the exact sample mean is used. This has the unpleasant property that the bias depends on $k$, while the bias of ${\widehat{C}}_{c,k}$ is independent of its argument, and smaller in magnitude. The advantage of $C_{c,k}^{\ast }$ is that its variance is $O(1/N)$ independent of $k$, thus it has a less noisy tail. Many authors prefer $C_{c,k}^{\ast }$ due to its smaller variance and to the fact that it strictly preserves properties of the correlation (like property 2 of § 2.3.1), which may not hold for ${\widehat{C}}_{c,k}$. Here we stick to ${\widehat{C}}_{c,k}$, as usual in the physics literature (e.g. [9, 12, 13]), so that we avoid worrying about possible distortions of the shape at small $k$. In practice however, it seems that as long as $N$ is greater than $\sim 10\tau$ (a necessary requirement in any case, see below), and for $k\le N/2$, there is little numerical difference between the two estimates.

4.2.1 Two properties of the estimator

We must mention two important features of the estimator that are highly relevant when attempting to compute numerically the time correlation. The first is that the difference between the non-connected and connected estimators is not ${\overline{a}}^2$, but

$${\widehat{C}}_k-{\widehat{C}}_{c,k}=\overline{a}\left[\frac{1}{N-k}\sum _j^{N-k}{a}_j+\frac{1}{N-k}\sum _j^{N-k}{a}_{j+k}-\overline{a}\right],$$ (39)

as it is not difficult to compute. The difference does tend to ${\overline{a}}^2$ for $N\to \mathrm{\infty }$, but in a finite sample fluctuations are important, especially at large $k$. Fluctuations are additionally amplified by a factor $\overline{a}$, so that when the signal mean is large with respect to its variance, the estimate ${\widehat{C}}_k$ is completely washed out by statistical fluctuations. In practice, this means that while it is true that $C(t)$ and $C_c(t)$ differ by a constant term (${⟨A⟩}^2$), in general it is a bad idea to compute ${\widehat{C}}_k$ and subtract ${\overline{a}}^2$ to obtain an estimate of the connected correlation. Instead, one computes the connected correlation directly (by computing first an estimate of the mean and using (37)) and then, if needed, adds ${\overline{a}}^2$ to obtain an estimation of $C(t)$.

Another important fact is that the estimator of the connected correlation will always have zero, whatever the length of the sample used, even if $N\ll \tau$. To see this, consider the quantity

$$B_i=(N-i){\widehat{C}}_{c,i}=\sum _{j=1}^{N-i}\delta a_j\delta a_{j+i},$$ (40)

and compute the sum

$$\sum _{i=0}^{N-1}{B}_i=\sum _{i=0}^{N-1}\sum _{j=1}^{N-i}\delta a_j\delta a_{j+i}=\sum _{i=0}^{N-1}\sum _{j=1}^N\sum _{k=1}^N\delta a_j\delta a_k{\delta }_{k,j+i}.$$ (41)

But ${\sum }_{i=0}^{N-1}{\delta }_{k,j+i}$ equals 1 if $k\ge j$ and 0, so

$$\begin{array}{c}\sum _{i=0}^{N-1}{B}_i=\sum _{j=1}^N\sum _{k=j}^N\delta a_j\delta a_k=\frac{1}{2}\sum _{j\ne k}^N\delta a_j\delta a_k+\sum _{j=1}^N{(\delta a_j)}^2=\hfill \\ \hfill \frac{1}{2}{\left[\sum _{j=1}^N\delta a_j\right]}^2+\frac{1}{2}\sum _{j=1}^N{(\delta a_j)}^2=\frac{1}{2}\sum _{j=1}^N{(\delta a_j)}^2>0,\end{array}$$ (42)

where the last equality follows because ${\sum }_j\delta a_j=0$. Now we can easily do the same sum starting from $i=1$:

(43)

This shows that at least some of the $B_i$ must be negative. But since $B_0>0$, the conclusion is that $B_k$, and hence ${\widehat{C}}_{c,k}$, which differs from it by a positive factor, must change sign at least once for $k\ge 1$.

The practical consequence of this is that when $N$ is of the order of $\tau$, or smaller, the estimate ${\widehat{C}}_{c,k}$ will suffer from strong finite-size effects, and its shape will be quite different from the actual $C_c(t)$. In particular, since ${\widehat{C}}_{c,k}$ will intersect the $x$-axis, it can mislead you into thinking that $\tau$ is smaller than $N$ when in fact it is several times larger. Be suspicious if ${\widehat{C}}_{c,k}$ changes sign once and stays very anticorrelated. Anticorrelation may be a feature of the actual $C_c(t)$, but if the sample is long enough, the estimate should decay until correlation is lost (noisily oscillating around zero). One must always perform tests estimating the correlation with different values of $N$: if the shape of the correlation at short times depends on $N$, then $N$ must be increased until one finds that estimates for different values of $N$ coincide for lags up to a few times the correlation time (see example below).

Once a sample-size-independent estimate has been obtained, the correlation time can be estimated (§3), and it must be checked that self-consistently $N$ is several times larger than $\tau$.

4.2.2 Example

To illustrate the above considerations, we generate a correlated sequence from the recursion

$$a_i=w{a}_{i-1}+(1-w){\xi }_i,$$ (44)

where $w\in [0,1]$ is a parameter and the ${\xi }_i$ are uncorrelated Gaussian random variables with mean $\mu$ and variance ${\sigma }^2$. Assuming the signal is stationary it is not difficult to find

$$⟨a⟩=\mu ,⟨{(a-\mu )}^2⟩=\frac{1-w}{1+w}{\sigma }^2,C_k^{(c)}=⟨(a_i-\mu )(a_{i+k}-\mu )⟩={\sigma }^2{w}^k,$$ (45)

so that the correlation time is $\tau =-1/\log w$.

We used the above recursion to generate artificial sequences and computed their time correlation functions with the estimates discussed above. Fig. 3 shows the problem that can face the non-connected estimate. When the average of the signal is smaller than or of the order of the noise amplitude (as in the top panels), one can get away with using (38). However if $\mu \gg \sigma$, the non-connected estimate is swamped by noise, while the connected estimate is essentially unaffected (bottom panels). Hence, if one is considering only one sequence, one should always use the connected estimator. The situation might be different if many sequences, corresponding to different realizations of the same experiment, are available (see the example after next).

Figure 3: Connected vs. nonconnected estimate. The estimate of $C(t)$ (equation (38), left panels) is much worse that the estimates of $C_c(t)$ (equation (37), right panels). The (artificial) signal was generated with (44). Top panels: $\mu =1$; bottom panels: $\mu =100$. In both cases, ${\sigma }^2=1$, $w=0.95$ ($\tau \approx 20$) and length $N=5\cdot {10}^4$.

In Fig. 4 we see how using samples that are too short affects the correlation estimates. The same artificial signal was generated with different lengths. For the shorter lengths, it is seen that the correlation estimate crosses the $t$ axis (as we have shown it must) but does not show signs of losing correlation. One might hope that $N=1000\approx 10\tau$ is enough (the estimate starts to show a flat tail), but comparing to the result of doubling the length shows that it is still suffering from finite-length effects. It is seen that a sequence at least $20\tau$ to $50\tau$ long is needed to get the initial part of the normalized connected correlation more or less right, while a length of about 1000$\tau$ is necessary to obtain a good estimate up to $t\sim 5\tau$. The unnormalized estimator suffers still more from finite size due to the increased error in the determination of the variance (left panel).

Figure 4: Finite size effects. Estimates of the connected correlation (top) and normalized connected correlation (bottom) for sequence (44) of different lengths as indicated, together with the analytical result $C_c(t)={\sigma }^2{w}^t$. Parameters are $\mu =0$, ${\sigma }^2=1$, $w=0.99$ ($\tau \approx 100$).

If the experimental situation is such that it is impossible to obtain sequences much longer than the correlation time, one can get some information on the time correlation if it is possible to repeat the experiment so as to obtain several independent and statistically equivalent sample sequences. In Fig. 5 take several (say $M$) sequences with the same parameters as in the previous example, but quite short ($N=200\approx 2\tau$). As we have shown, it is not possible to obtain a moderately reasonable estimate of $C_c(t)$ using one such sequence (as is also clear from the $M=1$ case of Fig. 5, where the analytical result is plotted for comparison). However, the figure shows how one may benefit from the multiple repetitions of the (short) experiment by averaging together all the estimates. The averaged connected estimates are always far from the actual correlation function, even for $M=500$ (the case which contains in total ${10}^5$ points, which proved quite satisfactory in the previous example): this is consequence of the fact that all connected estimates must become negative. Instead, the averaged non-connected estimates approach quite closely the analytical result even though not reaching the decorrelated region. Although it is tricky to try to extract a correlation time from this estimate (one may fit the initial decay but it is not possible to know whether the final decay will follow this trend or whether some very slow tail is present), this procedure at least offers a way to obtain some dynamical information in the face of experimental limitations.

Figure 5: Effect of averaging the estimates of many short sequences. We show the connected (left) and non-connected (right) estimates of $M$ different samples of sequence (44) as indicated in the legend, with $N=200$, $\mu =0$, ${\sigma }^2=1$, $w=0.99$ ($\tau \approx 100$).

4.3 Algorithms to compute the estimators

The estimators can be computed numerically by straightforward implementation of equations (35) or (37), although in the stationary case it is much more efficient to compute the connected correlation through relation (21) using a fast Fourier transform (FFT) algorithm. Let us focus on the stationary case and examine in some detail these algorithms.

Algorithm 1 presents the direct method. It is straightforward to translate the pseudo-code to an actual language of your choice. Apart from some missing variable declarations, the only thing to consider is that it is probably not convenient (or even illegal in some languages, as in classic C) to return a large array, and it is better to define $C$ as an output argument, using a pointer or reference (as e.g. FORTRAN or C do by default) to avoid copying large blocks of data. The advantages of this algorithm are that it is self-contained and simple to understand and implement. Its main disadvantage is that, due to the double loop of lines 8–11, it runs in a time that grows as $N^2$. For $N$ up to about ${10}^5$ this algorithm is perfectly fine: a good implementation in a compiled language should should run in a few seconds in a modern computer. But this time grows quickly; in the author’s computer $N=5\cdot {10}^5$ takes 35 seconds, for $N={10}^6$ the time is two and a half minutes. In contrast, the algorithm with FFT takes 1 second for $N={10}^6$ and 11 seconds for $N={10}^7$.

Algorithm 1 Compute the connected correlation of sequence $a$ (of length $N$) using the direct $O(N^2)$ method. The connected correlation is returned in vector $C$.

1:function timecorr($a$,$N$)

2:     $\mu \leftarrow 0$ $▷$ Compute average

3:     for $i=1,\mathrm{\dots },N$ do

4:         $\mu \leftarrow \mu +a_i$      

5:     $\mu \leftarrow \mu /N$

6:     for $i=1,\mathrm{\dots },N$ do $▷$ Clear $C$ vector

7:         $C_i\leftarrow 0$      

8:

9:     for $i=1,\mathrm{\dots },N$ do $▷$ Correlation loop

10:         $d\leftarrow a_i-\mu$

11:         for $k=0,\mathrm{\dots },N-i$ do

12:              $C_{k+1}\leftarrow C_{k+1}+d\ast (a_{i+k}-\mu )$               

13:

14:     for $i=1,\mathrm{\dots },N$ do $▷$ Normalize and return

15:         $C_i\leftarrow C_i/(N-i-1)$      

16:     return C

So, if you need the correlation of really long sequences, the FFT-based algorithm, though more difficult to get running, pays off with huge savings in CPU time at essentially the same numerical precision. The idea of the algorithm is to compute the Fourier transform of the signal, use (20) to obtain the Fourier transform of the connected correlation, then transform back to obtain $C_c(t)$. This is faster than algorithm 1 because the clever FFT algorithm can compute the Fourier transform in a time that is $O(N\log N)$.

Actually, we need discrete versions of the Fourier transform formulas (as we remarked before, the Fourier transform of the continuous time signal does not exist). The discrete Fourier transform (DFT) and its inverse operation are defined [11, §12.1] as (it is convenient to let the subindex of $a_i$ run from $0$ to $N-1$ to write the following two equations),

$${\tilde{a}}_k=\sum _{j=0}^{N-1}{e}^{2\pi ijk/N}{a}_j,a_j=\frac{1}{N}\sum _{k=0}^{N-1}{e}^{-2\pi ijk/N}{\tilde{a}}_k,$$ (46)

where we note that the inverse DFT effectively extends the sequence periodically (with period $N$). The discrete version of (20) is [11, §13.2]

${\tilde{D}}_k$ $={|{\tilde{a}}_k|}^2,$ where $D_j$ $={\displaystyle \sum _{k=0}^{N-1}}{a}_k{a}_{k+j},$ (47)

and where the definition of $D_j$ makes use of the (assumed) periodicity of $a_i$. $D_j$ is almost our estimate (37): we only need to take care of the normalization and of the fact that due to the assumed periodicity of $a_i$ some past times are regarded as future, e.g. for $k=10$, in the sum there appear the terms $a_0{a}_{10}$ up to $a_{N-11}{a}_{N-1}$ (which are fine), but also $a_{N-10}{a}_0$ through $a_{N-1}{a}_9$, which we do not want included. This is fixed by padding the original signal with $N$ zeros at the end, i.e. setting $a_k=0$ for $k=N,\mathrm{\dots },2N-1$ and ignoring the values of $D_j$ for $j\ge N$.

In summary, to compute the connected correlation using FFT the steps are i) estimate the mean and substract from the $a_i$, ii) add $N$ zeros to the end of the sequence, iii) compute the DFT of the sequence, iv) compute the squared modulus of the transform, iv) compute the inverse DFT of the squared modulus, v) multiply by the $1/(N-i)$ prefactor. Pseudocode for this algorithm is presented as algorithm 2.

Algorithm 2 Compute the connected correlation of sequence $a$ (of length $N$) using a fast Fourier transform. This algorithm is $O(N\log N)$.

1:function timecorr($a$,$N$)

2:     $\mu \leftarrow 0$ $▷$ Compute average

3:     for $i=1,\mathrm{\dots },N$ do

4:         $\mu \leftarrow \mu +a_i$      

5:     $\mu \leftarrow \mu /N$

6:     for $i=0,\mathrm{\dots },N$ do $▷$ Substract the average from signal

7:         $a_i\leftarrow a_i-\mu$      

8:     for $i=N,\mathrm{\dots },2N$ do$▷$ Pad with 0s at the end

9:         ${\text{a}}_i\leftarrow 0$      

10:

11:     $b\leftarrow$FFT($a$,$2N$) $▷$ Compute the FFT of a as a vector of length $2N$

12:     for $i=1,\mathrm{\dots },2N$ do $▷$ Compute squared modulus of $b$

13:         $b_i\leftarrow {|b_i|}^2$ $▷$ Note that the Fourier transform is complex      

14:     $C\leftarrow$IFFT($b$,$2N$) $▷$ Inverse FFT

15:

16:     $C\leftarrow$resize($C$,$N$) $▷$ Discard the last $N$ elements of $C$

17:     for $i=1,\mathrm{\dots },N$ do $▷$ Normalize and return

18:         ${\text{C}}_i\leftarrow {\text{C}}_i/(N-i-1)$      

19:     return C

To translate this into an actual programming language the comments made for algorithm 1 apply, and in addition some extra work is needed for lines 10–13. First, one needs to choose an FFT routine. If you are curious about the FFT algorithm, you can read for example [11, Ch. 12] or [14], but writing an FFT routine is not easy, and implementing a state-of-the-art FFT is stuff for professionals. Excellent free-software implementations of the FFT can be found on Internet. FFTW [15], at http://www.fftw.org deserves mention as particularly efficient, although it is a large library and a bit complex to use. Note that some simpler implementations require that $N$ be a power of two, failing or using a slow $O(N^2)$ algorithm if the requiriment is not fulfilled. Also pay attention to i) the difference between “inverse” and “backward” DFTs (the latter lacks the $1/N$ factor), ii) how the routine expects the data to be placed in the input array, iii) how it is returned, and iv) whether the transform is done “in place” (i.e. overwriting the original data) or not. If the routine is a “complex FFT” it will expect complex input data (so that for real sequences you will have to set to zero the imaginary part of the $a_i$), while if it is a “real FFT” routine it will typically arrange (“pack”) the output data in some way, making use of the discrete equivalent of the $\tilde{A}(-\omega )=A^{\ast }(\omega )$ symmetry so as to return $N$ real numbers instead of $2N$ (the real and imaginary parts of the complete DFT). This affects the way you must compute the squared modulus (lines 11–12). For example, for the packing used by the FFTW real routines, lines 11–12 translate to (in C)

⬇

b[0]*=b[0];

for (int i=1; i<N; i++) {

b[i] = b[i]*b[i] + b[2*N-i]*b[2*N-i];

b[2*N-i] = 0;

}

b[N]*=b[N];

5 Conclusion

I have tried to convey the basic notions about time correlation functions as well as some practical advice on how to compute them from actual data. I hope this account will be useful for students and researchers finding themselves in need to compute time correlations.

On closing, I wish to thank the colleagues with whom I have worked over the years, and which have helped shape my understanding of time correlations through discussions on concepts and practicalities, in particular A. Cavagna, I. Giardina, V. Martín-Mayor, G. Parisi and the late P. Verrocchio.

Last but not least, I thank D. Chialvo for contributing the title, for organizing and inviting me to the Complexity Weekend, for agreeing that a tutorial of this sort should be written, and for editing and getting this volume to the press.