1"""F2-smoke-01: OLS slope of y on x with firm-clustered SE. Reads data/ only; writes out/results.json."""2import json, hashlib, os3import numpy as np, pandas as pd4import statsmodels.api as sm5from scipy import stats67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8path = os.path.join(ROOT, "data", "panel.csv")9raw = pd.read_csv(path)10n_raw = len(raw)1112# Exclusion steps (audit found none needed, but rules are applied and logged)13excl = []14d = raw.copy()15def step(name, rule, mask, why):16 global d17 nb = len(d); d = d[mask(d)].copy(); excl.append([name, rule, nb, len(d), nb - len(d), why])18step("1", "drop rows with missing firm_id/x/y", lambda z: z[["firm_id","x","y"]].notna().all(axis=1), "OLS requires complete cases")19step("2", "drop exact duplicate rows", lambda z: ~z.duplicated(), "guard against duplicated records")20pd.DataFrame(excl, columns=["step","rule","n_before","n_after","n_removed","rationale"]).to_csv(21 os.path.join(ROOT, "out", "exclusions.csv"), index=False)2223X = sm.add_constant(d["x"].to_numpy()); y = d["y"].to_numpy(); g = d["firm_id"].to_numpy()24ols = sm.OLS(y, X).fit()25cl = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": g, "use_correction": True, "df_correction": True})26G = int(len(np.unique(g))); N = len(d); K = X.shape[1]2728# manual CR1 check29e = ols.resid; XtXi = np.linalg.inv(X.T @ X)30meat = sum(np.outer(X[g == k].T @ e[g == k], X[g == k].T @ e[g == k]) for k in np.unique(g))31V = (G/(G-1))*((N-1)/(N-K)) * XtXi @ meat @ XtXi32assert abs(np.sqrt(V[1,1]) - cl.bse[1]) < 1e-103334beta, alpha = float(ols.params[1]), float(ols.params[0])35se_cl = float(cl.bse[1]); tcrit = float(stats.t.ppf(0.975, G - 1))36res = {37 "task_id": "F2-smoke-01",38 "beta_hat": beta, "alpha_hat": alpha,39 "se_ols": float(ols.bse[1]), "se_cluster": se_cl,40 "beta_ci_low": beta - tcrit*se_cl, "beta_ci_high": beta + tcrit*se_cl,41 "n_analysed": int(N), "n_clusters": G,42 "feasibility": "complete", "feasibility_reason": None,43 "data_issues": [],44 "n_input_rows": {"panel.csv": int(n_raw)},45 "random_seed": None,46 # supplementary47 "se_ratio_cluster_over_ols": se_cl/float(ols.bse[1]),48 "ci_critical_value_t_df": G - 1, "ci_critical_value": tcrit,49 "baseline_ci_low_ols": beta - float(stats.t.ppf(0.975, N-K))*float(ols.bse[1]),50 "baseline_ci_high_ols": beta + float(stats.t.ppf(0.975, N-K))*float(ols.bse[1]),51 "obs_per_firm_min": int(d.groupby("firm_id").size().min()),52 "obs_per_firm_max": int(d.groupby("firm_id").size().max()),53 "input_sha256": hashlib.sha256(open(path,"rb").read()).hexdigest(),54}55with open(os.path.join(ROOT, "out", "results.json"), "w") as f:56 json.dump(res, f, indent=2)57print(json.dumps(res, indent=2))58