fix: convergence_advanced.heidelberger_welch_advanced non-stationary fallback ignores autocorrelation in MCSE

[AlexFiliakov]

Problem Statement

convergence_advanced.py heidelberger_welch_advanced (lines 334-342) has a non-stationary fallback path that computes MCSE as np.std(chain) / np.sqrt(n). This is the IID standard error formula, which ignores autocorrelation in the chain. For a non-stationary chain (which by definition has strong serial dependence), the true MCSE can be orders of magnitude larger.

Note: #783 covers the ddof=0 issue in this same line, but the deeper problem -- using an IID formula for a correlated chain -- is not addressed.

v1.0 Impact

High. When the Heidelberger-Welch stationarity test fails, the fallback computes a far-too-small MCSE, potentially leading downstream consumers to believe the mean estimate is precise when it is not. The halfwidth is set to 1.96 * mcse (line 339), and relative_halfwidth is set to inf (line 340), so the halfwidth test correctly fails. However, the returned mcse value itself is wrong and may be used by callers.

Affected Code

ergodic_insurance/convergence_advanced.py, lines 334-342:

else:
    # Non-stationary chain
    stationary_chain = chain
    mean_est = np.mean(chain)
    mcse = np.std(chain) / np.sqrt(n)  # IID formula -- ignores autocorrelation
    halfwidth = 1.96 * mcse
    relative_halfwidth = np.inf
    halfwidth_passed = False
    tau = np.nan

Current Behavior

For a non-stationary chain with autocorrelation rho = 0.9 and N = 1000:

• IID MCSE: std / sqrt(1000) = std * 0.032
• True MCSE: std * sqrt(tau / N) where tau ~ (1 + 2*0.9/(1-0.9)) = 19, giving std * sqrt(19/1000) = std * 0.138
• The IID MCSE underestimates the true MCSE by a factor of ~4.3x

Expected Behavior

Even for a non-stationary chain, the MCSE should account for autocorrelation. Use the spectral density at zero or batch means method:

else:
    stationary_chain = chain
    mean_est = np.mean(chain)
    # Use batch means MCSE even for non-stationary chains
    ess = self.calculate_ess_batch_means(chain)
    mcse = np.std(chain, ddof=1) / np.sqrt(max(1, ess))
    halfwidth = 1.96 * mcse
    relative_halfwidth = np.inf  # Still fail the halfwidth test
    halfwidth_passed = False
    tau = n / max(1, ess)

Alternative Solutions

1. Use batch means ESS for the non-stationary fallback MCSE (recommended)
2. Return mcse = np.inf for non-stationary chains (conservative but loses information)
3. Attempt to find a stationary sub-chain before giving up entirely

Recommended Approach

Option 1: use batch means to estimate ESS and compute a more realistic MCSE. The halfwidth test should still fail (relative_halfwidth = inf), but the returned mcse value will be meaningful.

Acceptance Criteria

• Non-stationary fallback MCSE accounts for autocorrelation
• halfwidth_passed remains False for non-stationary chains
• returned mcse value is consistent with batch means estimate
• Unit test with AR(1) chain (rho=0.95) showing correct MCSE

References

• Heidelberger, P. & Welch, P.D. (1983). 'Simulation Run Length Control in the Presence of an Initial Transient.' Operations Research, 31(6), 1109-1144.
• Flegal, J.M. & Jones, G.L. (2010). 'Batch means and spectral variance estimators in Markov chain Monte Carlo.' Annals of Statistics, 38(2), 1034-1070.
• Related: Math: heidelberger_welch_advanced non-stationary fallback uses ddof=0 for MCSE #783 (ddof=0 in this same line -- narrower fix)