1"""Replication script for F1-dev-01. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import pandas as pd67os.makedirs("out", exist_ok=True)8raw = pd.read_csv("data/firms.csv", dtype=str, keep_default_na=False)9n_rows = len(raw)10for c in raw.columns:11 raw[c] = raw[c].str.strip()1213# Baseline: raw revenue mean over non-empty cells, mixed units, before cleaning14rev_raw = pd.to_numeric(raw["revenue"].replace("", None), errors="coerce")15mean_revenue_raw = float(rev_raw.dropna().mean())1617ledger, issues = [], []18def step(rule, before, after, rationale):19 ledger.append(dict(step=len(ledger) + 1, rule=rule, n_before=before,20 n_after=after, n_removed=before - after, rationale=rationale))2122ledger_start = n_rows23step("read data/firms.csv", n_rows, n_rows, "All data rows read (header excluded).")2425# 1. exact duplicate rows26n_exact = int(raw.duplicated().sum())27d = raw.drop_duplicates().copy()28step("drop exact duplicate rows", n_rows, len(d),29 "Rows identical in every column are assembly duplicates of the same firm; keep first occurrence.")30issues.append({"issue": f"{n_exact} rows are exact duplicates of another row (same firm_id and all values)",31 "handling": "Dropped duplicates, keeping one row per firm."})32# check remaining id conflicts33n_conf = int(d["firm_id"].duplicated().sum())34assert n_conf == 0, "conflicting duplicate firm_ids present"35n_unique = int(d["firm_id"].nunique())3637# 2. unit recoding38units = d["unit"].value_counts().to_dict()39assert set(units) <= {"USD_millions", "USD_thousands"}, units40d["revenue_num"] = pd.to_numeric(d["revenue"].replace("", None), errors="coerce")41d["revenue_musd"] = d["revenue_num"].where(d["unit"] == "USD_millions", d["revenue_num"] / 1000.0)42step("convert revenue to millions USD", len(d), len(d),43 "Rows with unit USD_thousands divided by 1000; USD_millions unchanged. No rows removed.")44issues.append({"issue": f"Mixed revenue units: {units.get('USD_thousands',0)} distinct-firm rows in USD_thousands, "45 f"{units.get('USD_millions',0)} in USD_millions",46 "handling": "Converted USD_thousands to millions USD (divide by 1000)."})4748# employees: complete over all firms49d["employees"] = pd.to_numeric(d["employees"], errors="raise").astype(int)50mean_employees = float(d["employees"].mean())5152# 3. missing revenue53miss = d["revenue_musd"].isna()54n_missing = int(miss.sum())55missing_ids = sorted(d.loc[miss, "firm_id"])56before = len(d)57a = d.loc[~miss].copy()58step("exclude firms with unknown revenue", before, len(a),59 "Empty revenue cell; cannot contribute to revenue aggregates. Firms: " + ";".join(missing_ids))60issues.append({"issue": f"{n_missing} distinct firms have empty revenue ({', '.join(missing_ids)})",61 "handling": "Excluded from revenue aggregates and listed in the exclusion ledger; retained for mean_employees."})6263assert (a["revenue_musd"] > 0).all()64a = a.sort_values("firm_id")65a[["firm_id", "sector", "revenue_musd", "employees"]].to_csv("out/analytic.csv", index=False)66pd.DataFrame(ledger).to_csv("out/exclusions.csv", index=False)6768res = {69 "task_id": "F1-dev-01",70 "n_unique_firms": n_unique,71 "n_missing_revenue": n_missing,72 "n_analysed": int(len(a)),73 "mean_revenue_musd": round(float(a["revenue_musd"].mean()), 6),74 "median_revenue_musd": round(float(a["revenue_musd"].median()), 6),75 "mean_employees": round(mean_employees, 6),76 "mean_revenue_raw": round(mean_revenue_raw, 6),77 "feasibility": "complete",78 "feasibility_reason": None,79 "data_issues": issues,80 "n_input_rows": {"firms.csv": int(n_rows)},81 "random_seed": None,82}83with open("out/results.json", "w") as f:84 json.dump(res, f, indent=2)85print(json.dumps(res, indent=2))8687# Descriptive uncertainty (printed only, quoted in report; not part of the results schema)88import numpy as np89x = a["revenue_musd"].to_numpy()90se = x.std(ddof=1) / np.sqrt(len(x))91print(f"SE(mean revenue) = {se:.3f} MUSD; normal 95% CI = [{x.mean()-1.96*se:.3f}, {x.mean()+1.96*se:.3f}]")92print(f"SD revenue = {x.std(ddof=1):.3f}; min = {x.min():.3f}; max = {x.max():.3f}")93e = d["employees"].to_numpy(); print(f"SE(mean employees) = {e.std(ddof=1)/np.sqrt(len(e)):.3f}")94