1"""F2-smoke-01: OLS slope of y on x with firm-clustered SEs. Reads only data/panel.csv."""2import json, os3import numpy as np, pandas as pd4from scipy import stats56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7OUT = os.path.join(ROOT, "out")8df_raw = pd.read_csv(os.path.join(ROOT, "data", "panel.csv"))9n_raw = len(df_raw)1011# Exclusion ledger: drop rows with missing / non-finite values or missing firm_id12excl = []13df = df_raw.copy()14mask = np.isfinite(df[["x", "y"]]).all(axis=1) & df["firm_id"].notna()15excl.append(("1", "drop rows with missing/non-finite x, y or firm_id", len(df), int(mask.sum()),16 int((~mask).sum()), "OLS requires complete finite observations"))17df = df[mask]18pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(19 os.path.join(OUT, "exclusions.csv"), index=False)2021data_issues = []22if (~mask).sum() > 0:23 data_issues.append({"issue": f"{int((~mask).sum())} incomplete rows", "handling": "dropped"})24if df_raw.duplicated().sum() > 0:25 data_issues.append({"issue": f"{int(df_raw.duplicated().sum())} duplicate rows", "handling": "retained"})2627y = df["y"].to_numpy(float)28X = np.column_stack([np.ones(len(df)), df["x"].to_numpy(float)])29g = df["firm_id"].to_numpy()30N, K = X.shape31XtX_inv = np.linalg.inv(X.T @ X)32b = XtX_inv @ X.T @ y33e = y - X @ b # residuals3435# Baseline: conventional homoskedastic OLS SE36s2 = e @ e / (N - K)37se_ols = float(np.sqrt(s2 * XtX_inv[1, 1]))3839# Cluster-robust (CR1, Stata-style small-sample factor) by firm_id40groups = np.unique(g)41G = len(groups)42meat = np.zeros((K, K))43for gid in groups:44 idx = g == gid45 sg = X[idx].T @ e[idx]46 meat += np.outer(sg, sg)47c = G / (G - 1) * (N - 1) / (N - K)48V = c * XtX_inv @ meat @ XtX_inv49se_cl = float(np.sqrt(V[1, 1]))50tcrit = float(stats.t.ppf(0.975, G - 1))5152res = {53 "task_id": "F2-smoke-01",54 "beta_hat": float(b[1]),55 "alpha_hat": float(b[0]),56 "se_ols": se_ols,57 "se_cluster": se_cl,58 "beta_ci_low": float(b[1] - tcrit * se_cl),59 "beta_ci_high": float(b[1] + tcrit * se_cl),60 "n_analysed": int(N),61 "n_clusters": int(G),62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": data_issues,65 "n_input_rows": {"panel.csv": int(n_raw)},66 "random_seed": None,67 # supplementary (not required)68 "se_ratio_cluster_over_ols": se_cl / se_ols,69 "ci_t_critical_df_G_minus_1": tcrit,70 "se_cluster_alpha": float(np.sqrt(V[0, 0])),71 "baseline_ci_low_ols": float(b[1] - stats.t.ppf(0.975, N - K) * se_ols),72 "baseline_ci_high_ols": float(b[1] + stats.t.ppf(0.975, N - K) * se_ols),73 "share_between_firm_var_resid": float(pd.Series(e).groupby(g).mean().var() / pd.Series(e).var()),74 "share_between_firm_var_x": float(df.groupby("firm_id")["x"].mean().var() / df["x"].var()),75}76with open(os.path.join(OUT, "results.json"), "w") as f:77 json.dump(res, f, indent=2)78print(json.dumps(res, indent=2))79