1"""F2-dev-01: TWFE slope with firm-clustered SE. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4import statsmodels.formula.api as smf5from scipy import stats67os.makedirs("out", exist_ok=True)8raw = pd.read_csv("data/firm_panel.csv")9n_raw = len(raw)10issues, excl = [], []11d = raw.copy()1213def step(name, rule, mask, why):14 global d15 nb = len(d); d = d[~mask].copy()16 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=len(d), n_removed=nb - len(d), rationale=why))1718step("1", "drop rows with missing firm_id/year/x/y", d[["firm_id","year","x","y"]].isna().any(axis=1), "cannot be used in regression")19step("2", "drop exact duplicate rows", d.duplicated(), "duplicate records")20dupk = d.duplicated(["firm_id","year"], keep=False)21if dupk.any():22 issues.append({"issue": f"{int(dupk.sum())} rows with conflicting duplicate firm-year keys", "handling": "dropped"})23step("3", "drop conflicting duplicate firm-year keys", dupk, "ambiguous observation")24for e in excl:25 if e["n_removed"] and e["step"] != "3":26 issues.append({"issue": f"{e['n_removed']} rows removed by rule: {e['rule']}", "handling": "dropped"})27d = d.sort_values(["firm_id","year"]).reset_index(drop=True)28pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)2930# Baseline: pooled OLS31pooled = smf.ols("y ~ x", data=d).fit()32# TWFE via LSDV33fe = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit()34G = d.firm_id.nunique(); T = d.year.nunique(); N = len(d)35# Cluster-robust (CR1, small-sample factor G/(G-1)*(N-1)/(N-K), K counts all LSDV params)36fec = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit(37 cov_type="cluster", cov_kwds={"groups": d.firm_id.values, "use_correction": True, "df_correction": True})38b = float(fe.params["x"]); se_conv = float(fe.bse["x"]); se_cl = float(fec.bse["x"])39# Sensitivity: CR1 not counting firm dummies (nested in clusters; xtreg-fe convention, K = 1 + (T-1) + 1 for x, year dummies, const)40fec0 = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit(cov_type="cluster", cov_kwds={"groups": d.firm_id.values, "use_correction": False})41K_nf = 1 + (T - 1) + 142se_cl_nested = float(fec0.bse["x"]) * np.sqrt(G/(G-1) * (N-1)/(N-K_nf))43tcrit = float(stats.t.ppf(0.975, G - 1))4445# Cross-check by within-transformation (balanced panel => two-way demeaning exact)46def dm(v):47 return v - d.groupby("firm_id")[v.name].transform("mean") - d.groupby("year")[v.name].transform("mean") + v.mean()48xt, yt = dm(d.x), dm(d.y)49b_w = float((xt*yt).sum()/(xt**2).sum())50assert abs(b_w - b) < 1e-8, (b_w, b)5152res = {53 "task_id": "F2-dev-01",54 "beta_pooled": float(pooled.params["x"]),55 "beta_fe": b,56 "se_conventional": se_conv,57 "se_cluster": se_cl,58 "beta_ci_low": b - tcrit*se_cl,59 "beta_ci_high": b + tcrit*se_cl,60 "n_firms": int(G), "n_years": int(T), "n_analysed": int(N),61 "feasibility": "complete", "feasibility_reason": None,62 "data_issues": issues,63 "n_input_rows": {"firm_panel.csv": int(n_raw)},64 "random_seed": None,65 "baseline_se_pooled_conventional": float(pooled.bse["x"]),66 "ci_critical_value_t_Gminus1": tcrit,67 "se_ratio_cluster_to_conventional": se_cl/se_conv,68 "beta_fe_within_check": b_w,69 "se_cluster_sensitivity_firmFE_not_counted": se_cl_nested,70 "se_cluster_CR0_no_correction": float(fec0.bse["x"]),71 "n_params_lsdv": int(len(fe.params)),72}73with open("out/results.json", "w") as f:74 json.dump(res, f, indent=2)75print(json.dumps(res, indent=2))76