#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ rate_leak_montecarlo.py -- "What are the odds?" model for the price-stabilization value of the Okanogan PUD wireless network. ================================================================================ v1 (2026-06-22). WHY THIS MODEL EXISTS --------------------- The hardest thing to argue at the Board table is a NEGATIVE: "the towers prevent price increases that therefore never show up in anyone's bill." You cannot point to a harm that didn't happen. So instead of trying to prove the negative, this model asks the question the right way round: Across thousands of plausible 10-year futures, in what FRACTION does the network's price-discipline effect save the county MORE than the upgrade costs? That reframes "prove it won't happen" into "here is the probability distribution of what happens" -- which is exactly what a Monte Carlo is for. WHAT IT MODELS -------------- Two parallel worlds over a 10-year horizon (the Tarana gear's life), for the county's incorporated-town households: World A "No tower defense" -- the public wireless option is gone, so the private/ISP market is left to raise prices on its documented schedule. World B "With tower defense" -- a credible, already-operating public option disciplines the market, so prices rise more slowly (the yardstick / contestable-market effect). county LEAK = extra dollars households pay above today's bill in World A (the money that leaves the county when prices run). discipline BENEFIT = World A - World B (the leak the network holds back = the part attributable to the tower's presence). The decision question: P(BENEFIT > upgrade cost). INPUTS (every distribution is sourced; ranges are deliberately conservative) ----------------------------------------------------------------------------- HOUSEHOLDS H = 6,700 Okanogan incorporated-town households (WA OFM April-2024: 16,540 town residents / ~2.45 persons per household). CONSERVATIVE -- it EXCLUDES the ~11,000 rural households the wireless network serves most directly. BASE BILL B0 = $80/mo ~ Consumer Reports 2022 median internet bill ($74.99), rounded up modestly. HORIZON T = 10 years (Tarana's own planning horizon; wireless gear life). World A annual increase (% of the current avg bill), drawn EACH year: triangular(2%, 4.5%, 8%). Sources: monopoly markets raise rates ~5-8%/yr, ~2x competitive markets (ISP market-concentration analyses); 63% of households saw an increase last year averaging +$195/yr (CNET). Mode 4.5% blends magnitude with the ~37% who see no increase in a given year -- i.e. BELOW the headline monopoly rate. World B annual increase (% of current avg bill), drawn EACH year: triangular(0.5%, 2%, 4%). Competitive/disciplined markets rise roughly half as fast (Harvard Berkman Klein: community networks cheaper in 23 of 27 markets and price-stable; Consumer Reports: 3+ competitors ~ $5/mo less; competitive states 25-35% cheaper than monopoly states). PI_FAIL = 0.15 -- the probability that, in a given future, the public option delivers NO price discipline at all (World B = World A). This is an HONEST allowance against the mixed evidence (e.g. Landgraf, Telecommunications Policy 2023, found a single DSL incumbent did not improve on the mere THREAT of municipal entry). Including failures makes the headline LOWER and more credible -- it is not a sure thing in every world. No discounting and no inflation on the $1.2M cost: both are left nominal. This is conservative -- the upgrade is spent once, now, while the leak recurs and would compound with inflation. NOT A FORECAST. These are scenario distributions, not predictions of any one household's bill. The point is the SHAPE of the risk, not a point estimate. """ import json import numpy as np SEED = 20260622 # fixed -> anyone re-running gets identical numbers rng = np.random.default_rng(SEED) N = 200_000 # simulated futures T = 10 # years H = 6_700 # incorporated-town households (WA OFM) B0 = 80.0 # baseline monthly bill ($) UPGRADE_FULL = 1_180_000 # 7-site Tarana upgrade (District workbooks, ~$1.2M) UPGRADE_ONE = 118_000 # one mid-range tower (Number Hill ~$113K, Nortons ~$125K) PI_FAIL = 0.15 # P(no price discipline at all, in a given future) def run(pi_fail=PI_FAIL, a_mode=0.045, b_mode=0.02, seed=SEED): r = np.random.default_rng(seed) billA = np.full(N, B0) billB = np.full(N, B0) leak = np.zeros(N) # cumulative $ above today's bill, World A (county-wide) benefit = np.zeros(N) # cumulative $ saved by discipline = A - B (county-wide) fails = r.random(N) < pi_fail for _ in range(T): gA = r.triangular(0.02, a_mode, 0.08, N) gB = r.triangular(0.005, b_mode, 0.04, N) billA = billA * (1.0 + gA) billB = np.where(fails, billA, billB * (1.0 + gB)) leak += (billA - B0) * 12.0 * H benefit += (billA - billB) * 12.0 * H return leak, benefit leak, benefit = run() def pct(a, q): return float(np.percentile(a, q)) p_benefit_full = float((benefit > UPGRADE_FULL).mean()) p_benefit_one = float((benefit > UPGRADE_ONE).mean()) p_benefit_pos = float((benefit > 0).mean()) p_leak_full = float((leak > UPGRADE_FULL).mean()) # Histogram of the BENEFIT distribution ($M), clipped at the 99th pct for a readable chart. hi = pct(benefit, 99) / 1e6 counts, edges = np.histogram(benefit / 1e6, bins=32, range=(0.0, hi)) out = { "model": "rate_leak_montecarlo v1", "seed": SEED, "n_trials": N, "years": T, "households": H, "base_bill_mo": B0, "pi_fail": PI_FAIL, "upgrade_full": UPGRADE_FULL, "upgrade_one": UPGRADE_ONE, "p_benefit_gt_full_upgrade": round(p_benefit_full, 4), "p_benefit_gt_one_tower": round(p_benefit_one, 4), "p_benefit_positive": round(p_benefit_pos, 4), "p_leak_gt_full_upgrade": round(p_leak_full, 4), "benefit_M": {"p10": round(pct(benefit,10)/1e6,3), "p50": round(pct(benefit,50)/1e6,3), "p90": round(pct(benefit,90)/1e6,3), "mean": round(float(benefit.mean())/1e6,3)}, "leak_M": {"p10": round(pct(leak,10)/1e6,3), "p50": round(pct(leak,50)/1e6,3), "p90": round(pct(leak,90)/1e6,3), "mean": round(float(leak.mean())/1e6,3)}, "per_household_10yr_leak_p50": round(pct(leak,50)/H, 0), "benefit_hist": {"edges_M": [round(e,3) for e in edges.tolist()], "counts": counts.tolist()}, } # --- Sensitivity: is the headline robust to the failure rate and the World-A mode? --- sens = [] for pf in (0.10, 0.15, 0.25): for am in (0.035, 0.045, 0.06): _, ben = run(pi_fail=pf, a_mode=am, seed=SEED + int(pf*100) + int(am*1000)) sens.append({"pi_fail": pf, "a_mode": am, "p_benefit_gt_full": round(float((ben > UPGRADE_FULL).mean()), 3), "benefit_p50_M": round(float(np.percentile(ben,50))/1e6, 2)}) out["sensitivity"] = sens print(json.dumps(out, indent=2))