Wrong definition of auto-correlation function(ACF)

Question

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.

---

Update2

Currently using

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

test

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}*)

Update

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

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[\(ρ\), \(^\)]\)]"}}]

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;
SeedRandom[12345]
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}],
  CorrelationFunction[
    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"}]

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

(* Generate data *)
SeedRandom[123]
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 *)