One Card instrument, many synthetic moments

Mincer returns, PyFixest, T-Rex GMM, and many-instrument asymptotics

Author

Krabbs

Published

August 11, 2026

Question

Suppose a Mincer wage equation has one credible excluded instrument: the Card (1993) indicator for growing up near a four-year college. Can we manufacture more identifying information by interacting that one instrument with a growing collection of baseline covariates?

This notebook constructs

\[ Z_i,\quad Z_iW_{i1},\quad \ldots,\quad Z_iW_{iK}, \]

where only \(Z_i\) has a nonzero population first-stage coefficient. The interactions are valid under the deliberately favorable simulation assumptions, but they are synthetic moments, not new sources of exogenous variation. Each irrelevant interaction can fit a little endogenous first-stage noise in finite samples. When \(K\) grows in proportion to \(n\), those small overfits accumulate and conventional 2SLS drifts toward OLS.

The exercise uses PyFixest for formula-based IV estimation and Apoorva Lal’s T-Rex GMM library for the same linear moment problem.

Data-generating process

The structural equation is

\[ \log(w_i) = 1.5 + 0.08S_i + 0.035X_i - 0.0006X_i^2 + u_i, \]

and schooling is

\[ S_i = 12 + 0.75Z_i + 0.8A_i + 0.15\widetilde X_i + v_i. \]

Ability \(A_i\) is unobserved and enters \(u_i=0.20A_i+\varepsilon_i\), so schooling is endogenous. Proximity \(Z_i\) is independent of \((A_i,v_i,\varepsilon_i,W_i)\) and affects wages only through schooling. The auxiliary covariates \(W_{ij}\) are independent standard normals. Consequently every \(Z_iW_{ij}\) is a valid instrument, but its population first-stage coefficient is zero.

Show Python 
import platform
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyfixest as pf
import scipy
import torch

from trex.gmm import GMMEstimator

warnings.filterwarnings("ignore", category=FutureWarning)
plt.style.use("seaborn-v0_8-whitegrid")

BETA = 0.08

versions = pd.DataFrame(
    {
        "package": ["Python", "NumPy", "SciPy", "PyFixest", "PyTorch"],
        "version": [
            platform.python_version(),
            np.__version__,
            scipy.__version__,
            pf.__version__,
            torch.__version__,
        ],
    }
)
versions
package version
0 Python 3.12.9
1 NumPy 2.5.2
2 SciPy 1.18.0
3 PyFixest 0.60.0
4 PyTorch 2.13.0
Show Python 
def make_mincer(n=2_000, k_fake=20, seed=93_0811):
    rng = np.random.default_rng(seed)
    ability = rng.normal(size=n)
    near = rng.binomial(1, 0.45, size=n)
    experience = rng.uniform(5, 35, size=n)
    exp_c = (experience - 20) / 10
    W = rng.normal(size=(n, k_fake))

    schooling = (
        12
        + 0.75 * near
        + 0.80 * ability
        + 0.15 * exp_c
        + rng.normal(scale=1.20, size=n)
    )
    log_wage = (
        1.5
        + BETA * schooling
        + 0.35 * exp_c
        - 0.06 * exp_c**2
        + 0.20 * ability
        + rng.normal(scale=0.35, size=n)
    )

    data = pd.DataFrame(
        {
            "log_wage": log_wage,
            "schooling": schooling,
            "exp_c": exp_c,
            "exp2": exp_c**2,
            "near_college": near,
        }
    )
    for j in range(k_fake):
        # Centering W makes the base instrument and interactions nearly orthogonal.
        data[f"zW{j + 1}"] = near * (W[:, j] - W[:, j].mean())
    return data


d = make_mincer()
d.head()
log_wage schooling exp_c exp2 near_college zW1 zW2 zW3 zW4 zW5 ... zW11 zW12 zW13 zW14 zW15 zW16 zW17 zW18 zW19 zW20
0 1.859943 10.097746 0.580055 0.336464 1 1.242903 -2.701076 0.064349 0.320079 0.668261 ... -0.210513 -0.770781 -1.463736 0.915141 0.690822 0.016899 1.372583 -0.251429 0.348562 -0.248535
1 2.201412 10.711891 -0.062175 0.003866 1 -0.054603 0.006098 -1.015384 1.209893 1.932296 ... 0.204189 1.180709 -0.884827 -0.335250 -1.737977 -0.159704 -0.138648 -1.177560 0.232402 -1.243181
2 2.617197 15.381796 -1.141093 1.302094 0 0.000000 -0.000000 0.000000 0.000000 -0.000000 ... -0.000000 0.000000 0.000000 -0.000000 0.000000 -0.000000 -0.000000 0.000000 0.000000 -0.000000
3 1.803065 12.867565 -1.345927 1.811519 1 0.184060 0.692842 -0.214432 -0.071429 -0.636782 ... 1.187301 0.397638 1.636067 -0.607905 -1.143297 -0.311634 0.343588 -0.739469 0.369532 -0.594297
4 3.071465 14.821328 0.113232 0.012821 0 0.000000 -0.000000 -0.000000 0.000000 0.000000 ... 0.000000 0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.000000 0.000000

5 rows × 25 columns

The instrument list is nested: the one-instrument specification uses only near_college; the expanded specification adds 20 interactions. PyFixest automatically includes the exogenous Mincer controls as their own instruments.

PyFixest: one instrument versus 21 moments

Show Python 
controls = "exp_c + exp2"
fake_names = [f"zW{j}" for j in range(1, 21)]

ols = pf.feols(
    f"log_wage ~ schooling + {controls}", data=d, vcov="HC1"
)
iv_one = pf.feols(
    f"log_wage ~ {controls} | schooling ~ near_college",
    data=d,
    vcov="HC1",
)
iv_many = pf.feols(
    f"log_wage ~ {controls} | schooling ~ near_college + {' + '.join(fake_names)}",
    data=d,
    vcov="HC1",
)


def pick(model, term):
    tidy = model.tidy().loc[term]
    return float(tidy["Estimate"]), float(tidy["Std. Error"])


rows = []
for label, model, term in [
    ("OLS", ols, "schooling"),
    ("2SLS: solitary valid instrument", iv_one, "schooling"),
    ("2SLS: valid instrument + 20 fake interactions", iv_many, "schooling"),
]:
    estimate, se = pick(model, term)
    rows.append({"estimator": label, "estimate": estimate, "robust_se": se})

single_sample = pd.DataFrame(rows)
single_sample["bias_from_truth"] = single_sample["estimate"] - BETA
single_sample.round(4)
estimator estimate robust_se bias_from_truth
0 OLS 0.1518 0.0059 0.0718
1 2SLS: solitary valid instrument 0.0863 0.0230 0.0063
2 2SLS: valid instrument + 20 fake interactions 0.1048 0.0214 0.0248

Calling the interactions “fake” does not mean invalid here. It means they contain no population first-stage signal beyond the solitary instrument. Their sample first-stage coefficients are noise.

T-Rex: the same IV problem as GMM

T-Rex accepts moment functions directly. For linear IV,

\[ g_i(\theta)=Z_i^*(Y_i-X_i'\theta). \]

We whiten the instrument matrix so that identity-weighted GMM uses the 2SLS weight \((Z'Z/n)^{-1}\). This makes the T-Rex objective numerically equivalent to 2SLS rather than to an arbitrarily scaled identity-weighted criterion.

Show Python 
def iv_moment(z, y, x, beta):
    return z * (y - x @ beta)[:, None]


X = np.column_stack(
    [np.ones(len(d)), d["schooling"], d["exp_c"], d["exp2"]]
)
Z = np.column_stack(
    [np.ones(len(d)), d["exp_c"], d["exp2"], d["near_college"], d[fake_names]]
)
y = d["log_wage"].to_numpy()

# If Q = Z'Z/n and W = Q^{-1} = LL', then Z* = ZL gives
# g*(theta)'g*(theta) = g(theta)'Wg(theta).
Q = Z.T @ Z / len(Z)
L = np.linalg.cholesky(np.linalg.inv(Q))
Z_white = Z @ L

trex_iv = GMMEstimator(
    iv_moment,
    weighting_matrix="identity",
    backend="scipy",
)
np.random.seed(930811)  # T-Rex currently initializes the optimizer randomly.
trex_iv.fit(Z_white, y, X, two_step=False, fit_method="BFGS")

trex_comparison = pd.DataFrame(
    {
        "implementation": ["PyFixest 2SLS", "T-Rex identity-GMM on whitened moments"],
        "schooling_return": [pick(iv_many, "schooling")[0], trex_iv.theta_[1]],
    }
)
trex_comparison["difference"] = (
    trex_comparison["schooling_return"] - trex_comparison["schooling_return"].iloc[0]
)
trex_comparison.round(7)
implementation schooling_return difference
0 PyFixest 2SLS 0.104772 0.000000
1 T-Rex identity-GMM on whitened moments 0.104739 -0.000033

The two implementations solve the same sample moment problem. PyFixest is the convenient regression interface; T-Rex exposes the GMM criterion and makes clear that every added interaction is another sample moment.

A fast 2SLS kernel for the Monte Carlo

Repeated formula parsing and numerical GMM optimization would obscure the asymptotic exercise. The Monte Carlo therefore uses the closed-form 2SLS solution. The preceding cross-check verifies that the kernel targets the same estimator.

Show Python 
def residualize(a, controls):
    return a - controls @ np.linalg.lstsq(controls, a, rcond=None)[0]


def one_replication(n, k_fake, seed):
    rng = np.random.default_rng(seed)
    ability = rng.normal(size=n)
    near = rng.binomial(1, 0.45, size=n)
    experience = rng.uniform(5, 35, size=n)
    exp_c = (experience - 20) / 10
    controls = np.column_stack([np.ones(n), exp_c, exp_c**2])
    W = rng.normal(size=(n, k_fake))

    schooling = (
        12 + 0.75 * near + 0.80 * ability + 0.15 * exp_c
        + rng.normal(scale=1.20, size=n)
    )
    wage = (
        1.5 + BETA * schooling + 0.35 * exp_c - 0.06 * exp_c**2
        + 0.20 * ability + rng.normal(scale=0.35, size=n)
    )

    s = residualize(schooling, controls)
    y_resid = residualize(wage, controls)
    instruments = np.column_stack([near, near[:, None] * W])
    instruments = residualize(instruments, controls)

    s_hat = instruments @ np.linalg.lstsq(instruments, s, rcond=None)[0]
    iv = (s_hat @ y_resid) / (s_hat @ s)
    ols_value = (s @ y_resid) / (s @ s)
    return ols_value, iv


# Closed-form parity in the displayed sample.
def closed_form_on_frame(frame, fake_names):
    C = np.column_stack([np.ones(len(frame)), frame["exp_c"], frame["exp2"]])
    s = residualize(frame["schooling"].to_numpy(), C)
    y_resid = residualize(frame["log_wage"].to_numpy(), C)
    z = residualize(frame[["near_college", *fake_names]].to_numpy(), C)
    s_hat = z @ np.linalg.lstsq(z, s, rcond=None)[0]
    return (s_hat @ y_resid) / (s_hat @ s)


pd.DataFrame(
    {
        "implementation": ["PyFixest", "closed-form kernel"],
        "estimate": [
            pick(iv_many, "schooling")[0],
            closed_form_on_frame(d, fake_names),
        ],
    }
).round(8)
implementation estimate
0 PyFixest 0.104772
1 closed-form kernel 0.104772

Experiment 1: adding interactions at fixed \(n\)

At \(n=1{,}000\), add up to 200 irrelevant-but-valid interactions. Each point averages 150 replications.

Show Python 
def run_cell(n, k_fake, reps=150, seed_offset=0):
    estimates = np.array(
        [
            one_replication(n, k_fake, 10_000_000 + seed_offset + r)
            for r in range(reps)
        ]
    )
    return {
        "n": n,
        "K_fake": k_fake,
        "K_over_n": k_fake / n,
        "ols_mean": estimates[:, 0].mean(),
        "iv_mean": estimates[:, 1].mean(),
        "iv_median": np.median(estimates[:, 1]),
        "iv_sd": estimates[:, 1].std(ddof=1),
        "iv_rmse": np.sqrt(np.mean((estimates[:, 1] - BETA) ** 2)),
    }


k_grid = [0, 5, 10, 25, 50, 100, 200]
fixed_n = pd.DataFrame(
    [run_cell(1_000, k, seed_offset=10_000 * k) for k in k_grid]
)
fixed_n.round(4)
n K_fake K_over_n ols_mean iv_mean iv_median iv_sd iv_rmse
0 1000 0 0.000 0.1523 0.0788 0.0807 0.0358 0.0357
1 1000 5 0.005 0.1523 0.0841 0.0850 0.0347 0.0348
2 1000 10 0.010 0.1521 0.0874 0.0902 0.0279 0.0288
3 1000 25 0.025 0.1516 0.0999 0.1004 0.0294 0.0354
4 1000 50 0.050 0.1517 0.1141 0.1157 0.0240 0.0416
5 1000 100 0.100 0.1512 0.1272 0.1262 0.0210 0.0517
6 1000 200 0.200 0.1520 0.1371 0.1388 0.0161 0.0593
Show Python 
fig, ax = plt.subplots(figsize=(9, 5.4))
ax.plot(fixed_n["K_fake"], fixed_n["iv_mean"], "o-", lw=2, label="Mean 2SLS")
ax.plot(fixed_n["K_fake"], fixed_n["ols_mean"], "s--", lw=1.6, label="Mean OLS")
ax.axhline(BETA, color="black", ls=":", lw=2, label="Truth: 0.08")
ax.set(
    xlabel="Number of fake interactions, K",
    ylabel="Mean estimated schooling return",
    title="At fixed n, generated instruments pull 2SLS toward OLS",
)
ax.legend(frameon=True)
plt.show()

Experiment 2: which asymptotic sequence?

“More data” is incomplete without saying what happens to the number of moments. We compare:

  • Fixed-\(K\) asymptotics: keep 10 fake interactions as \(n\) grows, so \(K/n\to0\).
  • Bekker-style many-instrument asymptotics: set \(K=\lfloor0.10n\rfloor\), so \(K/n\to0.10\).

The solitary useful instrument is unchanged in both sequences.

Show Python 
n_grid = [250, 500, 1_000, 2_000]
asymptotic_rows = []

for n in n_grid:
    for regime, k in [
        ("Fixed K = 10", 10),
        ("Many IV: K/n = 0.10", int(0.10 * n)),
    ]:
        row = run_cell(n, k, reps=150, seed_offset=1_000_000 * n + 1_000 * k)
        row["regime"] = regime
        row["bias"] = row["iv_mean"] - BETA
        row["sqrt_n_bias"] = np.sqrt(n) * row["bias"]
        asymptotic_rows.append(row)

asymptotics = pd.DataFrame(asymptotic_rows)
asymptotics[
    ["regime", "n", "K_fake", "K_over_n", "iv_mean", "bias", "iv_sd", "iv_rmse"]
].round(4)
regime n K_fake K_over_n iv_mean bias iv_sd iv_rmse
0 Fixed K = 10 250 10 0.040 0.1145 0.0345 0.0513 0.0616
1 Many IV: K/n = 0.10 250 25 0.100 0.1312 0.0512 0.0452 0.0682
2 Fixed K = 10 500 10 0.020 0.1025 0.0225 0.0399 0.0457
3 Many IV: K/n = 0.10 500 50 0.100 0.1241 0.0441 0.0297 0.0531
4 Fixed K = 10 1000 10 0.010 0.0883 0.0083 0.0317 0.0327
5 Many IV: K/n = 0.10 1000 100 0.100 0.1291 0.0491 0.0254 0.0553
6 Fixed K = 10 2000 10 0.005 0.0876 0.0076 0.0211 0.0224
7 Many IV: K/n = 0.10 2000 200 0.100 0.1254 0.0454 0.0158 0.0481
Show Python 
fig, axes = plt.subplots(1, 2, figsize=(11, 5))
colors = {"Fixed K = 10": "#0072B2", "Many IV: K/n = 0.10": "#D55E00"}

for regime, group in asymptotics.groupby("regime"):
    group = group.sort_values("n")
    axes[0].plot(
        group["n"], group["iv_mean"], "o-", lw=2,
        color=colors[regime], label=regime,
    )
    axes[1].plot(
        group["n"], group["sqrt_n_bias"], "o-", lw=2,
        color=colors[regime], label=regime,
    )

axes[0].axhline(BETA, color="black", ls=":", lw=2, label="Truth")
axes[0].set(
    xlabel="Sample size, n", ylabel="Mean 2SLS estimate",
    title="Consistency depends on K/n",
)
axes[1].axhline(0, color="black", ls=":", lw=2)
axes[1].set(
    xlabel="Sample size, n", ylabel=r"$\sqrt{n}\,(E[\hat\beta]-\beta)$",
    title="Root-n centering fails in the many-IV sequence",
)
axes[0].legend(frameon=True)
fig.tight_layout()
plt.show()

Interpretation

  1. Moment validity is not instrument strength. Under the simulation’s strong conditional independence assumptions, the interacted moments are valid. Their population first-stage coefficients are nevertheless zero.
  2. With fixed \(K\), the nuisance disappears asymptotically. The true instrument’s concentration parameter grows with \(n\), while the finite collection of spurious sample correlations becomes negligible. The 2SLS mean returns toward 0.08.
  3. With \(K/n\to\alpha>0\), overfitting does not disappear. The projection matrix spends a non-vanishing fraction of the sample degrees of freedom fitting endogenous schooling noise. Conventional 2SLS retains bias toward OLS, and \(\sqrt n(\widehat\beta-\beta)\) is not centered at zero.
  4. Interactions require stronger identifying assumptions. \(E[Z_iu_i]=0\) alone does not imply \(E[Z_iW_{ij}u_i]=0\). The expanded moments need a conditional restriction such as \(E[u_i\mid Z_i,W_i]=E[u_i\mid W_i]\), plus correctly handled covariate main effects. This notebook grants those assumptions so that many-instrument distortion is isolated from outright invalidity.
  5. A larger first-stage \(F\) is not automatically better evidence. Adding many columns can increase in-sample fit mechanically. The economically credible source of variation remains college proximity.

The practical response is not to forbid interactions. It is to state the conditional exclusion argument, limit or regularize the instrument basis, report how results change with \(K\), and use methods designed for many or weak instruments—such as LIML, Fuller corrections, jackknife IV, or identification-robust inference—when the instrument count is not small relative to the sample.

References

  • Bekker, Paul A. 1994. “Alternative Approximations to the Distributions of Instrumental Variable Estimators.” Econometrica 62 (3): 657–681.
  • Bound, John, David A. Jaeger, and Regina M. Baker. 1995. “Problems with Instrumental Variables Estimation When the Correlation Between the Instruments and the Endogenous Explanatory Variable Is Weak.” Journal of the American Statistical Association 90 (430): 443–450.
  • Card, David. 1993. “Using Geographic Variation in College Proximity to Estimate the Return to Schooling.” NBER Working Paper 4483. https://doi.org/10.3386/w4483
  • Hansen, Christian, Jerry Hausman, and Whitney Newey. 2008. “Estimation with Many Instrumental Variables.” Journal of Business & Economic Statistics 26 (4): 398–422.