1"""F2-smoke-01: OLS slope of y on x with firm-clustered SE. Reads data/panel.csv only."""2import json, os3import numpy as np, pandas as pd4from scipy import stats56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7raw = pd.read_csv(os.path.join(ROOT, "data", "panel.csv"))8n_raw = len(raw)910# Exclusion ledger (rules applied in order)11ledger, d, issues = [], raw.copy(), []12def step(name, rule, mask, why):13 global d14 nb = len(d); d = d[mask(d)]; ledger.append(dict(step=name, rule=rule, n_before=nb, n_after=len(d), n_removed=nb - len(d), rationale=why))15step("1", "drop rows with missing firm_id/x/y", lambda z: z[["firm_id","x","y"]].notna().all(axis=1), "OLS needs complete cases")16step("2", "drop rows with non-finite x/y", lambda z: np.isfinite(z["x"]) & np.isfinite(z["y"]), "OLS needs finite values")17step("3", "drop exact duplicate rows", lambda z: ~z.duplicated(), "guard against duplicated records")18for e in ledger:19 if e["n_removed"] > 0:20 issues.append({"issue": f"{e['n_removed']} rows failed rule '{e['rule']}'", "handling": "excluded"})21pd.DataFrame(ledger).to_csv(os.path.join(ROOT, "out", "exclusions.csv"), index=False)2223d = d.sort_values(["firm_id", "x", "y"]).reset_index(drop=True) # deterministic order24y = d["y"].to_numpy(float); X = np.column_stack([np.ones(len(d)), d["x"].to_numpy(float)])25g = d["firm_id"].to_numpy()26n, k = X.shape27XtX_inv = np.linalg.inv(X.T @ X)28b = XtX_inv @ X.T @ y29e = y - X @ b30# Conventional homoskedastic SE31s2 = e @ e / (n - k)32se_ols = float(np.sqrt(s2 * XtX_inv[1, 1]))33# Cluster-robust (CR1, Stata-style small-sample factor) SE34groups = np.unique(g); G = len(groups)35meat = np.zeros((k, k))36for gg in groups:37 m = g == gg38 sg = X[m].T @ e[m]39 meat += np.outer(sg, sg)40c = G / (G - 1) * (n - 1) / (n - k)41V = c * XtX_inv @ meat @ XtX_inv42se_cl = float(np.sqrt(V[1, 1]))43tcrit = stats.t.ppf(0.975, G - 1)4445res = {46 "task_id": "F2-smoke-01",47 "beta_hat": float(b[1]), "alpha_hat": float(b[0]),48 "se_ols": se_ols, "se_cluster": se_cl,49 "beta_ci_low": float(b[1] - tcrit * se_cl), "beta_ci_high": float(b[1] + tcrit * se_cl),50 "n_analysed": int(n), "n_clusters": int(G),51 "feasibility": "complete", "feasibility_reason": None,52 "data_issues": issues,53 "n_input_rows": {"panel.csv": int(n_raw)},54 "random_seed": None,55 "ci_critical_value_t_df": int(G - 1), "ci_critical_value": float(tcrit),56 "cluster_se_correction": "CR1: G/(G-1)*(N-1)/(N-K)",57 "se_ratio_cluster_to_ols": se_cl / se_ols,58 "obs_per_cluster_min": int(d.groupby("firm_id").size().min()),59 "obs_per_cluster_max": int(d.groupby("firm_id").size().max()),60}61with open(os.path.join(ROOT, "out", "results.json"), "w") as f:62 json.dump(res, f, indent=2)63print(json.dumps(res, indent=2))64