1"""F2-dev-01: two-way FE slope with firm-clustered SE. Run from workspace root."""2import json, hashlib3import numpy as np, pandas as pd4from scipy import stats5import statsmodels.formula.api as smf67DATA = "data/firm_panel.csv"8raw = pd.read_csv(DATA)9n_raw = len(raw)10issues, excl = [], []1112def step(name, rule, before, after, why):13 excl.append(dict(step=name, rule=rule, n_before=before, n_after=after,14 n_removed=before - after, rationale=why))1516d = raw.copy()17b = len(d); d = d.dropna(subset=["firm_id", "year", "x", "y"])18step(1, "drop rows with missing firm_id/year/x/y", b, len(d), "FE estimator needs complete cases")19b = len(d); d = d.drop_duplicates()20step(2, "drop exact duplicate rows", b, len(d), "duplicates carry no new information")21b = len(d); d = d[np.isfinite(d.x) & np.isfinite(d.y)]22step(3, "drop non-finite x/y", b, len(d), "non-numeric values invalid")23if d.duplicated(["firm_id", "year"]).any():24 issues.append({"issue": "conflicting duplicate firm-year keys", "handling": "kept; reported"})25balanced = d.groupby("firm_id").size().nunique() == 1 and d.groupby("year").size().nunique() == 126if not balanced:27 issues.append({"issue": "panel unbalanced", "handling": "LSDV estimation handles imbalance"})2829N = len(d); G = d.firm_id.nunique(); T = d.year.nunique()3031# Baseline: pooled OLS32pooled = smf.ols("y ~ x", data=d).fit()33beta_pooled = float(pooled.params["x"])3435# Two-way FE via LSDV (exact, handles imbalance)36fe = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit()37beta_fe = float(fe.params["x"])38se_conv = float(fe.bse["x"])39K_all = int(fe.df_model) + 1 # all estimated params incl. FE4041# Firm-clustered (CR1) SE via Frisch-Waugh: residualise x on FE, use FE residuals.42fwl = smf.ols("x ~ C(firm_id) + C(year)", data=d).fit()43xt = fwl.resid.values; u = fe.resid.values44# check FWL45assert abs((xt @ d.y.values) / (xt @ xt) - beta_fe) < 1e-846scores = pd.Series(xt * u).groupby(d.firm_id.values).sum().values47V0 = (scores @ scores) / (xt @ xt) ** 248# Small-sample factor: firm FE are nested within clusters, so (as in Stata xtreg/reghdfe)49# they are not counted in K; K = slope + intercept + (T-1) year effects.50K_nest = 1 + 1 + (T - 1)51c = G / (G - 1) * (N - 1) / (N - K_nest)52se_cluster = float(np.sqrt(c * V0))53# Sensitivity: counting all FE parameters (statsmodels LSDV default)54c_all = G / (G - 1) * (N - 1) / (N - K_all)55se_cluster_allK = float(np.sqrt(c_all * V0))56lsdv_cl = fe.get_robustcov_results(cov_type="cluster", groups=d.firm_id.values)57se_cluster_sm = float(lsdv_cl.bse[list(fe.params.index).index("x")])5859tcrit = float(stats.t.ppf(0.975, G - 1))60lo, hi = beta_fe - tcrit * se_cluster, beta_fe + tcrit * se_cluster6162res = {63 "task_id": "F2-dev-01",64 "beta_pooled": beta_pooled,65 "beta_fe": beta_fe,66 "se_conventional": se_conv,67 "se_cluster": se_cluster,68 "beta_ci_low": float(lo),69 "beta_ci_high": float(hi),70 "n_firms": int(G),71 "n_years": int(T),72 "n_analysed": int(N),73 "feasibility": "complete",74 "feasibility_reason": None,75 "data_issues": issues,76 "n_input_rows": {"firm_panel.csv": int(n_raw)},77 "random_seed": None,78 # supplementary79 "baseline_se_pooled_conventional": float(pooled.bse["x"]),80 "se_cluster_unadjusted_CR0": float(np.sqrt(V0)),81 "se_cluster_allK_dfadj": se_cluster_allK,82 "se_cluster_statsmodels_lsdv": se_cluster_sm,83 "ci_t_critical_df_G_minus_1": tcrit,84 "cluster_small_sample_factor": float(c),85 "df_resid_conventional": int(fe.df_resid),86 "panel_balanced": bool(balanced),87 "input_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),88}89with open("out/results.json", "w") as f:90 json.dump(res, f, indent=2)91pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)92print(json.dumps(res, indent=2))93