#!/usr/bin/env python3
"""The rebuilt diagnostic, and the null tests that forced the rebuild.

v1 (levels regression of log cost on log revenue) failed P49 twice:
  - the naive "elasticity > 1" rule fires ~50% of the time on random data;
  - the t-stat version fires 98.8% of the time on a pure accounting
    reclassification -- a one-off step in cost with NO change in cost structure.

v2 works in FIRST DIFFERENCES. A one-off level shift contaminates a single
difference instead of tilting the whole slope, and the estimator is the median
quarter-over-quarter ratio, so that one contaminated difference cannot carry it.
"""
import json, os, math, random, statistics as st

OUT = os.path.dirname(os.path.abspath(__file__))
d = json.load(open(os.path.join(OUT, "panel_final.json")))
WIN, BAL, P = d["window"], d["balanced"], d["panel"]


def theil_sen(xs, ys):
    """Median of pairwise slopes -- resistant to a contaminated observation."""
    sl = [(ys[j]-ys[i])/(xs[j]-xs[i])
          for i in range(len(xs)) for j in range(i+1, len(xs))
          if abs(xs[j]-xs[i]) > 1e-9]
    return st.median(sl) if sl else float("nan")


def diff_elasticity(rev, cost):
    """v2: Theil-Sen slope of Dlog(cost) on Dlog(revenue)."""
    dr = [math.log(rev[i+1]/rev[i]) for i in range(len(rev)-1)]
    dc = [math.log(cost[i+1]/cost[i]) for i in range(len(cost)-1)]
    return theil_sen(dr, dc)


def boot_ci(rev, cost, reps=2000, seed=7):
    rnd = random.Random(seed)
    dr = [math.log(rev[i+1]/rev[i]) for i in range(len(rev)-1)]
    dc = [math.log(cost[i+1]/cost[i]) for i in range(len(cost)-1)]
    n = len(dr); out = []
    for _ in range(reps):
        idx = [rnd.randrange(n) for _ in range(n)]
        x = [dr[i] for i in idx]; y = [dc[i] for i in idx]
        if len(set(x)) < 3:
            continue
        s = theil_sen(x, y)
        if s == s:
            out.append(s)
    out.sort()
    return out[int(.025*len(out))], out[int(.975*len(out))]


# ------------------------------------------------------------- null tests
print("=== P49 NULL TESTS ON THE REBUILT DIAGNOSTIC (v2, first differences) ===")
random.seed(20260810)
REPS = 6000


def synth(step=False, noise=0.02, n=17, margin=0.75, g=1.04):
    rev, cost = [], []
    step_at = random.randint(4, n-4)
    for i in range(n):
        r = 100*(g**i)*math.exp(random.gauss(0, noise))
        c = r*(1-margin)*(1.08 if (step and i >= step_at) else 1.0)*math.exp(random.gauss(0, noise))
        rev.append(r); cost.append(c)
    return rev, cost


for label, kw in [("pure noise, constant margin", {}),
                  ("+ one-off cost reclassification", {"step": True})]:
    v1_fire, v2_fire, ests = 0, 0, []
    for _ in range(REPS):
        rev, cost = synth(**kw)
        # v1 for comparison
        xs = [math.log(r) for r in rev]; ys = [math.log(c) for c in cost]
        n = len(xs); mx, my = st.mean(xs), st.mean(ys)
        sxx = sum((x-mx)**2 for x in xs)
        b1 = sum((x-mx)*(y-my) for x, y in zip(xs, ys))/sxx
        a1 = my-b1*mx
        s2 = sum((y-(a1+b1*x))**2 for x, y in zip(xs, ys))/(n-2)
        if (b1-1)/math.sqrt(s2/sxx) > 2:
            v1_fire += 1
        e = diff_elasticity(rev, cost); ests.append(e)
        lo, hi = boot_ci(rev, cost, reps=400, seed=random.randrange(10**6))
        if lo > 1:
            v2_fire += 1
    print(f"\n{label}:")
    print(f"   v1 levels t-stat rule      fires {v1_fire/REPS:6.1%}")
    print(f"   v2 diff + bootstrap CI>1   fires {v2_fire/REPS:6.1%}   "
          f"(median estimate {st.median(ests):.3f})")

# ---- power: does v2 still detect a REAL structural change? --------------
print("\n=== POWER — can v2 still see a real effect? ===")
for extra in (0.10, 0.20, 0.30):
    hit = 0
    for _ in range(REPS):
        rev, cost = [], []
        for i in range(17):
            r = 100*(1.04**i)*math.exp(random.gauss(0, 0.02))
            share = extra*(i/16)                       # variable-cost share ramps in
            c = r*(0.25*(1-share) + 0.55*share)*math.exp(random.gauss(0, 0.02))
            rev.append(r); cost.append(c)
        lo, hi = boot_ci(rev, cost, reps=400, seed=random.randrange(10**6))
        if lo > 1:
            hit += 1
    print(f"   variable-cost mix reaching {extra:.0%} of revenue by the end: "
          f"detected {hit/REPS:5.1%} of the time")

# ------------------------------------------------------- the real panel
print("\n=== THE REBUILT DIAGNOSTIC ON THE 53-FIRM PANEL ===")
res = {}
for t in BAL:
    rev = [P[t]["quarters"][q]["revenue"] for q in WIN]
    cost = [P[t]["quarters"][q]["cost_of_revenue"] for q in WIN]
    if min(cost) <= 0:
        continue
    e = diff_elasticity(rev, cost)
    lo, hi = boot_ci(rev, cost)
    res[t] = (e, lo, hi)

flagged = sorted([t for t in res if res[t][1] > 1], key=lambda t: -res[t][0])
print(f"{'tick':6}{'elast':>8}{'95% CI':>18}  name")
for t in sorted(res, key=lambda t: -res[t][0])[:10]:
    e, lo, hi = res[t]
    mark = "  <-- CI excludes 1" if lo > 1 else ""
    print(f"{t:6}{e:8.3f}   [{lo:6.3f},{hi:6.3f}]  {P[t]['name'][:30]}{mark}")
vals = [v[0] for v in res.values()]
print(f"\nn={len(vals)}  median {st.median(vals):.3f}")
print(f"  cost growing FASTER than revenue (CI entirely above 1): "
      f"{len(flagged)}/{len(vals)}  {flagged}")
print(f"  cost growing SLOWER (CI entirely below 1): "
      f"{sum(1 for v in res.values() if v[2] < 1)}/{len(vals)}")

json.dump({t: {"elasticity": v[0], "ci_lo": v[1], "ci_hi": v[2]} for t, v in res.items()},
          open(os.path.join(OUT, "elasticity_v2.json"), "w"), indent=1)