1#!/usr/bin/env python32"""WP44 — structural-change toolkit for monthly publication series (RQ2, hypotheses H1/H3).34Functions (all pure, tested in tests/test_structural_change.py):5 its_fit(y, dates, event, lags) interrupted time series on log counts with month-of-year seasonality,6 level and slope change at `event`, Newey–West HAC standard errors7 placebo_breaks(y, dates, window) distribution of the ITS level/slope t-statistics over every candidate month of the8 baseline (placebo), returns the empirical rank/p-value of the true event9 bai_perron(y, dates, max_breaks) multiple structural breaks in mean + trend by dynamic programming (BIC), with10 minimum segment length h (Bai & Perron 1998/2003, pure-python DP on SSR)11 pelt(y, pen) PELT changepoints (ruptures, l2 cost) on log counts, penalty by BIC12 bocpd(y, hazard, mu0, kappa0, ...) Bayesian online changepoint detection (Adams & MacKay 2007), Student-t predictive,13 returns run-length posterior and P(changepoint at t)14Right-truncation: callers must stop the series at the last complete month (see bibliometrics/01_volume_series.py).15"""16from __future__ import annotations17import numpy as np, pandas as pd18import statsmodels.api as sm192021def _design(dates: pd.DatetimeIndex, event: pd.Timestamp, seasonal=True):22 t = np.arange(len(dates), dtype=float)23 post = (dates >= event).astype(float)24 t_post = np.where(post == 1, t - t[post == 1][0] if post.sum() else 0, 0.0)25 X = {"const": np.ones(len(t)), "t": t, "post": post, "t_post": t_post}26 if seasonal:27 for m in range(2, 13):28 X[f"m{m}"] = (dates.month == m).astype(float)29 return pd.DataFrame(X, index=dates)303132def its_fit(y, dates, event, lags=6, seasonal=True, log=True):33 dates = pd.DatetimeIndex(dates); event = pd.Timestamp(event)34 yy = np.log(np.asarray(y, dtype=float)) if log else np.asarray(y, dtype=float)35 X = _design(dates, event, seasonal)36 res = sm.OLS(yy, X).fit(cov_type="HAC", cov_kwds={"maxlags": lags})37 out = {"level_change": res.params["post"], "level_se": res.bse["post"], "level_t": res.tvalues["post"], "level_p": res.pvalues["post"],38 "slope_change": res.params["t_post"], "slope_se": res.bse["t_post"], "slope_t": res.tvalues["t_post"], "slope_p": res.pvalues["t_post"],39 "pre_slope": res.params["t"], "pre_slope_se": res.bse["t"], "n": int(res.nobs), "r2": res.rsquared, "event": str(event.date())}40 return out, res414243def placebo_breaks(y, dates, true_event, min_pre=24, min_post=12, lags=6, seasonal=True):44 dates = pd.DatetimeIndex(dates)45 rows = []46 for i in range(min_pre, len(dates) - min_post):47 ev = dates[i]48 o, _ = its_fit(y, dates, ev, lags, seasonal)49 rows.append({"event": ev, "level_t": o["level_t"], "slope_t": o["slope_t"], "level_change": o["level_change"], "slope_change": o["slope_change"]})50 df = pd.DataFrame(rows)51 tv = df.loc[df.event == pd.Timestamp(true_event)]52 out = {"n_placebo": len(df)}53 if len(tv):54 for k in ("level_t", "slope_t"):55 v = float(tv[k].iloc[0]); out[f"{k}_true"] = v56 out[f"{k}_placebo_p"] = float((np.abs(df[k]) >= abs(v)).mean())57 out[f"{k}_rank"] = int((np.abs(df[k]) >= abs(v)).sum())58 return out, df596061def bai_perron(y, dates=None, max_breaks=3, h=None, trend=True):62 """Global minimiser of SSR with m breaks (Bai & Perron 2003 dynamic programming) for a segment-wise linear trend63 (or mean) model; model selection by BIC and by the LWZ criterion. Returns dict with breaks (indices and dates)."""64 y = np.asarray(y, dtype=float); n = len(y); h = h or max(8, int(0.15 * n))65 t = np.arange(n, dtype=float)66 def ssr(i, j): # segment [i, j)67 yy = y[i:j]68 if trend:69 X = np.column_stack([np.ones(j - i), t[i:j]])70 beta, *_ = np.linalg.lstsq(X, yy, rcond=None); r = yy - X @ beta71 else:72 r = yy - yy.mean()73 return float(r @ r)74 S = np.full((n + 1, n + 1), np.inf)75 for i in range(n):76 for j in range(i + h, n + 1):77 S[i, j] = ssr(i, j)78 best = {}79 k_params = 2 if trend else 180 # DP: C[m][j] = min SSR of first j obs with m breaks81 C = {0: S[0, :].copy()}; B = {0: {j: [] for j in range(n + 1)}}82 for m in range(1, max_breaks + 1):83 C[m] = np.full(n + 1, np.inf); B[m] = {}84 for j in range((m + 1) * h, n + 1):85 cands = [(C[m - 1][i] + S[i, j], i) for i in range(m * h, j - h + 1) if np.isfinite(C[m - 1][i]) and np.isfinite(S[i, j])]86 if cands:87 v, i = min(cands); C[m][j] = v; B[m][j] = B[m - 1][i] + [i]88 for m in range(0, max_breaks + 1):89 tot = C[m][n]90 if not np.isfinite(tot):91 continue92 p = (m + 1) * k_params + m93 bic = n * np.log(tot / n) + p * np.log(n)94 lwz = n * np.log(tot / (n - p)) + p * 0.299 * np.log(n) ** 2.195 best[m] = {"ssr": tot, "bic": bic, "lwz": lwz, "breaks": B[m].get(n, [])}96 m_bic = min(best, key=lambda m: best[m]["bic"]); m_lwz = min(best, key=lambda m: best[m]["lwz"])97 res = {"by_m": best, "m_bic": m_bic, "m_lwz": m_lwz, "breaks_bic": best[m_bic]["breaks"], "h": h}98 if dates is not None:99 d = pd.DatetimeIndex(dates); res["break_dates_bic"] = [str(d[i].date()) for i in best[m_bic]["breaks"]]100 res["break_dates_lwz"] = [str(d[i].date()) for i in best[m_lwz]["breaks"]]101 return res102103104def pelt(y, pen=None, min_size=6, model="l2"):105 import ruptures as rpt106 y = np.asarray(y, dtype=float); n = len(y)107 pen = pen if pen is not None else 2 * np.log(n) * np.var(y) # BIC-type penalty for l2 cost108 algo = rpt.Pelt(model=model, min_size=min_size, jump=1).fit(y.reshape(-1, 1))109 return algo.predict(pen=pen)[:-1]110111112def bocpd(y, hazard=1 / 60.0, mu0=None, kappa0=1.0, alpha0=1.0, beta0=None):113 """Adams & MacKay (2007) with Normal–Inverse-Gamma conjugacy (Student-t predictive). Returns (R, cp_prob, map_rl):114 R[t, r] is the run-length posterior after observing y[:t]; with a constant hazard the mass on r = 0 equals the hazard115 by construction, so the detection statistics are cp_prob[t] = P(r_t <= 3 | y[:t+1]) (a changepoint in the last116 months) and map_rl[t] = argmax_r R[t+1, r] (the MAP run length, which resets after a change)."""117 from scipy.stats import t as student118 y = np.asarray(y, dtype=float); n = len(y)119 mu0 = y[:6].mean() if mu0 is None else mu0; beta0 = y[:6].var() + 1e-6 if beta0 is None else beta0120 R = np.zeros((n + 1, n + 1)); R[0, 0] = 1.0121 mu = np.array([mu0]); kappa = np.array([kappa0]); alpha = np.array([alpha0]); beta = np.array([beta0])122 cp = np.zeros(n); map_rl = np.zeros(n, dtype=int)123 for t_ in range(n):124 x = y[t_]125 df = 2 * alpha; scale = np.sqrt(beta * (kappa + 1) / (alpha * kappa))126 pred = student.pdf(x, df, loc=mu, scale=scale)127 growth = R[t_, : t_ + 1] * pred * (1 - hazard)128 change = (R[t_, : t_ + 1] * pred * hazard).sum()129 R[t_ + 1, 1 : t_ + 2] = growth; R[t_ + 1, 0] = change130 R[t_ + 1] /= R[t_ + 1].sum() + 1e-300131 cp[t_] = R[t_ + 1, :4].sum(); map_rl[t_] = int(np.argmax(R[t_ + 1, : t_ + 2]))132 mu_new = (kappa * mu + x) / (kappa + 1); kappa_new = kappa + 1; alpha_new = alpha + 0.5133 beta_new = beta + kappa * (x - mu) ** 2 / (2 * (kappa + 1))134 mu = np.concatenate([[mu0], mu_new]); kappa = np.concatenate([[kappa0], kappa_new]); alpha = np.concatenate([[alpha0], alpha_new]); beta = np.concatenate([[beta0], beta_new])135 return R, cp, map_rl136