Wild cluster bootstrap in Python: a boottest workflow for few clusters

Cluster-robust standard errors are justified by letting the number of clusters grow. Many studies have few of them: a policy adopted by some of a dozen states, an experiment run in eight schools, a panel of twenty countries. With few clusters, or clusters of very different sizes, the usual cluster-robust t test can reject a true null hypothesis far more often than its nominal level. The wild cluster bootstrap is the standard way to obtain more reliable p-values and confidence intervals in that situation.

What the wild cluster bootstrap does

The method builds many artificial samples that keep the dependence within clusters and the regressors exactly as observed. Cameron, Gelbach and Miller (2008) showed that it controls the rejection rate much better than conventional cluster-robust tests when there are few clusters. The restricted version imposes the null hypothesis when it generates the samples. Each artificial outcome is

y*_ig = x_ig′β̃ + ũ_ig · v*_g

where β̃ and ũ are the coefficients and residuals estimated with the null imposed, and v*_g is one random weight for cluster g, so that all residuals in a cluster are flipped or rescaled together. The procedure has four steps:

1. Estimate the model with the null hypothesis imposed, for example with the policy coefficient set to zero.
2. For each draw, build a new outcome y* with fresh cluster weights.
3. Refit the unrestricted model on y*, recompute the cluster-robust (CR1) standard error, and compute the bootstrap t or Wald statistic for the same hypothesis.
4. The p-value is the share of bootstrap statistics at least as extreme as the one from the original data.

The weights have mean zero and variance one. Rademacher weights are +1 or −1 with equal probability, so with G clusters there are only 2^G distinct sign patterns, which limits how small a p-value can be. Webb's six-point distribution allows 6^G patterns and was proposed for very few clusters (Webb 2023). Mammen's two-point distribution is a third option.

A confidence set comes from test inversion: it contains every value of the coefficient that the bootstrap test does not reject. Roodman, Nielsen, MacKinnon and Webb (2019) describe this workflow and its fast implementation in boottest, a community-contributed Stata command. MacKinnon, Nielsen and Webb (2023) review when the method still struggles: with very few treated clusters, for instance, no cluster-robust procedure is fully reliable.

In OpenEconometrics

Both functions take a fitted oe.ols result, or an absorbed or panel fixed-effects fit, together with the data it was fitted on: oe.wild_cluster_test(fit, data=df, null={"policy": 0}, wild="webb", reps=1999, seed=1) tests a hypothesis, and oe.wild_cluster_confidence_set(fit, data=df, term="policy", grid=grid, ...) inverts the test.

• null maps terms to their hypothesized values. One entry gives a test of one coefficient; several give a joint Wald test.
• cluster names the bootstrap cluster column. It defaults to the fit's own cluster column, so the fit itself may use any covariance.
• wild="rademacher" (the default), "webb" or "mammen", with reps draws and a seed. enumerate_all=True evaluates every Rademacher sign pattern instead of sampling, which is possible with up to 12 clusters.
• grid is an increasing list of candidate values for one term. The same draws are used at every grid point.

Sampled p-values use the plus-one convention, (1 + exceedances)/(B + 1). The test returns tables: out["test"] holds the statistic, p-value, degrees of freedom and the Monte Carlo standard error of the p-value; out["settings"] records the seed, weights, cluster sizes and conventions. The confidence set returns the p-value at each grid point in out["grid"] and the accepted components in out["components"]. For OLS clustered by state, the test corresponds to boottest policy, reps(1999) weighttype(webb) in Stata after regress y policy x, vce(cluster state). IV models and two-way clustering have their own functions, oe.iv_wild_test and oe.wild_multiway_test, described in the IV and multiway wild bootstrap guide.

Supported options

oe.wild_cluster_test

Stata command: boottest.

oe.wild_cluster_confidence_set

Stata command: boottest.

The inventory section inference_extensions.restricted_wild_cluster for oe.wild_cluster_test and oe.wild_cluster_confidence_set lists multipliers: rademacher, mammen, webb.

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 525 workers in 12 states whose sizes range from 10 to 90 workers. Five states adopt a policy with a true effect of 0.20 on the outcome. Each state adds a shared shock with standard deviation 0.3, so the errors are correlated within states, and a worker-level covariate x has coefficient 0.5. The question is what the data say about the policy effect, with only 12 clusters and 5 of them treated.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n_states=12, n_treated=5, seed=12):
    """Workers in 12 states of unequal size; a policy adopted by 5 states; a shared state shock."""
    rng = np.random.default_rng(seed)
    sizes = rng.integers(10, 91, size=n_states)
    state = np.repeat(np.arange(n_states), sizes)
    policy = np.repeat((np.arange(n_states) < n_treated).astype(int), sizes)
    state_shock = np.repeat(rng.normal(scale=0.3, size=n_states), sizes)
    x = rng.normal(size=state.size)
    y = 1.0 + 0.2 * policy + 0.5 * x + state_shock + rng.normal(size=state.size)
    return pd.DataFrame({"y": y, "policy": policy, "x": x, "state": state})


if __name__ == "__main__":
    df = simulate()
    fit = oe.ols(data=df, y="y", x=["policy", "x"], cluster="state")
    c = {k.term: k for k in fit.coefficients}["policy"]
    print(f"N = {len(df)}, clusters = {df['state'].nunique()}")
    print(f"policy: {c.estimate:.4f}, CR1 se {c.std_error:.4f}, t(11) p = {c.p_value:.4f}, "
          f"95% CI [{c.ci_low:.4f}, {c.ci_high:.4f}]")
    webb = oe.wild_cluster_test(fit, data=df, null={"policy": 0}, cluster="state",
                                wild="webb", reps=1999, seed=1)
    print(f"Wild cluster bootstrap, Webb, 1,999 draws: p = {webb['test']['p_value'][0]:.4f}")
    signs = oe.wild_cluster_test(fit, data=df, null={"policy": 0}, cluster="state",
                                 wild="rademacher", enumerate_all=True)
    print(f"Wild cluster bootstrap, all 4,096 Rademacher patterns: "
          f"p = {signs['test']['p_value'][0]:.4f}")
    grid = np.linspace(-1.0, 1.5, 51).round(2).tolist()
    cs = oe.wild_cluster_confidence_set(fit, data=df, term="policy", grid=grid,
                                        cluster="state", wild="webb", reps=1999, seed=1,
                                        max_work=10**9)
    for row in cs["components"].itertuples():
        print(f"Wild 95% confidence set: accepted {row.lower_grid:.2f} to {row.upper_grid:.2f}, "
              f"rejected at {row.left_rejected_grid:.2f} and {row.right_rejected_grid:.2f}")

With openecon 0.3.20a1 from PyPI the script prints:

N = 525, clusters = 12
policy: 0.0570, CR1 se 0.1944, t(11) p = 0.7748, 95% CI [-0.3709, 0.4849]
Wild cluster bootstrap, Webb, 1,999 draws: p = 0.8025
Wild cluster bootstrap, all 4,096 Rademacher patterns: p = 0.7915
Wild 95% confidence set: accepted -0.40 to 0.50, rejected at -0.45 and 0.55

Reading the output

• The estimate is imprecise, as it should be. The policy coefficient is 0.0570 against a true 0.20. With 12 states and a state shock that every worker in a state shares, the data contain little information about a state-level policy, and the CR1 standard error of 0.1944 says so.
• Here the bootstrap and the t test agree. The CR1 t test with 11 degrees of freedom gives p = 0.7748, the Webb bootstrap p = 0.8025 and the full Rademacher enumeration p = 0.7915. With 12 clusters, 5 of them treated, the three procedures give similar answers in this sample; the conventional test already uses a t distribution with G − 1 degrees of freedom. The guide to robust and clustered standard errors shows a case with four clusters where the procedures disagree by orders of magnitude.
• The confidence set is slightly wider. The bootstrap accepts every grid value from −0.40 to 0.50 and rejects −0.45 and 0.55, so its limits lie within those brackets. The CR1 interval is −0.3709 to 0.4849. Both contain the true 0.20, and both contain zero.
• Exact enumeration removes simulation noise. With 12 clusters there are 2^12 = 4,096 Rademacher sign patterns, so enumerate_all=True evaluates all of them and the p-value has no Monte Carlo error. The Webb p-value reports a Monte Carlo standard error in out["test"].
• The work budget is explicit. A confidence set repeats the bootstrap at every grid point; with 51 points and 1,999 draws the estimated work exceeds the default max_work, so the script raises it. It runs in a few seconds.

Limitations

• Supported domain of oe.wild_cluster_test and oe.wild_cluster_confidence_set: numeric unweighted OLS/absorbed FE/panel FE; independent one-way clusters; FE nested within clusters; full rank.
• Size limits of oe.wild_cluster_test and oe.wild_cluster_confidence_set: N<=8192; K_full<=256; Rademacher complete enumeration G<=12; bounded work and workspace.
• The bootstrap restricts linear coefficients of a model without weights or categorical regressors; IV models and two-way clustering need oe.iv_wild_test and oe.wild_multiway_test. In a fixed-effects model only reported slopes can be restricted.
• Sampled bootstraps use 49 to 100,000 draws and at most 2 million multiplier cells; max_work caps the estimated work before anything runs, and a confidence set multiplies it by the number of grid points.
• The confidence set is evaluated only on the declared grid of 2 to 201 increasing values. It reports accepted grid points and bracketing rejected points, without interpolation, and says nothing about values outside the grid.
• Complete Rademacher enumeration is conditional on the fitted sample. It removes simulation noise but is not an exact randomization test.
• A draw with a singular cluster-robust covariance stops the procedure with an error; no draws are dropped silently.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Cameron, A. C., Gelbach, J. B. and Miller, D. L. (2008). Bootstrap-based improvements for inference with clustered errors. Review of Economics and Statistics, 90(3), 414–427. doi:10.1162/rest.90.3.414
• MacKinnon, J. G., Nielsen, M. Ø. and Webb, M. D. (2023). Cluster-robust inference: A guide to empirical practice. Journal of Econometrics, 232(2), 272–299. doi:10.1016/j.jeconom.2022.04.001
• Roodman, D., Nielsen, M. Ø., MacKinnon, J. G. and Webb, M. D. (2019). Fast and wild: Bootstrap inference in Stata using boottest. Stata Journal, 19(1), 4–60. doi:10.1177/1536867X19830877
• Webb, M. D. (2023). Reworking wild bootstrap-based inference for clustered errors. Canadian Journal of Economics, 56(3), 839–858. doi:10.1111/caje.12661