Bootstrap Methods in Econometrics

Resampling, Inference, and Robust Standard Errors
using R, Python & Stata

Athanassios Stavrakoudis

astavrak@uoi.gr

Applied Informatics and Computational Economics Lab

2026-06-29

A Toy Example First

κεἰ σμικρόν ἐστι, σπέρμʼ ἰδεῖν βουλήσομαι.

even if it is small, I will want to see the seed

Σοφοκλῆς, Οἰδίπους Τύραννος 1077

Toy Example — The Mean of a Normal Sample

Setup: Draw a single sample of size \(N=50\) from \(\mathcal{N}(\mu=5,\, \sigma=2)\).

Quantity of interest: the population mean \(\mu\).

Estimator: the sample mean \(\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i\).

• Theory
• Bootstrap

For \(X_i \overset{iid}{\sim} \mathcal{N}(\mu, \sigma^2)\) the sampling distribution of the sample mean is known exactly:

\[\bar{X}_n \sim \mathcal{N}\!\left(\mu,\, \frac{\sigma^2}{n}\right)\]

\[SE(\bar{X}_n) = \frac{\sigma}{\sqrt{n}} = \frac{2}{\sqrt{50}} \approx 0.283\]

\[\text{95% CI} = \bar{X}_n \pm 1.96 \cdot \frac{\sigma}{\sqrt{n}}\]

In practice \(\sigma\) is unknown — replace with \(\hat\sigma\) (sample SD) and use \(t_{n-1}\) critical values:

\[\text{95% CI}_t = \bar{X}_n \pm t_{n-1,\,0.975} \cdot \frac{\hat\sigma}{\sqrt{n}}\]

We do not assume any distributional form. Instead, we let the data approximate \(F\) via the empirical distribution \(\hat{F}_n\):

1. Draw \(\mathbf{X}^*_b = (X_1^*, \ldots, X_n^*)\) with replacement from the observed sample
2. Compute \(\bar{X}^*_b = \frac{1}{n}\sum_{i=1}^n X_i^*\)
3. Repeat \(B\) times
4. Use \(\{\bar{X}^*_b\}_{b=1}^B\) to estimate SE and CI

\[\widehat{SE}_B = \sqrt{\frac{1}{B-1}\sum_{b=1}^B (\bar{X}^*_b - \overline{\bar{X}^*})^2}\]

\[\text{95% CI}_{\text{perc}} = \left(\bar{X}^*_{(0.025)},\; \bar{X}^*_{(0.975)}\right)\]

No formula for \(SE\) is needed — the algorithm is identical for any estimator (median, ratio, IV coefficient, etc.).

Toy Example — Code & Results

• R
• Python
• Stata

# Step 1: draw one sample of size N = 50 from N(5, 2²)
set.seed(14159)
N     <- 50L
mu    <- 5
sigma <- 2
x     <- rnorm(N, mu, sigma)

# Step 2: theoretical SE and CI (assumes σ known)
se_true <- sigma / sqrt(N)
ci_true <- mu + c(-1.96, 1.96) * se_true

# Step 3: classical t-based SE and CI (σ replaced by sample SD)
xbar <- mean(x)
se_t <- sd(x) / sqrt(N)
ci_t <- xbar + qt(c(0.025, 0.975), df = N - 1) * se_t

# Step 4: bootstrap — one replicate at a time via map_dbl
one_boot_mean <- function(b) mean(sample(x, N, replace = TRUE))

B_toy     <- 5000L
boot_mean <- purrr::map_dbl(seq_len(B_toy), one_boot_mean)

se_b <- sd(boot_mean)
ci_b <- quantile(boot_mean, c(0.025, 0.975))

# Step 5: tabulate the three methods side-by-side
tibble(
  Method     = c("Theoretical (σ known)",
                 "Classical t (σ̂)",
                 sprintf("Bootstrap (B = %d)", B_toy)),
  Estimate   = c(mu,         xbar,    mean(boot_mean)),
  SE         = c(se_true,    se_t,    se_b),
  `CI lower` = c(ci_true[1], ci_t[1], ci_b[1]),
  `CI upper` = c(ci_true[2], ci_t[2], ci_b[2])
) %>%
  mutate(across(where(is.numeric), \(x) round(x, 4))) %>%
  kbl(caption = "Toy example: SE & 95% CI for the sample mean (N = 50)") %>%
  kable_styling(font_size = 22, full_width = TRUE)

Toy example: SE & 95% CI for the sample mean (N = 50)

Method	Estimate	SE	CI lower	CI upper
Theoretical (σ known)	5.0000	0.2828	4.4456	5.5544
Classical t (σ̂)	5.0105	0.2551	4.4980	5.5231
Bootstrap (B = 5000)	5.0096	0.2559	4.4969	5.5081

import numpy as np
import pandas as pd
from scipy.stats import t as t_dist

# Step 1: draw one sample of size N = 50 from N(5, 2²)
rng           = np.random.default_rng(14159)
N, mu, sigma  = 50, 5.0, 2.0
x             = rng.normal(mu, sigma, N)

# Step 2: theoretical SE and CI (assumes σ known)
se_true = sigma / np.sqrt(N)
ci_true = (mu - 1.96*se_true, mu + 1.96*se_true)

# Step 3: classical t-based SE and CI (σ replaced by sample SD)
xbar    = x.mean()
se_t    = x.std(ddof=1) / np.sqrt(N)
tcrit   = t_dist.ppf(0.975, df=N-1)
ci_t    = (xbar - tcrit*se_t, xbar + tcrit*se_t)

# Step 4: bootstrap — one replicate per row
def one_boot_mean(b):
    return rng.choice(x, N, replace=True).mean()

B_toy     = 5000
boot_mean = np.array([one_boot_mean(b) for b in range(B_toy)])
se_b      = boot_mean.std(ddof=1)
ci_b      = np.percentile(boot_mean, [2.5, 97.5])

# Step 5: tabulate the three methods
results = pd.DataFrame({
    "Method"    : ["Theoretical (σ known)",
                   "Classical t (σ̂)",
                   f"Bootstrap (B = {B_toy})"],
    "Estimate"  : [mu,         xbar,    boot_mean.mean()],
    "SE"        : [se_true,    se_t,    se_b],
    "CI lower"  : [ci_true[0], ci_t[0], ci_b[0]],
    "CI upper"  : [ci_true[1], ci_t[1], ci_b[1]],
})
print(results.round(4).to_string(index=False))

               Method  Estimate     SE  CI lower  CI upper
Theoretical (σ known)    5.0000 0.2828    4.4456    5.5544
     Classical t (σ̂)    4.6936 0.3081    4.0744    5.3128
 Bootstrap (B = 5000)    4.7012 0.3022    4.1097    5.2909

# Output:
#                  Method  Estimate      SE  CI lower  CI upper
#   Theoretical (σ known)    5.0000  0.2828    4.4457    5.5543
#         Classical t (σ̂)    5.1742  0.2901    4.5914    5.7569
#      Bootstrap (B = 5000)   5.1731  0.2864    4.6157    5.7355

clear all
set seed 14159
set obs 50

* Step 1: draw the sample from N(5, 2^2)
gen x = rnormal(5, 2)

* Step 2: classical SE and 95% CI (uses t_{n-1} critical value)
mean x
* Output:
*    Mean estimation              Number of obs = 50
*    ---------------------------------------------------
*                |  Mean    Std. err.   [95% conf. interval]
*    ------------+--------------------------------------
*              x | 5.1742    0.2901       4.5914   5.7569

* Step 3: bootstrap SE and percentile CI
bootstrap _b[x], reps(5000) seed(14159) nodots: mean x

* Step 4: ask for all three CI types
estat bootstrap, all
* Output:
*                      Observed   Bootstrap                       
*                          mean    std. err.   [95% conf. interval]
*    -----------------------------------------------------------------
*    Mean: x          5.174165    0.286412     4.6155   5.7327   (N)
*                                              4.6157   5.7355   (P)
*                                              4.6109   5.7295   (BC)
*    (N) normal, (P) percentile, (BC) bias-corrected

Toy Example — Visual Comparison

What the plots show

The bootstrap distribution mimics the sampling distribution of \(\bar{X}\) — as if \(\bar{X}\) were the truth. Note that it is centred on \(\bar{X}\), not on \(\mu\). The CLT-style normal shape is visible: the bootstrap reproduces the theoretical \(\mathcal{N}(\mu, \sigma^2/n)\) density even though the algorithm never used it.

• R
• Python
• Stata

set.seed(14159)
N <- 50L; mu <- 5; sigma <- 2
x <- rnorm(N, mu, sigma)
xbar <- mean(x); se_t <- sd(x) / sqrt(N)

one_boot_mean <- function(b) mean(sample(x, N, replace = TRUE))
B_toy     <- 5000L
boot_mean <- purrr::map_dbl(seq_len(B_toy), one_boot_mean)
ci_b      <- quantile(boot_mean, c(0.025, 0.975))

p_a <- ggplot(tibble(x = x), aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), fill = col_main,
                 bins = 12, alpha = 0.55, colour = "white") +
  geom_vline(xintercept = mu,   colour = col_ok,    linetype = "dashed", linewidth = 1) +
  geom_vline(xintercept = xbar, colour = col_accent, linewidth = 1) +
  scale_x_continuous(breaks = scales::pretty_breaks(10)) +
  labs(title = "Observed sample (N = 50)", x = "x", y = "Density")

grid <- seq(xbar - 4*se_t, xbar + 4*se_t, length.out = 300)

p_b <- ggplot() +
  geom_histogram(data = tibble(b = boot_mean),
                 aes(x = b, y = after_stat(density)), fill = col_main,
                 bins = 40, alpha = 0.55, colour = "white") +
  geom_line(data = tibble(x = grid, y = dnorm(grid, mu, sigma/sqrt(N))),
            aes(x = x, y = y), colour = col_ok, linewidth = 1.1) +
  geom_line(data = tibble(x = grid, y = dnorm(grid, xbar, se_t)),
            aes(x = x, y = y), colour = col_accent, linetype = "dashed", linewidth = 1.1) +
  geom_vline(xintercept = ci_b, colour = col_warn, linetype = "dotted", linewidth = 1) +
  scale_x_continuous(breaks = scales::pretty_breaks(10)) +
  labs(title = sprintf("Bootstrap distribution of x̄  (B = %d)", B_toy),
       x = expression(bar(x)^"*"), y = "Density")

(p_a | p_b) & theme(text = element_text(size = 18),
                    plot.title = element_text(size = 18, face = "bold"),
                    legend.position = "none")

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from matplotlib.ticker import MaxNLocator

rng = np.random.default_rng(14159)
N, mu, sigma = 50, 5.0, 2.0
x = rng.normal(mu, sigma, N)
xbar = x.mean()
se_t = x.std(ddof=1) / np.sqrt(N)

B_toy = 5000
boot_mean = np.array([rng.choice(x, N, replace=True).mean() for _ in range(B_toy)])
ci_b = np.percentile(boot_mean, [2.5, 97.5])

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

ax1.hist(x, bins=12, density=True, alpha=0.55, color="#185FA5", edgecolor="white")
ax1.axvline(mu,   color="#1D9E75", linestyle="--", linewidth=1.5)
ax1.axvline(xbar, color="#D85A30", linewidth=1.5)
ax1.set_title("Observed sample (N = 50)", fontsize=18)
ax1.set_xlabel("x", fontsize=18)
ax1.set_ylabel("Density", fontsize=18)
ax1.tick_params(labelsize=16)
ax1.xaxis.set_major_locator(MaxNLocator(nbins=10))

grid = np.linspace(xbar - 4*se_t, xbar + 4*se_t, 300)
ax2.hist(boot_mean, bins=40, density=True, alpha=0.55, color="#185FA5", edgecolor="white")
ax2.plot(grid, norm.pdf(grid, mu,   sigma/np.sqrt(N)), color="#1D9E75", linewidth=1.5)
ax2.plot(grid, norm.pdf(grid, xbar, se_t),             color="#D85A30", linestyle="--", linewidth=1.5)
ax2.axvline(ci_b[0], color="#BA7517", linestyle="dotted", linewidth=1.5)
ax2.axvline(ci_b[1], color="#BA7517", linestyle="dotted", linewidth=1.5)
ax2.set_title(f"Bootstrap distribution of x̄  (B = {B_toy})", fontsize=18)
ax2.set_xlabel("x̄*", fontsize=18)
ax2.set_ylabel("Density", fontsize=18)
ax2.tick_params(labelsize=16)
ax2.xaxis.set_major_locator(MaxNLocator(nbins=10))

plt.tight_layout()
plt.show()

clear
set seed 14159
set obs 50
gen x = rnormal(5, 2)
quietly sum x
local xbar = r(mean)

* Panel 1: observed sample histogram with mean lines
histogram x, density bin(12) fcolor(%55) lcolor(white)     ///
    xline(`xbar', lcolor(orange_red) lwidth(medium))        ///
    xline(5, lpattern(dash) lcolor(emerald) lwidth(medium)) ///
    title("Observed sample (N = 50)", size(medlarge))       ///
    xtitle("x", size(medlarge)) ytitle("Density", size(medlarge)) ///
    xlabel(, labsize(medlarge)) ylabel(, labsize(medlarge)) ///
    legend(off) name(g1, replace)

* Panel 2: bootstrap distribution via bsample loop
set seed 14159
local B = 5000
tempfile bdraws
postfile handle double bmean using `bdraws', replace
forvalues b = 1/`B' {
    preserve
    bsample
    quietly sum x
    post handle (r(mean))
    restore
}
postclose handle

use `bdraws', clear
histogram bmean, density bin(40) fcolor(%55) lcolor(white) ///
    title("Bootstrap distribution of x̄  (B = 5000)", size(medlarge)) ///
    xtitle("x̄*", size(medlarge)) ytitle("Density", size(medlarge))   ///
    xlabel(, labsize(medlarge)) ylabel(, labsize(medlarge))           ///
    legend(off) name(g2, replace)

graph combine g1 g2, cols(2) xsize(13) ysize(5)

Toy Example — Why Does the Bootstrap Work?

Resampling techniques used:

Technique	What it does	When to use
Nonparametric (pairs)	Resample raw observations with replacement	No model assumed; the case we just demonstrated
Parametric	Fit \(\hat{F}_\theta\), draw new samples from it	Strong parametric assumption (e.g. \(\mathcal{N}\) known)
Smoothed	Kernel-smooth the EDF, then sample	Reduces discreteness when \(n\) is small
Residual	Resample fitted residuals \(\hat\varepsilon_i\)	Regression with homoskedastic errors
Wild	\(\hat\varepsilon_i \cdot w_i\), \(\mathbb{E}[w]=0\)	Regression with heteroskedasticity

Conditions for bootstrap consistency (the theoretical guarantee):

\[\sup_{x} \left| \Pr^*\!\left(\sqrt{n}(\bar{X}^* - \bar{X}) \le x\right) - \Pr\!\left(\sqrt{n}(\bar{X} - \mu) \le x\right) \right| \xrightarrow{p} 0\]

i.e. the bootstrap distribution converges (uniformly) to the true sampling distribution, conditional on the data. This is Bickel & Freedman’s (1981) consistency theorem for the sample mean.

Comparison summary — toy example ((N = 50), true (= 5), true (= 2))

The three SE estimates agree to within ~3%. The CIs are nearly indistinguishable. Why bother with the bootstrap when classical (t) works?

→ Because for non-trivial estimators (ratios, quantiles, IV coefficients, IRFs at horizon 5, …), no closed-form SE exists. The bootstrap delivers a usable SE/CI with the same algorithm, unchanged.

Outline

Background & Theory

•  Brief history of bootstrap
•  Why asymptotics sometimes fail
•  The bootstrap principle — plug-in & EDF
•  Bootstrap consistency conditions
•  CI types: normal, basic, percentile, BCa
•  Taxonomy of bootstrap methods
•  Bootstrap weight distributions

Implementation & Applications

•  Serial vs parallel bootstrap
•  Required libraries · DGP · Data preview
•  App 1 — SE and CI estimation
•  App 2 — Bootstrap hypothesis testing
•  App 3 — Wild cluster bootstrap (few clusters)
•  App 4 — IV / 2SLS with weak instruments
•  App 5 — Nonstandard distributions (unit roots)
•  TS 1 — Block bootstrap for AR models
•  TS 2 — VAR impulse-response CI
•  TS 3 — GARCH filtered bootstrap
•  Advantages, limitations & cutting-edge research

Introduction & History

παλαιότης γὰρ τῷ λόγῳ γʼ ἔνεστί τις.

there is a certain antiquity in the story

Εὐριπίδης, Ἑλένη 1056

Brief History of Bootstrap

• Timeline
• Note
• Tip

Year	Author(s)	Contribution
1979	Efron	Bootstrap introduced — doi:10.1214/aos/1176344552
1981	Efron	Bootstrap SE & CI — doi:10.1093/biomet/68.3.589
1986	Wu	Wild bootstrap for heteroskedastic regression — doi:10.1214/aos/1176350142
1988	Liu	Bootstrap under non-iid models — doi:10.1214/aos/1176351062
1989	Künsch	Block bootstrap — doi:10.1214/aos/1176347265
1992	Efron & Tibshirani	An Introduction to the Bootstrap (book, Chapman & Hall)
1993	Mammen	Wild bootstrap, high-dim — doi:10.1214/aos/1176349025
1994	Politis & Romano	Stationary bootstrap — doi:10.1080/01621459.1994.10476870
2001	Horowitz	Bootstrap in econometrics, Handbook of Econometrics Vol. 5 (ch. 52)
2008	Cameron, Gelbach & Miller	Wild cluster bootstrap — doi:10.1162/rest.90.3.414
2018	MacKinnon & Webb	WCB for few treated clusters, Econometrics Journal 21(2), 114–135
2019	Roodman et al.	boottest / fwildclusterboot — doi:10.1177/1536867X19830877
2023	MacKinnon, Nielsen & Webb	Cluster-robust guide — doi:10.1016/j.jeconom.2022.04.001

Note

Efron’s original insight (1979):

“Bootstrap” refers to the story of Baron Münchhausen pulling himself out of a swamp by his own bootstraps. The idea: use the data itself — not a parametric model — to approximate the sampling distribution of any statistic.

The term was deliberately provocative: you cannot literally lift yourself by your own bootstraps — yet the method works. The key insight is that the empirical distribution \(\hat{F}_n\) is the best available estimate of \(F\), and any quantity computable from \(F\) can be estimated by the same computation applied to \(\hat{F}_n\).

Tip

Recognition and impact:

The bootstrap is now one of the most-cited statistical ideas of the 20th century. Bradley Efron received the International Prize in Statistics in 2018 and the US National Medal of Science in 2023, in part for this contribution.

It has been estimated that the bootstrap (and its variants) appears in over 100 000 published articles across statistics, econometrics, biology, medicine, and machine learning.

A short and accessible introduction: Efron & Hastie (2016) Computer Age Statistical Inference, Ch. 10–11 (free PDF available at the link).

Why Do Asymptotics Sometimes Fail?

Classical inference assumes (unreliably) that for a statistic \(\hat\theta_n\): \[\sqrt{n}(\hat\theta_n - \theta_0) \xrightarrow{d} \mathcal{N}(0,\, V)\]

When is this unreliable?

Small sample	Normal approximation inaccurate
Heteroskedasticity	\(\widehat{SE}\) biased unless corrected
Few clusters (\(G < 30\))	CLT doesn’t kick in
Weak instruments (\(F < 10\))	Non-normal distribution of \(\hat\beta_{IV}\)
Near unit root (\(\rho \to 1\))	Dickey-Fuller, not normal
Nonlinear statistics	Variance formula complex or unknown
Boundary / inequality constraints	Asymptotic distribution non-standard

The bootstrap alternative

Instead of relying on ((0,V)), let the data itself estimate the sampling distribution.

No closed-form variance formula needed. No normality assumption. Often higher-order accurate than first-order asymptotics.

Limitations

Bootstrap is not a cure-all. It fails when the statistic is not asymptotically pivotal, at boundary cases, or when the sample size is so small that the EDF is a poor proxy for (F).

Mathematical Framework

ἡ γὰρ φύσις βέβαιος, οὐ τὰ χρήματα.

it is nature that is steadfast, not possessions

Εὐριπίδης, Ἠλέκτρα 941

The Bootstrap Principle

• EDF
• Plug-In
• Algorithm

Empirical distribution function (EDF) as a plug-in estimator of the true CDF \(F\):

\[\hat{F}_n(x) = \frac{1}{n}\sum_{i=1}^n \mathbf{1}\{X_i \le x\}\]

Plug-in principle: if \(\theta = T(F)\) is a functional of the CDF, then \(\hat\theta_n = T(\hat F_n)\).

The sampling distribution of \(\hat\theta_n - \theta_0\) is approximated by the bootstrap distribution of \(\hat\theta^* - \hat\theta_n\), where \(\hat\theta^*\) is computed on a resample \(\mathbf{X}^*\) drawn from \(\hat{F}_n\):

\[\underbrace{F \longrightarrow \mathbf{X} \longrightarrow \hat\theta_n}_{\text{true world}} \;\approx\; \underbrace{\hat{F}_n \longrightarrow \mathbf{X}^* \longrightarrow \hat\theta^*}_{\text{bootstrap world}}\]

Algorithm:

1. Draw \(\mathbf{X}^* = (X_1^*, \ldots, X_n^*)\) iid from \(\hat{F}_n\) (resample with replacement)
2. Compute \(\hat\theta^*_b = T(\hat{F}^*_{b,n})\)
3. Repeat \(B\) times
4. Use \(\{\hat\theta^*_b\}_{b=1}^B\) to approximate the sampling distribution of \(\hat\theta_n\)

The Bootstrap Principle — SE and Bias

From the \(B\) bootstrap replicates \(\{\hat\theta^*_b\}_{b=1}^B\) one can estimate:

Bootstrap standard error:	\[\widehat{SE}_B = \sqrt{\dfrac{1}{B-1}\sum_{b=1}^B\!\left(\hat\theta^*_b - \bar{\theta}^*\right)^2}, \quad \bar{\theta}^* = \dfrac{1}{B}\sum_{b=1}^B \hat\theta^*_b\]
Bootstrap bias estimate:	\[\widehat{\text{bias}}_B = \bar{\theta}^* - \hat\theta_n\]
Bias-corrected estimator:	\[\tilde\theta_n = 2\hat\theta_n - \bar{\theta}^* = \hat\theta_n - \widehat{\text{bias}}_B\]

Why prefer \(\widehat{SE}_B\) over the analytical SE?

When the analytical SE formula is derived under assumptions that may not hold (normality, homoskedasticity, iid), \(\widehat{SE}_B\) requires only that the EDF is a good approximation of \(F\) — a weaker condition.

Convergence

Convergence: the bootstrap approximation improves as \(B \to \infty\) and \(n \to \infty\). For a given \(n\), \(B = 999\) is usually more than sufficient for \(p\)-values and CIs. \(B = 99\) is acceptable for SE estimates only.

Bootstrap Confidence Intervals

Four main approaches, increasing accuracy (and complexity):

Type	Formula	Notes
Normal	\(\hat\theta \pm z_{1-\alpha/2}\cdot\widehat{SE}_B\)	Assumes symmetry; fastest
Basic	\((2\hat\theta - \theta^*_{1-\alpha/2},\; 2\hat\theta - \theta^*_{\alpha/2})\)	Inverts pivot; handles skew
Percentile	\((\theta^*_{\alpha/2},\; \theta^*_{1-\alpha/2})\)	Simplest; symmetric treatment
Studentized (BCa)	\(\left(\hat\theta - \hat{SE}\cdot t^*_{1-\alpha/2},\; \hat\theta - \hat{SE}\cdot t^*_{\alpha/2}\right)\)	Best coverage; second bootstrap

Bias-corrected and accelerated (BCa):

\[\text{CI}_{BCa} = \left(\theta^*_{\alpha_1},\; \theta^*_{\alpha_2}\right)\]

where: \[\alpha_1 = \Phi\!\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{\alpha/2})}\right), \quad \alpha_2 = \Phi\!\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{1-\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{1-\alpha/2})}\right)\]

\(\hat{z}_0\) — bias correction (median of bootstrap distribution relative to \(\hat\theta\)); \(\hat{a}\) — acceleration (jackknife skewness of the influence function)

Note

Recommendation (practical): - Default: percentile-\(t\) (studentized) for best coverage - Quick: BCa via boot::boot.ci(type="bca") - Avoid: normal when distribution is skewed - Avoid: basic for nonlinear statistics

Types of Bootstrap

καὶ πολλὰ καὶ παντοῖʼ ἀκουούσας κακά.

hearing evils many and of every kind

Ἀριστοφάνης, Θεσμοφοριάζουσαι 388

Taxonomy of Bootstrap Methods

Method	Resample unit	When to use	Key reference
Pairs (nonparametric)	\((y_i, \mathbf{x}_i)\) rows	Robust to heterosked.	Efron (1979)
Residual	\(\hat\varepsilon_i\)	Homoskedastic errors	Freedman (1981)
Wild	\(\hat\varepsilon_i \cdot w_i\)	Heterosked., regression	Wu (1986)
Wild cluster	\(\hat\varepsilon_g \cdot w_g\)	Clustered data, few \(G\)	Cameron et al. (2008)
Parametric	From fitted model \(\hat{F}_\theta\)	Parametric model known	—
Moving block (MBB)	Fixed-length blocks	Dependent time series	Künsch (1989)
Circular block (CBB)	Wrap-around blocks	Reduce end-block bias	Politis & Romano (1992)
Stationary	Geometric-length blocks	Second-order correct TS	Politis & Romano (1994)
Sieve bootstrap	AR model residuals	Parametric TS structure	Bühlmann (1997)
Jackknife	Leave-one-out subsamples	Quick SE / bias, small \(n\)	Efron & Tibshirani (1993)

Rule of Thumb

Rule of thumb: pairs bootstrap if the regression error structure is unknown; wild bootstrap if heteroskedasticity is suspected but errors are independent; block/stationary bootstrap if serial dependence is present.

The Jackknife

The jackknife (leave-one-out resampling) is a linear approximation to the bootstrap — cheaper but less accurate. It gets its own section later in the deck, with R/Python/Stata code.

Bootstrap Methods — How to Use Each One

Hands-on syntax for each method. The Pairs tab runs live against bs-cross.csv; remaining tabs are illustrative only.

• Pairs (nonparametric)
• Wild (regression)
• Wild cluster
• Block (time series)
• Sieve (AR-based)
• Filtered Historical Simulation (FHS for GARCH)

• R
• Python
• Stata

set.seed(14159)
stat_fn <- function(data, idx) coef(lm(y ~ x1 + x2, data = data[idx, ]))["x1"]
fit <- boot(data = bs_cross, statistic = stat_fn, R = 999,
            parallel = if (.Platform$OS.type == "unix") "multicore" else "snow",
            ncpus = 12)
ci <- quantile(fit$t, c(0.025, 0.975))
cat(sprintf("95%% percentile CI for beta_x1:  [%.4f,  %.4f]\n", ci[1], ci[2]))

95% percentile CI for beta_x1:  [1.9450,  2.4168]

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from joblib import Parallel, delayed

bs_cross_py = pd.read_csv("../data/bs-cross.csv")
y_arr = bs_cross_py["y"].values
X_arr = np.column_stack([np.ones(len(y_arr)),
                         bs_cross_py["x1"].values,
                         bs_cross_py["x2"].values])

def one_boot(b):
    rng = np.random.default_rng(14159 + b)
    idx = rng.integers(0, len(y_arr), len(y_arr))
    return float(np.linalg.lstsq(X_arr[idx], y_arr[idx], rcond=None)[0][1])

bt = np.array(Parallel(n_jobs=12)(delayed(one_boot)(b) for b in range(999)))
ci = np.percentile(bt, [2.5, 97.5])
print(f"95% percentile CI for beta_x1:  [{ci[0]:.4f},  {ci[1]:.4f}]")

95% percentile CI for beta_x1:  [1.9355,  2.4188]

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace

set seed 14159
bootstrap _b[x1], reps(999) nodots: quietly regress y x1 x2
estat bootstrap, percentile

Critical Parameters — Pairs

• R / reps: 999 or 9 999 — use \(B = 10\,000^x - 1\) so percentile indices fall on integers
• seed: always set for reproducibility
• Parallelism: multicore on Linux/macOS, snow on Windows
• Pairs bootstrap automatically handles heteroskedasticity — it resamples rows, so \(\text{Var}(\varepsilon_i \mid x_i)\) is preserved

• R
• Python
• Stata

set.seed(14159)
fit_w   <- lm(y ~ x1 + x2, data = bs_cross)
h_w     <- hatvalues(fit_w)
e_resc  <- residuals(fit_w) / sqrt(1 - h_w)   # HC3-style rescaling
y_fit_w <- fitted(fit_w)
n_w     <- nrow(bs_cross)

b_wild <- purrr::map_dbl(seq_len(999), function(b) {
  w      <- sample(c(-1, 1), n_w, replace = TRUE)
  y_star <- y_fit_w + e_resc * w
  coef(lm(y_star ~ x1 + x2, data = bs_cross))[2]
})
ci_wild <- quantile(b_wild, c(0.025, 0.975))
cat(sprintf("95%% percentile CI for beta_x1 (wild):  [%.4f,  %.4f]\n",
            ci_wild[1], ci_wild[2]))

95% percentile CI for beta_x1 (wild):  [1.9493,  2.3988]

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from joblib import Parallel, delayed

df_w = pd.read_csv("../data/bs-cross.csv")
y_w  = df_w["y"].values
X_w  = np.column_stack([np.ones(len(y_w)), df_w["x1"].values, df_w["x2"].values])
n_w  = len(y_w)

beta_w  = np.linalg.lstsq(X_w, y_w, rcond=None)[0]
y_fit_w = X_w @ beta_w
resid_w = y_w - y_fit_w
H_w     = X_w @ np.linalg.solve(X_w.T @ X_w, X_w.T)
h_w     = np.diag(H_w)
e_resc  = resid_w / np.sqrt(1 - h_w)

def wild_one(b):
    rng    = np.random.default_rng(14159 + b)
    w      = rng.choice([-1.0, 1.0], size=n_w)
    y_star = y_fit_w + e_resc * w
    return float(np.linalg.lstsq(X_w, y_star, rcond=None)[0][1])

bt_wild = np.array(Parallel(n_jobs=12)(delayed(wild_one)(b) for b in range(999)))
ci_wild = np.percentile(bt_wild, [2.5, 97.5])
print(f"95% percentile CI for beta_x1 (wild):  [{ci_wild[0]:.4f},  {ci_wild[1]:.4f}]")

95% percentile CI for beta_x1 (wild):  [1.9552,  2.4111]

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace
quietly regress y x1 x2
boottest x1, reps(999) weighttype(rademacher) nograph seed(14159)

Weight Distributions & Residual Rescaling

Weights (all satisfy \(\mathbb{E}[w] = 0\), \(\mathbb{E}[w^2] = 1\)):

• Rademacher \(\{-1, +1\}\): default; symmetric, \(\mathbb{E}[w^3] = 0\)
• Mammen \(\{-(\sqrt5-1)/2, (\sqrt5+1)/2\}\): \(\mathbb{E}[w^3] = 1\) (skewness correction)
• Webb (6-point) \(\{\pm\sqrt{3/2}, \pm 1, \pm\sqrt{1/2}\}\): recommended when \(G \le 10\)

Residual rescaling: divide residuals by \(\sqrt{1-h_t}\) (HC3-style) before applying weights. Rescaled residuals have constant variance, so the wild bootstrap mimics the true error variance more accurately than raw residuals. Davidson & Flachaire (2008) show Rademacher performs best when errors are roughly symmetric.

• R
• Python
• Stata

set.seed(14159)
fit_wc <- lm(y ~ treat + x1, data = bs_cluster)
res_wc <- boottest(fit_wc, clustid = "cluster", param = "treat", B = 999L)
print(res_wc)

boottest.lm(object = fit_wc, param = "treat", B = 999L, clustid = "cluster")
 
p value: 0.023 
confidence interval: 0.2415 3.0191 
test statistic 2.9472 

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from joblib import Parallel, delayed

df_c  = pd.read_csv("../data/bs-cluster.csv")
y_c   = df_c["y"].values
X_c   = np.column_stack([np.ones(len(y_c)), df_c["treat"].values, df_c["x1"].values])
cl    = df_c["cluster"].values
G_ids = np.unique(cl)
G     = len(G_ids)

beta_c  = np.linalg.lstsq(X_c, y_c, rcond=None)[0]
y_fit_c = X_c @ beta_c
e_hat_c = y_c - y_fit_c

def cse(r):
    scores = X_c * r[:, None]
    meat   = sum(np.outer(scores[cl==g].sum(0), scores[cl==g].sum(0))
                 for g in G_ids) * G / (G - 1)
    v = np.linalg.solve(X_c.T @ X_c, meat) @ np.linalg.solve(X_c.T @ X_c, np.eye(3))
    return np.sqrt(v[1, 1])

se_obs = cse(e_hat_c)
t_obs  = beta_c[1] / se_obs

def wcb_one(b):
    rng    = np.random.default_rng(14159 + b)
    w_g    = rng.choice([-1.0, 1.0], size=G)
    w      = w_g[np.searchsorted(G_ids, cl)]
    y_star = y_fit_c + e_hat_c * w
    b_j    = np.linalg.lstsq(X_c, y_star, rcond=None)[0][1]
    r_b    = y_star - X_c @ np.linalg.lstsq(X_c, y_star, rcond=None)[0]
    return (b_j - beta_c[1]) / cse(r_b)

t_boot = np.array(Parallel(n_jobs=12)(delayed(wcb_one)(b) for b in range(999)))
p_val  = float(np.mean(np.abs(t_boot) >= np.abs(t_obs)))
q      = np.quantile(t_boot, [0.025, 0.975])
ci_lo  = beta_c[1] - se_obs * q[1]
ci_hi  = beta_c[1] - se_obs * q[0]
print(f"Wild cluster bootstrap  (G={G}, B=999, Rademacher)")

Wild cluster bootstrap  (G=15, B=999, Rademacher)

print(f"  treat estimate: {beta_c[1]:.4f}   t = {t_obs:.4f}   p = {p_val:.4f}")

  treat estimate: 1.7085   t = 2.9521   p = 0.0210

print(f"  95% CI:  [{ci_lo:.4f},  {ci_hi:.4f}]")

  95% CI:  [0.1823,  3.0656]

quietly import delimited "../data/bs-cluster.csv", clear
quietly destring _all, replace
quietly regress y treat x1, vce(cluster cluster)
boottest treat, reps(999) weighttype(rademacher) cluster(cluster) nograph seed(14159)

Critical Parameters — Wild Cluster

• \(G\) (number of clusters): if \(G < 30\), asymptotic cluster-robust SEs over-reject — WCB is essential
• B: use ≥ 9 999 for tail \(p\)-values
• weighttype: webb if \(G \le 10\), rademacher otherwise
• Output: boottest does NOT return a single SE — it inverts a test. Report point estimate, \(t\), \(p\), and CI.

• R
• Python
• Stata

set.seed(14159)
ar1_stat <- function(ts, ...) {
  coef(lm(ts[-1] ~ ts[-length(ts)]))[2]
}
ts_fit <- tsboot(tseries   = bs_ts$y_ar,
                 statistic = ar1_stat,
                 R         = 999,
                 l         = 13,        # block length ≈ 1.75 · T^(1/3)
                 sim       = "fixed")
ci_mbb <- quantile(ts_fit$t, c(0.025, 0.975), na.rm = TRUE)
cat(sprintf("95%% percentile CI for AR(1) rho (MBB, l=13):  [%.4f,  %.4f]\n",
            ci_mbb[1], ci_mbb[2]))

95% percentile CI for AR(1) rho (MBB, l=13):  [0.6107,  0.7543]

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from joblib import Parallel, delayed

df_ts = pd.read_csv("../data/bs-ts.csv")
y_ar  = df_ts["y_ar"].values
T_ar  = len(y_ar)
l_blk = 13  # block length ≈ 1.75 * T^(1/3)

def ar1_coef(y):
    return float(np.linalg.lstsq(
        np.column_stack([np.ones(len(y)-1), y[:-1]]), y[1:], rcond=None)[0][1])

def mbb_one(b):
    rng    = np.random.default_rng(14159 + b)
    starts = rng.integers(0, T_ar - l_blk + 1, int(np.ceil(T_ar / l_blk)))
    y_star = np.concatenate([y_ar[s:s+l_blk] for s in starts])[:T_ar]
    return ar1_coef(y_star)

bt_mbb = np.array(Parallel(n_jobs=12)(delayed(mbb_one)(b) for b in range(999)))
ci_mbb = np.percentile(bt_mbb, [2.5, 97.5])
print(f"95% percentile CI for AR(1) rho (MBB, l={l_blk}):  [{ci_mbb[0]:.4f},  {ci_mbb[1]:.4f}]")

95% percentile CI for AR(1) rho (MBB, l=13):  [0.6194,  0.7548]

* No direct block-bootstrap command in Stata.
* Resample blocks manually using -bsample- in a loop,
* or use R / Python tsboot / MBB instead.

Critical: Block Length \(\ell\)

• Rule of thumb: \(\ell \approx 1.75 \cdot T^{1/3}\) (≈ 13 for \(T = 400\))
• Too short → fails to capture serial dependence
• Too long → too few distinct blocks, high variance
• Use np::b.star() in R or Politis–White–Patton automatic selector

• R
• Python
• Stata

set.seed(14159)
ar_fit  <- ar(bs_ts$y_ar, order.max = 5, method = "ols", aic = TRUE)
p_s     <- ar_fit$order
resid_w <- na.omit(ar_fit$resid) - mean(na.omit(ar_fit$resid))

sieve_one <- function(b) {
  set.seed(14159 + b)
  u_star <- sample(resid_w, length(bs_ts$y_ar), replace = TRUE)
  y_star <- as.numeric(stats::filter(u_star, ar_fit$ar, method = "recursive"))
  tryCatch(ar(y_star, order.max = p_s, method = "ols", aic = FALSE)$ar[1],
           error = \(e) NA_real_)
}
bt_sv <- purrr::map_dbl(seq_len(499), sieve_one)
ci_sv <- quantile(bt_sv, c(0.025, 0.975), na.rm = TRUE)
cat(sprintf("AR order selected: p = %d\n", p_s))
cat(sprintf("95%% percentile CI for AR(1) rho (sieve, B=499):  [%.4f,  %.4f]\n",
            ci_sv[1], ci_sv[2]))

AR order selected: p = 1

95% percentile CI for AR(1) rho (sieve, B=499):  [0.6741,  0.8104]

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from statsmodels.tsa.ar_model import AutoReg
from joblib import Parallel, delayed

df_sv = pd.read_csv("../data/bs-ts.csv")
y_sv  = df_sv["y_ar"].values

fit_sv = AutoReg(y_sv, lags=1, old_names=False).fit()
phi    = float(fit_sv.params[1])
u_w    = fit_sv.resid - fit_sv.resid.mean()

def sieve_one(b):
    rng    = np.random.default_rng(14159 + b)
    u_star = rng.choice(u_w, size=len(y_sv), replace=True)
    y_star = np.zeros(len(y_sv))
    y_star[0] = y_sv[0]
    for t in range(1, len(y_sv)):
        y_star[t] = phi * y_star[t-1] + u_star[t]
    try:
        return float(AutoReg(y_star, lags=1, old_names=False).fit().params[1])
    except Exception:
        return np.nan

bt_sv = np.array(Parallel(n_jobs=12)(delayed(sieve_one)(b) for b in range(499)))
ci_sv = np.nanpercentile(bt_sv, [2.5, 97.5])
print(f"AR(1) coefficient estimate: {phi:.4f}")

AR(1) coefficient estimate: 0.7625

print(f"95% percentile CI for rho (sieve, B=499):  [{ci_sv[0]:.4f},  {ci_sv[1]:.4f}]")

95% percentile CI for rho (sieve, B=499):  [0.6808,  0.8128]

* Sieve bootstrap requires AR model fitting per replicate.
* Implement manually with -bsample- and -arima-, or use R / Python.

Critical: AR Order \(p\)

Choose \(p\) by AIC or BIC. Sieve bootstrap is valid only if the AR(\(p\)) model is correctly specified — misspecification carries through to the bootstrap distribution. For long-memory or fractional series, block bootstrap is safer.

• R
• Python
• Stata

set.seed(14159)
r_g    <- bs_garch$r
fit_g  <- suppressWarnings(tseries::garch(r_g, order = c(1, 1), trace = FALSE))
s_hat  <- as.numeric(fitted(fit_g))   # conditional standard deviations σ_t
valid  <- which(!is.na(s_hat) & s_hat > 0)
z_hat  <- r_g[valid] / s_hat[valid];  z_hat <- z_hat - mean(z_hat)
params <- coef(fit_g)   # a0 = omega, a1 = alpha, b1 = beta

fhs_one <- function(b) {
  set.seed(14159 + b)
  z_star <- sample(z_hat, length(valid), replace = TRUE)
  r_star <- s_hat[valid] * z_star   # fitted vol path × resampled innovations
  g <- tryCatch({
    f  <- suppressWarnings(tseries::garch(r_star, order = c(1,1), trace = FALSE))
    cf <- coef(f)
    if (cf["a1"] + cf["b1"] >= 1) rep(NA_real_, 3) else cf
  }, error = \(e) rep(NA_real_, 3))
  g
}
bt <- do.call(rbind, lapply(seq_len(199), fhs_one))
ok <- !apply(bt, 1, anyNA)
cat(sprintf("Valid replicates: %d / 199\n", sum(ok)))
cat("FHS bootstrap (B=199) — GARCH(1,1)\n")
cat(sprintf("  alpha (ARCH)  est=%.4f  95%% CI [%.4f, %.4f]\n",
            params["a1"],
            quantile(bt[ok, 2], 0.025), quantile(bt[ok, 2], 0.975)))
cat(sprintf("  beta  (GARCH) est=%.4f  95%% CI [%.4f, %.4f]\n",
            params["b1"],
            quantile(bt[ok, 3], 0.025), quantile(bt[ok, 3], 0.975)))

Valid replicates: 199 / 199

FHS bootstrap (B=199) — GARCH(1,1)

  alpha (ARCH)  est=0.2209  95% CI [0.0000, 0.1375]

  beta  (GARCH) est=0.0000  95% CI [0.0000, 0.8543]

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from joblib import Parallel, delayed

df_g = pd.read_csv("../data/bs-garch.csv")
r_g  = df_g["r"].values
T_g  = len(r_g)

def garch11_nll(params, r):
    omega, alpha, beta = np.abs(params)
    h = np.full(len(r), max(float(np.var(r)), 1e-6))
    ll = 0.0
    for t in range(1, len(r)):
        h[t] = omega + alpha * r[t-1]**2 + beta * h[t-1]
        if h[t] <= 0:
            return 1e10
        ll += np.log(h[t]) + r[t]**2 / h[t]
    return 0.5 * ll

res0 = minimize(garch11_nll, [0.05, 0.15, 0.80], args=(r_g,),
                bounds=[(1e-6, 1), (1e-6, 0.999), (1e-6, 0.999)],
                method="L-BFGS-B")
w0, a0, b0 = np.abs(res0.x)

h_hat = np.full(T_g, float(np.var(r_g)))
for t in range(1, T_g):
    h_hat[t] = w0 + a0 * r_g[t-1]**2 + b0 * h_hat[t-1]
s_hat = np.sqrt(np.maximum(h_hat, 1e-8))
z_hat = r_g / s_hat;  z_hat -= z_hat.mean()

def fhs_one(b):
    rng    = np.random.default_rng(14159 + b)
    z_star = rng.choice(z_hat, size=T_g, replace=True)
    r_star = s_hat * z_star
    try:
        res_b = minimize(garch11_nll, [w0, a0, b0], args=(r_star,),
                         bounds=[(1e-6, 1), (1e-6, 0.999), (1e-6, 0.999)],
                         method="L-BFGS-B")
        ab, bb = np.abs(res_b.x[1]), np.abs(res_b.x[2])
        return (ab, bb) if ab + bb < 1 else (np.nan, np.nan)
    except Exception:
        return (np.nan, np.nan)

results = Parallel(n_jobs=12)(delayed(fhs_one)(b) for b in range(199))
bt_a   = np.array([r[0] for r in results])
bt_b   = np.array([r[1] for r in results])
ok     = ~np.isnan(bt_a)
ci_a   = np.nanpercentile(bt_a, [2.5, 97.5])
ci_b_v = np.nanpercentile(bt_b, [2.5, 97.5])
print(f"Valid replicates: {ok.sum()} / 199")

Valid replicates: 199 / 199

print(f"FHS bootstrap (B=199) — GARCH(1,1)")

FHS bootstrap (B=199) — GARCH(1,1)

print(f"  alpha (ARCH)  est={a0:.4f}  95% CI [{ci_a[0]:.4f}, {ci_a[1]:.4f}]")

  alpha (ARCH)  est=0.2748  95% CI [0.0703, 0.2431]

print(f"  beta  (GARCH) est={b0:.4f}  95% CI [{ci_b_v[0]:.4f}, {ci_b_v[1]:.4f}]")

  beta  (GARCH) est=0.6695  95% CI [0.6198, 0.8784]

* Stata's -arch- has a bootstrap option but no built-in FHS.
* Use R (tseries::garch) or Python (scipy) for filtered historical simulation.

Why Filtered, Not Raw?

iid resampling of \(r_t\) destroys volatility clustering, biasing \(\hat\alpha^*\) toward zero (Gonçalves & Kilian, 2004). FHS resamples standardised \(\hat z_t = r_t/\hat\sigma_t\) and reconstructs using the fitted volatility path, preserving the dynamics.

Stationarity check: drop replicates with \(\hat\alpha^* + \hat\beta^* \ge 1\).

Pairs, Residual, and Wild Bootstrap

Pairs bootstrap (nonparametric):	\[\{(y_i^*, \mathbf{x}_i^*)\}_{i=1}^n \overset{iid}{\sim} \hat{F}_n \;\Rightarrow\; \hat\beta^* = (\mathbf{X}^{*\top}\mathbf{X}^*)^{-1}\mathbf{X}^{*\top}\mathbf{y}^*\]
Residual bootstrap (homoskedastic):	\[y_i^* = \mathbf{x}_i^\top\hat\beta + \varepsilon_i^*, \quad \varepsilon_i^* \overset{iid}{\sim} \hat{F}_{\hat\varepsilon} \;\Rightarrow\; \hat\beta^* = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}^*\]
Wild bootstrap (heteroskedasticity-robust):	\[y_i^* = \mathbf{x}_i^\top\hat\beta + \hat\varepsilon_i \cdot w_i, \quad w_i \overset{iid}{\sim} F_w, \quad \mathbb{E}[w_i]=0,\;\mathbb{E}[w_i^2]=1\]

Wild bootstrap weight distributions \(F_w\):

Distribution	Values	Probabilities	\(\mathbb{E}[w^3]\)
Rademacher	\(\{-1,\, +1\}\)	\(\{1/2,\, 1/2\}\)	0
Mammen	\(\{-(\sqrt5-1)/2,\,(\sqrt5+1)/2\}\)	\(\{p^*,\,1-p^*\}\)	1
Webb (6-point)	\(\{\pm\sqrt{3/2},\pm 1,\pm\sqrt{1/2}\}\)	\(\{1/6,\ldots,1/6\}\)	0

where \(p^* = (\sqrt5+1)/(2\sqrt5)\)

Mammen vs Webb Weights

The Mammen weight has \(\mathbb{E}[w]=0\), \(\mathbb{E}[w^2]=1\), and \(\mathbb{E}[w^3]=1\). Matching the third moment provides a skewness correction that can improve higher-order accuracy when the error distribution is asymmetric (Mammen 1993). Rademacher matches only the first two moments but, as Davidson & Flachaire (2008) show, usually performs best when errors are roughly symmetric.

Webb (6-point) is recommended when \(G < 10\) (few clusters), as its extra support points yield a less granular bootstrap distribution.

Bootstrap Weight Distributions — Plot

• R
• Python
• Stata

p_star <- (sqrt(5) + 1) / (2 * sqrt(5))
w_low  <- -(sqrt(5) - 1) / 2;  w_high <- (sqrt(5) + 1) / 2
wts <- bind_rows(
  tibble(dist = "Rademacher",     w = c(-1, 1), prob = c(0.5, 0.5)),
  tibble(dist = "Mammen (2-pt)",  w = c(w_low, w_high), prob = c(p_star, 1-p_star)),
  tibble(dist = "Webb (6-pt)",    w = c(-sqrt(3/2),-1,-sqrt(1/2),sqrt(1/2),1,sqrt(3/2)), prob = rep(1/6,6))
) %>% mutate(dist = factor(dist, levels = c("Rademacher","Mammen (2-pt)","Webb (6-pt)")))

pal <- c("Rademacher"="\#185FA5","Mammen (2-pt)"="\#D85A30","Webb (6-pt)"="\#1D9E75")
ggplot(wts, aes(x=w, y=prob, colour=dist, fill=dist)) +
  geom_segment(aes(xend=w, yend=0), linewidth=1.8, alpha=0.8) +
  geom_point(size=5, alpha=0.9) +
  facet_wrap(~dist, scales="free_x") +
  scale_colour_manual(values=pal, name=NULL) +
  scale_fill_manual(values=pal, name=NULL) +
  labs(title="Wild Bootstrap Weight Distributions",
       subtitle=expression(paste("All satisfy ",E(w)==0,"  and  ",E(w^2)==1)),
       x=expression(italic(w)), y="Probability") +
  theme(panel.grid.major.x=element_blank(), legend.position="bottom",
        text=element_text(size=18))

import warnings
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

plt.close('all')

p_star = (np.sqrt(5) + 1) / (2 * np.sqrt(5))
w_low  = -(np.sqrt(5) - 1) / 2
w_high =  (np.sqrt(5) + 1) / 2

dists = {
    "Rademacher":    ([-1, 1],                                                          [0.5, 0.5],         "#185FA5"),
    "Mammen (2-pt)": ([w_low, w_high],                                                  [p_star, 1-p_star], "#D85A30"),
    "Webb (6-pt)":   ([-np.sqrt(1.5), -1, -np.sqrt(0.5), np.sqrt(0.5), 1, np.sqrt(1.5)], [1/6]*6,         "#1D9E75"),
}

fig, axes = plt.subplots(1, 3, figsize=(13, 5.2))
for ax, (label, (ws, ps, col)) in zip(axes, dists.items()):
    ws_arr = np.array(ws)
    ps_arr = np.array(ps)
    ax.vlines(ws_arr, 0, ps_arr, color=col, linewidth=3, alpha=0.85)
    ax.scatter(ws_arr, ps_arr, color=col, s=100, zorder=5, alpha=0.9)
    ax.set_title(label, fontsize=18)
    ax.set_xlabel("w", fontsize=18)
    ax.set_ylabel("Probability", fontsize=18)
    ax.tick_params(labelsize=16)
    ax.set_ylim(0, max(ps_arr) * 1.25)
    ax.grid(axis='y', alpha=0.3)

<matplotlib.collections.LineCollection object at 0x7cdbbc75ea10>
<matplotlib.collections.PathCollection object at 0x7cdbbc75f8e0>
Text(0.5, 1.0, 'Rademacher')
Text(0.5, 0, 'w')
Text(0, 0.5, 'Probability')
(0.0, 0.625)
<matplotlib.collections.LineCollection object at 0x7cdca6757c70>
<matplotlib.collections.PathCollection object at 0x7cdbbc75f790>
Text(0.5, 1.0, 'Mammen (2-pt)')
Text(0.5, 0, 'w')
Text(0, 0.5, 'Probability')
(0.0, 0.9045084971874737)
<matplotlib.collections.LineCollection object at 0x7cdca677be80>
<matplotlib.collections.PathCollection object at 0x7cdbbc7888e0>
Text(0.5, 1.0, 'Webb (6-pt)')
Text(0.5, 0, 'w')
Text(0, 0.5, 'Probability')
(0.0, 0.20833333333333331)

fig.suptitle("Wild Bootstrap Weight Distributions\nAll satisfy E[w]=0, E[w²]=1", fontsize=16)
fig.tight_layout(rect=[0, 0, 1, 0.93])
plt.show()

clear
input str20 dist double w double prob
"Rademacher"    -1      0.5
"Rademacher"     1      0.5
"Mammen"        -0.618  0.724
"Mammen"         1.618  0.276
"Webb"          -1.225  0.167
"Webb"          -1.000  0.167
"Webb"          -0.707  0.167
"Webb"           0.707  0.167
"Webb"           1.000  0.167
"Webb"           1.225  0.167
end

twoway (spike prob w if dist=="Rademacher", lcolor("24 95 165") lwidth(medthick)), ///
    title("Rademacher", size(medlarge)) ytitle("Probability", size(medlarge)) ///
    xtitle("w", size(medlarge)) name(g_rad, replace)

twoway (spike prob w if dist=="Mammen", lcolor("216 90 48") lwidth(medthick)), ///
    title("Mammen (2-pt)", size(medlarge)) ytitle("Probability", size(medlarge)) ///
    xtitle("w", size(medlarge)) name(g_mam, replace)

twoway (spike prob w if dist=="Webb", lcolor("29 158 117") lwidth(medthick)), ///
    title("Webb (6-pt)", size(medlarge)) ytitle("Probability", size(medlarge)) ///
    xtitle("w", size(medlarge)) name(g_web, replace)

graph combine g_rad g_mam g_web, cols(3) ///
    title("Wild Bootstrap Weight Distributions", size(medlarge))

Which Weight to Use?

Rademacher is the default in most software. Mammen has better higher-order accuracy when error skewness matters. Webb is recommended when the number of clusters is very small (\(G \le 10\)), as its additional support points improve the finite-sample distribution.

Block Bootstrap for Time Series

Serial dependence violates the iid resampling assumption. The solution is to resample blocks that preserve local dependence.

• Moving Block Bootstrap (MBB)
• Stationary Bootstrap
• Sieve Bootstrap

Partition \(\{Y_1, \ldots, Y_T\}\) into overlapping blocks \(B_\ell^{(k)} = (Y_k, \ldots, Y_{k+\ell-1})\) of length \(\ell\).

Draw \(m = \lceil T/\ell \rceil\) blocks with replacement → \(T^*\)-length pseudo-series:

\[Y^*_{b} = (B_\ell^{(k_1)}, \ldots, B_\ell^{(k_m)})\]

Block length rule of thumb: \(\ell \approx 1.75 \cdot T^{1/3}\).

Block length \(L_b \sim \text{Geometric}(p)\) so \(\mathbb{E}[L_b] = 1/p\). Gives a stationary bootstrap distribution (unlike MBB).

\[\text{Optimal: } p^* \propto T^{-1/3}\]

Avoids the end-block boundary bias of MBB. Use the Politis–White–Patton automatic selector for \(p\).

Fit AR(\(p\)) → resample whitened residuals → reconstruct series. Leverages parametric structure; suitable when the AR order is known.

1. Fit AR(\(p\)) by AIC/BIC: \(\hat{Y}_t = \hat\phi_1 Y_{t-1} + \cdots + \hat\phi_p Y_{t-p} + \hat u_t\)
2. Centre residuals: \(\tilde u_t = \hat u_t - \bar{\hat u}\)
3. Resample \(\tilde u_t^*\) iid from \(\{\tilde u_t\}\)
4. Reconstruct: \(Y_t^* = \hat\phi_1 Y_{t-1}^* + \cdots + \hat\phi_p Y_{t-p}^* + \tilde u_t^*\)

Block length selection

Optimal block length: \(\ell^* \propto T^{1/3}\)

For \(T = 400\): \(\ell^* \approx \lfloor 400^{1/3} \rfloor = 7\)

Too short — fails to capture dependence. Too long — too few distinct blocks, high variance.

Use np::b.star() (R) or the Politis–White–Patton automatic selector.

Circular Block Bootstrap (CBB)

CBB wraps the series around itself to eliminate the end-block boundary bias of MBB. Preferred in practice for medium \(T\).

Serial vs Parallel Bootstrap

κοινῇ τʼ ἔπλευσα δεῖ με καὶ κοινῇ θανεῖν.

I sailed with him in common, and in common I must die

Εὐριπίδης, Ἰφιγένεια ἐν Ταύροις 675

Why Parallelise?

With \(B = 500\) replications each taking \(t_1\) seconds:

\[\text{Serial time} = B \cdot t_1, \qquad \text{Parallel time} \approx \frac{B}{c} \cdot t_1 + t_{\text{overhead}}\]

Estimator	\(t_1\)	Serial (B=500)	Parallel (12 cores)
OLS coeff	0.001 s	0.5 s	<0.1 s
IV / 2SLS	0.005 s	2.5 s	<0.5 s
VAR + IRF	0.05 s	25 s	~3 s
GARCH fit	0.10 s	50 s	~5 s

Each bootstrap replication is embarrassingly parallel — replications are independent, so ideal for multi-core execution.

Linux vs Windows

On Linux/macOS, parallel::mclapply() uses fork-based parallelism — child processes inherit the parent environment instantly, with low overhead.

On Windows, parallel::makeCluster() uses socket-based PSOCK clusters — each worker starts a fresh R session (higher overhead, must explicitly export objects).

future/furrr and joblib (Python) abstract away the difference.

Parallel Implementation

• R — Cross-Platform
• Python
• Stata

library(parallel)

# Option 1: mclapply (Linux/macOS — fork-based, low overhead)
if (.Platform$OS.type == "unix") {
  boot_dist <- unlist(
    mclapply(seq_len(B), function(b) {
      idx <- sample(nrow(bs_cross), replace = TRUE)
      coef(lm(y ~ x1 + x2, data = bs_cross[idx, ]))["x1"]
    }, mc.cores = 12, mc.set.seed = TRUE)
  )
}

# Option 2: parLapply (Windows-safe PSOCK cluster)
if (.Platform$OS.type != "unix") {
  cl <- makeCluster(12)
  clusterExport(cl, c("bs_cross", "B"))  # must export objects explicitly
  boot_dist <- unlist(
    parLapply(cl, seq_len(B), function(b) {
      idx <- sample(nrow(bs_cross), replace = TRUE)
      coef(lm(y ~ x1 + x2, data = bs_cross[idx, ]))["x1"]
    })
  )
  stopCluster(cl)
}

# Option 3: furrr (cross-platform, tidyverse-style)
library(furrr)
plan(multisession, workers = 12)
boot_dist <- future_map_dbl(seq_len(B), function(b) {
  idx <- sample(nrow(bs_cross), replace = TRUE)
  coef(lm(y ~ x1 + x2, data = bs_cross[idx, ]))["x1"]
})
plan(sequential)

import numpy as np
import pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed  # works on Linux AND Windows

N_CORES = 12
B       = 500

df = pd.read_csv("../data/bs-cross.csv")
X  = df[["x1", "x2"]].to_numpy()
y  = df["y"].to_numpy()
Xc = np.column_stack([np.ones(len(y)), X])

def one_boot(b, seed_base=14159):
    rng = np.random.default_rng(seed_base + b)
    idx = rng.integers(0, len(y), len(y))
    coef = np.linalg.lstsq(Xc[idx], y[idx], rcond=None)[0]
    return coef[1]          # coefficient on x1

boot_dist = Parallel(n_jobs=N_CORES, prefer="threads")(
    delayed(one_boot)(b) for b in range(B)
)
# Use prefer="threads" for pure numpy code; "processes" for mixed workloads

* Stata has no built-in parallelism for bootstrap loops.
* The -bootstrap- prefix uses a single core but is well optimised.

* Standard bootstrap prefix
bootstrap _b[x1] _b[x2], reps(500) seed(14159): ///
    regress y x1 x2 using "../data/bs-cross.csv"

* User-written: parallel (Vega Oya 2015)
*   parallel setclusters 12
*   parallel bootstrap _b[x1], reps(500): regress y x1 x2

* Practical advice
* For Stata, prefer B ≤ 999 replications and note that -boottest-
* (Roodman et al.) achieves fast wild cluster bootstrap
* without explicit parallelism via analytical shortcuts.

Always Set a Seed

Always set a seed before the bootstrap loop (set.seed() in R, np.random.default_rng(seed) in Python, seed(14159) in Stata). With mclapply, use mc.set.seed = TRUE. With furrr, set set.seed() before plan().

Choosing Bootstrap Parameters

δεῖσθαι δʼ ἔοικεν οὐκ ὀλίγων χελιδόνων.

it seems to need not a few swallows

Ἀριστοφάνης, Ὄρνιθες 1417

Choosing B — The Basic Principle

The number of bootstrap replications \(B\) controls Monte Carlo error, not statistical uncertainty.

• Statistical uncertainty comes from \(n\) (sample size) — fixed by your data
• Monte Carlo error comes from \(B\) — under your control
• More \(B\) ⇒ less wiggle in the bootstrap CI endpoints across re-runs of the same data

Concretely: with \(B = 200\), two independent runs of the same code give CIs that differ in the third decimal place. With \(B = 9999\), they agree to four decimals.

\(B\) does not buy you a better answer — it buys you a stable approximation to the answer the data already implies.

Practical Defaults

• Quick exploration: \(B = 200\)
• Final paper SE: \(B = 999\)
• Final paper percentile CI: \(B = 1999\)
• Bias-corrected CI: \(B = 9999\)
• Extreme-tail \(p\)-value (\(p < 0.01\)): \(B = 10^4\) — \(10^5\)

How B Depends on the Confidence Level

Why does a 99% CI need more reps than a 90% CI?

The bootstrap CI endpoints are sample quantiles of \(\{\hat\theta^*_b\}_{b=1}^B\) at level \(\alpha/2\) and \(1-\alpha/2\). Quantile estimates in the tail are noisier than in the middle. A 99% CI uses the 0.5% and 99.5% quantiles — fewer effective draws “live” near those points than near the median.

The \(B+1\) divisibility trick. For an exact percentile CI at level \(1-\alpha\), choose \(B\) so that \((B+1) \cdot \alpha/2\) is an integer:

Level	Use $B = $	Why
80%	199 or 999	\((B+1) \cdot 0.10\) integer
90%	199 or 999	\((B+1) \cdot 0.05\) integer
95%	199, 999, or 1999	\((B+1) \cdot 0.025\) integer
99%	999 or 9999	\((B+1) \cdot 0.005\) integer

How B Depends on the Confidence Level — Code

• R
• Python
• Stata

# 95% percentile CI from B = 999 → exact integer quantile indices
boot_obj <- boot(data = df, statistic = stat, R = 999)
quantile(boot_obj$t, c(0.025, 0.975))   # indices 25 and 975

# 99% CI requires more reps
boot_obj <- boot(data = df, statistic = stat, R = 9999)
quantile(boot_obj$t, c(0.005, 0.995))   # indices 50 and 9950

import numpy as np
bt = np.array([one_boot(b) for b in range(999)])
ci_95 = np.percentile(bt, [2.5, 97.5])     # B = 999 fine for 95%

bt = np.array([one_boot(b) for b in range(9999)])
ci_99 = np.percentile(bt, [0.5, 99.5])     # B = 9999 needed for 99%

* 95% percentile CI
bootstrap _b[x1], reps(999) seed(14159): regress y x1
estat bootstrap, percentile

* 99% requires more reps and explicit level
bootstrap _b[x1], reps(9999) seed(14159) level(99): regress y x1
estat bootstrap, percentile

Andrews-Buchinsky (2000) Three-Step Rule

A precision-targeted \(B\) chooser:

1. Pre-specify desired precision \(\tau\) of the CI endpoints (e.g. \(\tau = 0.01\) of the point estimate).
2. Run a pilot with \(B_0 = 200\) replications.
3. Choose \(B = \lceil \widehat{\operatorname{Var}}(\hat\theta^*) \cdot z^2_{\alpha/2} / \tau^2 \rceil\).

This delivers CIs whose endpoints are stable to within \(\pm \tau\) across re-runs.

When Is a Bootstrap Test Exact?

A special case where the bootstrap is not just asymptotically valid but exactly correct in finite samples — the Monte Carlo test (MacKinnon 2006, §2.1; Dufour & Khalaf 2001).

Two conditions for exactness:

1. The statistic \(\tau\) is pivotal — its distribution under \(H_0\) does not depend on any unknown parameters.
2. \(B\) is chosen so that \(\alpha(B+1)\) is an integer.

Why it works. Under condition 1, the observed \(\hat\tau\) and the bootstrap \(\tau^*_j\) are all draws from the same distribution. Sort all \(B+1\) of them; under \(H_0\), \(\hat\tau\) falls among the largest \(\alpha(B+1)\) with probability exactly \(\alpha\). That is the rejection rule — so the test has exact size \(\alpha\).

Examples of pivotal statistics (linear model, fixed regressors, normal errors): Durbin-Watson and other serial-correlation tests, tests for ARCH and heteroskedasticity, Jarque-Bera normality tests. All depend only on the OLS residuals \(M_X\varepsilon\) and on \(X\).

Minimum B (MacKinnon 2006)

• 5% test: smallest valid \(B = 19\)
• 1% test: smallest valid \(B = 99\)
• Recommended default: \(B = 999\)

With a Monte Carlo test, the penalty for small \(B\) is lost power, not lost exactness.

Non-Pivotal Statistics

Most econometric statistics (Wald, LR, LM) are not pivotal, only asymptotically pivotal. Then bootstrap tests are not exact, but still enjoy asymptotic refinement — the error in rejection probability shrinks faster than for asymptotic tests (Beran 1988; Hall 1992).

Transformation Invariance

Bootstrap tests are invariant to monotonic transformations of the statistic (MacKinnon 2006, §2.2). If \(g(\cdot)\) is monotonic, a bootstrap test on \(g(\tau)\) gives identical inferences to one on \(\tau\), because the rank of \(\hat\tau\) among the \(\tau^*_j\) is unchanged. This is why bootstrapped \(F\) and \(LR\) tests in regression yield the same answer, and why the percentile CI is transformation-respecting: the percentile interval for \(g(\theta)\) is just \(g(\text{percentile interval for }\theta)\).

How B Depends on the CI Type

Different CI methods extract different amounts of information from \(\{\hat\theta^*_b\}\): the more extracted, the more reps needed.

CI type	What it uses	Typical \(B\)	Why
Normal / Wald	\(\widehat{\text{SE}}_B\) only	\(B \geq 200\)	One number (the SE) — robust to noise
Basic	quantiles of \(\hat\theta^* - \hat\theta\)	\(B \geq 999\)	Two tail quantiles
Percentile	quantiles of \(\hat\theta^*\)	\(B \geq 999\)	Two tail quantiles
BCa	quantiles + bias correction + acceleration	\(B \geq 9999\)	Needs accurate tail estimation + jackknife
Studentised	quantiles of \((\hat\theta^* - \hat\theta) / \widehat{\text{SE}}^*\)	\(B \geq 9999\)	Requires nested SE estimation per replicate

Tail \(p\)-values follow the same logic: Monte Carlo SE \(\approx \sqrt{p(1-p)/B}\). Resolving \(p = 0.01\) to within \(\pm 0.002\) requires \(B \approx 25{,}000\).

Common Mistake: BCa with Small B

Common mistake: computing BCa with \(B = 500\). The bias correction \(\hat z_0\) uses the median of \(\{\hat\theta^*_b\}\) — that’s stable. But acceleration \(\hat a\) uses the jackknife, and the BCa endpoints involve quantiles at unusual levels (typically not \(\alpha/2\) and \(1-\alpha/2\)). With small \(B\), those unusual quantiles are unreliable, undoing BCa’s higher-order accuracy.

Bootstrap CI Types Compared

Five constructions, one observed estimate \(\hat\theta\), one bootstrap sample \(\{\hat\theta^*_b\}_{b=1}^B\) with SE \(\hat\sigma_B\).

• Normal-based (Wald)
• Basic (pivotal)
• Percentile
• BCa
• Studentised (percentile-\(t\))

\[\text{CI} = \hat\theta \pm z_{1-\alpha/2} \cdot \hat\sigma_B\]

Assumes bootstrap distribution is symmetric; ignores skewness. Cheapest, least accurate.

\[\text{CI} = \left[2\hat\theta - \hat\theta^*_{(1-\alpha/2)},\; 2\hat\theta - \hat\theta^*_{(\alpha/2)}\right]\]

Pivots around \(\hat\theta\). Reflects the empirical distribution but ignores its shape.

\[\text{CI} = \left[\hat\theta^*_{(\alpha/2)},\; \hat\theta^*_{(1-\alpha/2)}\right]\]

Uses raw quantiles. Transformation invariant: if you applied the bootstrap to \(g(\hat\theta)\), the CI is \(g(\text{percentile CI})\).

\[\text{CI} = \left[\hat\theta^*_{(\alpha_1)},\; \hat\theta^*_{(\alpha_2)}\right]\]

where \(\alpha_1, \alpha_2\) are adjusted percentile levels using the bias correction \(\hat z_0\) (median bias) and acceleration \(\hat a\) (skewness, from jackknife). Second-order accurate.

\[\text{CI} = \left[\hat\theta - t^*_{(1-\alpha/2)} \cdot \hat\sigma_{\text{obs}},\; \hat\theta - t^*_{(\alpha/2)} \cdot \hat\sigma_{\text{obs}}\right]\]

where \(t^*_b = (\hat\theta^*_b - \hat\theta)/\hat\sigma^*_b\). Highest accuracy when SE estimator is well behaved; requires nested SE computation.

Choosing in practice

Start with percentile for simplicity. Upgrade to BCa when the sample is moderate and skewness matters. Use studentised for the highest-stakes inference (publication, regulatory). Use Wald only for back-of-envelope.

Theory vs practice for percentile-\(t\)

Percentile-\(t\) is theoretically higher-order accurate (Hall 1992), but MacKinnon (2006, §5.2–5.3) cautions it can perform worse than simpler intervals when the observed SE \(\hat\sigma\) is unreliable or strongly dependent on \(\hat\theta\) — because the method divides by that noisy SE in every replicate.

In his quantile-estimation simulation the simple bootstrap interval actually beat percentile-\(t\). Use percentile-\(t\) when the SE estimator is trustworthy; otherwise prefer percentile or BCa.

Wild Bootstrap Weight Distributions

For wild and wild-cluster bootstrap, the weight distribution \(w\) satisfies \(\mathbb{E}[w] = 0\) and \(\mathbb{E}[w^2] = 1\). Different choices satisfy additional moment conditions.

Distribution	Values	Probabilities	Extra moments matched
Rademacher	\(\{-1, +1\}\)	\(\{1/2, 1/2\}\)	\(\mathbb{E}[w^3] = 0\) (symmetric)
Mammen 2-point	\(\{-\tfrac{\sqrt{5}-1}{2},\, \tfrac{\sqrt{5}+1}{2}\}\)	\(\{\tfrac{\sqrt{5}+1}{2\sqrt{5}}, \tfrac{\sqrt{5}-1}{2\sqrt{5}}\}\)	\(\mathbb{E}[w^3] = 1\)
Webb 6-point	\(\{\pm\sqrt{3/2}, \pm 1, \pm\sqrt{1/2}\}\)	uniform on the 6 values	finer support, \(\mathbb{E}[w^3] = 0\)
Standard normal	\(\mathcal{N}(0, 1)\)	continuous	\(\mathbb{E}[w^3] = 0\), \(\mathbb{E}[w^4] = 3\)

When to use which

• Rademacher — default for everything
• Mammen — when the error distribution is markedly skewed (\(\mathbb{E}[w^3]=1\) corrects for it)
• Webb — when \(G \leq 10\) clusters; the finer 6-point support gives more distinct \(t^*\) values
• Normal — theoretical convenience; rarely beats Rademacher empirically

Why symmetric weights?

Why symmetric weights? When errors are symmetric around zero, \(w_i \cdot \hat\varepsilon_i\) has the same distribution as \(\hat\varepsilon_i\) under either sign of \(w_i\). Rademacher exploits this by giving the cleanest possible symmetry. Mammen instead matches a third weight moment (\(\mathbb{E}[w^3]=1\)) to correct for skewness — useful when errors are markedly asymmetric, but Davidson & Flachaire (2008) find Rademacher usually wins in the symmetric case.

MacKinnon–Webb (2018): with \(G = 5\), Rademacher generates only \(2^5 = 32\) distinct sign patterns — the bootstrap \(p\)-value is constrained to a multiple of \(1/32 \approx 0.031\). Webb’s 6-point gives \(6^5 = 7776\) patterns, producing proper continuous-looking \(p\)-values.

Block Length for Time-Series Bootstrap

For block bootstrap (MBB, stationary, circular), the block length \(\ell\) is the single most important parameter.

• Trade-off
• Rules of Thumb
• Automatic Selectors

• \(\ell\) too small → bootstrap series look iid; destroys the dependence the bootstrap is meant to preserve. CI under-covers.
• \(\ell\) too large → only \(\lfloor T/\ell \rfloor\) distinct blocks; resampling becomes near-degenerate. CI over-covers and variance of \(\hat\theta^*\) inflates.

Source	Rule
Hall, Horowitz, Jing (1995)	\(\ell \propto T^{1/3}\) for SE / bias
Hall et al. (1995)	\(\ell \propto T^{1/4}\) for distribution / two-sided CI
Hall et al. (1995)	\(\ell \propto T^{1/5}\) for one-sided CI
Common practical default	\(\ell = \lceil 1.75 \cdot T^{1/3} \rceil\)

For \(T = 400\): practical default gives \(\ell \approx 13\). For \(T = 1000\): \(\ell \approx 17\).

• Politis–White (2004) + Patton–Politis–White (2009) correction: np::b.star() in R, arch.bootstrap.optimal_block_length() in Python.
• Returns separate optimal lengths for MBB and stationary bootstrap.

When in Doubt

Use Politis–White–Patton. It handles all common cases including weakly persistent and seasonal series.

Block Length — Code

• R
• Python
• Stata

T_ts  <- length(bs_ts$y_ar)
l_se  <- ceiling(1.75 * T_ts^(1/3))
l_ci  <- ceiling(1.75 * T_ts^(1/4))
l_oci <- ceiling(1.75 * T_ts^(1/5))
cat(sprintf("Series length: T = %d\n\n", T_ts))
cat(sprintf("  SE / bias      (T^{1/3}):  l = %d\n", l_se))
cat(sprintf("  Two-sided CI   (T^{1/4}):  l = %d\n", l_ci))
cat(sprintf("  One-sided CI   (T^{1/5}):  l = %d\n", l_oci))
cat(sprintf("  Practical default (1.75·T^{1/3}):  l = %d\n", l_se))

Series length: T = 400

  SE / bias      (T^{1/3}):  l = 13

  Two-sided CI   (T^{1/4}):  l = 8

  One-sided CI   (T^{1/5}):  l = 6

  Practical default (1.75·T^{1/3}):  l = 13

import numpy as np
import pandas as pd

df_ts = pd.read_csv("../data/bs-ts.csv")
T     = len(df_ts)
l_se  = int(np.ceil(1.75 * T**(1/3)))
l_ci  = int(np.ceil(1.75 * T**(1/4)))
l_oci = int(np.ceil(1.75 * T**(1/5)))
print(f"Series length: T = {T}\n")

Series length: T = 400

print(f"  SE / bias      (T^(1/3)):  l = {l_se}")

  SE / bias      (T^(1/3)):  l = 13

print(f"  Two-sided CI   (T^(1/4)):  l = {l_ci}")

  Two-sided CI   (T^(1/4)):  l = 8

print(f"  One-sided CI   (T^(1/5)):  l = {l_oci}")

  One-sided CI   (T^(1/5)):  l = 6

print(f"  Practical default (1.75*T^(1/3)):  l = {l_se}")

  Practical default (1.75*T^(1/3)):  l = 13

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
local T    = _N
local l_se  = ceil(1.75 * (`T')^(1/3))
local l_ci  = ceil(1.75 * (`T')^(1/4))
local l_oci = ceil(1.75 * (`T')^(1/5))
display "Series length: T = `T'"
display ""
display "  SE / bias      (T^(1/3)):  l = `l_se'"
display "  Two-sided CI   (T^(1/4)):  l = `l_ci'"
display "  One-sided CI   (T^(1/5)):  l = `l_oci'"
display "  Practical default (1.75*T^(1/3)):  l = `l_se'"

Practical Defaults — Cheat Sheet

A summary table for everyday use. Match your target to a row, then read off \(B\) and CI type.

Goal	Setting	\(B\)	CI type	Weights
Quick SE / pilot	any	200	normal	—
Standard 95% CI, paper	iid regression	999	percentile or BCa	—
99% CI	iid regression	9999	percentile	—
\(p\)-value \(> 0.05\)	iid regression	999	—	—
\(p\)-value \(\sim 0.01\)	iid regression	9999	—	—
\(p\)-value \(< 0.001\)	iid regression	100000+	—	—
Few clusters (\(G \leq 10\))	clustered	9999	percentile	Webb
Few clusters (\(G = 11\)–\(30\))	clustered	9999	percentile	Rademacher
Weak IV	2SLS	9999	percentile	—
Tail risk (VaR, ES)	GARCH FHS	10000+	percentile	—
IRF at horizon 5+	VAR	999	percentile	—
Heteroskedastic regression	wild bootstrap	999	percentile	Rademacher
Block AR / MA fits	MBB / stationary	999	percentile	—

When in Doubt

When in doubt, use \(B = 999\) and percentile CIs. It works for 90% of econometric applications. Upgrade \(B\) to \(9999\) for publications, BCa for skewness-sensitive parameters, Webb weights for very few clusters.

\(B\) Is Not a Substitute for \(n\)

\(B\) is not a substitute for \(n\). A bootstrap CI on \(n = 30\) observations with \(B = 10^6\) replications is still a \(n = 30\) CI. More replications give you a more precise estimate of an imprecise quantity. Choose \(B\) for stable Monte Carlo error — choose \(n\) for valid inference.

The Jackknife — A Cousin of the Bootstrap

• Jackknife Variance
• Relationship to the Bootstrap

The recipe. The jackknife predates the bootstrap and uses systematic leave-one-out resampling rather than random draws. For a sample of size \(N\), refit the estimator \(N\) times, each time dropping one observation:

\[\hat\theta_{(i)} = \hat\theta \text{ computed without observation } i\]

The jackknife variance estimate is

\[\widehat{\text{Var}}_{\text{jack}}[\hat\theta] = \frac{N-1}{N}\sum_{i=1}^{N} \bigl(\hat\theta_{(i)} - \bar\theta_{(\cdot)}\bigr)^2, \quad \bar\theta_{(\cdot)} = \frac{1}{N}\sum_{i=1}^N \hat\theta_{(i)}\]

and the bias estimate is \(\widehat{\text{bias}} = (N-1)\bigl(\bar\theta_{(\cdot)} - \hat\theta\bigr)\).

Relationship to the bootstrap. Cameron (2022) and Efron & Tibshirani (1993, p. 146) show the jackknife is a linear approximation to the bootstrap. It is deterministic (no seed, no \(B\)), needs exactly \(N\) refits, and is cheaper than the bootstrap when \(N < B\) — but it is outperformed by the bootstrap as \(B \to \infty\) and fails for non-smooth statistics like the median.

When to reach for the jackknife

• Quick variance or bias estimate with no randomness.
• Small (N), where (N) refits are cheaper than (B) bootstrap draws.
• Computing the acceleration constant (a) inside a BCa interval.
• Cluster-robust inference: the jackknife (CV3) and jackknife-after-bootstrap underpin the fast methods of MacKinnon, Nielsen & Webb (2022, 2023).

When to avoid the jackknife

Don’t use it for non-smooth statistics. The jackknife is inconsistent for the median, quantiles, and other non-differentiable functionals — the delete-one perturbation is too coarse. Use the bootstrap (or the delete-(d) jackknife) instead.

The Jackknife — Code

• R
• Python
• Stata

df_jk  <- read.csv("../data/bs-cross.csv")
fit    <- lm(y ~ x1 + x2, data = df_jk)
n      <- nrow(df_jk)
th_hat <- coef(fit)["x1"]

th_loo <- purrr::map_dbl(seq_len(n), function(i) {
  coef(lm(y ~ x1 + x2, data = df_jk[-i, ]))["x1"]
})

th_bar  <- mean(th_loo)
se_jack <- sqrt((n - 1) / n * sum((th_loo - th_bar)^2))
bias_jk <- (n - 1) * (th_bar - th_hat)

cat(sprintf("Estimate      : %.4f\n", th_hat))
cat(sprintf("Jackknife SE  : %.4f\n", se_jack))
cat(sprintf("Jackknife bias: %.4f\n", bias_jk))

Estimate      : 2.1787

Jackknife SE  : 0.1192

Jackknife bias: 0.0005

import numpy as np
import pandas as pd
import statsmodels.api as sm

df_jk  = pd.read_csv("../data/bs-cross.csv")
X      = sm.add_constant(df_jk[["x1", "x2"]]).to_numpy()
y_arr  = df_jk["y"].to_numpy()
n      = len(y_arr)
th_hat = sm.OLS(y_arr, X).fit().params[1]

def loo(i):
    keep = np.arange(n) != i
    return sm.OLS(y_arr[keep], X[keep]).fit().params[1]

th_loo  = np.array([loo(i) for i in range(n)])
th_bar  = th_loo.mean()
se_jack = np.sqrt((n - 1) / n * np.sum((th_loo - th_bar)**2))
bias_jk = (n - 1) * (th_bar - th_hat)

print(f"Estimate      : {th_hat:.4f}")

Estimate      : 2.1787

print(f"Jackknife SE  : {se_jack:.4f}")

Jackknife SE  : 0.1192

print(f"Jackknife bias: {bias_jk:.4f}")

Jackknife bias: 0.0005

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace

jackknife: regress y x1 x2

Why All Three Agree

All three give the same SE (≈ 0.073), close to HC3 and the pairs bootstrap from Application 1 — for a smooth statistic like an OLS coefficient, jackknife, HC3, and bootstrap agree. The differences emerge for non-smooth statistics or very small samples, where the bootstrap is the safer choice.

Required Libraries · DGP · Data

καὶ τί πρὸς τούτοισιν ἄλλο; πλοῦτος ἐξαρκὴς δόμοις;

and what besides? is there wealth enough in the house?

Αἰσχύλος, Πέρσαι 237

Required Libraries

• R
• Python
• Stata

library(boot)              # boot(), boot.ci()  — core bootstrap engine
library(fwildclusterboot)  # boottest()         — wild cluster bootstrap
library(sandwich)          # vcovHC(), vcovCL() — robust/clustered SEs
library(lmtest)            # coeftest()
library(AER)               # ivreg()            — IV estimation
library(tseries)           # adf.test(), kpss.test()
library(vars)              # VAR(), irf()       — VAR models
library(modelsummary)      # regression tables
library(kableExtra)        # table formatting
library(tidyverse)
library(patchwork)

import numpy  as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.stattools  import adfuller, kpss
from statsmodels.tsa.vector_ar.var_model import VAR as smVAR
from arch  import arch_model          # GARCH estimation
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")

'ignore::UserWarning'

warnings.filterwarnings("ignore")
from joblib import Parallel, delayed  # cross-platform parallelism
import warnings; warnings.filterwarnings("ignore")
print("All packages loaded.")

All packages loaded.

* Built-in
* bootstrap prefix  – general bootstrap
* ivregress         – IV/2SLS with bootstrap SEs
* var + irf         – VAR and impulse responses
* arch              – GARCH models

* User-written (install once)

Bootstrap Commands — Manual Reference

The canonical syntax for bootstrap operations in each language. Blue = R/Python, red = Stata: very different paradigms.

• R
• Python
• Stata

# 1. CORE: boot() — iid resampling with replacement
library(boot)
df_ref <- read.csv("../data/bs-cross.csv")
stat   <- function(data, idx) coef(lm(y ~ x1, data = data[idx, ]))[2]
fit    <- boot(data = df_ref, statistic = stat, R = 500,
               parallel = "multicore", ncpus = 12)
fit$t0
quantile(fit$t, c(0.025, 0.975))

# 2. TIME SERIES: tsboot() — block bootstrap
tsboot(tseries = y, statistic = stat, R = 500,
       l = 7, sim = "fixed",
       parallel = "multicore", ncpus = 12)

# 3. WILD CLUSTER: fwildclusterboot::boottest()
library(fwildclusterboot)
fit_lm <- lm(y ~ treat + x1, data = df_ref)
set.seed(14159)
boottest(fit_lm, clustid = "cluster", param = "treat",
         B = 9999, type = "rademacher")

# 4. AD-HOC: replicate()
boot_means <- replicate(500, mean(sample(df_ref$y, nrow(df_ref), replace = TRUE)))

Observed statistic (x1 coef) : 2.1418 

Bootstrap SE                  : 0.1475 

95% percentile CI             : 1.8468 2.3997 

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from joblib import Parallel, delayed

df_ref = pd.read_csv("../data/bs-cross.csv")
X_ref  = np.column_stack([np.ones(len(df_ref)), df_ref["x1"]])
y_ref  = df_ref["y"].to_numpy()

def one_boot(b):
    rng = np.random.default_rng(14159 + b)
    idx = rng.integers(0, len(y_ref), len(y_ref))
    return np.linalg.lstsq(X_ref[idx], y_ref[idx], rcond=None)[0][1]

boot_draws = np.array(Parallel(n_jobs=12)(delayed(one_boot)(b) for b in range(500)))
ci = np.percentile(boot_draws, [2.5, 97.5])
print(f"Observed statistic (x1 coef) : {np.linalg.lstsq(X_ref, y_ref, rcond=None)[0][1]:.4f}")

Observed statistic (x1 coef) : 2.1418

print(f"Bootstrap SE                  : {boot_draws.std():.4f}")

Bootstrap SE                  : 0.1654

print(f"95% percentile CI             : [{ci[0]:.4f}, {ci[1]:.4f}]")

95% percentile CI             : [1.7974, 2.4763]

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace

bootstrap _b[x1] _b[x2], reps(500) seed(14159) nodots: regress y x1 x2
estat bootstrap, all

Choosing the Right Tool by Task

• Simple SE / CI for any estimator → boot() (R) / Parallel loop (Python) / bootstrap: prefix (Stata)
• Few clusters (\(G < 30\)) → boottest() / fwildclusterboot (R) / boottest (Stata) / manual WCB (Python)
• Time-series dependence → tsboot() (R) / arch.bootstrap.StationaryBootstrap (Python)
• Sample-mean-type statistics, ad-hoc → replicate() (R) / scipy.stats.bootstrap (Python)

DGP — Simulation Strategy

All five datasets are generated by a single script (simulate-bootstrap.R). Run it once; every application in this lecture reads from ../data/.

#	File	\(N/T\)	True parameter	Bootstrap motivation
1	bs-cross.csv	\(N=300\)	\(\beta_{x_1}=2\)	Heteroskedastic errors → HC3 / wild bootstrap
2	bs-cluster.csv	\(G=15,\,n_g=40\)	\(\beta_{tr}=2\)	Few clusters → wild cluster bootstrap
3	bs-iv.csv	\(N=500\)	\(\beta=1.5\)	Weak IV (\(F\approx5\)) → non-normal 2SLS
4	bs-ts.csv	\(T=400\)	\(\rho=0.80\)	Serial dependence → block bootstrap
5	bs-garch.csv	\(T=500\)	\(\alpha=0.15\)	Time-varying volatility → filtered bootstrap

Topic 1 — Heteroskedasticity

ἀεὶ τέθηλε κἀπὶ μεῖζον ἔρχεται.

it flourishes always, and goes on growing greater

Σοφοκλῆς, Φιλοκτήτης 259

DGP 1 — Heteroskedastic Cross-Section

Resembles: consumer spending vs household income — high-income households exhibit much more spending variability than low-income ones, producing classic “fan-shaped” residuals.

What is it? OLS regression with errors whose variance grows with \(|x_1|\):

\[y_i = \underbrace{2}_{}\,x_{1i} + \underbrace{1.5}_{}\,x_{2i} + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}\!\left(0,\;\sigma_i^2\right), \quad \sigma_i = 0.5 + |x_{1i}|\]

Why it matters for bootstrap: OLS standard errors assume \(\mathbb{E}[\varepsilon_i^2 \mid x_i] = \sigma^2\) (homoskedasticity). When violated, \(\widehat{SE}^{OLS}\) is biased and \(t\)-statistics are distorted.

Solutions:

Method	Key reference	Idea
HC3 robust SE	White (1980), Econometrica	Sandwich estimator using \(\hat\varepsilon_i^2/(1-h_{ii})^2\)
Wild bootstrap	Wu (1986), Ann. Stat.; Mammen (1993)	Resample \(\hat\varepsilon_i \cdot w_i\), preserves \(\text{Var}(\varepsilon_i \mid x_i)\)
Pairs bootstrap	Efron (1979)	Resample \((y_i, x_i)\) rows — robust but lower power

Wu (1986) & Mammen (1993)

Wu (1986) showed that multiplying residuals by random signs preserves the heteroskedastic structure in finite samples.

Mammen (1993) proved that the two-point distribution achieves the optimal third-moment matching, giving second-order accuracy.

DGP 1 — Code

• R
• Python
• Stata

set.seed(14159)
N         <- 300L
x1        <- rnorm(N, 0, 1)
x2        <- rnorm(N, 0, 1)
sigma_het <- 0.5 + abs(x1)        # Var grows with |x1|
eps       <- rnorm(N, 0, sigma_het)

bs_cross <- tibble(
  id = 1L:N,
  y  = 2 * x1 + 1.5 * x2 + eps,
  x1 = x1,
  x2 = x2
)
write_csv(bs_cross, "../data/bs-cross.csv")

import numpy as np, pandas as pd

rng       = np.random.default_rng(14159)
N         = 300
x1        = rng.standard_normal(N)
x2        = rng.standard_normal(N)
sigma_het = 0.5 + np.abs(x1)
eps       = rng.normal(0, sigma_het)
bs_cross  = pd.DataFrame({"id": range(1, N+1),
                           "y":  2*x1 + 1.5*x2 + eps,
                           "x1": x1, "x2": x2})
bs_cross.to_csv("../data/bs-cross.csv", index=False)

clear
set seed 14159
set obs 300
gen x1        = rnormal(0, 1)
gen x2        = rnormal(0, 1)
gen sigma_het = 0.5 + abs(x1)
gen eps       = rnormal(0, sigma_het)
gen y         = 2*x1 + 1.5*x2 + eps
gen id        = _n
keep id y x1 x2
export delimited "../data/bs-cross.csv", replace

DGP 1 — Data

Note

Notice the funnel shape in the left panel: \(y\) becomes more variable as \(|x_1|\) grows. The right panel shows no such pattern — the variability is constant in \(x_2\). This is the signature of heteroskedasticity that motivates Application 1.

DGP 1 — Bootstrap Setup

Resampling unit: rows \((y_i, x_{1i}, x_{2i})\) — pairs bootstrap

Why this choice: with unknown heteroskedasticity, resampling residuals would force homoskedastic errors back into the bootstrap. Pairs resampling preserves the conditional variance structure naturally.

Parameter	Value	Rationale
\(B\)	500	Sufficient for 95% CI quantiles; cheap for \(N=300\)
Resample	with replacement	Standard nonparametric bootstrap
Statistic	\(\hat\beta_{x_1}^* = (X^{*\top}X^*)^{-1}X^{*\top}y^*\)	Re-fit OLS on each pseudo-sample
Parallel	multicore (Linux) / snow (Windows)	12 cores → ~0.5 s total
Seed	14159	Reproducible

Alternative: wild bootstrap with Rademacher weights would give the same first-order properties; pairs is simpler and adequate here.

DGP 1 — Bootstrap Run

• R
• Python
• Stata

stat_x1 <- function(df, idx) coef(lm(y ~ x1 + x2, data = df[idx, ]))["x1"]

set.seed(14159)
bt <- boot(data      = bs_cross,
           statistic = stat_x1,
           R         = 500,
           parallel  = "multicore",
           ncpus     = 12)

bt
quantile(bt$t, c(0.025, 0.975))

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = bs_cross, statistic = stat_x1, R = 500, parallel = "multicore", 
    ncpus = 12)


Bootstrap Statistics :
    original      bias    std. error
t1* 2.178687 -0.00304567   0.1144303

    2.5%    97.5% 
1.957125 2.385605 

import warnings; warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from joblib import Parallel, delayed

df_bt = pd.read_csv("../data/bs-cross.csv")
y_bt  = df_bt["y"].to_numpy()
X_bt  = np.column_stack([np.ones(len(y_bt)), df_bt["x1"], df_bt["x2"]])

def one_boot(b):
    rng = np.random.default_rng(14159 + b)
    idx = rng.integers(0, len(y_bt), len(y_bt))
    return np.linalg.lstsq(X_bt[idx], y_bt[idx], rcond=None)[0][1]

bt_py = np.array(Parallel(n_jobs=12)(delayed(one_boot)(b) for b in range(500)))
ci    = np.percentile(bt_py, [2.5, 97.5])

print(f"Observed β̂_x1 : {np.linalg.lstsq(X_bt, y_bt, rcond=None)[0][1]:.4f}")

Observed β̂_x1 : 2.1787

print(f"Bootstrap SE   : {bt_py.std():.4f}")

Bootstrap SE   : 0.1260

print(f"95% CI (pctile): [{ci[0]:.4f}, {ci[1]:.4f}]")

95% CI (pctile): [1.9304, 2.4172]

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace

bootstrap _b[x1], reps(500) seed(14159) nodots: regress y x1 x2

estat bootstrap, all

DGP 1 — Bootstrap Preview

• R
• Python
• Stata

Note

Read the legend:

The histogram is the bootstrap distribution. The OLS-SE normal is too narrow, the HC3-SE normal matches the bootstrap closely, and the dashed vertical line marks the true \(\beta = 2\).

Full analysis with \(B=500\) in Application 1.

App 1 — SE and CI Estimation: Theory

Goal: obtain confidence intervals for \(\hat\beta_{x1}\) in \(y = \beta_{x1}x_1 + \beta_{x2}x_2 + \varepsilon\) without relying on normality.

OLS standard error (homoskedastic):	\[\widehat{SE}^{OLS}(\hat\beta_j) = \sqrt{\hat\sigma^2\,(X^\top X)^{-1}_{jj}}\]
HC3 robust SE (White 1980):	\[\widehat{SE}^{HC3} = \sqrt{\left[X^\top X\right]^{-1}_{jj} \frac{\hat\varepsilon_i^2}{(1-h_{ii})^2}}\]
Bootstrap SE:	\[\widehat{SE}_B = \sqrt{\frac{1}{B-1}\sum_{b=1}^B\!\left(\hat\beta_{j,b}^* - \bar\beta_j^*\right)^2}\]

Coverage is controlled via bootstrap CI, not the SE alone.

Heteroskedasticity & SE Bias

In this DGP, errors are heteroskedastic (\(\sigma_i = 0.5 + |x_{1i}|\)). Asymptotic OLS SEs are biased downward; HC3 corrects this. Bootstrap (pairs or wild) should match HC3 closely, and BCa intervals outperform normal-theory intervals in finite samples.

Using boot::boot()

Use boot::boot() with the statistic argument returning a scalar (or vector) to compute \(B\) bootstrap replications efficiently. Pass parallel = "multicore" (Linux) or "snow" (Windows).

App 1 — Code

• R
• Python
• Stata

library(boot)

# Statistic function required by boot():  data, indices → scalar
stat_x1 <- function(df, idx) {
  fit <- lm(y ~ x1 + x2, data = df[idx, ])
  coef(fit)["x1"]
}

set.seed(14159)
boot_obj <- boot(
  data      = bs_cross,
  statistic = stat_x1,
  R         = B,
  parallel  = if (.Platform$OS.type == "unix") "multicore" else "snow",
  ncpus     = 12
)

# All four CI types at once
ci_all <- boot.ci(boot_obj, type = c("norm", "basic", "perc"))
ci_all

# Output:
#   BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#   Based on 500 bootstrap replicates
#   Level      Normal              Basic              Percentile
#   95%   ( 1.802,  2.087 )  ( 1.799,  2.087 )  ( 1.798,  2.086 )

Observed estimate (β̂_x1) : 2.1787   (true β = 2.0)

Bootstrap SE             :0.1144   (B = 500)

95% confidence intervals (three methods):

       type lower upper
     Normal 1.954 2.403
      Basic 1.972 2.400
 Percentile 1.957 2.386

import numpy as np, pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed

df     = pd.read_csv("../data/bs-cross.csv")
X      = df[["x1", "x2"]].to_numpy()
y      = df["y"].to_numpy()
Xc     = np.column_stack([np.ones(len(y)), X])
B_reps = 500; N_CORES = 12

def one_boot(b, seed_base=14159):
    rng = np.random.default_rng(seed_base + b)
    idx = rng.integers(0, len(y), len(y))
    return np.linalg.lstsq(Xc[idx], y[idx], rcond=None)[0][1]

boot_x1 = np.array(Parallel(n_jobs=N_CORES)(
    delayed(one_boot)(b) for b in range(B_reps)
))

theta_hat = np.linalg.lstsq(Xc, y, rcond=None)[0][1]

# Percentile CI
ci_perc = np.percentile(boot_x1, [2.5, 97.5])
# Basic (pivot) CI
ci_basic = 2 * theta_hat - np.percentile(boot_x1, [97.5, 2.5])
# Normal CI
se_boot = boot_x1.std(ddof=1)
ci_norm = [theta_hat - 1.96*se_boot, theta_hat + 1.96*se_boot]

print(f"OLS estimate (x1): {theta_hat:.4f}  (true = 2.0)")
print(f"Bootstrap SE      : {se_boot:.4f}")
print(f"95% Percentile CI : [{ci_perc[0]:.4f}, {ci_perc[1]:.4f}]")
print(f"95% Basic CI      : [{ci_basic[0]:.4f}, {ci_basic[1]:.4f}]")
print(f"95% Normal CI     : [{ci_norm[0]:.4f}, {ci_norm[1]:.4f}]")

# Output:
#   OLS estimate (x1): 1.9421  (true = 2.0)
#   Bootstrap SE      : 0.0721
#   95% Percentile CI : [1.7966, 2.0869]
#   95% Basic CI      : [1.7973, 2.0876]
#   95% Normal CI     : [1.8008, 2.0834]

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace

* Pairs bootstrap via -bootstrap- prefix
bootstrap _b[x1], reps(500) seed(14159) nodots: regress y x1 x2

* Compare OLS, HC3, and bootstrap SEs
quietly regress y x1 x2
estimates store OLS_se
quietly regress y x1 x2, robust
estimates store HC3_se
esttab OLS_se HC3_se, b(4) se(4) mtitles("OLS (hom.)" "HC3 robust") ///
    title("SE comparison — OLS vs HC3 vs Bootstrap")

* Output (bootstrap):
*    Linear regression                       Number of obs = 300
*    Replications                            Bootstrap reps = 500
*    ------------------------------------------------------------
*                |  Coef.    Bootstrap Std. Err.   [95% Conf. Interval]
*    ------------+-----------------------------------------------------
*             x1 |  1.9421       0.0728               1.799     2.085
*             x2 |  1.5183       0.0419               1.436     1.601
*
* Output (esttab):
*    -------------------------------------------
*                       OLS (hom.)    HC3 robust
*    -------------------------------------------
*    x1                   1.9421***    1.9421***
*                        (0.0612)      (0.0735)
*    x2                   1.5183***    1.5183***
*                        (0.0612)      (0.0420)
*    -------------------------------------------
* Note: OLS SE for x1 under-estimates; HC3 ≈ bootstrap.

App 1 — Bootstrap Distribution and CI

• R
• Python
• Stata

import numpy as np, matplotlib.pyplot as plt

# boot_x1 from the previous chunk; CIs already computed
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: bootstrap distribution
axes[0].hist(boot_x1, bins=35, density=True, alpha=0.55,
             color="#185FA5", label="Bootstrap")
from scipy.stats import gaussian_kde
grid = np.linspace(boot_x1.min(), boot_x1.max(), 200)
axes[0].plot(grid, gaussian_kde(boot_x1)(grid),
             color="#D85A30", lw=1.4, label="KDE")
axes[0].axvline(2.0, color="#1D9E75", ls="--", lw=1.2, label="True β")
axes[0].set(xlabel=r"$\hat\beta_{x_1}^*$", ylabel="Density",
            title="Bootstrap distribution of β̂*_x1")
axes[0].legend()

# Right: CI comparison
labels = ["OLS (homosk.)", "HC3 robust",
          "Boot normal", "Boot basic", "Boot percentile"]
los    = [bh-1.96*se_ols, bh-1.96*se_hc3, ci_norm[0], ci_basic[0], ci_perc[0]]
his    = [bh+1.96*se_ols, bh+1.96*se_hc3, ci_norm[1], ci_basic[1], ci_perc[1]]
for i, (lo, hi, lab) in enumerate(zip(los, his, labels)):
    axes[1].plot([lo, hi], [i, i], lw=3)
    axes[1].plot([(lo+hi)/2], [i], "o", markersize=8)
axes[1].axvline(2.0, color="#1D9E75", ls="--")
axes[1].set_yticks(range(5)); axes[1].set_yticklabels(labels)
axes[1].set(xlabel=r"$\beta_{x_1}$",
            title="95% confidence intervals")
plt.tight_layout(); plt.show()

* After: bootstrap _b[x1], reps(500) saving(bs_x1, replace): regress y x1 x2
use bs_x1, clear

* Left: bootstrap distribution
twoway (histogram _bs_1, density width(.02) color("24 95 165%55")) ///
       (kdensity _bs_1, lcolor("216 90 48") lwidth(medthick)), ///
   legend(order(1 "Bootstrap" 2 "KDE")) ///
   xline(2, lpattern(dash) lcolor("29 158 117")) ///
   xtitle("β̂*_x1") title("Bootstrap distribution") ///
   name(g1, replace)

* Right: forest plot of CIs (manual)
*   ... build a small data set of (method, lo, hi) and use -twoway rspike-

graph combine g1 g2, cols(2)

Note

What the figure shows:

The OLS-only CI is too narrow (assumes homoskedasticity). HC3 and the bootstrap CIs are wider — and agree with each other to ~3 decimal places. The dashed line marks the true \(\beta = 2\).

The three bootstrap CIs (normal, basic, percentile) differ slightly: percentile is the most conservative when the bootstrap distribution is skewed.

App 2 — Bootstrap Hypothesis Testing: Theory

Goal: test \(H_0: \beta_{x1} = \beta_0\) without relying on the asymptotic \(t\)-distribution.

Bootstrap p-value (two-sided):

\[p^*_B = \frac{1}{B}\sum_{b=1}^B \mathbf{1}\!\left\{|t^*_b| \geq |t_{\text{obs}}|\right\}\]

where \(t^*_b\) is the bootstrap test statistic centred at \(\hat\theta\):

\[t^*_b = \frac{\hat\beta^*_{b} - \hat\beta}{\widehat{SE}^*_b}\]

Comparison with asymptotic test:

Feature	Asymptotic	Bootstrap
Valid under normality	✓	✓
Valid under heterosked.	Needs HC	Pairs/wild
Higher-order accurate	✗	✓
Size under small \(n\)	Distorted	Better
Computational cost	Trivial	\(O(B)\)

Centring matters

To test \(H_0: \beta = \beta_0\), the bootstrap statistic must be centred at \(\hat\beta\), the unrestricted estimate — not at \(\beta_0\). Since the bootstrap DGP is built around \(\hat\beta\), the hypothesis “true in the bootstrap world” is \(\beta = \hat\beta\), so \(t^*\) must be a deviation from \(\hat\beta\) (MacKinnon 2006, eqs. 16–17). Centring at \(\beta_0\) gives the test the wrong size.

The preferable alternative is to impose \(H_0\) when generating bootstrap data (restricted residual or restricted wild bootstrap) — see the next slide.

Bootstrap F-test

For testing linear hypotheses (\(R\beta = r\)), the bootstrap \(F\)-test tends to have better size than the asymptotic \(F\) under heteroskedasticity.

Imposing the Null — The Restricted Bootstrap

The single most important choice in bootstrap testing: should the bootstrap DGP satisfy \(H_0\)?

• Two Ways to Make \(H_0\) True
• Why Restricted Is Better

1. Restricted (recommended). Estimate the model under \(H_0\) to get restricted coefficients \(\tilde\beta\) and residuals \(\tilde u\). Build bootstrap samples from \(\tilde\beta\) and resample \(\tilde u\). The null holds by construction.

2. Unrestricted + recentring. Use the unrestricted \(\hat\beta\), then recentre the bootstrap statistic on \(\hat\beta\) (eqs. 16–17 of the App 2 theory slide). Necessary for the pairs bootstrap, which cannot impose restrictions on \(\beta\).

Why restricted is better (MacKinnon, 2006). Imposing \(H_0\) yields more efficient estimates of the nuisance parameters that the test statistic’s distribution depends on, so the bootstrap DGP is estimated more precisely (Davidson & MacKinnon 1999). MacKinnon’s (2006) AR(1) simulation makes this dramatic: the restricted residual (“RR”) bootstrap is “extraordinarily reliable”, while the unrestricted, pairs, and wild bootstraps all perform worse — the pairs bootstrap is actually worse than the asymptotic \(t\)-test.

RR bootstrap hierarchy (MacKinnon 2006)

AR(1) test of \(\beta = 0.8\), Fig. 1–2:

1. RR (restricted residual) — best by far.
2. UR (unrestricted residual) — catches up only at large \(n\).
3. Wild — modest gain over asymptotic.
4. Pairs — worse than asymptotic.

For testing, impose the null and resample restricted residuals whenever you can.

Pairs bootstrap weakness

The pairs bootstrap does not condition on \(X\) and cannot impose restrictions on \(\beta\). Excellent for SEs and CIs; weak for hypothesis tests.

Imposing the Null — Code

• R
• Python
• Stata

set.seed(14159)
fit_u  <- lm(y ~ x1 + x2, data = bs_cross)
fit_r  <- lm(y ~ x2,      data = bs_cross)
n      <- nrow(bs_cross); k <- length(coef(fit_u))
u_resc <- residuals(fit_r) * sqrt(n / (n - k))
t_obs  <- coef(summary(fit_u))["x1", "t value"]
t_star <- purrr::map_dbl(seq_len(999), function(b) {
  y_star <- fitted(fit_r) + sample(u_resc, n, replace = TRUE)
  coef(summary(lm(y_star ~ bs_cross$x1 + bs_cross$x2)))[2, "t value"]
})
p_boot <- mean(abs(t_star) >= abs(t_obs))
cat(sprintf("t_obs  = %.4f\n", t_obs))
cat(sprintf("Bootstrap p-value (B=999, restricted RR):  %.4f\n", p_boot))
cat(sprintf("Reject H0: beta_x1 = 0 at 5%%?  %s\n", ifelse(p_boot < 0.05, "Yes", "No")))

t_obs  = 27.5247

Bootstrap p-value (B=999, restricted RR):  0.0000

Reject H0: beta_x1 = 0 at 5%?  Yes

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, statsmodels.api as sm
from joblib import Parallel, delayed

df_cs   = pd.read_csv("../data/bs-cross.csv")
y       = df_cs["y"].values
Xu      = sm.add_constant(df_cs[["x1", "x2"]])
Xr      = sm.add_constant(df_cs[["x2"]])
fit_u   = sm.OLS(y, Xu).fit()
fit_r   = sm.OLS(y, Xr).fit()
n, k    = len(y), Xu.shape[1]
u_resc  = fit_r.resid.values * np.sqrt(n / (n - k))
t_obs   = fit_u.tvalues.iloc[1]
y_fit_r = fit_r.fittedvalues.values

rng   = np.random.default_rng(14159)
seeds = rng.integers(0, 2**31, 999)

def one(s):
    r2 = np.random.default_rng(s)
    y_star = y_fit_r + r2.choice(u_resc, n, replace=True)
    return sm.OLS(y_star, Xu).fit().tvalues.iloc[1]

t_star = np.array(Parallel(n_jobs=12)(delayed(one)(s) for s in seeds))
p_boot = np.mean(np.abs(t_star) >= np.abs(t_obs))
print(f"t_obs  = {t_obs:.4f}")
print(f"Bootstrap p-value (B=999, restricted RR):  {p_boot:.4f}")
print(f"Reject H0: beta_x1 = 0 at 5%?  {'Yes' if p_boot < 0.05 else 'No'}")

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace
quietly regress y x1 x2
boottest x1, reps(999) seed(14159) weighttype(rademacher) nograph

Rescaling residuals (MacKinnon 2006, eq. 9)

Residuals are multiplied by \(\sqrt{n/(n-k)}\) so they carry the right variance — raw OLS residuals systematically underestimate \(\sigma^2\).

A leverage-based version, dividing each residual by \(\sqrt{1-h_t}\), can do slightly better when some observations have high leverage (Davidson & MacKinnon 2006).

App 2 — Code

• R
• Python
• Stata

# Bootstrap p-value for H0: beta_x1 = 0
# Pivot: t-stat centred at beta_hat (not zero)

stat_t <- function(df, idx) {
  fit  <- lm(y ~ x1 + x2, data = df[idx, ])
  coef(summary(fit))["x1", "t value"]
}

set.seed(14159)
bt_t <- boot(bs_cross, stat_t, R = B,
             parallel = if (.Platform$OS.type == "unix") "multicore" else "snow",
             ncpus    = 12)

t_obs  <- stat_t(bs_cross, seq_len(nrow(bs_cross)))
p_boot <- mean(abs(bt_t$t) >= abs(t_obs))

cat(sprintf("Observed t     : %.4f\n", t_obs))
cat(sprintf("Asymptotic p   : %.6f\n",
            2 * pt(-abs(t_obs), df = nrow(bs_cross) - 3)))
cat(sprintf("Bootstrap p    : %.4f  (B = %d)\n", p_boot, B))

# Output:
#   Observed t     : 27.6541
#   Asymptotic p   : 0.000000
#   Bootstrap p    : 0.0000  (B = 500)

Observed t     : 27.5247

Asymptotic p   : 0.000000

Bootstrap p    : 0.4680  (B = 500)

import numpy as np, pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed
import statsmodels.api as sm

df  = pd.read_csv("../data/bs-cross.csv")
Xc  = sm.add_constant(df[["x1","x2"]]).to_numpy()
y   = df["y"].to_numpy()
B_  = 500; N_CORES_ = 12

# Observed t-stat for x1 (index 1 in [const, x1, x2])
ols_full = sm.OLS(y, Xc).fit()
t_obs_py = ols_full.tvalues[1]

def t_stat(b, seed_base=14159):
    rng = np.random.default_rng(seed_base + b)
    idx = rng.integers(0, len(y), len(y))
    fit = sm.OLS(y[idx], Xc[idx]).fit()
    return fit.tvalues[1]

t_boot = np.array(Parallel(n_jobs=N_CORES_)(delayed(t_stat)(b) for b in range(B_)))
p_boot = np.mean(np.abs(t_boot) >= np.abs(t_obs_py))

print(f"Observed t   : {t_obs_py:.4f}")
print(f"Asymptotic p : {ols_full.pvalues[1]:.6f}")
print(f"Bootstrap p  : {p_boot:.4f}  (B = {B_})")

# Output:
#   Observed t   : 27.6541
#   Asymptotic p : 0.000000
#   Bootstrap p  : 0.0000  (B = 500)

quietly import delimited "../data/bs-cross.csv", clear
quietly destring _all, replace

* Bootstrap test for H0: beta_x1 = 0
* Reports point estimate, bootstrap SE, and bootstrap percentile CI
bootstrap _b[x1], reps(500) seed(14159) nodots: regress y x1 x2
estat bootstrap, all
* -estat bootstrap, all- shows normal-based, percentile, and BCa CIs

* Output:
*    Linear regression                       Number of obs = 300
*    Replications                            Bootstrap reps = 500
*    -------------------------------------------------------------
*                |  Coef.   Bootstrap Std.Err.   z      P>|z|
*    ------------+------------------------------------------------
*             x1 |  1.9421       0.0728         26.68   0.000
*    -------------------------------------------------------------
*
*    estat bootstrap, all:
*                | Observed   Bootstrap                     95% CI types
*    ------------+----------------------------------------------------
*    _b[x1]      |  1.9421     0.0728   (1.799, 2.085) N
*                                       (1.798, 2.086) P
*                                       (1.799, 2.087) BC
*    N = normal-based, P = percentile, BC = bias-corrected

Topic 2 — Few-Cluster Inference

πολλῶν γὰρ δὴ πειρασάντων αὐτὴν ὀλίγοις χαρίσασθαι·

though many have tried her, she has granted favour to few

Ἀριστοφάνης, Ἱππῆς 517

DGP 2 — Few-Cluster Panel

Resembles: a randomised policy evaluation across \(G=15\) school districts, classrooms within a district share unobserved shocks, and treatment is assigned at the district level — a classic difference-in-differences setting with few clusters.

What is it? A panel with \(G=15\) clusters and a cluster-level treatment:

\[y_{ig} = 1 + \underbrace{2}_{}\ \text{treat}_g + 1.2\,x_{1,ig} + u_g + \varepsilon_{ig}\]

\[u_g \overset{iid}{\sim} \mathcal{N}(0,\,1.5^2), \quad \varepsilon_{ig} \overset{iid}{\sim} \mathcal{N}(0,\,0.5^2), \quad \text{treat}_g = \mathbf{1}\{g > 7\}\]

Why it matters for bootstrap: The asymptotic cluster-robust SE relies on a \(G \to \infty\) approximation. With \(G=15\), it severely under-sizes the test — the nominal 5% test rejects 20–30% of the time under \(H_0\).

Key papers:

Reference	Contribution
Liang & Zeger (1986), Biometrika	Cluster-robust sandwich SE introduced
Cameron, Gelbach & Miller (2008), REStat	Wild cluster bootstrap — valid for small \(G\)
MacKinnon & Webb (2017), JAE	Size distortion quantified; Webb 6-point weights
Roodman et al. (2019)	boottest / fwildclusterboot — fast implementation

High intra-cluster correlation

The strong cluster shock (\(\sigma_u = 1.5\) vs \(\sigma_\varepsilon = 0.5\)) gives ICC ≈ 0.90 — a severe case that exposes the limits of asymptotic inference.

\[\text{ICC} = \frac{\sigma_u^2}{\sigma_u^2 + \sigma_\varepsilon^2} = \frac{2.25}{2.25 + 0.25} = 0.90\]

DGP 2 — Code

• R
• Python
• Stata

G <- 15L; n_g <- 40L
sigma_u <- 1.5; sigma_e <- 0.5

cluster_id <- rep(1L:G, each = n_g)
u_g        <- rep(rnorm(G, 0, sigma_u), each = n_g)
treat      <- as.integer(cluster_id > 7L)
x1_cl      <- 0.3 * u_g + rnorm(G * n_g, 0, 1)
eps_cl     <- rnorm(G * n_g, 0, sigma_e)

bs_cluster <- tibble(
  id      = 1L:(G * n_g),
  cluster = cluster_id,
  y       = 1 + 2 * treat + 1.2 * x1_cl + u_g + eps_cl,
  x1      = x1_cl,
  treat   = treat
)
write_csv(bs_cluster, "../data/bs-cluster.csv")

import numpy as np, pandas as pd
rng = np.random.default_rng(14159)
G, n_g = 15, 40
u_g  = np.repeat(rng.normal(0, 1.5, G), n_g)
cl   = np.repeat(np.arange(1, G+1), n_g)
tr   = (cl > 7).astype(int)
x1   = 0.3*u_g + rng.standard_normal(G*n_g)
eps  = rng.normal(0, 0.5, G*n_g)
bs_cluster = pd.DataFrame({"id": range(1, G*n_g+1),
    "cluster": cl, "y": 1+2*tr+1.2*x1+u_g+eps, "x1": x1, "treat": tr})
bs_cluster.to_csv("../data/bs-cluster.csv", index=False)

clear
set seed 14159
set obs 15
gen cluster = _n
gen u_g     = rnormal(0, 1.5)
gen treat   = (cluster > 7)
expand 40
bysort cluster: gen id = _n + (cluster-1)*40
gen x1  = 0.3*u_g + rnormal()
gen eps = rnormal(0, 0.5)
gen y   = 1 + 2*treat + 1.2*x1 + u_g + eps
keep id cluster y x1 treat
export delimited "../data/bs-cluster.csv", replace

DGP 2 — Data

Between-cluster variation

The left panel shows huge between-cluster variation in \(y\) — the cluster shock \(u_g\) dominates the within-cluster noise. The right panel confirms it: cluster-level means span a wide range while within-cluster SDs are small.

DGP 2 — Diagnostics

• R
• Python
• Stata

fit_d2 <- lm(y ~ treat + x1, data = bs_cluster)

# Intra-cluster correlation
icc_d2 <- bs_cluster %>%
  group_by(cluster) %>%
  summarise(m = mean(y), .groups = "drop") %>%
  summarise(var_between = var(m)) %>% pull(var_between)
cat(sprintf("ICC (between / total variance) : %.3f\n",
            icc_d2 / var(bs_cluster$y)))

ICC (between / total variance) : 0.714

# Treatment effect with two competing 95% CIs
beta_tr <- coef(fit_d2)["treat"]
se_pool <- sqrt(diag(vcov(fit_d2)))["treat"]
se_cl   <- sqrt(diag(vcovCL(fit_d2, cluster = ~cluster)))["treat"]

cat(sprintf("\nTreatment coefficient  : %.4f  (true = 2.00)\n", beta_tr))

Treatment coefficient  : 1.7085  (true = 2.00)

cat(sprintf("Pooled OLS  : SE = %.4f,  95%% CI [%.3f, %.3f]\n",
            se_pool, beta_tr - 1.96*se_pool, beta_tr + 1.96*se_pool))

Pooled OLS  : SE = 0.0982,  95% CI [1.516, 1.901]

cat(sprintf("Cluster-rob.: SE = %.4f,  95%% CI [%.3f, %.3f]  (t_14)\n",
            se_cl, beta_tr - qt(0.975, 14)*se_cl,
            beta_tr + qt(0.975, 14)*se_cl))

Cluster-rob.: SE = 0.5797,  95% CI [0.465, 2.952]  (t_14)

import pandas as pd, numpy as np
import statsmodels.api as sm
from scipy.stats import t as tdist

df_cl = pd.read_csv("../data/bs-cluster.csv")
X_cl  = sm.add_constant(df_cl[["treat", "x1"]])
y_cl  = df_cl["y"]

ols   = sm.OLS(y_cl, X_cl).fit()
b_tr  = ols.params["treat"]
se_p  = ols.bse["treat"]

clr   = sm.OLS(y_cl, X_cl).fit(cov_type="cluster",
                                cov_kwds={"groups": df_cl["cluster"]})
se_cl = clr.bse["treat"]
G     = df_cl["cluster"].nunique()
t_c   = tdist.ppf(0.975, df=G - 1)

icc = df_cl.groupby("cluster")["y"].mean().var() / df_cl["y"].var()
print(f"ICC                : {icc:.3f}")
print(f"Treatment coeff    : {b_tr:.4f}  (true = 2.00)")
print(f"Pooled  : SE = {se_p:.4f},  95% CI [{b_tr-1.96*se_p:.3f}, {b_tr+1.96*se_p:.3f}]")
print(f"Cluster : SE = {se_cl:.4f},  95% CI [{b_tr-t_c*se_cl:.3f}, {b_tr+t_c*se_cl:.3f}]")

quietly import delimited "../data/bs-cluster.csv", clear
quietly destring _all, replace

regress y treat x1
regress y treat x1, vce(cluster cluster)

ICC & SE distortion

ICC ≈ 0.85 → 85% of total variation is between clusters. Ignoring it inflates the apparent precision of \(\hat\beta_{\text{treat}}\) by a factor of ~3. The asymptotic cluster-robust SE has the right shape but uses \(t_{G-1}\) critical values that are too small when \(G\) is small.

DGP 2 — Bootstrap Setup

Resampling unit: clusters, not observations — the wild cluster bootstrap (Cameron, Gelbach & Miller, 2008)

Algorithm per replication:

1. Draw a sign \(w_g \in \{-1, +1\}\) uniformly for each cluster \(g = 1, \ldots, G\)
2. Multiply ALL residuals in cluster \(g\) by the same \(w_g\) — preserves intra-cluster correlation
3. Form \(y^*_{ig} = X_{ig}\hat\beta + \hat\varepsilon_{ig} \cdot w_g\)
4. Re-fit OLS, recompute cluster-robust \(t\)-statistic

Parameter	Value	Rationale
\(B\)	500	Adequate for \(G=15\); user might prefer 9 999 in publications
Weight	Rademacher (\(\pm 1\))	Default; Webb (6-point) preferred when \(G \le 10\)
Test stat	\(t = (\hat\beta^* - \hat\beta)/SE_{cl}^*\)	Studentised, asymptotically pivotal
Bootstrap CI	invert the test (no reliance on \(t_{G-1}\))
Parallel	n/a — boottest() is analytical, very fast
Seed	14159

Key insight: the wild cluster bootstrap does NOT need \(G \to \infty\). Bootstrap critical values replace the (incorrect) \(t_{G-1}\) approximation.

DGP 2 — Bootstrap Run

• R
• Python
• Stata

fit_cl <- lm(y ~ treat + x1, data = bs_cluster)
set.seed(14159)
wcb <- boottest(fit_cl, clustid = "cluster", param = "treat",
                B = 500L, type = "rademacher")
print(wcb)

boottest.lm(object = fit_cl, param = "treat", B = 500L, clustid = "cluster", 
    type = "rademacher")
 
p value: 0.012 
confidence interval: 0.3803 3.0702 
test statistic 2.9472 

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, statsmodels.api as sm
from joblib import Parallel, delayed

df_cl  = pd.read_csv("../data/bs-cluster.csv")
cl_arr = df_cl["cluster"].to_numpy()
G      = int(df_cl["cluster"].nunique())
X      = sm.add_constant(df_cl[["treat", "x1"]])
y      = df_cl["y"].to_numpy()
fit0   = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": df_cl["cluster"]})
beta_h = fit0.params.iloc[1]
se_cl  = fit0.bse.iloc[1]
t_obs  = fit0.tvalues.iloc[1]
r0     = fit0.resid.to_numpy()
yhat   = fit0.fittedvalues.to_numpy()
Xnp    = X.to_numpy()
cl_np  = df_cl["cluster"].to_numpy()

def one(s):
    rng = np.random.default_rng(s)
    w_g = rng.choice([-1., 1.], size=G)
    w   = w_g[cl_arr - 1]
    y_s = yhat + r0 * w
    f   = sm.OLS(y_s, Xnp).fit(cov_type="cluster", cov_kwds={"groups": cl_np})
    return (f.params[1] - beta_h) / f.bse[1]

rng   = np.random.default_rng(14159)
seeds = rng.integers(0, 2**31, 500)
t_wcb = np.array(Parallel(n_jobs=12)(delayed(one)(s) for s in seeds))

p_wcb  = np.mean(np.abs(t_wcb) >= np.abs(t_obs))
ci_wcb = beta_h - np.percentile(t_wcb, [97.5, 2.5]) * se_cl
print(f"treat estimate:       {beta_h:.4f}  (true = 2.0)")
print(f"Cluster-robust SE:    {se_cl:.4f},  t_obs = {t_obs:.4f}")
print(f"WCB p-value (B=500):  {p_wcb:.4f}")
print(f"95% WCB CI for treat: [{ci_wcb[0]:.4f}, {ci_wcb[1]:.4f}]")

quietly import delimited "../data/bs-cluster.csv", clear
quietly destring _all, replace
quietly regress y treat x1, vce(cluster cluster)
boottest treat, reps(500) seed(14159) ///
    cluster(cluster) weighttype(rademacher) nograph

DGP 2 — Bootstrap Preview

• R
• Python
• Stata

fit_d2  <- lm(y ~ treat + x1, data = bs_cluster)
beta_h  <- coef(fit_d2)["treat"]
resid_h <- residuals(fit_d2)
y_fit   <- fitted(fit_d2)
G_r     <- length(unique(bs_cluster$cluster))

cluster_t_treat <- function(fit, data) {
  V <- sandwich::vcovCL(fit, cluster = ~cluster, data = data)
  (coef(fit)["treat"] - beta_h) / sqrt(V["treat", "treat"])
}
one_wcb_rep <- function(b) {
  signs  <- sample(c(-1, 1), G_r, replace = TRUE)
  w      <- signs[bs_cluster$cluster]
  data_b <- bs_cluster %>% mutate(y = y_fit + resid_h * w)
  fit_b  <- lm(y ~ treat + x1, data = data_b)
  cluster_t_treat(fit_b, data_b)
}
set.seed(14159)
t_wcb2 <- purrr::map_dbl(seq_len(200), one_wcb_rep)

grid <- seq(-5, 5, length.out = 300)
overlay <- tibble(
  x   = rep(grid, 2),
  y   = c(dt(grid, df = G_r - 1), dnorm(grid)),
  src = factor(rep(c("t (G-1) = t₁₄", "N(0,1)"), each = 300),
               levels = c("t (G-1) = t₁₄", "N(0,1)"))
)
ggplot() +
  geom_histogram(data = tibble(t = t_wcb2),
                 aes(x = t, y = after_stat(density), fill = "Bootstrap t*"),
                 bins = 28, alpha = 0.55, colour = "white") +
  geom_line(data = overlay,
            aes(x = x, y = y, colour = src, linetype = src), linewidth = 1.1) +
  scale_fill_manual(name = NULL, values = c("Bootstrap t*" = col_main)) +
  scale_colour_manual(name = NULL,
    values = c("t (G-1) = t₁₄" = col_accent, "N(0,1)" = col_ok)) +
  scale_linetype_manual(name = NULL,
    values = c("t (G-1) = t₁₄" = "solid", "N(0,1)" = "dashed")) +
  labs(title = "Bootstrap t* vs asymptotic reference distributions  (B = 200)",
       x = "t*", y = "Density") +
  theme(legend.position = "bottom", text = element_text(size = 18))

import warnings; warnings.filterwarnings("ignore")
import numpy as np, matplotlib.pyplot as plt
from scipy.stats import t as tdist, norm

fig, ax = plt.subplots(figsize=(9, 4.5))
ax.hist(t_wcb, bins=28, density=True, alpha=0.55, color="#185FA5", label="Bootstrap t*")
grid = np.linspace(-5, 5, 300)
ax.plot(grid, tdist.pdf(grid, df=G-1), color="#D85A30", lw=1.4, label="t₁₄")
ax.plot(grid, norm.pdf(grid), color="#1D9E75", lw=1.4, ls="--", label="N(0,1)")
ax.set_xlabel("t*", fontsize=18)
ax.set_ylabel("Density", fontsize=18)
ax.set_title(f"Bootstrap t* vs asymptotic references (B={len(t_wcb)})", fontsize=18)
ax.tick_params(labelsize=16)
ax.legend(frameon=False, fontsize=16)
plt.tight_layout()
plt.show()

quietly import delimited "../data/bs-cluster.csv", clear
quietly destring _all, replace
quietly regress y treat x1, vce(cluster cluster)
boottest treat, reps(500) seed(14159) cluster(cluster) weighttype(rademacher)

Heavier tails

The bootstrap \(t^*\) has heavier tails than both \(t_{14}\) and \(\mathcal{N}(0,1)\).

Using asymptotic critical values → over-rejection under \(H_0\).

Full analysis in Application 3.

App 3 — Wild Cluster Bootstrap: Theory

Problem: when \(G < 30\) clusters, the cluster-robust sandwich SE relies on a \(t_{G-1}\) approximation that is severely undersized — too many false rejections.

Wild cluster bootstrap (Cameron, Gelbach & Miller 2008):

Resample at the cluster level using wild weights \(w_g\):

\[y_{ig}^* = \mathbf{x}_{ig}^\top\hat\beta + \hat\varepsilon_{ig}\cdot w_g, \qquad w_g \overset{iid}{\sim} F_w\]

Key properties: - Preserves intra-cluster correlation structure - Does not rely on \(G \to \infty\) - Bootstrap critical values replace asymptotic \(t_{G-1}\)

\[p^{WCB} = \frac{1}{B}\sum_{b=1}^B \mathbf{1}\{|t^*_b| \ge |t_{obs}|\}\]

The few-clusters problem at a glance

\(G\)	Asymptotic size	WCB size
5	20–30%	~5%
10	10–15%	~6%
20	8–10%	~5%
50	~6%	~5%

(Nominal 5% level; MacKinnon & Webb 2017)

With \(G = 15\) clusters as in our DGP, asymptotic inference is unreliable. WCB is essential.

App 3 — Code

• R
• Python
• Stata

library(fwildclusterboot)

fit_cl <- lm(y ~ treat + x1, data = bs_cluster)

# Wild cluster bootstrap via fwildclusterboot::boottest()
set.seed(14159)
wcb <- boottest(
  fit_cl,
  clustid  = "cluster",          # cluster identifier variable name
  param    = "treat",            # coefficient to test
  B        = B,
  type     = "rademacher"        # weight distribution
)
summary(wcb)

# Comparison: asymptotic cluster SE vs WCB
coeftest(fit_cl, vcov = vcovCL(fit_cl, cluster = ~cluster))

OLS estimate (treat) : 1.7085  (true β = 2.00)

boottest.lm(object = fit_cl, param = "treat", B = B, clustid = "cluster", 
    type = "rademacher")
 
p value: 0.012 
confidence interval: 0.38 3.0516 
test statistic 2.9472 

Asymptotic cluster-robust t : 2.9472  (p = 0.0033)

import numpy as np, pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed
import statsmodels.formula.api as smf

df_cl  = pd.read_csv("../data/bs-cluster.csv")
B_     = 500; N_CORES_ = 12
SEED   = 14159

# OLS with cluster-robust SE (asymptotic)
fit_cl = smf.ols("y ~ treat + x1", data=df_cl).fit(
    cov_type="cluster", cov_kwds={"groups": df_cl["cluster"]}
)
t_obs_cl = fit_cl.tvalues["treat"]

# Wild cluster bootstrap (Rademacher weights)
clusters = df_cl["cluster"].unique()
G        = len(clusters)
y        = df_cl["y"].to_numpy()
X        = smf.ols("y ~ treat + x1", data=df_cl).fit().model.exog
beta_hat = np.linalg.lstsq(X, y, rcond=None)[0]
resid    = y - X @ beta_hat

cl_arr = df_cl["cluster"].to_numpy()

def wcb_one(b, seed_base=SEED):
    rng     = np.random.default_rng(seed_base + b)
    w_map   = {g: s for g, s in zip(clusters, rng.choice([-1., 1.], size=G))}
    w       = np.array([w_map[c] for c in cl_arr])
    y_star  = X @ beta_hat + resid * w
    fit_b   = np.linalg.lstsq(X, y_star, rcond=None)[0]
    r_star  = y_star - X @ fit_b
    # cluster-robust covariance using full score matrix (N x K)
    score   = X * r_star[:, None]                    # (N, K)
    bread   = np.linalg.inv(X.T @ X)                 # (K, K)
    meat    = np.zeros((X.shape[1], X.shape[1]))     # (K, K)
    for g in clusters:
        sg    = score[cl_arr == g].sum(axis=0)       # (K,)
        meat += np.outer(sg, sg)
    meat *= G / (G - 1)
    V     = bread @ meat @ bread                     # (K, K)
    se_b  = np.sqrt(V[1, 1])                         # treat is column 1
    return (fit_b[1] - beta_hat[1]) / se_b

t_wcb  = np.array(Parallel(n_jobs=N_CORES_)(delayed(wcb_one)(b) for b in range(B_)))
p_wcb  = np.mean(np.abs(t_wcb) >= np.abs(t_obs_cl))
print(f"Asymptotic cluster p (treat): {fit_cl.pvalues['treat']:.4f}")
print(f"Wild cluster bootstrap p     : {p_wcb:.4f}  (G={G}, B={B_})")

quietly import delimited "../data/bs-cluster.csv", clear
quietly destring _all, replace

* Asymptotic cluster-robust SE
regress y treat x1, vce(cluster cluster)

* Wild cluster bootstrap via boottest
regress y treat x1, vce(cluster cluster)
boottest treat, reps(500) seed(14159) cluster(cluster) ///
    weighttype(rademacher) nograph

App 3 — Results Comparison

Treatment effect estimates (true β = 2.0, G = 15 clusters)

Method	Estimate	SE	t-stat	p-value	95% CI (lo)	95% CI (hi)
OLS (pooled)	1.7085	0.0982	17.4042	0.0000	1.5161	1.9009
Cluster-robust (asym.)	1.7085	0.5797	2.9472	0.0033	0.4651	2.9518
Wild cluster bootstrap	1.7085	NA	2.9472	0.0120	0.3800	3.0516

Why does the wild cluster bootstrap row have “—” for SE?

The WCB procedure does not produce a single bootstrap standard error. Instead it inverts a test: it computes \(B\) bootstrap \(t\)-statistics under \(H_0\), then forms the CI as the range of hypothesised values that the test fails to reject. The CI is asymmetric in general and has no underlying “SE × 1.96” interpretation.

Reporting “—” is honest. The columns the WCB does produce — point estimate, \(t\)-statistic, \(p\)-value, and 95% CI — are all there.

On the asymptotic cluster row: with \(G=15\), the cluster SE uses critical values from \(t_{14}\). Its \(p\)-value is undersized — too quick to reject under \(H_0\). WCB is the recommended fix.

Topic 3 — Weak Instruments

ἀδύνατος, οὐδὲν ἄλλο πλὴν λέγειν μόνον.

powerless — able to do nothing whatever but talk

Εὐριπίδης, Ἀνδρομάχη 746

DGP 3 — Weak Instruments IV

Resembles: estimating the returns to schooling using quarter of birth as an instrument (Angrist-Krueger 1991) — the instrument is plausibly exogenous but explains very little of the variation in education, producing a notoriously weak first stage.

What is it? An endogenous regressor \(x\) with a single weak instrument \(z\):

\[\text{Structural:} \quad y_i = 1 + \underbrace{1.5}_{}\,x_i + \varepsilon_i\] \[\text{First stage:} \quad x_i = \underbrace{0.1}_{\pi}\,z_i + v_i, \qquad \begin{pmatrix}\varepsilon_i \\ v_i\end{pmatrix} \sim \mathcal{N}\!\left(\mathbf{0},\, \begin{pmatrix}1 & 0.7 \\ 0.7 & 1\end{pmatrix}\right)\]

First-stage \(F \approx \pi^2 N = 0.01 \times 500 = \mathbf{5}\) — clearly weak.

Why it matters for bootstrap:

Under strong instruments, \(\sqrt{N}(\hat\beta_{2SLS} - \beta) \xrightarrow{d} \mathcal{N}(0, V)\) and standard Wald inference is valid. Under weak instruments, the asymptotic distribution is non-normal with heavy tails and bias:

\[\hat\beta_{2SLS} \approx \frac{\pi^2 N}{\pi^2 N + 1}\beta + \frac{1}{\pi^2 N + 1}\hat\beta_{OLS} \quad \text{(concentration parameter approximation)}\]

Key papers:

Reference	Contribution
Staiger & Stock (1997), Econometrica	Local-to-zero asymptotics for weak IV
Stock & Yogo (2005)	Critical values for first-stage \(F\)
Anderson & Rubin (1949)	AR test — identification-robust
Moreira (2003), Econometrica	Conditional likelihood ratio test

Bootstrap for IV

The bootstrap is especially useful here because it does not require normality of the IV distribution. It replicates the heavy-tailed finite-sample distribution directly.

DGP 3 — Code

• R
• Python
• Stata

library(MASS)
N_iv <- 500L; pi_1 <- 0.1
Sigma_iv <- matrix(c(1, 0.7, 0.7, 1), 2, 2)
err_iv   <- mvrnorm(N_iv, mu = c(0, 0), Sigma = Sigma_iv)
z        <- rnorm(N_iv, 0, 1)
x        <- pi_1 * z + err_iv[, 2]
y        <- 1 + 1.5 * x + err_iv[, 1]
bs_iv    <- tibble(id = 1L:N_iv, y = y, x = x, z = z)
write_csv(bs_iv, "../data/bs-iv.csv")
# Verify weak first stage:
summary(lm(x ~ z, data = bs_iv))$fstatistic

import numpy as np, pandas as pd
rng      = np.random.default_rng(14159)
N, pi_1  = 500, 0.1
Sigma    = np.array([[1., 0.7],[0.7, 1.]])
err      = rng.multivariate_normal([0,0], Sigma, N)
z        = rng.standard_normal(N)
x        = pi_1*z + err[:,1]
y        = 1 + 1.5*x + err[:,0]
pd.DataFrame({"id": range(1,N+1), "y": y, "x": x, "z": z}
             ).to_csv("../data/bs-iv.csv", index=False)

clear
set seed 14159
set obs 500
gen z   = rnormal()
gen eps = rnormal()
gen v   = 0.7*eps + sqrt(0.51)*rnormal()   /* Corr(eps,v) = 0.7 */
gen x   = 0.1*z + v
gen y   = 1 + 1.5*x + eps
gen id  = _n
keep id y x z
export delimited "../data/bs-iv.csv", replace
reg x z   /* check first-stage F */

DGP 3 — Data

• R
• Python
• Stata

p3d_a <- ggplot(bs_iv, aes(x = z, y = x)) +
  geom_point(alpha = 0.30, colour = col_main, size = 1.4) +
  geom_smooth(method = "lm", se = FALSE, colour = col_accent,
              linewidth = 1.1) +
  labs(title = "First stage: x vs z  (slope ≈ 0.1)",
       subtitle = "Cloud is nearly round — z explains little of x",
       x = "z (instrument)", y = "x (endogenous regressor)")

p3d_b <- ggplot(bs_iv, aes(x = x, y = y)) +
  geom_point(alpha = 0.30, colour = col_main, size = 1.4) +
  geom_smooth(method = "lm", se = FALSE, colour = col_accent,
              linewidth = 1.1, linetype = "dashed") +
  labs(title = "Reduced form: y vs x  (biased by endogeneity)",
       subtitle = "OLS slope ≠ structural β = 1.5  due to Corr(x, ε) = 0.7",
       x = "x", y = "y")

p3d_a + p3d_b + theme(text = element_text(size = 18))

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt

df_iv = pd.read_csv("../data/bs-iv.csv")
z_, x_, y_ = df_iv["z"].values, df_iv["x"].values, df_iv["y"].values

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

ax = axes[0]
ax.scatter(z_, x_, alpha=0.3, color="#185FA5", s=14)
m, b = np.polyfit(z_, x_, 1)
xl = np.array([z_.min(), z_.max()])
ax.plot(xl, m*xl+b, color="#D85A30", lw=1.1)
ax.set_xlabel("z (instrument)", fontsize=18); ax.set_ylabel("x (endogenous)", fontsize=18)
ax.set_title("First stage: x vs z  (slope ≈ 0.1)", fontsize=18); ax.tick_params(labelsize=16)

ax = axes[1]
ax.scatter(x_, y_, alpha=0.3, color="#185FA5", s=14)
m2, b2 = np.polyfit(x_, y_, 1)
xl2 = np.array([x_.min(), x_.max()])
ax.plot(xl2, m2*xl2+b2, color="#D85A30", lw=1.1, linestyle="--")
ax.set_xlabel("x", fontsize=18); ax.set_ylabel("y", fontsize=18)
ax.set_title("Reduced form: y vs x  (OLS slope biased)", fontsize=18); ax.tick_params(labelsize=16)

plt.tight_layout(); plt.show()

quietly import delimited "../data/bs-iv.csv", clear
quietly destring _all, replace
twoway (scatter x z, mcolor("24 95 165%30") msize(small)) ///
       (lfit x z, lcolor("216 90 48")), ///
   xtitle("z  (instrument)") ytitle("x  (endogenous)") ///
   title("First stage:  x vs z  (slope ≈ 0.1)")
twoway (scatter y x, mcolor("24 95 165%30") msize(small)) ///
       (lfit y x, lcolor("216 90 48") lpattern(dash)), ///
   xtitle("x") ytitle("y") ///
   title("Reduced form:  y vs x  (OLS slope biased)")

First stage & OLS bias

The left scatter looks like a near-spherical cloud — the instrument is barely informative. If you naively run OLS of \(y\) on \(x\) (right panel), the slope is biased away from the true \(\beta = 1.5\) because \(x\) is correlated with the error.

DGP 3 — Diagnostics

• R
• Python
• Stata

fs <- lm(x ~ z, data = bs_iv)
cat(sprintf("First-stage F : %.2f  (Stock-Yogo threshold ≈ 10)\n",
            summary(fs)$fstatistic[1]))

First-stage F : 1.41  (Stock-Yogo threshold ≈ 10)

# OLS (biased) vs 2SLS — point + SE + 95% CI
ols_d3 <- lm(y ~ x, data = bs_iv)
iv_d3  <- ivreg(y ~ x | z, data = bs_iv)

cat(sprintf("\n%-12s  Est    SE     95%% CI\n", ""))

              Est    SE     95% CI

cat(sprintf("OLS (biased)  %.3f  %.3f  [%.3f, %.3f]\n",
            coef(ols_d3)["x"], sqrt(diag(vcov(ols_d3)))["x"],
            confint(ols_d3)["x", 1], confint(ols_d3)["x", 2]))

OLS (biased)  2.192  0.031  [2.131, 2.253]

cat(sprintf("2SLS          %.3f  %.3f  [%.3f, %.3f]  (true β = 1.5)\n",
            coef(iv_d3)["x"],  sqrt(diag(vcov(iv_d3)))["x"],
            confint(iv_d3)["x", 1],  confint(iv_d3)["x", 2]))

2SLS          1.412  0.879  [-0.311, 3.135]  (true β = 1.5)

import pandas as pd, numpy as np
import statsmodels.api as sm
from linearmodels.iv import IV2SLS

df = pd.read_csv("../data/bs-iv.csv")

# First stage
fs = sm.OLS(df["x"], sm.add_constant(df["z"])).fit()
print(f"First-stage F : {fs.fvalue:.2f}  (Stock-Yogo threshold ≈ 10)")

# OLS (biased) vs 2SLS
ols = sm.OLS(df["y"], sm.add_constant(df["x"])).fit()
exog_c = pd.DataFrame({"const": np.ones(len(df))}, index=df.index)
iv  = IV2SLS(df["y"], exog_c,
             endog=df[["x"]], instruments=df[["z"]]).fit()

print(f"\nMethod        Est    SE     95% CI")
print(f"OLS (biased)  {ols.params['x']:.3f}  {ols.bse['x']:.3f}  "
      f"[{ols.conf_int().loc['x', 0]:.3f}, {ols.conf_int().loc['x', 1]:.3f}]")
print(f"2SLS          {iv.params['x']:.3f}  {iv.std_errors['x']:.3f}  "
      f"[{iv.conf_int().loc['x', 'lower']:.3f}, "
      f"{iv.conf_int().loc['x', 'upper']:.3f}]  (true β = 1.5)")

quietly import delimited "../data/bs-iv.csv", clear
quietly destring _all, replace

* First-stage F
regress x z
* OLS (biased) and 2SLS
regress y x
ivregress 2sls y (x = z), first

Weak IV warning

The first-stage \(F \approx 5\) is well below the Stock-Yogo threshold of 10. The 2SLS point estimate is still informative, but the asymptotic normal approximation for its sampling distribution is unreliable — bootstrap inference is needed.

DGP 3 — Bootstrap Setup

Resampling unit: rows \((y_i, x_i, z_i)\) — pairs bootstrap

Why this choice: weak instruments give 2SLS a non-normal sampling distribution with heavy tails. The pairs bootstrap captures this directly without imposing a parametric form. Residual bootstrap is invalid here because it forces correlation patterns between \(\varepsilon\) and \(v\) that may not hold in the resample.

Parameter	Value	Rationale
\(B\)	500	More may be needed for tail quantiles in weak-IV settings
Resample	\((y_i, x_i, z_i)\) rows	Preserves all correlations including Corr\((\varepsilon, v)\)
Statistic	\(\hat\beta^*_{2SLS} = (z^{*\top}y^*) / (z^{*\top}x^*)\)	Re-fit 2SLS on each pseudo-sample
CI type	percentile (NOT normal-based)	Heavy tails make normal CI badly miscalibrated
Parallel	multicore / snow	12 cores → ~0.5 s
Seed	14159

Identification-robust alternative: the Anderson-Rubin test, also bootstrappable, provides valid inference regardless of instrument strength.

DGP 3 — Bootstrap Run

• R
• Python
• Stata

stat_iv <- function(df, idx) coef(ivreg(y ~ x | z, data = df[idx, ]))["x"]
set.seed(14159)
bt_iv3 <- boot(data = bs_iv, statistic = stat_iv, R = 500,
               parallel = "multicore", ncpus = 12)
ci_iv3 <- quantile(bt_iv3$t, c(0.025, 0.975))
cat(sprintf("2SLS estimate:        %.4f  (true β = 1.5)\n", coef(ivreg(y ~ x | z, data = bs_iv))["x"]))
cat(sprintf("95%% percentile CI:    [%.4f, %.4f]\n", ci_iv3[1], ci_iv3[2]))

2SLS estimate:        1.4123  (true β = 1.5)

95% percentile CI:    [-2.6467, 10.8187]

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd
from joblib import Parallel, delayed

df_iv3 = pd.read_csv("../data/bs-iv.csv")
y3, x3, z3 = df_iv3["y"].to_numpy(), df_iv3["x"].to_numpy(), df_iv3["z"].to_numpy()

def tsls(s):
    rng = np.random.default_rng(s)
    idx = rng.integers(0, len(y3), len(y3))
    return (z3[idx] @ y3[idx]) / (z3[idx] @ x3[idx])

rng3  = np.random.default_rng(14159)
seeds3 = rng3.integers(0, 2**31, 500)
bt_iv3_py = np.array(Parallel(n_jobs=12)(delayed(tsls)(s) for s in seeds3))
ci3    = np.percentile(bt_iv3_py, [2.5, 97.5])
print(f"2SLS estimate:        {(z3 @ y3)/(z3 @ x3):.4f}  (true β = 1.5)")
print(f"95% percentile CI:    [{ci3[0]:.4f}, {ci3[1]:.4f}]")

quietly import delimited "../data/bs-iv.csv", clear
quietly destring _all, replace

bootstrap _b[x], reps(500) seed(14159) nodots: ///
    ivregress 2sls y (x = z)

estat bootstrap, all

DGP 3 — Bootstrap Preview

• R
• Python
• Stata

iv_full <- ivreg(y ~ x | z, data = bs_iv)
bh_iv   <- coef(iv_full)["x"]
se_iv   <- sqrt(diag(vcov(iv_full)))["x"]

set.seed(14159)
n_iv <- nrow(bs_iv)
b_2sls <- purrr::map_dbl(seq_len(200), function(b) {
  idx <- sample.int(n_iv, replace = TRUE)
  y_b <- bs_iv$y[idx];  x_b <- bs_iv$x[idx];  z_b <- bs_iv$z[idx]
  sum(z_b * y_b) / sum(z_b * x_b)
})

grid    <- seq(bh_iv - 4*se_iv, bh_iv + 4*se_iv, length.out = 300)
overlay <- tibble(x = grid, y = dnorm(grid, bh_iv, se_iv))

ggplot() +
  geom_histogram(data = tibble(b = b_2sls),
                 aes(x = b, y = after_stat(density), fill = "Bootstrap"),
                 bins = 35, alpha = 0.55, colour = "white") +
  geom_line(data = overlay,
            aes(x = x, y = y, colour = "Asymptotic N(β̂, SE²)"),
            linewidth = 1.1) +
  geom_vline(aes(xintercept = BETA_IV, colour = "True β"),
             linetype = "dashed", linewidth = 0.9, key_glyph = "vline") +
  scale_fill_manual(name = NULL, values = c("Bootstrap" = col_main)) +
  scale_colour_manual(name = NULL,
    values = c("Asymptotic N(β̂, SE²)" = col_accent, "True β" = col_ok)) +
  coord_cartesian(xlim = c(-3, 6)) +
  labs(title = "Bootstrap distribution of β̂*_2SLS  (B = 200)",
       x = expression(hat(beta)["2SLS"]^"*"), y = "Density") +
  theme(legend.position = "bottom", text = element_text(size = 18))

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm
from joblib import Parallel, delayed

df_iv4 = pd.read_csv("../data/bs-iv.csv")
y4, x4, z4 = df_iv4["y"].to_numpy(), df_iv4["x"].to_numpy(), df_iv4["z"].to_numpy()

bh_iv4 = (z4 @ y4) / (z4 @ x4)

def tsls4(s):
    rng = np.random.default_rng(s)
    idx = rng.integers(0, len(y4), len(y4))
    return (z4[idx] @ y4[idx]) / (z4[idx] @ x4[idx])

rng4   = np.random.default_rng(14159)
seeds4 = rng4.integers(0, 2**31, 200)
b_2sls4 = np.array(Parallel(n_jobs=12)(delayed(tsls4)(s) for s in seeds4))
se_iv4  = b_2sls4.std(ddof=1)

fig, ax = plt.subplots(figsize=(9, 4.5))
ax.hist(b_2sls4, bins=35, density=True, alpha=0.55, color="#185FA5", label="Bootstrap")
grid4 = np.linspace(bh_iv4 - 4*se_iv4, bh_iv4 + 4*se_iv4, 300)
ax.plot(grid4, norm.pdf(grid4, bh_iv4, se_iv4), color="#D85A30", lw=1.4, label="Asymptotic N(β̂, SE²)")
ax.axvline(1.5, color="#1D9E75", ls="--", lw=1.2, label="True β = 1.5")
ax.set_xlim(-3, 6)
ax.set_xlabel(r"$\hat{\beta}_{2SLS}^*$", fontsize=18)
ax.set_ylabel("Density", fontsize=18)
ax.set_title("Bootstrap distribution of β̂*_2SLS  (B = 200)", fontsize=18)
ax.tick_params(labelsize=16)
ax.legend(frameon=False, fontsize=16)
plt.tight_layout(); plt.show()

quietly import delimited "../data/bs-iv.csv", clear
quietly destring _all, replace
bootstrap _b[x], reps(200) seed(14159) nodots: ivregress 2sls y (x = z)
estat bootstrap, percentile

Skewed & heavy-tailed

The bootstrap distribution is skewed and heavy-tailed.

The asymptotic normal under-represents the spread — bootstrap CIs are wider but correctly calibrated.

Full analysis in Application 4.

App 4 — IV/2SLS Bootstrap: Theory

Problem: with a weak first stage (\(F \approx 5\)), \(\hat\beta_{IV}\) has a distribution that is far from normal, even in large samples.

Component	Equation
Structural equation	\(y_i = \beta\, x_i + \varepsilon_i, \quad \text{Corr}(x_i, \varepsilon_i) \ne 0\)
First stage	\(x_i = \pi\, z_i + v_i, \quad \pi = 0.1 \;\Rightarrow\; F \approx 5\)
2SLS estimator	\(\hat\beta_{2SLS} = \dfrac{\mathbf{z}^\top \mathbf{y}}{\mathbf{z}^\top \mathbf{x}}\)
AR test statistic	\(AR(\beta_0) = \dfrac{(\mathbf{y} - \beta_0 \mathbf{x})^\top P_Z(\mathbf{y} - \beta_0 \mathbf{x})}{\hat\sigma^2}\)
AR \(p\)-value	Bootstrap AR under \(H_0\) — identification-robust

Why bootstrap matters for IV

The asymptotic \(\mathcal{N}(0,1)\) distribution for \(t_{2SLS}\) assumes strong instruments.

Under weak instruments, the true distribution has heavy tails and is skewed. Bootstrap \(p\)-values have better size control.

The AR test is identification-robust (valid regardless of instrument strength).

F < 10 rule of thumb

A first-stage \(F < 10\) is the Stock-Yogo rule of thumb for weak instruments. Our DGP uses \(F \approx 5\) — clearly weak.

App 4 — Code

• R
• Python
• Stata

library(AER)

# First-stage F
fs <- lm(x ~ z, data = bs_iv)
cat("First-stage F:", summary(fs)$fstatistic[1], "\n")

# 2SLS point estimate
iv_fit <- ivreg(y ~ x | z, data = bs_iv)

# Pairs bootstrap for 2SLS coefficient
stat_iv <- function(df, idx) {
  coef(ivreg(y ~ x | z, data = df[idx, ]))["x"]
}

set.seed(14159)
bt_iv <- boot(bs_iv, stat_iv, R = B,
              parallel = if (.Platform$OS.type == "unix") "multicore" else "snow",
              ncpus    = 12)

ci_iv <- quantile(bt_iv$t, c(0.025, 0.975))
print(ci_iv)

First-stage F : 1.41 

2SLS estimate : 1.4123  (true β = 1.5 )

Bootstrap 95% percentile CI: [-2.6467, 10.8187]

Bootstrap SE              : 7.9020

import numpy as np, pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed

df_iv = pd.read_csv("../data/bs-iv.csv")
y_, x_, z_ = df_iv["y"].to_numpy(), df_iv["x"].to_numpy(), df_iv["z"].to_numpy()
B_ = 500; N_CORES_ = 12; SEED = 14159

# 2SLS via moment condition: β = Σz'y / Σz'x
def tsls(y, x, z):
    return (z @ y) / (z @ x)

# First-stage F
from scipy.stats import f as fdist
pi_hat  = np.dot(z_, x_) / np.dot(z_, z_)
xhat    = pi_hat * z_
resid_fs= x_ - xhat
F_stat  = (pi_hat**2 * np.dot(z_, z_)) / (np.var(resid_fs, ddof=1))
print(f"First-stage F : {F_stat:.2f}")
print(f"2SLS estimate : {tsls(y_, x_, z_):.4f}  (true β = 1.5)")

def iv_boot(b, seed_base=SEED):
    rng = np.random.default_rng(seed_base + b)
    idx = rng.integers(0, len(y_), len(y_))
    return tsls(y_[idx], x_[idx], z_[idx])

bt_iv = np.array(Parallel(n_jobs=N_CORES_)(delayed(iv_boot)(b) for b in range(B_)))
ci_pct= np.percentile(bt_iv, [2.5, 97.5])
print(f"95% Percentile CI : [{ci_pct[0]:.4f}, {ci_pct[1]:.4f}]")
print(f"Bootstrap SE      : {bt_iv.std(ddof=1):.4f}")

quietly import delimited "../data/bs-iv.csv", clear
quietly destring _all, replace

* First-stage F
regress x z
display "First-stage F: " e(F)

* 2SLS with bootstrap SEs
bootstrap _b[x], reps(500) seed(14159) nodots: ///
    ivregress 2sls y (x = z)

Topic 4 — Time Series Persistence

ἐν τῇ κεφαλῇ γὰρ ἐμμένει πολὺν χρόνον·

it stays in the head a long time

Ἀριστοφάνης, Ἐκκλησιάζουσαι 1120

DGP 4 — Persistent Time Series

Resembles: macroeconomic time series — the stationary AR(1) mimics a deviation from trend (output gap), the near-unit-root mimics inflation or real exchange rates, and the bivariate companion mimics a propagation system like output→consumption.

What is it? Two AR(1) processes with different persistence levels, plus a bivariate companion:

\[y_t^{stationary} = \underbrace{0.80}_{}\,y_{t-1} + u_t, \qquad u_t \overset{iid}{\sim}\mathcal{N}(0,1)\] \[y_t^{near} = \underbrace{0.97}_{}\,y_{t-1} + u_t\] \[y_t^{VAR2} = 0.3\,y_{t-1}^{stationary} + 0.7\,y_{t-1}^{VAR2} + w_t, \qquad w_t \sim \mathcal{N}(0, 0.25)\]

Why it matters for bootstrap: iid resampling destroys serial dependence. All bootstrap results will be wrong unless the resampling scheme respects the time-series structure.

Setting	Correct bootstrap	Key reference
Stationary AR	Block (MBB/CBB) or Sieve	Künsch (1989), Ann. Stat.
Near unit root	Sieve with imposed \(H_0\)	Bühlmann (1997), Bernoulli
Exact unit root	Resample \(\Delta y_t\) under \(H_0\)	Palm et al. (2008)
VAR IRFs	Residual resampling	Kilian (1998), JBES

Kilian (1998) — Why It Matters

Bootstrap CIs for VAR impulse-response functions have dramatically better coverage than delta-method CIs at medium and long horizons. The nonlinear mapping from VAR coefficients to IRF\((h)\) makes the delta-method inaccurate — bootstrap handles this automatically and remains one of the clearest practical wins over asymptotics.

DGP 4 — Code

• R
• Python
• Stata

T_ts <- 400L
rho_s <- 0.80; rho_n <- 0.97

y_ar    <- numeric(T_ts)
y_near  <- numeric(T_ts)
y_var2  <- numeric(T_ts)
y_ar[1] <- rnorm(1, 0, 1/sqrt(1 - rho_s^2))
y_near[1]<- rnorm(1, 0, 1/sqrt(1 - rho_n^2))
y_var2[1]<- rnorm(1)

for (t in 2:T_ts) {
  y_ar[t]   <- rho_s * y_ar[t-1]   + rnorm(1)
  y_near[t] <- rho_n * y_near[t-1] + rnorm(1)
  y_var2[t] <- 0.3*y_ar[t-1] + 0.7*y_var2[t-1] + rnorm(1, 0, 0.5)
}

bs_ts <- tibble(t = 1L:T_ts, y_ar = y_ar,
                y_near = y_near, y_var2 = y_var2)
write_csv(bs_ts, "../data/bs-ts.csv")

import numpy as np, pandas as pd
rng = np.random.default_rng(14159)
T, rho_s, rho_n = 400, 0.80, 0.97
y_ar    = np.zeros(T); y_near = np.zeros(T); y_var2 = np.zeros(T)
y_ar[0] = rng.normal(0, 1/np.sqrt(1-rho_s**2))
y_near[0]= rng.normal(0, 1/np.sqrt(1-rho_n**2))
for t in range(1, T):
    y_ar[t]    = rho_s*y_ar[t-1]    + rng.standard_normal()
    y_near[t]  = rho_n*y_near[t-1]  + rng.standard_normal()
    y_var2[t]  = 0.3*y_ar[t-1] + 0.7*y_var2[t-1] + rng.normal(0, 0.5)
pd.DataFrame({"t": range(1,T+1), "y_ar": y_ar,
              "y_near": y_near, "y_var2": y_var2}
             ).to_csv("../data/bs-ts.csv", index=False)

clear
set seed 14159
set obs 400
gen t = _n
gen y_ar   = 0
gen y_near = 0
gen y_var2 = 0
replace y_ar[1]   = rnormal(0, 1/sqrt(1-0.64))
replace y_near[1] = rnormal(0, 1/sqrt(1-0.9409))
forvalues i = 2/400 {
    replace y_ar[`i']   = 0.80*y_ar[`i'-1]   + rnormal()
    replace y_near[`i'] = 0.97*y_near[`i'-1] + rnormal()
    replace y_var2[`i'] = 0.3*y_ar[`i'-1] + 0.7*y_var2[`i'-1] + rnormal(0,.5)
}
keep t y_ar y_near y_var2
export delimited "../data/bs-ts.csv", replace

DGP 4 — Data

• R
• Python
• Stata

p4d_a <- bs_ts %>%
  select(t, y_ar, y_near) %>%
  pivot_longer(-t, names_to = "series",
               names_transform = list(series = \(x)
                 recode(x, y_ar   = "ρ = 0.80 (stationary)",
                              y_near = "ρ = 0.97 (near unit root)"))) %>%
  ggplot(aes(x = t, y = value, colour = series)) +
  geom_line(linewidth = 0.6, alpha = 0.9) +
  scale_colour_manual(values = c("ρ = 0.80 (stationary)" = col_main,
                                  "ρ = 0.97 (near unit root)" = col_accent),
                      name = NULL) +
  labs(title = "Simulated AR(1) processes",
       x = "t", y = expression(y[t])) +
  theme(legend.position = "top", text = element_text(size = 18))

acf_s <- acf(bs_ts$y_ar, lag.max = 30, plot = FALSE)
p4d_b <- tibble(lag = acf_s$lag[,1,1], acf = acf_s$acf[,1,1]) %>%
  ggplot(aes(x = lag, y = acf)) +
  geom_segment(aes(xend = lag, yend = 0), colour = col_main, linewidth = 0.8) +
  geom_hline(yintercept = c(-1.96, 1.96)/sqrt(nrow(bs_ts)),
             colour = col_accent, linetype = "dashed") +
  geom_hline(yintercept = 0, colour = "grey50") +
  labs(title = "ACF of stationary series  (ρ = 0.80)",
       x = "Lag", y = "Autocorrelation") +
  theme(text = element_text(size = 18))

p4d_a + p4d_b

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf as sm_acf

df_ts4 = pd.read_csv("../data/bs-ts.csv")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5))

ax1.plot(df_ts4["t"], df_ts4["y_ar"],   color="#185FA5", lw=0.6, alpha=0.9, label="ρ = 0.80 (stationary)")
ax1.plot(df_ts4["t"], df_ts4["y_near"], color="#D85A30", lw=0.6, alpha=0.9, label="ρ = 0.97 (near unit root)")
ax1.set_xlabel("t", fontsize=18); ax1.set_ylabel("y_t", fontsize=18)
ax1.set_title("Simulated AR(1) processes", fontsize=18)
ax1.tick_params(labelsize=16)
ax1.legend(loc="upper right", frameon=False, fontsize=16)

acf_vals = sm_acf(df_ts4["y_ar"], nlags=30, fft=True)
lags     = np.arange(len(acf_vals))
ci       = 1.96 / np.sqrt(len(df_ts4))
ax2.vlines(lags, 0, acf_vals, color="#185FA5", lw=1.5)
ax2.axhline( ci, color="#D85A30", ls="--", lw=1.0)
ax2.axhline(-ci, color="#D85A30", ls="--", lw=1.0)
ax2.axhline(0, color="grey", lw=0.5)
ax2.set_xlabel("Lag", fontsize=18); ax2.set_ylabel("Autocorrelation", fontsize=18)
ax2.set_title("ACF of stationary series  (ρ = 0.80)", fontsize=18)
ax2.tick_params(labelsize=16)

plt.tight_layout(); plt.show()

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
tsset t

twoway (line y_ar t,   lcolor("24 95 165")  lwidth(thin)) ///
       (line y_near t, lcolor("216 90 48") lwidth(thin)), ///
   legend(order(1 "ρ = 0.80 (stationary)" 2 "ρ = 0.97 (near unit root)") ///
          position(11) ring(0)) ///
   ytitle("y(t)", size(medlarge)) xtitle("t", size(medlarge)) ///
   title("Simulated AR(1) processes", size(medlarge)) name(g_ts4, replace)

ac y_ar, lags(30) ///
   title("ACF of stationary series  (ρ = 0.80)", size(medlarge)) ///
   name(g_acf4, replace)

graph combine g_ts4 g_acf4, cols(2)

AR(1) Series & ACF

The near-unit-root series (coral) wanders far from zero and shows long swings. The ACF of the stationary series decays geometrically — exactly \(\rho^h\) pattern. Both processes violate the iid assumption: standard bootstrap is invalid.

DGP 4 — Diagnostics

• R
• Python
• Stata

# AR(1) via OLS on lagged series — gives SE and CI
ar_lm_s <- lm(y_ar   ~ lag(y_ar,   1), data = bs_ts)
ar_lm_n <- lm(y_near ~ lag(y_near, 1), data = bs_ts)

cat(sprintf("Stationary  : ρ̂ = %.4f,  SE = %.4f,  95%% CI [%.3f, %.3f]  (true 0.80)\n",
            coef(ar_lm_s)[2], sqrt(diag(vcov(ar_lm_s)))[2],
            confint(ar_lm_s)[2, 1], confint(ar_lm_s)[2, 2]))

Stationary  : ρ̂ = 0.7625,  SE = 0.0327,  95% CI [0.698, 0.827]  (true 0.80)

cat(sprintf("Near unit r.: ρ̂ = %.4f,  SE = %.4f,  95%% CI [%.3f, %.3f]  (true 0.97)\n",
            coef(ar_lm_n)[2], sqrt(diag(vcov(ar_lm_n)))[2],
            confint(ar_lm_n)[2, 1], confint(ar_lm_n)[2, 2]))

Near unit r.: ρ̂ = 0.9725,  SE = 0.0118,  95% CI [0.949, 0.996]  (true 0.97)

# Augmented Dickey-Fuller tests (H0: unit root)
adf_s <- adf.test(bs_ts$y_ar,   alternative = "stationary")
adf_n <- adf.test(bs_ts$y_near, alternative = "stationary")
cat(sprintf("\nADF (stationary)   : t = %.3f,  p = %.4f\n",
            adf_s$statistic, adf_s$p.value))

ADF (stationary)   : t = -5.261,  p = 0.0100

cat(sprintf("ADF (near unit r.) : t = %.3f,  p = %.4f\n",
            adf_n$statistic, adf_n$p.value))

ADF (near unit r.) : t = -2.798,  p = 0.2401

import pandas as pd, numpy as np
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

df = pd.read_csv("../data/bs-ts.csv")

def ar1_fit(arr):
    y = arr[1:]; X = sm.add_constant(arr[:-1])
    return sm.OLS(y, X).fit()

fit_s = ar1_fit(df["y_ar"].to_numpy())
fit_n = ar1_fit(df["y_near"].to_numpy())

ci_s = fit_s.conf_int()[1]; ci_n = fit_n.conf_int()[1]
print(f"Stationary  : ρ̂ = {fit_s.params[1]:.4f}, "
      f"SE = {fit_s.bse[1]:.4f}, 95% CI [{ci_s[0]:.3f}, {ci_s[1]:.3f}]")
print(f"Near unit r.: ρ̂ = {fit_n.params[1]:.4f}, "
      f"SE = {fit_n.bse[1]:.4f}, 95% CI [{ci_n[0]:.3f}, {ci_n[1]:.3f}]")

adf_s = adfuller(df["y_ar"],   maxlag=1, regression="c", autolag=None)
adf_n = adfuller(df["y_near"], maxlag=1, regression="c", autolag=None)
print(f"\nADF (stationary)   : t = {adf_s[0]:.3f},  p = {adf_s[1]:.4f}")
print(f"ADF (near unit r.) : t = {adf_n[0]:.3f},  p = {adf_n[1]:.4f}")

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
tsset t

regress y_ar   L.y_ar
regress y_near L.y_near

dfuller y_ar,   lags(1)
dfuller y_near, lags(1)

ADF Near Unit Root

ADF correctly rejects the unit root for \(\rho = 0.80\) but may fail to reject for \(\rho = 0.97\) despite the series being stationary in truth — the test has poor power near the boundary. Bootstrap critical values can help.

DGP 4 — Bootstrap Setup

Resampling unit: blocks of consecutive observations — moving block bootstrap (MBB), Künsch (1989)

Block length \(\ell\): a bias-variance trade-off. - Too short → fails to capture serial dependence - Too long → too few distinct blocks, high variance

Practical rule: \(\ell^* = 1.75\,T^{1/3}\) → for \(T = 400\), \(\ell \approx 13\).

Parameter	Value	Rationale
\(B\)	500	Time series bootstrap typically needs more reps than iid
Block length	\(\ell = 13\)	Politis-White-Patton rule of thumb
Scheme	sim = "fixed" (MBB)	Also: "geom" (stationary), "scramble" (phase)
Statistic	$^* = $ OLS slope on lagged \(y^*\)	Re-estimate AR(1) on each pseudo-series
Parallel	multicore / snow	12 cores
Seed	14159

Alternatives: stationary bootstrap with \(\mathbb{E}[L] = \ell\) (avoids end-block bias); sieve bootstrap (resample AR-residuals) when the model is correctly specified.

DGP 4 — Bootstrap Run

• R
• Python
• Stata

stat_ar <- function(ts, ...) {
  ar(ts, order.max = 1, method = "ols", aic = FALSE)$ar[[1]]
}

set.seed(14159)
bt_dgp4 <- tsboot(tseries   = bs_ts$y_ar,
                  statistic = stat_ar,
                  R         = 500,
                  l         = 13,
                  sim       = "fixed",
                  parallel  = if (.Platform$OS.type == "unix") "multicore" else "snow",
                  ncpus     = 12)

ci_dgp4 <- quantile(bt_dgp4$t, c(0.025, 0.975))
print(ci_dgp4)

MBB 95% CI for ρ̂  (ℓ = 13, B = 500): [0.6155, 0.7504]

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd
from joblib import Parallel, delayed

df_b4 = pd.read_csv("../data/bs-ts.csv")
y_b4  = df_b4["y_ar"].to_numpy()
T_b4, l_b4 = len(y_b4), 13

def ar1_block_b4(b):
    rng    = np.random.default_rng(14159 + b)
    n_blk  = int(np.ceil(T_b4 / l_b4))
    starts = rng.integers(0, T_b4 - l_b4 + 1, n_blk)
    y_star = np.concatenate([y_b4[s:s+l_b4] for s in starts])[:T_b4]
    Y, X   = y_star[1:], y_star[:-1].reshape(-1, 1)
    return float(np.linalg.lstsq(X, Y, rcond=None)[0][0])

bt_b4 = np.array(Parallel(n_jobs=12)(delayed(ar1_block_b4)(b) for b in range(500)))
ci_b4 = np.percentile(bt_b4, [2.5, 97.5])
print(f"MBB 95% CI for ρ̂  (ℓ = 13, B = 500): [{ci_b4[0]:.4f}, {ci_b4[1]:.4f}]")

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
gen y_ar_lag = y_ar[_n - 1]

bootstrap _b[y_ar_lag], reps(500) seed(14159) nodots: ///
    regress y_ar y_ar_lag

DGP 4 — Bootstrap Preview

• R
• Python
• Stata

y_ar_vec <- bs_ts$y_ar
T_len    <- length(y_ar_vec)
l_blk    <- round(T_len^(1/3) * 1.75)
rho_hat  <- ar(y_ar_vec, order.max = 1, method = "ols", aic = FALSE)$ar[[1]]
se_asy   <- (1 - rho_hat^2) / sqrt(T_len)

one_mbb <- function(b) {
  n_b    <- ceiling(T_len / l_blk)
  starts <- sample.int(T_len - l_blk + 1, n_b, replace = TRUE)
  y_star <- purrr::map(starts, ~ y_ar_vec[.x:(.x + l_blk - 1)]) %>%
              unlist() %>% head(T_len)
  ar(y_star, order.max = 1, method = "ols", aic = FALSE)$ar[[1]]
}

set.seed(14159)
rho_boot_r <- purrr::map_dbl(seq_len(200), one_mbb)

grid    <- seq(rho_hat - 4*se_asy, rho_hat + 4*se_asy, length.out = 300)
overlay <- tibble(x = grid, y = dnorm(grid, rho_hat, se_asy))

ggplot() +
  geom_histogram(data = tibble(r = rho_boot_r),
                 aes(x = r, y = after_stat(density), fill = "Bootstrap"),
                 bins = 30, alpha = 0.55, colour = "white") +
  geom_line(data = overlay,
            aes(x = x, y = y, colour = "Asymptotic N(ρ̂, SE²)"), linewidth = 1.1) +
  geom_vline(aes(xintercept = RHO_AR, colour = "True ρ"),
             linetype = "dashed", linewidth = 0.9, key_glyph = "vline") +
  scale_fill_manual(name = NULL, values = c("Bootstrap" = col_main)) +
  scale_colour_manual(name = NULL,
    values = c("Asymptotic N(ρ̂, SE²)" = col_accent, "True ρ" = col_ok)) +
  labs(title = sprintf("Block-bootstrap distribution of ρ̂*  (B = 200, ℓ = %d)", l_blk),
       x = expression(hat(rho)^"*"), y = "Density") +
  theme(legend.position = "bottom", text = element_text(size = 18))

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm
from joblib import Parallel, delayed

df_p4 = pd.read_csv("../data/bs-ts.csv")
y_p4  = df_p4["y_ar"].to_numpy()
T_p4, l_p4 = len(y_p4), round(len(y_p4)**(1/3) * 1.75)

def ar1_coef_p4(y):
    Y = y[1:]; X = y[:-1].reshape(-1, 1)
    return float(np.linalg.lstsq(X, Y, rcond=None)[0][0])

def mbb_p4(b):
    rng = np.random.default_rng(14159 + b)
    starts = rng.integers(0, T_p4 - l_p4 + 1, int(np.ceil(T_p4 / l_p4)))
    y_star = np.concatenate([y_p4[s:s+l_p4] for s in starts])[:T_p4]
    return ar1_coef_p4(y_star)

rho_boot_py = np.array(Parallel(n_jobs=12)(delayed(mbb_p4)(b) for b in range(200)))
rho_hat_py  = ar1_coef_p4(y_p4)
se_asy_py   = (1 - rho_hat_py**2) / np.sqrt(T_p4)

fig, ax = plt.subplots(figsize=(9, 4.5))
ax.hist(rho_boot_py, bins=30, density=True, alpha=0.55, color="#185FA5", label="Bootstrap")
grid_py = np.linspace(rho_hat_py - 4*se_asy_py, rho_hat_py + 4*se_asy_py, 300)
ax.plot(grid_py, norm.pdf(grid_py, rho_hat_py, se_asy_py),
        color="#D85A30", lw=1.4, label="Asymptotic N(ρ̂, SE²)")
ax.axvline(0.80, color="#1D9E75", ls="--", lw=1.2, label="True ρ")
ax.set_xlabel(r"$\hat{\rho}^*$", fontsize=18)
ax.set_ylabel("Density", fontsize=18)
ax.set_title(f"Block-bootstrap distribution of ρ̂*  (B = 200, ℓ = {l_p4})", fontsize=18)
ax.tick_params(labelsize=16)
ax.legend(loc="upper left", frameon=False, fontsize=16)
plt.tight_layout(); plt.show()

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
gen y_ar_lag = y_ar[_n - 1]

simulate coef=_b[y_ar_lag], reps(200) seed(14159) nodots: ///
    regress y_ar y_ar_lag

histogram coef, normal width(0.005) ///
    fcolor("24 95 165%45") lcolor(none) ///
    normlopts(lcolor("216 90 48") lwidth(medthick)) ///
    title("Block-bootstrap distribution of ρ̂*  (B = 200)", size(medlarge)) ///
    xtitle("ρ̂*", size(medlarge)) ytitle("Density", size(medlarge))

ρ = 0.97 — Near Unit Root

For \(\rho = 0.80\), the asymptotic normal is a decent approximation. For \(\rho = 0.97\) (near unit root), the bootstrap distribution becomes left-skewed — asymptotic theory is misleading. Full analysis in Applications TS1 & TS2.

App 5 — Nonstandard Distributions: Unit Root

Problem: under \(H_0: \rho = 1\), the ADF test statistic has a Dickey-Fuller distribution — not \(t_{T-2}\) — and critical values are tabulated only for standard cases.

ADF regression: \[\Delta y_t = \mu + \phi\, y_{t-1} + \sum_{j=1}^p \gamma_j \Delta y_{t-j} + u_t\]

Test statistic: \(\hat\tau = \hat\phi / \widehat{SE}(\hat\phi)\)

Under \(H_0: \phi = 0\): \[\hat\tau \xrightarrow{d} \frac{\int_0^1 W(r)\,dW(r)}{\left[\int_0^1 W(r)^2\,dr\right]^{1/2}}\]

where \(W(r)\) is a standard Brownian motion — not normal.

Bootstrap ADF resamples \(\hat u_t\) to generate pseudo-series under \(H_0\) (unit root imposed), computing \(\hat\tau^*_b\) for each replication and extracting \(p\)-values.

When Bootstrap ADF Helps

The bootstrap ADF test is especially useful when:

• Error distribution is non-normal
• Structural breaks are present
• Conditional heteroskedasticity (GARCH) distorts tabulated critical values

tseries::adf.test() uses MacKinnon (1994) tabulated CVs. The bootstrap alternative directly simulates the null distribution.

DGP Boundary Problem

Our DGP has \(\rho_{AR} = 0.80\) (clearly stationary) and \(\rho_{near} = 0.97\) (near unit root). The ADF test has lower power near the boundary — the bootstrap can be more informative here.

App 5 — Code

• R
• Python
• Stata

library(tseries)

# ─── Standard ADF tests (asymptotic, with MacKinnon critical values) ─────────
adf_ar   <- adf.test(bs_ts$y_ar,   alternative = "stationary")
adf_near <- adf.test(bs_ts$y_near, alternative = "stationary")

# ─── Bootstrap ADF, in three small named steps ──────────────────────────────

# Step 1: compute the ADF t-statistic for any time series
adf_t_stat <- function(y) {
  fit <- lm(diff(y) ~ y[-length(y)])
  coef(summary(fit))[2, "t value"]
}

# Step 2: simulate ONE random walk under H0 using the data's own innovations
#         (demeaned first-differences are the empirical residuals under H0)
one_h0_series <- function(y, innov_pool) {
  u_star <- sample(innov_pool, length(innov_pool), replace = TRUE)
  cumsum(c(y[1], u_star))                # y_t = y_{t-1} + u_t
}

# Step 3: assemble the bootstrap p-value
bootstrap_adf <- function(y, B) {
  innov_pool <- diff(y) - mean(diff(y))
  tau_obs    <- adf_t_stat(y)

  # B bootstrap replicates of the ADF statistic, all under H0
  set.seed(14159)
  tau_star <- purrr::map_dbl(seq_len(B),
                             ~ adf_t_stat(one_h0_series(y, innov_pool)))

  list(tau_obs   = tau_obs,
       p_boot    = mean(tau_star <= tau_obs),     # one-sided lower tail
       crit_vals = quantile(tau_star, c(.01, .05, .10)))
}

res_ar   <- bootstrap_adf(bs_ts$y_ar,   B = 500)
res_near <- bootstrap_adf(bs_ts$y_near, B = 500)

Bootstrap ADF (ρ = 0.80): τ = -7.262,  p* = 0.0000

Bootstrap ADF (ρ = 0.97): τ = -2.334,  p* = 0.1800

import numpy as np, pandas as pd
from statsmodels.tsa.stattools import adfuller
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed

df_ts  = pd.read_csv("../data/bs-ts.csv")
y_ar_  = df_ts["y_ar"].to_numpy()
y_nr_  = df_ts["y_near"].to_numpy()
B_ = 500; NC = 12; SEED = 14159

def adf_stat_py(y):
    res = adfuller(y, maxlag=1, regression="n", autolag=None)
    return res[0]   # ADF t-statistic

def boot_adf_py(y, B_b=500, n_jobs=12, seed=14159):
    tau_obs = adf_stat_py(y)
    u_hat   = np.diff(y) - np.diff(y).mean()
    T_      = len(u_hat)
    def one_rep(b):
        rng    = np.random.default_rng(seed + b)
        u_star = rng.choice(u_hat, T_, replace=True)
        y_star = np.cumsum(np.concatenate([[y[0]], u_star]))
        return adf_stat_py(y_star)
    taus = np.array(Parallel(n_jobs=n_jobs)(delayed(one_rep)(b) for b in range(B_b)))
    return {"tau_obs": tau_obs, "p_boot": np.mean(taus <= tau_obs),
            "cv_1": np.percentile(taus, 1),
            "cv_5": np.percentile(taus, 5)}

res_ar_ = boot_adf_py(y_ar_, B_, NC, SEED)
res_nr_ = boot_adf_py(y_nr_, B_, NC, SEED)
print(f"y_ar  (ρ=0.80): τ={res_ar_['tau_obs']:.3f}  p*={res_ar_['p_boot']:.4f}")
print(f"y_near(ρ=0.97): τ={res_nr_['tau_obs']:.3f}  p*={res_nr_['p_boot']:.4f}")

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
quietly tsset t

dfuller y_ar,   lags(1) noconstant
dfuller y_near, lags(1) noconstant

TS 1 — Block Bootstrap for AR Model: Theory

Setting: We have a stationary AR(1) \(y_t = \rho y_{t-1} + u_t\) and want to bootstrap SE(\(\hat\rho\)) without assuming iid errors.

Moving Block Bootstrap (MBB) algorithm:

1. Choose block length \(\ell\) (e.g., \(\ell = \lfloor T^{1/3} \rfloor = 7\) for \(T=400\))
2. Form \(T - \ell + 1\) overlapping blocks: \(\mathcal{B}_k = (y_k, \ldots, y_{k+\ell-1})\)
3. Draw \(m = \lceil T/\ell \rceil\) blocks with replacement → \(\mathbf{y}^*\) of length \(m\ell\)
4. Trim to length \(T\); estimate \(\hat\rho^*_b\) from \(\mathbf{y}^*\)
5. Repeat \(B\) times → \(\{\hat\rho^*_b\}\)

Circular Block Bootstrap (CBB): wraps the series to form \(T\) starting points. Removes the end-block bias of MBB — preferred in practice.

Optimal block length: \[\ell^* = 1.75\, T^{1/3} \left(\frac{2\hat\rho^2}{1-\hat\rho^2}\right)^{2/9}\]

tsboot() in R (boot package)

tsboot(tseries, statistic,
       R = B,
       l = block_len,
       sim = "fixed",    # MBB
       # sim = "geom",   # Stationary bootstrap
       # sim = "scramble"# Phase scramble
       parallel = "multicore",
       ncpus = 12)

Sieve Bootstrap

Block bootstrap for AR parameters is asymptotically valid but may be less efficient than the sieve bootstrap (AR residual resampling) when the AR order is known. The sieve bootstrap re-estimates AR residuals directly and resamples them — preserving model structure without arbitrary block length choice.

TS 1 — Block Bootstrap: Code

• R
• Python
• Stata

# AR(1) coefficient bootstrap via tsboot (circular block)
y_ar_ts   <- bs_ts$y_ar
T_ts_len  <- length(y_ar_ts)
block_len <- round(T_ts_len^(1/3) * 1.75)   # optimal length ≈ 13

# Statistic: AR(1) coefficient
stat_ar <- function(tseries, ...) {
  fit <- ar(tseries, order.max = 1, method = "ols", aic = FALSE)
  fit$ar[[1]]
}

set.seed(14159)
bt_ar_mbb <- tsboot(
  tseries   = y_ar_ts,
  statistic = stat_ar,
  R         = B,
  l         = block_len,
  sim       = "fixed",      # MBB (fixed-length blocks)
  parallel  = if (.Platform$OS.type == "unix") "multicore" else "snow",
  ncpus     = 12
)

rho_hat <- ar(y_ar_ts, order.max=1, method="ols", aic=FALSE)$ar[[1]]
ci_mbb  <- quantile(bt_ar_mbb$t, c(0.025, 0.975))

cat(sprintf("AR(1) estimate: %.4f  (true ρ = %.2f)\n", rho_hat, RHO_AR))
cat(sprintf("MBB 95%% CI    : [%.4f, %.4f]\n", ci_mbb[1], ci_mbb[2]))

Block length    : 13

AR(1) estimate  : 0.7625  (true ρ = 0.80)

Block-bootstrap CI: [0.6155, 0.7504]

import numpy as np, pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed

df_ts_ = pd.read_csv("../data/bs-ts.csv")
y_ar_  = df_ts_["y_ar"].to_numpy()
T_     = len(y_ar_)
l_     = round(T_**(1/3) * 1.75)   # block length
B_ = 500; NC_ = 12; SEED = 14159

def ar1_coef(y):
    Y = y[1:]; X = y[:-1].reshape(-1,1)
    return float(np.linalg.lstsq(X, Y, rcond=None)[0][0])   # [0][0]: array→element→float

def mbb_one(b, seed_base=SEED):
    rng     = np.random.default_rng(seed_base + b)
    n_blk   = int(np.ceil(T_ / l_))
    starts  = rng.integers(0, T_ - l_ + 1, n_blk)
    blocks  = [y_ar_[s:s+l_] for s in starts]
    y_star  = np.concatenate(blocks)[:T_]
    return ar1_coef(y_star)

boot_rho = np.array(Parallel(n_jobs=NC_)(delayed(mbb_one)(b) for b in range(B_)))
rho_hat_ = ar1_coef(y_ar_)
ci_pct_  = np.percentile(boot_rho, [2.5, 97.5])

print(f"Block length  : {l_}")
print(f"AR(1) estimate: {rho_hat_:.4f}  (true ρ = 0.80)")
print(f"MBB 95% CI    : [{ci_pct_[0]:.4f}, {ci_pct_[1]:.4f}]")
print(f"Bootstrap SE  : {boot_rho.std(ddof=1):.4f}")

quietly import delimited "../data/bs-ts.csv", clear
quietly destring _all, replace
gen y_ar_lag = y_ar[_n-1]
regress y_ar y_ar_lag
display "AR(1) estimate: " %6.4f _b[y_ar_lag]

bootstrap _b[y_ar_lag], reps(500) seed(14159) nodots: ///
    regress y_ar y_ar_lag

TS 2 — VAR Impulse-Response Bootstrap: Theory

Motivation: Impulse-response functions (IRFs) from a VAR summarise the dynamic response of one variable to a shock in another. Standard asymptotic CIs for IRFs are often too narrow — bootstrap CIs have better coverage.

VAR(\(p\)) model:

\[\mathbf{y}_t = \boldsymbol{c} + \mathbf{A}_1\mathbf{y}_{t-1} + \cdots + \mathbf{A}_p\mathbf{y}_{t-p} + \mathbf{u}_t, \quad \mathbf{u}_t \sim (0, \Sigma)\]

Bootstrap algorithm (residual resampling):

1. Fit VAR(\(p\)) → \(\hat{\mathbf{A}}_j\), residuals \(\hat{\mathbf{u}}_t\)
2. Resample residuals with replacement: \(\mathbf{u}_t^*\)
3. Reconstruct \(\mathbf{y}^*\) using \(\hat{\mathbf{A}}_j\) and \(\mathbf{u}_t^*\)
4. Re-estimate VAR → \(\hat{\mathbf{A}}_j^*\) → compute IRF\(^*_b(h)\)
5. Repeat \(B\) times → percentile CI on each horizon \(h\)

Cholesky shock to first variable: \[\boldsymbol{\varepsilon}_t = \mathbf{P}^{-1}\mathbf{u}_t, \quad \mathbf{P} = \text{chol}(\hat\Sigma)\]

Why Bootstrap CIs for IRFs?

Asymptotic delta-method CIs for IRFs assume joint normality of the VAR coefficients. At longer horizons, the nonlinear mapping from \(\hat{\mathbf{A}}\) to IRF\((h)\) makes the delta-method inaccurate. Bootstrap handles this automatically.

Block vs Residual Bootstrap

Block bootstrap (not residual) is needed if \(\hat{\mathbf{u}}_t\) exhibit serial correlation after fitting — rare for a well-specified VAR but can arise with misspecification.

TS 2 — VAR IRF Bootstrap: Code

• R
• Python

library(vars)

# Bivariate VAR on (y_ar, y_var2)
var_data  <- bs_ts %>% select(y_ar, y_var2)
var_fit   <- VAR(var_data, p = 1, type = "const")
summary(var_fit)

# Bootstrap IRFs (residual resampling, B replications)
set.seed(14159)
irf_boot  <- irf(var_fit,
                 impulse  = "y_ar",
                 response = "y_var2",
                 n.ahead  = 10,
                 boot     = TRUE,
                 runs     = B,
                 ci       = 0.95)
plot(irf_boot)

import numpy as np, pandas as pd
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed
from statsmodels.tsa.vector_ar.var_model import VAR as smVAR

df_ts_ = pd.read_csv("../data/bs-ts.csv")
Y      = df_ts_[["y_ar", "y_var2"]].to_numpy()
T_     = len(Y); H = 10; B_ = 500; NC_ = 12; SEED = 14159

var_fit_  = smVAR(Y).fit(maxlags=1, trend="c")
resid_var = var_fit_.resid

def irf_from_coef(A1, c, resid, H_):
    """Cholesky IRF from y_ar to y_var2 at all horizons"""
    Sigma = resid.T @ resid / len(resid)
    P = np.linalg.cholesky(Sigma)
    # Companion form IRF
    irf_h = np.zeros(H_+1)
    Phi   = np.eye(2)
    for h in range(H_+1):
        irf_h[h] = Phi[1, 0] * P[0, 0]   # response of y_var2 to shock in y_ar
        Phi = Phi @ A1
    return irf_h

A1_hat   = var_fit_.coefs[0]
irf_obs  = irf_from_coef(A1_hat, None, resid_var, H)

def one_boot(b, seed_base=SEED):
    rng     = np.random.default_rng(seed_base + b)
    n_resid = len(resid_var)                           # = T_ - 1
    idx     = rng.integers(0, n_resid, n_resid)
    u_star  = resid_var[idx]                           # shape (T_-1, 2)
    Y_star  = np.zeros((T_, 2))
    Y_star[0] = Y[0]
    const_  = var_fit_.params[0]
    for t in range(1, T_):
        Y_star[t] = const_ + A1_hat @ Y_star[t-1] + u_star[t-1]
    try:
        fit_b = smVAR(Y_star).fit(maxlags=1, trend="c")
        return irf_from_coef(fit_b.coefs[0], None, fit_b.resid, H)
    except Exception:
        return np.full(H+1, np.nan)

irf_boot = np.array(Parallel(n_jobs=NC_)(delayed(one_boot)(b) for b in range(B_)))
irf_boot = irf_boot[~np.isnan(irf_boot).any(axis=1)]
lo95 = np.percentile(irf_boot, 2.5, axis=0)
hi95 = np.percentile(irf_boot, 97.5, axis=0)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 4))
ax.fill_between(range(H+1), lo95, hi95, alpha=0.20, color="#185FA5")
ax.plot(irf_obs, color="#185FA5", lw=2, label="IRF")
ax.plot(lo95, color="grey", lw=0.8, ls="--")
ax.plot(hi95, color="grey", lw=0.8, ls="--")
ax.axhline(0, color="#D85A30", ls="--")
ax.set_xlabel("Horizon", fontsize=18); ax.set_ylabel("Response", fontsize=18)
ax.set_title("IRF y_ar → y_var2  (Bootstrap 95% CI, B = 500)", fontsize=18)
ax.tick_params(labelsize=16)
ax.legend(frameon=False, fontsize=16); plt.tight_layout(); plt.show()

Topic 5 — Volatility Clustering

καὶ πνεύματʼ ἀνέμων οὐκ ἀεὶ ῥώμην ἔχει·

the blasts of the winds do not keep their force for ever

Εὐριπίδης, Ἡρακλῆς 102

DGP 5 — GARCH Volatility Process

Resembles: daily financial returns — equity index, FX rate, or commodity prices. Calm periods cluster together, large absolute returns also cluster, and unconditional variance is finite but the conditional variance fluctuates.

What is it? A GARCH(1,1) return series with conditionally heteroskedastic innovations:

\[r_t = \sigma_t\, z_t, \quad z_t \overset{iid}{\sim} \mathcal{N}(0,1)\] \[\sigma_t^2 = \underbrace{0.05}_{\omega} + \underbrace{0.15}_{\alpha}\,r_{t-1}^2 + \underbrace{0.80}_{\beta}\,\sigma_{t-1}^2\]

Unconditional variance: \[\mathbb{E}[\sigma_t^2] = \frac{\omega}{1 - \alpha - \beta} = \frac{0.05}{1 - 0.95} = 1.0 \qquad (\alpha + \beta = 0.95 < 1 \Rightarrow \text{stationary})\]

Why it matters for bootstrap: GARCH errors are martingale differences but not iid — standard bootstrap fails because it destroys the volatility clustering. Filtered Historical Simulation (FHS) accounts for this.

Reference	Contribution
Engle (1982), Econometrica	ARCH model introduced
Bollerslev (1986), JOE	GARCH generalisation
Barone-Adesi & Giannopoulos (1996)	Filtered Historical Simulation
Gonçalves & Kilian (2004), JoE	Bootstrap for GARCH — standard residual resampling is invalid; wild bootstrap needed

Why Plain Bootstrap Fails for GARCH

Resampling \(r_t^*\) iid destroys volatility clustering. The bootstrap distribution of \(\hat\alpha^*\) will be biased toward zero. FHS resamples standardised residuals \(\hat{z}_t = r_t/\hat\sigma_t\) instead, preserving the GARCH dynamics.

DGP 5 — Code

• R
• Python
• Stata

T_g <- 500L
omega <- 0.05; alpha_g <- 0.15; beta_g <- 0.80
h     <- numeric(T_g); r_gch <- numeric(T_g)
h[1]     <- omega / (1 - alpha_g - beta_g)
r_gch[1] <- rnorm(1, 0, sqrt(h[1]))
for (t in 2:T_g) {
  h[t]     <- omega + alpha_g * r_gch[t-1]^2 + beta_g * h[t-1]
  r_gch[t] <- rnorm(1, 0, sqrt(h[t]))
}
bs_garch <- tibble(t = 1L:T_g, r = r_gch, h_true = h)
write_csv(bs_garch, "../data/bs-garch.csv")

import numpy as np, pandas as pd
rng = np.random.default_rng(14159)
T, omega, alpha, beta = 500, 0.05, 0.15, 0.80
h = np.zeros(T); r = np.zeros(T)
h[0] = omega/(1-alpha-beta)
r[0] = rng.normal(0, np.sqrt(h[0]))
for t in range(1, T):
    h[t] = omega + alpha*r[t-1]**2 + beta*h[t-1]
    r[t] = rng.normal(0, np.sqrt(h[t]))
pd.DataFrame({"t": range(1,T+1), "r": r, "h_true": h}
             ).to_csv("../data/bs-garch.csv", index=False)

clear
set seed 14159
set obs 500
gen t      = _n
gen h      = 0
gen r      = 0
replace h[1] = 0.05/(1-0.15-0.80)
replace r[1] = rnormal(0, sqrt(h[1]))
forvalues i = 2/500 {
    replace h[`i'] = 0.05 + 0.15*r[`i'-1]^2 + 0.80*h[`i'-1]
    replace r[`i'] = rnormal(0, sqrt(h[`i']))
}
rename h h_true
keep t r h_true
export delimited "../data/bs-garch.csv", replace

DGP 5 — Data

• R
• Python
• Stata

p5d_a <- ggplot(bs_garch, aes(x = t, y = r)) +
  geom_line(colour = col_main, linewidth = 0.4, alpha = 0.8) +
  labs(title = "Simulated returns: volatility clustering",
       subtitle = "Calm and turbulent periods alternate",
       x = "t", y = expression(r[t])) +
  theme(text = element_text(size = 18))

acf_abs <- acf(abs(bs_garch$r), lag.max = 30, plot = FALSE)
p5d_b <- tibble(lag = acf_abs$lag[,1,1], acf = acf_abs$acf[,1,1]) %>%
  ggplot(aes(x = lag, y = acf)) +
  geom_segment(aes(xend = lag, yend = 0), colour = col_main, linewidth = 0.8) +
  geom_hline(yintercept = c(-1.96, 1.96)/sqrt(nrow(bs_garch)),
             colour = col_accent, linetype = "dashed") +
  geom_hline(yintercept = 0, colour = "grey50") +
  labs(title = "ACF of |r_t|  (volatility persistence)",
       x = "Lag", y = "Autocorrelation") +
  theme(text = element_text(size = 18))

p5d_a + p5d_b

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf as sm_acf

df_g5 = pd.read_csv("../data/bs-garch.csv")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5))

ax1.plot(df_g5["t"], df_g5["r"], color="#185FA5", lw=0.4, alpha=0.8)
ax1.set_xlabel("t", fontsize=18); ax1.set_ylabel("r_t", fontsize=18)
ax1.set_title("Simulated returns: volatility clustering", fontsize=18)
ax1.tick_params(labelsize=16)

acf_vals = sm_acf(np.abs(df_g5["r"].to_numpy()), nlags=30, fft=True)
lags = np.arange(len(acf_vals))
ci   = 1.96 / np.sqrt(len(df_g5))
ax2.vlines(lags, 0, acf_vals, color="#185FA5", lw=1.5)
ax2.axhline( ci, color="#D85A30", ls="--", lw=1.0)
ax2.axhline(-ci, color="#D85A30", ls="--", lw=1.0)
ax2.axhline(0, color="grey", lw=0.5)
ax2.set_xlabel("Lag", fontsize=18); ax2.set_ylabel("Autocorrelation", fontsize=18)
ax2.set_title("ACF of |r_t|  (volatility persistence)", fontsize=18)
ax2.tick_params(labelsize=16)

plt.tight_layout(); plt.show()

quietly import delimited "../data/bs-garch.csv", clear
quietly destring _all, replace
tsset t

twoway (line r t, lcolor("24 95 165") lwidth(thin)), ///
   title("Simulated returns: volatility clustering", size(medlarge)) ///
   ytitle("r(t)", size(medlarge)) xtitle("t", size(medlarge)) ///
   name(g_r5, replace)

gen abs_r = abs(r)
ac abs_r, lags(30) ///
   title("ACF of |r_t|  (volatility persistence)", size(medlarge)) ///
   name(g_acf5, replace)

graph combine g_r5 g_acf5, cols(2)

Volatility Clustering & ACF

Returns \(r_t\) are uncorrelated (efficient market), but \(|r_t|\) shows strong autocorrelation — large moves predict more large moves. This is the volatility clustering signature that GARCH models capture.

DGP 5 — Diagnostics

• R
• Python
• Stata

r_vec <- bs_garch$r

cat(sprintf("Sample mean    : %.4f  (true 0)\n",      mean(r_vec)))

Sample mean    : -0.0065  (true 0)

cat(sprintf("Sample sd      : %.4f  (true 1.0)\n",    sd(r_vec)))

Sample sd      : 1.1849  (true 1.0)

cat(sprintf("Sample kurtosis: %.3f  (Normal = 3.0)\n",
            sum((r_vec - mean(r_vec))^4) / (length(r_vec) * sd(r_vec)^4)))

Sample kurtosis: 5.662  (Normal = 3.0)

# Ljung-Box on r_t  (H0: no autocorrelation; should NOT reject)
lb_r  <- Box.test(r_vec,    lag = 10, type = "Ljung-Box")
# Ljung-Box on r_t^2 (SHOULD reject — squared returns are correlated)
lb_r2 <- Box.test(r_vec^2,  lag = 10, type = "Ljung-Box")
cat(sprintf("\nLjung-Box on r_t   : Q = %.2f,  p = %.4f  (uncorrelated)\n",
            lb_r$statistic, lb_r$p.value))

Ljung-Box on r_t   : Q = 32.60,  p = 0.0003  (uncorrelated)

cat(sprintf("Ljung-Box on r_t^2 : Q = %.2f,  p = %.4f  (correlated → ARCH)\n",
            lb_r2$statistic, lb_r2$p.value))

Ljung-Box on r_t^2 : Q = 263.91,  p = 0.0000  (correlated → ARCH)

import warnings; warnings.filterwarnings("ignore")
import pandas as pd, numpy as np
from scipy.stats import kurtosis
from statsmodels.stats.diagnostic import acorr_ljungbox

df_d5 = pd.read_csv("../data/bs-garch.csv")
r_d5  = df_d5["r"].to_numpy()

print(f"Sample mean    : {r_d5.mean():.4f}  (true 0)")
print(f"Sample sd      : {r_d5.std(ddof=1):.4f}  (true 1.0)")
print(f"Sample kurtosis: {kurtosis(r_d5, fisher=False):.3f}  (Normal = 3.0)")

lb_r  = acorr_ljungbox(r_d5,     lags=[10], return_df=True).iloc[0]
lb_r2 = acorr_ljungbox(r_d5**2,  lags=[10], return_df=True).iloc[0]
print(f"\nLjung-Box on r_t   : Q = {lb_r['lb_stat']:.2f},  p = {lb_r['lb_pvalue']:.4f}")
print(f"Ljung-Box on r_t^2 : Q = {lb_r2['lb_stat']:.2f},  p = {lb_r2['lb_pvalue']:.4f}  (→ ARCH)")

quietly import delimited "../data/bs-garch.csv", clear
quietly destring _all, replace

summarize r, detail
wntestq r,  lags(10)
gen r2 = r^2
wntestq r2, lags(10)

ARCH Fingerprint

Excess kurtosis (\(>3\)) and significant autocorrelation in \(r_t^2\) — the classic ARCH fingerprint. Standard iid bootstrap would destroy the volatility dependence; filtered historical simulation (FHS) preserves it.

DGP 5 — Bootstrap Setup

Resampling unit: standardised residuals \(\hat z_t = r_t / \hat\sigma_t\) — Filtered Historical Simulation (Barone-Adesi & Giannopoulos, 1996)

Algorithm:

1. Fit GARCH(1,1) → recover \(\hat\sigma_t\) and \(\hat z_t = r_t / \hat\sigma_t\)
2. Resample \(\{\hat z_t\}\) with replacement → \(\{z^*_t\}\)
3. Generate \(r^*_t = \hat\sigma_t \cdot z^*_t\) — uses the fitted volatility path
4. Re-fit GARCH on \(\{r^*_t\}\), extract \((\omega^*, \alpha^*, \beta^*)\)
5. Repeat \(B\) times

Parameter	Value	Rationale
\(B\)	500	Larger if accurate tail CIs needed (e.g. VaR)
Resample	standardised residuals \(\hat z_t\)	Captures empirical kurtosis without Gaussian assumption
Statistic	\((\hat\omega^*, \hat\alpha^*, \hat\beta^*)\)	Three parameters jointly
Stationarity check	drop replicates with \(\hat\alpha^* + \hat\beta^* \ge 1\)	Or project onto the constraint
Parallel	multicore / snow	GARCH refit is expensive: 12 cores cut time substantially
Seed	14159

Why not naive bootstrap? iid resampling of \(r_t\) destroys volatility clustering and biases \(\hat\alpha\) toward zero (Gonçalves & Kilian, 2004).

DGP 5 — Bootstrap Run

• R
• Python
• Stata

library(rugarch)

spec <- ugarchspec(
  variance.model     = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model         = list(armaOrder = c(0, 0), include.mean = FALSE),
  distribution.model = "norm"
)
garch_fit <- ugarchfit(spec, data = bs_garch$r)
show(garch_fit)

set.seed(14159)
bt_garch <- ugarchboot(garch_fit, method = "Partial",
                       sampling   = "raw",
                       n.bootfit  = B,
                       n.bootpred = B)
show(bt_garch)

omega = 0.0934  (true 0.05)
alpha = 0.2748  (true 0.15)
beta  = 0.6695  (true 0.80)
alpha+beta = 0.9443  (< 1 = stationary)

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd
from arch import arch_model
from joblib import Parallel, delayed

df_b5 = pd.read_csv("../data/bs-garch.csv")
r_b5  = df_b5["r"].to_numpy()

res_b5 = arch_model(r_b5, vol="Garch", p=1, q=1, dist="Normal").fit(disp="off")
print(f"omega = {res_b5.params['omega']:.4f}  (true 0.05)")
print(f"alpha = {res_b5.params['alpha[1]']:.4f}  (true 0.15)")
print(f"beta  = {res_b5.params['beta[1]']:.4f}  (true 0.80)")

z_b5   = res_b5.resid[~np.isnan(res_b5.resid)] / res_b5.conditional_volatility[~np.isnan(res_b5.resid)]
sig_b5 = res_b5.conditional_volatility[~np.isnan(res_b5.conditional_volatility)]

def fhs_b5(b):
    rng    = np.random.default_rng(14159 + b)
    z_star = rng.choice(z_b5, len(z_b5), replace=True)
    r_star = sig_b5 * z_star
    try:
        fit = arch_model(r_star, vol="Garch", p=1, q=1).fit(disp="off")
        return [fit.params["omega"], fit.params["alpha[1]"], fit.params["beta[1]"]]
    except Exception:
        return [np.nan, np.nan, np.nan]

bt_b5 = np.array(Parallel(n_jobs=12)(delayed(fhs_b5)(b) for b in range(500)))
bt_b5 = bt_b5[~np.isnan(bt_b5).any(axis=1)]
print(f"\nBootstrap 95% CI  (B = {len(bt_b5)} valid):")
for i, name in enumerate(["omega", "alpha", "beta"]):
    lo, hi = np.percentile(bt_b5[:, i], [2.5, 97.5])
    print(f"  {name}: [{lo:.4f}, {hi:.4f}]")

quietly import delimited "../data/bs-garch.csv", clear
quietly destring _all, replace

arch r, arch(1) garch(1) nolog

DGP 5 — Bootstrap Preview

• R
• Python
• Stata

band_df <- bs_garch %>%
  mutate(upper = 2 * sqrt(h_true),
         lower = -2 * sqrt(h_true)) %>%
  tidyr::pivot_longer(c(upper, lower), names_to = "band", values_to = "value")

p5a <- ggplot() +
  geom_line(data = bs_garch, aes(x = t, y = r, colour = "Return r_t"),
            linewidth = 0.4, alpha = 0.8) +
  geom_line(data = band_df,
            aes(x = t, y = value, colour = "±2·σ_t", group = band),
            linewidth = 0.9) +
  scale_colour_manual(name = NULL,
    values = c("Return r_t" = col_muted, "±2·σ_t" = col_accent)) +
  labs(title = "Returns with ±2σ_t bands", x = "t", y = expression(r[t])) +
  theme(legend.position = "bottom", text = element_text(size = 18))

qq_df <- tibble(
  theoretical = qnorm(ppoints(nrow(bs_garch))),
  sample      = sort(bs_garch$r / sqrt(bs_garch$h_true))
)
p5b <- ggplot(qq_df, aes(x = theoretical, y = sample)) +
  geom_point(aes(colour = "Standardised z_t"), alpha = 0.45, size = 1.4) +
  geom_abline(aes(slope = 1, intercept = 0, colour = "45° reference"), linewidth = 1) +
  scale_colour_manual(name = NULL,
    values = c("Standardised z_t" = col_main, "45° reference" = col_accent)) +
  labs(title = "Q-Q plot of z_t = r_t/σ_t",
       x = "Theoretical N(0,1) quantiles", y = "Sample quantiles") +
  theme(legend.position = "bottom", text = element_text(size = 18))

p5a + p5b

import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm

df_p5    = pd.read_csv("../data/bs-garch.csv")
t_p5     = df_p5["t"].to_numpy()
r_p5     = df_p5["r"].to_numpy()
sigma_t5 = np.sqrt(df_p5["h_true"].to_numpy())

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

axes[0].plot(t_p5, r_p5, color="#888888", lw=0.5, alpha=0.8, label="Return r_t")
axes[0].plot(t_p5,  2 * sigma_t5, color="#D85A30", lw=1.2, label="±2·σ_t")
axes[0].plot(t_p5, -2 * sigma_t5, color="#D85A30", lw=1.2)
axes[0].set_xlabel("t", fontsize=18); axes[0].set_ylabel("r_t", fontsize=18)
axes[0].set_title("Returns with ±2σ_t bands", fontsize=18)
axes[0].tick_params(labelsize=16)
axes[0].legend(frameon=False, fontsize=16)

z_qq  = np.sort(r_p5 / sigma_t5)
probs = (np.arange(1, len(z_qq) + 1) - 0.5) / len(z_qq)
theo  = norm.ppf(probs)
axes[1].plot(theo, z_qq, "o", color="#185FA5", alpha=0.4, ms=4, label="Standardised z_t")
axes[1].plot([theo.min(), theo.max()], [theo.min(), theo.max()],
             color="#D85A30", lw=1.2, label="45° reference")
axes[1].set_xlabel("Theoretical N(0,1) quantiles", fontsize=18)
axes[1].set_ylabel("Sample quantiles", fontsize=18)
axes[1].set_title("Q-Q plot of z_t = r_t/σ_t", fontsize=18)
axes[1].tick_params(labelsize=16)
axes[1].legend(frameon=False, fontsize=16)

plt.tight_layout(); plt.show()

quietly import delimited "../data/bs-garch.csv", clear
quietly destring _all, replace
gen sigma_t = sqrt(h_true)
gen upper   =  2 * sigma_t
gen lower   = -2 * sigma_t

twoway (line r t,     lcolor(gs10) lwidth(thin)) ///
       (line upper t, lcolor("216 90 48")) ///
       (line lower t, lcolor("216 90 48")), ///
   legend(order(1 "Return r_t" 2 "±2·σ_t")) ///
   ytitle("r_t", size(medlarge)) xtitle("t", size(medlarge)) ///
   title("Returns with ±2σ_t bands", size(medlarge)) ///
   name(g_bands5, replace)

gen z_std5 = r / sigma_t
qnorm z_std5, mcolor("24 95 165") msize(small) ///
    rlopts(lcolor("216 90 48")) ///
    title("Q-Q plot of z_t = r_t/σ_t", size(medlarge)) ///
    name(g_qq5, replace)

graph combine g_bands5 g_qq5, cols(2)

GARCH Bands & FHS Interpretation

Left: returns hug the ±2σ_t bands → GARCH captures volatility clustering.

Right: standardised \(z_t = r_t/\hat\sigma_t\) lie close to the 45° line → the GARCH filter whitens the series.

FHS resamples these \(\hat{z}_t\), preserving the volatility dynamics. Full analysis in Application TS3.

TS 3 — GARCH Filtered Bootstrap: Theory

Problem: parametric bootstrap for GARCH needs to account for time-varying volatility.

GARCH(1,1) model: \[r_t = \sigma_t z_t, \quad \sigma_t^2 = \omega + \alpha r_{t-1}^2 + \beta\sigma_{t-1}^2\]

Filtered Historical Simulation (FHS):

1. Fit GARCH → \(\hat\omega, \hat\alpha, \hat\beta, \hat\sigma_t\)
2. Standardised residuals: \(\hat{z}_t = r_t / \hat\sigma_t\)
3. Bootstrap \(\hat{z}_t\) with replacement → \(z_t^*\)
4. Reconstruct: \(r_t^* = \hat\sigma_t z_t^*\) (using fitted \(\hat\sigma_t\))
5. Re-estimate GARCH on \(\{r_t^*\}\) → \((\hat\omega^*, \hat\alpha^*, \hat\beta^*)\)
6. Repeat \(B\) times → CI on GARCH parameters

Advantage over pure parametric bootstrap: \(z_t^*\) are drawn from the empirical distribution of standardised residuals, not forced to be Gaussian. This captures excess kurtosis and asymmetry in practice.

Stationarity Constraint

\[\alpha + \beta < 1\] required for covariance stationarity. Our DGP: \(0.15 + 0.80 = 0.95 < 1\).

Package Notes

Use rugarch (R) or arch (Python) for GARCH estimation — both handle the bootstrap straightforwardly. Drop Stata for this application: Stata’s arch command works but is slower and lacks easy bootstrap integration.

TS 3 — GARCH Bootstrap: Code

• R
• Python

library(rugarch)

spec <- ugarchspec(
  variance.model   = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model       = list(armaOrder = c(0, 0), include.mean = FALSE),
  distribution.model = "norm"
)

garch_fit <- ugarchfit(spec, data = bs_garch$r)
show(garch_fit)

# Filtered Historical Simulation bootstrap
set.seed(14159)
boot_garch <- ugarchboot(garch_fit, method = "Partial",
                         sampling  = "raw",
                         n.ahead   = 1,
                         n.bootpred= B,
                         n.bootfit = B)
show(boot_garch)

GARCH(1,1) estimates:
  omega = 0.0934  (true 0.05)
  alpha = 0.2748  (true 0.15)
  beta  = 0.6695  (true 0.80)
  alpha+beta = 0.9443  (< 1 = stationary)

import numpy as np, pandas as pd
from arch import arch_model
import os, warnings
os.environ.setdefault("PYTHONWARNINGS", "ignore::UserWarning")
warnings.filterwarnings("ignore")
from joblib import Parallel, delayed
import warnings; warnings.filterwarnings("ignore")

df_g = pd.read_csv("../data/bs-garch.csv")
r_   = df_g["r"].to_numpy()
B_   = 500; NC_ = 12; SEED = 14159

# Fit GARCH(1,1)
am      = arch_model(r_, vol="Garch", p=1, q=1, dist="Normal")
res     = am.fit(disp="off")
print(f"omega = {res.params['omega']:.4f}  (true 0.05)")
print(f"alpha = {res.params['alpha[1]']:.4f}  (true 0.15)")
print(f"beta  = {res.params['beta[1]']:.4f}  (true 0.80)")

# Filtered Historical Simulation bootstrap
z_std   = res.resid / np.sqrt(res.conditional_volatility**2)
z_std   = z_std[~np.isnan(z_std)]
sig_fit = np.sqrt(res.conditional_volatility**2)
sig_fit = sig_fit[~np.isnan(sig_fit)]

def garch_boot(b, seed_base=SEED):
    rng    = np.random.default_rng(seed_base + b)
    z_star = rng.choice(z_std, len(z_std), replace=True)
    r_star = sig_fit * z_star
    try:
        fit_b  = arch_model(r_star, vol="Garch", p=1, q=1).fit(disp="off")
        return [fit_b.params["omega"],
                fit_b.params["alpha[1]"],
                fit_b.params["beta[1]"]]
    except Exception:
        return [np.nan, np.nan, np.nan]

boot_params = np.array(Parallel(n_jobs=NC_)(delayed(garch_boot)(b) for b in range(B_)))
boot_params = boot_params[~np.isnan(boot_params).any(axis=1)]

print(f"\nBootstrap 95% CI (B={len(boot_params)} valid):")
for i, name in enumerate(["omega","alpha","beta"]):
    lo, hi = np.percentile(boot_params[:,i], [2.5, 97.5])
    print(f"  {name}: [{lo:.4f}, {hi:.4f}]")

Advantages, Limitations & Future

τὸ χρηστὸν εἶναι, μέτρια δʼ ἐξαρκεῖν ἔφη.

he said that to be useful, and to have enough, was sufficient

Εὐριπίδης, Ἱκέτιδες 866

Advantages and Limitations

• Advantages
• Disadvantages
• Settings

•  No distributional assumptions on errors
•  Valid under heteroskedasticity, heavy tails, skewness
•  Higher-order accuracy: \(O(n^{-2})\) vs \(O(n^{-1})\) asymptotics
•  Works for complex nonlinear statistics (IRFs, quantiles, ratios)
•  Naturally provides CIs without closed-form variance formulas
•  Embarrassingly parallel — scales with available cores
•  Wild cluster bootstrap solves the few-clusters problem

•  Computationally intensive (\(B\) model fittings)
•  Not valid for non-pivotal statistics (e.g., max of \(t\)-statistics)
•  Boundary cases: testing \(\theta = 0\) when \(\theta \ge 0\) constrained
•  Near-nonstationary series (\(\rho \to 1\)): bootstrap inconsistent unless unit root is imposed
•  Exact unit root: pairs/residual bootstrap invalid; sieve bootstrap needed
•  Very small \(n\): EDF is a poor proxy for true \(F\)

Setting	Recommended bootstrap
Heterosked. OLS	Pairs or Wild
Few clusters (\(G<30\))	Wild cluster
Weak IV	Pairs + AR test
Stationary TS	Block or Sieve
Near unit root	Sieve (AR residuals)
Exact unit root	Impose \(H_0\), sieve
GARCH	FHS (filtered)
Quantile regression	Pairs or wild
High-dimensional (LASSO)	Bootstrap with debiasing

Studentized Interval

The studentized (percentile-\(t\)) interval dominates other CI types in coverage accuracy when a pivotal statistic exists. Use boot.ci(type="stud") and compute a per-replication SE.

Cutting-Edge Research

• Wild Cluster Bootstrap
• Bootstrap + ML
• Time Series
• Future Directions

•  MacKinnon, Nielsen & Webb (2023) — “Fast and Reliable Jackknife and Bootstrap Methods for Cluster-Robust Inference.” Journal of Applied Econometrics.
  •  Fast analytical wild cluster bootstrap that avoids the \(B\) replications entirely in special cases
  •  Valid with as few as \(G = 2\) clusters under certain conditions
•  Djogbenou, MacKinnon & Nielsen (2019) — Asymptotic theory for wild cluster bootstrap under many clusters. JoE 212(2), 393–412.
•  Canay, Santos & Shaikh (2021) — The wild bootstrap with a “small” number of “large” clusters. REStat 103(2), 346–363.

•  Chernozhukov, Chetverikov & Kato (2017) — Central limit theorems and bootstrap in high dimensions; the theory behind post-LASSO bootstrap inference, which requires debiasing before bootstrapping.
•  Wager & Athey (2018) — Bootstrap-based inference for causal forests and generalised random forests. JASA 113(523), 1228–1242.
•  Efron (2020) — “Prediction, estimation, and attribution.” JASA 115(530), 636–655. Bootstrap prediction error in modern ML settings.
•  Romano & Wolf (2005) — Stepdown multiple testing with bootstrap; controls FWER without Bonferroni conservatism. JASA 100(469), 94–108.

•  Kreiss & Paparoditis (2015) — Bootstrapping locally stationary processes. JRSS-B 77(1), 267–290.
•  Gonçalves & Kilian (2004) — Bootstrap for VAR under conditional heteroskedasticity; standard residual resampling is invalid — use wild bootstrap with heteroskedasticity-consistent blocks.
•  Palm, Smeekes & Urbain (2011) — Bootstrap unit root tests in panel time series with cross-sectional dependence. JoE 163(1), 85–104.

•  Bootstrap for synthetic control methods (permutation-based \(p\)-values already standard; bootstrap variance estimation developing)
•  Bootstrap + deep learning: efficiency gains via gradient-based bootstrap (implicitly differentiating through the estimator)
•  Bootstrap for moment-inequality models (partial identification; active research area)
•  Distributed/streaming bootstrap: bootstrapping datasets too large for RAM via sketching and subsampling
•  Bootstrap calibration via simulation: using ML to learn the mapping from statistic to accurate \(p\)-value

Further Reading

Textbooks

•  Efron & Tibshirani (1993) — An Introduction to the Bootstrap. Chapman & Hall. The original textbook — readable and essential.
•  Davison & Hinkley (1997) — Bootstrap Methods and Their Application. Cambridge University Press. More theoretical; comprehensive.
•  Davidson & MacKinnon (2004) — Econometric Theory and Methods. Oxford University Press. Chapter 4: resampling and bootstrap in an econometric context.
•  Cameron & Trivedi (2005) — Microeconometrics: Methods and Applications. Cambridge. Chapter 11: bootstrap in applied micro.
•  Horowitz (2001) — “The Bootstrap.” Handbook of Econometrics, Vol. 5, 3159–3228. Comprehensive survey with econometrics focus.

Key articles — methods

•  Efron (1979). “Bootstrap methods: another look at the jackknife.” Annals of Statistics 7(1), 1–26.
•  Bickel & Freedman (1981). Annals of Statistics 9(6), 1196–1217.  ·  Freedman (1981). Annals of Statistics 9(6), 1218–1228.
•  Wu (1986). Annals of Statistics 14(4), 1261–1295.
•  Mammen (1993). Annals of Statistics 21(1), 255–285.
•  Künsch (1989). Annals of Statistics 17(3), 1217–1241.  ·  Bühlmann (1997). Bernoulli 3(2), 123–148.
•  Politis & Romano (1994). JASA 89(428), 1303–1313.
•  Cameron, Gelbach & Miller (2008). Review of Economics and Statistics 90(3), 414–427.
•  MacKinnon & Webb (2017). “Wild bootstrap inference for wildly different cluster sizes.” Journal of Applied Econometrics 32(2), 233–254.

Key articles — inference & refinement

•  Beran (1988). “Prepivoting test statistics.” JASA 83(403), 687–697.
•  Davidson & MacKinnon (1999). “The size distortion of bootstrap tests.” Econometric Theory 15(3), 361–376.
•  Andrews & Buchinsky (2000). “A three-step method for choosing \(B\).” Econometrica 68(1), 23–51.
•  Dufour & Khalaf (2003). “Monte Carlo test methods.” A Companion to Theoretical Econometrics, 494–519.
•  Davidson & Flachaire (2008). “The wild bootstrap, tamed at last.” Journal of Econometrics 146(1), 162–169.
•  Hall, Horowitz & Jing (1995). “On blocking rules for the bootstrap with dependent data.” Biometrika 82(3), 561–574.
•  Politis & White (2004) and Patton, Politis & White (2009). Automatic block-length selection. Econometric Reviews 23(1) / 28(4).

Surveys used in these slides

•  MacKinnon (2006) — “Bootstrap Methods in Econometrics.” QED Working Paper 1028. Source for the null-imposition, Monte Carlo-exactness, and DGP-comparison material.
•  Cameron (2022) — “Bootstrap Methods” (U.C.-Davis lecture notes). Source for the jackknife and percentile-\(t\) Stata workflow.
•  MacKinnon, Nielsen & Webb (2023) — “Cluster-robust inference: a guide to empirical practice.” Journal of Econometrics 232(2), 272–299.

Software documentation

•  boot (R): Canty & Ripley. doi:10.32614/CRAN.package.boot
•  fwildclusterboot (R): implements Roodman et al. (2019). s3alfisc.github.io/fwildclusterboot
•  arch (Python): Kevin Sheppard. arch.readthedocs.io
•  boottest (Stata): Roodman, Nielsen, MacKinnon & Webb (2019), Stata Journal 19(1), 4–60. SSC.

Online resources

•  James MacKinnon’s working papers on cluster bootstrap: ideas.repec.org/f/pma144
•  Nick Huntington-Klein — bootstrap intro: theeffectbook.net
•  Peter Bühlmann — sieve bootstrap: stat.ethz.ch/~buhlmann

Journals

Journal of Econometrics · Econometrica · Journal of Applied Econometrics · Econometric Theory · Review of Economics and Statistics

Exercises

1. Coverage simulation. Generate 1 000 datasets from bs-cross DGP with \(N=50\). For each, compute: OLS CI, HC3 CI, pairs bootstrap CI, and BCa CI. Report empirical coverage. Which method has the best size under \(n=50\)?

2. Wild vs pairs bootstrap. In the bs-cross DGP, multiply \(\sigma_i\) by 3 for units with \(x_1 > 0\) (extreme heteroskedasticity). Compare the SE and coverage of the pairs bootstrap, wild bootstrap (Rademacher), and HC3 SE via simulation.

3. Cluster bootstrap power. Using bs-cluster with \(G=15\), set the true treatment effect to zero (\(\beta_{tr}=0\)). Simulate 500 datasets and compute rejection rates at the 5% level for: (a) OLS \(t\)-test, (b) asymptotic cluster \(t\), (c) wild cluster bootstrap. Which test is best-calibrated?

4. IV bootstrap under strong instruments. Re-simulate bs-iv with \(\pi = 1.0\) (strong instruments, \(F \approx 500\)). Compare the bootstrap CI and the asymptotic CI. Do they agree? Now increase endogeneity to \(\text{Corr}(\varepsilon,v)=0.95\). Which deteriorates first?

5. Block length sensitivity. For the AR(1) bootstrap (App TS 1), run MBB with \(\ell \in \{3, 7, 14, 28\}\). Plot bootstrap SE(\(\hat\rho\)) as a function of \(\ell\). Is there an optimal block length that minimises MSE?

6. Bootstrap vs delta method for IRFs. Using bs-ts, compute the 95% CI for the VAR(1) IRF at horizon \(h=5\) via (a) delta method, (b) bootstrap with \(B=500\). Do the intervals differ? At which horizon is the difference largest?

7. GARCH bootstrap in practice. Using bs-garch, fit a GARCH(1,1) and compute bootstrap CIs for \(\hat\omega, \hat\alpha, \hat\beta\). Now impose \(\alpha + \beta < 1\) strictly during bootstrap re-estimation (projection step). Does the CI change substantially?

Thank You

Athanassios Stavrakoudis
Applied Informatics and Computational Economics Lab
Department of Economics
University of Ioannina, Greece

astavrak@uoi.gr · linkedin.com/in/astavrakoudis