feat(3.9): variance-conventions inventory + component-aware absorbed-FE rank (D3)

The measured map behind the 3.9 variance-consolidation program, plus the
first of its fixes.

Inventory (docs/methodology/variance-conventions.md, repo-internal):
- Measured matrix of every surface reaching the shared clustered-CR1
  denominator (linalg._compute_robust_vcov_numpy, vcov_type="hc1" only -
  a clustered call in any other family fails its row) and every tail-df
  convention passed to safe_inference/safe_inference_batch. The table is
  generated from expected literals in tests/test_variance_conventions.py;
  the parametrized tests assert those literals against live
  instrumentation (fast subset by default, full sweep under -m slow) and
  a byte-equality test keeps the committed doc in sync.
- Classifies the defects scheduled for the consolidation PRs - D1
  (absorb= vs fixed_effects=: 10.35% SE split on the same model, ratio
  exactly sqrt((360-2)/(360-66))), D2 (clustered CR1 never counts
  absorbed FE not nested in the cluster; anti-conservative, up to 5.51%
  at G=60), D4 (SunAbraham reports residual df per cohort-period cell but
  normal theory on aggregates) - and the legitimate exceptions L1-L4
  (StackedDiD CR1S, LPDiD G-1, influence-function paths, hc2/hc2_bm +
  survey TSL) with reasons and external anchors.

D3 fix (this PR's behavior change):
- New diff_diff.utils.absorbed_fe_rank: two-way absorbed-FE df from the
  connected components of the bipartite level graph (sum(levels) - C,
  minus 1 when the visible design carries an intercept column), exact per
  Abowd-Creecy-Kramarz. Levels and edges use positive-weight rows only,
  per the REGISTRY zero-weight-padding guarantee; NaN group keys raise an
  actionable ValueError. N>=3 keeps sum(levels-1) with the over-count
  limitation documented and tracked in TODO.md.
- Consumers: DiD(absorb=) and MultiPeriodDiD(absorb=) (intercept form,
  computed on the pre-transform frame), within-transform
  TwoWayFixedEffects (intercept form), SunAbraham (no-intercept form -
  its saturated design has coef_offset=0). Connected independent panels
  are bit-identical to the old count; disconnected and hierarchical
  designs (absorb=["state","state_year"]: true rank 29, old count 34) get
  corrected residual df and non-clustered classical/HC1 SE scaling, with
  the fail-closed NaN boundary moving consistently in both directions.
- External anchor: matches fixest::ssc(K.exact = TRUE) at <=1e-12 on a
  committed hierarchical R golden (benchmarks/data/fixest_kexact_golden.json,
  fixest 0.14.2, generator script included); fixest's DEFAULT
  K.exact=FALSE reproduces the old approximate count, recorded as a
  labeled deviation from the R default in REGISTRY.md.
- demean_by_groups/demean_by_group's second return value is now
  documented as the raw level count, NOT a valid df adjustment, pointing
  to absorbed_fe_rank (no callers consume it as df any more; returning
  the rank from the demeaner instead was rejected because within_transform
  delegates there and would compute the graph twice per TWFE/SA fit).
- REGISTRY.md: the TwoWayFixedEffects absorbed-FE df note is now
  component-aware, with cross-references from the DiD, MultiPeriodDiD,
  and SunAbraham sections; both sides of the D1 split and the Wooldridge
  hc1 mechanism status are documented. doc-deps.yaml registers the
  inventory doc against its 11 source modules; performance-plan.md
  records the helper's measured 1.9 ms / 186k-row cost with the
  codes-reuse optimization tracked as its own TODO row.

Author: igerber

CHANGELOG.md

@@ -8,6 +8,29 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Fixed
+- **Absorbed-FE degrees of freedom over-counted on disconnected and hierarchical
+  panels.** The df adjustment for absorbed fixed effects (`TwoWayFixedEffects`,
+  `SunAbraham`, `DifferenceInDifferences(absorb=)`, `MultiPeriodDiD(absorb=)`)
+  used `sum_d (levels_d − 1)`, which assumes the FE dimensions are mutually
+  independent and the unit×time incidence graph is connected. On disconnected
+  panels, and on nested specs like `absorb=["state", "state_year"]` (measured:
+  6 states × 5 years → true absorbed rank 29, old count 34), it over-stated the
+  absorbed rank, making the reported residual df too small and the non-clustered
+  classical/hc1 SEs too large. The new `diff_diff.utils.absorbed_fe_rank`
+  computes the two-way rank from the connected components of the bipartite level
+  graph; connected panels are bit-identical (`rtol=0, atol=1e-14` across all ten
+  estimator surfaces). Levels and connectivity are now evaluated over
+  positive-weight rows only, restoring the documented zero-weight-padding
+  inference-invariance on weighted fits with inert rows.
+- **New CI-enforced variance-convention audit matrix**
+  (`tests/test_variance_conventions.py` +
+  `docs/methodology/variance-conventions.md`): pins, per estimator surface, the
+  visible `k` reaching the shared clustered CR1 denominator and the tail-df
+  convention passed to inference, with each cell classified as a documented
+  defect (scheduled for the 3.9 variance-consolidation program) or a declared
+  legitimate exception with its reason. This is the map for the follow-up PRs
+  that converge the clustered CR1 `k` on the reghdfe/fixest nested convention
+  and the tail df on `t(G−1)`.
 - **`WooldridgeDiD` silently dropped genuine post-treatment effects
   (issue #724).** With `control_group="never_treated"`, the ETWFE design emitted
   every cohort×time indicator *including* the reference period — which is

DEFERRED.md

@@ -50,8 +50,8 @@ exists but parity can't be verified without a local toolchain.
 | Extend `WooldridgeDiD` `method ∈ {logit, poisson}` with `vcov_type ∈ {classical, hc2, hc2_bm}`: composing HC2 leverage + Bell-McCaffrey DOF with the QMLE pseudo-residual sandwich needs derivation + R parity vs `clubSandwich::vcovCR(glm, type="CR2")`. Rejected at `__init__`. | `wooldridge.py` | follow-up | Medium |
 | Multi-constraint CR2 parallel-trends test (AHT/HTZ) for `hc2_bm` fits: DiagnosticReport's PT check routes `vcov_type="hc2_bm"` sources to Bonferroni over the BM-adjusted per-row p-values because the generic chi-square joint Wald would discard the CR2 small-sample correction (see REPORTING.md "hc2_bm parallel-trends policy"). The proper joint test is the AHT/HTZ Wald with a Satterthwaite-style denominator df over the pre-period contrast block; needs derivation for the stacked/pooled WLS-CR2 layout + parity vs `clubSandwich::Wald_test(..., test="HTZ")`. | `diagnostic_report.py`, `linalg.py` | vcov/df round-trip PR | Low |
 | `PreTrendsPower` CS/SA `anticipation=1` R-parity fixture: R `pretrends` has no anticipation parameter, so the Python `_extract_pre_period_params` anticipation filter isn't R-parity-locked. Build a synthetic CS/SA result with `anticipation=1` and assert γ_p matches R's `slope_for_power()`. (Mechanism already covered by MC + full-VCV tests.) | `tests/test_methodology_pretrends.py`, `generate_pretrends_golden.R` | PR-C | Low |
-| Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. | `sun_abraham.py`, `linalg.py` | follow-up | Low |
-| Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. | `linalg.py`, `estimators.py` | SE-audit D4 | Low |
+| Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. Mechanism + measured inventory now in `docs/methodology/variance-conventions.md` (defect D2 family; 3.9 consolidation PR B). | `sun_abraham.py`, `linalg.py` | follow-up | Low |
+| Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. Closed form now derived and externally verified (Stata ~1e-15 / fixest ~1e-12) — see `docs/methodology/variance-conventions.md` defect D2; scheduled as PR B of the 3.9 consolidation program. | `linalg.py`, `estimators.py` | SE-audit D4 | Low |
 | Rust multiplier-bootstrap weight RNG (`generate_bootstrap_weights_batch`) seeds `Xoshiro256PlusPlus::seed_from_u64(seed+i)` per row; audit Python callers (`sdid.py`, `efficient_did_bootstrap.py`, `bootstrap_utils.py`) for parity-test gaps and, where a numpy-canonical equivalent exists, pre-generate in Python and pass through PyO3 (same fix shape as TROP RNG parity #354). | `rust/src/bootstrap.rs`, `bootstrap_utils.py` | follow-up | Medium |
 | `SyntheticDiD` bootstrap cross-language parity anchor vs R `synthdid::vcov(method="bootstrap")` or Julia `Synthdid.jl` (refit-native). Same-library validation is in place; Julia is the cleanest target. Tolerance ~1e-6 (BLAS+RNG paths preclude 1e-10). | `benchmarks/R/`, `benchmarks/julia/`, `tests/` | follow-up | Low |
 | CS R helpers hard-code `xformla = ~1`; no covariate-adjusted R benchmark for the IRLS path. | `tests/test_methodology_callaway.py` | #202 | Low |

TODO.md

@@ -21,9 +21,10 @@ Related tracking surfaces:
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
+| `absorbed_fe_rank` N>=3 general rank: the helper keeps `sum(levels-1)` for 3+ absorbed dims, exact for independent connected dims but an over-count for duplicated/nested triples (measured `a(5),b(4),c==b(4)`: true 7 vs formula 10) and for disconnected N-way graphs. Two-way is component-exact. Deriving general N-way FE rank is a hypergraph problem; do it with a reference (fixest's `fixef.rm`/reghdfe df_a) rather than a guess. See `docs/methodology/variance-conventions.md` D3. | `diff_diff/utils.py` | #variance-inventory | Mid | Low |
 | `SyntheticControl` conformal (CWZ 2021) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies — the residual-permutation shortcut is only valid for time-permutation-invariant proxies (SC/Lasso/DiD); an AR proxy needs innovation permutation. | `diff_diff/conformal.py`, `diff_diff/synthetic_control_results.py` | CWZ-2021 | Heavy | Low |
 | Make the post-fit `results.aggregate("event_study")` container consumable downstream. `EventStudyResults` is rejected by all THREE consumers that read a CS event study — `compute_honest_did` (`honest_did.py`, dispatches on `CallawaySantAnnaResults` and raises `TypeError`), `compute_pretrends_power` (`pretrends.py`, same), and `plot_event_study` (`visualization`, same) — so `fit(aggregate="event_study")` is still the only route for them and their error messages say so explicitly. Needs an `EventStudyResults` branch in each extraction path (consuming `event_time` / `is_reference` / `vcov` / `vcov_index` / per-row `df`) PLUS `base_period` and `anticipation` provenance, which the unified container does not carry and HonestDiD needs for its universal-base-period warning and pre-period classification. Gate with end-to-end tests: `compute_honest_did(res.aggregate("event_study"))` at `base_period="universal"`, and `compute_pretrends_power(...)` at `anticipation=1`. | `diff_diff/honest_did.py`, `diff_diff/pretrends.py`, `diff_diff/results_base.py` | #726 | Mid | Medium |
-| Derive and fix the `WooldridgeDiD` `hc1` SE gap vs Stata `jwdid`. Measured: every SE is uniformly SMALLER than `jwdid`'s, by 1.0280 at G=20 / 1.0132 at G=40 / 1.00264 at G=191 / 1.0010 at G=500 (ATT(g,t) points match exactly). The gap tracks `sqrt(G/(G-1))` but sits consistently above it, and `solve_ols` already applies the full CR1 `(G/(G-1))*((n-1)/(n-k))`, so a missing cluster factor is ruled out -- the likely source is the within-transform `k` accounting vs `hdfe`/`reghdfe`'s. Derive the exact factor FIRST; a `sqrt(G/(G-1))` patch would match only approximately and would force a loose tolerance on the parity test, defeating its purpose. Move the REGISTRY note and the pinned ratio in `tests/test_etwfe_cs_stata_parity.py` together with the fix. **Required artifact:** a committed subsample LADDER in the Stata golden (rosters = first N units per `first_treat` by sorted `countyreal`; rungs spanning G≈20..500, storing Stata `G`/`n`/`df_a`/`rank`/`df_r` and per-cell `att`/`se`), with parameterized ratio assertions. Only G=500 (full panel) and G=191 (all-eventually-treated arm) are pinned today, so the few-cluster behavior — where the gap is materially largest (~2.8%) — is ungated, and the ladder is also the instrument for comparing Stata's `df_a`/`rank` against the library's within-transform `k`, which is the leading hypothesis for the factor. | `diff_diff/wooldridge.py`, `diff_diff/linalg.py` | #723-followup | Mid | Medium |
+| Derive and fix the `WooldridgeDiD` `hc1` SE gap vs Stata `jwdid`. **MECHANISM NOW DERIVED** (see `docs/methodology/variance-conventions.md`, defect D2): the clustered CR1 denominator uses `k_visible` and never counts absorbed FE not nested in the cluster; the closed form `K_reference = explicit cols + (1 if no intercept col) + rank(non-nested FE | nested)` reproduces jwdid/reghdfe to ~1e-15 on three arms and fixest to ~1e-12 on two, retrodicting the G=20/G=40 rungs. Remaining work (PR B of the 3.9 consolidation program) is the threading + Rust port + Stata ladder, not derivation. Measured: every SE is uniformly SMALLER than `jwdid`'s, by 1.0280 at G=20 / 1.0132 at G=40 / 1.00264 at G=191 / 1.0010 at G=500 (ATT(g,t) points match exactly). The gap tracks `sqrt(G/(G-1))` but sits consistently above it, and `solve_ols` already applies the full CR1 `(G/(G-1))*((n-1)/(n-k))`, so a missing cluster factor is ruled out -- the likely source is the within-transform `k` accounting vs `hdfe`/`reghdfe`'s. Derive the exact factor FIRST; a `sqrt(G/(G-1))` patch would match only approximately and would force a loose tolerance on the parity test, defeating its purpose. Move the REGISTRY note and the pinned ratio in `tests/test_etwfe_cs_stata_parity.py` together with the fix. **Required artifact:** a committed subsample LADDER in the Stata golden (rosters = first N units per `first_treat` by sorted `countyreal`; rungs spanning G≈20..500, storing Stata `G`/`n`/`df_a`/`rank`/`df_r` and per-cell `att`/`se`), with parameterized ratio assertions. Only G=500 (full panel) and G=191 (all-eventually-treated arm) are pinned today, so the few-cluster behavior — where the gap is materially largest (~2.8%) — is ungated, and the ladder is also the instrument for comparing Stata's `df_a`/`rank` against the library's within-transform `k`, which is the leading hypothesis for the factor. | `diff_diff/wooldridge.py`, `diff_diff/linalg.py` | #723-followup | Mid | Medium |
 | `SunAbraham`: a cohort not observed at its own reference relative period (`e = -1 - anticipation`) makes that cohort's block collinear, so QR drops an unnamed column (`dropping 1 of 12 columns (column 9)`) and `overall_att` comes back **NaN**. Found by auditing the sibling estimator while fixing the ETWFE analogue (#724); PRE-EXISTING, not introduced there. Lower severity than #724 — that returned a silently WRONG finite number, this returns NaN with a rank warning — but the event-study surface still looks complete, so a user may not notice the loss. SA already omits its reference explicitly and tracks `_reference_observed`, so the fix is per-cohort support for that flag rather than the ETWFE-style redesign. | `diff_diff/sun_abraham.py` | #724-audit | Mid | Low |
 | Define `N_g` (W2025 Eqs. 7.4/7.6) for UNBALANCED panels where comparison-support filtering removes every observation of some units in an estimated cohort, then replace the fail-closed guard with the defined behavior. `_n_g_per_cohort` is read off the final sample, so those units vanish from the cohort-share weights; measured on a cohort supplied with 100 units of which 90 appear only at a dropped period, `aggregate(weights="cohort_share")` moves 1.8078 -> 3.8157. The paper assumes a balanced panel and does not say whether `N_g` counts the supplied cohort or the surviving units, and the two disagree materially, so `aggregate` currently raises naming the cohorts and counts ([M-125]); `weights="cell"` is unaffected and balanced panels never trip it. Settle the estimand (likely: count the supplied cohort, since ATT(g,t) is a cohort-level quantity, but that weights units with no retained observation) and gate with a test computing Eq. 7.4 by hand on unequal cohort sizes. | `diff_diff/wooldridge_results.py`, `diff_diff/wooldridge.py` | #729-followup | Mid | Medium |
 | `WooldridgeDiD` + `survey_design=` does not support DOMAIN ESTIMATION, so BOTH row-deleting paths are currently REFUSED (`NotImplementedError`, all three methods) rather than performed: unidentified-cohort exclusion ([M-123]) and comparison-support period filtering ([M-125]). One fix unblocks both. Implementing it properly means zero-padding the excluded rows' weights while retaining strata/PSU/FPC, per REGISTRY *Subpopulation Analysis (Phase 6)* / Lumley (2004) 3.4, so TSL variance and `df_survey = n_PSU - n_strata` use the full design (naive deletion measured 22 -> 14 on a two-stratum panel). `SurveyDesign.subpopulation()` already implements the contract and SpilloverDiD Wave E.3 is the in-repo precedent; the blocker is that the weighted within-transform rejects zero-weight units, shared machinery behind 7 estimators. Landing it would turn both refusals back into supported fits. Gate with a `SurveyDesign.subpopulation()` parity test on ATT, TSL SE and survey df where the excluded cohort exhausts a PSU. | `diff_diff/wooldridge.py`, `diff_diff/utils.py` | #724-codex-R4/R5 | Heavy | Medium |
@@ -48,6 +49,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
+| Reuse the demeaner's factorized codes in `absorbed_fe_rank` instead of re-factorizing: the helper adds ~1.9 ms per absorbed fit at 186k rows (7.7% of the fastest Rust-served TWFE fit; see `docs/performance-plan.md` "Component-aware absorbed-FE rank"), and both the helper and `demean_by_groups` factorize the same group columns. Threading the codes through the four call sites halves the factorize work; the `connected_components` call itself is ~1.1 ms. Deliberately not done in the correctness PR. | `diff_diff/utils.py` | #variance-inventory | Quick | Low |
 | `EfficientDiD` conditional path: the largest remaining O(n) stage is the sieve/nuisance construction outside the tiled pass (~9s at 10k). (The `_ridge_solve_weights` Python-prep shave landed 2026-07-07 — the `omega_stack[rest]` fancy-index copy and tail scatter are skipped when no row is zero-masked, byte-identical outputs; the `zero_mask` abs scan itself remains, needed for correctness.) | `efficient_did_covariates.py` | CS-scaling | Mid | Low |
 | `_rq_fit` LP assembly is dense (`A_eq = [X, I, -I]` with dense identity blocks, rebuilt per cell fit): a `scipy.sparse` construction would cut memory and likely HiGHS time for large cells / bootstrap-heavy covariate CiC/QDiD fits. CAVEAT before doing it: a different matrix representation can change HiGHS's vertex selection at degenerate/tied QR optima - end-to-end covariate goldens are tie-selection-gated (fine), but the `qr_cases` tight coefficient matches may shift to the equal-loss branch; re-run the parity suite and re-calibrate if needed. | `diff_diff/changes_in_changes.py::_rq_fit` | covariates PR | Quick | Low |
 | Evaluate flipping `DIFF_DIFF_SOLVE_OLS_FASTPATH` default-ON after an opt-in soak (the 2026-07 certified normal-equations Cholesky fast path, both backends). A flip needs: golden/parity-suite recapture at the tol-bounded posture (fitted ~1e-8 abs / SE ~1e-6 rel — the default today is byte-pinned in several benchmark conventions), certification-rate telemetry across real workloads (any decline is silent-correct but forfeits the speedup), and the staged default-flip protocol used for `df_convention` (v4-class change). Lifecycle tracked in docs/v4-deprecations.yaml (M-008). | `diff_diff/linalg.py::_resolve_solve_ols_fastpath`, `rust/src/linalg.rs::solve_ols_chol` | CS-scaling | Mid | Low |

benchmarks/R/generate_fixest_kexact_golden.R

@@ -0,0 +1,59 @@
+# Golden: fixest exact-vs-default FE counting on a hierarchical two-way design.
+#
+# absorb = [state, state_year] with state_year nested in state splits the
+# bipartite level graph into one component per state (C = 6), so the absorbed
+# dummy-space rank is 29 beyond the intercept -- not the naive
+# sum(levels - 1) = 34. fixest's default ssc(K.exact = FALSE) uses the naive
+# count (df.K = 36 here); ssc(K.exact = TRUE) computes the exact rank
+# (df.K = 31). diff-diff's component-aware absorbed_fe_rank matches the EXACT
+# side at machine precision (a documented deviation from the R *default*).
+#
+# Regenerate: Rscript benchmarks/R/generate_fixest_kexact_golden.R
+suppressMessages(library(fixest))
+suppressMessages(library(jsonlite))
+
+set.seed(7)
+d <- expand.grid(s = 0:5, y = 0:4, r = 1:4)
+d$state <- d$s
+d$state_year <- d$s * 100 + d$y
+d$x <- rnorm(nrow(d))
+d$out <- rnorm(nrow(d)) + 0.5 * d$x + 0.3 * d$s
+
+m_default <- feols(out ~ x | state + state_year, data = d, vcov = "iid")
+m_exact <- feols(out ~ x | state + state_year, data = d, vcov = "iid",
+                 ssc = ssc(K.exact = TRUE))
+
+golden <- list(
+  meta = list(
+    generator = "benchmarks/R/generate_fixest_kexact_golden.R",
+    r_version = paste(R.version$major, R.version$minor, sep = "."),
+    fixest_version = as.character(packageVersion("fixest")),
+    description = paste(
+      "Hierarchical two-way FE (state_year nested in state, C=6):",
+      "fixest default ssc(K.exact=FALSE) vs exact FE-rank counting.",
+      "diff-diff absorbed_fe_rank matches the K.exact=TRUE side."
+    )
+  ),
+  data = list(
+    state = d$state,
+    state_year = d$state_year,
+    x = d$x,
+    out = d$out
+  ),
+  n_obs = nrow(d),
+  coef = unname(coef(m_default)[["x"]]),
+  iid_default = list(
+    se = unname(se(m_default)[["x"]]),
+    df_k = degrees_freedom(m_default, "k")
+  ),
+  iid_k_exact = list(
+    se = unname(se(m_exact)[["x"]]),
+    df_k = degrees_freedom(m_exact, "k")
+  )
+)
+
+path <- "benchmarks/data/fixest_kexact_golden.json"
+write_json(golden, path, digits = NA, auto_unbox = TRUE, pretty = TRUE)
+cat("wrote", path, "\n")
+cat(sprintf("coef=%.15f default_se=%.15f exact_se=%.15f\n",
+            golden$coef, golden$iid_default$se, golden$iid_k_exact$se))