1"""F2-dev-01: two-way FE slope with firm-clustered SE. Run from workspace root.2No randomness is used (random_seed = null)."""3import json, hashlib4import numpy as np, pandas as pd5from scipy import stats6import statsmodels.formula.api as smf78PATH = "data/firm_panel.csv"9raw = pd.read_csv(PATH)10n_raw = len(raw)11issues = []12excl = []1314d = raw.copy()15# Step 1: missing values16nb = len(d); d = d.dropna(subset=["firm_id", "year", "x", "y"])17excl.append(("1", "drop rows with missing firm_id/year/x/y", nb, len(d), "FE estimator needs complete cases"))18# Step 2: duplicate firm-year keys19nb = len(d); d = d.drop_duplicates(subset=["firm_id", "year"], keep=False)20excl.append(("2", "drop all rows of duplicated firm_id-year keys", nb, len(d), "panel must have one row per firm-year"))21if n_raw - len(d) > 0:22 issues.append({"issue": f"{n_raw-len(d)} rows missing or duplicated", "handling": "dropped"})23# Step 3: singleton firms (no within-firm variation)24nb = len(d); cnt = d.groupby("firm_id")["y"].transform("size"); d = d[cnt > 1]25excl.append(("3", "drop singleton firms", nb, len(d), "singletons carry no within-firm information"))26d = d.sort_values(["firm_id", "year"]).reset_index(drop=True)2728N = len(d); G = d.firm_id.nunique(); T = d.year.nunique()29balanced = bool((d.groupby("firm_id").size() == T).all())3031# Baseline: pooled OLS with intercept32pooled = smf.ols("y ~ x", data=d).fit()33beta_pooled = float(pooled.params["x"])3435# Two-way FE via LSDV (exact)36lsdv = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit()37beta_fe = float(lsdv.params["x"])38se_conv = float(lsdv.bse["x"]) # homoskedastic, resid df = N - G - T39k_lsdv = int(lsdv.df_model) + 140resid = lsdv.resid.values4142# Cross-check with two-way within transformation (exact for balanced panel)43def dm(v):44 return v - v.groupby(d.firm_id).transform("mean") - v.groupby(d.year).transform("mean") + v.mean()45xt, yt = dm(d.x), dm(d.y)46beta_within = float((xt * yt).sum() / (xt ** 2).sum())4748# Firm-clustered sandwich for beta via FWL (x residualised on both FE sets).49xr = smf.ols("x ~ C(firm_id) + C(year)", data=d).fit().resid.values50scores = pd.Series(xr * resid).groupby(d.firm_id.values).sum().values51v_raw = float((scores ** 2).sum()) / float((xr ** 2).sum()) ** 252# Small-sample factor G/(G-1)*(N-1)/(N-K). Firm FE are nested within clusters and53# are not counted in K (Stata xtreg-fe / reghdfe convention): K = slope + (T-1) year effects.54K_nested = 1 + (T - 1)55c_nested = G / (G - 1) * (N - 1) / (N - K_nested)56se_cluster = float(np.sqrt(c_nested * v_raw))57# Sensitivity: all LSDV dummies counted in K (statsmodels default on LSDV)58c_lsdv = G / (G - 1) * (N - 1) / (N - k_lsdv)59se_cluster_lsdv = float(np.sqrt(c_lsdv * v_raw))60sm_cl = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit(cov_type="cluster", cov_kwds={"groups": d.firm_id.values})61se_cluster_sm = float(sm_cl.bse["x"])6263tcrit = float(stats.t.ppf(0.975, G - 1))64lo, hi = beta_fe - tcrit * se_cluster, beta_fe + tcrit * se_cluster6566res = {67 "task_id": "F2-dev-01",68 "beta_pooled": beta_pooled,69 "beta_fe": beta_fe,70 "se_conventional": se_conv,71 "se_cluster": se_cluster,72 "beta_ci_low": float(lo),73 "beta_ci_high": float(hi),74 "n_firms": int(G),75 "n_years": int(T),76 "n_analysed": int(N),77 "feasibility": "complete",78 "feasibility_reason": None,79 "data_issues": issues,80 "n_input_rows": {"firm_panel.csv": int(n_raw)},81 "random_seed": None,82 # supplementary keys83 "se_pooled_conventional": float(pooled.bse["x"]),84 "beta_fe_within_check": beta_within,85 "se_cluster_lsdv_all_dummies_counted": se_cluster_lsdv,86 "se_cluster_statsmodels_lsdv": se_cluster_sm,87 "se_ratio_cluster_to_conventional": se_cluster / se_conv,88 "cluster_small_sample_factor": c_nested,89 "ci_critical_value_t": tcrit,90 "ci_df": int(G - 1),91 "panel_balanced": balanced,92 "input_sha256": hashlib.sha256(open(PATH, "rb").read()).hexdigest(),93}94with open("out/results.json", "w") as f:95 json.dump(res, f, indent=2)96pd.DataFrame([dict(step=s, rule=r, n_before=a, n_after=b, n_removed=a - b, rationale=w) for s, r, a, b, w in excl]97 ).to_csv("out/exclusions.csv", index=False)98print(json.dumps(res, indent=2))99