Consider $X$ to be a random variate with defined expectation $\mu$ and finite variance $\sigma^2$ . Consider a ton of points $X_1,X_2,..X_n$ generated via some process. The points need not be all independent of each other i.e there exists some auto-correlation amongst them.

We are interested in calculating the auto-correlation for some function $f$, which maps over these points, defined as $$ \rho_h(f)=R_h(f)/R_0(f)$$ where $$ R_h(f)=\frac{1}{n-h}\sum_{i=1}^{n-h}\left(f(x_i)-\langle f\rangle\right)(f(x_{i+h})-\langle f\rangle)\\ \langle f\rangle=\frac{1}{n}\sum_{i=1}^{n}f(x_i)$$ In particular consider $f(x)=x$. Then $$ \begin{align} \rho_h&=\frac{\frac{1}{n-h}\sum_{i=1}^{n-h}(x_i-\mu)(x_{i+h}-\mu)}{\frac{1}{n}\sum_{i=1}^{n}(x_i-\mu)^2}\\ &\\ &=\frac{\langle (X_i-\mu)(X_{i+h}-\mu)\rangle}{\sigma^2}\tag1 \end{align}$$ However the ACF given by mathematica,CorrelationFunction, is $$ \rho_h=\frac{\sum_{i=1}^{n-h}(x_i-\mu)(x_{i+h}-\mu)}{\sum_{i=1}^{n}(x_i-\mu)^2}\tag 2$$ 1. How to get the correct value (1) from (2)?
2. Is there some inbuilt option to get this correction instead of having to multiply (2) by $\frac{n}{n-h}$? The problem with this is I wish to calculate $\sum_h\rho_h$ using CorrelationFunction[..,{1,10^4}] and so explicit casewise correction is awkward.


Currently using

ACFTotal[data_, {a_, b_}] := 
 With[{n = Length@data}, 
  CorrelationFunction[data, {a, b}].Array[1.0/(1.0 - #/n) &, 
    b - a + 1, a]]


rn = 10^4;
rPts = RandomFunction[ARMAProcess[1, {-.7, .1}, {.8}, 1], {0, rn}];

AbsoluteTiming[ACFTotal[rpts, {1, rn - 1}]]
(*{0.0157095, -4.44417}*)

compared to @JimB's

AbsoluteTiming[Sum[CorrelationFunction[rpts, h] rn/(rn - h), {h, rn - 1}]]
(* {95.271, -4.44417}*)


Currently I am using
ACFTotal[data_, {a_, b_}] := With[{n = Length@data},CorrelationFunction[data,{a, b}].Array[1.0/(1.0 - #/n) &,b - a + 1, a]] but its slow

  • $\begingroup$ Please add in the definition of $\mu$. I'm assuming that you mean $\bar{x}$ (the sample mean). You seem to be mixing parameters with sample moments. Did you also mean $\hat{\sigma}$ (the sample variance)? Have you considered using CovarianceFunction rather than CorrelationFunction? Please consider changing the title. The definition isn't wrong. It just isn't the one you want. $\endgroup$
    – JimB
    Jan 23, 2020 at 15:55
  • $\begingroup$ @JImB $\mu$ is the expectation of the underlying probability distribution that all $X_i$ follow. Sample mean,as I recall, is the arithmetic mean, right? Thats not what $\mu$ is. Similarly for variance, its not the sample variance but the expectation of the deviation squared of every $X_i$ wrt its underlying probability distrib. $\endgroup$
    – lineage
    Jan 23, 2020 at 16:08
  • $\begingroup$ OK. So you do mean the parameters rather than the sample moments. So where do you plug in $\mu$ and $\sigma$ into your function ACFTotal? Yes, I'm easily confused. For example your two definitions of $\rho_h$ use different notations. And the CorrelationFunction uses $\bar{x}$ rather than $\mu$. $\endgroup$
    – JimB
    Jan 23, 2020 at 16:14
  • $\begingroup$ @JImB As for the definition for ACF--its autocovariance normalized wrt autocovarinace at lag 0. For eg. as said on para 1, pg 30, arxiv.org/abs/hep-lat/0702020 or as perhaps given in the wiki at en.wikipedia.org/wiki/Autocorrelation. In either case, the definition involves expectations, which necessarily require appropriate division by the number of quantities. Its in this sense that I consider the inbuilt definition to be wrong. However, pardon my ignorance towards the possibility of a completely plausible reason behind defining it that way in Mathematica. $\endgroup$
    – lineage
    Jan 23, 2020 at 16:17
  • $\begingroup$ @JImB see (1) where both $\mu,\sigma$ find use. Its just that when calculating for some $f(x)=x$,$\langle f\rangle \to \mu$ and $R_0(f)\to R_0(x)=\sigma^2$ $\endgroup$
    – lineage
    Jan 23, 2020 at 16:22

1 Answer 1


Update: My original answer which constructed the estimate of the sum of all of the correlations one-at-a-time was much, much slower than the method provided by the OP. (I should have followed the RTFM advice I keep giving others.) But below I give what I hope are relevant comments:

  1. Parameters and estimates should have separate notations such as using $\mu$ (as a parameter) and $\hat{\mu}$ as an estimator.

  2. There is another reason not to use the sum of the unadjusted serial correlations: that sum is always -1/2 no matter what the data looks like, no matter what the process was to generate the data.

  3. Multiplying each correlation by the associated $n/(n-h)$ might remove much of the bias but greatly increases the variability in the estimate of the sum of the correlations. There is undoubtedly a better estimator especially if some form of the generating process is known (ARMA, ARIMA, etc.).

  4. If you do know the form of the process, then you can use EstimatedProcess and then obtain the correlation coefficients for that estimated process. That approach will get you much better precision. In fact, the brute force approach of adding up the individual correlations (no matter how fast) will almost certainly be completely unacceptable in terms of precision.

Here's an approach that I hope is convincing about (2). The sum of the (unadjusted) correlations can be found in general for any set of $n$ data points:

TableForm[Table[{n, Simplify[
    Total[CorrelationFunction[Table[z[i], {i, n}], {1, n - 1}]],
    Assumptions -> Table[z[i] ∈ Reals, {i, n}]]}, {n, 4, 15}],
 TableHeadings -> {None, {"n", "Sum[\!\(\*OverscriptBox[\(ρ\), \(^\)]\)]"}}]

Sum of correlations for various sample sizes

Here is an example addressing (3): removing bias in the manner proposed greatly increases the variability of the sum of the estimated correlations because of the increased variability for the correlations of large lags (which have few samples to use).

rn = 100;
rPts = RandomFunction[ARProcess[{3/4}, 1], {0, rn}][[2, 1, 1]];
ListPlot[{Table[(3/4)^h, {h, 1, rn - 1}],
  CorrelationFunction[rPts, {1, rn - 1}],
    rPts, {1, rn - 1}] ((rn/(rn - #)) & /@ Range[rn - 1])},
 Joined -> True, PlotRange -> All,
 PlotStyle -> {{Thickness[0.01], Red}, {Thickness[0.01], Green}, Blue},
 PlotLegends -> {"True", "Unadjusted", "Adjusted"}]

Variability of estimates of correlation

For (4) here's an example showing how to use the form of the data generation process to obtain an estimate:

(* Generate data *)
rn = 10^4
rPts = RandomFunction[ARMAProcess[1, {-.7, .1}, {.8}, 1], {0, rn}][[2, 1, 1]];

(* Estimate parameters *)
ep = EstimatedProcess[rPts, ARMAProcess[a, {b, c}, {d}, e]]

(* Estimated total of correlation coefficients *)
Total[CorrelationFunction[ep, {1, rn - 1}]]

(* True total of correlation coefficients *)
Total[CorrelationFunction[ARMAProcess[1, {-.7, .1}, {.8}, 1], {1, rn - 1}]]
(* 0.125702 *)

(* Estimated total of correlation coefficients *)
Total[CorrelationFunction[ep, {1, rn - 1}]]
(* 0.127698 *)
  • $\begingroup$ the $\hat{\mu} \to \mu$ is an OK allowance as I am taking large n The factor,$n/(n-h)$ on the other hand, affects heavily $\Sigma_{h}\rho_h$ $\endgroup$
    – lineage
    Jan 23, 2020 at 16:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.