From b0c46cb9949e8ef5e069d2e983b793c1bf14b21d Mon Sep 17 00:00:00 2001 From: igerber Date: Sat, 29 Aug 2026 08:30:18 -0400 Subject: [PATCH 1/2] docs(dml): tutorial 32 - double ML DiD (numbers-locked) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fulfills the CONTRIBUTING new-estimator-checklist tutorial obligation for DMLDiD (TODO row DML PR-B1, removed here) now that survey/cluster support (#796) completed the estimator surface. docs/tutorials/32_dml_did.ipynb (16 md + 15 code cells, committed outputs; executed once via nbconvert on the pure-Python backend): - A 600-unit staggered panel DGP with nonlinear confounding in BOTH nuisances (expit-of-quadratic assignment; nonlinear covariate-dependent trend), engineered so the bias demo is robust: across 30 seeds the linear outcome learner's mean bias is +0.284 (t = 16.9) while sieve and a hand-rolled duck-typed PolynomialRidge learner are unbiased (t = 0.3). At the locked display seed: linear 2.5909 vs truth 2.2388 (4.7 SE off), sieve 2.2804 / PolynomialRidge 2.2818 (within ~1 SE). - The comparison estimand is overall_att with the DGP-implied truth computed from the fit's own per-cell n_treated weights (the panel lane's _aggregate_simple recipe). The narrative presents Chang's moment as a DISTINCT score family from the Sant'Anna-Zhao DR score (REGISTRY DR-score families note; conditional g-hat(X) vs unconditional treated share p-hat) and carries the Theorem-1 caveat: recovery works via double robustness with the propensity deliberately misspecified in every arm, so the table's SEs/CIs are illustrative, not theory-backed. - Cross-fit diagnostics (p_hat / n_clipped / fold losses + the rate-condition caveat; n_clipped=0 framed as no-trimming-required, not established overlap), event-study/group aggregation, HonestDiD (universal-base main fit; relative-magnitude restriction stated in consecutive-first-difference units; robust CI [2.1393, 2.5677] at M=1.0, qualified as an API illustration), and a varying-base bootstrap fit for sup-t bands (universal-base reference cells have structurally NaN bootstrap SEs, so the sup-t demo keeps the default base; point estimates bit-identical; nominal-coverage caveat carried). - panel=False RCS lane (A2.3 warning as the teaching beat) + a brief survey_design= mention (PSU-cohesive folds, df = n_PSU - n_strata = 16), and seed semantics via two explicit seeds (no committed seed=None cell). Registration: index.rst Advanced Methods card + toctree (meta count set to the post-add directory total 33, correcting a pre-existing off-by-one), tutorials README catalog entries for 31 (previously missing) and 32, doc-deps tutorial entries on dml_did.py / dml_did_results.py / _learners.py, CHANGELOG, and a dev-status kernelspec-note correction (10 of 33 committed notebooks carry the plain python3 kernelspec the t31/t32 hygiene tests pin). tests/test_t32_dml_did_drift.py (20 tests, 0.3s): rendered-surface quote pins (including semantic pins for the score-family distinction, RM units, and coverage qualifications), full DGP + fit re-derivation (sections 3/4/6/7 including the survey lane; the sup-t cband golden forces the NumPy multiplier-weight backend, since Rust draws a different equally-valid matrix from the same seed per the REGISTRY weight-backend identity Note — verified under both backends), ALL_CODE_CELL_HASHES source pins with a mutation negative control, source-fragment sync pins for the duplicated DGP/learner, an RCS warning gate that requires exactly one Assumption 2.3 UserWarning and fails on any other warning, and the kernelspec/no-error hygiene test. Verified: nbmake run green on DIFF_DIFF_BACKEND=python; drift suite green under both the Rust and pure-Python backends; docs-IA + doc-deps-integrity green; Sphinx -W build green. --- CHANGELOG.md | 8 + TODO.md | 1 - docs/dev-status.md | 4 +- docs/doc-deps.yaml | 6 + docs/tutorials/32_dml_did.ipynb | 1481 +++++++++++++++++++++++++++++++ docs/tutorials/README.md | 13 + docs/tutorials/index.rst | 10 +- tests/test_t32_dml_did_drift.py | 452 ++++++++++ 8 files changed, 1972 insertions(+), 3 deletions(-) create mode 100644 docs/tutorials/32_dml_did.ipynb create mode 100644 tests/test_t32_dml_did_drift.py diff --git a/CHANGELOG.md b/CHANGELOG.md index e4991162e..79b5e31bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [3.11.0] - 2026-08-29 ### Added +- **DMLDiD tutorial notebook** (`docs/tutorials/32_dml_did.ipynb`): authored via the + numbers-locked workflow (companion drift test `tests/test_t32_dml_did_drift.py`). + A staggered panel with engineered nonlinear confounding motivates the estimator: + the misspecified linear outcome learner misses by ~5 SE while the sieve and a + hand-rolled duck-typed `PolynomialRidge` learner recover the truth via double + robustness; also covers cross-fitting diagnostics, event-study/group aggregation, + HonestDiD, bootstrap sup-t bands, the `panel=False` repeated-cross-section lane + with a `SurveyDesign`, and seed semantics. - **`DMLDiD` survey-design and cluster support (both lanes)**: `fit()` gains `survey_design=` (pweight-only full-design TSL — weights/strata/PSU/FPC) and the constructor gains `cluster=` (coarser-than-unit clustering). A diff --git a/TODO.md b/TODO.md index c4124e13a..348dfa896 100644 --- a/TODO.md +++ b/TODO.md @@ -82,7 +82,6 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| DMLDiD tutorial notebook (CONTRIBUTING's new-estimator checklist requires a tutorial; deferred to its own PR per the numbers-locked notebook protocol — prototype in scripts, lock numbers, assemble + execute once, register in `docs/tutorials/index.rst` with a toctree short label + group card): staggered DGP with nonlinear covariate confounding, learner comparison (linear/ridge/sieve/sklearn object), post-fit aggregation + HonestDiD via the event-study container, seed/reproducibility note, and a `panel=False` repeated-cross-section example | `docs/tutorials/`, `docs/tutorials/index.rst` | DML PR-B1 | Mid | Medium | | Replicate Chang (2020) §4's own RCS simulation DGPs (pp. 17-21, "fully specified" per the paper review) as recovery/coverage fixtures for the `DMLDiD(panel=False)` lane — the shipped tests use a library-authored RCS design (documented in the REGISTRY checklist caveat); needs the paper PDF to extract the parameterization | `tests/test_methodology_dml_did.py`, `docs/methodology/papers/chang-2020-review.md` | DML PR-B2 | Mid | Low | | Optional scheduled end-to-end execution gate for the MMM tutorials (29/30): a cron-only workflow (or extension of `mmm-interop.yml`) that executes both notebooks in isolated exact-pin environments, so a stale/invalid committed posterior cannot stay green indefinitely - today the hybrid posture (deliberate: notebooks execute locally with committed outputs; CI smoke-tests the exporters without sampling; drift tests pin source + committed-output needles) leaves the MCMC claims un-re-executed in CI | `.github/workflows/mmm-interop.yml`, `docs/tutorials/29_mmm_calibration_pymc.ipynb`, `docs/tutorials/30_mmm_calibration_meridian.ipynb` | mmm-interop | Mid | Low | | Committed `fixest::feols` event-study golden for TWFE `event_study=True` (within + pooled specs, unbalanced + covariate panels, matched CR1 cluster convention, per-period effects + vcov block) - the in-suite gates are shared-core cross-checks (TWFE-within == MPD-absorb, pooled == MPD bit-exact), so a defect common to the shared core would pass; the live-R harness (`benchmarks/R/benchmark_multiperiod.R`, `feols(y ~ treated * time_f \| unit)`) validated the within design in `docs/benchmarks.rst` but is not a committed regression test - follow the `fixest_did_twfe_golden.json` committed-golden pattern (pytest.skip when absent) | `tests/test_fixest_did_twfe_parity.py`, `benchmarks/R/` | 3(a) R2 | Mid | Medium | diff --git a/docs/dev-status.md b/docs/dev-status.md index 444f4c6c4..96add53e5 100644 --- a/docs/dev-status.md +++ b/docs/dev-status.md @@ -170,7 +170,9 @@ Jupyter kernel (TensorFlow does not support the 3.14 dev environment): Both use the CI jobs' `diff_diff_dev.pth` shim (the repo root written into site-packages) instead of `pip install -e .` (the maturin build backend would -demand a Rust build). Committed notebooks stay kernelspec-free; every local +demand a Rust build). Committed notebooks carry either no kernelspec or the +plain `python3` one (10 of 33 carry it; the t31/t32 drift tests pin `python3` +where present) - never a machine-local venv kernel name; every local execution names the kernel explicitly (`--nbmake-kernel=...` / `nbconvert --ExecutePreprocessor.kernel_name=...`) - never the default `python3` kernel. Fragile edge: google-meridian pins an exact `tfp-nightly` diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index b4f64a023..7beb3e6c3 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -988,6 +988,8 @@ sources: - path: docs/index.rst section: "Supported Estimators (one-line catalog row)" type: user_guide + - path: docs/tutorials/32_dml_did.ipynb + type: tutorial diff_diff/dml_did_results.py: drift_risk: low @@ -998,6 +1000,8 @@ sources: section: "DMLDiD" type: methodology note: "DMLDiDResults subclasses CallawaySantAnnaResults (kit-based post-fit aggregate(), bootstrap replay); adds learner/fold/cross-fit provenance and the seed/n_bootstrap/bootstrap_weights/cband inference-provenance fields." + - path: docs/tutorials/32_dml_did.ipynb + type: tutorial # ── TROP (trop group) ────────────────────────────────────────────── @@ -1454,6 +1458,8 @@ sources: section: "Cross-fitting, DR-score, and ridge infrastructure (DML)" type: methodology note: "Duck-typed learner protocol (RegressorLearner/ClassifierLearner Protocols, validate_learner, _validate_predictions) + native learners (LinearLearner/RidgeLearner/LogitLearner/SieveLearner) wrapping linalg solvers. Contracts documented in REGISTRY: raw-X-no-intercept input, fit-reset semantics (documented limitation for stateful user learners), identified-columns prediction under rank deficiency." + - path: docs/tutorials/32_dml_did.ipynb + type: tutorial diff_diff/_dr_scores.py: drift_risk: low diff --git a/docs/tutorials/32_dml_did.ipynb b/docs/tutorials/32_dml_did.ipynb new file mode 100644 index 000000000..cbf2d59ad --- /dev/null +++ b/docs/tutorials/32_dml_did.ipynb @@ -0,0 +1,1481 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b42dc4aa", + "metadata": {}, + "source": [ + "# Double Machine Learning DiD - When Covariate Adjustment Is Too Hard to Get Right by Hand\n", + "\n", + "Staggered difference-in-differences with covariates rests on modeling two nuisance\n", + "functions: the probability of treatment given covariates (the propensity score) and the\n", + "untreated outcome trend given covariates. Callaway-Sant'Anna estimates both with fixed\n", + "parametric forms - logit and OLS. When the true relationships are nonlinear, those\n", + "models are misspecified and the \"doubly robust\" estimate quietly inherits the bias.\n", + "\n", + "**DMLDiD** (Chang 2020) keeps the same cell architecture and aggregation surface as\n", + "Callaway-Sant'Anna but plugs *cross-fitted machine-learning learners* into the\n", + "nuisances. Neyman orthogonality makes the treatment-effect estimate first-order\n", + "insensitive to nuisance estimation error, and cross-fitting (each observation scored by\n", + "learners trained on the *other* folds) removes the own-observation overfitting that\n", + "would otherwise invalidate inference.\n", + "\n", + "This tutorial:\n", + "\n", + "1. Explains when to reach for DMLDiD instead of Callaway-Sant'Anna.\n", + "2. Builds a staggered panel with *nonlinear confounding* in both nuisances.\n", + "3. Runs the first fit and reads the provenance-rich `summary()`.\n", + "4. Compares outcome learners - and watches the misspecified linear model fail.\n", + "5. Reads the cross-fitting diagnostics.\n", + "6. Aggregates to event-study and group effects, runs HonestDiD, and adds\n", + " bootstrap uniform bands.\n", + "7. Runs the repeated-cross-section lane (`panel=False`), with a survey design.\n", + "8. Shows what `seed=` does and does not pin.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "ae5eb843", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:53.564132Z", + "iopub.status.busy": "2026-08-29T11:56:53.563786Z", + "iopub.status.idle": "2026-08-29T11:56:54.417315Z", + "shell.execute_reply": "2026-08-29T11:56:54.416898Z" + } + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "from diff_diff import (\n", + " DMLDiD,\n", + " SurveyDesign,\n", + " compute_honest_did,\n", + " plot_event_study,\n", + " practitioner_next_steps,\n", + ")\n", + "\n", + "pd.set_option(\"display.precision\", 4)\n" + ] + }, + { + "cell_type": "markdown", + "id": "214fdaad", + "metadata": {}, + "source": [ + "## 1. Why ML Nuisances - and Why Cross-Fitting\n", + "\n", + "Chang's estimator targets the same cohort-period estimand and keeps the\n", + "Callaway-Sant'Anna cell architecture, but scores each `(cohort g, period t)` cell\n", + "with **Chang's own Neyman-orthogonal moment** - a distinct score family from the\n", + "Sant'Anna-Zhao doubly-robust score CS uses (Chang normalizes by the unconditional\n", + "treated share $\\hat p_g$ and carries its associated augmentation term; the\n", + "repeated-cross-section case adds further $\\lambda$-related terms). The two\n", + "plug-in nuisances\n", + "\n", + "- $\\hat g_{g}(X)$ - the conditional propensity of belonging to cohort $g$ given covariates, and\n", + "- $\\hat\\ell_{gt}(X) = E[\\Delta Y \\mid X, \\text{control}]$ - the control outcome-change\n", + " regression,\n", + "\n", + "come from **cross-fitted learners**: the sample is split into $K$ folds, and each\n", + "observation's nuisance predictions come from learners trained on the other $K-1$\n", + "folds. Two properties carry the theory (Chang 2020, Thm. 1-2):\n", + "\n", + "- **Neyman orthogonality** - the moment's derivative with respect to each nuisance is\n", + " zero at the truth, so small nuisance errors have only second-order effect on the ATT;\n", + "- **sample splitting** - predictions never come from a learner that saw the\n", + " observation, so flexible learners cannot overfit their own scores.\n", + "\n", + "*When to prefer it over `CallawaySantAnna`*: covariate relationships that are\n", + "nonlinear or high-dimensional. When a logit + linear-OLS specification is plausible,\n", + "prefer CS - fewer moving parts, replicate-weight support, and no learner tuning (see\n", + "*Choosing an estimator* in the docs). DMLDiD **requires** covariates: without them the\n", + "learners have nothing to do, and `fit()` points you back to CS.\n" + ] + }, + { + "cell_type": "markdown", + "id": "911d9e80", + "metadata": {}, + "source": [ + "## 2. A Staggered Panel with Nonlinear Confounding\n", + "\n", + "We build a 600-unit, 6-period panel with cohorts treated at $t=4$ and $t=5$ plus a\n", + "large never-treated pool, and we bury a *quadratic* confounder in both nuisances:\n", + "\n", + "- **assignment**: treatment probability is an expit of $0.9(x_1^2 - 1) + 0.7 x_1 x_2$\n", + " - so a logistic model in *raw* $(x_1, x_2)$ is misspecified;\n", + "- **trend**: untreated outcomes drift along $f(X)\\, t/6$ with\n", + " $f(X) = 1.6(x_1^2 - 1) + x_1 x_2$ - a covariate-dependent trend that survives\n", + " differencing (levels difference out in DiD; only the trend confounds), and is\n", + " nonlinear in exactly the way a linear outcome learner cannot absorb.\n", + "\n", + "The true effect is dynamic and known by construction:\n", + "$\\tau(e) = 2.0 + 0.3\\,e$ at event time $e = t - g$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "41f5b35e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.418612Z", + "iopub.status.busy": "2026-08-29T11:56:54.418512Z", + "iopub.status.idle": "2026-08-29T11:56:54.428801Z", + "shell.execute_reply": "2026-08-29T11:56:54.428404Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "first_treat\n", + "0 354\n", + "4 121\n", + "5 125\n", + "Name: count, dtype: int64" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rng = np.random.default_rng(7)\n", + "n_units, periods = 600, [1, 2, 3, 4, 5, 6]\n", + "\n", + "x1 = rng.normal(size=n_units)\n", + "x2 = rng.normal(size=n_units)\n", + "\n", + "# Assignment: expit of a QUADRATIC score -> logit-on-raw-X is misspecified.\n", + "# Capped at 0.60 so a large never-treated pool always survives (the default\n", + "# control_group=\"never_treated\" needs one).\n", + "score = 0.9 * (x1**2 - 1.0) + 0.7 * x1 * x2\n", + "p_any = np.clip(1.0 / (1.0 + np.exp(-score)), 0.05, 0.60)\n", + "treated = rng.uniform(size=n_units) < p_any\n", + "early = rng.uniform(size=n_units) < 0.5\n", + "cohort = np.where(treated, np.where(early, 4, 5), 0)\n", + "\n", + "# Untreated outcome: nonlinear covariate-dependent TREND (this is what\n", + "# confounds a DiD; time-invariant levels difference out).\n", + "f_nl = 1.6 * (x1**2 - 1.0) + 1.0 * x1 * x2\n", + "rows = []\n", + "for i in range(n_units):\n", + " alpha_i = 0.5 * x1[i] - 0.3 * x2[i] + rng.normal(scale=0.3)\n", + " for t in periods:\n", + " y = 1.0 + 0.25 * t + alpha_i + f_nl[i] * (t / 6) + rng.normal(scale=0.5)\n", + " if cohort[i] > 0 and t >= cohort[i]:\n", + " y += 2.0 + 0.3 * (t - cohort[i]) # true tau(e) = 2.0 + 0.3 e\n", + " rows.append((i, t, y, cohort[i], x1[i], x2[i]))\n", + "df = pd.DataFrame(rows, columns=[\"unit\", \"time\", \"y\", \"first_treat\", \"x1\", \"x2\"])\n", + "\n", + "df.groupby(\"unit\")[\"first_treat\"].first().value_counts().sort_index()\n" + ] + }, + { + "cell_type": "markdown", + "id": "e7b974a6", + "metadata": {}, + "source": [ + "354 never-treated units against cohorts of 121 ($g=4$) and 125 ($g=5$).\n", + "\n", + "## 3. First Fit\n", + "\n", + "The propensity learner has exactly one built-in spec (`\"logit\"`); the outcome learner\n", + "is where the flexibility lives - `\"sieve\"` is an adaptive polynomial whose degree is\n", + "selected by BIC per cell. `seed=` pins the fold draws, and\n", + "`base_period=\"universal\"` benchmarks every pre-period against $t = g-1$ (the default\n", + "is `\"varying\"`; universal is what HonestDiD in section 6 needs).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e9d02484", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.429813Z", + "iopub.status.busy": "2026-08-29T11:56:54.429757Z", + "iopub.status.idle": "2026-08-29T11:56:54.474313Z", + "shell.execute_reply": "2026-08-29T11:56:54.473986Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "=====================================================================================\n", + " DML DiD (Chang 2020) Staggered Difference-in-Differences Results \n", + "=====================================================================================\n", + "\n", + "Propensity learner: 'logit'\n", + "Outcome learner: 'sieve'\n", + "Cross-fitting folds (K): 5\n", + "Seed: 42\n", + "\n", + "Total observations: 3600\n", + "Treated units: 246\n", + "Never-treated units: 354\n", + "Treatment cohorts: 2\n", + "Time periods: 6\n", + "Control group: never_treated\n", + "Base period: universal\n", + "\n", + "-------------------------------------------------------------------------------------\n", + " Overall Average Treatment Effect on the Treated \n", + "-------------------------------------------------------------------------------------\n", + "Parameter Estimate Std. Err. z-stat P>|z| Sig.\n", + "-------------------------------------------------------------------------------------\n", + "ATT 2.2804 0.0438 52.053 0.0000 ***\n", + "-------------------------------------------------------------------------------------\n", + "\n", + "95% Confidence Interval: [2.1945, 2.3662]\n", + "CV (SE/abs(ATT)): 0.0192\n", + "\n", + "Signif. codes: '***' 0.001, '**' 0.01, '*' 0.05, '.' 0.1\n", + "=====================================================================================\n" + ] + } + ], + "source": [ + "fit_kw = dict(\n", + " outcome=\"y\", unit=\"unit\", time=\"time\", first_treat=\"first_treat\",\n", + " covariates=[\"x1\", \"x2\"],\n", + ")\n", + "\n", + "est = DMLDiD(\n", + " propensity_learner=\"logit\", # the only built-in propensity spec\n", + " outcome_learner=\"sieve\", # adaptive polynomial, BIC-selected degree\n", + " n_folds=5,\n", + " seed=42, # pins the fold draws\n", + " base_period=\"universal\",\n", + ")\n", + "res = est.fit(df, **fit_kw)\n", + "print(res.summary())\n" + ] + }, + { + "cell_type": "markdown", + "id": "6de2565e", + "metadata": {}, + "source": [ + "The header records the full nuisance provenance - learners, folds, seed - so a\n", + "serialized result is self-describing.\n", + "\n", + "**What should this number be?** `overall_att` is the simple aggregation: an average of\n", + "the post-treatment cell effects $ATT(g,t)$ weighted by each cell's treated count. The\n", + "DGP knows every cell's true effect, so we can compute the exact target by applying\n", + "the *fit's own* aggregation weights to the true per-cell effects:\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "829cfa1e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.475460Z", + "iopub.status.busy": "2026-08-29T11:56:54.475396Z", + "iopub.status.idle": "2026-08-29T11:56:54.477740Z", + "shell.execute_reply": "2026-08-29T11:56:54.477340Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DGP-implied overall ATT: 2.2388\n", + "DMLDiD (sieve) estimate: 2.2804 +/- 0.0438\n" + ] + } + ], + "source": [ + "BASE, SLOPE = 2.0, 0.3\n", + "num = den = 0.0\n", + "for (g, t), cell in res.group_time_effects.items():\n", + " if t < g or cell.get(\"is_reference\"):\n", + " continue # overall_att averages post-treatment cells only\n", + " w = cell[\"n_treated\"] # the fit's own per-cell aggregation weight\n", + " num += w * (BASE + SLOPE * (t - g))\n", + " den += w\n", + "true_att = num / den\n", + "print(f\"DGP-implied overall ATT: {true_att:.4f}\")\n", + "print(f\"DMLDiD (sieve) estimate: {res.overall_att:.4f} +/- {res.overall_se:.4f}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "46da3845", + "metadata": {}, + "source": [ + "The sieve fit lands on **2.2804 ± 0.0438** against a true **2.2388** - within one\n", + "standard error.\n", + "\n", + "## 4. Learner Comparison - Watching Misspecification Fail\n", + "\n", + "Now the punchline. We refit with four outcome learners, keeping everything else\n", + "identical. One of them is a **custom learner object**: DMLDiD accepts any object with\n", + "`fit(X, y, sample_weight=None)` and `predict(X)` (classifiers: `predict_proba`).\n", + "Here is a numpy-only ridge regression on degree-2 polynomial features - flexible\n", + "enough to represent the DGP's quadratic trend:\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b6496708", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.478753Z", + "iopub.status.busy": "2026-08-29T11:56:54.478685Z", + "iopub.status.idle": "2026-08-29T11:56:54.481343Z", + "shell.execute_reply": "2026-08-29T11:56:54.481007Z" + } + }, + "outputs": [], + "source": [ + "class PolynomialRidge:\n", + " \"\"\"Degree-2 polynomial features + ridge normal equations (numpy only).\"\"\"\n", + "\n", + " def __init__(self, alpha=1.0):\n", + " self.alpha = alpha\n", + "\n", + " def _expand(self, X):\n", + " n, d = X.shape\n", + " cols = [np.ones(n)] + [X[:, j] for j in range(d)]\n", + " for j in range(d):\n", + " for k in range(j, d):\n", + " cols.append(X[:, j] * X[:, k])\n", + " return np.column_stack(cols)\n", + "\n", + " def fit(self, X, y, sample_weight=None):\n", + " Z = self._expand(np.asarray(X, dtype=np.float64))\n", + " y = np.asarray(y, dtype=np.float64)\n", + " w = np.ones(len(y)) if sample_weight is None else np.asarray(sample_weight, float)\n", + " ZtW = Z.T * w\n", + " A = ZtW @ Z + self.alpha * np.eye(Z.shape[1])\n", + " A[0, 0] -= self.alpha # leave the intercept unpenalized\n", + " self.coef_ = np.linalg.solve(A, ZtW @ y)\n", + " return self\n", + "\n", + " def predict(self, X):\n", + " return self._expand(np.asarray(X, dtype=np.float64)) @ self.coef_\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "8e3f77d0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.482238Z", + "iopub.status.busy": "2026-08-29T11:56:54.482178Z", + "iopub.status.idle": "2026-08-29T11:56:54.610519Z", + "shell.execute_reply": "2026-08-29T11:56:54.610180Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
outcome learnerATTSE|bias| / SE
0linear2.59090.07434.7392
1ridge2.58980.07364.7720
2sieve2.28040.04380.9482
3PolynomialRidge (custom)2.28180.04380.9823
\n", + "
" + ], + "text/plain": [ + " outcome learner ATT SE |bias| / SE\n", + "0 linear 2.5909 0.0743 4.7392\n", + "1 ridge 2.5898 0.0736 4.7720\n", + "2 sieve 2.2804 0.0438 0.9482\n", + "3 PolynomialRidge (custom) 2.2818 0.0438 0.9823" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "comparison = []\n", + "for name, learner in [\n", + " (\"linear\", \"linear\"),\n", + " (\"ridge\", \"ridge\"),\n", + " (\"sieve\", \"sieve\"),\n", + " (\"PolynomialRidge (custom)\", PolynomialRidge(alpha=1.0)),\n", + "]:\n", + " r = DMLDiD(outcome_learner=learner, n_folds=5, seed=42,\n", + " base_period=\"universal\").fit(df, **fit_kw)\n", + " comparison.append((name, r.overall_att, r.overall_se,\n", + " abs(r.overall_att - true_att) / r.overall_se))\n", + "pd.DataFrame(comparison, columns=[\"outcome learner\", \"ATT\", \"SE\", \"|bias| / SE\"])\n" + ] + }, + { + "cell_type": "markdown", + "id": "84774620", + "metadata": {}, + "source": [ + "Both misspecified learners miss high: **linear 2.5909** and **ridge 2.5898** sit\n", + "almost five standard errors above the true 2.2388 (a ridge penalty on the wrong\n", + "functional form is still the wrong functional form). The flexible learners recover\n", + "it: **sieve 2.2804** and the custom **PolynomialRidge 2.2818**, each within about\n", + "one standard error.\n", + "\n", + "Why does fixing only the *outcome* learner rescue the estimate when the propensity\n", + "model (`\"logit\"` on raw $X$) is misspecified in **every** arm? That is **double\n", + "robustness** at work: Chang's moment is consistent when *either* nuisance is right,\n", + "and the flexible outcome learners get theirs right. Two honest caveats:\n", + "\n", + "- Double robustness protects the *point estimate*. Chang's Theorem 1 inference needs\n", + " **both** nuisances converging at $o(N^{-1/4})$ - \"a fast learner cannot compensate\n", + " a slow one\" - so with a deliberately misspecified propensity, the SEs and CIs in\n", + " this table are illustrative rather than theory-backed.\n", + "- Custom learner objects are **deep-copied, never-fit, once per fold** - so they must\n", + " survive `copy.deepcopy`; under `survey_design=` (section 7) the `sample_weight`\n", + " keyword becomes mandatory, and stochastic learners need their own internal seeding\n", + " (`seed=` pins folds, not your learner's RNG).\n" + ] + }, + { + "cell_type": "markdown", + "id": "15d70cb8", + "metadata": {}, + "source": [ + "## 5. Cross-Fitting Diagnostics\n", + "\n", + "Every fit records per-cell cross-fitting health in `cross_fit_diagnostics`: the\n", + "cell's treated share `p_hat`, the number of propensity scores clipped by the trim\n", + "(`pscore_trim=0.01` by default), and out-of-fold losses per fold (log-loss for the\n", + "propensity learner, MSE for the outcome learner).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "c2ddc9f9", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.611538Z", + "iopub.status.busy": "2026-08-29T11:56:54.611478Z", + "iopub.status.idle": "2026-08-29T11:56:54.615289Z", + "shell.execute_reply": "2026-08-29T11:56:54.614916Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
gtp_hatn_clippedpropensity log-lossoutcome MSE
0410.254700.57670.5145
1420.254700.57720.5026
2440.254700.57641.6779
3450.254700.57591.7839
4460.254700.57802.3770
5510.261000.57060.4929
6520.261000.57580.4782
7530.261000.57540.5367
8550.261000.57201.5103
9560.261000.57191.8739
\n", + "
" + ], + "text/plain": [ + " g t p_hat n_clipped propensity log-loss outcome MSE\n", + "0 4 1 0.2547 0 0.5767 0.5145\n", + "1 4 2 0.2547 0 0.5772 0.5026\n", + "2 4 4 0.2547 0 0.5764 1.6779\n", + "3 4 5 0.2547 0 0.5759 1.7839\n", + "4 4 6 0.2547 0 0.5780 2.3770\n", + "5 5 1 0.2610 0 0.5706 0.4929\n", + "6 5 2 0.2610 0 0.5758 0.4782\n", + "7 5 3 0.2610 0 0.5754 0.5367\n", + "8 5 5 0.2610 0 0.5720 1.5103\n", + "9 5 6 0.2610 0 0.5719 1.8739" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "diag_rows = []\n", + "for (g, t), d in sorted(res.cross_fit_diagnostics.items()):\n", + " if \"p_hat\" not in d:\n", + " continue # universal-base reference cells carry no learner fits\n", + " diag_rows.append((g, t, d[\"p_hat\"], d[\"n_clipped_ps\"],\n", + " np.nanmean(d[\"propensity\"][\"fold_losses\"]),\n", + " np.nanmean(d[\"outcome\"][\"fold_losses\"])))\n", + "pd.DataFrame(diag_rows, columns=[\n", + " \"g\", \"t\", \"p_hat\", \"n_clipped\", \"propensity log-loss\", \"outcome MSE\"])\n" + ] + }, + { + "cell_type": "markdown", + "id": "64fea9dd", + "metadata": {}, + "source": [ + "`n_clipped = 0` in every cell means no trimming was required - that alone does not\n", + "establish overlap or valid nuisance rates, but overlap trouble does announce itself\n", + "loudly: a nonzero clip count **plus** a `UserWarning` at fit time means some units'\n", + "propensities were pushed against the trim boundary (clip-never-drop: DMLDiD warns\n", + "and clips rather than silently discarding units).\n", + "\n", + "One caution the numbers cannot give you: Chang's theory requires both nuisances to\n", + "converge at rate $o(N^{-1/4})$, and **the fold losses do not verify that** - they\n", + "measure predictive fit, not convergence rates. A slow learner silently invalidates\n", + "normal inference while still posting respectable losses.\n", + "\n", + "## 6. Aggregation, HonestDiD, and Uniform Bands\n", + "\n", + "The staggered surface is post-fit, exactly like Callaway-Sant'Anna:\n", + "`aggregate('event_study')`, `aggregate('group')`, `aggregate('simple')`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "359ea48e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.616211Z", + "iopub.status.busy": "2026-08-29T11:56:54.616154Z", + "iopub.status.idle": "2026-08-29T11:56:54.709486Z", + "shell.execute_reply": "2026-08-29T11:56:54.709180Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAJOCAYAAABBfN/cAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAa19JREFUeJzt3QmcTfX/+PH37IuZYRhMWYuKkC0kbaQiRItKpZGQvr7t9a2+9a1vv3+LFqWi1VckSnyLUlSkJFGyFIWyZClmjBlmzD5z/4/3p++dZrkz5o573XPueT0f3T73nnvumc9ZnHPe57OFuFwulwAAAAAAAL8I9c9iAQAAAAAAgTcAAAAAAH5GiTcAAAAAAH5E4A0AAAAAgB8ReAMAAAAA4EcE3gAAAAAA+BGBNwAAAAAAfkTgDQAAAACAHxF4AwAAAADgRwTeAAC/2L9/v4SEhMjEiRPZwjYSHR0tt99+u0+XmZqaKgkJCfL+++979TunHkNOXW+1Z88eiYmJkUWLFgU6KwDgUwTeAGAxRUVF5qbb/QoLC5N69epJ+/btZcSIEbJ48eJqfzN06FCPy50xY0bpPPPmzav026uvvtqnefLncuxg9+7dZh0nTZpk6WUeCw8++KCcdNJJMmTIkEBnBRbXpEkTGTt2rNx1113mfAEAwYLAGwAs6qqrrhKXyyXFxcWyc+dOmTp1qsTHx8tFF11kguvCwsJKv6lTp458+OGHcuDAgUrfTZs2zXx/rPKUlJRk5vVUelqbdYM9/fHHH/LGG2/I3//+d/PQwBvVHUPBzKnr7TZu3Dj56aefZP78+YHOCgD4DIE3ANiAVtPt3r27vPjii/LWW2/J3Llz5b777qs0X//+/SU0NFTefvvtctN/++03Wbp0qVxxxRXHPE/Hajmwpv/85z8SHh7u02MPwa1169bSo0cPeeWVVwKdFQDwGQJvALCZYcOGmUB18uTJkp2dXe67unXrmuq8Wrpdln7WUrQBAwYcszzVpp1qdetWkZaKP/HEE3Lqqaeadsn169c3Jek7duww3x8+fNgE9VqFvaKCggJp2LChmb+my1O//vqrWacpU6bIf//7X2nbtq1ERUXJaaedJp988knpfMuXL5dmzZqZ97fccktp1Xot9a3O888/Lx06dJC4uDg5/vjjzb5cu3ZtjZb5xRdfmM+e2sZq4Hv33XeXm7Zv3z657rrrTFV/fV1//fVy6NChcvPoNtTvdL6K8vLyzDF15ZVXVrtOWgOjW7dupkaDN+tb3TEUiH3vSU2PB2/3T1XrfaTt5V6W+9iIjIyUE088Ue655x6zTWqTb/f2efLJJ83f1vbXTZs2lRtvvNG0x/Z2G9ZkHVSfPn3MNsvKyqp2HwCAXRB4A4ANXXjhhZKfny8rV66s9N0NN9wgq1evlo0bN5rPWmV1+vTpcu2110pERERA8uTr5ZSUlMill14qzz77rDz88MOyd+9e+fbbb81N+plnnmk689Jq9XrjP2fOnEoB5QcffGCCGw0earq8sj799FNZtmyZfPbZZyawOOGEE8zvNZhVZ511luzatcu815J83Qf6qq5t9quvvmoCJG0PrdWzf/zxRxk1apQ888wztV5mVTRoPv/882XFihUmMNbq/pdffrlpW1uWbkMNXrUWQlpaWrnvZs+eLenp6aXb0JPc3FxZs2aNnH766V6vr1X3vSdHOh58oabbS9tFu48N3T/6EGvmzJmm+nZt8q3NQQYNGiTjx483Dwf0WPnuu+/k7LPPltdff92rbejNPtcHcLoueowCQFBwAQAspbCw0KWn56uuuqrKeV555RUzz/Tp08v95sYbb3QVFxe7mjVr5rr77rvNd59//rn5bv369a7333/fvNfUm79XmzylpaWZz88999xRLceT2bNnm3nmzJlTbnpmZqarXr16rnvvvdd8XrlypZnv1VdfLTdfv379XM2bNzfbypvl/fLLL2a+7t27l5vvt99+M9MnTJhQOm3Xrl1m2osvvuiqiauvvtrVtm3bauepbplLly413y1cuLDSd2FhYa677rqr9LNuD51Xf+Np2992222l0zZv3uwKCQlxjR8/vty8PXr0MMeZext68uuvv5rlPfXUU7VaX0/HUKD2vSfeHA/e7B9P612T7VWVSZMmuUJDQ105OTle5/utt94y02bNmlXl8mu6Db1Zh6r2HwDYFSXeAGBDWpqlPHVWpW28hw8fbtpLa2mVdmzVpUsXU400UHny9XK0lFar0V5yySWVqtpr6eqXX35pPms7Ue0xXTtvK9szuJb0ac0A3VbeLM+tYpX95s2bm6rU27Ztq/V6d+zYUX7++WfTodb3339v9p2/LFmyxOT3vPPOKzfdU6/jJ598sikd19JKLdlUWjV41apV5bahJ5mZmSb1VM28tusb6H3viT+Oh9puL10/rTXSoEEDs47u5gi67yrmpyb5/vjjj0319apGS/BmG3qzz7WpQNljCADsjsAbAGxIAwilbSQ90erBWt1Tq9pqG05PbV2PdZ58uRxdN22rGxsba4ICHZbMHWTokGRaxdZNqxRrkOiueq/t3TW41+CrNstTxx13nMdA4WiCBB0+6d///rfpyVmDFQ2ctEMyrap9tNwPM9x0fRo3blxpvkaNGnkMpP/2t7/J9u3bS9snv/TSS2bblN2Gnmj7cFWxuvfRrG+g970nR3s8VNw/ntRke2m1ce1gsWXLlqaphlb112W7+3yoOFpATfKt1cT1WNFtU5WabkNv9rn7mElMTDzitgEAOyDwBgAb0lIt7cDojDPO8Pi9jpmsbSs1YNJ2ktdcc03A8+TL5WinXhog6M2+rp+WnGmJnrtt65YtW0rn1dJ/LY3Tkk93ENK3b19p0aJFrZbni1J9T7T9vbaP1QBXX88995xs2LBBzjnnnNKHEdXR0kVVsTMqDYrcJdVuGvB4arus0yrOq7QkUzt204D74MGDMmvWLFMKrgFedfThiQZimgdfrW+g970nNTkevNk/ntRke2ktF+38THsD13OAdpimdP7a5ls7otM239WNqV3TbejNPtc24KrsvgIAOyPwBgCb0aHCtHMj7dm6unG5tVQvIyPDdIykgZYV8uSr5eg6aYmYpx6iK9J11yrUM2bMMEH91q1bK3UI5s3yasqdf+0ozlsa0Or+006ntDdqd4lgdcvUjrG0lFEDmbIWLFjgscdoXd+K1airGjdZSzDHjBkjCxcuNCWWOTk51Xaq5qZBYNeuXU1nf7VZX7vue0+82T9HUt320gcNZQNqDYL1QUltDRw40ATTWnumKrXZhkfa59o5mz606dmzZ63zDgBWQuANADagQ2tpQHrrrbeaUjwdwumxxx6r9jfaW7CWNmmP1FbJk6+Wc/XVV5uAQKvQv/nmm6bUUEsS9eZdxwCv2EuyBlvaK/fo0aPNMEcV2zJ7u7ya0CqyWuL7+eefm1LiI0lJSTG9lW/evNkE1jr2upbQ6rBL7l7Bq1umVuvWartaKv3VV1+Z/L/33ntm22rgXJYOHaZDSOl2+frrr8282k5XA3F3KWlFuu00ENIhrnQbai/WNaFBmeahYklvTdbXrvveE2/2jyc12V5aM0Eftv3rX/8yx4cOG6ZD9Gnb6trS3uEvuugiU3tGR0fQbanbSP+2PoTxZht6s8/1GNc+CDz1DwAAthTo3t0AAOW5e/52v7Q34vj4eNepp57qSklJcS1evLjSJivbq3l1quvVvKrX4MGDa5Wn6no192Y5VSkqKnJNnDjR1aVLF1dsbKyrbt26rtNPP9319NNPm96UyyopKXG1bNnS/M1bb7211stz9wb9+uuvV/p9kyZNXNdee225aYsWLXK1b9/eFRkZaX43bty4Ktdn27ZtrjvuuMPVpk0bV3R0tOv44493DR061LV27doaL1O3uf4mLi7O9CY9ZswY05N1xV6z1R9//GF6mU5ISDCva665xqxnVFRUuV7Nyxo2bJj5m7fccourpvTvREREuKZOner1+no6hgK17z3x9nio6f7xtN41PT60F/CTTz7ZzKM9iL/55puut99+2yzPPa+3+c7Pz3c99thjZnl63DVt2tSca/bs2ePVNqzpOrjzN3fu3Cq3PQDYTYj+L9DBPwAAsD4d51t7N1+3bp1Xpaj6Oy3Z1Z6sgSO54447TNOA9evXV9upGwDYCYE3AAA4Im3nq71ga1tlbX/rDe20TTv70qHtLrvsMrY2qvT7779Lq1at5P3335d+/fqxpQAEDR4jAgCAIwbdEyZMkP3798uUKVO83lo6TFlN2rkD2oeBDoMGAMGGwBsAAFTpnXfeMR106ZBR//d//yeDBw9mawEA4CWqmgMAAAAAEMwl3tpxhnagocPJtGvXzgxPEhERUe1v9u3bZ8ak1LRDhw5mqAs63wAAAAAAWFFAx/G+++67TU+nmZmZ5vNDDz0k3bp1k0OHDlX5m19++cUE24sWLTIB+sMPP2w63yguLj6GOQcAAAAAwAZVzTdt2iRt2rQp/aydtjRu3Fjeeust057ME+0NVef74osvJDQ0VHbu3CmtW7eW//znPzJ8+PAa/d2SkhLTa2Z8fLyEhIT4bH0AAAAAAM7gcrkkKyvLdAypsallq5qXDbqV9mKpmU9MTPQ4f2FhoXz88ccyceLE0hVr3ry5nHfeeTJv3rwaB94adDdr1swHawAAAAAAcLJdu3ZJ06ZNrd3GW0u9X3nlFVO9fMWKFfLYY49VOW6jlm7n5+eb8R3L0s9ff/11lX9Df6MvN3ch/2+//SYJCQliRVoqn5eXZ95TKg8AcJqCggLTl4vWhIuMjBSncd+rREdHH7EUBQAQGBrDtmjRwtSkPpKAB94xMTHSsmVLU31cg00NoLWjtbi4uErz5uTkmLTiimnw7P7OkyeeeEIeeeSRStM1GHcHt1aj20JL+BWBNwDAafTal5ycbN67r4dOUrYlIIE3AFiTu3C3JvGapYYTO3jwoJx88slyyy23yIMPPljp+x07dsgJJ5wgCxcuLFcqPmbMGPnuu+9k7dq1NSrx1icTWtU8IyPD8iXeuhMJvAEATqMP1Ddu3GhGPImNjRWn0dszfVHiDQDWpXGlNpPWOPZIcWXAS7zLqlu3rmn3rdXPPdH23FoSrt+XDbz186mnnlrlcqOiosyrIn2CbOWnyJo3Am8AgBMdPnzYNEHTWnGeasE5JfC2+r0KADhZqBfn59BAtt1avnx5uWnbt283pdZdu3YtnTZ79mx58sknS1fs8ssvl2nTppVWEf/hhx9M9fQrr7zyGK8BAADwF23bfccdd5gUAAC7C1hVc22vdfHFF5uezLW0Wsfy1irk/fv3lxkzZpSWUI8aNUpWrlwpGzZsMJ/37t0r5557rulopXPnzqaX8wEDBsj06dO9qhKgpes1qRIQKFQ1BwDAuahqDgDW501cGfA23mvWrDGl3HXq1JFOnTpVGmJs8eLFJti+7rrrSqdpsK5BuvZ22qFDBznrrLO8+psE3gAAWFt6enrpA/kGDRqI0xB4A4D1eRNXBryNd5cuXcyrKn379vXYE/pll13m55wBAIBACQ8PNx3WaAoAgN1xNQMAAJajJQjalAwAgGBAN5kAAMByiouLzZBimgIAYHcE3gAAwHLS0tLkxRdfNCkAAHZH4A0AACynXr16pj8XTQEAsDvaeAMAAMuJjo6Wk046KdDZAADAJyjxBgAAlqPtu9etW2dSAADsjsAbAABYcmzUTz/91KQAANhdiMvlconDeDPQeaCUlJRIXl6ehISEmBcAAHAOvT3Tl1a5Dw2lnAQA7B5XciYHAAAAAMCPCLwBAIDlHDhwQGbPnm1SAADsjsAbAABYjlavjoqKopo1ACAoMJwYAACwHB2/e8iQIYHOBgAAPkHgDQAALNnJaFFRkYSHh1PqDQBBLCffu76+Y6Ps2fE0gTcAALCc1NRUmT59uqSkpEhycnKgswMA8JNO49K9mn/LlCRb7gvaeAMAAMvR4VkGDRpkUgAA7I4SbwAAYDkxMTFy6qmnBjobAAA/Wze5QbnPufku6XnnnyNafPNsfYmxadXyigi8AQCA5eTm5sq2bdvkxBNPNEE4ACA4xVYTWGvQbdc23RVR1RwAAFjOwYMHZcGCBSYFAMDuKPEGAACW07hxY7nrrrskLCws0FkBAOCoEXgDAADLCQkJMUOJAQAQDKhqDgAALCczM1Pef/99kwIAYHcE3gAAwHJcLpcUFRWZFAAAu6MOFwAAsJzExEQZOnRooLMBAIBPUOINAAAAAIAfEXgDAADL2bt3rzz55JMmBQDA7gi8AQCA5SQkJEi/fv1MCgCA3dHGGwAAWE5sbKx07Ngx0NkAAMAnKPEGAACWk5eXJ5s3bzYpAAB2R+ANAAAsR8fvnjdvHuN4AwCCAlXNAQCA5TRq1EhuvfVWiYqKCnRWAAA4agTeAADAckJDQyUmJibQ2QAAwCeoag4AACxZ1fzDDz+kqjkAOEh+oUsWrMov/XzX61ky75s8M93uKPEGAACWU1JSIllZWSYFAAS/Jevy5d6p2XIo568ge+n6AlmyrkAeffuwPDUyTvp0sm/zI0q8AQCA5dSvX1+uueYakwIAgj/o/tvkLMkqE3Srkv991Ok3T84y89kVgTcAAAAAICDyC12mpFtc5j+PzHSXmPnsWu2cwBsAAFjOvn375NlnnzUpACB4LVydb6qXHymc1u91vkWr7VnqTeANAAAsJy4uTs455xyTAgCC1+K1BRIaUrN5db7P1haIHdG5GgAAsJw6derI6aefHuhsAAD8LDO7pLQt95HofJmH7dnpJiXeAADAcvLz82Xbtm0mBQAEr3pxoV6VeNerY88Q1p65BgAAQS0jI0PmzJljUgBA8OrbOdKrEu8LOkeKHRF4AwAAy2nYsKHcfPPNJgUABK/+p0dJQmyIHKnQW7/X+fqdbs+xvAm8AQCA5YSFhUlCQoJJAQDBKyoiRJ4aGWci66qCbzM9RMx8Or8dEXgDAADLOXTokCxatMikAIDg1qdTlLw0Ll7iY8sH1e623zr95XHxZj67oldzAABgOYWFhWYMb00BAMHv/E5R8vWESPlgZb48MD3bTOvdMVL6dY001cvtWtLtFuJyuWrYlD146NPzunXrysGDB001NisqKSmRvLw8CQkJMS8AAOAcenumr+joaAkNpYIiAOfIyXdJp3Hp5v26yQ0kNiokKOJKzuQAAAAAAPgRgTcAALCc1NRUeeGFF0wKAIDdEXgDAADLiY2NlW7dupkUAAC7o3M1AABgOXFxcdKzZ89AZwMAAJ+gxBsAAFhOQUGB7N6926QAANgdgTcAALCcAwcOyMyZM00KAIDdUdUcAABYTlJSkowaNcoM0wIAgN0ReAMAAMsJDw+XBg0aBDobAAD4BFXNAQCA5Rw6dEiWLFliUgAA7I7AGwAAWI52qrZ9+3Y6VwMABAWqmgMAAMu28QYAIBhQ4g0AAAAAgB8ReAMAAMtJS0uTl19+2aQAANgdVc0BAIDlREdHS7t27UwKAAheOfmucp9zy3wu+94tNipE7CjE5XJVXpsgpz2k6rigBw8elISEBLGikpISycvLk5CQEPMCAADOobdn+tIHD6GhVFAEELxOHrXfq/m3TEkSO8aVnMkBAIDlFBYWyt69e00KAIDdUdUcAABYTnp6ukyfPl1SUlIkOTk50NkBAPjJuskNHLFtCbwBAIDlNGjQwATdmgIAglesTdtse4vAGwAAWE5ERAQl3QCAoBHwwHvVqlWyYsUKCQ8Pl7POOks6d+5c7fwffPCBfPvtt+WmNW7cWG655RY/5xQAABwrWVlZsnbtWnNfEB8fz4YHANhaaCB77T777LPl9ttvl127dsnGjRvN53/84x/V/u7jjz+WefPmmV4+3a+oqKhjlm8AAOB/OrKH3htoCgCA3QVsODH9s1rS3atXr9JpCxYskEGDBpkL7amnnurxd2PHjpX9+/fL3Llza/23GU4MAABYGcOJAYD12WI4MR2bumzQrbp162bS3377rdrf6vePP/64TJ48WdavX+/XfAIAAAAAYOs23mXNnDnTVBvv2rVrtfOFhoZKZmamKRm/88475Y477pDx48dXOX9+fr55lX0y4a7uri8rcudNH1DoCwAApw0nNn/+fBk8eLAjezZ3l3hb9T4FACBenaMtE3h/88038sADD8ijjz4qjRo1qnI+DbRPPvnk0s/XXnutDBgwQPr37y/nnnuux9888cQT8sgjj1SanpaWZtm2Y7oTCwsLzXsCbwCA0+Tk5EiLFi3k8OHDJgB1Gvc6a+/uWuAAALBmR6CWb+Nd1po1a+T888+XG264QZ599lmvf9+0aVMZM2aMPPTQQzUu8W7WrJlkZGQcsS5+IANvfShAiTcAAM5DG28AsD6NKxMTE2vUxjvgJd46VEjfvn0lJSWlVkG3O0jNzc2t8nutvu6p53N9gmzlp8iaNwJvAIATFRUVmRsZ7bRGhxx1auBt9XsVAHCyUC/OzwE9k69bt84E3ddff71MnDjR4zw6dNiLL75YehH+7rvvKn3/xx9/mBJzAAAQHHQEkylTppgUAAC7C1hVc22z1bJlS/OUYPTo0eW+u/TSS0s7WBs1apSsXLlSNmzYIMXFxaYdd2xsrLRr10527txphiA7UudqFTGcGAAA1lZQUCCpqamm35fIyEhxGqqaA4D1eRNXBqzulgbct912m8fvwsLCygXh3bt3L52+fPlyWbZsmami3rlzZxNwn3TSSccs3wAAwP802NY+XAAACAaW6FztWKPEGwAAa8vOzpYff/xROnToIHFxceI0lHgDQHDFlfTWAQAALDmcmPbroikAAHbnvG5CAQCA5Wnb7ltvvTXQ2QAAwCco8QYAAAAAwI8IvAEAgOWkp6fL9OnTTQoAgN0ReAMAAMuJiIiQxo0bmxQAALujjTcAALAc7R22X79+gc4GAAA+QeANAAAsp7i4WA4fPix16tSRsLCwQGcHgB/l5Hs3unFsVIjf8gL4C4E3AACwnLS0NNPGOyUlRZKTkwOdHQB+1Gmcd305bJmS5Le8AP5CG28AAGA5iYmJMnToUJMCAGB3lHgDAADLiYqKkhNPPDHQ2QBwDKyb3KDc59x8l/S884B5/82z9SWGquUIAgTeAADAcrR9988//yxt27Y17bwBBK/q2mxr0E2bbgQDqpoDAADLyc7OlmXLlpkUAAC7o8QbAABYjo7hfeeddwY6GwAA+AQl3gAAAAAA+BGBNwAAsJwDBw7IrFmzTAoAgN0ReAMAAMsJDQ2V+Ph4kwIAYHe08QYAAJZTr149GTRoUKCzAQCAT/AYGQAAWE5JSYnk5uaaFIBz5Be6ZMGq/NLPd72eJfO+yTPTATsj8AYAAJaTmpoqL7zwgkkBOMOSdfnS664D8uCbfw0juHR9gfzjP9lm+ufr/grIAbsh8AYAAJasaj5kyBCTAnBG0P23yVmSlVO+ZLvkfx91+s2Ts8x8gB0ReAMAAMuJjo6WU045xaQAgptWI793araIy/znkZnuEjMf1c5hRwTeAADAcnJycmT9+vUmBRDcFq7Ol0M5riqDbjf9XudbtJpSb9gPgTcAALCcQ4cOyaJFi0wKILgtXlsgoSE1m1fn+2xtgb+zBPgcw4kBAADLSU5OlnvvvTfQ2QBwDGRml5S25T4SnS/zMKMdwH4o8QYAAAAQMPXiQr0q8a5XhxAG9sNRCwAALCcjI0PmzJljUgDBrW/nSK9KvC/oHOnvLAE+R+ANAAAsJyQkRMLDw00KILj1Pz1KEmJD5Ej/2vV7na/f6VHHKGeA7xB4AwAAy9Hxuy+99FLG8QYcICoiRJ4aGWci66qCbzM9RMx8Oj9gNwTeAADAclwulxQVFZkUQPDr0ylKXhoXL/Gx5YNqd9tvnf7yuHgzH2BHBN4AAMBy9u3bJxMmTDApAGc4v1OUfD2hvjyWElc6rXfHSHn6xjgznaAbdkbgDQAALKdu3boycOBAkwJwDq1GPqD7X6XaE0bHy+Ce0VQvh+0xjjcAALCcmJgYadeuXaCzAQCAT1DiDQAALCc3N1d++uknkwIAYHcE3gAAwHIOHjwoH374oUkBALA7qpoDAADLadSokdxxxx1mLG8AAOyOqxkAALCc0NBQiYyMDHQ2AADwCaqaAwAAy8nMzJR58+aZFEBwy8l3lXvl5rtKv9P3Fb8H7IgSbwAAYDklJSWSn59vUgDBrdO49Cq/63nngUrTtkxJ8nOOAN8j8AYAAJZTv359ueqqqwKdDQAAfILAGwAAAEDArJvcgK2PoEcbbwAAYDl79+6Vp556yqQAgltsVIhXL8COCLwBAIDlJCQkyIUXXmhSAADsjqrmAADAcmJjY6VTp06BzgYAAD5BiTcAALCcvLw8+eWXX0wKAIDdEXgDAADL0fG733vvPcbxBgAEBaqaAwAAy2nYsKHccsstEhUVFeisAABw1Ai8AQCA5YSFhZl23gh+Ofkur+anV2sAdkTgDQAALOfgwYOyfPlyOeuss6Ru3bqBzg78qNO4dK/m3zIlyW95AQB/oY03AACwnKKiIsnIyDApAAB2F+Jyubyr3xMEDh06ZJ6e69N0q44PWlJSYnpyDQkJMS8AAOAcenumr+joaAkNDXVUVfPcfJf0vPOAef/Ns/UlJqr8fRBVzQHYMa6kqjkAAAACprpAWoNuAm0AwSC4H6ECAABb2rdvnzz33HMmBQDA7gi8AQCA5cTFxUmvXr1MCgCA3VHVHAAAWE6dOnWke/fugc4GAAA+QYk3AACwnPz8fNmxY4dJAQCwOwJvAABgOTqU2OzZs00KAIDdUdUcAABYTlJSkowdO9ZUOQcAwO4o8QYAAJYTHh5uxkbVFM6RX+iSBav+al5w1+tZMu+bPDMdAOyMwBsAAFjOoUOH5NNPPzUpnGHJunzpddcBefDN7NJpS9cXyD/+k22mf76O9v4A7IvAGwAAWE5BQYH8/vvvJoUzgu6/Tc6SrJzyJdsl//uo02+enGXmAwA7CnG5XI6ru6NPz7X62sGDByUhIUGsqKSkRPLy8iQkJMS8AACAc+jtmb6io6MlNDS4y0m0GrmWaGtwXd1Nqd4NxceGyNcT6ktUBPdGAOwVVwb3mRwAAACWtnB1vhw6QtCt9Hudb9FqSr0B2A+BNwAAsJzU1FSZNGmSSRHcFq8tkNAaFmDrfJ+tpfkBAPsh8AYAAJYTGxsrXbp0MSmCW2Z2SWlb7iPR+TIPl/g7SwDgcwEfo2P9+vWyYsUKM1zImWeeKe3atTvib7KysmTevHmyb98+6dChg1x00UXHJK8AAODYiIuLM/cFCH714kJNSXZNgm+dr14dyo0A2E/AzlzaYcgFF1wgI0aMkB9//FG+/vpr6datmzz00EPV/m737t0m2H7xxRdl69atcsMNN8jQoUPN8gAAQHDQ3sz37NlDr+YO0LdzpFcl3hd0jvR3lgAgeHo11167P//8c+nbt2/ptPfff18uu+wy+fnnn6VNmzYef3fttdfKli1bTCl5RESEbN682ZSSv/POO3LFFVfU6G/TqzkAANa2d+9emT59uqSkpEhycrI4Db2aV0av5gCsxha9muvQGGWDbuWuUrZ9+3aPvykuLjZVzK+//noTdKtTTjlFzj77bJk7d+4xyDUAADgWGjRoICNHjjQpgpsODfbUyDgTWVfVx5qZHiJmPoYSA2BHAW/jXda7775rAmrtTMWTnTt3Sk5Ojpx88snlpuvnVatWVbnc/Px88yr7ZMJd6q4vK3LnjXG8AQBOFBYWVhp0W/VafSxKvJ2y7uedFiGTb46T+6YdNkOGubnbfuv43eNvqGPmc8o2AWB93pyPLBN4r169Wu699155+OGHpXHjxh7nyc7ONqkW55dVr1690u88eeKJJ+SRRx6pND0tLU3y8vLEqjuxsLDQvNfgGwAAJ9EH7dqXS6tWrRzZs7m7JaAWSGgtQSdof7zIf+8R+Wx9iDw1L8xM63lKifRu75Jz27kkKiJTGF0OgJVop9+2Cry1c7V+/fqZjtIeeOCBKuerU6dOuRJrN61T7/7Ok/vvv1/uvPPO0s/6+2bNmknDhg2PWBc/kIG3PhSgxBsA4ET79++Xbdu2mQ5Vk5KSxGmc1Ma7oquTXPLUvAzz/sVxDSQ2igIIANak52jbBN4bNmyQPn36mJ7JJ02aVO28zZs3l5iYGPn111/lwgsvLJ3+yy+/mLbeVYmKijKvivRCZuWLmeaNwBsA4ESNGjWS0aNHi1O5A2+r36v4Q2homarmZv0JvAFYkzfn54CeyTdu3GiCbu2N/KWXXvJYpfqjjz6S1157zbzXsb4HDRokM2bMkKKiIjNNq6EtW7bM9IYOAAAAAIDVhAey7db5559vgu3jjz9eHnvssdLvBg4cKJ06dSodYmzlypUyZswY8/nJJ5+UXr16Se/evc2439qbuVZTv/LKKwO1KgAAwMe0H5b33nvPPFjXpmEAANhZQKuajxo1yqRlexx3DxtWNgjv2LFj6eeWLVua6ulz5syRffv2yQsvvCCXXHKJ46phAQAQ7O3m2rRp41X7OQAArCrE5e4200G8Geg8UOhcDQAA53JS52o5+eVvRXPzXdLzzgPm/TfP1peYCp2r0dkaADvGlQHvXA0AAKAiHVIzIyNDEhMTzZBaCF6dxqVX+Z07AC9ryxTn9XIPwP6C+xEqAACwpfT0dHnjjTdMCgCA3VHiDQAALKd+/fpy/fXXmxTBbd3kBoHOAgD4HYE3AACwnMjISDnuuOMCnQ0cA7TZBuAEVDUHAACWk52dLcuXLzcpAAB2R+ANAAAsJzc3V9avX29SAADsjqrmAADAcho2bCjjxo0LdDYAAPAJSrwBAAAAAPAjAm8AAGA5+/fvl6lTp5oUAAC7I/AGAACW7NW8efPmJgUAwO5o4w0AACwnISFB+vbtG+hsAADgE5R4AwAAyykqKpKMjAyTAgBgdwTeAADAcrRt92uvvUYbbwBAUCDwBgAAlpOYmCjDhg0zKQAAdkcbbwAAYDlRUVGmczUAAIIBJd4AAMByDh8+LKtWrTIpAAB2R+ANAAAsRwPulStXEngDAIICVc0BAIDlNGrUSG677bZAZwMAAJ+gxBsAAAAAAD8i8AYAAJaTnp4uM2bMMCkAAHZH4A0AACwnPDxckpKSTAoAgN1xNQMAAJZTt25d6d+/f6CzAQCAT1DiDQAALKe4uFiys7NNCgCA3RF4AwAAy0lLS5PJkyebFAAAuyPwBgAAlpOYmChXXHGFSQEAsDvaeAMAAMuJioqSVq1aBTobAAD4BCXeAADAcnJycmTNmjUmBQDA7gi8AQCA5WRlZcnnn39uUgAA7I6q5gAAwHIaN24sd999d6CzAQCAT1DiDQAAAACAHxF4AwAAyzlw4IC88847JgUAwO4IvAEAgOWEhoZKbGysSQEAsLtaXc0yMzNr9R0AAEBN1KtXTy655BKTAgDgyMA7MTGxVt8BAADURElJieTn55sUAAC782n9rby8PImOjvblIgEAgAOlpqbKxIkTTQoAgKOGExs/frzH90qfSK9evVrat2/vu9wBAABH0irmgwcPpqo5AMB5gffcuXM9vlcRERHSsmVLmTp1qu9yBwAAHElr0LVp0ybQ2QAA4NgH3lqirfr16yeLFi3yTQ4AAAAqyM3NlV9//VVat24tMTExbB8AgPPaeBN0AwAAfzp48KB8/PHHJgUAwJGB9/r16+Xuu++uNF2n6XcAAABHo3HjxnLPPfeYFAAARwbet956q+nwpCIdb/P222/3Rb4AAICDhYSESGhoqEkBAHBk4P3tt99K586dK03XafodAADA0cjIyJD//ve/JgUAwJGBt1b7+u677ypN16A7KSnJF/kCAAAAAMC5gff1118vN9xwg3zwwQeSmZlpnkbPnz/fTNPvAAAAjkZiYqJcfvnlJgUAwFHDibn961//kj179sill14qJSUlZpq2wxoxYoQ89NBDvs4jAABwGJfLZV7axpt23gAAuwtx6VWtljT4XrdunbkgduzYUZo0aSJ2cOjQIalbt64ZoiQhIUGsSB9o5OXlccMBAHCkvXv3yvTp0yUlJUWSk5PFqQ8eoqOjTeEGAMDecWWtSrzdNNC2S7ANAADsQ29kLr74YpMCAGB3tX6EunTpUhk+fLj06tWrdNrkyZNNtA8AAHA0YmJipEOHDiYFAMCRgfecOXNk0KBBEh8fLytWrCidnpOTI08//bQv8wcAABxIm1tt2rTJpAAAODLwfvTRR2X27Nny0ksvlZuuna29+eabvsobAABwKB01RUdM0RQAAEd2rqbVvg4cOGBS7fDD3bN5VlaWNGjQQAoKCsTK6FwNAABr03uLwsJCiYiIcGTnYnSuBgDW501cWasrWVJSkmzdutW8LzvExxdffCEtWrSozSIBAABKabAdFRXlyKAbABB8anU10/G6b7rpJtmwYYMJvNPT02XmzJkyatQoufHGG32fSwAA4ChaxfyDDz6gqjkAICjUajixhx9+WFJTU83Y3VoVTEvA9Yn06NGj5Z577vF9LgEAgKPo/YV22upuzgYAgOPaeLvt3btX1q5day6KnTp1ss2Y3rTxBgAAVkYbbwBwaBtvDazL9mqukpOTpX///jJgwADbBN0AAAAAABxLNQ68f/75Z9O7qPrXv/7lzzwBAACH27dvnzzzzDMmBQDAMW2827dvL2PGjJHu3bubz6+88kqV844dO9Y3uQMAAI4UHx8vffr0MSkAAI5p471+/Xp54IEH5Ndff5XNmzdLq1atqpxX57Ey2ngDAAAro403AFifN3FljUu8GzRoIAsWLDDvdQgxqwfXAADAvvLz82X37t3StGlTM543AACOaOPdrFkz/+YEAADgfzIyMmTu3LkmBQDAMYG3trFKT0/3b24AAABEpGHDhjJu3DiTAgBgdzWuan7BBRdIt27dpG3btubzwIEDq5zXXSUdAACgNsLCwiQuLo6NBwBwVuA9Y8YMefPNN03b7o8//lhat27t35wBAADH0o5qVqxYIWeeeabpuAYAAEcE3rGxsaXDhK1cuVImTpzokwwUFxfLwoULZdOmTXLVVVcdsS354sWLZd26deWmJSUlyYgRI3ySHwAAEHhFRUWyf/9+kwIA4JjAu6zly5f75I9rpyn33HOPCba/+uorOf30048YeOtvli5dKoMGDSqdpr2sAwCA4KGjqQwfPjzQ2QAA4NgH3n379jUlzm4PPvigPProo+Xm0SroNR1qTDtM+eKLL0w7Lm96Te/QoYM888wzXuQcAAAAAACL92qulixZUu7zY489VmmerVu31nh55557rrRo0UK89fvvv8tLL70kM2fO9OrvAQAAe0hNTZXnn3/epAAAOLKqeaDpmJ7aznvPnj0ycuRI+fe//y33339/lfPn5+ebl9uhQ4dMWlJSYl5W5M6bVqOnKj0AwGm0b5kePXqY1KrXan9yuVzm5cR1BwC78OYcbbvA++abb5aXX365NBidPXu2DBs2THr37i1nnHGGx9888cQT8sgjj1SanpaWJnl5eWLVnVhYWGjeE3gDAJzoxBNPlJycHPNyGg26VUREhISGelVBEQBwjGRlZdV43hCX+8xek5lDQkovBJ4+VzXtSHbv3m3aeGunaeedd55467jjjpO///3v8sADD9S4xFv/npacJyQkiFUDb30oQIk3AMCJ9Lqt1cwbNWokUVFR4tQS7+joaAJvALAojSsTExPNEJhHiiu9LvG++uqrq/0cCOHh4aXVxz3RC7ani7Y+QbbyU2TNG4E3AMCJ9CbmnXfekZSUFElOThanBt5Wv1cBACcL9eL87FXgfdFFF0lmZmaVn93TfOnTTz81nanpON065vcvv/wibdq0Kf1ee1nXEvNzzjnHp38XAAAETlJSkowZM0bi4+PZDQAA2/Mq8F60aJFP//jGjRtl4cKF5qm2u7326tWr5cwzzzQv9e6778rKlStN4K1Pfq+55hozZFm7du1k586dMmvWLLnppptkwIABPs0bAAAIbG02rb4HAEAwCGjnarm5ubJ3717z/q677jKpfs7Ozi5Xgt62bdvSi7AG5gsWLJC1a9dKx44dZdy4cdKlS5cArQEAAPAHbUL27bffSvfu3S3bHwsAAH7pXC2YLuZ169atUSP4QKFzNQCAk+3fv18++OADueSSS0y1c6ehczUACK640nbDiQEAgOCnwfbIkSMDnQ0AAHyiVt1kVuxQrabfAQAAAADgNLUKvKvr7ISOUAAAwNFKS0uTyZMnmxQAALvz6cCQeXl5Eh0d7ctFAgAAB4qJiTGdqGoKAIDdedXGe/z48R7fuzsD0x7H27dv77vcAQAAR4qLi5Ozzjor0NkAAODYB95z5871+F5FRERIy5YtZerUqewaAABwVAoKCiQ9PV0aNGggkZGRbE0AgHMCby3RVv369ZNFixb5K08AAMDhDhw4IG+++aakpKRIcnJyoLMDAMBRYRxvxvEGAMByCgsLJSMjw3TaqrXqnIZxvAHA+o7JON4fffSRfP311+aJdEWvvPJKbRcLAABggu1GjRqxJQAAQaFWgfeDDz4oTz31lJx99tkMHwYAAHwuKytL1qxZI126dJH4+Hi2MADAeYH3lClT5JNPPpHevXv7PkcAAMDxdIjSTZs2yamnnkrgDQBwZuBdVFQkPXr08H1uAAAARKRhw4Zy0003sS0AAEEhtDY/0irmn332me9zAwAAAABAkKlViXerVq3kmmuukREjRkjr1q0lJCSk3Pe33367r/IHAAAcKC0tTebNmydDhgwxpd8AADhuOLH27dtX+/2GDRskWLp9D5SSkhLTvk0falR8sAEAQLDTa/X3338vXbt2tey12p8YTgwArM/vw4lZPbAGAAD2pjcwdOIKAHB0G28AAAB/KiwsNNXNNQUAwLGB99KlS2X48OHSq1ev0mmTJ082xewAAABHIz09XaZOnWpSAAAcGXjPmTNHBg0aZMbVXLFiRen0nJwcefrpp32ZPwAA4ED169eX6667zqQAADiyc7WOHTvK448/LgMGDDAdf7kX8euvv0qfPn1k586dYmV0rgYAAKyMztUAwPq8iStrVeK9ZcsWE2Crsj1uN27cWPbu3VubRQIAAJTKzs42teo0BQDA7moVeCclJcnWrVsrBd5ffPGFtGjRwne5AwAAjqTN19asWWNSAAAcGXiPGDFCbrrpJjOsmAbe2vHJzJkzZdSoUXLjjTf6PpcAAMBRGjVqJH//+99NCgCA3dVqHO+HH35YUlNTTVvvkpISUwIeGhoqo0ePlnvuucf3uQQAAAAAwEmdq7lpe+61a9ea4LtTp07SpEkTsQM6VwMAwNr2798vCxYskIEDB5oH/E5D52oAEFxxZa1KvN2Sk5Olf//+R7MIAACASiIjI+X44483KQAAdhd+NKXd33//vWRkZFT6TsfdBAAAqC0tObjwwgvZgAAA5wbes2bNMp2oacdqnorUCbwBAMDRKCoqksOHD0udOnUkPPyoKugBAGDPXs3vv/9+mTBhgrkgasl3xRcAAMDRtvF+5ZVXTAoAgCMD77S0NElJSSk3hjcAAICvJCYmylVXXWVSAAAcGXj36tVLli9f7vvcAAAAiEhUVJS0bNnSpAAA2F2tGk1NnjzZdHhyxRVXSKtWrSqVfI8dO9ZX+QMAAA6kzdk2btwo7dq1M+28AQBwXOA9e/Zs+e2332TatGlSr169St8TeAMAgKORnZ0tX3/9tbRo0YLAGwBgeyEul8vl7Y+SkpJM52razjvYBzoPlJKSEsnLyzO1CWhLDwCAs+jtmb6io6MlNLRWLQMBABaKK2t1Ji8uLjbVzAEAAAAAgPg+8O7WrZssWbKkNj8FAAA4ovT0dHnrrbdMCgCAI9t4t27dWoYNGyYjRoww7ytWhb799tt9lT8AAOBA4eHhZigxTQEAcGQb7/bt21f7/YYNG8TKaOMNAACsjDbeAGB93sSVtXqMrGN4e+rNXGVmZtZmkQAAAOX6k8nPzzfjeIeFhbFlAADOa+OtVb9q8x0AAEBNpKWlyYsvvmhSAADszqfjU+jwVzrsBQAAwNHQmnWXXXZZlTXsAACwE6+qmo8fP97je/e406tXrz5i+28AAIAj0Qf5J510EhsKAOC8wHvu3Lke36uIiAhp2bKlTJ061Xe5AwAAjpSTkyNbtmyRk08+WWJjYwOdHQAAjl3grSXaql+/frJo0aKj+8sAAADV9BT76aefSnJyMoE3AMCZw4nZHcOJAQAAK2M4MQCwPr8PJ+aWmpoqO3fulKKionLTzzjjjKNZLAAAAAAAQaNWgfeePXvkmmuukWXLlnn83oGF6AAAwIcOHDggn332mVxwwQVSv359ti0AwHnDid1+++1y3HHHya5du8xnHWPzo48+klatWsmkSZN8nUcAAOAwoaGhEhUVZVIAABzZxrtRo0ayZs0aadq0qYSEhJiq5mFhYbJq1SpJSUmRTZs2iZXRxhsAAFgZbbwBwPq8iStr9RhZS7g16FZa/UvbeqsOHTrI9u3ba7NIAACAUiUlJVJQUGBSAADs7qjrb3Xq1Elefvllyc7OlldffVWaN2/um5wBAADH0of6zz33XOnDfQAAHNe52uDBg0vfP/roozJgwAD5f//v/0l0dLTMmDHDl/kDAAAOpFX3Bg0aZFIAAOzOJ+N4a2m3tutu2bKlJCUlidXRxhsAAFgZbbwBwPr83sa7ori4ODn99NNtEXQDAADry83NlY0bN5oUAAC7q3XgvXTpUhk+fLj06tWrdNrkyZNNtA8AAHA09H5iwYIF3FcAAJwbeM+ZM8e0u4qPj5cVK1aUTs/JyZGnn37al/kDAAAO1LhxY7nrrrtMCgCAIwNv7VBt9uzZ8tJLL5Wbfumll8qbb77pq7wBAACHCgkJkfDwcJMCAODIwHvLli3Sp08f877sBVGfSu/du9d3uQMAAI6UmZkp77//vkkBAHBk4K2dqG3durVS4P3FF19IixYtfJc7AADg2F69i4qKTAoAgCMD7xEjRshNN90kGzZsMIF3enq6zJw5U0aNGiU33nij73MJAAAcJTExUYYOHWpSAADsLrw2P3r44YclNTVVOnbsKCUlJaYEPDQ0VEaPHi333HOP73MJAAAAAIBNhbiOog6Xtudeu3atCb47deokTZo0kWAb6DxQdJvm5eWZGgV0LAMAcBq9x5g+fbqkpKRIcnKyOI3enukrOjraFG4AAOwdV9aqxDsmJkZyc3PNhbB///61zScAAIBHegPTr18/yz4gBwDAG7V6hKrjd+/fv782PwUAADii2NhY06RNUwAAHBl433DDDfLII49IYWHhUWdAq1EtXrxYJk2aJHv27KnRb/Lz8+XDDz+UKVOmyKpVq446DwAAwFq0udXmzZtNCgCA3dWqqvny5ctlxYoV8s4778hJJ50kkZGR5b7XYcVqYv78+aYzNu2x9Ntvv5X27dsfsZ24dup23nnnmff6JPwf//iHXHbZZSYIBwAAwUHH7543b55j23gDABwceLdu3Vp+/fVX6d27t3kdrTp16shHH31k2ow3a9asRr+57777JCIiQlauXGl+t27dOunatasMHjxYBg0adNR5AgAAgdeoUSO59dZbJSoqKtBZAQDg2AbeW7duNemjjz569H9ZRPr27WvS3bt317in77lz55pq7hp0K+1N/cwzz5TZs2cTeAMAECS0J2/3tR4AAEdWNQ+UXbt2SVZWlrRt27bcdP28evXqatuE66tst+/uQF5fVuTOG8OJAQCcSIdm0aZtZ511lhmqxanDiVn1PgUAIF6do70OvPUieCR6kfQHDbpVvXr1yk3XNuLuYNqTJ554wpSSV5SWlmbZTlt0J7o7r2McbwCA02RnZ5trdEZGhk86c7UbDbqVNq9jHG8AsCZ3fOqXwPvss8+u8cXC19xVziquoAbd1Q03cv/998udd95Zbn5tU96wYUPLjg+qgbfecFDiDQBwoqSkJGnZsqU4lbvEOzo6msAbACxKz9F+C7x//vlnCZTmzZubHtS3b99ebvq2bdtM7+pV0Y5ZPHXOok+QrfwUWfNG4A0AgHMDb6vfqwCAk4V6cX72OvBu06aNHEtLly6VvXv3yrBhw0x1q/79+8usWbNk9OjRJijVjtl0+LLXXnvtmOYLAAD4z759+2TmzJly7bXXSuPGjdnUAABbC2jnaps2bZLFixebsTrV+++/Lxs2bJDu3bubl9KLrg4dpoG3evLJJ00v5gMHDpQePXrIW2+9ZT7rhRkAAASHuLg4Oeecc0wKAICjAm8NdH3dY6kG32rcuHFSXFxsPp944oml8/Tp06fc51NOOcUE5zNmzDBPw//5z3+aoDs83FYdtAMAgGrUqVNHTj/9dLYRACAohLj81ROahWnnajo0iQb+dK4GAID16DCge/bskSZNmnjspyXY0bkaAARXXElvHQAAwHJ0GLE5c+aYFAAAu6N+NgAAsBwd8vPmm282Vc4BALA7Am8AAGA5YWFhlm0OBgCAt6hqDgAALNlubtGiRSYFAMDuCLwBAIDlFBYWmtFLNAUAwO6oag4AACynQYMGkpKSEuhsAADgE5R4AwAAAADgRwTeAADAclJTU+WFF14wKQAAdkfgDQAALCc2Nla6detmUgAA7I423gAAwHLi4uKkZ8+egc4GAAA+QYk3AACwnIKCAtm9e7dJAQCwOwJvAABgOQcOHJCZM2eaFAAAu6OqOQAAsJykpCQZNWqU1K1bN9BZAQDgqBF4AwAAywkPDzdjeQMAEAyoag4AACzn0KFDsmTJEpMCAGB3BN4AAMBytFO17du307kaACAoUNUcAABYto03AADBgBJvAAAAAAD8iMAbAABYTlpamrz88ssmBQDA7gi8AQCA5URHR0u7du1MCgCA3dHGGwAAWE58fLycc845gc4GAAA+QYk3AACwnMLCQtm7d69JAQCwOwJvAABgOenp6TJ9+nSTAgBgdwTeAADAcho0aCApKSkmBQDA7mjjDQAALCciIkKSk5MDnQ0AAHyCEm8AAGA5WVlZsmzZMpMCAGB3BN4AAMBy8vLyZOPGjSYFAMDuqGoOAAAsp2HDhnLzzTcHOhsAAPgEJd4AAAAAAPgRgTcAALCc/fv3y5QpU0wKAIDdEXgDAADLiYyMlBNOOMGkAADYHW28AQCA5SQkJMj5558f6GwAAOATlHgDAADLKSoqkvT0dJMCAGB3BN4AAMByaOMNAAgmBN4AAMBy6tevL9dee61JAQCwO9p4AwAAy9FO1Zo2bRrobAAA4BOUeAMAAMvJzs6Wb775xqQAANgdgTcAALCcnJwc+e6770wKAIDdUdUcAABYTqNGjeTWW28NdDYAAPAJSrwBAAAAAPAjAm8AAGA5Oob39OnTTQoAgN0ReAMAAMuJiIiQxo0bmxQAALujjTcAALCchIQE6devX6CzAQCAT1DiDQAALKe4uFgOHTpkUgAA7I7AGwAAWE5aWpq8/PLLJgUAwO4IvAEAgOUkJibK0KFDTQoAgN3RxhsAAFhOVFSUnHjiiYHOBgAAPkGJNwAAsJzDhw/L6tWrTQoAgN0ReAMAAMvJzs6WZcuWmRQAALujqjkAALAcHcP7zjvvDHQ2AADwCUq8AQAAAADwIwJvAABgOQcOHJBZs2aZFAAAuyPwBgAAlhMaGirx8fEmBQDA7mjjDQAALKdevXoyaNCgQGcDAACf4DEyAACwnJKSEsnNzTUpAAB2R+ANAAAsJzU1VV544QWTAgBgdwTeAADAklXNhwwZYlIAAOyONt4AAMByoqOj5ZRTTgl0NgAA8AlKvAEAgOXk5OTI+vXrTQoAgN0ReAMAAMs5dOiQLFq0yKQAANhdiMvlconD6EW8bt26cvDgQUlISBAr0l5c8/LyJCQkxLwAAIBz6O2ZvrTKPWOZA4D940pKvAEAAAAA8CMCbwAAYDkZGRkyZ84ckwIAYHcE3gAAwHK0mVV4eDjNrQAAQSHgw4kVFRXJV199Jfv27ZMOHTpIu3btqp3/m2++kc2bN5eblpiYKIMHD/ZzTgEAwLGi43dfeumlbHAAQFAIaOB94MABueCCC0zavn17+fLLL2XkyJEyceLEKn8zffp008vpeeedVzqtadOmBN4AAAQR7VisuLhYwsLCKPUGANheQAPvf/7zn5Kbmys//vijxMXFyapVq6Rnz57Sv39/ueiii6r83emnny7Tpk07pnkFAADHjtaE04ftKSkpkpyczKYHANhaaCCfZL/zzjumhFuDbtWjRw/zmjVrVrW/TU9Pl3fffVc++eQTSUtLO0Y5BgAAx4oOzzJw4ECTAgBgdwEr8d61a5cZ76xim26tcr5mzZpqf6ttvDVo37Nnj2zcuFEmTJggN910U5Xz5+fnm1fZ8dbcY2Xry4rceWMcbwCAE0VFRUnbtm3Ne6teq4/FON5OXHcAsAtvztEBC7w16HZ3nlJW/fr1S7/zZPjw4fLCCy9IZGSk+fzyyy/LuHHjTEl5p06dPP7miSeekEceeaTSdC0tz8vLE6vuxMLCQvNeg28AAJykoKDAVDdv3Lhx6TXfSTToVhERERIayiA0AGBFWVlZ1g+8Y2JiTHr48OFKmXd/50mvXr3Kfb755pvloYceMh2uVRV433///XLnnXeWK/Fu1qyZNGzYUBISEsSqgbc+FKDEGwDgRBp0L168WK6//npJSkoSp5Z4R0dHE3gDgEXpOdrygXfz5s3NU9wdO3aUm66fW7Vq5dWyYmNjTbvv6qqr6asifYJs5afImjcCbwCAE2lJ9x133GHG8rbytdrfgbfV71UAwMlCvTg/B+xMrtXGdCix2bNnl3u6/fnnn5vOVNxWrFgh8+fPLy0F1nbdZWlP6Npe/IwzzjiGuQcAAP6+mdF7BYJOAEAwCHG5GxEFwIYNG0zVcR06TAPnN954Q+rUqSNfffWVKQ1Xo0aNkpUrV5p5tc1zx44dzRje2inbzp075ZVXXpELL7zQBPA1vThrVXPtJVXbklPVHAAA68nMzJQvvvjCXPMr9gfjBFQ1BwDr8yauDGjdJe3BfP369dKmTRvTU7n2TK4XWXfQrTQwHzJkiHmv07XHcx3He9OmTab6+Ny5c2XOnDk8EQcAIIhoLTcdkYRevQEAwSCgJd6BQok3AACwMkq8AcD6bFPiDQAAAABAsCPwBgAAlrN371556qmnTAoAgN0ReAMAAMvRKnvaeapVO0EFAMAbARvHGwAAoCqxsbHSqVMnNhAAIChQ4g0AACwnLy9PfvnlF5MCAGB3BN4AAMCS43i/9957JgUAwO6oag4AACynYcOGcsstt0hUVFSgswIAwFEj8AYAAJYTFhZm2nkDABAMqGoOAAAs5+DBg/LRRx+ZFAAAuyPwBgAAllNUVCQZGRkmBQDA7qhqDgAALKdBgwZy3XXXBTobAAD4BCXeAAAAAAD4EYE3AACwnH379slzzz1nUgAA7I7AGwAAWE5cXJz06tXLpAAA2B1tvAEAgOXUqVNHunfvHuhsAADgE5R4AwAAy8nPz5cdO3aYFAAAuyPwBgAAlqNDic2ePdukAADYHVXNAQCA5SQlJcnYsWNNlXMAAOyOwBsAAFhOeHi41K1bN9DZAADAJ6hqDgAALOfQoUPy6aefmhQAALsj8AYAAJZTUFAgv//+u0kBALA7qpoDAABLtvEeMWJEoLMBAIBPUOINAAAAAIAfEXgDAADLSU1NlUmTJpkUAAC7I/AGAACWExsbK126dDEpAAB2RxtvAABgOXFxcXLmmWcGOhsAAPgEJd4AAMBytDfzPXv20Ks5ACAoEHgDAADLOXDggLz11lsmBQDA7qhqDgAALKdBgwYycuRIqVevXqCzAgDAUSPwBgAAlhMRESENGzYMdDYAAPAJqpoDAADLOXTokCxdutSkAADYHYE3AACwnPz8fPn1119NCgCA3VHVHAAAWI5WMx89enSgswEAgE9Q4g0AAAAAgB8ReAMAAMtJS0uTV1991aQAANgdgTcAALCc6OhoadOmjUkBALA72ngDAADLiY+Pl3PPPTfQ2QAAwCco8QYAAJZTWFgoqampJgUAwO4IvAEAgOWkp6fLG2+8YVIAAOyOwBsAAFhO/fr15frrrzcpAAB2RxtvAABgOZGRkXLccccFOhsAAPgEJd4AAMBysrOzZfny5SYFAMDuCLwBAIDl5Obmyvr1600KAIDdUdUcABBQOfkur+aPjQrxW15gHQ0bNpRx48YFOhsAAPgEgTcAIKDOuOuwV/P/MCnOb3kBAADwB6qaAwAAy9m/f79MnTrVpAAA2B0l3gCAgFo5oU65z7n5Lun9zxzzfunjsRJD1XLH9mrevHlzkwIAYHcE3gCAgKquzbYG3bTpdqaEhATp27dvoLMBAIBPEHjDEuhcCQBQVlFRkWRlZUl8fLyEh3O7AgCwN65ksAQ6VwIAlKVtu6dPny4pKSmSnJzMxgEA2BqdqwEAAMtJTEyUYcOGmRQAALujxBuWQOdKAICyoqKiTOdqAAAEAwJvWAKdKzkbbfwBVHT48GHZsGGDtG/fXurUKd/zPQAAdkPgDSDgaOMPwFPgvXLlSjnhhBMIvAEAtkcbbwCAZeQXumTh6sLSz/dNy5MPvy000+EsjRo1kttuu82kAADYHSXeAAKONv5QS38okn/NyJNDuX9tj2UbiuWLH4vlyTn58uj10XJeBy5bAADAfijxBmCJNv5lXzFRIaXf6fuK3yM4g+7bX8+TrDJBtyr5X0G3Tr/ttTwzH5whPT1dZsyYYVIAAOyOwBsAEFBajVxLusVl/vPITHeJmY9q584QHh4uSUlJJgUAwO4IvAEAAfXp2iJTvfxIrbj1e53vs7WUejtB3bp1pX///iYFAMDuCLwBAAH1+foiCa1hCwKdb8l6Am8nKC4uluzsbJMCAGB3BN4AgIA6mOMqbct9JDqfzo/gl5aWJpMnTzYpAAB2R+ANAAiourEhXpV46/wIfomJiXLFFVeYFAAAuyPwBgAEVJ+O4V6VeJ/fkc62gt2e9BJZtDZUvtnZzKT6GQAAO+PuBYClaI/VC1cXln6+b1qeXNA5XC7sHC5REZR0BiPdtzpOtw4ZVl38rXs/PkbM8YDgtOG3Ynnl4wL56qdicZU5GEJCRM4+NUzGXhwp7VuEBTKLAADYs8S7pKREvvnmG5k3b578+uuvfvsNAOvTMZrP/+dheeTtgtJpyzYUywNv5pvpX/xIp1rBSB+oPHp9tImsq3q0YqaHiJmPBzDBafG6Ikl5LleWbSwfdCv9rNP1e50PAAC7CWjgffDgQenVq5dpw6UdqHTs2FHuu+8+n/8GgD2C7ttfzzOlnmW5qyDr9NteyzPzIfic1yFcJo6ONiXaZbnbfuv058dEm/kQnCXd907Lk8L//fM+rn6o/H1QjEwYHWfS5MQ/b1f0e51P5wcAwE4CegfzwAMPyIEDB+Snn34y43QuX75czj77bLngggvk/PPP99lvYC9UNXbmPv/XjDxTz7iqqsY6PcQlZr4lj9eh1DMI9T4t3OzbBd8VySOz8s20c9qHmaroWr2cku7gpdXL3UH3kJ5R8mhKnESG/1X/YeyAWHlgerbM/ybfzPfKwgKZNLbCUxoAACwsYCXeLpdLZs6cKTfeeKMJoNVZZ50l3bt3l7feestnv4G9UNXYmT5dWySHjtC+V+n3Ot9nayn1DlYaXPfv+tcz4fEjomVg9wiC7iCmHadpm253SXfFoFvp58dS4kpLvr/aWCy/H6DDNQCAfQSsxHv37t2SmZkp7du3Lze9Q4cOsm7dOp/9RuXn55uX26FDh0z6xx9/SHZ2dun06OhoqVevnhQVFcn+/fsrLSc5Odmk6enpUlj4V+dPSh8ExMTEyOHDhyUrK6vcd5GRkVK/fn3TNj01NbXSchs2bChhYWGSkZFRmk+dt6CgQOLj46VOnTpmuq57WeHh4dKgQQPzXperDybK0r8ZERFh1jc3t3z9XV1mXFycx+VqXpKSksx73Q7FxeWr9Ok2ioqKMttO17cs3QYJCQlm+2jNhLJCQkKkUaNGpdtQt3NZP+ypK/+YXlIp+qpY1fjhK3Kk58nFlbZhaGio2Ya63crS/Gi+8vLyTFOFsnT76HZS+/btK/+HRcz21e2sv9Pfl6Xbr6p9U3Yb6hi0uj/L0uFx9LjQYyUnJ6fcd7GxsWa/12YbuveN7peyx7b7+NbjVPenp+Nbl6vLr24b6nHk/vfjpuuh66PHn6fjW7eDbo/qtuFna/IlJMQlLteRO0/T+T5alS0XdIw1+6+6bajroetTlh4nerwc6fiubhvqdtftX1Hjxo1NqvvN0zlCf1/dNtTjxNOYxe5tqMdZ2XOZ0vXU9Q2mc4SrTBfn+/bulejI8svV/V3xPKvT9fuqjm/OEdY9R6zcpt/9+f7yXlGVgu7S5YSHyOVnRcnkD3PN/J9/nynXnF/viMe3Xc8R7vXR7afTvb2PcNN/x/rS5Va8Buo5wv1vee/evZWWq+cP3Yb6u4rnCN2Gut/173nahu5rlebX0zVQt6GuZ8V9o8eRbgvdPp62ofteTPdNxWug+15M90vFfaN/T/+u7k9P51nNr+Zb903F49t9L+ZpG+r2cZ9nPW1D3b66nfU8WvEaqMvUZXvahro/3cehp22ox4MeG95uQ/037D4OPW1DPb71ePO0Dd33yVVtQ12uLt/TNtRjRY8Zzaunc4Sujx7znu7F3Me3p23oPr51uqdroPv41uVWPEfU5PjW9ah4rSp7fOt2qHiOcB/feu6peA20UqxRcRtyjoiq1Tmi4r6wZODtPnG5g56y/wgq/sM5mt+oJ554Qh555JFK09977z3zD8CtRYsWcuaZZ5oNuGDBgkrzDxs2zKSfffZZpX8sPXv2lJYtW8qWLVvk+++/r/SPqHfv3uYf0Ny5cyst99JLLzX5WLZsmezZs6fSQ4WTTz7ZPHRYtWpVue/0H667ev38+fMrnZT79u1r/pFqfnbs2FHuu1NOOcU8wNADS/9uWZqXAQMGmPcfffRRpZPcOeecY/4Bb9iwQTZv3lzuO90GXbt2Nftq8eLF5b7Tg1XXVS1ZsqTcPisqCZWZW68WcYWKq4rulf6sauySJ/4bJtef8r6Eh/61vpdccok5eX711VeVTjidOnWSVq1ayc6dO+W7774r950eS7pvlHbWV9FFF11kTkbffvut7Nq1q9x3bdu2lVNPPdVcaL/++utKF9N+/fqZ93osVbwAnXfeeea4Xb9+faUOAk888UTp3LmzOTl+/vnnlS4igwcPNu8//fTTSv/Y9Tg8/vjjZdOmTbJx48Zy3zVp0kTOOOMMc9FbuHBhpXUdMmSIOSl/+eWXlY7vLl26yAknnCDbt2+XNWvWlPtOL2rnnnuuOVF52ob9+/c3F7CVK1dWOr7btWsnbdq0kb3pWqMlVmpCg/Ptu9Jlx45Uc4Jcu3atbNu2rdw8rVu3Nv0/6IXriy++qHRxGjRokHm/aNGiShdF7UNC/81qc5aff/653HfNmjUzNWz0ZuSTTz6plLfLL7/cpEuXLq10oe7WrZs0b95ctm7dWulBoZ7MtcmMniM++OCDSssdOHCguYivWLHCPDAM9nNErvnn8ufx8OGCBRIR+ueNYY8ePaRp06bmPPvjjz+W++1xxx1nzt96Q+Hp/M05wrrniMJGel2IM+9PSK6+x/ITGv/1/ffrNsmFHVv/eZx8+GGlAKJPnz5Bc444mvsIvZ7oeVavgRWvVbp93NeqOXPmVDpH6PlbzyN6bqm4DfUaqNdXDWgqXqs0ONDjRek+rxjU6L7RIE3Xs+I21Gug/lvXe4SKx6GeI6666qrSfVMxWNV9o9tQj2/d7xWPbz036fnq/fffr7QNte8gvY/QfVMxgNZzlp5n9TypnftWPL61uaOaPXu2x/O3Btd6/v7tt9/KfafnWD2H63m94nGo9x7u41DvVysGS3r+1vOs/nureJ7V41uPJz3GKh6Heh8xdOhQ8/7jjz+u9CBBr0V6ntXzww8//FDp+NaapnqO0GtKRVdeeaU5R+g9XsV7Mf13ofdienzrPVXFa6Beq/Qc8e6771Zart736H2ENi+teC922mmnmXsJvQbqPWBZen27+OKLzXv9d1PxHKH3eHofqPeGFe/F9Bqo5zW9Bla8Vun1+LLLLis991R8QKH3eHpN0uuUXgfLsmKswTni6M4RFe+BqhPiqvj45xj55ZdfzElMDyw9ebiNGzfOXNQrHqi1/U1VJd568tATs54MrfgUymkl3ot/DJenP6h5e71/XJIr53f46/eUZtm7xPv213Lky401L/E+8+QimTiGEu9gLfHOzXdJz3v+zM/8e7Io8Q7yWjErtiXKI2//mU/tSO3WwXWkKs/PP2xKvNU/BhdQ4v0/lGb9iRLvP1Hi/SdKvP/COcJ/5wiN+/Thpt7n6rXQkoG33szpTd2rr74qI0eOLPdkUFfKU6lPbX7jid4U6IavyQYKFA28NUjRmxx9Bbs7Xs+VpT8Ul1Yrr472ctz7tDB5bjQd6wSLD78tNEOG1dTj10eZdr8ITjn5Ljnjrj8D9pUT6khsVPCfA53exvvif+eY6uPaxvuzxxM9VjcvKHJJ3/szZG9GiRnXe+EjsXJ8/YCPiuo3enumL30govc4AADr8SauDNiZXJ8qaLUSrdrkpk/YtbqSuwqj0uofWp3Im9/Afg7muGoUdCudT+dH8NBeqxNiqh7D2U2/1/m0h2sAwaFJg1A5+9Q/q5D/caDE9F6uQXZZ+lmna9Ctzm4XFtRBNwAg+AT07nX8+PGmrcjw4cNNu4UpU6aYNkMjRowonUdLt7VtqLsNUk1+A/upGxtiSrJrWuKt8yO4erJ+9Ppo03meDhnm6TAwezxEzHwMKxV8JdxlaVVzT+/dKAEPPmMvjpRvNueaocJ0yLBVmwpNR2rapnv7vmL57/L80qA7IlxkbP//9bgHAIBNBKyquZt2EKDBs3bOoR1M3HzzzaZtodtrr71m2nY//fTTNf7NkVDV3Hqoagz3cHI6TrcOGebmfiCjJd0adJ/XgdLuYHPa38u3szqSHyb92REXgsvidUVy77S80vG8PdGg+8kR0dK3U/CfB6hqDgDW501cGfDAOxAIvK0nv9Al5//zsBkyrLoDUks942NEljxeh1LPID4WFnxXJI/M+rPN93kdwkxVdK1eTkl3cCLwhtuG34rllYUFZpzusncn2qZbq5drSXf7FtX3fB4sCLwBwPoIvH24gQLFaZ2rqS9+LDJVjeUIVY2fH0OpZ7Cjcy1nVzU/EqqaO6PDtdW/FElWrkviY0Kk28nhjmvTTeANAMEVVwZ/XS3YhlYhnjg6usqqxlrSTVVjIPgQSMNTh2tNGtCOGwAQPJz1+BiW1/u0cFON/OFrokqnndM+zAwfpdNp3wsAzqDjhf/3v/81KQAAdkfgDcvRdrz9u/5VGWP8iGgzZjPtewEAAADYEVXNAQCA5SQmJsrll18e6GwAAOATlHgDAABLdi6mHY06cPAVAEAQIvAGAACWs2/fPnn66adNCgCA3VHVHIDlhpPKLfO57Hs3esEGgp8Oz3LxxRebFAAAuyPwBhBwZ9x1uMrvev8zp9K0HybF+TlHAAItJiZGOnToEOhsAADgE1Q1BwAAlpOXlyebNm0yKQAAdkeJN4CAWzmhTqCzAMBiMjMzZf78+ZKSkiLJycmBzg4AAEeFwBtAwNFmG0BFjRo1kttvv10iIiLYOAAA2yPwBgAAlhMaGipRUVGBzgYAAD5BG28AAGDJquYffPCBSQEAsDsCbwAAYDklJSWSk5NjUgAA7I6q5gAAwHLq168vV199daCzAQCATxB4wxJy8l3lPueW+Vz2vRudcQEAAACwCwJvWMIZdx2u8rve/8ypNO2HSXF+zhEAIJD27dsnM2bMkOHDh0vjxo3ZGQAAW6ONNwAAsJz4+Hjp06ePSQEAsDtKvGEJKyfUCXQWAAAWEhsbK126dAl0NgAA8AkCb1gCbbYBAGXl5+fL7t27pWnTpoznDQCwPaqaAwAAy8nIyJC5c+eaFAAAu6PEGwAAWE7Dhg1l3LhxEhMTE+isAABw1Ai8AQCA5YSFhUlcHCNYAACCA1XNAQCA5Rw8eFAWLlxoUgAA7I7AGwAAWE5RUZHs37/fpAAA2B1VzQEAgOU0aNBAhg8fHuhsAADgE5R4AwAAAADgRwTeAADAclJTU+X55583KQAAdkfgDQAALKdOnTpyxhlnmBQAALujjTcAALAcDbh79OgR6GwAAOATlHgDAADLyc/Pl507d5oUAAC7I/AGAACWk5GRIW+//bZJAQCwO6qaAwAAy0lKSpIxY8ZIfHx8oLMCAMBRI/AGAACWEx4eLomJiYHOBgAAPkFVcwAAYDmHDh2SxYsXmxQAALsj8AYAAJZTUFBgOlfTFAAAu6OqOQAAsGQb75EjRwY6GwAA+AQl3gAAAAAA+JEjS7xdLpdJrdxurKSkRPLy8sz7kJCQQGcHAIBjav/+/TJv3jwZMmSIKf126r2KVrUPDaWcBACsyB1Pus/Z1Qlx1WSuILN7925p1qxZoLMBAAAAALC5Xbt2SdOmTaudx5GBt5Ym//7772ZsUEqTrf0ESR+Q6IGckJAQ6OzgGGP/Oxv7HxwDzsb+B8eAsx2ySRygoXRWVpYcf/zxR6yd5Miq5rpRjvREAtah/9is/A8O/sX+dzb2PzgGnI39D44BZ0uwQRxQt27dGs1HoyEAAAAAAPyIwBsAAAAAAD8i8IZlRUVFycMPP2xSOA/739nY/+AYcDb2PzgGnC0qCOMAR3auBgAAAADAsUKJNwAAAAAAfkTgDQAAAACAHxF4AwAAAABA4A2n27ZtmyxfvlwOHjwY6KzgGPv9999l3bp1kpmZybZ3IP03v3btWklNTQ10VhAg+/btM+d/zgHO8ttvv8nq1aslKysr0FlBABQVFcmqVatk06ZNbH8HKigokB9//FF27NghJSUlEiwo8YYtAq8zzzxTzj77bHMDDmdYsWKF9OjRQ3r27CkjRoyQ4447TkaNGmUuxgh+W7ZskSFDhsgJJ5wgI0eOlFatWsmgQYMkIyMj0FnDMaIP3K666io57bTTzPl/5cqVbHsHyMnJMf/WTz31VBk+fLgkJyfLa6+9Fuhs4Rju/3//+9/mnH/BBRfIgw8+yLZ32P6/6667zD3ftddea+4B27dvL99++60EAwJvWJo+5bruuutM4AXnlXZMnTrVpHoDrg9d3n77bZkyZUqgs4ZjYOvWrXLDDTdIenq62ffbt2+XX375RW677Ta2v0Ns3LhRLr30Uh64OswDDzwgGzZsMP/mf/75Z/nPf/4jY8eONdcBBD8956uvv/5a+vTpE+jsIAD7v0mTJrJr1y754YcfTNq9e3fzIL64uNj2+4PAG5b26KOPSp06dcxFF84ybNgwadeuXennNm3amCege/bsCWi+cGz0799fBg8eLCEhIeZzUlKSDB061FQ5hjNoacfVV18tkZGRgc4KjuHD9unTp5trfqNGjcw0PQZat24t06ZNYz84QLNmzUyJd9OmTQOdFQRo/995550SGxtrPoeHh8vo0aPljz/+MA/k7S480BkAqqI32Fq9bM2aNabqCZwnLy/PtPHT/T9//nxzU6YnYDjTd999Z27AAQQnbc+pzUm6du1abnq3bt2o+QA4+NofGRlpgnK7I/CGJR04cMCUdrz++uvmqbdejOHMKkf33Xef6WBLqx3+85//5Cm4Q2lp12effSaff/55oLMCwI/XftWgQYNy0/WzVjsF4CybN282NSDuvvtuiYmJEbsj8EbAaVsud2+1YWFhpiMF/QemJVvx8fGm5Hvv3r3me+3hUANx7XQFwUH3vR4DbrrftTMdpe183FWLta3fOeecIy6Xy7QBRPAoW6tFq5d16dKl3Pcffvih3HTTTTJp0iQ599xzA5RL+LPXcm2/76ZNTBITE9ngDhQREVFa26ms3NxcmhwADrNz50656KKLpG/fvvJ///d/EgwIvBFw2nGKViNx33R/+umnUrduXcnPzzelnUrfK616vnv3bnnyyScDmmf49mmmez+re+65x7Ttraht27YycOBAWbhwIYF3kJkwYYLpRE+1aNFCZs6cWfrdRx99ZNp2P/PMM3LzzTcHMJfwFz3/jx8/vvTzU089ZUaygPPov39VsS8P/dy8efMA5QpAIILu8847Tzp37mw61tWCuWAQ4tLiI8DitKq5Diu0dOlS8w8Rwe/w4cOmY72ytMRbS8K0vTeCnz5k0V6tNRC79dZbA50dBMj+/fulYcOG5njo168f+yHIaXtufdD65ptvms/a1Eg71tSHb3/7298CnT0cQ9qTtXauNXfuXLa7g+zevdvUbtOhJN99993SmjDBgBJvAJak43cOGDDAVDvWGg/vvPOOKRmjja8zLFu2zATdl112mTkG3E0OQkNDKQ11UMC9adMmE3i5hxeLi4szHey4S0YRfB5//HG5+OKLpWXLlqaTteeee86UduvwgnCGb775xgwdpW3+taRTz/8afPXo0SPQWcMx6Nund+/epmmJPnBftWpV6XcdOnQwNWLtjBJv2IK28b7iiivkxRdfNNVOEPz0Znvy5MnmAqwXXm3Xr0PMUN3QGXQMd31VFBUVJUuWLAlInnBs6X5++OGHK02/7rrrGGIyyH355Zfy8ssvS1pamnTq1Mk0R9JaD3CGCy+8sNJoNvXq1ZMFCxYELE84Nn766ScZM2aMx++CIQYg8AYAAAAAwI9C/blwAAAAAACcjsAbAAAAAAA/IvAGAAAAAMCPCLwBAAAAAPAjAm8AAAAAAPyIwBsAAAAAAD8i8AYAAAAAwI8IvAEAjrBy5UpZtWrVUS9n9erV8s0334i/7du3T9555x3LLAcAANReiMvlch3F7wEA8JmCggJ57733Sj/Xq1dP2rZtKy1atDjqZV933XUSHh4u06ZNO6rljB07VjIzM6sMZtesWSNbtmwx7/XvNWnSRDp37izR0dFe/Z3FixfLBRdcIN5cpvfu3SvLli2TK6+88qiW4y1Pf9db+lBk+/bt1c5z1VVXSUhIiPjrwYwuu0ePHmIVVswTAKB2wmv5OwAAfO7QoUMybNgwOfvss+X444+X9PR0E9CNHj1aJk2adFTL7tmzp4SG+r+i19SpU2XWrFly4YUXSlFRkaxfv15yc3Pl/fffl27duvn1b69bt06uueaacgFwcnKyCViP9d/1lj6w+PLLL0s/z54922yvE088sXSaLt9fgbceX/qgxEpBrhXzBACoHQJvAIDl/OMf/5CBAwea9x999JF5P2jQILnooovMtLVr18pvv/1mSsI7depULhhbvny5xMbGmoBNSwzDwsJMiW/Xrl0rBW0aGGu18f3795uS9TZt2lTKS05Ojgn+Y2JiTMl1TejfdpeIFxcXy/nnny9///vfK1V115Li7777TuLi4qRLly5St27dKpd58OBBWbhwoXkfFRUlrVq1ktNOO630+4yMDBO4asm2+2+ffPLJpsR9yJAh5nNeXp7MmzdP+vbtK0lJSaW/1QcD8+fPLze9pnmr6u/qb5SW/m/cuFHq168vZ555pkRERHhczs0332xeZQPvMWPGyKhRo8xn3d8ffvihORZ0O+7Zs0cuvvhiqVOnTqVjomPHjuUesvzwww/y008/mfeJiYlmux133HGl37t/q79xr4MeM5s2bTLbWtdHHy5kZ2fLueeea/6m1npYsWKFOdZ0vSIjIyutU3V5+vrrr82y9ZjT+XTfaICdkJBQbZ4aNGjgcfsBAKyNwBsAYGkaXGlQo4GPloBeeumlsnPnThM8bdiwwQSWGjRqQKWeeeYZ+eOPP0wwfcopp5hgWQOWiqWHGtRoIK/V2zWw0iBq6NCh8p///Kf0b2/dulV69+5tqomfcMIJJhBr1KiRCXprSgN/Lf3+f//v/5Wb/sQTT8iTTz4pZ5xxhgl8NTB86623Sh8uVJSVlWWCZqVBmgaf7du3lwULFpgATgNgfYigAbB7Pl3WgQMHTC2Cq6++2qzHQw89JDt27JD77ruvdNk6vwa52h7c27xV9Xc18NZAWn+ngaluS33w8cknn5Qrxa6pr776SsaNG2cCWH1g0rx5cznvvPPM/rvssstMNXX9ToP8xo0bywcffFAapOo0PUaUHhea3/Hjx8stt9xipv3444+ya9cukz/3Opx++uny3HPPmXxrzYsOHTqY/a9/T/flI488IqeeeqpZth4TGki7HyroNjlSnnTZehxrAN+6dWvzXn+nD470WKsqTwTeAGBT2sYbAAArSEtL04bIrg8//LB02vbt28206dOnu6688krXVVdd5SoqKjLfFRQUuC644ALXuHHjSucfPHiwKyYmxrVly5Zyy7722mtdKSkppZ+HDBniOvfcc115eXnm88aNG11RUVGu2bNnl84zaNAg10UXXeQqLCw0n1esWOEKCQkxeaiK5qVr167lpo0cOdLVrFmz0s8ff/yxKykpyayb27Rp01yNGjVyHT582Hz+7LPPzHpXRedr166da+LEiaXTFi5c6AoLCys3X8XlPPzww6727duXm2fAgAGu4cOH1zhvFXn6u/PmzXNFRka6fvjhB/M5Pz/fdf7557v69+/vqgnN8+uvv176ecaMGWbac889V26+a665xnXFFVeU7iNN+/Xr57rpppuqXPby5cvNvt6zZ0+Vx4e6/PLLXQkJCa4dO3aYz7r+uh3q169f+tvMzEwzT9njpiZ5ci9727Zt5nNxcbGrR48erltvvbXaPAEA7IkSbwCA5Wjpplbr1dLayZMnmxJBLbUeOXKkPPjgg6a9tMZm+tJqvEuXLi33+379+slJJ51U5fK1xFhLQPWlpcVKSy+1NF2r9WpbYi3p1dJkrd6tJeXuduLa/vxItORSl6PVzL///nuZOXOmvPbaa6Xfv/HGG+bvaQ/pWp1b10NLxlNTU03paFVtwUtKSkxb6N27d5t10FLfb7/9VrztZE5La7X6tdYa0BLgTz/91Kzr0eStIl1/rRauJcVKay3cc889Zt9oW353lWpvaD7+9re/lX4+fPiwvPvuu6b0XkuF3ceEbpeKx4SWLGutCV0P3Y5aOq3t77Uvgepoft2d+2m1ct1m+hv377QKvm4vd4d63uRJawZo6bbSKuV6bGlJNwAg+BB4AwAsR6sCa1VwDWpuuukmueGGG0y1Ww1kNRjUKr9l9erVq9znsu13PdFqvRoMVazyrFXI3e2o9e/pPC1btiw3jwZKGvQeKfDWoEurJWsVdg2orrjiitLvtaq3PliYO3duud9pJ2juIL8iDba1DbY+ENAq5vHx8WY5hYWF4g19iNG9e3fzMECDSA0StfqytkOvbd480f1Xcb+4q+jrd+6A3Buaz7JtqbWdt1Y714cbv/zyS7l5zzrrrNL3Wt1dA3ZtT920aVOzDD2WNAg/EncTBjd9UONpmvuYqGmelLZ7r2o5AIDgQuANALB052pu7hLS2267zbSZrs6Rer5u2LChSbVEvSz97O5czN2WVoPosvSzdrRW087VdH5tV67tid3tx3VdNADW3s9r6vHHHzdB42effVa6ftrxmAbK3tJS76efftq0c9YAXNuAa2lybfPmiW5HT9vX/V1tVNyv7mNCO67TvgA80Ycn2jb8+eefNw9w3NO0gzR/DLFWkzwBAJzH/+OqAADgAxp0tmvXTl555ZVK32kpoze0xFJ7MS87Znh+fr7pNdtdKqmBt5aQuju2cgfRFasL1+Rvvfjii6YKt5bWu6sv69+qmO/q1kN7GddO4NzBp1ZpXrRoUbl5tAdyLck9Uim4ll5rB3TTp083JfIaiLvVJm+e/q5ux48//rhcCe6cOXNMqfeRaiTUlA6Vpp2XVXdMaA0BrdquHe25aU/5Or3iOviitLkmeaopX+UJABB4lHgDAGxjypQpJjDs37+/DB482FSJ1uBTe99+9NFHvVrWxIkTTam6Bota7XnGjBmmqu8dd9xROo/27H355ZeboFyD3tdff93jsFFHom15tXd0LcnXwF1LQzW41ZJwLY3VIF+Hj1qyZElpW+GKdEiwsWPHml7cNZjXbaHBd1na1ljbId97772mOrnm2RPthVvbzGspvD6AcA/9pWqTN09/V5etY5rreqekpMjPP/8sL7/8crmHHb6gbed1++pL2+jrMaFt1nX4OO2dXfOl66o9rOu6/f7776bfAO3hvSztMfz++++XV1991TRx0N/4K0815SlP9GoOAPZE4A0AsAwNfLU0VoNLTzTA1qGtpk2bZtqBa+midramw0q5aXvqim1n3R2jlR1HWaur65BcGnDrUFA6TrgOqaWljG6XXHKJLF682FS73rx5szz22GMmcKsY8JalwZWnMa8nTJhgqnZr1XBtN66BrI5VrWOEa5tnDX6fffbZ0vl13XRbuF1//fWmXbe2Qdc26hrkaimzDlnlpuut+X377bdNx3G6jtoZWtnluOkDBq0WretYlgakR8pbRZ7+rv5Gt68GobqvdB4tXddgsiY0z2WHbdNt5h6PvCwN9MseEzpsl3Zs1qdPn9J5NNjXYFuH6tKHDp9//rl5CODu2EyNGDHCdLqmndVpoKz51FJ73eZlnXPOOWbflKV/q+wY8DXJk6dla5v7sg92POWJwBsA7ClEuzYPdCYAAAAAAAhWtPEGAAAAAMCPCLwBAAAAAPAjAm8AAAAAAPyIwBsAAAAAAD8i8AYAAAAAwI8IvAEAAAAA8CMCbwAAAAAA/IjAGwAAAAAAPyLwBgAAAADAjwi8AQAAAADwIwJvAAAAAAD8iMAbAAAAAADxn/8PtifsWPPAPw0AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
event_timeattseconf_int_lowerconf_int_upper
0-4-0.01560.0731-0.15880.1276
1-3-0.00040.0520-0.10230.1014
2-20.06180.0501-0.03640.1599
3-10.0000NaNNaNNaN
402.04340.04921.94692.1399
512.30500.04962.20782.4022
622.71210.07602.56312.8611
\n", + "
" + ], + "text/plain": [ + " event_time att se conf_int_lower conf_int_upper\n", + "0 -4 -0.0156 0.0731 -0.1588 0.1276\n", + "1 -3 -0.0004 0.0520 -0.1023 0.1014\n", + "2 -2 0.0618 0.0501 -0.0364 0.1599\n", + "3 -1 0.0000 NaN NaN NaN\n", + "4 0 2.0434 0.0492 1.9469 2.1399\n", + "5 1 2.3050 0.0496 2.2078 2.4022\n", + "6 2 2.7121 0.0760 2.5631 2.8611" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "es = res.aggregate(\"event_study\")\n", + "fig = plot_event_study(es, title=\"DMLDiD event study (sieve nuisances)\", show=False)\n", + "plt.show()\n", + "es.to_dataframe()[[\"event_time\", \"att\", \"se\", \"conf_int_lower\", \"conf_int_upper\"]]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "344ecd88", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.710646Z", + "iopub.status.busy": "2026-08-29T11:56:54.710584Z", + "iopub.status.idle": "2026-08-29T11:56:54.714353Z", + "shell.execute_reply": "2026-08-29T11:56:54.714087Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
labelattseconf_int_lowerconf_int_upper
042.37580.06002.25812.4934
152.14180.06422.01602.2677
\n", + "
" + ], + "text/plain": [ + " label att se conf_int_lower conf_int_upper\n", + "0 4 2.3758 0.0600 2.2581 2.4934\n", + "1 5 2.1418 0.0642 2.0160 2.2677" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "res.aggregate(\"group\").to_dataframe()[[\"label\", \"att\", \"se\",\n", + " \"conf_int_lower\", \"conf_int_upper\"]]\n" + ] + }, + { + "cell_type": "markdown", + "id": "4fe239d2", + "metadata": {}, + "source": [ + "The event-study estimates track the true dynamic path $\\tau(e) = 2.0 + 0.3e$\n", + "(**2.04, 2.30, 2.71** against true 2.0, 2.3, 2.6), pre-treatment effects sit near\n", + "zero, and the group effects order correctly (cohort 4, treated longer under a growing\n", + "effect, averages higher: **2.38** vs **2.14**).\n", + "\n", + "**HonestDiD** (Rambachan & Roth 2023) bounds the effect under relaxations of\n", + "parallel trends: the relative-magnitude restriction bounds each post-treatment\n", + "*consecutive first difference* of the violation path by `M` times the largest\n", + "pre-treatment first difference (it does not cap violation *levels*, which can\n", + "accumulate over horizons). The event-study container\n", + "feeds it directly - this is why the fit used `base_period=\"universal\"` (the\n", + "varying-base container warns: consecutive-period pre-treatment comparisons are not\n", + "the level deviations HonestDiD's $\\delta$ bounds describe).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "a72978e5", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.715414Z", + "iopub.status.busy": "2026-08-29T11:56:54.715350Z", + "iopub.status.idle": "2026-08-29T11:56:54.718811Z", + "shell.execute_reply": "2026-08-29T11:56:54.718521Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "======================================================================\n", + " Honest DiD Sensitivity Analysis Results \n", + " (Rambachan & Roth 2023) \n", + "======================================================================\n", + "\n", + "Method: Relative Magnitudes (Delta^RM)\n", + "Target: Equal-weight avg over post horizons\n", + "Restriction parameter (M): 1.0000\n", + "CI method: FLCI\n", + "\n", + "----------------------------------------------------------------------\n", + " Original Estimate (under parallel trends) \n", + "----------------------------------------------------------------------\n", + "Point estimate: 2.3535\n", + "Standard error: 0.0458\n", + "\n", + "----------------------------------------------------------------------\n", + " Robust Results (allowing for violations) \n", + "----------------------------------------------------------------------\n", + "Identified set: [2.2291, 2.4779]\n", + "95% Robust CI: [2.1393, 2.5677]\n", + "\n", + "Effect robust to violations: Yes\n", + "\n", + "Pre horizons used: [-4, -3, -2]\n", + "Post horizons used: [0, 1, 2]\n", + "\n", + "----------------------------------------------------------------------\n", + " Interpretation \n", + "----------------------------------------------------------------------\n", + "Post-treatment first differences bounded at 1.0x max pre-period first difference.\n", + "Effect remains POSITIVE even with violations up to M=1.0.\n", + "\n", + "======================================================================\n" + ] + } + ], + "source": [ + "honest = compute_honest_did(es, method=\"relative_magnitude\", M=1.0)\n", + "print(honest.summary())\n" + ] + }, + { + "cell_type": "markdown", + "id": "cc0e1fa9", + "metadata": {}, + "source": [ + "Even allowing each post-treatment *change* in the violation path to be as large\n", + "as the worst pre-treatment change (`M=1.0`), the robust CI **[2.14, 2.57]** stays\n", + "far from zero. One inheritance to keep in mind: this robust CI is built from the\n", + "section-3 fit, whose propensity model this DGP deliberately misspecifies - so\n", + "like every interval in this tutorial, read it as an API illustration rather than\n", + "a coverage guarantee (section 4's rate-condition caveat).\n", + "\n", + "**Uniform (sup-t) bands** need a multiplier bootstrap: refit with `n_bootstrap=199`.\n", + "Point estimates are bit-identical to the analytical fit - the bootstrap only replaces\n", + "inference (the per-cell p-values become `Boot. p` in `summary()`). This fit keeps the\n", + "default `base_period=\"varying\"`: a universal-base fit materializes zero-effect\n", + "reference cells whose bootstrap SE is undefined by construction (the fit would warn),\n", + "and this fit does not feed HonestDiD. `cband=True` (the default) emits the sup-t\n", + "band on the post-fit event-study replay - the band designed to cover the *entire*\n", + "dynamic path at 95%, not each point separately (nominal coverage, under the rate\n", + "conditions this DGP deliberately violates - the section-4 caveat applies here too).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "a32a5862", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.719954Z", + "iopub.status.busy": "2026-08-29T11:56:54.719888Z", + "iopub.status.idle": "2026-08-29T11:56:54.768639Z", + "shell.execute_reply": "2026-08-29T11:56:54.768279Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "analytical ATT 2.2804 | bootstrap ATT 2.2804\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
event_timeattconf_int_lowerconf_int_uppercband_lowercband_upper
0-30.0092-0.13850.1567-0.19680.2152
1-20.0624-0.03560.1735-0.07850.2033
2-1-0.0582-0.15600.0344-0.19470.0784
302.04341.96072.11871.92452.1622
412.30502.22392.40422.18102.4290
522.71212.55502.83042.51662.9075
\n", + "
" + ], + "text/plain": [ + " event_time att conf_int_lower conf_int_upper cband_lower \\\n", + "0 -3 0.0092 -0.1385 0.1567 -0.1968 \n", + "1 -2 0.0624 -0.0356 0.1735 -0.0785 \n", + "2 -1 -0.0582 -0.1560 0.0344 -0.1947 \n", + "3 0 2.0434 1.9607 2.1187 1.9245 \n", + "4 1 2.3050 2.2239 2.4042 2.1810 \n", + "5 2 2.7121 2.5550 2.8304 2.5166 \n", + "\n", + " cband_upper \n", + "0 0.2152 \n", + "1 0.2033 \n", + "2 0.0784 \n", + "3 2.1622 \n", + "4 2.4290 \n", + "5 2.9075 " + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "res_boot = DMLDiD(outcome_learner=\"sieve\", n_folds=5, seed=42,\n", + " n_bootstrap=199, cband=True).fit(df, **fit_kw)\n", + "print(f\"analytical ATT {res.overall_att:.4f} | bootstrap ATT {res_boot.overall_att:.4f}\")\n", + "res_boot.aggregate(\"event_study\").to_dataframe()[\n", + " [\"event_time\", \"att\", \"conf_int_lower\", \"conf_int_upper\",\n", + " \"cband_lower\", \"cband_upper\"]]\n" + ] + }, + { + "cell_type": "markdown", + "id": "a5e45b6c", + "metadata": {}, + "source": [ + "## 7. Repeated Cross-Sections (`panel=False`)\n", + "\n", + "Chang's Case 2 handles *repeated cross-sections* - fresh units sampled each wave\n", + "(think BRFSS, ACS, CPS supplements). Rows must be unique units (one observation per\n", + "row), and the estimator warns about the assumption you cannot test from data:\n", + "stationary cross-sectional sampling (Chang's Assumption 2.3 - each wave draws from\n", + "the same target population, so the composition of $(D, X)$ is stable across waves).\n", + "That warning is a caveat about the *lane*, not this dataset; it applies equally to\n", + "every `panel=False` fit in this section.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "96e28ee9", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.769757Z", + "iopub.status.busy": "2026-08-29T11:56:54.769692Z", + "iopub.status.idle": "2026-08-29T11:56:54.821274Z", + "shell.execute_reply": "2026-08-29T11:56:54.820941Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "RCS ATT: 2.4428 +/- 0.2373\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/bh/mzf05nq92hs6t7vn2ssvfhpr0000gn/T/ipykernel_41788/2521273073.py:27: UserWarning: panel=False uses Chang (2020) Case 2 repeated-cross-section scores, which assume stationary cross-sectional sampling (Assumption 2.3): each wave samples the SAME target population — conditional on the period, rows are i.i.d. draws from the distribution of (Y(0), D, X) (pre) or (Y(1), D, X) (post), so the composition of (D, X) is stable across waves while outcomes are the period-specific potential outcomes (trends and treatment effects are expected, not violations). This assumption is not data-checkable.\n", + " res_rcs = DMLDiD(outcome_learner=\"sieve\", n_folds=5, seed=42,\n" + ] + } + ], + "source": [ + "rng2 = np.random.default_rng(11)\n", + "n_rows = 5000\n", + "rx1 = rng2.normal(size=n_rows)\n", + "rx2 = rng2.normal(size=n_rows)\n", + "rscore = 0.9 * (rx1**2 - 1.0) + 0.7 * rx1 * rx2\n", + "rp = np.clip(1.0 / (1.0 + np.exp(-rscore)), 0.05, 0.60)\n", + "rtreated = rng2.uniform(size=n_rows) < rp\n", + "rearly = rng2.uniform(size=n_rows) < 0.5\n", + "rcohort = np.where(rtreated, np.where(rearly, 3, 4), 0)\n", + "rt = rng2.integers(1, 6, size=n_rows) # waves 1..5\n", + "rf = 1.6 * (rx1**2 - 1.0) + 1.0 * rx1 * rx2\n", + "ry = (1.0 + 0.25 * rt + 0.5 * rx1 - 0.3 * rx2 + rf * (rt / 5.0)\n", + " + rng2.normal(scale=0.6, size=n_rows))\n", + "rpost = (rcohort > 0) & (rt >= rcohort)\n", + "ry = ry + np.where(rpost, 2.0 + 0.3 * (rt - rcohort), 0.0)\n", + "\n", + "# Survey structure for the closing example: 20 PSUs in 4 strata (5 each),\n", + "# mean-1 sampling weights.\n", + "rpsu = rng2.integers(0, 20, size=n_rows)\n", + "rw = rng2.uniform(0.5, 1.5, size=n_rows)\n", + "rcs = pd.DataFrame({\n", + " \"unit\": np.arange(n_rows), \"time\": rt, \"y\": ry, \"first_treat\": rcohort,\n", + " \"x1\": rx1, \"x2\": rx2,\n", + " \"psu\": rpsu, \"stratum\": rpsu % 4, \"w\": rw / rw.mean(),\n", + "})\n", + "\n", + "res_rcs = DMLDiD(outcome_learner=\"sieve\", n_folds=5, seed=42,\n", + " panel=False).fit(rcs, **fit_kw)\n", + "print(f\"RCS ATT: {res_rcs.overall_att:.4f} +/- {res_rcs.overall_se:.4f}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "55db2545", + "metadata": {}, + "source": [ + "Repeated cross-sections are where complex survey data usually lives, and DMLDiD\n", + "accepts a `SurveyDesign` at `fit()` time (it is a `fit()` keyword, not a constructor\n", + "parameter). A declared design weights the moments and the learners (`sample_weight`),\n", + "makes cross-fitting folds **PSU-cohesive** (whole PSUs enter a fold together, so\n", + "within-cluster dependence cannot leak across the train/score split), and switches\n", + "inference to a design-based $t$ with $df = n_{PSU} - n_{strata}$:\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "47179b76", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.822373Z", + "iopub.status.busy": "2026-08-29T11:56:54.822307Z", + "iopub.status.idle": "2026-08-29T11:56:54.881982Z", + "shell.execute_reply": "2026-08-29T11:56:54.881610Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "survey ATT: 2.2886 +/- 0.1509\n", + "design df: 16 (20 PSUs - 4 strata)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/bh/mzf05nq92hs6t7vn2ssvfhpr0000gn/T/ipykernel_41788/4091951892.py:1: UserWarning: panel=False uses Chang (2020) Case 2 repeated-cross-section scores, which assume stationary cross-sectional sampling (Assumption 2.3): each wave samples the SAME target population — conditional on the period, rows are i.i.d. draws from the distribution of (Y(0), D, X) (pre) or (Y(1), D, X) (post), so the composition of (D, X) is stable across waves while outcomes are the period-specific potential outcomes (trends and treatment effects are expected, not violations). This assumption is not data-checkable.\n", + " res_svy = DMLDiD(outcome_learner=\"sieve\", n_folds=5, seed=42, panel=False).fit(\n" + ] + } + ], + "source": [ + "res_svy = DMLDiD(outcome_learner=\"sieve\", n_folds=5, seed=42, panel=False).fit(\n", + " rcs, **fit_kw,\n", + " survey_design=SurveyDesign(weights=\"w\", strata=\"stratum\", psu=\"psu\"),\n", + ")\n", + "print(f\"survey ATT: {res_svy.overall_att:.4f} +/- {res_svy.overall_se:.4f}\")\n", + "print(f\"design df: {res_svy.survey_metadata.df_survey} \"\n", + " f\"({res_svy.survey_metadata.n_psu} PSUs - \"\n", + " f\"{res_svy.survey_metadata.n_strata} strata)\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "bcc70076", + "metadata": {}, + "source": [ + "Both fits bracket the design's true overall effect (about 2.24 under the same\n", + "$\\tau(e) = 2.0 + 0.3e$): the unweighted fit gives **2.44 ± 0.24**, the\n", + "survey-weighted fit **2.29 ± 0.15** on **16** design degrees of freedom. Survey\n", + "tutorials proper - weight types, replicate designs, and the rest of the survey\n", + "toolkit - live in the *Survey-weighted DiD* tutorial.\n", + "\n", + "## 8. Seeds and Reproducibility\n", + "\n", + "`seed=` pins the cross-fitting fold draws, and therefore the *point estimates* - an\n", + "unseeded fit re-randomizes folds, so even the ATT moves between runs (`summary()`\n", + "prints the entropy an unseeded fit drew, so it can be reported after the fact). It\n", + "does **not** reach inside a custom learner: a stochastic learner needs its own\n", + "`random_state`-style seeding.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "e92fffe6", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:54.883033Z", + "iopub.status.busy": "2026-08-29T11:56:54.882967Z", + "iopub.status.idle": "2026-08-29T11:56:55.007685Z", + "shell.execute_reply": "2026-08-29T11:56:55.007322Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "seed=1: 2.271637\n", + "seed=2: 2.273391 (different folds -> different estimate)\n", + "seed=1 again: 2.271637 (bit-identical)\n" + ] + } + ], + "source": [ + "a = DMLDiD(outcome_learner=\"sieve\", seed=1, base_period=\"universal\").fit(df, **fit_kw)\n", + "b = DMLDiD(outcome_learner=\"sieve\", seed=2, base_period=\"universal\").fit(df, **fit_kw)\n", + "a2 = DMLDiD(outcome_learner=\"sieve\", seed=1, base_period=\"universal\").fit(df, **fit_kw)\n", + "print(f\"seed=1: {a.overall_att:.6f}\")\n", + "print(f\"seed=2: {b.overall_att:.6f} (different folds -> different estimate)\")\n", + "print(f\"seed=1 again: {a2.overall_att:.6f} (bit-identical)\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "d81d0c7c", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "| Situation | Reach for |\n", + "|-----------|-----------|\n", + "| Covariates with nonlinear / high-dimensional relationships to treatment or trends | `DMLDiD` |\n", + "| A logit + linear-OLS nuisance specification is plausible | `CallawaySantAnna` (fewer moving parts, replicate-weight support) |\n", + "| No covariates at all | `CallawaySantAnna` (DMLDiD requires covariates) |\n", + "| Repeated cross-sections, optionally with a complex survey design | `DMLDiD(panel=False)` + `survey_design=` |\n", + "| Conclusion must survive a learner swap | Refit with `outcome_learner=` alternatives (section 4; this is Baker et al.'s step 8) |\n", + "\n", + "DMLDiD shares Callaway-Sant'Anna's cell architecture and post-fit aggregation\n", + "surface, so everything downstream - event studies, group effects, HonestDiD,\n", + "`DiagnosticReport`, `practitioner_next_steps` - works identically. The methodology\n", + "details (score functions, survey extension, trimming and skip conventions) are in\n", + "`docs/methodology/REGISTRY.md`.\n", + "\n", + "**References**: Chang, N.-C. (2020). Double/debiased machine learning for\n", + "difference-in-differences models. *The Econometrics Journal*, 23(2), 177-191.\n", + "https://doi.org/10.1093/ectj/utaa001 | Rambachan, A. & Roth, J. (2023). A more\n", + "credible approach to parallel trends. *Review of Economic Studies*, 90(5).\n" + ] + }, + { + "cell_type": "markdown", + "id": "cda2d2a9", + "metadata": {}, + "source": [ + "## What next\n", + "\n", + "Fitted results know their own follow-up work: `practitioner_next_steps()` returns\n", + "the Baker, Callaway, Cunningham, Goodman-Bacon & Sant'Anna (2025) practitioner\n", + "workflow steps as generic runnable templates - they carry the fit's design shape\n", + "(panel vs RCS, cluster, survey) but use placeholder column names and their own\n", + "seeds, so adapt them to your specification.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "c45bcda7", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-29T11:56:55.008963Z", + "iopub.status.busy": "2026-08-29T11:56:55.008893Z", + "iopub.status.idle": "2026-08-29T11:56:55.010605Z", + "shell.execute_reply": "2026-08-29T11:56:55.010260Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "============================================================\n", + "Practitioner Guidance — DMLDiD (Chang 2020 double/debiased ML)\n", + "Baker et al. (2025) 8-Step Workflow\n", + "============================================================\n", + "\n", + "Recommended next steps (7 remaining):\n", + "\n", + " * [HIGH] Step 1: Define target parameter\n", + " Why: State explicitly what causal effect you are estimating (ATT, ATT(g,t), weighted/unweighted) and what policy question it answers.\n", + " >>> # What is the target parameter? ATT? Weighted or unweighted?\n", + "\n", + " * [HIGH] Step 2: State identification assumptions\n", + " Why: Name the parallel trends variant you are invoking (unconditional, conditional, PT-GT-NYT, etc.), the no-anticipation assumption, and any overlap conditions.\n", + " >>> # Which PT variant? No-anticipation? Overlap?\n", + "\n", + " * [HIGH] Step 3: Test parallel trends (event-study pre-periods)\n", + " Why: For staggered designs, inspect event-study pre-period coefficients rather than the generic check_parallel_trends() which assumes a single binary treatment with universal pre-periods. Pre-treatment ATTs should be near zero. Use CS post-fit results.aggregate('event_study') or check the estimator's event-study output directly.\n", + " >>> # Inspect pre-treatment event-study coefficients:\n", + " >>> # (available after fitting with event-study aggregation)\n", + " >>> # Pre-period effects should be near zero and insignificant.\n", + "\n", + " * [HIGH] Step 6: Run HonestDiD sensitivity analysis\n", + " Why: Bounds the treatment effect under plausible violations of (conditional) parallel trends. Aggregate the event study post-fit — no refit needed.\n", + " >>> from diff_diff import compute_honest_did\n", + " >>> # Aggregate post-fit; the container feeds HonestDiD directly:\n", + " >>> es = results.aggregate('event_study')\n", + " >>> honest = compute_honest_did(es, method='relative_magnitude', M=1.0)\n", + " >>> print(honest.summary())\n", + "\n", + " * [HIGH] Step 7: Examine group and event study effects\n", + " Why: Aggregate ATT may mask heterogeneity across cohorts or dynamic effects over time. Inspect group and event study aggregations.\n", + " >>> # Aggregate post-fit - no refit needed:\n", + " >>> print(results.aggregate('group').to_dataframe()) # Per-cohort ATTs\n", + " >>> print(results.aggregate('event_study').to_dataframe()) # Dynamic effects\n", + "\n", + " * [HIGH] Step 8: Compare with alternative estimators (CallawaySantAnna, SunAbraham, or ImputationDiD)\n", + " Why: Agreement across estimators with different assumptions strengthens conclusions. Disagreement reveals sensitivity.\n", + " >>> # Re-estimate with CallawaySantAnna, SunAbraham, or ImputationDiD and compare ATT, SE, CI\n", + " >>> # If results agree, confidence increases.\n", + " >>> # If they disagree, investigate which assumptions differ.\n", + "\n", + " * [HIGH] Step 8: Check learner sensitivity\n", + " Why: DML estimates plug in cross-fitted nuisance functions; a conclusion that survives a different outcome learner (e.g. the adaptive polynomial sieve, or a penalized ridge) is more credible than one that depends on the default linear specification.\n", + " >>> # Refit with alternative nuisance learners (carry the fit's\n", + " >>> # design forward so the comparison isolates the learner):\n", + " >>> alt = DMLDiD(outcome_learner='sieve', seed=0).fit(\n", + " >>> df, outcome=..., unit=..., time=..., first_treat=...,\n", + " >>> covariates=[...])\n", + " >>> print(alt.att, results.att) # should be close\n", + "\n", + "============================================================\n", + "\n" + ] + } + ], + "source": [ + "guidance = practitioner_next_steps(res)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials/README.md b/docs/tutorials/README.md index 9ee2ebb1f..d36b8ad04 100644 --- a/docs/tutorials/README.md +++ b/docs/tutorials/README.md @@ -173,6 +173,19 @@ Executed end-to-end against the real framework: a geo-holdout launch experiment - Requires `pip install "google-meridian==1.8.0"` on Python 3.12 (this tutorial only); cannot share an environment with the PyMC-Marketing tutorial (conflicting `arviz` pins); fully self-contained data - Companion drift-test file (`tests/test_t30_mmm_calibration_meridian_drift.py`) +### 31. LWDiD Rolling Transformations (`31_lwdid.ipynb`) +Lee & Wooldridge (2026) rolling-transformation DiD: demeaning (Procedure 2.1) and detrending (Procedure 3.1) turn the panel into a cross-section with exact small-sample t inference: +- Prop 99 (common timing, one treated unit) and Walmart entry (staggered) replications +- Cluster-robust and wild-bootstrap inference, randomization tests, transformation diagnostics +- Companion drift-test file (`tests/test_t31_lwdid_drift.py`) + +### 32. Double ML DiD (`32_dml_did.ipynb`) +Chang (2020) double/debiased machine learning DiD: cross-fitted ML nuisances on a staggered panel with engineered nonlinear confounding: +- The bias demo: a misspecified linear outcome learner misses by almost five standard errors while the adaptive sieve and a hand-rolled duck-typed `PolynomialRidge` learner recover the truth (double robustness in action) +- Cross-fitting diagnostics, event-study/group aggregation, HonestDiD sensitivity, bootstrap sup-t bands +- The repeated-cross-section lane (`panel=False`) with a `SurveyDesign` (PSU-cohesive folds, design-df t inference) +- Companion drift-test file (`tests/test_t32_dml_did_drift.py`) + ## Running the Notebooks 1. Install diff-diff with dependencies: diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 5ce156166..63467834c 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -1,5 +1,5 @@ .. meta:: - :description: Hands-on diff-diff tutorials — 31 Jupyter notebooks covering basic 2x2 DiD, staggered adoption, synthetic DiD, power analysis, and business applications. + :description: Hands-on diff-diff tutorials — 33 Jupyter notebooks covering basic 2x2 DiD, staggered adoption, synthetic DiD, power analysis, and business applications. :keywords: DiD tutorial, difference-in-differences examples, causal inference notebooks Tutorials @@ -245,6 +245,13 @@ Modern estimators for designs the basic toolkit cannot handle. Turn the panel into a cross-section, get exact small-sample inference, and replicate the Prop 99 and Walmart applications. + .. grid-item-card:: Double ML DiD (Chang 2020) + :link: 32_dml_did + :link-type: doc + + Cross-fitted ML nuisances for nonlinear confounding - watch a + misspecified linear model fail where the sieve recovers the truth. + .. toctree:: :maxdepth: 1 @@ -262,6 +269,7 @@ Modern estimators for designs the basic toolkit cannot handle. Synthetic Control for Policy <25_synthetic_control_policy> Regression Discontinuity (RDD) <28_rdd_scholarship_illusion> LWDiD Rolling Transformations <31_lwdid> + Double ML DiD (Chang 2020) <32_dml_did> Study Design ------------ diff --git a/tests/test_t32_dml_did_drift.py b/tests/test_t32_dml_did_drift.py new file mode 100644 index 000000000..c66b72264 --- /dev/null +++ b/tests/test_t32_dml_did_drift.py @@ -0,0 +1,452 @@ +"""Drift detection for Tutorial 32 (``docs/tutorials/32_dml_did.ipynb``). + +The tutorial narrative quotes locked, seed-specific numbers (the DGP-implied +truth, the learner-comparison bias table, the HonestDiD robust CI, bootstrap +sup-t bands, and the RCS/survey fits). ``pytest --nbmake`` only checks that +cells *execute*; it does not check the prose or the committed outputs +(``nbsphinx_execute = "never"`` renders the committed outputs verbatim). +Three layers here: + +1. ``assert_quotes_in_rendered`` pins the load-bearing quoted values against + the committed rendered surface (markdown + outputs). +2. Full re-derivation: the panel and RCS DGPs are rebuilt from the locked + seeds and the section 3/4/6/7 estimates re-checked against the quoted + values - sections 6-7 exercise the survey/cluster lane, so a library + numerics change surfaces here even without re-executing the notebook. +3. ``ALL_CODE_CELL_HASHES`` pins every code cell's normalized source, and + source-fragment pins keep the duplicated DGP/learner definitions below in + sync with the notebook cells they mirror. +""" + +from __future__ import annotations + +import json +import warnings +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from diff_diff import DMLDiD, SurveyDesign, compute_honest_did + +from ._tutorial_drift import assert_quotes_in_rendered, notebook_markdown + +NB = "docs/tutorials/32_dml_did.ipynb" + +# sha256[:16] of EVERY code cell's normalized source, in notebook order - +# the complete stale-output contract (see test_all_code_cells_hash_pinned) +ALL_CODE_CELL_HASHES = [ + "9405e25ece47ce8b", + "68f3f1e50a2a3c60", + "128c08960b3dfd45", + "c48f19ab5c5b76ad", + "d4eed65cbf569b74", + "8ecbf97a1a817414", + "98695413c281f252", + "84d167bbdd9af3c8", + "9efefccaf9379ba7", + "c3573fc2030e099a", + "bbd6bc001d063ff8", + "ffe1220f7721cb2f", + "6b746ea09fd9e3da", + "ed87cf8c4a4cd42f", + "9ee221cd97948c35", +] + +FIT_KW = dict( + outcome="y", unit="unit", time="time", first_treat="first_treat", covariates=["x1", "x2"] +) + + +def _load_nb(): + nb_path = Path(__file__).resolve().parents[1] / NB + if not nb_path.exists(): + pytest.skip("notebook not available in this CI environment") + return json.loads(nb_path.read_text()) + + +def _code_cell_hashes(): + import hashlib + + hashes = [] + for c in _load_nb()["cells"]: + if c["cell_type"] != "code": + continue + src = "".join(c["source"]) + normalized = "\n".join(ln.rstrip() for ln in src.strip().splitlines()) + hashes.append(hashlib.sha256(normalized.encode()).hexdigest()[:16]) + return hashes + + +class TestRenderedSurface: + def test_dgp_and_first_fit_quotes(self): + assert_quotes_in_rendered( + NB, + [ + "DGP-implied overall ATT: 2.2388", + "DMLDiD (sieve) estimate: 2.2804 +/- 0.0438", + ], + surface="output", + ) + md = notebook_markdown(NB) + assert "**2.2804 ± 0.0438**" in md and "**2.2388**" in md + assert "354 never-treated units" in md + + def test_learner_table_quotes(self): + assert_quotes_in_rendered( + NB, + ["2.5909", "2.5898", "2.2804", "2.2818"], + surface="output", + ) + md = notebook_markdown(NB) + for quoted in ( + "**linear 2.5909**", + "**ridge 2.5898**", + "**sieve 2.2804**", + "**PolynomialRidge 2.2818**", + ): + assert quoted in md, f"prose quote missing: {quoted}" + + def test_score_family_and_rm_units_wording(self): + """Review-round pins: Chang's score is a DISTINCT family from the + Sant'Anna-Zhao DR score (REGISTRY DR-score families note), the + relative-magnitude restriction is stated in consecutive + first-difference units, and the downstream HonestDiD/sup-t claims + carry the no-nominal-coverage qualification.""" + md = notebook_markdown(NB) + assert "distinct score family" in md + assert "same doubly-robust moment" not in md + assert "consecutive first difference" in md + assert "coverage guarantee" in md + assert "deliberately violates" in md # sup-t nominal-coverage caveat + + def test_dr_inference_caveat_present(self): + """The bias demo deliberately misspecifies the propensity in every + arm, so the prose MUST carry the Theorem-1 both-nuisances caveat + (plan-review round 3): SEs/CIs in the table are illustrative.""" + md = notebook_markdown(NB) + assert "double robustness" in md.lower() + assert "cannot compensate" in md + assert "illustrative rather than theory-backed" in md + + def test_diagnostics_quotes(self): + assert_quotes_in_rendered(NB, ["0.2547", "0.2610"], surface="output") + md = notebook_markdown(NB) + assert "fold losses do not verify that" in md + + def test_aggregation_and_honest_quotes(self): + assert_quotes_in_rendered( + NB, + [ + "Identified set: [2.2291, 2.4779]", + "95% Robust CI: [2.1393, 2.5677]", + "analytical ATT 2.2804 | bootstrap ATT 2.2804", + ], + surface="output", + ) + md = notebook_markdown(NB) + assert "**[2.14, 2.57]**" in md + assert "**2.04, 2.30, 2.71**" in md + assert "**2.38**" in md and "**2.14**" in md + + def test_rcs_and_survey_quotes(self): + assert_quotes_in_rendered( + NB, + [ + "RCS ATT: 2.4428 +/- 0.2373", + "survey ATT: 2.2886 +/- 0.1509", + "design df: 16 (20 PSUs - 4 strata)", + ], + surface="output", + ) + md = notebook_markdown(NB) + assert "**2.44 ± 0.24**" in md and "**2.29 ± 0.15**" in md + + def test_seed_quotes(self): + assert_quotes_in_rendered( + NB, + ["seed=1: 2.271637", "seed=2: 2.273391", "seed=1 again: 2.271637"], + surface="output", + ) + + def test_source_cells_match_rederived_dgps(self): + """The quote pins read committed OUTPUTS and the rederivation tests + below duplicate the DGPs and the PolynomialRidge learner, so a + source cell edited without re-execution could leave both layers + green. Pin the load-bearing source fragments (seeds + DGP structure + + learner math) so a source change that diverges from the + rederivations fails here.""" + nb = _load_nb() + src = "\n".join("".join(c["source"]) for c in nb["cells"] if c["cell_type"] == "code") + for fragment in [ + "default_rng(7)", # panel DGP seed + "default_rng(11)", # RCS DGP seed + "0.9 * (x1**2 - 1.0) + 0.7 * x1 * x2", # nonlinear assignment score + "1.6 * (x1**2 - 1.0) + 1.0 * x1 * x2", # nonlinear trend + "2.0 + 0.3 * (t - cohort[i])", # true dynamic effect + "np.clip(1.0 / (1.0 + np.exp(-score)), 0.05, 0.60)", + 'base_period="universal"', + "n_bootstrap=199, cband=True", + "A[0, 0] -= self.alpha", # PolynomialRidge intercept exemption + 'SurveyDesign(weights="w", strata="stratum", psu="psu")', + "seed=42,", + ]: + assert fragment in src, f"source fragment missing: {fragment!r}" + + def test_all_code_cells_hash_pinned(self): + """Complete source/output contract: EVERY code cell's normalized + source is hash-pinned, so ANY source edit fails here until the pins + are re-locked together with a fresh execution and updated + rederivation constants.""" + cells = _code_cell_hashes() + assert cells == ALL_CODE_CELL_HASHES, ( + "notebook code cells changed - re-execute the notebook and " + "re-lock ALL_CODE_CELL_HASHES plus any affected rederivation " + f"constants. Got: {cells}" + ) + + def test_hash_guard_detects_mutation(self): + """Negative control: a one-character estimator-argument mutation + must change the cell hash.""" + import hashlib + + nb = _load_nb() + src = next( + "".join(c["source"]) + for c in nb["cells"] + if c["cell_type"] == "code" and 'outcome_learner="sieve"' in "".join(c["source"]) + ) + mutated = src.replace('outcome_learner="sieve"', 'outcome_learner="ridge"', 1) + norm = lambda x: "\n".join(ln.rstrip() for ln in x.strip().splitlines()) # noqa: E731 + h = lambda x: hashlib.sha256(norm(x).encode()).hexdigest()[:16] # noqa: E731 + assert h(mutated) != h(src) + + def test_paper_reference_present(self): + assert_quotes_in_rendered( + NB, + ["The Econometrics Journal", "10.1093/ectj/utaa001"], + surface="markdown", + ) + + def test_notebook_hygiene(self): + nb = _load_nb() + assert nb["metadata"]["kernelspec"]["name"] == "python3" + errors = [ + out + for cell in nb["cells"] + if cell["cell_type"] == "code" + for out in cell.get("outputs", []) + if out.get("output_type") == "error" + ] + assert not errors + + +# --------------------------------------------------------------------------- +# Re-derivation (DGPs mirror the notebook cells exactly; source-fragment +# pins above keep the copies honest) +# --------------------------------------------------------------------------- + + +class _PolynomialRidge: + """Mirror of the notebook's custom learner (degree-2 features + ridge).""" + + def __init__(self, alpha=1.0): + self.alpha = alpha + + def _expand(self, X): + n, d = X.shape + cols = [np.ones(n)] + [X[:, j] for j in range(d)] + for j in range(d): + for k in range(j, d): + cols.append(X[:, j] * X[:, k]) + return np.column_stack(cols) + + def fit(self, X, y, sample_weight=None): + Z = self._expand(np.asarray(X, dtype=np.float64)) + y = np.asarray(y, dtype=np.float64) + w = np.ones(len(y)) if sample_weight is None else np.asarray(sample_weight, float) + ZtW = Z.T * w + A = ZtW @ Z + self.alpha * np.eye(Z.shape[1]) + A[0, 0] -= self.alpha + self.coef_ = np.linalg.solve(A, ZtW @ y) + return self + + def predict(self, X): + return self._expand(np.asarray(X, dtype=np.float64)) @ self.coef_ + + +@pytest.fixture(scope="module") +def panel_df(): + rng = np.random.default_rng(7) + n_units, periods = 600, [1, 2, 3, 4, 5, 6] + x1 = rng.normal(size=n_units) + x2 = rng.normal(size=n_units) + score = 0.9 * (x1**2 - 1.0) + 0.7 * x1 * x2 + p_any = np.clip(1.0 / (1.0 + np.exp(-score)), 0.05, 0.60) + treated = rng.uniform(size=n_units) < p_any + early = rng.uniform(size=n_units) < 0.5 + cohort = np.where(treated, np.where(early, 4, 5), 0) + f_nl = 1.6 * (x1**2 - 1.0) + 1.0 * x1 * x2 + rows = [] + for i in range(n_units): + alpha_i = 0.5 * x1[i] - 0.3 * x2[i] + rng.normal(scale=0.3) + for t in periods: + y = 1.0 + 0.25 * t + alpha_i + f_nl[i] * (t / 6) + rng.normal(scale=0.5) + if cohort[i] > 0 and t >= cohort[i]: + y += 2.0 + 0.3 * (t - cohort[i]) + rows.append((i, t, y, cohort[i], x1[i], x2[i])) + return pd.DataFrame(rows, columns=["unit", "time", "y", "first_treat", "x1", "x2"]) + + +@pytest.fixture(scope="module") +def sieve_fit(panel_df): + return DMLDiD( + propensity_learner="logit", + outcome_learner="sieve", + n_folds=5, + seed=42, + base_period="universal", + ).fit(panel_df, **FIT_KW) + + +class TestPanelRederivation: + def test_cohort_split(self, panel_df): + counts = panel_df.groupby("unit")["first_treat"].first().value_counts().to_dict() + assert counts == {0: 354, 5: 125, 4: 121} + + def test_truth_and_sieve_fit(self, sieve_fit): + num = den = 0.0 + for (g, t), cell in sieve_fit.group_time_effects.items(): + if t < g or cell.get("is_reference"): + continue + w = cell["n_treated"] + num += w * (2.0 + 0.3 * (t - g)) + den += w + np.testing.assert_allclose(num / den, 2.2388, atol=5e-4) + np.testing.assert_allclose(sieve_fit.overall_att, 2.2804, atol=5e-4) + np.testing.assert_allclose(sieve_fit.overall_se, 0.0438, atol=5e-4) + + def test_learner_table(self, panel_df): + expected = { + "linear": (2.5909, 0.0743), + "ridge": (2.5898, 0.0736), + "PolynomialRidge": (2.2818, 0.0438), + } + for name, (att, se) in expected.items(): + learner = _PolynomialRidge(alpha=1.0) if name == "PolynomialRidge" else name + r = DMLDiD(outcome_learner=learner, n_folds=5, seed=42, base_period="universal").fit( + panel_df, **FIT_KW + ) + np.testing.assert_allclose(r.overall_att, att, atol=5e-4, err_msg=name) + np.testing.assert_allclose(r.overall_se, se, atol=5e-4, err_msg=name) + + def test_honest_bounds(self, sieve_fit): + honest = compute_honest_did( + sieve_fit.aggregate("event_study"), method="relative_magnitude", M=1.0 + ) + np.testing.assert_allclose(honest.lb, 2.2291, atol=1e-3) + np.testing.assert_allclose(honest.ub, 2.4779, atol=1e-3) + np.testing.assert_allclose(honest.ci_lb, 2.1393, atol=1e-3) + np.testing.assert_allclose(honest.ci_ub, 2.5677, atol=1e-3) + + def test_bootstrap_supt_bands(self, panel_df, sieve_fit, monkeypatch): + # The notebook was executed on the NumPy weight backend, and the two + # backends draw DIFFERENT (equally valid) multiplier matrices from + # the same seed (REGISTRY weight-backend identity Note) - so force + # the NumPy generator here to match the committed cband goldens on + # Rust-enabled installs (pattern: test_bootstrap_chunking.py). + from diff_diff import bootstrap_chunking + + monkeypatch.setattr(bootstrap_chunking, "_rust_bootstrap_weights", None) + assert bootstrap_chunking.effective_weight_backend() == "numpy" + r = DMLDiD(outcome_learner="sieve", n_folds=5, seed=42, n_bootstrap=199, cband=True).fit( + panel_df, **FIT_KW + ) + # bootstrap overrides inference only; point estimate matches the + # analytical fit across base_period (post cells share the g-1 base) + np.testing.assert_allclose(r.overall_att, sieve_fit.overall_att, rtol=1e-12) + es = r.aggregate("event_study").to_dataframe() + row0 = es[es["event_time"] == 0].iloc[0] + np.testing.assert_allclose(row0["cband_lower"], 1.9245, atol=5e-4) + np.testing.assert_allclose(row0["cband_upper"], 2.1622, atol=5e-4) + + +@pytest.fixture(scope="module") +def rcs_df(): + rng2 = np.random.default_rng(11) + n_rows = 5000 + rx1 = rng2.normal(size=n_rows) + rx2 = rng2.normal(size=n_rows) + rscore = 0.9 * (rx1**2 - 1.0) + 0.7 * rx1 * rx2 + rp = np.clip(1.0 / (1.0 + np.exp(-rscore)), 0.05, 0.60) + rtreated = rng2.uniform(size=n_rows) < rp + rearly = rng2.uniform(size=n_rows) < 0.5 + rcohort = np.where(rtreated, np.where(rearly, 3, 4), 0) + rt = rng2.integers(1, 6, size=n_rows) + rf = 1.6 * (rx1**2 - 1.0) + 1.0 * rx1 * rx2 + ry = ( + 1.0 + + 0.25 * rt + + 0.5 * rx1 + - 0.3 * rx2 + + rf * (rt / 5.0) + + rng2.normal(scale=0.6, size=n_rows) + ) + rpost = (rcohort > 0) & (rt >= rcohort) + ry = ry + np.where(rpost, 2.0 + 0.3 * (rt - rcohort), 0.0) + rpsu = rng2.integers(0, 20, size=n_rows) + rw = rng2.uniform(0.5, 1.5, size=n_rows) + return pd.DataFrame( + { + "unit": np.arange(n_rows), + "time": rt, + "y": ry, + "first_treat": rcohort, + "x1": rx1, + "x2": rx2, + "psu": rpsu, + "stratum": rpsu % 4, + "w": rw / rw.mean(), + } + ) + + +class TestRCSRederivation: + @staticmethod + def _fit_expecting_only_a23(est, rcs_df, **kw): + """Record everything; the ONLY tolerated warning is the deliberate + Assumption 2.3 lane warning - anything else (clipping, degenerate + cells, fold reduction, lonely PSU, weight normalization) fails.""" + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + r = est.fit(rcs_df, **kw) + expected = [ + w for w in caught if w.category is UserWarning and "Assumption 2.3" in str(w.message) + ] + unexpected = [f"{w.category.__name__}: {w.message}" for w in caught if w not in expected] + assert not unexpected, f"unexpected warnings: {unexpected}" + # ... and the documented lane UserWarning must actually FIRE (the + # lane contract, not just be tolerated). + assert len(expected) == 1, [str(w.message) for w in caught] + return r + + def test_rcs_fit(self, rcs_df): + r = self._fit_expecting_only_a23( + DMLDiD(outcome_learner="sieve", n_folds=5, seed=42, panel=False), rcs_df, **FIT_KW + ) + np.testing.assert_allclose(r.overall_att, 2.4428, atol=5e-4) + np.testing.assert_allclose(r.overall_se, 0.2373, atol=5e-4) + + def test_survey_fit(self, rcs_df): + r = self._fit_expecting_only_a23( + DMLDiD(outcome_learner="sieve", n_folds=5, seed=42, panel=False), + rcs_df, + **FIT_KW, + survey_design=SurveyDesign(weights="w", strata="stratum", psu="psu"), + ) + np.testing.assert_allclose(r.overall_att, 2.2886, atol=5e-4) + np.testing.assert_allclose(r.overall_se, 0.1509, atol=5e-4) + assert r.survey_metadata.df_survey == 16 + assert r.survey_metadata.n_psu == 20 + assert r.effective_n_folds is None # 20 PSUs >= 5 folds: no reduction From 545c31fcfbaf26f51be9c499488bf701b849c513 Mon Sep 17 00:00:00 2001 From: igerber Date: Sat, 29 Aug 2026 08:43:50 -0400 Subject: [PATCH 2/2] docs(dml): tutorial 32 - scope the bootstrap bit-identity claim, document the non-copyable-learner fallback PR #800 review P2s: (1) the bootstrap fit keeps the default varying base while the analytical fit is universal-base, so only the overall ATT and post-treatment cells are bit-identical - the prose now says so and flags the differing pre-treatment rows as a base-regime change, not a bootstrap effect; (2) a learner that cannot be deep-copied is not rejected - it is reused across folds with a loud UserWarning relying on fit-reset behavior, so the prose now documents the fallback and the leak risk for stateful non-copyable learners instead of claiming deepcopy is mandatory (_crossfit._fresh_learner contract). Markdown-only; code cells, outputs, and hashes untouched. --- docs/tutorials/32_dml_did.ipynb | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/docs/tutorials/32_dml_did.ipynb b/docs/tutorials/32_dml_did.ipynb index cbf2d59ad..7ac5de43f 100644 --- a/docs/tutorials/32_dml_did.ipynb +++ b/docs/tutorials/32_dml_did.ipynb @@ -491,10 +491,14 @@ " **both** nuisances converging at $o(N^{-1/4})$ - \"a fast learner cannot compensate\n", " a slow one\" - so with a deliberately misspecified propensity, the SEs and CIs in\n", " this table are illustrative rather than theory-backed.\n", - "- Custom learner objects are **deep-copied, never-fit, once per fold** - so they must\n", - " survive `copy.deepcopy`; under `survey_design=` (section 7) the `sample_weight`\n", - " keyword becomes mandatory, and stochastic learners need their own internal seeding\n", - " (`seed=` pins folds, not your learner's RNG).\n" + "- Custom learner objects are **deep-copied, never-fit, once per fold**. A learner\n", + " that cannot be deep-copied still runs, but with a loud `UserWarning`: the same\n", + " instance is REUSED across folds relying on its fit-reset behavior, so a\n", + " warm-start/stateful non-copyable learner can leak data across folds - make such\n", + " a learner fully re-initialize on every `fit()`. Under `survey_design=`\n", + " (section 7) the `sample_weight` keyword becomes mandatory, and stochastic\n", + " learners need their own internal seeding (`seed=` pins folds, not your\n", + " learner's RNG).\n" ] }, { @@ -1005,9 +1009,11 @@ "a coverage guarantee (section 4's rate-condition caveat).\n", "\n", "**Uniform (sup-t) bands** need a multiplier bootstrap: refit with `n_bootstrap=199`.\n", - "Point estimates are bit-identical to the analytical fit - the bootstrap only replaces\n", - "inference (the per-cell p-values become `Boot. p` in `summary()`). This fit keeps the\n", - "default `base_period=\"varying\"`: a universal-base fit materializes zero-effect\n", + "The overall ATT and every post-treatment cell are bit-identical to the analytical\n", + "fit - the bootstrap only replaces inference (the per-cell p-values become\n", + "`Boot. p` in `summary()`). This fit keeps the default `base_period=\"varying\"`,\n", + "so the *pre-treatment* event-study rows below differ from section 6's\n", + "universal-base rows (a different comparison regime, not a bootstrap effect): a universal-base fit materializes zero-effect\n", "reference cells whose bootstrap SE is undefined by construction (the fit would warn),\n", "and this fit does not feed HonestDiD. `cband=True` (the default) emits the sup-t\n", "band on the post-fit event-study replay - the band designed to cover the *entire*\n",