1"""Replication entry point for F6-conf-04. Run from workspace root."""2import json, os3import numpy as np4import pandas as pd5import statsmodels.api as sm6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910SEED = None # no randomness used11OUT = "out"12os.makedirs(OUT, exist_ok=True)1314raw = pd.read_csv("data/firms.csv")15n_raw = len(raw)16excl = []17issues = []1819# Step 1: missing values (none expected)20d = raw.dropna()21excl.append(("1", "drop rows with any missing value", n_raw, len(d), n_raw - len(d),22 "complete cases required for OLS"))23# Step 2: duplicate firm_id (none expected)24n0 = len(d); d = d.drop_duplicates("firm_id")25excl.append(("2", "drop duplicate firm_id", n0, len(d), n0 - len(d), "one row per firm"))26# Step 3: negative training hours are physically impossible27n0 = len(d)28n_neg = int((d.training_hours < 0).sum())29d = d[d.training_hours >= 0].copy()30excl.append(("3", "drop training_hours < 0", n0, len(d), n0 - len(d),31 "hours per employee cannot be negative; treated as invalid records"))32if n_neg:33 issues.append({"issue": f"{n_neg} rows have negative training_hours (min {raw.training_hours.min():.3f}), physically impossible",34 "handling": "excluded from main analysis (exclusion step 3); sensitivity re-estimate including them reported in table1 and report"})35pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(36 f"{OUT}/exclusions.csv", index=False)373839def fit(df, cols):40 X = sm.add_constant(df[cols])41 return sm.OLS(df["productivity_growth"], X).fit(cov_type="HC1")4243naive = fit(d, ["training_hours"])44adj = fit(d, ["training_hours", "skill_index"])45adj_all = fit(raw, ["training_hours", "skill_index"]) # sensitivity: keep negative-hour rows46naive_all = fit(raw, ["training_hours"])4748# specification check: quadratic + interaction terms (joint robust Wald test)49_d2 = d.assign(th2=d.training_hours ** 2, sk2=d.skill_index ** 2, thsk=d.training_hours * d.skill_index)50flex = fit(_d2, ["training_hours", "skill_index", "th2", "sk2", "thsk"])51flex_p = float(flex.wald_test("th2 = 0, sk2 = 0, thsk = 0", scalar=True).pvalue)5253rmse = lambda m: float(np.sqrt(np.mean(m.resid ** 2)))54ci = lambda m: m.conf_int(alpha=0.05).loc["training_hours"].astype(float).tolist()55nl, nh = ci(naive); el, eh = ci(adj)56if el > 0: concl = "positive_effect"57elif eh < 0: concl = "negative_effect"58else: concl = "no_detectable_effect"5960res = {61 "task_id": "F6-conf-04",62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": issues,65 "n_input_rows": {"firms.csv": int(n_raw)},66 "n_analysed": int(adj.nobs),67 "random_seed": SEED,68 "naive_estimate": float(naive.params["training_hours"]),69 "naive_ci_low": nl, "naive_ci_high": nh,70 "baseline_metric": rmse(naive),71 "model_metric": rmse(adj),72 "effect_estimate": float(adj.params["training_hours"]),73 "effect_se": float(adj.bse["training_hours"]),74 "effect_ci_low": el, "effect_ci_high": eh,75 "skill_index_coefficient": float(adj.params["skill_index"]),76 "conclusion": concl,77 "se_type": "HC1 heteroskedasticity-robust",78 "naive_se": float(naive.bse["training_hours"]),79 "skill_index_se": float(adj.bse["skill_index"]),80 "sensitivity_incl_negative_effect_estimate": float(adj_all.params["training_hours"]),81 "sensitivity_incl_negative_effect_se": float(adj_all.bse["training_hours"]),82 "sensitivity_incl_negative_naive_estimate": float(naive_all.params["training_hours"]),83 "spec_check_nonlinear_terms_joint_pvalue": flex_p,84 "corr_skill_training": float(d[["skill_index", "training_hours"]].corr().iloc[0, 1]),85}86with open(f"{OUT}/results.json", "w") as f:87 json.dump(res, f, indent=2)8889rows = []90for name, m in [("unadjusted_ols (main sample)", naive), ("adjusted_ols_skill (main sample, target)", adj),91 ("unadjusted_ols (incl. negative hours)", naive_all), ("adjusted_ols_skill (incl. negative hours)", adj_all)]:92 lo, hi = ci(m)93 rows.append({"model": name, "estimate": m.params["training_hours"], "se": m.bse["training_hours"],94 "ci_low": lo, "ci_high": hi, "n": int(m.nobs)})95pd.DataFrame(rows).to_csv(f"{OUT}/table1.csv", index=False)9697fig, ax = plt.subplots(figsize=(7, 5))98ax.scatter(d.training_hours, d.productivity_growth, s=6, alpha=0.4, color="grey", label="firms")99xs = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)100ax.plot(xs, naive.params["const"] + naive.params["training_hours"] * xs, "b-",101 label=f"unadjusted OLS (slope {naive.params['training_hours']:.3f})")102sk = d.skill_index.mean()103ax.plot(xs, adj.params["const"] + adj.params["skill_index"] * sk + adj.params["training_hours"] * xs, "r-",104 label=f"adjusted for skill_index at mean (slope {adj.params['training_hours']:.3f})")105ax.set_xlabel("training hours per employee"); ax.set_ylabel("productivity growth (%)")106ax.set_title("Productivity growth vs training hours (simulated)"); ax.legend()107fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)108109print(json.dumps({k: res[k] for k in ["naive_estimate", "effect_estimate", "effect_se", "effect_ci_low",110 "effect_ci_high", "conclusion", "n_analysed"]}, indent=1))111