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 stats56PATH = "data/firm_panel.csv"7raw = pd.read_csv(PATH)8n_raw = len(raw)9issues = []10excl = []1112# --- audit / exclusions (each step logged even if nothing removed)13d = raw.copy()14def step(name, rule, mask, why):15 global d16 nb = len(d); d = d[mask(d)].copy()17 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=len(d), n_removed=nb-len(d), rationale=why))18step("1", "drop rows with missing/non-finite firm_id, year, x or y",19 lambda t: np.isfinite(t[["firm_id","year","x","y"]].astype(float)).all(axis=1), "FE estimator needs complete cases")20step("2", "drop exact duplicate rows", lambda t: ~t.duplicated(), "duplicates would double-weight observations")21step("3", "drop duplicated firm_id-year keys", lambda t: ~t.duplicated(["firm_id","year"], keep=False), "panel must have one row per firm-year")22pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)23counts = d.groupby("firm_id").size()24if counts.nunique() != 1 or counts.iloc[0] != d.year.nunique():25 issues.append({"issue": "unbalanced panel", "handling": "FE estimator applied to available obs"})26d = d.sort_values(["firm_id","year"]).reset_index(drop=True)2728# --- baseline: pooled OLS with intercept29X0 = np.column_stack([np.ones(len(d)), d.x.values])30beta_pooled = float(np.linalg.lstsq(X0, d.y.values, rcond=None)[0][1])3132# --- two-way FE via LSDV (firm + year dummies), exact for any panel33firm = pd.get_dummies(d.firm_id, drop_first=False).values.astype(float)34yr = pd.get_dummies(d.year, drop_first=True).values.astype(float)35X = np.column_stack([d.x.values, firm, yr])36y = d.y.values37N, K = X.shape38XtX_inv = np.linalg.inv(X.T @ X)39b = XtX_inv @ X.T @ y40u = y - X @ b41beta_fe = float(b[0])42G = d.firm_id.nunique(); T = d.year.nunique()43rank = np.linalg.matrix_rank(X)44assert rank == K4546# conventional SE: s^2 with N-K dof (K counts all firm and year effects)47s2 = u @ u / (N - K)48se_conv = float(np.sqrt(s2 * XtX_inv[0, 0]))4950# cluster-robust (CR1) by firm. Only the x-row of the sandwich is needed.51# Firm effects are nested within clusters, so (as in reghdfe / xtreg,fe) they52# are not counted in the small-sample K: K_cr = 1 (x) + (T-1) year dummies.53a = XtX_inv[0, :] @ X.T # influence weights for beta_fe54sc = pd.Series(a * u).groupby(d.firm_id.values).sum().values55V0 = float(sc @ sc)56K_cr = 1 + (T - 1)57c_nested = G / (G - 1) * (N - 1) / (N - K_cr)58c_full = G / (G - 1) * (N - 1) / (N - K)59se_cluster = float(np.sqrt(c_nested * V0))60se_cluster_allK = float(np.sqrt(c_full * V0))61se_cluster_noadj = float(np.sqrt(V0))62tcrit = float(stats.t.ppf(0.975, G - 1))63lo, hi = beta_fe - tcrit * se_cluster, beta_fe + tcrit * se_cluster6465res = {66 "task_id": "F2-dev-01",67 "beta_pooled": beta_pooled,68 "beta_fe": beta_fe,69 "se_conventional": se_conv,70 "se_cluster": se_cluster,71 "beta_ci_low": float(lo),72 "beta_ci_high": float(hi),73 "n_firms": int(G),74 "n_years": int(T),75 "n_analysed": int(N),76 "feasibility": "complete",77 "feasibility_reason": None,78 "data_issues": issues,79 "n_input_rows": {"firm_panel.csv": int(n_raw)},80 "random_seed": None,81 # supplementary (not required)82 "ci_critical_value_t_Gminus1": tcrit,83 "ci_dof": int(G - 1),84 "se_cluster_ratio_to_conventional": se_cluster / se_conv,85 "se_cluster_alt_K_all_dummies": se_cluster_allK,86 "se_cluster_alt_no_small_sample_adj": se_cluster_noadj,87 "ci_normal_low": beta_fe - 1.959963984540054 * se_cluster,88 "ci_normal_high": beta_fe + 1.959963984540054 * se_cluster,89 "n_parameters_lsdv": int(K),90 "input_sha256": hashlib.sha256(open(PATH, "rb").read()).hexdigest(),91}92with open("out/results.json", "w") as f:93 json.dump(res, f, indent=2)94print(json.dumps({k: res[k] for k in ["beta_pooled","beta_fe","se_conventional","se_cluster","beta_ci_low","beta_ci_high"]}, indent=2))95