1"""Run-level inference. The unit of observation is one simulation run (cell × seed); trades and ticks2within a run are never treated as independent replicates.34* paired_effects : matched-seed differences treatment − control with t-based 95 % CI, Cohen's d (paired),5 and a Wilcoxon signed-rank p-value.6* ols_seed_fe : y = α_seed + Σ_k β_k 1[cell = k] + ε with cluster-robust (by seed) standard errors.7* dose_response : effect at each treatment level vs control (paired), with Holm-adjusted p-values.8* local_projection: β_h from (gap_{t+h} − gap_{t-1}) regressed on treatment, h = 0..H, around a shock.9"""10from __future__ import annotations1112import numpy as np13import pandas as pd14from scipy import stats151617def holm(p: np.ndarray) -> np.ndarray:18 p = np.asarray(p, dtype=float)19 m = np.isfinite(p).sum()20 order = np.argsort(np.where(np.isfinite(p), p, np.inf))21 adj = np.full_like(p, np.nan)22 running = 0.023 for rank, idx in enumerate(order[:m]):24 val = min(1.0, (m - rank) * p[idx])25 running = max(running, val)26 adj[idx] = running27 return adj282930def paired_effects(df: pd.DataFrame, outcome: str, treat_col: str, treat_val, control_val, seed_col: str = "seed") -> dict:31 a = df[df[treat_col] == treat_val].set_index(seed_col)[outcome]32 b = df[df[treat_col] == control_val].set_index(seed_col)[outcome]33 idx = a.index.intersection(b.index)34 d = (a.loc[idx] - b.loc[idx]).astype(float).dropna()35 n = len(d)36 if n < 2:37 return {"outcome": outcome, "level": treat_val, "n": n, "diff": float(d.mean()) if n else np.nan}38 se = d.std(ddof=1) / np.sqrt(n)39 tcrit = stats.t.ppf(0.975, n - 1)40 t = d.mean() / se if se > 0 else np.nan41 p = 2 * stats.t.sf(abs(t), n - 1) if np.isfinite(t) else np.nan42 try:43 pw = stats.wilcoxon(d).pvalue if (d != 0).any() else 1.044 except ValueError:45 pw = np.nan46 return {"outcome": outcome, "level": treat_val, "n": n, "control_mean": float(b.loc[idx].mean()), "treat_mean": float(a.loc[idx].mean()),47 "diff": float(d.mean()), "se": float(se), "ci_lo": float(d.mean() - tcrit * se), "ci_hi": float(d.mean() + tcrit * se),48 "t": float(t), "p": float(p), "p_wilcoxon": float(pw), "cohen_d": float(d.mean() / d.std(ddof=1)) if d.std(ddof=1) > 0 else np.nan,49 "pct_change": float(100 * d.mean() / abs(b.loc[idx].mean())) if b.loc[idx].mean() != 0 else np.nan}505152def dose_response(df: pd.DataFrame, outcome: str, treat_col: str, control_val, levels: list | None = None, seed_col: str = "seed") -> pd.DataFrame:53 levels = levels or sorted(v for v in df[treat_col].unique() if v != control_val)54 rows = [paired_effects(df, outcome, treat_col, lv, control_val, seed_col) for lv in levels]55 out = pd.DataFrame(rows)56 if "p" in out:57 out["p_holm"] = holm(out["p"].to_numpy())58 return out596061def ols_seed_fe(df: pd.DataFrame, outcome: str, cell_col: str, control_val, seed_col: str = "seed") -> pd.DataFrame:62 """Cluster-robust (by seed) OLS with seed fixed effects; returns one row per non-control cell."""63 import statsmodels.formula.api as smf64 d = df[[outcome, cell_col, seed_col]].dropna().copy()65 d["cell"] = d[cell_col].astype(str)66 ctrl = str(control_val)67 cats = [ctrl] + sorted(c for c in d["cell"].unique() if c != ctrl)68 d["cell"] = pd.Categorical(d["cell"], categories=cats)69 d["seed_f"] = d[seed_col].astype(str)70 d = d.rename(columns={outcome: "y"})71 m = smf.ols("y ~ C(cell) + C(seed_f)", data=d).fit(cov_type="cluster", cov_kwds={"groups": d[seed_col]})72 rows = []73 for c in cats[1:]:74 name = f"C(cell)[T.{c}]"75 if name in m.params:76 ci = m.conf_int().loc[name]77 rows.append({"outcome": outcome, "level": c, "beta": m.params[name], "se": m.bse[name], "ci_lo": ci[0], "ci_hi": ci[1],78 "p": m.pvalues[name], "n": int(m.nobs), "r2": m.rsquared})79 out = pd.DataFrame(rows)80 if len(out):81 out["p_holm"] = holm(out["p"].to_numpy())82 return out838485def local_projection(ticks_by_run: dict, treat: dict, shock_t: int, H: int = 40, control_val=0.0) -> pd.DataFrame:86 """ticks_by_run: {run_id: ticks DataFrame}; treat: {run_id: treatment value}. Gap = log(mid) - log(full_value).87 For each horizon h, regress (gap_{t+h} - gap_{t-1}) on 1[treated] across runs (paired by seed not required)."""88 rows = []89 gaps = {}90 for rid, tk in ticks_by_run.items():91 d = tk[tk.asset == 0].set_index("t")92 g = np.log(d["mid"]) - np.log(d["full_value"])93 gaps[rid] = g94 for h in range(0, H + 1):95 y, x = [], []96 for rid, g in gaps.items():97 if shock_t - 1 in g.index and shock_t + h in g.index:98 y.append(g.loc[shock_t + h] - g.loc[shock_t - 1])99 x.append(1.0 if treat[rid] != control_val else 0.0)100 y, x = np.array(y), np.array(x)101 if len(y) < 4 or x.var() == 0:102 continue103 X = np.column_stack([np.ones_like(x), x])104 beta, *_ = np.linalg.lstsq(X, y, rcond=None)105 resid = y - X @ beta106 # HC1 robust SE107 XtX_inv = np.linalg.inv(X.T @ X)108 V = XtX_inv @ (X.T * resid ** 2) @ X @ XtX_inv * len(y) / (len(y) - 2)109 rows.append({"h": h, "beta": beta[1], "se": np.sqrt(V[1, 1]), "control_path": float(y[x == 0].mean()), "treated_path": float(y[x == 1].mean())})110 return pd.DataFrame(rows)111112113def power_seeds(d: float, power: float = 0.8, alpha: float = 0.05) -> int:114 """Paired design: seeds needed to detect a standardized effect d (one-sample t approximation)."""115 z = stats.norm.ppf(1 - alpha / 2) + stats.norm.ppf(power)116 return int(np.ceil((z / d) ** 2))117118119def summarize_cells(df: pd.DataFrame, outcomes: list[str], cell_col: str = "name") -> pd.DataFrame:120 g = df.groupby(cell_col)[outcomes]121 m, s, n = g.mean(), g.std(ddof=1), g.size()122 out = m.copy()123 for c in outcomes:124 out[c + "_se"] = s[c] / np.sqrt(n)125 out["n_runs"] = n126 return out127