When your errors aren’t equal

I built a time series whose errors I controlled completely, because I wrote the noise myself: an AR(1) with slope 0.5 (each value is half the previous value plus a fresh random shock, so the 0.5 is a slope you can fit like any regression), six hundred points long, where the shocks in the second half are three times louder than in the first, nine times the variance. Then I fit the model, asked a residual bootstrap for the standard error of that slope (the typical miss between the fitted 0.5 and the truth) and compared it to the truth, which I can compute to four decimal places by brute force because I own the data-generating process. The truth is 0.0453. The bootstrap said 0.034.

That is not a rounding disagreement. It is a quarter of the uncertainty, gone, on a series where nothing was misspecified and the bootstrap did exactly what it was designed to do. In the first piece in this series, the ordinary bootstrap lied by shredding autocorrelation, and the block bootstrap fixed it by resampling contiguous chunks. This is the second failure mode, and it lives one layer deeper: not in how you resample the observations, but in how you resample the model’s errors.

The bootstrap under the bootstrap¶

The block methods from the first article resample the observed series directly. There is a second family, the residual bootstrap, that goes through a model instead: fit the model, collect the leftover errors it could not explain (the residuals), resample those, and rebuild synthetic series by running the fitted model forward on the resampled errors. It is the natural choice when you trust the model’s structure and want uncertainty on its parameters, which is precisely the situation of anyone bootstrapping a regression slope, an autoregression coefficient, or a fitted forecast.

And how does it resample the residuals? The default answer is the same default as always: independently and uniformly, with replacement. Pick a residual at random from anywhere in the series, put it back, pick another. Which is a mathematical statement that the residuals are exchangeable, that the error the model made at time 12 could just as well have happened at time 580, because they all came from one common distribution.

There is a name for the property that assumption requires: homoskedasticity, from the Greek for “same scatter”: the error variance is constant over the whole series. Its violation, heteroskedasticity (“different scatter”), shows up everywhere: financial returns in calm and turbulent regimes, sensor noise at different temperatures, forecast errors between weekdays and holidays. If your data was collected over any interesting stretch of time, the odds are good its noise level did not hold still. The check is cheap: plot the squared or absolute residuals in time order and look for a level that drifts or steps.

The tight number I believed¶

I did not go looking for this failure. It turned up in a mean-reversion estimate I once made on a spread whose noise regime changed partway through the sample: a venue migration, roughly doubling the tick-to-tick noise from that date on. I knew better than to trust a textbook formula, so I bootstrapped the reversion coefficient’s standard error from the fitted model’s residuals and got a satisfyingly tight number. What I did not register was that my resampler was lifting residuals from the loud months and planting them in the quiet months, and vice versa, building synthetic histories in which the noise regime never changed at all. The tight standard error was a measurement of those counterfeit series. The strategy sized off that estimate was, in hindsight, sized off roughly three quarters of the real uncertainty; it survived, but by luck rather than arithmetic, and I only worked out why long after the trade was closed.

As in the first article, nothing looked wrong: a residual bootstrap on a fitted model is the responsible choice. The failure has no error message.

Shuffling averages the variance away¶

The mechanism fits in one sentence: shuffling residuals averages the variance away; flipping their signs keeps it in place.

In my built series, the second-half residuals are three times larger than the first-half ones, and, this is the load-bearing detail, the second half is also where the series itself swings widest, since louder shocks push the process further from its mean. The uncertainty of a regression slope is dominated by the errors that sit on the extreme values of the regressor (the input on the x-axis of the fit, which for an AR(1) is simply yesterday’s value). A big error on a big swing tilts the fitted line hard. Heteroskedasticity that is correlated with the regressor level concentrates the largest errors exactly where they do the most damage, and an honest standard error has to account for that.

The i.i.d. shuffle (independent and identically distributed) destroys the correlation. It takes the loud residuals and scatters them uniformly across time, so in the synthetic series a big error is no likelier to land on a big swing than a quiet one is. Averaged over replicates, every time position experiences the average error variance: the profile is flattened and the slope wobbles less than it truly would. The bootstrap then reports that reduced wobble as your standard error. On my series: 0.034 reported, 0.0453 true.

Flip the signs instead¶

Now the second half of the sentence. The wild bootstrap refuses to move anything. Each synthetic error is the original residual, at its original time position, multiplied by a fresh random draw:

where the multiplier $v_t$ is drawn independently at each $t$ from a distribution satisfying two moment conditions:

Mean zero makes the synthetic errors centered, like errors should be. Variance one is the trick: it means $e_t^*$ has exactly the variance of $\hat{e}_t$, at every single $t$ (the multiplier is drawn independently of the residual, so $\operatorname{Var}(v_t \, \hat{e}_t) = \mathbb{E}[v^2] \, \hat{e}_t^2 = \hat{e}_t^2$). A loud residual stays loud, a quiet one stays quiet, and the whole variance profile, whatever shape it has, known or unknown, is preserved by construction. You never estimate the heteroskedasticity; you just decline to destroy it.

The simplest distribution meeting both conditions is the Rademacher: $v = +1$ or $-1$, a fair coin. Under it, the wild bootstrap literally just flips residual signs, which sounds too crude to work, except that a sign flip is exactly as much randomness as a resampler needs: the refit cares only about which way each error pushes the line and how hard, and the flip randomizes the direction while the residual keeps the how-hard. Some multipliers add a third condition, $\mathbb{E}[v^3] = 1$, which additionally preserves each residual’s skewness (the asymmetry of the errors, not just their size). That is why the Mammen two-point distribution exists, with its odd-looking golden-ratio values: it is the two-point law that satisfies all three conditions at once. The working answer: use the coin, and reach for Mammen only when the errors are strongly skewed and the asymmetry of the interval matters.

Five hundred series, one verdict¶

One series proves nothing, so here is the experiment in full, self-contained enough to check. Generate 500 fresh series from the heteroskedastic AR(1) above (slope 0.5, n = 600, error variance 9x larger in the second half). For each, fit the model and bootstrap the standard error of the lag-1 slope four ways (the i.i.d.-residual shuffle and the wild bootstrap under each of the three multipliers) with 999 replicates each. Divide every bootstrap standard error by the true one, 0.045326, computed from 200,000 independent Monte-Carlo draws of the process. A perfect resampler averages 1.0. All of it ran on tsbootstrap 0.6.1 on one 8-core Linux workstation, measured 2026-07-06.

The i.i.d.-residual bootstrap recovered on average 78% of the true standard error (mean ratio 0.783). The wild bootstrap with Rademacher multipliers recovered 100% (mean ratio 1.001); Mammen and Gaussian multipliers landed at 0.997 and 0.990. And the miss was unanimous: the i.i.d. standard error came out below the wild one in 500 of 500 series: the wild estimate ran 1.28x larger on average, and the gap is exactly the uncertainty the shuffle deleted.

Two scope notes. The 78% belongs to this process, where the loud errors sit on the big regressor values; heteroskedasticity that is uncorrelated with the regressor bites far less. And on a genuinely homoskedastic, well-specified model the i.i.d.-residual bootstrap is fine: the shuffle only lies when there is a variance profile to flatten. The wild bootstrap removes the question: flipping works whether or not the shuffle would have failed.

What a collapsed standard error costs you downstream¶

A standard error is an intermediate quantity; the thing you act on is the interval. So: 2,000 more fresh series from the same process, a nominal 90% percentile confidence interval (read straight off the replicates’ percentiles) for the slope from each method, 999 replicates, and a count of how often the interval contains the true 0.5.

The i.i.d.-residual intervals covered the truth 79.1% of the time. The wild (Rademacher) intervals covered 88.8%, running 1.27x wider on average, width the shuffle had been withholding.

And no, 88.8% is not 90%. Part of the residual gap has nothing to do with variance at all: the least-squares slope of an AR(1) is biased slightly downward in finite samples (across the 200,000 truth draws it averaged 0.4947 against a true 0.5), and percentile intervals at 999 replicates carry their own small-sample wobble. Both leftovers belong to the estimator and the percentile rule, not to the resampling. The wild bootstrap recovers the standard error and most of the coverage; it does not repair the estimator it is wrapped around. If the first article taught suspicion of intervals that are too tight, the habit worth keeping from both is asking which assumption each resampler restores, and checking coverage when the answer matters.

The fix is one argument¶

In tsbootstrap the residual bootstrap is a method spec, and how the residuals are resampled is its innovation, a typed argument:

import numpy as np
from tsbootstrap import (
    AR, IID, Wild,
    ResidualBootstrap, bootstrap)

# hetero AR(1): slope 0.5, n=600,
# shock sd 1 then 3 (9x variance)
def hetero_ar1(n, phi, seed):
    rng = np.random.default_rng(seed)
    z = rng.standard_normal(n)
    late = np.arange(n) >= n // 2
    sd = np.where(late, 3.0, 1.0)
    x = np.empty(n)
    x[0] = z[0]/np.sqrt(1 - phi**2)
    for t in range(1, n):
        x[t] = phi*x[t-1]+sd[t]*z[t]
    return x

x = hetero_ar1(600, 0.5, seed=42)

def slope_se(paths):
    # refit the slope per replicate
    slopes = [
      np.polyfit(p[:-1], p[1:], 1)[0]
      for p in paths]
    return np.std(slopes, ddof=1)

# the lie: shuffle residuals in time
res = bootstrap(
    x,
    method=ResidualBootstrap(
        model=AR(order=1),
        innovation=IID()),
    n_bootstraps=999,
    random_state=0,
)
se = slope_se(res.values())
print(round(se, 3))
# slope SE: 0.034   (truth 0.0453)

# the fix: keep each residual in
# place, flip its sign at random
res = bootstrap(
    x,
    method=ResidualBootstrap(
        model=AR(order=1),
        innovation=Wild()),
    n_bootstraps=999,
    random_state=0,
)
se = slope_se(res.values())
print(round(se, 3))
# slope SE: 0.045   (truth 0.0453)

Wild() defaults to Rademacher multipliers; Wild(distribution="mammen") and "gaussian" are one keyword away. Everything else (model, data, replicate count, seed) is identical between the lie and the fix. As of this writing (July 2026), none of arch, darts, MAPIE, or TorchCP ships a wild-bootstrap innovation spec or a block-wild variant as a first-class resampling option.

When the residuals themselves are dependent¶

One assumption is still hiding in the wild bootstrap: drawing each $v_t$ independently asserts that the residuals, whatever their variances, are at least uncorrelated with each other. For a well-specified model that is what residuals should be. But when the model is a little wrong (an underfit autoregression, structure at a horizon the model does not see) the leftover errors are themselves locally dependent, and independent sign flips would shred that dependence the same way the i.i.d. bootstrap shredded the wave in the first article. The concrete picture: fit an AR(1) to a series that is truly AR(3), and the two lags the model never sees leave ripples behind: each residual is correlated with its own recent past, and the residual autocorrelations show it at short lags.

The extension is the block-wild bootstrap: partition the time axis into contiguous blocks and give every residual in a block the same multiplier. Within a block, relative signs (and therefore local correlation) are preserved exactly; across blocks, the flips are independent. You can measure the machinery doing its job directly, on the multiplier process itself: at block length 10, neighboring multipliers correlate at 0.904, multipliers five steps apart at 0.513, and twelve steps apart (past the block edge) at 0.007, indistinguishable from zero. The measured points sit on the triangular law $\max(0,\, 1 - h/L)$ (correlation at lag h for block length L) that the block construction implies.

Block-wild is not a free upgrade over plain wild, though. On this article’s process the residuals are independent by construction, so there is no local dependence to preserve, and block-wild duly buys nothing. Its automatic block-length rule resolved to a length of 2 on my showcase series, and its standard-error recovery there, 0.856 of the truth, actually sat below plain wild’s. Blocking the multipliers spends resampling variety to protect a structure, so when the structure is absent it costs variety for no gain. So the order of operations is explicit. First, fix the model: dependent residuals are a misspecification symptom before they are a resampling problem, and a refit that absorbs the missing lags beats any bootstrap patch. When you cannot refit, block-wild is the patch that respects what the residuals still carry. And when the residuals pass the checks (no drift in their spread, no short-lag correlation) plain wild is all you need.

Which resampler, then?¶

The decision this piece adds to the last one’s, in one table:

The first article’s thesis holds here unchanged: match the resampler to the assumption your data violates. Only the assumption moved. That failure lived in the observations; this one lives in the errors of a model you already trusted, which is a more dangerous address, because trusting the model is what made the analysis feel careful. My tight little 0.034 came from a fitted model and 999 replicates of honest-looking work. The truth, 0.0453, was never hidden. It was flattened, averaged into the quiet half of a history that never happened. The wild bootstrap’s whole contribution is refusing to build that history: every residual stays home, and the only randomness is a sign.

The next failure in this series has nothing to do with statistics: when you need a million replicates, you have to be able to build them at all.