#!/usr/bin/env python3
"""What the filings say about software gross margin, 2022Q1-2026Q1."""
import json, os, math, 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"]

print(f"BALANCED PANEL: {len(BAL)} firms x {len(WIN)} quarters "
      f"= {len(BAL)*len(WIN)} firm-quarters\n")

# ---- 1. the aggregate series -------------------------------------------
print("quarter   agg(rev-wtd)  equal-wtd    median      p25     p75   n")
agg_series, ew_series = [], []
for q in WIN:
    ms = [P[t]["quarters"][q]["gross_margin"] for t in BAL]
    R = sum(P[t]["quarters"][q]["revenue"] for t in BAL)
    G = sum(P[t]["quarters"][q]["gross_profit"] for t in BAL)
    agg, ew = G / R, st.mean(ms)
    agg_series.append(agg); ew_series.append(ew)
    qs = st.quantiles(ms, n=4)
    print(f"{q}      {agg*100:6.2f}%     {ew*100:6.2f}%   {st.median(ms)*100:6.2f}%  "
          f"{qs[0]*100:6.2f}% {qs[2]*100:6.2f}%  {len(ms)}")

print(f"\nAggregate  {WIN[0]} -> {WIN[-1]}: {agg_series[0]*100:.2f}% -> "
      f"{agg_series[-1]*100:.2f}%   change {(agg_series[-1]-agg_series[0])*100:+.2f} pp")
print(f"Equal-wtd  {WIN[0]} -> {WIN[-1]}: {ew_series[0]*100:.2f}% -> "
      f"{ew_series[-1]*100:.2f}%   change {(ew_series[-1]-ew_series[0])*100:+.2f} pp")

# ---- 2. firm-level change, 4-quarter blocks to kill seasonality ---------
PRE = WIN[:4]                 # 2022Q1-Q4
POST = WIN[-4:]               # 2025Q2-2026Q1
print(f"\nPRE  = {PRE}\nPOST = {POST}")

deltas = {}
for t in BAL:
    pre = st.mean(P[t]["quarters"][q]["gross_margin"] for q in PRE)
    post = st.mean(P[t]["quarters"][q]["gross_margin"] for q in POST)
    deltas[t] = (pre, post, post - pre)

ds = [v[2] for v in deltas.values()]
up = sum(1 for x in ds if x > 0)
print(f"\nFirm-level change (pp): mean {st.mean(ds)*100:+.2f}  median "
      f"{st.median(ds)*100:+.2f}  sd {st.stdev(ds)*100:.2f}")
print(f"  improved: {up}/{len(ds)}   declined: {len(ds)-up}/{len(ds)}")

# paired t-test
n = len(ds); mean_d = st.mean(ds); se = st.stdev(ds)/math.sqrt(n)
tstat = mean_d/se
def t_p(t, df):                       # two-sided p via incomplete beta
    x = df/(df+t*t)
    def betacf(a,b,x):
        MAXIT,EPS,FPMIN=200,3e-16,1e-300
        qab,qap,qam=a+b,a+1,a-1
        c=1.0; dd=1-qab*x/qap
        if abs(dd)<FPMIN: dd=FPMIN
        dd=1/dd; h=dd
        for m in range(1,MAXIT+1):
            m2=2*m
            aa=m*(b-m)*x/((qam+m2)*(a+m2))
            dd=1+aa*dd; c=1+aa/c
            if abs(dd)<FPMIN: dd=FPMIN
            if abs(c)<FPMIN: c=FPMIN
            dd=1/dd; h*=dd*c
            aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
            dd=1+aa*dd; c=1+aa/c
            if abs(dd)<FPMIN: dd=FPMIN
            if abs(c)<FPMIN: c=FPMIN
            dd=1/dd; de=dd*c; h*=de
            if abs(de-1)<EPS: break
        return h
    a,b=df/2,0.5
    lbeta=math.lgamma(a)+math.lgamma(b)-math.lgamma(a+b)
    if x < (a+1)/(a+b+2):
        ib=math.exp(a*math.log(x)+b*math.log(1-x)-lbeta)*betacf(a,b,x)/a
    else:
        ib=1-math.exp(b*math.log(1-x)+a*math.log(x)-lbeta)*betacf(b,a,1-x)/b
    return ib
print(f"  paired t = {tstat:+.3f} (df {n-1}), p = {t_p(tstat, n-1):.4f}")

# Wilcoxon signed-rank, normal approximation
nz = [x for x in ds if x != 0]
ranks = sorted(range(len(nz)), key=lambda i: abs(nz[i]))
rk = [0]*len(nz)
i = 0
while i < len(ranks):
    j = i
    while j+1 < len(ranks) and abs(nz[ranks[j+1]]) == abs(nz[ranks[i]]): j += 1
    avg = (i+j)/2 + 1
    for k in range(i, j+1): rk[ranks[k]] = avg
    i = j+1
Wp = sum(rk[i] for i in range(len(nz)) if nz[i] > 0)
N = len(nz); mu = N*(N+1)/4; sd = math.sqrt(N*(N+1)*(2*N+1)/24)
z = (Wp-mu)/sd
print(f"  Wilcoxon W+ = {Wp:.0f}, z = {z:+.3f}, p = "
      f"{2*(1-0.5*(1+math.erf(abs(z)/math.sqrt(2)))):.4f}")

# ---- 3. biggest movers -------------------------------------------------
print("\nLargest declines:")
for t, v in sorted(deltas.items(), key=lambda kv: kv[1][2])[:8]:
    print(f"  {t:6}{P[t]['name'][:30]:32}{v[0]*100:6.1f}% -> {v[1]*100:6.1f}%  {v[2]*100:+6.1f} pp")
print("Largest gains:")
for t, v in sorted(deltas.items(), key=lambda kv: -kv[1][2])[:8]:
    print(f"  {t:6}{P[t]['name'][:30]:32}{v[0]*100:6.1f}% -> {v[1]*100:6.1f}%  {v[2]*100:+6.1f} pp")

# ---- 4. dispersion over time -------------------------------------------
print("\nDispersion (sd of firm margins, pp):")
for q in [WIN[0], WIN[4], WIN[8], WIN[12], WIN[-1]]:
    ms = [P[t]["quarters"][q]["gross_margin"] for t in BAL]
    print(f"  {q}  sd={st.stdev(ms)*100:5.2f}  IQR="
          f"{(st.quantiles(ms,n=4)[2]-st.quantiles(ms,n=4)[0])*100:5.2f}")

# ---- 5. what would it take to move the blended line? -------------------
print("\nSENSITIVITY — blended margin change from an AI mix shift")
print("base margin 80%; AI revenue at margin m and share s of total revenue")
print(f"{'share':>7}", "".join(f"{m:>9.0%}" for m in (0.7,0.6,0.5,0.4,0.3)))
for s in (0.01, 0.03, 0.05, 0.10, 0.20):
    row = "".join(f"{((0.80*(1-s)+m*s)-0.80)*100:>8.2f}p" for m in (0.7,0.6,0.5,0.4,0.3))
    print(f"{s:>7.0%}", row)

json.dump({"agg": agg_series, "ew": ew_series, "window": WIN,
           "deltas": {k: list(v) for k, v in deltas.items()}},
          open(os.path.join(OUT, "results.json"), "w"), indent=1)