"""Step 2: turn the paper's capital estimate into a revenue requirement and a bill impact.

A transparent rate-base model, not EPRI's. Capital in constant 2021 dollars is spent
2026-2050, enters service the year it is spent, depreciates straight-line, and earns a
pre-tax return on average net plant. Revenue requirement = depreciation + return +
other carrying costs. Deferred taxes (ADIT) are ignored, which overstates slightly.

Verified inputs are marked V (with source); judgment inputs are marked J.
Run:  python -I step2_revenue_requirement.py
"""
import json

YEARS = list(range(2026, 2051))          # 25 years, matches the paper's /25 average (p.29)
MULTS = (1.0, 1.70, 1.95, 2.33)   # sales vs 2024. 1.0 = no growth (reference only). 1.70 / 1.96 / 2.33 = Pipeline Gas / All Options / Limited Efficiency: Roadmap Energy Pathways Fig.11 (p.43) final electricity demand 2050 less rooftop PV 8.3 TWh (Table 22), over 2024 EIA sales 49.419 TWh

# ---- capital cases ($B, from Figure 14 bar labels on p.5) ----------------------
CAP_HIGH = 22.84 + 23.04                 # V: no mitigation, 45.88
CAP_LOW = 27.19                          # V: all mitigations (implied by $1,088M/yr x 25)
HOST_GROSS = 23.04                       # V: solar-hosting bar, p.5
CAP_LOW_RP = CAP_LOW - 0.5 * HOST_GROSS  # low case with developers paying 50% of hosting (p.30 assumption)

# ---- financial parameters -------------------------------------------------------
ROE = 0.0935       # V: National Grid MECo, D.P.U. 23-150 (Oct 2024); Eversource NSTAR Electric 9.80% (verified: D.P.U. 22-22 final order, 30 Nov 2022, p.421 of the PDF)
EQUITY = 0.5283    # V: same MECo order
KD = 0.05          # J: cost of debt 4-6%
TAX = 1 - (1 - 0.21) * (1 - 0.08)        # J/V: 21% federal, 8% MA corporate rate -> 27.3%
OTHER = 0.01       # J: property tax and other carrying cost, share of net plant
LIFE = 40          # J: distribution plant book life, 30-50
INFL = 0.025       # J: used only in the "real" variant

# ---- sales and bill (V unless noted) --------------------------------------------
SALES_TWH = 49.419                       # V: EIA Table 2.8, MA all sectors 2024
HH_KWH_MONTH = 600                       # V: repo's R-1 bill basis (MASTER_DATA.md)
R1_BILL = 225.75                         # V: repo, Eversource R-1 at 600 kWh, Oct 2026
DIST_C = 9.837                           # V: Eversource R-1 Greater Boston, filed rate effective 1 Sep 2026 (was 9.443 until 31 Aug)
ALLIN_C = 30.49                          # V: repo, MA residential all-in cents/kWh (EIA, Jul 2026)


def weights(profile):
    n = len(YEARS)
    if profile == "even":
        w = [1.0] * n
    elif profile == "ramp":                              # back-loaded
        w = [i + 1.0 for i in range(n)]
    elif profile == "mid":                               # peak ~2040 (paper: fastest growth 2030-40)
        w = [max(0.2, 1 - abs(y - 2040) / 16) for y in YEARS]
    s = sum(w)
    return [x / s for x in w]


def rev_req(cap, profile="even", life=LIFE, roe=ROE, eq=EQUITY, kd=KD, tax=TAX,
            other=OTHER, infl=0.0):
    w = weights(profile)
    capex = {y: cap * wi for y, wi in zip(YEARS, w)}
    rate = kd * (1 - eq) + eq * roe / (1 - tax) - infl
    out, prev_net = {}, 0.0
    for t in YEARS:
        dep = sum(c / life for v, c in capex.items() if v <= t and (t - v) < life)
        accdep = sum(c * min(t - v + 1, life) / life for v, c in capex.items() if v <= t)
        gross = sum(c for v, c in capex.items() if v <= t)
        net = gross - accdep
        avg = (prev_net + net) / 2
        out[t] = dep + avg * (rate + other)
        prev_net = net
    return out


def cents(rr_b, mult):                                   # $B / TWh = cents/kWh x 10 -> /10
    return rr_b / (SALES_TWH * mult) * 100


res = {}
cases = {"A High, ratepayers pay all": CAP_HIGH,
         "B Low, ratepayers pay all": CAP_LOW,
         "C Low, developers pay 50% of hosting": CAP_LOW_RP}

print(f"Capital: A {CAP_HIGH:.2f}  B {CAP_LOW:.2f}  C {CAP_LOW_RP:.2f} ($B, 2021 dollars)")
print(f"Pre-tax return rate (nominal) = {KD*(1-EQUITY)+EQUITY*ROE/(1-TAX):.4f}; real variant subtracts {INFL}")

print("\n== Reconciliation with EPRI (p.30): ~$4B/yr high, ~$1.4B/yr low by 2050 ==")
for label, infl in (("nominal return on constant-$ plant", 0.0), ("real-adjusted (return - 2.5% inflation)", INFL)):
    a = rev_req(CAP_HIGH, infl=infl)[2050]
    c = rev_req(CAP_LOW_RP, infl=infl)[2050]
    b = rev_req(CAP_LOW, infl=infl)[2050]
    print(f"{label:<44} A ${a:.2f}B (EPRI 4.0)   B ${b:.2f}B   C ${c:.2f}B (EPRI 1.4)")
    res.setdefault("reconcile", {})[label] = {"A": round(a, 2), "B": round(b, 2), "C": round(c, 2)}

print("\n== Annual revenue requirement, even spend, nominal variant ($B) ==")
for k, cap in cases.items():
    r = rev_req(cap)
    print(f"{k:<40} 2030 {r[2030]:.2f}   2040 {r[2040]:.2f}   2050 {r[2050]:.2f}")
    res.setdefault("rr_even_nominal", {})[k] = {y: round(r[y], 2) for y in (2030, 2040, 2050)}

print("\n== 2050 bill impact (nominal variant), household fixed at 600 kWh/month ==")
print(f"{'Case':<40}{'sales x':>8}{'cents/kWh':>11}{'$/month':>9}{'% of R-1 bill':>15}{'% of dist. charge':>19}")
bill = {}
for k, cap in cases.items():
    r50 = rev_req(cap)[2050]
    for m in MULTS:
        c = cents(r50, m)
        d = c / 100 * HH_KWH_MONTH
        print(f"{k:<40}{m:>8.1f}{c:>11.2f}{d:>9.2f}{d / R1_BILL * 100:>14.0f}%{c / DIST_C * 100:>18.0f}%")
        bill[f"{k} | x{m}"] = {"cents_per_kwh": round(c, 2), "usd_per_month": round(d, 2)}
res["bill_2050_nominal"] = bill

print("\n== Same, real-adjusted variant ==")
for k, cap in cases.items():
    r50 = rev_req(cap, infl=INFL)[2050]
    print(f"{k:<40}" + "   ".join(f"x{m}: {cents(r50, m):.2f}c = ${cents(r50, m) / 100 * HH_KWH_MONTH:.0f}/mo" for m in MULTS))

print("\n== Baseline context (paper fn 38, 2020 dollars) ==")
for y, v in (("2020", 2.8), ("2050 (CECP)", 6.3)):
    print(f"Distribution revenue requirement {y}: ${v}B  -> {v / SALES_TWH * 100:.1f} c/kWh at 2024 sales (all sectors)")
print("CECP increase to 2050: +$3.5B; EPRI capacity-only increment: +$1.4B to +$4.0B")

print("\n== Sensitivity of 2050 revenue requirement (case A even-spend nominal, and case C) ==")
tests = [
    ("ROE 9.0% / 10.5%", dict(roe=0.09), dict(roe=0.105)),
    ("Cost of debt 4% / 6%", dict(kd=0.04), dict(kd=0.06)),
    ("Equity share 48% / 53%", dict(eq=0.48), dict(eq=0.53)),
    ("Book life 50 / 30 yrs", dict(life=50), dict(life=30)),
    ("Other carrying 0% / 2%", dict(other=0.0), dict(other=0.02)),
    ("Spend profile mid-peak / ramp", dict(profile="mid"), dict(profile="ramp")),
    ("Inflation-adjusted / not", dict(infl=INFL), dict()),
]
print(f"{'Parameter (low / high)':<34}{'A low':>8}{'A high':>8}{'C low':>8}{'C high':>8}")
sens = []
for name, lo, hi in tests:
    row = [rev_req(CAP_HIGH, **lo)[2050], rev_req(CAP_HIGH, **hi)[2050],
           rev_req(CAP_LOW_RP, **lo)[2050], rev_req(CAP_LOW_RP, **hi)[2050]]
    print(f"{name:<34}" + "".join(f"{v:>8.2f}" for v in row))
    sens.append({"param": name, "A": [round(row[0], 2), round(row[1], 2)], "C": [round(row[2], 2), round(row[3], 2)]})
res["sensitivity_2050"] = sens

with open("step2-results.json", "w") as f:
    json.dump(res, f, indent=2)
print("\nwrote step2-results.json")