1"""F2-smoke-01: OLS slope of y on x with firm-clustered SE. Run from workspace root."""2import hashlib, json, os3import numpy as np, pandas as pd4from scipy import stats5import statsmodels.api as sm67DATA = "data/panel.csv"8raw = pd.read_csv(DATA)9n_raw = len(raw)10issues, excl = [], []1112# Exclusion steps (documented even if they remove nothing)13d = raw.copy()14n0 = len(d); d = d.dropna(subset=["firm_id", "x", "y"])15excl.append(("1", "drop rows with missing firm_id/x/y", n0, len(d), "OLS needs complete cases"))16n0 = len(d); d = d[np.isfinite(d[["x", "y"]]).all(axis=1)]17excl.append(("2", "drop rows with non-finite x/y", n0, len(d), "OLS needs finite values"))18ndup = int(d.duplicated().sum())19if ndup:20 issues.append({"issue": f"{ndup} exact duplicate rows", "handling": "kept (cannot distinguish from genuine repeats)"})21for s, r, a, b, why in excl:22 if a != b:23 issues.append({"issue": f"{a-b} rows removed by rule: {r}", "handling": why})2425X = sm.add_constant(d["x"].to_numpy())26yv = d["y"].to_numpy()27g = d["firm_id"].to_numpy()28ols = sm.OLS(yv, X).fit() # conventional homoskedastic SE (baseline)29cl = sm.OLS(yv, X).fit(cov_type="cluster", cov_kwds={"groups": g, "use_correction": True}) # CR130G = int(len(np.unique(g)))31beta, alpha = float(cl.params[1]), float(cl.params[0])32se_c = float(cl.bse[1])33tcrit = float(stats.t.ppf(0.975, G - 1)) # small-G t(G-1) reference distribution34res = {35 "task_id": "F2-smoke-01",36 "beta_hat": beta,37 "alpha_hat": alpha,38 "se_ols": float(ols.bse[1]),39 "se_cluster": se_c,40 "beta_ci_low": beta - tcrit * se_c,41 "beta_ci_high": beta + tcrit * se_c,42 "n_analysed": int(len(d)),43 "n_clusters": G,44 "feasibility": "complete",45 "feasibility_reason": None,46 "data_issues": issues,47 "n_input_rows": {"panel.csv": int(n_raw)},48 "random_seed": None,49 # supplementary50 "t_crit_cluster": tcrit,51 "ci_df": G - 1,52 "se_ratio_cluster_over_ols": se_c / float(ols.bse[1]),53 "baseline_ci_low_ols": float(ols.conf_int()[1][0]),54 "baseline_ci_high_ols": float(ols.conf_int()[1][1]),55 "obs_per_cluster_min": int(d.groupby("firm_id").size().min()),56 "obs_per_cluster_max": int(d.groupby("firm_id").size().max()),57 "input_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),58}59os.makedirs("out", exist_ok=True)60with open("out/results.json", "w") as f:61 json.dump(res, f, indent=2)62pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"]).assign(63 n_removed=lambda t: t.n_before - t.n_after)[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]].to_csv("out/exclusions.csv", index=False)64print(json.dumps(res, indent=2))65