|
| 1 | +"""DoubleML parity spike for the Chang (2020) Case 1 (repeated outcomes) score. |
| 2 | +
|
| 3 | +Hand-rolled DML2 cross-fitted Chang estimator vs DoubleMLDID under identical |
| 4 | +folds, learners, and clipping. This is the reproducible anchor cited by the |
| 5 | +REGISTRY "Cross-fitting, DR-score, and ridge infrastructure (DML)" section's |
| 6 | +global-p-hat Note. |
| 7 | +
|
| 8 | +Environment (side venv; doubleml/sklearn are NEVER diff-diff dependencies): |
| 9 | +
|
| 10 | + python -m venv .venv-doubleml |
| 11 | + .venv-doubleml/bin/pip install "doubleml==0.11.4" scikit-learn |
| 12 | + .venv-doubleml/bin/python benchmarks/doubleml/chang_case1_parity.py |
| 13 | +
|
| 14 | +Observed transcript (2026-08-22, doubleml 0.11.4, sklearn 1.9.0, macOS arm64): |
| 15 | +
|
| 16 | + DoubleML ATT = 3.260530717619 SE = 0.173029660048 |
| 17 | + hand Chang (global p) ATT = 3.260530717619 SE = 0.173029660048 |
| 18 | + diff vs DoubleML: ATT -4.441e-16 SE +5.551e-17 |
| 19 | + hand Chang (fold-mean p) ATT = 3.262788670382 (finite-sample gap +2.258e-03) |
| 20 | + PARITY OK (global-p ATT and SE within 1e-10 of DoubleML) |
| 21 | +
|
| 22 | +The global-p-hat convention (p-hat = full-sample treated share) matches |
| 23 | +DoubleML exactly; the fold-mean convention from the paper's proofs differs |
| 24 | +only in finite samples. The SE parity uses the augmented score |
| 25 | +psi_bar_i = summand_i - D_i * theta / p_hat. |
| 26 | +""" |
| 27 | + |
| 28 | +import numpy as np |
| 29 | +from doubleml import DoubleMLDID |
| 30 | +from doubleml.data import DoubleMLDIDData |
| 31 | +from sklearn.linear_model import LinearRegression, LogisticRegression |
| 32 | + |
| 33 | +rng = np.random.default_rng(42) |
| 34 | +N, d = 500, 5 |
| 35 | +X = rng.standard_normal((N, d)) |
| 36 | +g0 = 1 / (1 + np.exp(-(X[:, 0] - 0.5 * X[:, 1]))) |
| 37 | +D = (rng.uniform(size=N) < g0).astype(float) |
| 38 | +# Delta Y = ell(X) + theta*D + noise, theta = 3 |
| 39 | +ell0 = X[:, 0] + 0.5 * X[:, 2] ** 2 |
| 40 | +dY = ell0 + 3.0 * D + rng.standard_normal(N) |
| 41 | + |
| 42 | +K = 5 |
| 43 | +TRIM = 1e-2 # match DoubleML clipping default behavior |
| 44 | + |
| 45 | +# Fixed fold assignment shared by both implementations. |
| 46 | +perm = rng.permutation(N) |
| 47 | +folds = [np.sort(perm[i::K]) for i in range(K)] |
| 48 | +smpls = [(np.setdiff1d(np.arange(N), te), te) for te in folds] # (train, test) |
| 49 | + |
| 50 | + |
| 51 | +def _nuisances(tr, te): |
| 52 | + lg = LogisticRegression(penalty=None, solver="lbfgs", max_iter=1000).fit(X[tr], D[tr]) |
| 53 | + g_hat = np.clip(lg.predict_proba(X[te])[:, 1], TRIM, 1 - TRIM) |
| 54 | + ctrl = tr[D[tr] == 0] # outcome nuisance fit on untreated complement only |
| 55 | + lr = LinearRegression().fit(X[ctrl], dY[ctrl]) |
| 56 | + return g_hat, lr.predict(X[te]) |
| 57 | + |
| 58 | + |
| 59 | +# --- DoubleML reference ---------------------------------------------------- |
| 60 | +data = DoubleMLDIDData.from_arrays(X, dY, D) |
| 61 | +m = DoubleMLDID( |
| 62 | + data, |
| 63 | + ml_g=LinearRegression(), |
| 64 | + ml_m=LogisticRegression(penalty=None, solver="lbfgs", max_iter=1000), |
| 65 | + n_folds=K, |
| 66 | + n_rep=1, |
| 67 | + score="observational", |
| 68 | + in_sample_normalization=False, |
| 69 | + clipping_threshold=TRIM, |
| 70 | + draw_sample_splitting=False, |
| 71 | +) |
| 72 | +m.set_sample_splitting([smpls]) |
| 73 | +m.fit() |
| 74 | +att_dml, se_dml = float(m.coef[0]), float(m.se[0]) |
| 75 | + |
| 76 | +# --- hand-rolled Chang Case 1, DML2 cross-fitting -------------------------- |
| 77 | +p_glob = D.mean() # global treated share (the library convention) |
| 78 | +theta_k_glob = np.empty(K) |
| 79 | +theta_k_fold = np.empty(K) |
| 80 | +for k, (tr, te) in enumerate(smpls): |
| 81 | + g_hat, ell_hat = _nuisances(tr, te) |
| 82 | + w_num = (D[te] - g_hat) / (1 - g_hat) |
| 83 | + resid = dY[te] - ell_hat |
| 84 | + theta_k_glob[k] = np.mean(w_num * resid / p_glob) |
| 85 | + theta_k_fold[k] = np.mean(w_num * resid / D[te].mean()) |
| 86 | +theta_glob = theta_k_glob.mean() |
| 87 | +theta_fold = theta_k_fold.mean() |
| 88 | + |
| 89 | +# Variance from the augmented score psi_bar = summand - D*theta/p (global p). |
| 90 | +sig2_k = np.empty(K) |
| 91 | +for k, (tr, te) in enumerate(smpls): |
| 92 | + g_hat, ell_hat = _nuisances(tr, te) |
| 93 | + summand = (D[te] - g_hat) / (p_glob * (1 - g_hat)) * (dY[te] - ell_hat) |
| 94 | + psi_bar = summand - D[te] * theta_glob / p_glob |
| 95 | + sig2_k[k] = np.mean(psi_bar**2) |
| 96 | +se_glob = np.sqrt(sig2_k.mean() / N) |
| 97 | + |
| 98 | +print(f"DoubleML ATT = {att_dml:.12f} SE = {se_dml:.12f}") |
| 99 | +print(f"hand Chang (global p) ATT = {theta_glob:.12f} SE = {se_glob:.12f}") |
| 100 | +print(f" diff vs DoubleML: ATT {theta_glob - att_dml:+.3e} SE {se_glob - se_dml:+.3e}") |
| 101 | +print( |
| 102 | + f"hand Chang (fold-mean p) ATT = {theta_fold:.12f} " |
| 103 | + f"(finite-sample gap {theta_fold - theta_glob:+.3e})" |
| 104 | +) |
| 105 | + |
| 106 | +assert abs(theta_glob - att_dml) < 1e-10, "global-p ATT parity broken" |
| 107 | +assert abs(se_glob - se_dml) < 1e-10, "global-p SE parity broken" |
| 108 | +print("PARITY OK (global-p ATT and SE within 1e-10 of DoubleML)") |
0 commit comments