import matplotlib
matplotlib.use("Agg") # headless, non-interactive backend
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from linearmodels.iv import IV2SLS # linearmodels 7.0 import path
rng = np.random.default_rng(20260528) # one seeded generator for the whole notebook
import linearmodels
print("numpy", np.__version__, "| pandas", pd.__version__, "| linearmodels", linearmodels.__version__)nb4.5 — Bartik / Shift-Share Instrument Decomposition
Companion to Ch 4.5, “Bartik / Shift-Share Designs.” This notebook closes the causal-inference arc of Weeks 3–4.
The chapter made a bold claim: when the world refuses to hand you a clean coin flip, you can build an instrument out of pieces you already have — national industry shifts g_j (growth rates no single town caused) and local baseline shares s_{rj} (how much each town rode on each industry) — combined as
B_r \;=\; \sum_j s_{rj}\, g_j.
This is the Bartik / shift-share instrument. The chapter also promised two things that sound almost contradictory and one thing that sounds like magic. The magic: with an unobserved regional confounder poisoning the local regressor, OLS lies, but instrumenting with B_r recovers the truth. The two views: Goldsmith-Pinkham–Sorkin–Swift (GPSS, 2020) say identification comes from the shares; Borusyak–Hull–Jaravel (BHJ, 2022) say it comes from the shifts. This notebook turns all of that into running code.
Maya is our anchor, as in the chapter. She wants to know whether a negative shock to local labor demand — a downturn in the industries that happen to employ a region’s workers — pushes households in that region into credit-card delinquency and mortgage default. The realized local shock is endogenous (local optimism, housing bubbles, labor-supply moves all hide in the error), so OLS is hopeless. The Bartik instrument is her escape.
We do six things, in order.
- Build a region × industry world with shares, shifts, an endogenous local shock, and a confounder that lives in the outcome’s error — engineered so OLS is guaranteed to lie.
- Show OLS is biased (\approx -1.06 against a true -1.5).
- Instrument with B_r in 2SLS (
linearmodels.iv.IV2SLS), recover the true -1.5, and read a first-stage F in the hundreds. - GPSS decomposition: compute the exact Rotemberg weights \hat\alpha_j, confirm \hat\beta^{\text{Bartik}} = \sum_j \hat\alpha_j \hat\beta_{1,j} reassembles the headline number to many decimals, and see which industries drive it.
- The stress demo: tie a few industry shares to the confounder and watch the most-suspect shares earn the biggest weights while the headline drifts from the truth — even though the first-stage F still looks reassuring.
- BHJ view + AKM inference: an effective-number-of-shocks diagnostic, and a shock-level look at why ordinary region-clustered standard errors are too small.
The chapter’s headline numbers — OLS \approx -1.06, 2SLS \approx -1.5, first-stage F \approx 571, Rotemberg weights reassembling to -1.5 — are reproduced exactly below.
Setup
We fix one seeded generator so the notebook reproduces bit-for-bit (the CONVENTIONS.md rule), and force matplotlib’s non-interactive Agg backend so it runs headless — laptop, server, or CI — with no display attached.
The seed is 20260528, the same one the chapter’s code block used, so our numbers match the prose exactly. The one import to notice is the last: in linearmodels 7.0 the 2SLS estimator lives at linearmodels.iv.IV2SLS, and that is the path we use throughout.
1. A region × industry world built to break OLS
We simulate R = 400 regions (think commuting zones) and J = 20 industries. The data-generating process mirrors the chapter exactly, with every channel labeled. The true causal effect is \beta_1 = -1.5: a worse local labor-demand shock raises the household-delinquency change (coded so the sign works out).
The four ingredients:
- Confounder
confounder\sim \mathcal N(0,1) — latent regional fragility (loose lending standards, weak balance sheets). It is unobserved; we simulate it but never let any estimator see it, and it lives in the outcome’s error. - Baseline local shares
S— a Dirichlet-style draw, one row per region, each row summing to 1. We userng.gamma(shape=0.5, ...)then normalize, giving a sparse-ish industry mix (each region leans on a few industries). In this clean baseline the shares are drawn independently of the confounder — the GPSS exogenous-shares ideal. §5 deliberately breaks this. - National shifts
g\sim \mathcal N(0, 0.05) — the industry growth rates, as-good-as-random across the 20 industries (the BHJ ideal). A national number: the same auto shock hits every region, scaled by exposure. - The Bartik instrument
B = S @ g, standardized for clean scaling.
Then the realized local shock blends the clean Bartik part with the endogenous confounder, and the outcome gets the true effect plus the confounder in its error.
The spec in CONVENTIONS §4 discipline: outcome = regional delinquency-rate change; key regressor = realized local employment growth (shock); controls = constant only in this pass; instrument = the Bartik shift-share B_r; identifying assumption = B_r affects delinquency only through the realized local shock and is uncorrelated with regional fragility (defended via exogenous shares in §4–5, or exogenous shifts in §6).
R, J = 400, 20 # 400 regions, 20 industries
beta1_true = -1.5 # the TRUE causal effect we will try to recover
# Unobserved confounder: latent regional fragility (lives in the outcome's error)
confounder = rng.normal(size=R)
# Baseline LOCAL shares s_rj: a Dirichlet-style draw per region, rows sum to 1.
# Clean baseline: shares are drawn INDEPENDENTLY of the confounder (we STRESS this in section 5).
alpha0 = rng.gamma(shape=0.5, size=(R, J)) # sparse-ish industry mix per region
S = alpha0 / alpha0.sum(axis=1, keepdims=True) # (R x J), each row sums to 1
# NATIONAL shifts g_j: as-good-as-random across the 20 industries (the BHJ ideal)
g = rng.normal(0.0, 0.05, size=J)
# Shift-share (Bartik) instrument: B_r = sum_j s_rj * g_j (standardized for clean scaling)
B = S @ g
B = (B - B.mean()) / B.std() # (R,)
# Realized local shock = clean part (driven by B) + endogenous local part (the confounder)
local_shock = 1.0 * B + 0.8 * confounder + rng.normal(0, 0.3, size=R)
# Outcome: delinquency change. True effect beta1_true; the confounder ALSO raises delinquency.
eps = 1.0 * confounder + rng.normal(0, 0.3, size=R) # confounder lives in the error
Y = 0.0 + beta1_true * local_shock + eps
df = pd.DataFrame({"Y": Y, "shock": local_shock, "B": B, "const": 1.0})
print(f"R = {R} regions, J = {J} industries, true beta1 = {beta1_true}")
print(f"each region's shares sum to 1? {np.allclose(S.sum(axis=1), 1.0)}")
print(f"corr(B, confounder) = {np.corrcoef(B, confounder)[0, 1]:+.3f} <- clean: instrument is exogenous")
print(f"corr(shock, confounder) = {np.corrcoef(local_shock, confounder)[0, 1]:+.3f} <- the endogeneity in the regressor")Looking at the world before we estimate
Two pictures make the construction concrete. On the left, the shares: a heatmap of the first few regions’ exposure to each industry — each row sums to 1, and the sparse-ish draw means each region leans heavily on a handful of industries (the dark cells). On the right, the instrument B_r against the confounder: the cloud is a shapeless blob (correlation near zero), which is exactly what makes B_r a candidate instrument — it carries the national-tide signal but not the local fragility.
This is the chapter’s §2 arithmetic in pictures: B_r turns one shared set of national shifts into a different predicted shock per region, and the only thing distinguishing regions is their shares.
fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))
# Left: shares heatmap for the first 12 regions
ax = axes[0]
im = ax.imshow(S[:12], aspect="auto", cmap="viridis")
ax.set_xlabel("industry j"); ax.set_ylabel("region r (first 12)")
ax.set_title("A. Baseline shares s_rj (each row sums to 1)")
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="share")
# Right: instrument vs confounder -- a shapeless blob (clean exogeneity)
ax = axes[1]
ax.scatter(confounder, B, s=14, alpha=0.5, color="#1f77b4")
ax.set_xlabel("unobserved regional fragility (confounder)")
ax.set_ylabel("Bartik instrument B_r (standardized)")
ax.set_title(f"B. Instrument is exogenous: corr = {np.corrcoef(B, confounder)[0,1]:+.3f}")
fig.tight_layout()
fig.savefig("nb4.5-world.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("saved figure: nb4.5-world.png")
print("Left: each region leans on a few industries (sparse shares). Right: B carries no fragility -> exogenous.")2. OLS is biased — and we can see whyRun the regression Maya is tempted to report: the delinquency change Y on a constant and the realized local shock. With no endogenous regressor declared, IV2SLS simply runs OLS, giving a clean apples-to-apples comparison with the 2SLS call to come.The bias has a predictable sign. Fragile regions (high confounder) both have a higher realized shock (the +0.8 * confounder channel) and higher delinquency on their own (the +1.0 * confounder in Y’s error). Both channels enter with a positive sign, so together they manufacture a spurious positive association between the shock and delinquency — which works against the true negative effect and pulls \hat\beta_1 toward zero. We therefore expect \approx -1.06, a number whose magnitude understates the true -1.5. The point is what matters: OLS is wrong, and we know it is wrong — and in which direction — by construction.
# Naive OLS via IV2SLS: dependent, exog=[const, shock], no endog, no instruments.
ols = IV2SLS(df["Y"], df[["const", "shock"]], None, None).fit()
beta_ols = ols.params["shock"]
print(f"OLS beta_hat(shock) = {beta_ols:+.3f} (true = {beta1_true})")
print(f"OLS bias = {beta_ols - beta1_true:+.3f} "
f"(toward zero: the positive confounder channel masks the true negative effect)")3. 2SLS with the Bartik instrument recovers the truth
Now the instrument earns its keep. We declare the local shock endogenous and hand IV2SLS the Bartik instrument B. The call signature is
IV2SLS(dependent, exog=[const + controls], endog=[shock], instruments=[B])We fit with cov_type="robust" (heteroskedasticity-robust SEs — the Ch 2.4 default). The point estimate should land right on the true -1.5, and the first-stage F should be enormous because we built a real first stage (1.0 * B feeds directly into shock).
linearmodels hands you the first-stage diagnostics in iv.first_stage.diagnostics; we read off the partial f.stat. With a single instrument it equals the squared t-statistic on B in the first stage. The chapter’s value is F \approx 571.
# 2SLS: instrument the endogenous local shock with the Bartik instrument B.
iv = IV2SLS(df["Y"], df[["const"]], df[["shock"]], df[["B"]]).fit(cov_type="robust")
beta_2sls = iv.params["shock"]
diag = iv.first_stage.diagnostics
F_stage = float(diag["f.stat"].iloc[0])
print(f"OLS beta_hat(shock) = {beta_ols:+.3f} (biased toward zero)")
print(f"2SLS beta_hat(shock) = {beta_2sls:+.3f} (true = {beta1_true}) SE = {iv.std_errors['shock']:.3f}")
print()
print("First-stage diagnostics (linearmodels):")
print(diag[["rsquared", "f.stat", "f.pval"]].to_string())
print(f"\nfirst-stage partial F = {F_stage:.1f} (>> 10: a STRONG lever)")
print()
print("The OLS-vs-2SLS gap IS the regional endogeneity, made visible:")
print(f" gap = {beta_ols - beta_2sls:+.3f}")6. The BHJ view and AKM inference: shifts, effective shocks, and honest standard errors
GPSS put the exogeneity burden on the shares. Borusyak–Hull–Jaravel flip the lens: treat the national shifts g_j as the random object, as-good-as-randomly drawn across many industries, with shares demoted to mere exposure weights. The shift-share IV moment can be rewritten as an exposure-weighted regression at the level of the shocks (industries). The catch BHJ stress: “many shocks” must be real, not an illusion. If a few industries carry all the weight, the law-of-large-numbers logic fails and you cannot lean on the BHJ story — you are back to defending shares.
The natural diagnostic is the effective number of shocks: the inverse Herfindahl of the (normalized) Rotemberg weights,
n_{\text{eff}} = \frac{1}{\sum_j p_j^2}, \qquad p_j = \frac{|\hat\alpha_j|}{\sum_k |\hat\alpha_k|}.
If n_{\text{eff}} \approx J, exposure is well spread and “many shocks” is honest; if n_{\text{eff}} is small, a handful of shocks run the show. We compute it for the clean baseline, the stressed design, and a deliberately concentrated few-shock world to see the contrast.
def effective_shocks(weights):
p = np.abs(weights) / np.abs(weights).sum()
return 1.0 / np.sum(p ** 2)
# A deliberately CONCENTRATED few-shock world: J=6, very sparse shares (a la Maya's bank case).
conc = np.random.default_rng(20260528)
Rc, Jc = 400, 6
conf_c = conc.normal(size=Rc)
a0_c = conc.gamma(shape=0.2, size=(Rc, Jc)) # shape=0.2 -> highly concentrated shares
S_c = a0_c / a0_c.sum(axis=1, keepdims=True)
g_c = conc.normal(0.0, 0.05, size=Jc)
B_c = S_c @ g_c; B_c = (B_c - B_c.mean()) / B_c.std()
shock_c = 1.0 * B_c + 0.8 * conf_c + conc.normal(0, 0.3, size=Rc)
D_dm_c = shock_c - shock_c.mean()
S_dm_c = S_c - S_c.mean(axis=0, keepdims=True)
alpha_c = g_c * (S_dm_c.T @ D_dm_c); alpha_c = alpha_c / alpha_c.sum()
print("Effective number of shocks (1 / HHI of |Rotemberg weights|):")
print(f" clean baseline (J={J}, dispersed shares): n_eff = {effective_shocks(alpha):5.2f} of {J}")
print(f" stressed (J={J}, suspect-loaded): n_eff = {effective_shocks(alpha_s):5.2f} of {J}")
print(f" concentrated (J={Jc}, sparse shares): n_eff = {effective_shocks(alpha_c):5.2f} of {Jc}")
print()
print("Many dispersed shocks -> n_eff near J -> BHJ 'quasi-random shifts' is honest.")
print("Few concentrated shocks -> tiny n_eff -> 'many shocks' is an illusion; lean on GPSS shares instead.")Adão–Kolesár–Morales: ordinary clustering lies about precision
AKM (2019) showed that regions with similar shares are exposed to the same shocks, so their regression errors are correlated even if they share no geographic cluster — two auto-heavy commuting zones on opposite coasts move together through the common auto shock. Ordinary region-clustered standard errors miss this entirely and come out too small, so you over-reject.
We make the point with a shock-level resampling experiment, which is exactly the BHJ “shocks are the random object” view. We perturb the national shifts g_j (the thing AKM say drives the cross-regional dependence) and recompute the Bartik 2SLS estimate many times; the spread of those estimates is the honest sampling variability. We compare its standard deviation to the vanilla heteroskedasticity-robust SE that IV2SLS reports. We do this in two regimes: the dispersed baseline (many shocks — the gap should be mild) and a concentrated few-shock world (the gap should be large — vanilla SEs are badly too small).
def bartik_beta(S_mat, gvec, shock_vec, y_vec):
# Closed-form single-instrument Bartik 2SLS coefficient (demeaned).
Bx = S_mat @ gvec
Bx = (Bx - Bx.mean()) / Bx.std()
d = shock_vec - shock_vec.mean()
y = y_vec - y_vec.mean()
bx = Bx - Bx.mean()
return (bx @ y) / (bx @ d)
def akm_vs_vanilla(S_mat, gvec, shock_vec, y_vec, label, n_boot=3000, seed=1):
# vanilla heteroskedasticity-robust SE from IV2SLS
Bx = S_mat @ gvec; Bx = (Bx - Bx.mean()) / Bx.std()
d = pd.DataFrame({"Y": y_vec, "shock": shock_vec, "B": Bx, "const": 1.0})
fit = IV2SLS(d["Y"], d[["const"]], d[["shock"]], d[["B"]]).fit(cov_type="robust")
se_vanilla = float(fit.std_errors["shock"])
# shock-level resampling: perturb the shifts (the source of cross-regional dependence)
br = np.random.default_rng(seed)
sd_g = gvec.std()
boot = np.array([bartik_beta(S_mat, gvec + br.normal(0, sd_g, size=gvec.size),
shock_vec, y_vec) for _ in range(n_boot)])
se_shock = boot.std()
print(f"{label}")
print(f" vanilla robust (region) SE = {se_vanilla:.4f}")
print(f" AKM/shock-level resample SD = {se_shock:.4f}")
print(f" ratio (shock / vanilla) = {se_shock / se_vanilla:.2f} "
f"({'OK: mild' if se_shock/se_vanilla < 1.5 else 'DANGER: vanilla SE far too small'})")
return se_vanilla, se_shock
print("How badly does ordinary clustering understate uncertainty?\n")
akm_vs_vanilla(S, g, local_shock, Y, "MANY dispersed shocks (clean baseline, J=20):")
print()
akm_vs_vanilla(S_c, g_c, shock_c,
beta1_true * shock_c + conc.normal(0, 0.3, size=Rc) + conf_c,
"FEW concentrated shocks (J=6):")
print()
print("AKM's lesson: with concentrated exposure the honest SD dwarfs the vanilla SE, so vanilla")
print("region-clustered t-stats over-reject. Report AKM-robust (or BHJ shock-level) SEs instead.")Your Turn
You have watched the clean design succeed and the confounded-shares design fail. Now drive it yourself.
Correlate a confounder with the shares and watch identification break. In §5 we tied the three biggest-shift industries to the confounder. Instead, tie a single, low-shift industry (say the one with the smallest |g_j|) to the confounder with a strong loading (
np.exp(2.0 * confounder)). (a) Does the headline 2SLS drift as much? (b) Look at that industry’s Rotemberg weight — is it large or small? (c) Explain the lesson: a violated share only damages the estimate to the extent it earns weight, and weight comes from the shift. This is why GPSS say to aim your skepticism at the high-weight shares.Sweep the contamination strength. Loop the stress loading from
0.0to2.5on the suspect industries, recording at each step the headline 2SLS estimate and the first-stage F. Plot both against the loading. Confirm in your own picture that the F stays flat and reassuring while the estimate marches away from -1.5 — the “F can’t see exclusion failure” lesson, made into a curve.Drop-the-suspect robustness check. The standard GPSS remedy: re-estimate the Bartik 2SLS excluding the top-weighted industry from the instrument (rebuild B_r from the other J-1 industries, re-standardize, re-fit). In the clean design the estimate should barely move; in the stressed design it should snap back toward -1.5. What does “the estimate barely moves when I drop the suspect” tell you, and what does “it snaps back” tell you?
Effective shocks vs. the AKM gap. Build three worlds with J \in \{5, 20, 80\} industries (same R, same DGP). For each, compute the effective number of shocks and the AKM-vs-vanilla SE ratio from §6. Plot the ratio against n_{\text{eff}}. Does the AKM problem (ratio \gg 1) shrink as the number of effective shocks grows? Connect your finding to the chapter’s claim that BHJ’s “many shocks” assumption is exactly what makes the LLN — and honest inference — kick in.