1"""F2-smoke-01: OLS slope of y on x with firm-clustered SE. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4from scipy import stats56d_raw = pd.read_csv("data/panel.csv")7n_raw = len(d_raw)8d = d_raw.dropna(subset=["firm_id", "x", "y"])9issues = []10if len(d) < n_raw:11 issues.append({"issue": f"{n_raw-len(d)} rows with missing values", "handling": "dropped (complete-case)"})12if d.duplicated().sum() > 0:13 issues.append({"issue": "exact duplicate rows", "handling": "retained (no basis to treat as errors)"})1415y = d["y"].to_numpy(float); x = d["x"].to_numpy(float); g = d["firm_id"].to_numpy()16n = len(y); X = np.column_stack([np.ones(n), x]); k = 217XtX_inv = np.linalg.inv(X.T @ X)18b = XtX_inv @ X.T @ y19e = y - X @ b20# conventional homoskedastic SE (naive baseline)21s2 = e @ e / (n - k)22se_ols = float(np.sqrt(s2 * XtX_inv[1, 1]))23# cluster-robust (CR1, Stata-style small-sample factor), clustered by firm_id24groups = np.unique(g); G = len(groups)25meat = np.zeros((k, k))26for gg in groups:27 m = g == gg28 s = X[m].T @ e[m]29 meat += np.outer(s, s)30c = G / (G - 1) * (n - 1) / (n - k)31V = c * XtX_inv @ meat @ XtX_inv32se_cl = float(np.sqrt(V[1, 1]))33tcrit = float(stats.t.ppf(0.975, G - 1))34res = {35 "task_id": "F2-smoke-01",36 "beta_hat": float(b[1]), "alpha_hat": float(b[0]),37 "se_ols": se_ols, "se_cluster": se_cl,38 "beta_ci_low": float(b[1] - tcrit * se_cl), "beta_ci_high": float(b[1] + tcrit * se_cl),39 "n_analysed": int(n), "n_clusters": int(G),40 "feasibility": "complete", "feasibility_reason": None,41 "data_issues": issues, "n_input_rows": {"panel.csv": int(n_raw)},42 "random_seed": None,43 # supplementary44 "ci_t_critical_df_Gminus1": tcrit,45 "se_ratio_cluster_over_ols": se_cl / se_ols,46 "baseline_ci_ols_low": float(b[1] - stats.t.ppf(0.975, n - k) * se_ols),47 "baseline_ci_ols_high": float(b[1] + stats.t.ppf(0.975, n - k) * se_ols),48}49os.makedirs("out", exist_ok=True)50with open("out/results.json", "w") as f:51 json.dump(res, f, indent=2)52print(json.dumps(res, indent=2))53