import numpy as np, importlib.util, sys, os
spec = importlib.util.spec_from_file_location("mc", os.path.join(os.path.dirname(__file__),"montecarlo_robust.py"))
mc = importlib.util.module_from_spec(spec); spec.loader.exec_module(mc)

print("="*70); print("1) SEED SENSITIVITY of the DISTRICT cost-recovery headline (default strict cliff)")
for tag,scn in [("flat_rate_floor",{}),("recovers_half",{'rate_track':0.5}),("recovers_costs",{'rate_track':1.0})]:
    vals=[]
    for s in range(20260629,20260629+8):
        o=mc.cost_recovery_mc(seed=s,scn=scn,verbose=False)
        vals.append(o['p_payback'])
    vals=np.array(vals)
    print(f"  {tag:16s} mean={vals.mean():.2f}  sd={vals.std(ddof=1):.3f}  min={vals.min():.1f} max={vals.max():.1f}  range={vals.max()-vals.min():.2f}  (reported SE~0.1)")

print("="*70); print("2) SEED SENSITIVITY of county p_county_ahead (10yr) equal-anchor")
vals=[]
for s in range(100,108):
    r=mc.county_value_mc(10,seed=s); vals.append(r['p_county_ahead'])
vals=np.array(vals); print(f"  county10 ahead mean={vals.mean():.2f} sd={vals.std(ddof=1):.3f} range={vals.max()-vals.min():.2f}")

print("="*70); print("3) CONVERGENCE: p_payback vs n (recovers_costs)")
for n in [2000,20000,200000,1000000]:
    o=mc.cost_recovery_mc(n=n,scn={'rate_track':1.0},verbose=False)
    print(f"  n={n:>8}  p={o['p_payback']:.2f}  SE={o['se']}")

print("="*70); print("4) TRY TO REPRODUCE scenarios_output.json from the JS calculator scen configs")
# JS SCEN mapped to python scn keys (Sticky/Comm div 100, Star/Shock div 100, Cov div 1)
scen_cost = {
 "pessimistic": dict(scn={'starlink':0.02,'shock':0.20,'cliff':True,'rate_track':0.0}),
 "balanced":    dict(scn={'rate_track':0.6,'cliff':False,'sticky':0.002,'community':0.002}),
 "optimistic":  dict(scn={'rate_track':1.0,'cliff':False,'coverage':180,'sticky':0.006,'community':0.009}),
}
for k,kw in scen_cost.items():
    o=mc.cost_recovery_mc(headline_arpu=mc.ARPU_YR,verbose=False,**kw)
    print(f"  {k:12s} district p_payback={o['p_payback']:.1f}   (json says pess26.5/bal85.9/opt95.3)")

scen_cty = {
 "pessimistic": {'starlink':0.02,'shock':0.20,'premium':0.06,'rate_track':0.0},
 "balanced":    {'rate_track':0.6,'premium':0.12,'sticky':0.002,'community':0.002},
 "optimistic":  {'rate_track':1.0,'premium':0.22,'coverage':180,'sticky':0.006,'community':0.009},
}
for k,scn in scen_cty.items():
    r10=mc.county_value_mc(10,scn=scn); r15=mc.county_value_mc(15,scn=scn)
    print(f"  {k:12s} county10={r10['p_county_ahead']:.1f} county15={r15['p_county_ahead']:.1f} cty_P50_15={r15['county_P50']} hh15={r15['household_savings_mean'] if False else r15['household_savings_P50']}")

print("="*70); print("5) Is 'macro-stress only loads the downside' actually conservative? distribution of Zp")
rng=np.random.default_rng(1)
Z=rng.standard_normal(200000); Zp=np.clip(Z,0,None)
print(f"  E[Z]=0 by construction; E[Zp]={Zp.mean():.3f} (>0). So clipping ADDS a positive mean shock to churn/overbuild and subtracts from growth -> net pessimistic. Mean churn add = {0.0015*Zp.mean():.5f}/yr")