#!/usr/bin/env python3
"""The diagnostic the essay will prescribe, and P49's null test of it.

Margin LEVEL answers "how profitable", not "is cost becoming variable". The second
question is an elasticity: regress log(cost of revenue) on log(revenue) across
quarters. Under a fixed margin the elasticity is exactly 1 by construction. Above 1,
cost is growing superlinearly in revenue -- which is what "COGS became variable"
means mechanically.

P49 requires this be run under the null before it is published: on synthetic data
with NO effect present, at the reader's realistic n, against the artefacts the
design invites.
"""
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 ols(xs, ys):
    """slope, intercept, se(slope)"""
    n = len(xs); mx, my = st.mean(xs), st.mean(ys)
    sxx = sum((x-mx)**2 for x in xs)
    sxy = sum((x-mx)*(y-my) for x, y in zip(xs, ys))
    b = sxy/sxx; a = my - b*mx
    resid = [y-(a+b*x) for x, y in zip(xs, ys)]
    s2 = sum(r*r for r in resid)/(n-2)
    return b, a, math.sqrt(s2/sxx)


print("=== ELASTICITY OF COST OF REVENUE TO REVENUE, 2022Q1-2026Q1 ===")
print("(1.00 = fixed margin by construction; >1 = cost superlinear in revenue)\n")
els = {}
for t in BAL:
    xs = [math.log(P[t]["quarters"][q]["revenue"]) for q in WIN]
    ys = [math.log(P[t]["quarters"][q]["cost_of_revenue"]) for q in WIN]
    if min(ys) <= 0 or any(P[t]["quarters"][q]["cost_of_revenue"] <= 0 for q in WIN):
        continue
    b, a, se = ols(xs, ys)
    els[t] = (b, se, (b-1)/se)

print(f"{'tick':6}{'elast':>8}{'se':>7}{'t vs 1':>8}  name")
for t, (b, se, tv) in sorted(els.items(), key=lambda kv: -kv[1][0])[:12]:
    print(f"{t:6}{b:8.3f}{se:7.3f}{tv:8.2f}  {P[t]['name'][:34]}")
print("  ...")
for t, (b, se, tv) in sorted(els.items(), key=lambda kv: -kv[1][0])[-6:]:
    print(f"{t:6}{b:8.3f}{se:7.3f}{tv:8.2f}  {P[t]['name'][:34]}")

vals = [v[0] for v in els.values()]
print(f"\nn={len(vals)}  median elasticity {st.median(vals):.3f}  mean {st.mean(vals):.3f}")
print(f"  above 1.0: {sum(1 for v in vals if v>1)}/{len(vals)}")
print(f"  significantly above 1 (t>2): {sum(1 for v in els.values() if v[2]>2)}/{len(vals)}")
print(f"  significantly below 1 (t<-2): {sum(1 for v in els.values() if v[2]<-2)}/{len(vals)}")

# ---------------------------------------------------------------- P49 null
print("\n=== P49 NULL TEST ===")
print("Synthetic firm, TRUE margin constant, only reporting noise. n=17 quarters.")
print("If the diagnostic fires here, it is not a diagnostic.\n")
random.seed(20260810)
REPS = 20000

for noise in (0.01, 0.02, 0.04):
    fired_t2, fired_raw, elastics = 0, 0, []
    for _ in range(REPS):
        rev0, g = 100.0, 1.04           # 4%/qtr growth, ~17% p.a.
        margin = 0.75
        xs, ys = [], []
        for i in range(17):
            r = rev0*(g**i)*math.exp(random.gauss(0, noise))
            c = r*(1-margin)*math.exp(random.gauss(0, noise))
            xs.append(math.log(r)); ys.append(math.log(c))
        b, a, se = ols(xs, ys)
        elastics.append(b)
        if (b-1)/se > 2: fired_t2 += 1
        if b > 1: fired_raw += 1
    print(f"noise sd={noise:.0%}  elasticity: median {st.median(elastics):.3f} "
          f"sd {st.stdev(elastics):.3f}")
    print(f"   naive rule 'elasticity > 1'        fires {fired_raw/REPS:6.1%} of the time  <-- coin flip")
    print(f"   rule 't-stat vs 1 exceeds 2'       fires {fired_t2/REPS:6.1%} of the time")

# the artefact the design invites: revenue growth alone can induce spurious elasticity
# if cost is measured with error correlated to scale -- test a second null
print("\nSecond null: TRUE margin constant but cost carries a one-off step change")
print("(an accounting reclassification, not an economics change)")
fired = 0
for _ in range(REPS):
    xs, ys = [], []
    step_at = random.randint(4, 13)
    for i in range(17):
        r = 100*(1.04**i)*math.exp(random.gauss(0, 0.02))
        c = r*0.25*(1.08 if i >= step_at else 1.0)*math.exp(random.gauss(0, 0.02))
        xs.append(math.log(r)); ys.append(math.log(c))
    b, a, se = ols(xs, ys)
    if (b-1)/se > 2: fired += 1
print(f"   rule 't-stat vs 1 exceeds 2' fires {fired/REPS:.1%} of the time on a pure "
      f"reclassification\n   -> the test CANNOT separate a cost-structure change from a "
      f"one-off restatement.")

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