1"""F2-smoke-01: OLS slope of y on x with firm-clustered SE. Reads data/panel.csv only."""2import hashlib, json, os3import numpy as np, pandas as pd4from scipy import stats5import statsmodels.api as sm67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8DATA = os.path.join(ROOT, "data", "panel.csv")9OUT = os.path.join(ROOT, "out")1011raw = pd.read_csv(DATA)12n_raw = len(raw)13issues, excl = [], []1415def step(name, rule, df_before, df_after, why):16 excl.append(dict(step=name, rule=rule, n_before=len(df_before), n_after=len(df_after),17 n_removed=len(df_before) - len(df_after), rationale=why))1819d0 = raw20d1 = d0.dropna(subset=["firm_id", "x", "y"])21step("1", "drop rows with missing firm_id/x/y", d0, d1, "OLS requires complete cases")22if len(d1) < len(d0):23 issues.append({"issue": f"{len(d0)-len(d1)} rows with missing values", "handling": "dropped"})24d2 = d1[np.isfinite(d1["x"]) & np.isfinite(d1["y"])]25step("2", "drop rows with non-finite x/y", d1, d2, "OLS requires finite values")26if len(d2) < len(d1):27 issues.append({"issue": f"{len(d1)-len(d2)} non-finite rows", "handling": "dropped"})28d3 = d2.drop_duplicates()29step("3", "drop exact duplicate rows", d2, d3, "exact duplicates treated as data-entry errors")30if len(d3) < len(d2):31 issues.append({"issue": f"{len(d2)-len(d3)} exact duplicate rows", "handling": "dropped"})32df = d3.reset_index(drop=True)3334X = sm.add_constant(df["x"].to_numpy())35y = df["y"].to_numpy()36g = df["firm_id"].to_numpy()37ols = sm.OLS(y, X).fit()38cl = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": g, "use_correction": True})39G = int(pd.Series(g).nunique())40beta, alpha = float(ols.params[1]), float(ols.params[0])41se_ols, se_cl = float(ols.bse[1]), float(cl.bse[1])42tcrit = float(stats.t.ppf(0.975, G - 1)) # t with G-1 df, standard for cluster-robust inference4344# manual CR1 cross-check45e = ols.resid46XtXi = np.linalg.inv(X.T @ X)47meat = sum(np.outer(X[g == k].T @ e[g == k], X[g == k].T @ e[g == k]) for k in np.unique(g))48N, K = X.shape49V = (G / (G - 1)) * ((N - 1) / (N - K)) * XtXi @ meat @ XtXi50assert abs(np.sqrt(V[1, 1]) - se_cl) < 1e-10 * max(1, se_cl)5152res = {53 "task_id": "F2-smoke-01",54 "beta_hat": beta,55 "alpha_hat": alpha,56 "se_ols": se_ols,57 "se_cluster": se_cl,58 "beta_ci_low": beta - tcrit * se_cl,59 "beta_ci_high": beta + tcrit * se_cl,60 "n_analysed": int(N),61 "n_clusters": G,62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": issues,65 "n_input_rows": {"panel.csv": int(n_raw)},66 "random_seed": None,67 "se_ratio_cluster_to_ols": se_cl / se_ols,68 "t_crit_cluster_df": tcrit,69 "cluster_df": G - 1,70 "baseline_ci_ols_low": beta - float(stats.t.ppf(0.975, N - K)) * se_ols,71 "baseline_ci_ols_high": beta + float(stats.t.ppf(0.975, N - K)) * se_ols,72 "obs_per_cluster_min": int(pd.Series(g).value_counts().min()),73 "obs_per_cluster_max": int(pd.Series(g).value_counts().max()),74 "input_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),75}76assert n_raw - sum(r["n_removed"] for r in excl) == N77os.makedirs(OUT, exist_ok=True)78pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)79with open(os.path.join(OUT, "results.json"), "w") as f:80 json.dump(res, f, indent=2)81print(json.dumps(res, indent=2))82