1#!/usr/bin/env python32"""WP44 — RQ2/H1/H3 driver: interrupted time series, placebo breaks, Bai–Perron, PELT and BOCPD on monthly series.34Series: total eligible works (month-precision records only), by area and by category (results/graph/category_counts_monthly.csv,5n_month_precision column when present; before classification only the total is run). Window: 2018-01 → --last-complete6(default 2026-06; sensitivity 2026-08 and 2026-03). Event dates: 2022-11-30 (ChatGPT, primary) and alternatives72020-05 (GPT-3), 2022-04 (DALL·E 2), 2023-03 (GPT-4), 2024-09 (o1). Outputs results/econometrics/its_results.json,8placebo_<series>.csv, breaks.csv, and a per-series CSV of fitted values."""9import argparse, json, os, sys, numpy as np, pandas as pd10sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))11from econometrics.structural_change import its_fit, placebo_breaks, bai_perron, pelt, bocpd1213EVENTS = {"chatgpt": "2022-12-01", "gpt3": "2020-06-01", "dalle2": "2022-05-01", "gpt4": "2023-04-01", "o1": "2024-10-01"}141516def series_total(R):17 d = pd.read_csv(f"{R}/results/bibliometrics/volume_monthly.csv")18 if "date_precision" in d.columns:19 d = d[d.date_precision == "month"]20 m = d.groupby("first_month")["n"].sum(); return {"total": m}212223def series_categories(R):24 p = f"{R}/results/graph/category_counts_monthly.csv"25 if not os.path.exists(p):26 return {}27 d = pd.read_csv(p); col = "n_month_precision" if "n_month_precision" in d.columns else "n_multilabel"28 return {f"cat_{c}": g.groupby("first_month")[col].sum() for c, g in d.groupby("category_id")}293031def run_one(name, s, first, last, out, placebo_baseline_only=True, donut=0):32 s = s[(s.index >= first) & (s.index <= last)].astype(float)33 idx = pd.period_range(first, last, freq="M").strftime("%Y-%m"); s = s.reindex(idx).fillna(0.0)34 if (s > 0).sum() < 48 or s.sum() < 500:35 return None36 y = s.values + 1.0; dates = pd.to_datetime(idx + "-01"); res = {"series": name, "n_months": len(y), "total": float(s.sum())}37 for ev, d in EVENTS.items():38 if pd.Timestamp(d) <= dates[12] or pd.Timestamp(d) >= dates[-12]:39 continue40 o, fit = its_fit(y, dates, d); res[f"its_{ev}"] = o41 if ev == "chatgpt":42 pl, df = placebo_breaks(y, dates, d); res["placebo_chatgpt"] = pl; df.to_csv(f"{out}/placebo_{name}.csv", index=False)43 if placebo_baseline_only:44 dfb = df[df.event < pd.Timestamp(d)]; tv = df[df.event == pd.Timestamp(d)]45 if len(tv) and len(dfb):46 res["placebo_chatgpt_baseline_only"] = {"n_placebo": int(len(dfb)), **{f"{k}_placebo_p": float((np.abs(dfb[k]) >= abs(float(tv[k].iloc[0]))).mean()) for k in ("level_t", "slope_t")}}47 if donut:48 keep = np.abs((dates - pd.Timestamp(d)).days) > donut * 3149 od, _ = its_fit(y[keep], dates[keep], d); res["its_chatgpt_donut"] = {"months_excluded_each_side": donut, **od}50 pd.DataFrame({"month": idx, "y": s.values, "fitted": np.exp(fit.fittedvalues) - 1}).to_csv(f"{out}/fit_{name}.csv", index=False)51 ly = np.log(y); months_of_year = dates.month; base = dates < pd.Timestamp("2022-12-01")52 seas = {m: ly[base & (months_of_year == m)].mean() - ly[base].mean() for m in range(1, 13)}53 ly_sa = ly - np.array([seas[m] for m in months_of_year]) # seasonal adjustment from baseline month effects54 bp = bai_perron(ly_sa, dates, max_breaks=3); res["bai_perron"] = {k: v for k, v in bp.items() if k != "by_m"}; res["bai_perron"]["input"] = "log, seasonally adjusted (baseline month means)"55 res["bai_perron"]["bic_by_m"] = {m: round(v["bic"], 2) for m, v in bp["by_m"].items()}56 res["pelt"] = [idx[i] for i in pelt(ly_sa) if i < len(idx)]57 _, cp, mrl = bocpd(ly_sa); resets = [idx[t] for t in range(12, len(idx)) if mrl[t] <= 2 and mrl[t - 1] > 12]58 res["bocpd_resets"] = resets; res["bocpd_cp_prob_2022_12_to_2023_03"] = float(cp[[i for i, m in enumerate(idx) if "2022-12" <= m <= "2023-03"]].max()) if any("2022-12" <= m <= "2023-03" for m in idx) else None59 return res606162def main():63 ap = argparse.ArgumentParser(); ap.add_argument("--root", required=True); ap.add_argument("--first", default="2018-02"); ap.add_argument("--last-complete", default="2026-06")64 ap.add_argument("--sensitivity", default="2026-03,2026-08"); ap.add_argument("--donut", type=int, default=3); a = ap.parse_args(); R = a.root65 out = f"{R}/results/econometrics"; os.makedirs(out, exist_ok=True)66 S = series_total(R); S.update(series_categories(R)); all_res = {"last_complete": a.last_complete, "first": a.first, "events": EVENTS, "series": {}}67 for name, s in S.items():68 r = run_one(name, s, a.first, a.last_complete, out, donut=a.donut)69 if r:70 all_res["series"][name] = r; print(name, "ITS chatgpt level", round(r.get("its_chatgpt", {}).get("level_change", float("nan")), 3), "p", round(r.get("its_chatgpt", {}).get("level_p", float("nan")), 3), "placebo p", r.get("placebo_chatgpt", {}).get("level_t_placebo_p"), "BP", r["bai_perron"].get("break_dates_bic"), flush=True)71 sens = {}72 for lc in a.sensitivity.split(","):73 r = run_one("total", S["total"], a.first, lc, out + "/sens_tmp" if os.makedirs(out + "/sens_tmp", exist_ok=True) is None else out)74 if r:75 sens[lc] = {"its_chatgpt": r.get("its_chatgpt"), "bai_perron": r["bai_perron"].get("break_dates_bic")}76 all_res["sensitivity_last_complete_total"] = sens77 json.dump(all_res, open(f"{out}/its_results.json", "w"), indent=1, default=float)78 rows = [{"series": n, **{f"{ev}_level": r.get(f"its_{ev}", {}).get("level_change") for ev in EVENTS}, **{f"{ev}_level_p": r.get(f"its_{ev}", {}).get("level_p") for ev in EVENTS},79 "chatgpt_slope": r.get("its_chatgpt", {}).get("slope_change"), "chatgpt_slope_p": r.get("its_chatgpt", {}).get("slope_p"), "placebo_level_p": r.get("placebo_chatgpt", {}).get("level_t_placebo_p"),80 "placebo_slope_p": r.get("placebo_chatgpt", {}).get("slope_t_placebo_p"), "bp_breaks_bic": ";".join(r["bai_perron"].get("break_dates_bic", [])), "pelt": ";".join(r["pelt"]), "bocpd_resets": ";".join(r["bocpd_resets"])} for n, r in all_res["series"].items()]81 pd.DataFrame(rows).to_csv(f"{out}/breaks.csv", index=False); print("series analysed:", len(rows))828384if __name__ == "__main__":85 main()86