---
title: "Village officials and rural subsidies"
subtitle: "He and Wang (AEJ Applied, 2017) · Table 5, Panel A, Column 1"
author: "Replication comparison"
date: "2026-10-10"
format:
  html:
    theme: cosmo
    toc: true
    embed-resources: true
    self-contained-math: true
    html-math-method: katex
    code-copy: true
    include-in-header:
      text: |
        <style>
        body{font-family:system-ui,sans-serif;line-height:1.6}
        main,.cell,.cell-output-display{min-width:0;max-width:100%}
        .cell-output-display{overflow-x:auto}
        table{display:block;overflow-x:auto;font-size:.85rem}
        img{max-width:100%;height:auto!important}
        p code,li code{white-space:normal;overflow-wrap:anywhere}
        .katex-display{overflow-x:auto;max-width:100%}
        </style>
jupyter: python3
execute:
  warning: false
  error: false
  cache: false
  fig-dpi: 140
---

[He and Wang (2017)](https://doi.org/10.1257/app.20160079) study college graduates
serving as village officials in rural China. This example uses the deposited
village-year panel and the subsidy outcome from Table 5, Panel A, Column 1.
It applies the same estimator and reporting code as the
[existing replication collection](https://lalten.org/pages/did-replication/).

**The original coefficient reproduces.** The adoption-template comparison uses
a smaller, explicitly defined sample. Its regression, CS-DID and matrix-completion
fits succeed; the other counterfactual adapters encounter missing-data or
convergence limits, recorded below.

[Results JSON](he-wang-village-officials.json) ·
[Source replication check](he-wang-source-check.json) ·
[Executable QMD](he-wang-village-officials.qmd) ·
[Code and configurations](he-wang-example-source.zip)

```{python}
#| tags: [parameters]
study_file = "configs/new-library/he-wang.json"
comparison_file = "configs/new-library/he-wang-comparison.json"
reuse_results = True
```

```{python}
#| label: setup
#| output: false
import json
import sys
from pathlib import Path

ROOT = next(
    p
    for p in [Path.cwd(), *Path.cwd().parents]
    if (p / "templates/expensive/runner.py").exists()
)
sys.path.insert(0, str(ROOT / "templates/expensive"))
from runner import Comparison, markdown_table

settings = json.loads((ROOT / study_file).read_text())
options = json.loads((ROOT / comparison_file).read_text())
run = Comparison(settings, options, ROOT)
source = json.loads(
    (ROOT / "analyses/new-library/he-wang/source-check.json").read_text()
)
sample = run.prepare()
```

```{python}
#| label: estimate-models
#| output: false
run.estimate(reuse=reuse_results)
for spec in run.base.results["specifications"].values():
    for result in spec["estimates"]:
        reason = spec.get("method_status", {}).get(result["method"])
        if reason:
            for field in ["estimate", "se", "pvalue", "lower", "upper"]:
                result[field] = None
            result.update(status="unvalidated", reason=reason)
            spec["event_studies"].pop(result["method"], None)
```

## Original specification

The deposited command is
`reg l_subsidy_rate cgvo v_fe* t1-t11, cluster(sm)`:
village and year fixed effects, no additional controls, with standard errors
clustered by town/subdistrict. The benchmark below is the library's saved
author-code run, not a rounded number transcribed from the article.

```{python}
check = source["source_fit"]
markdown_table(
    [
        {
            "Fit": "Saved author Stata",
            "Coefficient": check["saved_author_stata_coefficient"],
            "SE": check["saved_author_stata_se"],
            "Observations": 2102,
        },
        {
            "Fit": "Independent PyFixest",
            "Coefficient": check["coefficient"],
            "SE": check["se"],
            "Observations": check["n"],
        },
    ]
)
```

These coefficients are in the original log-outcome units. The PyFixest check
counts all fixed effects in the finite-sample covariance adjustment, matching
the deposited dummy-variable regression.

## Data and adoption sample

```{python}
run.base.describe()
run.base.df.head()
```

The source contains 2,809 records for 255 villages in 2000–2011; 2,102 have the
subsidy outcome. Treatment records are used to construct first adoption before
dropping missing outcomes. Some officials leave, so treatment is not strictly
absorbing: the template excludes untreated records after first adoption rather
than relabeling those records as treated. It also removes initially treated
villages and villages without an observed adoption year and year −1 baseline.
No missing outcomes are filled with zero.

```{python}
s = sample.iloc[0]
markdown_table(
    [
        {"Quantity": "Original regression observations", "Value": check["n"]},
        {"Quantity": "Template observations", "Value": int(s["rows"])},
        {"Quantity": "Villages", "Value": int(s["units"])},
        {
            "Quantity": "Treated / never-treated villages",
            "Value": f"{s['treated_units']} / {s['never_treated_units']}",
        },
        {"Quantity": "Adoption cohorts", "Value": int(s["treated_cohorts"])},
        {"Quantity": "Town clusters", "Value": int(s["clusters"])},
        {"Quantity": "Outcome SD before template exclusions", "Value": s["source_sd"]},
    ]
)
```

All comparisons below use $Y/s_Y$, where $s_Y$ is the unweighted standard
deviation of observed source outcomes (`ddof=1`). The source check and the
template results therefore differ in both sample and scale. The excluded
records and missing cells remain documented in the JSON.

```{python}
#| label: adoption
#| fig-alt: "Treatment timing for a seeded sample of 100 villages. White cells are unobserved."
run.base.panelview()
```

## Estimates

Static TWFE is the treatment regression coefficient. Dynamic TWFE, BJS/did2s,
Sun–Abraham, matched CS-DID and matrix completion average over the observed
treated post-adoption cells. CS-DID also reports its equal-event-time average;
long horizons receive more relative weight in that summary.

```{python}
rows = []
for row in run.rows():
    if row["Status"] != "estimated":
        continue
    rows.append(
        {
            "Estimator": row["Estimator"],
            "Estimate (SDs)": row["Estimate (SDs)"],
            "SE": f"{row['SE']:.5f}" if row["SE"] is not None else "—",
            "Target": "Treated cells"
            if row["Aggregation"] == "observations"
            else row["Aggregation"],
        }
    )
markdown_table(rows)
```

Regression and CS-DID intervals use the 19 town clusters. Matrix completion
is a point estimate here; no refit bootstrap was run.

```{python}
#| label: forest
#| fig-alt: "Estimator comparison in outcome standard deviations, with available 95 percent intervals."
sys.path.insert(0, str(ROOT / "templates/reporting"))
from comparison_report import export_report, forest_plots

forest_files = forest_plots(run)
```

### Adapter coverage

The same input sample was attempted for every enabled estimator. Multisynth's
current adapter requires a complete panel. Horizontal ridge lacks donor-history
overlap for an early cohort. Masked SDID did not converge with either 1,000 or
10,000 iterations; its acceptance criteria were unchanged.

```{python}
markdown_table(
    [
        {"Estimator": r["Estimator"], "Status": r["Status"], "Reason": r["Note"]}
        for r in run.rows()
        if r["Status"] != "estimated"
    ]
)
```

### Tuning and fit checks

Matrix completion retains the template's nuclear-norm penalty fraction of
0.25, with tolerance $10^{-12}$ and a 1,000-iteration limit. It converged in
45 iterations. These settings were not selected using the treatment effect.
CS-DID uses not-yet-treated controls, including never-treated villages.
For this unbalanced panel, the existing CS-DID adapter uses its
repeated-cross-section implementation on the available observations, with
town-clustered influence functions; it does not silently balance the panel.

```{python}
spec = next(iter(run.base.results["specifications"].values()))
markdown_table(
    [
        {"Check": name.replace("_", " "), "Absolute/relative difference": value}
        for name, value in spec["checks"].items()
    ]
)
```

The checks compare did2s with direct untreated-outcome imputation,
Sun–Abraham aggregation with explicit cell weights, and the saturated
regression with the imported deviations design.

## Cohort homogeneity

The design comes from `TestingInEventStudies`. As in the existing collection,
the reported statistic is the **homoskedastic nested-model F-test** of
cohort deviations at positive event times; adoption-year deviations remain
unrestricted. This is separate from the town-clustered estimator intervals.

```{python}
#| label: homogeneity
run.base.f_report()
```

## Event studies

All available paths are shown without folding. The year −1 normalization is
used by the regression event studies; the matrix-completion path is its fitted
outcome gap. Shading shows available 95% intervals. Cohort composition changes
with event time.

```{python}
#| label: event-studies
run.event_plots()
```

## Reproduction

The full library is read through
`~/Desktop/big-did-replication/DiD-Replication-Library`, a symlink to
`/Volumes/LAL_BAK/DiD-Replication-Library`. Its original package ZIP is never
modified. Only the required Stata member is read into the local analysis workspace.

From the repository root:

```bash
.venv/bin/python scripts/prepare_he_wang.py
QUARTO_PYTHON="$PWD/.venv/bin/python" quarto render \
  reports/new-library/he-wang-village-officials.qmd
```

Preparation verifies the author-code coefficient and saves the source/package
hashes. The two JSON configurations specify columns, sample rules, targets and
solver settings. `Comparison` calls the existing PyFixest, CS-DID,
crabbymetrics and augsynth adapters; the F-test uses the existing reference
environment. Cached fits are validated against their inputs and code hashes.
Set `reuse_results=False` in the QMD to refit.

```{python}
#| label: export
#| output: false
results_path = export_report(run)
```