When your errors aren’t equal

{/* Assets BUILT 2026-07-06 (see docs/plans/2026-07-06_tsbootstrap-wb-media-plan.md):
/assets/terminal/tsb-wb-signflip.{mp4,webm,frame1.png} (HyperFrames triplet; source at
docs/media-src/tsb-wb-signflip/gen.mjs); the phone cut (panels stacked vertically,
served below 620px via portraitSrc) is tsb-wb-signflip.portrait.{mp4,webm,frame1.png}
from gen-portrait.mjs (same seeded data and beats, 1080x1920). Plus SVG twins in
src/figures/charts/ for
tsb-wb-se-recovery, tsb-wb-multipliers (tier="text"), tsb-wb-blockwild-acf.
The .png src paths are twin references resolved to inline SVG by Figure.astro at build
(no flat PNGs are published, matching the A1 se-collapse convention). Do NOT regenerate
these assets blind — they are tuned; see the media plan before touching. */}

TL;DR: The residual bootstrap resamples your model's errors by shuffling them across time, which silently assumes every error was drawn from the same distribution. When the error variance shifts over time (larger in some regimes than others) that shuffle averages the variance profile away, and the price is a standard error that is too small exactly where it matters. On a heteroskedastic autoregression I built, the shuffle recovered on average only 78% of the true standard error of the slope, and its nominal 90% confidence interval covered the truth 79.1% of the time. The wild bootstrap (keep each residual at its own time position, multiply it by a random sign) recovered 100% of the true standard error and lifted coverage to 88.8%. Same model and data; one argument changed.

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.

The same fitted residuals, resampled two ways. The wild bootstrap keeps every residual at its own time position and only randomizes its sign, so the loud half of the series stays loud in every replicate, and the refitted slopes spread out to the honest width. The i.i.d. shuffle relocates the residuals, flattening the variance profile, and reports a standard error a quarter too small.

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:

$$
e_t^* = v_t \, \hat{e}_t
$$

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

$$
\mathbb{E}[v] = 0, \qquad \mathbb{E}[v^2] = 1 .
$$

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.

The coin flip is the actual literature recommendation, crude as it looks. Davidson and Flachaire (2008) found the Rademacher version generally performs best in practice even on skewed errors, which is why it is the default here and in most implementations. Mammen (1993) supplied the theory for the skewness-preserving alternative.

Three multiplier distributions, one contract. All satisfy the two moment conditions; Mammen alone adds the third, preserving skewness. On this article’s data all three are statistically indistinguishable, so read this as a glossary of mechanisms; on this data there is no leaderboard.

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.

Mean recovered fraction of the true slope standard error, 500 heteroskedastic AR(1) series, 999 replicates each. The i.i.d.-residual shuffle (orange) stalls at 78% of the truth; every wild variant (blue) sits at the grey truth line. The three multipliers are indistinguishable here; the whole gap sits between shuffling and flipping.

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:

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.

This block-constant scheme is the piecewise-constant case of Shao's (2010) dependent wild bootstrap, and the same construction as the wild cluster bootstrap of Cameron, Gelbach and Miller (2008), with time blocks playing the role of clusters.

What a block multiplier preserves. With block length 10 (Rademacher multipliers, 2,000 replicates), the sign flips applied to nearby residuals are almost perfectly correlated (0.904 at lag 1), half-correlated at lag 5, and independent past the block edge; the measured correlations track the triangular law from the construction. This is what lets block-wild keep local residual dependence that independent flips would destroy.

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:

your data's situation	the assumption at risk	reach for
dependent observations, no model you trust	independence of observations	block bootstrap (moving, circular, stationary)
trusted model, steady error variance	none (defaults hold)	i.i.d.-residual bootstrap
trusted model, error variance shifts over time	homoskedasticity of errors	wild bootstrap, Rademacher default
shaky model, residuals locally dependent	independence of errors	block-wild, then question the model

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.

{ q: "What is the wild bootstrap?", a: "The wild bootstrap is a resampling method for regression and time series models with heteroskedastic errors, meaning error variance that shifts over the sample. Instead of shuffling the fitted residuals across time, it keeps each residual at its own position and multiplies it by an independent random draw with mean zero and variance one (by default a Rademacher sign flip: +1 or −1 with equal probability). That preserves each residual’s variance and position exactly, so the bootstrap reproduces the true variance profile without ever estimating it. In tsbootstrap you select it by passing innovation=Wild() to ResidualBootstrap." },
{ q: "When does the ordinary residual bootstrap fail?", a: "When the errors are heteroskedastic (their variance shifts over time), and especially when the loud errors coincide with extreme regressor values. Shuffling residuals across time averages the variance profile away, so the bootstrap understates the uncertainty. On a heteroskedastic AR(1) (slope 0.5, *n* = 600, error variance 9x larger in the second half), the i.i.d.-residual bootstrap recovered on average only 78% of the true standard error of the slope across 500 simulated series, and its nominal 90% interval covered the truth 79.1% of the time (tsbootstrap 0.6.1, 999 replicates per bootstrap, measured 2026-07-06). The wild bootstrap recovered 100% of the true standard error on the same series, with 88.8% coverage." },
{ q: "Which wild bootstrap multiplier should I use: Rademacher, Mammen, or Gaussian?", a: "Rademacher (random sign flips) is the standard default, following Davidson and Flachaire (2008), who found it generally performs best in practice. Mammen’s two-point distribution additionally preserves the skewness of each residual (its third moment equals 1), which theory favors for strongly asymmetric errors; the Gaussian multiplier is a smooth alternative. In the measured experiment here all three recovered 99 to 100% of the true standard error; the choice that matters is wild versus shuffle, not the multiplier." },
{ q: "What is the block-wild bootstrap and when is it needed?", a: "The block-wild bootstrap gives every residual within a contiguous time block the same multiplier, preserving local correlation among residuals, useful when the model is misspecified enough that its residuals are still dependent. It is the piecewise-constant case of Shao’s (2010) dependent wild bootstrap. It is not a free upgrade: on a well-specified model with independent residuals, plain wild is the better choice, and in tsbootstrap BlockWild(block_length=\"auto\") will resolve to a very short block in that case." },
]} />

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.