Your ACF Can Be Right While Its Confidence Interval Is Wrong

The autocorrelation coefficients agreed down to the last stored bit. Their uncertainty did not.

In the statsmodels behavior covered by a contribution to improve missing-data handling, calling acf with missing="drop" on an array containing NaNs produced the same point estimates as removing the missing values by hand. On the affected code, the confidence intervals were narrower. The archived lag-one cases discussed below also returned larger Ljung-Box statistics. One function call was describing two sample sizes.

This is an awkward failure mode because the familiar diagnostic - the ACF plot - still looks plausible. If a test checks only the coefficients, it passes. The inconsistency appears in the quantities used to decide which spikes matter.

Here is the smallest version of the comparison.

Reproduce the invariant before inspecting the formula

import numpy as np
from statsmodels.tsa.stattools import acf

rng = np.random.default_rng(0)
kept = rng.standard_normal(200)
with_nans = np.concatenate([kept, np.full(200, np.nan)])
options = dict(nlags=10, qstat=True, alpha=0.05, fft=False)

drop = acf(with_nans, missing="drop", **options)
manual = acf(kept, missing="none", **options)

for label, index in [
    ("ACF", 0), ("interval", 1), ("Q statistic", 2), ("p-value", 3)
]:
    print(label, np.allclose(drop[index], manual[index]))

On the pre-fix checkout, the four lines read True, False, False, False. With the merged correction, all four read True. The comparison uses the documented meaning of missing="drop": remove missing observations, then estimate autocovariances while treating the survivors as contiguous.

This snippet uses the first draw from its seeded generator to expose tuple equality. The archived fixture behind the numerical examples used a separately frozen input. Any fixture observation below is labelled as such; its exact output should not be attributed to this snippet.

That last phrase matters. Manual deletion is a reference for this specific estimator; it is not a general recipe for irregular time series. If a missing reading represents elapsed time, joining its neighbors changes the time axis.

Figure 1. The point estimator and its downstream uncertainty calculations must receive the same effective sample size within this drop case.

The diagram names the two counts I will use throughout:

• n_raw: the length of the supplied array, including missing entries;
• n_eff: the number of non-missing observations that survives the chosen missing-data rule.

In the archived reproduction, n_raw = 400 and n_eff = 200. Before the patch, the ACF coefficients followed the effective series while the interval and Q-statistic paths still received 400.

Why the coefficients could still agree

acf delegates the autocovariance calculation to acovf. In the drop branch, those missing values are removed before the autocovariances are estimated. The point ACF therefore matched the manually compressed series.

The caller also retained x.shape[0] as its own observation count. That stale count was reused later for two calculations:

1. the Bartlett or non-Bartlett confidence interval;
2. q_stat(acf[1:], nobs=...).

So the point estimate and the uncertainty attached to it had silently diverged. A point-value assertion could not catch this branch; assertions on the returned interval and Q statistic could.

The source already had a comment asking whether nobs should shrink for missing="drop", and the tests contained a disabled FIXME for the corresponding Q-statistic case.

The confidence-band error is an exact square-root ratio

For a positive lag k, the Bartlett half-width implemented by acf can be written as

h_k(n) = z * sqrt((1 + 2 * sum_{j=1}^{k-1} r_j^2) / n)

Here, z is the normal quantile for the requested confidence level and r_j is the estimated ACF at lag j. Hold those coefficients and the quantile fixed. Everything above n is then a common factor, so the old-to-corrected half-width ratio is

h_k(n_raw) / h_k(n_eff) = sqrt(n_eff / n_raw)

This identity is exact for the implemented positive-lag Bartlett formula under those conditions. It also holds for the non-Bartlett branch, whose half-width is z / sqrt(n). It does not establish the empirical coverage of either interval.

At n_raw = 400 and n_eff = 200, the ratio is sqrt(200 / 400) = 0.7071. The affected result returned a lag-one half-width of 0.097998; the manually dropped reference returned 0.138590. I use “half-width” deliberately: it is the distance from the ACF estimate to one endpoint. The full interval width is twice that value.

Figure 2. In this seeded 50%-missing drop case, changing the count widens the interval while leaving every plotted ACF coefficient fixed.

The plot isolates the bug from the estimator. There is only one set of point markers because the old and corrected ACF arrays are bit-identical in the frozen experiment. The band changes because its denominator changes from 40 to 20.

The Q-statistic ratio is lag-dependent

The returned Q statistic needs a separate derivation. Statsmodels uses the Ljung-Box expression

Q_m(n) = n(n+2) * sum_{k=1}^m r_k^2 / (n-k)

The leading factor suggests that doubling n might double Q. The lag-specific denominators prevent that from being an identity. For fixed ACF coefficients, the exact cumulative old-to-corrected ratio is

R_m = [n_raw(n_raw+2) / n_eff(n_eff+2)]
      * [sum r_k^2/(n_raw-k) / sum r_k^2/(n_eff-k)]

For a single lag-k contribution, this simplifies to

n_raw(n_raw+2)(n_eff-k) / [n_eff(n_eff+2)(n_raw-k)]

Only when n is large relative to the lags does n_raw / n_eff become a useful large-n, small-lag approximation. In the archived 400-versus-200 lag-one fixture, the approximation is 2; the exact single-term formula gives 1.985111293. For a nonzero lag-one coefficient, that ratio follows from the two counts because the common r_1^2 cancels.

Figure 3. The confidence ratio and computed lag-one Q ratio are exact for each frozen case; only the n_raw / n_eff line is an approximation to the Q-statistic ratio.

The five points in Figure 3 come from a seeded AR(1) series with phi = 0.7, n_raw = 40, and nested missing indices at 0%, 10%, 25%, 50%, and 75%. Each fraction is a different sample, so its ACF coefficients can differ from the other fractions. Within any one fraction, the old and corrected calculations share that case’s coefficients; only their count changes.

In the archived 400-versus-200 lag-one fixture and every nonzero-missingness lag-one case plotted in Figure 3, the stale larger count produced the larger Q statistic. That direction is case-specific. The (n-k) terms can reverse the direction for some lag-and-count combinations. For a fixed chi-square reference distribution, the p-value falls only when Q rises; if Q falls, the p-value rises. None of these calculations measures a false-positive rate. Estimating one would require simulations under a null model across explicit missingness mechanisms.

The raw qstat=True result also has a narrower interpretation than “a universal test for autocorrelation.” It is a returned Ljung-Box portmanteau statistic. When testing fitted-model residuals, the dedicated acorr_ljungbox interface supports a model-degrees-of-freedom adjustment.

What changed upstream

The related statsmodels PR #10017 aligned the non-missing count across the affected paths. The implementation supplies the non-missing count to the Bartlett and non-Bartlett interval paths and to q_stat when the missing mode is drop or conservative. Inputs using none or raise, and inputs without NaNs, retain their previous count.

The patch added four fixture-backed regression checks: confidence intervals and Q statistics for both missing-data modes. A captured focused run of test_stattools.py at the fix tip reported 223 passing tests and 2 expected failures. Those checks establish agreement with the checked-in references covered by the patch.

The stable documentation linked in this article identifies version 0.14.6 and predates the merged implementation. It describes the API contract, not the release containing this change. Check the version or commit you run when reproducing the before-and-after behavior.

Limits

The repaired count makes the calculations internally consistent with the documented convention. Several statistical questions remain outside that statement.

• drop changes the clock. Survivors become contiguous. This is defensible only when that interpretation fits the data-generating process.
• conservative is a different estimator. It preserves positions. With interior gaps and adjusted=True, autocovariances can use lag-specific valid-pair denominators while the interval and Q paths receive one global non-missing count. The patch aligns the stated count contract; it does not prove exact missing-data inference for this case.
• The default lag choice has another edge. The default nlags is computed from raw length before the effective count is known; under heavy missingness, returned arrays can be truncated. That behavior was outside this patch.

Operational checklist

When a time-series diagnostic accepts missing values, a practical review uses five checks:

1. Write down the estimator’s treatment of gaps: delete, preserve, impute, or reject.
2. Track the sample count used by the point estimate, standard error, interval, and test statistic separately.
3. Compare the missing-data path with an explicit reference transformation when the documented contract permits it.
4. Derive the expected ratio before trusting rounded output; distinguish identities from approximations and observed values.
5. Test the whole return tuple. Matching coefficients do not validate their intervals, Q statistics, or p-values.

For this bug, the decisive invariant was simple: under missing="drop", the library call and manual deletion should agree across the returned ACF, confidence interval, Q statistic, and p-value. Once that invariant covered the entire result, the stale count had nowhere left to hide.

Primary sources

• statsmodels documentation for acf, acovf, q_stat, and acorr_ljungbox.
• M. S. Bartlett, “On the Theoretical Specification and Sampling Properties of Autocorrelated Time-Series” (1946), DOI 10.2307/2983611, with a 1948 corrigendum. This is the source for the large-sample sampling-covariance basis and 1/n scaling, not evidence about statsmodels’ missing-data implementation.
• G. M. Ljung and G. E. P. Box, “On a measure of lack of fit in time series models” (1978), DOI 10.1093/biomet/65.2.297. This is the source for the modified portmanteau statistic and its model-diagnostic framing.