1"""F2-dev-01: two-way FE slope with firm-clustered SE. Run from workspace root."""2import json3import numpy as np4import pandas as pd5import statsmodels.formula.api as smf6from scipy import stats78d0 = pd.read_csv("data/firm_panel.csv")9n_input = len(d0)10issues = []11d = d0.copy()12na = d[["firm_id", "year", "x", "y"]].isna().any(axis=1)13if na.any():14 issues.append({"issue": f"{int(na.sum())} rows with missing values", "handling": "dropped"})15 d = d[~na]16dup = d.duplicated(["firm_id", "year"], keep=False)17if dup.any():18 issues.append({"issue": f"{int(dup.sum())} rows with duplicated firm-year keys", "handling": "exact duplicates dropped"})19 d = d.drop_duplicates()20sizes = d.groupby("firm_id").size()21balanced = sizes.nunique() == 1 and sizes.iloc[0] == d.year.nunique()22if not balanced:23 issues.append({"issue": "panel is unbalanced", "handling": "kept all observations; FE estimator valid for unbalanced panels"})24within_sd = d.groupby("firm_id").x.std()25if (within_sd == 0).any():26 issues.append({"issue": f"{int((within_sd == 0).sum())} firms with no within variation in x", "handling": "retained (contribute nothing to beta)"})2728d = d.sort_values(["firm_id", "year"]).reset_index(drop=True)29pooled = smf.ols("y ~ x", data=d).fit()30fe = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit()31G = d.firm_id.nunique(); T = d.year.nunique(); N = len(d)32# Cluster-robust (CR1) with statsmodels default small-sample factor G/(G-1)*(N-1)/(N-K), K incl. all dummies33fe_cl = smf.ols("y ~ x + C(firm_id) + C(year)", data=d).fit(34 cov_type="cluster", cov_kwds={"groups": d.firm_id.values, "use_correction": True, "df_correction": True})35beta = float(fe.params["x"])36se_conv = float(fe.bse["x"])37se_cl = float(fe_cl.bse["x"])38tcrit = stats.t.ppf(0.975, G - 1)3940# Sensitivity: within-transformed regression (K=1 slope; like xtreg fe), CR1 with G/(G-1)*(N-1)/(N-1)41dm = d.copy()42for v in ["x", "y"]:43 dm[v + "_w"] = dm[v] - dm.groupby("firm_id")[v].transform("mean") - dm.groupby("year")[v].transform("mean") + dm[v].mean()44e = dm.y_w - beta * dm.x_w45sxx = (dm.x_w ** 2).sum()46score = (dm.x_w * e).groupby(dm.firm_id).sum()47v0 = (score ** 2).sum() / sxx ** 248se_cl_nocorr = float(np.sqrt(v0))49se_cl_xtreg = float(np.sqrt(v0 * G / (G - 1) * (N - 1) / (N - 2)))5051res = {52 "task_id": "F2-dev-01",53 "beta_pooled": float(pooled.params["x"]),54 "beta_fe": beta,55 "se_conventional": se_conv,56 "se_cluster": se_cl,57 "beta_ci_low": beta - tcrit * se_cl,58 "beta_ci_high": beta + tcrit * se_cl,59 "n_firms": int(G),60 "n_years": int(T),61 "n_analysed": int(N),62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": issues,65 "n_input_rows": {"firm_panel.csv": int(n_input)},66 "random_seed": None,67 "sensitivity": {68 "se_pooled_conventional": float(pooled.bse["x"]),69 "se_cluster_no_small_sample_correction": se_cl_nocorr,70 "se_cluster_K_excluding_absorbed_FE": se_cl_xtreg,71 "t_critical_df_G_minus_1": float(tcrit),72 "ci_normal_low": beta - 1.959964 * se_cl,73 "ci_normal_high": beta + 1.959964 * se_cl,74 "fe_resid_df": float(fe.df_resid),75 },76}77with open("out/results.json", "w") as f:78 json.dump(res, f, indent=2)79print(json.dumps(res, indent=2))80