The Arbitrage Value of Snowy 2.0

A backtest of storage arbitrage on NSW spot prices, 2022 to 2025

Run this yourself.

Produced with help from Claude Code. All method choices, data, and conclusions are the author’s own.

TipSummary

In this report, we backtest energy arbitrage for a plant with Snowy 2.0’s published specifications against four years of real half-hourly NSW spot prices (2022 to 2025). We consider three dispatch schemes that provide bounds on the profitability of the strategy. Firstly, we consider a perfect-foresight case where the plant knows the future prices and can dispatch optimally. Secondly, we consider a rule of thumb where the plant dispatches on two fixed price triggers that say when to charge and discharge. Finally, we consider a forecaster which dispatches on a statistical model of prices, trained on the first three years and tested on the fourth. Our results show that knowing the distribution of prices is almost as good as knowing the actual prices, capturing 84 to 93% of the perfect-foresight profits in the three years after the 2022 fuel crisis. Note, however, that this analysis assumes that the plant is a price taker, and does not account for the fact that a plant of this size will move prices against itself. Hence, the profit estimates should be read as optimistic rather than fair. A financial appendix then prices the backtest as a project: at the official $12 billion capital cost the venture returns about 4.85% real after tax, so the discount rate, not the operating strategy, decides whether it pays.

Introduction

Snowy 2.0 is a 2,200 megawatt pumped-hydro plant under construction in New South Wales. Like a giant battery, it can earn money through arbitrage: store energy when power is cheap, sell it back to the grid when power is expensive. This report estimates what that trade is worth.

We answer this with a backtesting approach using historical prices and asking how the plant would have traded against them. Specifically, we feed four years of real half-hourly NSW spot prices (2022 to 2025) through a Snowy-sized plant, year by year, and measure what arbitrage could have paid under three operating schemes:

1. a perfect foresight scheme where the operator knows every future price;
2. a simple rule of thumb scheme with two price triggers (buy when cheap, sell when dear);
3. a forecaster scheme where the operator knows the statistical pattern of prices but never the actual value.

Each case answers the same question with a different amount of knowledge.

The perfect foresight operator is the theoretical best case: for each year it knows the exact sequence of prices and can time every half-hour against it, even holding charge for a spike it sees coming. So, this scheme sets the ceiling on the arbitrage profit that year could have allowed.

The rule of thumb operator knows only price levels, not their timing: it pumps whenever the price is among the year’s cheapest and generates whenever it is among the most expensive. We allow this operator to use the data from the first three years to set its two price triggers, but it never sees the actual 2025 prices. This scheme is a simple, mechanical benchmark: it needs no forecasting at all, and it sets the floor any smarter operator should beat.

The forecaster is probably the most realistic operator. It learns the typical time-of-day shape of prices from 2022 to 2024, commits to a single standing policy, and is then judged on the untouched 2025; that out-of-sample result is the report’s headline.

As such, these three cases let us explore where the money comes from. The difference in profits between the rule of thumb and the forecaster scheme pins down the worth of skill, and that between the forecaster and the perfect foresight scheme pins down what a perfect price forecast would have been worth (the perfect-foresight gap). The forecaster’s profit is the main result of the report.

Note

Profit throughout means arbitrage profit: revenue from selling electricity, minus the cost of buying it, minus the variable cost of running the machines. Fixed operating costs, capital costs, and financing sit outside the dispatch problem and are not modelled.

Measuring only energy-arbitrage profit leaves the estimate incomplete in two opposite ways. It under-counts revenue, because frequency-control (FCAS) markets and any capacity or reliability payments are separate income streams we leave out entirely. And it over-counts trading skill, because we assume no outages and no network limits, and, above all, we treat the plant as a price taker: able to buy and sell any amount without moving the price. That is false for a plant this large. When 2,200 MW sells into an expensive evening it pushes that price down, eating into the very gap it is chasing, so every number here over-states what the plant would really earn, and the results should be read as bounds rather than forecasts.

The rest of the report is ordered as follows: we start off with a short background of the plant and the market in which it trades, then discuss the data, illustrate the arbitrage trade on one real day, and finally set out the dispatch model and its results. The appendix contains a glossary, the full assumption audit, a hindsight benchmark for the fixed rule, a first pass at price impact, the financial analysis, and data notes.

Background and Data

Snowy 2.0

Snowy 2.0 links two existing reservoirs in the Snowy Mountains, Tantangara (upper) and Talbingo (lower), through about 27 kilometres of tunnels and a new underground power station. When power is cheap it pumps water uphill; when power is expensive the water runs back down through the turbines. Snowy Hydro rates the plant at 2,200 MW of capacity with a headline storage figure of about 350 GWh, roughly 160 hours of generation at full power (Snowy Hydro, About Snowy 2.0). It is Australia’s largest renewable-energy project, budgeted at about $12 billion after several revisions, with full operation expected around 2028 (project overview; SBS News).

The 350 GWh figure is contested. Critics argue it describes the energy released by draining the upper reservoir once, while the storage that can actually be cycled (pumped back up and generated again, repeatedly) may be somewhere between about 40 and 230 GWh (The Conversation). Snowy Hydro and others dispute that reading (SBS News). We take the 350 GWh at face value, and show in the Analysis that the dispute barely moves the answer: arbitrage value comes from daily cycling, not from the size of the store.

NSW spot market

Once operational, Snowy 2.0 will sell into the National Electricity Market (NEM), where a wholesale spot price is set every five minutes for each region. There is a consistent intra-day price swing in that market: renewable generation floods the grid during the middle of the day, pushing prices down, while demand peaks after sunset just as solar disappears, pushing prices up. Those daily swings are the whole business case for storage, as they create the opportunity to buy low in the afternoon and sell high in the evening.

Note

The NEM has a price cap and a price floor. The floor is fixed at -$1,000/MWh. The cap is indexed each July and rose across our window: $15,100/MWh until June 2022, $15,500 from July 2022, $16,600 from July 2023, $17,500 from July 2024, and $20,300 from July 2025.

Data: NSW spot prices, 2022 to 2025

We load four full years of the New South Wales spot price from the market operator (AEMO), averaged from its five-minute settlements to half-hourly steps. A saved snapshot of each year keeps every run identical (details in the appendix under Data notes).

# On a fresh cloud machine (for example Colab) the notebook arrives on its own: no pinned
# packages and no price snapshots. This cell makes the notebook self-bootstrapping in that case,
# and does nothing when you already have the repository checked out locally. Everything the report
# needs is defined in this one notebook, so there is no companion module to fetch.
import importlib.util, subprocess, sys, urllib.request
from pathlib import Path

RAW = "https://raw.githubusercontent.com/yobin-tim/snowy-arbitrage/main/python"   # the published repository

def _fetch(relpath: str, dest: Path):
    """Download python/<relpath> from the repository to dest, unless dest already exists."""
    if dest.exists():
        return
    dest.parent.mkdir(parents=True, exist_ok=True)
    urllib.request.urlretrieve(f"{RAW}/{relpath}", str(dest))

# 1. Packages. On a fresh cloud machine (for example Colab) a few packages are missing; install
#    ONLY those, at the pinned versions from requirements.txt. We deliberately leave the machine's
#    preinstalled scientific stack (numpy, pandas, scipy, matplotlib) untouched: the running kernel
#    has already loaded those libraries, and replacing them on disk mid-session leaves the kernel
#    half old and half new until a restart. The result checks throughout the report verify that the
#    numbers come out the same either way.
_missing = [p for p in ("cvxpy", "highspy", "pyarrow", "nemosis") if importlib.util.find_spec(p) is None]
if _missing:
    _req = Path("requirements.txt")
    if not _req.exists():
        try:
            _fetch("requirements.txt", _req)          # pull the pinned list onto a bare cloud machine
        except Exception:
            pass                                      # no list available; fall back to unpinned names
    _pins = {}                                        # package name -> its "name==version" pin
    if _req.exists():
        for _line in _req.read_text().splitlines():
            _spec = _line.split("#")[0].strip()       # drop comments and whitespace
            if "==" in _spec:
                _pins[_spec.split("==")[0].strip()] = _spec
    _install = [_pins.get(p, p) for p in _missing]    # pin each missing package where known
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", *_install])

from dataclasses import dataclass, replace
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cvxpy as cp
from IPython.display import Markdown, display   # lets us print sentences filled in with live numbers

DATA_DIR = Path("data") if Path("data").exists() else Path("python/data")   # works from python/ or repo root
CACHE    = DATA_DIR.parent / ".nemosis_cache"       # raw AEMO downloads (git-ignored)
FIG_DIR  = DATA_DIR.parent / "figures"              # every figure is saved here as a PNG for reuse
FIG_DIR.mkdir(parents=True, exist_ok=True)

# 2. Price snapshots. Fetch the committed parquet for each year so the report reproduces its exact
#    numbers offline. Without them a cloud run would fall back to a slow live AEMO download instead.
for _year in (2022, 2023, 2024, 2025):
    _snap = DATA_DIR / f"nsw1_prices_{_year}_30min.parquet"
    if not _snap.exists():
        try:
            _fetch(f"data/{_snap.name}", _snap)
        except Exception:
            pass          # not fatal: the loader below will pull that year live from AEMO if needed

# One shared plotting style so every figure reads as a set: bold black data lines, light fills,
# a light grid, and a full border around each plot. Colours come from a validated, colour-vision-
# safe palette, and each colour keeps one fixed meaning throughout.
INK, INK2, MUTED = "#0b0b0b", "#52514e", "#898781"
GRID, SURFACE    = "#e1e0d9", "#fcfcfb"
BUY   = "#2a78d6"   # blue   -> pumping / buying / charging
SELL  = "#eb6834"   # orange -> generating / selling / discharging
STORE = "#1baf7a"   # green  -> state of charge (how full the store is)
SPIKE = "#d03b3b"   # red    -> extreme price spikes (the fragile part of the revenue)

plt.rcParams.update({
    "figure.figsize": (9, 4.2), "figure.dpi": 120,
    "axes.facecolor": SURFACE, "figure.facecolor": SURFACE,
    "axes.edgecolor": MUTED, "axes.linewidth": 1.0,   # a full border on every plot
    "axes.labelcolor": INK, "text.color": INK, "xtick.color": INK2, "ytick.color": INK2,
    "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.8,
    "axes.spines.top": True, "axes.spines.right": True,   # keep the graph bordered
    "font.size": 11, "axes.titlesize": 13, "axes.titleweight": "bold", "axes.titlepad": 10,
})

def dollars(x, _=None):
    '''Format an axis tick as a whole-dollar label: 1500 -> $1,500.'''
    return f"${x:,.0f}"

def clock(h):
    '''Format an hour-of-day number as clock time: 0 -> 12am, 15 -> 3pm, 24 -> 12am.'''
    h = int(round(h)) % 24
    return "12am" if h == 0 else "12pm" if h == 12 else (f"{h}am" if h < 12 else f"{h-12}pm")

def clock_axis(ax, last=24, step=3):
    '''Label an hour-of-day x axis in clock time, midnight to midnight.'''
    ticks = list(range(0, last + 1, step))
    ax.set_xticks(ticks)
    ax.set_xticklabels([clock(t) for t in ticks])
    ax.set_xlim(0, last)

def finish(fig, name):
    '''Save the figure as a PNG under python/figures/, then show it here.'''
    fig.savefig(FIG_DIR / name, dpi=200, bbox_inches="tight", facecolor=SURFACE)
    plt.show()

def ensure_year_prices(year: int) -> pd.Series:
    '''Return the 30-minute NSW1 price series for `year`. Uses the saved snapshot if it exists,
    otherwise downloads that year from AEMO, tidies it, saves it, and returns it. Idempotent:
    safe to call every run, it downloads at most once per year.'''
    snapshot = DATA_DIR / f"nsw1_prices_{year}_30min.parquet"
    if snapshot.exists():
        return pd.read_parquet(snapshot)["price"]

    # ---- first-time download path ----
    from nemosis import dynamic_data_compiler
    DATA_DIR.mkdir(parents=True, exist_ok=True)
    CACHE.mkdir(parents=True, exist_ok=True)
    raw = dynamic_data_compiler(
        f"{year}/01/01 00:00:00", f"{year+1}/01/01 00:00:00", "DISPATCHPRICE",
        raw_data_location=str(CACHE),
        select_columns=["SETTLEMENTDATE", "REGIONID", "INTERVENTION", "RRP"],
        filter_cols=["REGIONID"], filter_values=[["NSW1"]], fformat="parquet")
    raw = raw[raw["INTERVENTION"] == 0].copy()                 # drop duplicate intervention-pricing rows
    raw["SETTLEMENTDATE"] = pd.to_datetime(raw["SETTLEMENTDATE"])
    keep = (raw["SETTLEMENTDATE"] > f"{year}-01-01") & (raw["SETTLEMENTDATE"] <= f"{year+1}-01-01")
    raw = raw.loc[keep]                                        # clip to exactly one calendar year
    raw["interval_start"] = raw["SETTLEMENTDATE"] - pd.Timedelta(minutes=5)   # label by interval start
    # Average the six 5-minute settlements inside each half hour into one 30-minute price.
    prices = (raw.sort_values("interval_start").set_index("interval_start")["RRP"]
                 .astype(float).resample("30min").mean().rename("price"))
    prices.to_frame().to_parquet(snapshot)
    return prices

YEARS = [2022, 2023, 2024, 2025]
prices_by_year = {year: ensure_year_prices(year) for year in YEARS}
CLOSE_YEAR = 2025                                             # close-up year: the latest, and later the held-out test year
price_close = prices_by_year[CLOSE_YEAR]
print(f"Loaded {', '.join(map(str, YEARS))}: "
      f"{sum(len(s) for s in prices_by_year.values()):,} half-hourly prices in total.")

Loaded 2022, 2023, 2024, 2025: 70,128 half-hourly prices in total.

p = price_close                      # short alias for the close-up year, reused by the figure cells below
share_negative = (p < 0).mean() * 100
display(Markdown(
    f"In {CLOSE_YEAR}, the NSW spot price averaged **\\${p.mean():,.0f}/MWh**, ranged from "
    f"**\\${p.min():,.0f}** to **\\${p.max():,.0f}/MWh**, and was "
    f"**below zero {share_negative:.0f}%** of the time."
))

In 2025, the NSW spot price averaged $103/MWh, ranged from $-588 to $13,849/MWh, and was below zero 13% of the time.

Before we set up the dispatch model, we look at the raw prices in three ways: the average daily cycle, the distribution of prices within a year, and whether the opportunity persists from year to year. The close-ups use 2025, the most recent year of data and, later in the report, the held-out test year; the final figure compares all four years.

Daily price cycle

The average 2025 day makes the cycle visible: a small peak in the morning, followed by prices sagging through the middle of the day and climbing into the evening. The shaded blocks mark the cheapest and most expensive four-hour windows, which are roughly when a storage plant would buy and sell.

# Step 1. The average day: the close-up year's mean price at each of the 48 half-hour slots of the day.
average_by_time = p.groupby(p.index.hour + p.index.minute/60).mean()

# Step 2. Find the two four-hour windows to shade. We use CONTIGUOUS blocks (8 half-hours) so
# the shading reads as one window each: the cheapest consecutive block anywhere in the day, and
# the most expensive consecutive block after 2 pm (the evening peak, where a storage plant sells).
window = 8
running = np.convolve(average_by_time.values, np.ones(window), "valid")   # sum of each 8-slot block
cheap_start = running.argmin()
cheapest_hours = average_by_time.index[cheap_start:cheap_start + window]
afternoon = average_by_time[average_by_time.index >= 14]
peak_start = np.convolve(afternoon.values, np.ones(window), "valid").argmax()
evening_peak = afternoon.index[peak_start:peak_start + window]

# Step 3. Draw: the average-day curve, the two shaded windows, and one annotation for each.
fig, ax = plt.subplots()
ax.plot(average_by_time.index, average_by_time.values, color=INK, lw=2.2)
ax.fill_between(average_by_time.index, average_by_time.values,
                where=average_by_time.index.isin(cheapest_hours),
                color=BUY, alpha=0.25, label="Cheapest 4 hours")
ax.fill_between(average_by_time.index, average_by_time.values,
                where=average_by_time.index.isin(evening_peak),
                color=SELL, alpha=0.30, label="Most expensive 4 hours")
trough_centre = cheapest_hours[window // 2]            # point the arrow into the shaded trough
ax.annotate("Renewables flood midday",
            xy=(trough_centre, average_by_time.loc[trough_centre]),
            xytext=(7.5, average_by_time.max()*0.45),
            color=BUY, fontsize=9, arrowprops=dict(arrowstyle="->", color=BUY))
ax.annotate("Evening peak",
            xy=(average_by_time.idxmax(), average_by_time.max()),   # the true daily peak
            xytext=(11.5, average_by_time.max()*0.9),
            color=SELL, fontsize=9, arrowprops=dict(arrowstyle="->", color=SELL))
ax.set(title=f"Average Daily Price, NSW ({CLOSE_YEAR})", xlabel="Time of day", ylabel="Price ($/MWh)")
clock_axis(ax)
ax.yaxis.set_major_formatter(dollars)
ax.legend(frameon=False)
finish(fig, "01-average-day.png")

Price distribution

We sort every half-hour of the year from cheapest to most expensive and plot the same curve twice below. Sorted this way, the curve is the year’s cumulative price distribution: pick any share of the year on the horizontal axis and the curve gives the price below which that share of the year traded. (Formally it is the CDF drawn with its axes swapped, which keeps the extreme prices readable.)

The top panel shows the full range of prices in 2025; a handful of extreme half-hours at the far right dominate the view. The bottom panel redraws the same curve with the axis clipped at $400/MWh, and the shape of the ordinary year becomes visible: prices below zero for the first 13% of the year, half the year below $79/MWh, and a gentle climb through the modest band where almost the whole year sits.

As such, prices are modest for most of the year and extreme only occasionally; the Analysis section shows how much of the profit rides on those few extreme hours.

# Step 1. Sort the year cheapest-first. Against the share of the year on the x axis, the sorted
# curve IS the cumulative price distribution (the CDF with its axes swapped): at x% of the year,
# the curve's height is the price below which x% of the half-hours traded.
sorted_prices = np.sort(p.values)
share_of_year = np.arange(1, len(sorted_prices) + 1) / len(sorted_prices) * 100
top_1pct = int(0.01 * len(sorted_prices))               # how many half-hours make up the top 1%
neg_share = (sorted_prices < 0).mean() * 100            # where the curve crosses zero
median_price = float(np.median(sorted_prices))

# Step 2. Draw the same curve twice; the panels differ only in the range of the price axis.
fig, (full, zoom) = plt.subplots(2, 1, figsize=(9, 5.8), height_ratios=[1, 1.25])

# Top panel: the full range. The top 1% is shaded and the single most expensive half hour
# marked, because at this scale the extremes are the only feature the eye can read.
full.plot(share_of_year, sorted_prices, color=INK, lw=2.2)
full.fill_between(share_of_year[-top_1pct:], sorted_prices[-top_1pct:], color=SPIKE, alpha=0.35)
full.plot(share_of_year[-1], sorted_prices[-1], marker="o", color=SPIKE, ms=7)
full.annotate(f"The most expensive half-hour\nreaches ${sorted_prices[-1]:,.0f}/MWh;\nthe shaded band is the top 1%",
              xy=(share_of_year[-1], sorted_prices[-1]), xytext=(58, sorted_prices[-1]*0.5),
              color=SPIKE, fontsize=9, arrowprops=dict(arrowstyle="->", color=SPIKE))
full.set(title=f"The {CLOSE_YEAR} Prices as a Cumulative Distribution", ylabel="Price ($/MWh)")
full.yaxis.set_major_formatter(dollars)

# Bottom panel: the same curve with the axis clipped to the everyday band, so its shape shows.
# The spikes run off the top on purpose; the two dotted markers are the zero crossing and the
# median, each labelled with its cumulative reading.
zoom.plot(share_of_year, sorted_prices, color=INK, lw=2.2)
zoom.axhline(0, color=MUTED, lw=1, ls=":")              # the zero-price line
below = sorted_prices < 0
zoom.fill_between(share_of_year[below], sorted_prices[below], 0, color=BUY, alpha=0.30)
zoom.annotate(f"{neg_share:.1f}% of the year\nis priced below zero", xy=(neg_share, 0),
              xytext=(neg_share + 6, -72), color=BUY, fontsize=9,
              arrowprops=dict(arrowstyle="->", color=BUY))
zoom.plot([50, 50], [-100, median_price], color=MUTED, lw=1, ls=":")
zoom.plot(50, median_price, marker="o", color=INK, ms=5)
zoom.annotate(f"half the year is below ${median_price:,.0f}/MWh", xy=(50, median_price),
              xytext=(52, 190), color=INK2, fontsize=9, arrowprops=dict(arrowstyle="->", color=MUTED))
zoom.set(ylim=(-100, 400), xlabel="Share of the year (%)", ylabel="Price ($/MWh)")
zoom.yaxis.set_major_formatter(dollars)
zoom.text(30, 320, "Almost the whole year sits in this modest band", color=INK2, fontsize=9.5, ha="center")
finish(fig, "02-sorted-half-hours.png")

The opportunity, year by year

The following graphs show that there is a market opportunity for storage arbitrage.

The left panel shows the spread of daily price ranges (each day’s high minus its low), not just the median: the median day’s range is $210 to $245/MWh in every year, and even the 10th-percentile day clears about $80 to $130/MWh, several times what the round-trip loss costs on a typical cycle. In other words, cycling pays on almost every day of the year, not only on the rare spectacular ones.

The right panel shows the share of half-hours priced below zero, which more than quadruples across the window as solar capacity grows. So the opportunity is structural and increasing over the four years that we have considered in this report.

# Two summaries per year: the DISTRIBUTION of daily price ranges (the size of the everyday
# opportunity, shown as a median with two bands rather than a single bar), and the share of
# half-hours priced below zero (the cheapening buy side), labelled with its value.
spread_q, negative_share = {}, []
for year in YEARS:
    series = prices_by_year[year]
    by_day = series.groupby(series.index.date)
    gaps = by_day.max() - by_day.min()                            # one range (high minus low) per day
    spread_q[year] = gaps.quantile([0.10, 0.25, 0.50, 0.75, 0.90])
    negative_share.append((series < 0).mean() * 100)              # how often the plant is paid to charge

fig, (spread_ax, negative_ax) = plt.subplots(1, 2, figsize=(9, 3.9))
for i, year in enumerate(YEARS):
    q = spread_q[year]
    spread_ax.plot([i, i], [q[0.10], q[0.90]], color=MUTED, lw=7, alpha=0.35,
                   solid_capstyle="butt", zorder=1)               # middle 80% of days
    spread_ax.plot([i, i], [q[0.25], q[0.75]], color=MUTED, lw=7, alpha=0.8,
                   solid_capstyle="butt", zorder=2)               # middle 50% of days
    spread_ax.plot(i, q[0.50], marker="_", color=INK, ms=15, mew=3, zorder=3)
    spread_ax.annotate(f"${q[0.50]:,.0f}", (i, q[0.50]), xytext=(10, -3),
                       textcoords="offset points", fontsize=8.5, color=INK)
spread_ax.set_xticks(range(len(YEARS)))
spread_ax.set_xticklabels([str(y) for y in YEARS])
spread_ax.set_xlim(-0.6, len(YEARS) - 0.4)
spread_ax.set(title="Daily Price Range", ylabel="Range ($/MWh)")
spread_ax.yaxis.set_major_formatter(dollars)
spread_ax.legend(handles=[plt.Line2D([], [], color=MUTED, lw=7, alpha=0.35, label="Middle 80% of days"),
                          plt.Line2D([], [], color=MUTED, lw=7, alpha=0.8, label="Middle 50% of days"),
                          plt.Line2D([], [], marker="_", color=INK, ls="none", ms=12, mew=3, label="Median day")],
                 frameon=False, loc="upper left", fontsize=8.5)
negative_ax.bar([str(y) for y in YEARS], negative_share, color=BUY, alpha=0.85, edgecolor=INK, linewidth=1.2)
for i, share in enumerate(negative_share):                        # the actual percentage on each bar
    negative_ax.text(i, share + 0.3, f"{share:.1f}%", ha="center", fontsize=9.5, color=INK)
negative_ax.set(title="Prices Below Zero", ylabel="Share of the year (%)",
                ylim=(0, max(negative_share) * 1.18))
finish(fig, "03-opportunity-by-year.png")

Taken together, this pins down the key reasons for an arbitrage opportunity: a daily cycle that reliably pairs cheap afternoons with expensive evenings, a thin set of extreme hours carrying outsized prices, and a buy side that gets cheaper each year.

What they do not show is how a plant converts that opportunity into profit, or what the round-trip loss takes on the way. The next section works through both on a single winter day.

Storage Arbitrage: An Illustration

Before setting out the model, we walk through the trade itself on one real day. The section fixes the two ideas everything later rests on: the mechanics, where buying cheap half-hours fills the store and selling expensive ones drains it, with the store’s level tying each decision to the ones around it; and the economics, where a price gap becomes profit only once it clears the round-trip loss.

One day’s trade (24 July 2025)

We take one real winter day and rank its 48 half-hours by price: the plant buys during the ten cheapest (blue) and sells during the ten most expensive (orange). The ranking uses the whole day at once, so it is hindsight rather than one of the three schemes; it is the trade itself at its simplest, the thing every scheme in the Analysis is trying to find. On this day the ten most expensive half-hours all sit in the evening peak, and the ten cheapest in the solar-flooded middle of the day, where the price actually goes negative: the plant is briefly paid to charge. (The smaller morning peak goes unsold; the evening beats it.) The lower panel tracks the state of charge this produces: the store starts the day empty, fills across the cheap midday hours, and drains through the evening peak, ending the day exactly as charged as it began (the same start-equals-end convention the dispatch model later imposes on each year).

# Step 1. Pick the day and rank its 48 half-hours by price, using the whole day at once.
# That is hindsight, not one of the three schemes: the point is to show the trade itself.
one_day = p.loc["2025-07-24"]
window = 10                                        # trade 10 half-hours on each side
cheap_slots = one_day.nsmallest(window).index      # when the plant buys (store fills)
pricey_slots = one_day.nlargest(window).index      # when the plant sells (store drains)

# Step 2. Track the store at each half-hour boundary: up one notch after a buy, down one
# after a sell. Lifting the whole path so it never dips below empty is the same as saying
# the store began the day partly charged; with ten buys and ten sells it also ends the day
# exactly where it started.
boundary_level = [0.0]
for t in one_day.index:
    boundary_level.append(boundary_level[-1] + (1 if t in cheap_slots else -1 if t in pricey_slots else 0))
boundary_level = np.array(boundary_level) - min(boundary_level)

# Step 3. Draw. Each half hour has one settled price, so the price is drawn as a step. The
# buy/sell columns are shaded over the FULL panel height (not just up to the price line), so
# a traded half-hour reads as a column even where the price is negative; the store line below
# then rises and falls exactly across those columns (the level is plotted at each half-hour
# boundary, so a buy interval's ramp spans that interval, not the one after it).
hour_of_day = one_day.index.hour + one_day.index.minute/60
edge_hours = np.append(hour_of_day, 24.0)          # the 49 half-hour boundaries of the day
fig, (price_ax, store_ax) = plt.subplots(2, 1, sharex=True, figsize=(9, 5.2), height_ratios=[2, 1])
price_ax.plot(hour_of_day, one_day.values, color=INK, lw=2.2, drawstyle="steps-post", zorder=3)
price_ax.fill_between(hour_of_day, 0, 1, where=one_day.index.isin(cheap_slots), step="post",
                      transform=price_ax.get_xaxis_transform(), color=BUY, alpha=0.22,
                      label="Buy (cheap)")
price_ax.fill_between(hour_of_day, 0, 1, where=one_day.index.isin(pricey_slots), step="post",
                      transform=price_ax.get_xaxis_transform(), color=SELL, alpha=0.25,
                      label="Sell (expensive)")
price_ax.set(title="Illustration: Arbitrage Action (24 July 2025)", ylabel="Price ($/MWh)")
price_ax.yaxis.set_major_formatter(dollars)
price_ax.legend(frameon=False, ncol=1, loc="upper left")
store_ax.fill_between(edge_hours, boundary_level, color=STORE, alpha=0.5)
store_ax.plot(edge_hours, boundary_level, color=STORE, lw=2.2)
store_ax.set(xlabel="Time of day", ylabel="State of charge")
store_ax.set_yticks([])
clock_axis(store_ax)
# the caveats travel with the PNG, as a note under the figure rather than inside the plot
fig.text(0.005, -0.01, "Sketch of the idea, not model output. Buys and sells are simply placed in "
         "the day's ten cheapest and ten most expensive\nhalf-hours; this shows the mechanics, "
         "not an optimised schedule.", ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "04-one-day-trade.png")

Losses and the minimum profitable spread

We sort the day’s prices from cheap to expensive. We buy across the cheap block on the left and sell across the expensive block on the right; the vertical gap between the blocks is the price gap, the raw spread per megawatt-hour traded. Storage is leaky: Snowy’s round-trip efficiency is about 77%, so of every megawatt-hour bought, only 0.77 comes back out. Breaking even therefore needs the sell price to beat the buy price by about 30% (one divided by 0.77 is about 1.30), not merely by the 23% lost on the way through. Only the spread beyond that hurdle is profit, and when the gap falls short of it, the rational move is to do nothing.

# Step 1. Sort the day from cheapest to most expensive half hour; the traded blocks sit at the
# two ends, and the gap between their average price levels is the raw spread, BEFORE losses.
day_sorted = np.sort(one_day.values)
buy_price, sell_price = day_sorted[:window].mean(), day_sorted[-window:].mean()

# Step 2. Draw: the sorted curve, the two traded blocks, and the price gap between them.
fig, ax = plt.subplots()
ax.plot(range(len(day_sorted)), day_sorted, color=INK, lw=2.2)
ax.fill_between(range(window), day_sorted[:window], color=BUY, alpha=0.30, label="Buy (cheapest hours)")
ax.fill_between(range(len(day_sorted)-window, len(day_sorted)), day_sorted[-window:],
                color=SELL, alpha=0.35, label="Sell (most expensive hours)")
ax.annotate("", xy=(len(day_sorted)-window/2, sell_price), xytext=(len(day_sorted)-window/2, buy_price),
            arrowprops=dict(arrowstyle="<->", color=SPIKE, lw=2))
ax.text(len(day_sorted)-window, (buy_price+sell_price)/2,
        f"  Price gap\n  ~${sell_price-buy_price:,.0f}/MWh",
        color=SPIKE, fontsize=10, va="center")
ax.set(title="Illustration: The Same Day, Sorted (24 July 2025)",
       xlabel="Half-hours of the day, sorted by price", ylabel="Price ($/MWh)")
ax.yaxis.set_major_formatter(dollars)
ax.legend(frameon=False, loc="upper left")
# same caveat as the previous figure, travelling with the PNG
fig.text(0.005, -0.01, "Sketch of the idea, not model output; the shaded blocks are the same ten "
         "cheapest and ten most expensive half-hours as above.",
         ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "05-gap-vs-losses.png")

That is the whole trade in miniature: fill the store when the price is low, drain it when the price is high, and act only when the spread clears the loss hurdle. What the sketch cannot say is how an operator should pick those hours without seeing the day in advance, or how the plant’s power rating and reservoir size limit the trade. The Model takes up both.

The Model

This section sets the whole problem out in one place: the plant’s parameters, the dispatch problem every scheme solves, the assumptions behind it, and the one economic idea (the water value of stored energy) that turns the physics into decisions.

Plant parameters

The model needs five primitives: the power rating \(P\), the dischargeable storage \(E_{\max}\), the round-trip efficiency \(\eta\), the variable cost \(c_{\text{om}}\), and the time step \(\Delta t\). Physical figures are Snowy 2.0’s public specifications; the cost figure is set at the high end of published estimates.

Parameter	Value	Note
Power rating \(P\)	2,200 MW	Snowy 2.0 nameplate (rated) capacity, assumed the same for pumping and generating (A6).
Dischargeable storage \(E_{\max}\)	350 GWh	Snowy’s headline figure; the cyclable capacity is contested, so we test it in Scheme 1’s storage-capacity sensitivity (A4).
Round-trip efficiency \(\eta\)	77%	Reported for the plant; split evenly across the two legs (\(\eta_{\text{pump}} = \eta_{\text{gen}} = \sqrt{\eta}\)), which does not affect revenue (A3).
Variable cost \(c_{\text{om}}\)	$4 / MWh generated	High end of the published range for pumped hydro, so conservative here.
Time step \(\Delta t\)	30 minutes	We average the 5-minute market price to half-hourly steps, which smooths sub-half-hour spikes and mildly under-states the ceiling.

Fixed operating costs and the up-front capital cost are left out on purpose: they are the same whatever the plant dispatches, so they matter for whether the project pays off overall, but not for the best operating schedule, which is all this model chooses.

@dataclass(frozen=True)
class StoragePlant:
    '''Physical description of a storage plant. Defaults are Snowy 2.0.'''
    power_mw: float = 2200.0            # pump and generate power rating (symmetric here), MW
    dischargeable_mwh: float = 350_000  # usable generation energy (Snowy's 350 GWh, taken at face value)
    eta_pump: float = 0.77 ** 0.5       # pumping-leg efficiency  (square root of the 0.77 round trip)
    eta_gen:  float = 0.77 ** 0.5       # generating-leg efficiency (product of the two legs = 0.77)
    var_om_per_mwh: float = 4.0         # variable O&M per MWh generated, $/MWh
    interval_hr: float = 0.5            # half-hourly dispatch

    @property
    def eta_roundtrip(self) -> float:
        return self.eta_pump * self.eta_gen

    @property
    def max_stored_energy(self) -> float:
        # We track stored *potential* energy (the intuitive two-term balance in The dispatch problem below). The reservoir
        # holds up to `dischargeable_mwh` of generatable energy, i.e. dischargeable / eta_gen here.
        return self.dischargeable_mwh / self.eta_gen

snowy = StoragePlant()
hours_full_power = round(snowy.dischargeable_mwh / snowy.power_mw, -1)   # ~159, rounded to the nearest ten
print(f"Snowy 2.0: {snowy.power_mw:,.0f} MW, {snowy.dischargeable_mwh/1e3:.0f} GWh, "
      f"round trip {snowy.eta_roundtrip:.0%}, about {hours_full_power:.0f} hours at full power.")

Snowy 2.0: 2,200 MW, 350 GWh, round trip 77%, about 160 hours at full power.

The dispatch problem

Every scheme chooses the same two quantities for every half hour of the year: how hard to pump, \(u_t\) (megawatts bought), and how hard to generate, \(g_t\) (megawatts sold). The objective is profit. As always,

\[ \text{profit} = \text{revenue} - \text{cost}. \]

Revenue. Only generating earns money. Selling \(g_t\) megawatts for \(\Delta t\) hours at the spot price \(p_t\) brings in \(p_t\, g_t\, \Delta t\), so over the year

\[ \text{revenue} = \sum_t p_t\, g_t\, \Delta t . \]

Cost. Two costs change with the schedule: the electricity bought to pump (\(u_t\)), and variable operations and maintenance on what is generated:

\[ \text{cost} = \sum_t p_t\, u_t\, \Delta t \;+\; c_{\text{om}} \sum_t g_t\, \Delta t . \]

Putting them together, the operator maximises

\[ \text{profit} = \sum_t p_t\,(g_t - u_t)\,\Delta t \;-\; c_{\text{om}} \sum_t g_t\, \Delta t . \]

The constraints. Four constraints complete the problem; the first three are ordinary physics, and the fourth is an accounting convention.

1. Power. Pump or generate at most the rating: \(0 \le u_t \le P\) and \(0 \le g_t \le P\).
2. The store fills and drains. Let \(s_t\) be the energy in the reservoir at the start of interval \(t\). It changes by what flows in minus what flows out: \[ s_{t+1} = s_t + \eta_{\text{pump}}\, u_t\, \Delta t \;-\; \frac{g_t}{\eta_{\text{gen}}}\, \Delta t . \] Pumping adds only the fraction \(\eta_{\text{pump}}\) of the energy bought; generating draws \(1/\eta_{\text{gen}}\) from the store per unit sold. Pump in then generate out and you keep \(\eta_{\text{pump}}\eta_{\text{gen}} = 0.77\), the round trip. (One wording note: \(s_t\) tracks potential energy, the energy sitting in the reservoir; the generatable energy that could be sold from it is smaller by the factor \(\eta_{\text{gen}}\).)
3. Capacity. The store never overflows or goes negative: \(0 \le s_t \le E_{\max}/\eta_{\text{gen}}\). Since \(s_t\) is potential energy, the cap is the level at which the store holds exactly \(E_{\max} = 350\) GWh of generatable energy.
4. A repeating year. Require \(s_1 = s_{T+1}\): end the year with the store as full as it began.

All three schemes in the Analysis are evaluated by the same profit measure and face the same physical and year-end constraints. They differ along two dimensions: the price information available when dispatch is chosen, and the flexibility of the policy that turns that information and the reservoir state into an action (a continuous feasible schedule in Scheme 1, two fixed triggers in Scheme 2, a learned state-dependent action map in Scheme 3).

Assumptions

The assumptions that genuinely move the answer are stated up front, each with the reason we make it, the direction it pushes the estimate, and whether we test it. The bracketed labels key each one to the full audit of all fourteen (A1 to A14) in the appendix, which later sections cite by label.

• Snowy is a price taker (A1). The plant can buy and sell any amount without moving the price: the standard first pass, and it keeps every scheme fast and clean. Over-states profit, probably by a lot; a first-pass reduced-form price-impact band in the appendix sizes the effect.
• Storage is 350 GWh, taken at face value (A4). Snowy’s headline figure, though the cyclable capacity is disputed. Possibly over-states profit; tested in Scheme 1’s storage-capacity sensitivity, and shown non-binding there.
• A frictionless plant (A8, A11, A12). No ramp, minimum-load, or start-up limits; no planned or forced outages; no network limits on clearing at the NSW regional price. Each over-states profit; none is tested here.

Perfect foresight (A2) is deliberately not on this list: it is not an assumption of the analysis but the definition of Scheme 1, and Scheme 3 exists to remove it. Every assumption above pushes the same way, which is why even the rule of thumb here leans in the project’s favour. A1 does most of the work, and we return to it at the end.

Water value

The constraints say what the plant can do. Deciding what it should do comes down to one question: what is a stored megawatt-hour worth? Hydro traders call that number the water value, and it answers every decision at once: sell only when the price beats it, buy only when the price sits far enough below it to survive the round-trip loss, and otherwise wait. The forecaster in Scheme 3 is this calculation industrialised, with the water value shifting with the time of day, the current price, and how full the store is. The folded example below works the number out once, in the simplest possible market.

NoteA two-price example

Suppose one dischargeable megawatt-hour held now must be sold tonight, when the price will be either $50 or $300 per MWh, with equal odds. Its expected sale margin, after the $4/MWh variable operating cost, is \(\tfrac{1}{2}(50 - 4) + \tfrac{1}{2}(300 - 4) = \$171/\text{MWh}\). That expected margin is the water value, here in dischargeable-output units, and it answers both decisions. Sell now only if the current sale margin beats $171/MWh, that is, only if the price beats about $175/MWh, since that is what the energy earns on average by waiting. Pump now only when it pays: buying one megawatt-hour delivers just 0.77 of one for later sale after the round-trip loss, worth 0.77 × $171 ≈ $132, so buy only below about $132/MWh. Between roughly $132 and $175/MWh, do nothing.

(One unit note: the constraints above track potential energy \(s_t\) rather than dischargeable output, so in those units the water value is smaller by the factor \(\eta_{\text{gen}}\); the economic decision is unchanged.)

Analysis

The three dispatch schemes

Every scheme below trades the same plant against the same prices; they differ in what the operator is allowed to know, and in how freely its policy can act on that knowledge, the two dimensions set out at the end of The Model. Those differences are what we are measuring.

	Scheme	The operator knows	Sees the year it trades?	Role here
S1	Perfect foresight	Every future price, exactly	Yes, perfectly	The ceiling: no strategy can beat it
S2	Rule of thumb	Two trigger prices, frozen on 2022 to 2024	No	The deployable baseline: honest to beat
S3	The forecaster	The statistical pattern of prices, never the actual path	No: 2025 fully held out	The central estimate

The “Sees the year it trades?” column is about dispatch decisions: only S1 uses the scored year’s actual prices to choose when to trade, so its profit flatters it. S2 and S3 both fix their decision rule before 2025, S2 as two frozen price triggers and S3 as a learned pattern, and react only to information available at each half-hour. All three do use the realised annual path, ex post, for one accounting step: choosing a start = end reservoir level so no year banks uncredited inventory. That choice can shift the starting stock, and so the realised reservoir path and trades, but it never changes the frozen rules or supplies future prices when a trade is chosen.

Two gaps between the schemes carry the meaning. The gap from S2 up to S3 is the value of trading intelligently on the pattern rather than mechanically. The gap from S3 up to S1 is the realised value of a perfect forecast (the perfect-foresight gap): how much money knowing the future would actually have been worth. If that second gap turns out small, foresight is not the binding constraint on this project, and the bigger worries live elsewhere (Assumption A1 above).

Scheme 1: Perfect foresight

We start from the impossible best case. It takes the dispatch problem above and adds one extraordinary information assumption (A2): the optimiser receives the entire realised price path before choosing the year’s schedule, so it can time every half hour against prices no real operator could know. Every term in the objective and the constraints is a straight line in the decisions, so this is a linear program, which a solver handles exactly in about a second. Within the common constraints, this scheme may also choose any feasible pumping or generation level in every half hour; it is not restricted to fixed triggers or a learned action map, so its advantage over the later schemes is one of policy freedom as well as information.

NoteWhy no “do not pump and generate at once” rule is needed

While the price is positive it never pays to pump and generate in the same interval: cutting both by the same amount leaves revenue unchanged but adds energy to the store (the round-trip loss saved), which only helps later. When the price is negative that argument breaks down: a plant that is paid to consume could in principle pump and generate at once purely to dissipate energy, and the optimiser would take that deal. Rather than complicate the program with an either-or rule (which would destroy its linearity), we verify the schedules it actually produces: in all four solved years the plant never pumps and generates in the same half hour, and the solver code below asserts that on every run.

The solver below is a direct transcription of the dispatch problem in The Model.

def solve_perfect_foresight(prices: np.ndarray, plant: StoragePlant, cyclic: bool = True) -> dict:
    '''Best arbitrage schedule with every price known in advance (the ceiling). Returns total
    profit and the half-hourly pump / generate / stored-energy paths. Mirrors the equations above.'''
    n  = len(prices)
    dt = plant.interval_hr
    P  = plant.power_mw

    # Decisions the solver chooses, one value per half hour:
    pump_power = cp.Variable(n, nonneg=True)      # buy / charge, MW
    gen_power  = cp.Variable(n, nonneg=True)      # sell / discharge, MW
    stored     = cp.Variable(n + 1, nonneg=True)  # stored potential energy at each boundary, MWh

    constraints = [
        pump_power <= P, gen_power <= P,          # (1) power rating
        stored <= plant.max_stored_energy,        # (3) reservoir capacity
        # (2) energy balance, exactly the two-term equation: inflow adds eta_pump of what we buy,
        #     outflow draws 1/eta_gen of the store per unit generated.
        stored[1:] == stored[:-1] + plant.eta_pump * pump_power * dt - (gen_power / plant.eta_gen) * dt,
    ]
    # (4) a repeating year: end as full as it began (or pin both ends empty if cyclic=False)
    constraints += [stored[0] == stored[n]] if cyclic else [stored[0] == 0, stored[n] == 0]

    # profit = sales - purchases - variable O&M on generation
    sales     = cp.sum(cp.multiply(prices, gen_power) * dt)
    purchases = cp.sum(cp.multiply(prices, pump_power) * dt)
    variable_om = plant.var_om_per_mwh * cp.sum(gen_power * dt)
    problem = cp.Problem(cp.Maximize(sales - purchases - variable_om), constraints)
    problem.solve(solver=cp.HIGHS)

    g, u, s = gen_power.value, pump_power.value, stored.value
    # The LP has no explicit rule against pumping and generating at once (see the fold above);
    # assert on every solve that the optimum never actually does it.
    assert float(np.minimum(u, g).max()) < 1e-6, "the LP pumped and generated simultaneously"
    return {"status": problem.status, "profit": float(problem.value),
            "gen": g, "pump": u, "stored": s,
            "net_dollars": prices * (g - u) * dt - plant.var_om_per_mwh * g * dt,
            "generation_twh": g.sum() * dt / 1e6,
            "capacity_factor": g.sum() * dt / (P * len(prices) * dt),
            # volume-weighted average prices; NaN rather than a crash if a leg never runs
            "mean_sell": (g * prices).sum() / g.sum() if g.sum() > 0 else float("nan"),
            "mean_buy":  (u * prices).sum() / u.sum() if u.sum() > 0 else float("nan")}

check = solve_perfect_foresight(price_close.to_numpy(), snowy)
print(f"{CLOSE_YEAR} solver status:", check["status"])   # expect 'optimal'

2025 solver status: optimal

The ceiling, 2022 to 2025

Rather than trust one year, we run the model on four. We report the answer as a band rather than a single number, because a large slice of the top of it depends on a few extreme-price hours that Snowy itself would erase.

Two reading notes for the table. Perfect foresight and Ex-spikes are the pair to compare: the difference between them is how much of the ceiling rides on a handful of extreme hours. The last three columns (the capacity factor and the average buy and sell prices) return at the end of the report as the backtest’s summary statistics.

# Solve every year twice: with all prices, and with the top 1% of spikes trimmed to the
# 99th-percentile price (the part where the price-taker assumption is least believable).
ceiling_by_year, rows = {}, []
for year in YEARS:
    year_prices = prices_by_year[year].to_numpy()
    full = solve_perfect_foresight(year_prices, snowy)
    trimmed = solve_perfect_foresight(np.minimum(year_prices, np.percentile(year_prices, 99)), snowy)
    ceiling_by_year[year] = full                                   # kept for the scheme comparison later
    rows.append({"Year": year, "Perfect foresight ($bn)": full["profit"]/1e9, "Ex-spikes ($bn)": trimmed["profit"]/1e9,
                 "Capacity factor": full["capacity_factor"], "Mean buy ($/MWh)": full["mean_buy"],
                 "Mean sell ($/MWh)": full["mean_sell"]})
panel = pd.DataFrame(rows).set_index("Year")
panel.round(2)

	Perfect foresight ($bn)	Ex-spikes ($bn)	Capacity factor	Mean buy ($/MWh)	Mean sell ($/MWh)
Year
2022	0.89	0.75	0.27	101.39	307.74
2023	0.71	0.57	0.30	38.84	178.80
2024	1.43	0.79	0.34	40.46	276.11
2025	1.13	0.77	0.35	28.72	210.16

low, high = panel["Perfect foresight ($bn)"].min(), panel["Perfect foresight ($bn)"].max()
trim_low, trim_high = panel["Ex-spikes ($bn)"].min(), panel["Ex-spikes ($bn)"].max()
display(Markdown(
    f"**The ceiling: \\${low:.1f} bn to \\${high:.1f} bn of arbitrage profit per year** including "
    f"every extreme spike, or **\\${trim_low:.1f} bn to \\${trim_high:.1f} bn** with prices in the most "
    f"extreme 1% of half-hours capped at the 99th-percentile level. No operating strategy can beat this; "
    f"the next two schemes measure how "
    f"much of it an operator can capture without seeing the future."
))

The ceiling: $0.7 bn to $1.4 bn of arbitrage profit per year including every extreme spike, or $0.6 bn to $0.8 bn with prices in the most extreme 1% of half-hours capped at the 99th-percentile level. No operating strategy can beat this; the next two schemes measure how much of it an operator can capture without seeing the future.

Profit concentration

Sort 2025’s half-hours from most to least profitable and add them up; the curve below carries that running total across the whole year, as a share of the year’s final profit. It tells the story in three phases. First a near-vertical climb: the best ~330 hours (about 4% of the year) earn half the profit, and those are exactly the extreme-price hours where the price-taker assumption breaks down worst. The climb continues until every money-making hour has been counted, at about 126% of the final profit. The curve then runs flat across the hours the plant sat idle, and finally descends back to 100%: the give-back on the right is the pumping bill, the cheap hours in which the plant pays to charge. Selling high only nets the year’s profit after buying cheap has been paid for.

# Rank the close-up year's half-hours by the profit each contributed, then accumulate: the
# curve's height at x is the share of the year's final profit earned by its x most profitable
# hours. Generating hours push it up, idle hours leave it flat, and pumping hours (which cost
# money) pull it back down to exactly 100% at the year's end.
result_close = ceiling_by_year[CLOSE_YEAR]
profit_sorted = np.sort(result_close["net_dollars"])[::-1]
total_year_hours = len(profit_sorted) * 0.5                     # leap-safe length of the year in hours
cumulative_share = np.cumsum(profit_sorted) / profit_sorted.sum() * 100
hours_ranked = np.arange(1, len(profit_sorted)+1) * 0.5
hours_for_half = hours_ranked[np.searchsorted(cumulative_share, 50)]   # hours that earn half the profit
peak_ix = cumulative_share.argmax()                             # the last hour that still adds money

fig, ax = plt.subplots()
ax.plot(hours_ranked, cumulative_share, color=INK, lw=2.5)                       # bold black data line
ax.axhline(100, color=MUTED, lw=1, ls=":")                      # where the year must end
ax.axhline(50, color=MUTED, lw=1, ls=":")
ax.axvline(hours_for_half, color=SPIKE, lw=1, ls=":")
ax.fill_between(hours_ranked, cumulative_share, where=(hours_ranked <= hours_for_half), color=SPIKE, alpha=0.2)
ax.annotate(f"The best ~{hours_for_half:,.0f} hours\n(~{hours_for_half/total_year_hours*100:.0f}% of the year)\nearn half the profit",
            xy=(hours_for_half, 50), xytext=(hours_for_half + 380, 22), color=SPIKE, fontsize=10,
            arrowprops=dict(arrowstyle="->", color=SPIKE))
ax.annotate(f"every money-making hour counted: {cumulative_share[peak_ix]:.0f}%",
            xy=(hours_ranked[peak_ix], cumulative_share[peak_ix]), xytext=(2500, 135),
            color=INK2, fontsize=9, arrowprops=dict(arrowstyle="->", color=MUTED))
ax.annotate("the descent is the pumping bill:\nhours the plant pays to charge",
            xy=(7900, float(cumulative_share[int(7900 / 0.5) - 1])), xytext=(5300, 60),
            color=BUY, fontsize=9, arrowprops=dict(arrowstyle="->", color=BUY))
ax.set(title=f"Profit Distribution, Perfect Foresight ({CLOSE_YEAR})",
       ylabel="Cumulative share of the year's profit (%)",
       xlim=(0, total_year_hours), ylim=(0, 142))
# the x label as one line with reading arrows at each end, in place of a plain set_xlabel
ax.text(0.5, -0.13, "Hours of the year, ranked by profit", transform=ax.transAxes, ha="center", color=INK)
ax.text(0, -0.13, "← most profitable", transform=ax.transAxes, ha="left", fontsize=9, color=INK2)
ax.text(1, -0.13, "least profitable →", transform=ax.transAxes, ha="right", fontsize=9, color=INK2)
finish(fig, "06-profit-concentration.png")

Sensitivity to storage capacity

Snowy’s headline 350 GWh is contested: the cyclable capacity may be far smaller (see Background and Data). Re-running 2025 with much smaller stores barely changes the result, because the value is in daily cycling, not in holding energy for weeks. The headline figure can therefore be taken at face value without the answer leaning on it.

# Re-solve the close-up year's ceiling with progressively smaller reservoirs and report each
# as a share of the 350 GWh answer: if the share stays high, the contested size does not matter.
full_ceiling = ceiling_by_year[CLOSE_YEAR]["profit"]
storage_rows = []
for gwh in (50, 175, 350):
    smaller = solve_perfect_foresight(price_close.to_numpy(), replace(snowy, dischargeable_mwh=gwh*1000))
    storage_rows.append({"Storage (GWh)": gwh, "Perfect foresight ($bn)": round(smaller["profit"]/1e9, 2),
                         "Share of 350 GWh result": f"{round(smaller['profit']/full_ceiling*100)}%"})
pd.DataFrame(storage_rows).set_index("Storage (GWh)")

	Perfect foresight ($bn)	Share of 350 GWh result
Storage (GWh)
50	1.08	95%
175	1.12	99%
350	1.13	100%

The optimal schedule over one week

Over one real week, everything the optimal plan does sits in one panel: the price (black line, left axis), the store level (green line, right axis), and the plant’s actions as shaded columns (blue while pumping, orange while generating). The model fills the store when prices dip and empties it into the peaks, and the green line ties each action to the next: every climb happens inside blue columns, every drain inside orange ones.

# One winter week of the optimal plan, merged into a single panel: price on the left axis, the
# reservoir level on the right, and the actions as light full-height columns (blue = pumping,
# orange = generating) so the timing reads directly against the price line.
week = slice("2025-06-30", "2025-07-06")
week_index = price_close.loc[week].index
pos = price_close.index.get_indexer(week_index)
store_gwh = result_close["stored"][pos] * snowy.eta_gen / 1e3       # show the store in generatable GWh (Snowy's units)

fig, price_ax = plt.subplots(figsize=(9, 4.8))
price_ax.fill_between(week_index, 0, 1, where=result_close["pump"][pos] > 0, step="mid",
                      transform=price_ax.get_xaxis_transform(), color=BUY, alpha=0.15)
price_ax.fill_between(week_index, 0, 1, where=result_close["gen"][pos] > 0, step="mid",
                      transform=price_ax.get_xaxis_transform(), color=SELL, alpha=0.15)
price_ax.plot(week_index, price_close.loc[week].values, color=INK, lw=1.6, zorder=3)
price_ax.set(title="Optimal Strategy, Perfect Foresight (30 June to 6 July 2025)", ylabel="Price ($/MWh)")
price_ax.yaxis.set_major_formatter(dollars)
store_ax = price_ax.twinx()                                         # state of charge as a bold green line
store_ax.plot(week_index, store_gwh, color=STORE, lw=2.4)
store_ax.set_ylabel("State of charge (GWh)", color=STORE)
store_ax.set_ylim(bottom=0)
store_ax.tick_params(axis="y", colors=STORE)
store_ax.spines["right"].set_color(STORE)
# in-figure legend, because this PNG is reused outside the report without its caption; placed
# under the plot so it cannot collide with the title
from matplotlib.patches import Patch
price_ax.legend(handles=[plt.Line2D([], [], color=INK, lw=1.6, label="Price"),
                         plt.Line2D([], [], color=STORE, lw=2.4, label="State of charge"),
                         Patch(color=BUY, alpha=0.3, label="Pumping"),
                         Patch(color=SELL, alpha=0.3, label="Generating")],
                frameon=False, loc="upper center", ncol=4, fontsize=8.5, bbox_to_anchor=(0.5, -0.22))
fig.autofmt_xdate()
finish(fig, "07-optimal-week.png")

Scheme 2: The rule of thumb

The opposite extreme is an operator with two fixed trigger prices and no discretion. Buy at full power whenever the price sits at or below the cheaper trigger, sell whenever it sits at or above the more expensive one, do nothing in between. The triggers are set the way a real operator would have to set them: chosen once on the training years (2022 to 2024), frozen, and then applied unchanged. So this is not a rule tuned with hindsight to the year it trades; it commits to its prices in advance and is judged out of sample on the held-out 2025, exactly as the forecaster is. It reacts only to the current price against its two fixed levels: it never learns the time-of-day pattern and never forecasts the coming sequence. Its role is to show what plain mechanical trading at sensible price levels would have earned, so the forecaster that follows has an honest, deployable baseline to beat.

# ---- Part 1. Simulate one year of the two-trigger rule ----
def threshold_dispatch(prices, plant, buy_below, sell_above, start_soc=0.0):
    '''A fixed-threshold rule with explicit DOLLAR triggers, run from a given start level: buy at
    full power whenever the price is at or below `buy_below`, sell whenever it is at or above
    `sell_above`, idle in between. The triggers are fixed dollars, so they can be frozen on the
    training years rather than read off the year being scored, and returning the reservoir path
    (with a `start_soc`) lets us impose the same start = end boundary the other two schemes use.'''
    dt, P = plant.interval_hr, plant.power_mw
    max_store, eta_pump, eta_gen = plant.max_stored_energy, plant.eta_pump, plant.eta_gen
    stored, profit = float(start_soc), 0.0
    path = np.empty(len(prices) + 1)                               # reservoir level at each boundary
    path[0] = stored
    for i, price in enumerate(prices):
        if price <= buy_below and stored < max_store:              # at or below buy trigger, room left: buy
            mw = min(P, (max_store - stored) / (eta_pump * dt))    # full power, or just enough to fill up
            stored += eta_pump * mw * dt                           # only eta_pump of the buy arrives
            profit -= price * mw * dt                              # pay for the electricity bought
        elif price >= sell_above and stored > 0:                   # at or above sell trigger, energy held: sell
            mw = min(P, stored * eta_gen / dt)                     # full power, or just what is left
            stored -= (mw / eta_gen) * dt                          # each MWh sold draws 1/eta_gen from store
            profit += price * mw * dt - plant.var_om_per_mwh * mw * dt   # earn the price, net of O&M
        path[i + 1] = stored
    return {"profit": profit, "stored": path}


# ---- Part 2. The start = end boundary search, shared by Schemes 2 and 3 ----
def find_cyclic_boundary(run_year, plant, tol_dis_mwh=100.0,
                         start_fracs=(0.0, 0.25, 0.5, 0.75, 1.0), max_iters=40):
    '''Impose the REPEATED-YEAR (start = end reservoir) boundary on any dispatch rule, so every
    scheme is compared on identical accounting. `run_year(start_soc)` simulates one year from
    that starting level and returns at least the reservoir path ("stored") and the year's
    "profit". From several starting levels, iterate start = end until the year's end level
    matches its start within `tol_dis_mwh` dischargeable MWh, falling back to the
    smallest-residual run if none settles. A dispatch rule can in principle have several
    self-consistent start = end levels, so a single small residual proves a fixed point exists
    but not that the cyclic profit is unique; we therefore report whether every starting level
    settled to the same economic outcome (`unique_fixed_point`), so a non-unique fixed point
    cannot pass silently.'''
    eg = plant.eta_gen
    runs = []                                                       # one record per starting level
    for frac in start_fracs:
        start = frac * plant.max_stored_energy
        out = run_year(start)
        end = float(out["stored"][-1])
        residual_dis = abs(end - start) * eg                        # start-vs-end gap, dischargeable MWh
        it = 1
        while residual_dis >= tol_dis_mwh and it < max_iters:
            start = end                                             # feed the end level back in as the start
            out = run_year(start)
            end = float(out["stored"][-1])
            residual_dis = abs(end - start) * eg
            it += 1
        runs.append({"start": start, "end": end, "outcome": out,
                     "residual_dis_mwh": residual_dis, "iterations": it})
    best = min(runs, key=lambda r: r["residual_dis_mwh"])           # smallest-residual fallback if none converged
    # A unique fixed point means every starting level settled to the same shallow basin: all the
    # settled levels within a few convergence tolerances (5 * tol_dis_mwh) of each other AND
    # their profits agreeing to 0.1%.
    basin_width = (max(r["start"] for r in runs) - min(r["start"] for r in runs)) * eg
    profits = [r["outcome"]["profit"] for r in runs]
    unique = (basin_width < 5 * tol_dis_mwh
              and (max(profits) - min(profits)) < 0.001 * abs(best["outcome"]["profit"]))
    return {**best["outcome"], "start_soc": best["start"], "end_soc": best["end"],
            "residual_dis_mwh": best["residual_dis_mwh"], "iterations": best["iterations"],
            "converged": best["residual_dis_mwh"] < tol_dis_mwh,
            "unique_fixed_point": unique,
            "basin_dis_mwh": basin_width,                           # spread of the settled levels
            "profit_spread": (max(profits) - min(profits)) / abs(best["outcome"]["profit"]),
            "all_starts": [(r["start"], r["end"], r["outcome"]["profit"]) for r in runs]}


def cyclic_threshold_score(prices, plant, buy_below, sell_above):
    '''Score a dollar-threshold rule under the repeated-year convention: run the fixed triggers
    through the start = end search above, with all of its convergence and uniqueness checks.'''
    return find_cyclic_boundary(
        lambda start: threshold_dispatch(prices, plant, buy_below, sell_above, start_soc=start),
        plant)

# ---- Part 3. Choose the triggers on 2022 to 2024, freeze them, score every year ----
# Train the rule like a deployable strategy, not one fitted to the year it is scored on. Its two
# dollar triggers are chosen on the training years alone (2022 to 2024): take a small grid of buy
# and sell percentiles of the pooled training prices and convert each to a fixed dollar level.
TRAIN_YEARS = [2022, 2023, 2024]                                  # 2025 stays unseen until scoring
train_prices = np.concatenate([prices_by_year[y].to_numpy() for y in TRAIN_YEARS])
candidates = [(np.percentile(train_prices, low), np.percentile(train_prices, high))
              for low in (10, 20, 30) for high in (70, 80, 90)]

# Score every candidate on every training year up front, and assert each cyclic solution is a
# converged, UNIQUE fixed point (a bang-bang rule can in principle settle at more than one
# self-consistent level). The training objective is then built only from validated results, so an
# ill-defined candidate cannot silently enter the maximisation.
candidate_scores = {pair: {y: cyclic_threshold_score(prices_by_year[y].to_numpy(), snowy, *pair)
                           for y in TRAIN_YEARS} for pair in candidates}
assert all(r["converged"] and r["unique_fixed_point"]
           for by_year in candidate_scores.values() for r in by_year.values()), \
    "a candidate rule's cyclic boundary is not a unique fixed point on some training year"
# Keep the pair with the highest total profit summed across the three training years.
buy_trigger, sell_trigger = max(
    candidates, key=lambda pair: sum(candidate_scores[pair][y]["profit"] for y in TRAIN_YEARS))

# Score the frozen triggers on every year under the same start = end boundary (2025 out of sample),
# again asserting a converged, unique fixed point so each reported annual profit is well defined.
s2_scored = {year: cyclic_threshold_score(prices_by_year[year].to_numpy(), snowy,
                                          buy_trigger, sell_trigger) for year in YEARS}
assert all(s2_scored[y]["converged"] and s2_scored[y]["unique_fixed_point"] for y in YEARS), \
    "the frozen rule's cyclic boundary is not a unique fixed point in some year"
s2_by_year = {year: s2_scored[year]["profit"] for year in YEARS}
assert all(s2_by_year[y] <= ceiling_by_year[y]["profit"] + 1e-6 for y in YEARS), \
    "a rule of thumb should never beat perfect foresight"
held_out_s2 = s2_by_year[2025] / ceiling_by_year[2025]["profit"] * 100
display(Markdown(
    f"Frozen on 2022 to 2024, the rule buys at or below **\\${buy_trigger:.0f}/MWh** and sells at "
    f"or above **\\${sell_trigger:.0f}/MWh**. Applied unchanged to the held-out 2025, it earns "
    f"**\\${s2_by_year[2025]/1e9:.2f} bn**, or **{held_out_s2:.0f}%** of that year's ceiling. The "
    f"gap between this mechanical rule and the ceiling is the room in which forecasting and "
    f"smarter policies operate; Scheme 3 shows how much of it the forecaster captures."
))

Frozen on 2022 to 2024, the rule buys at or below $68/MWh and sells at or above $133/MWh. Applied unchanged to the held-out 2025, it earns $0.90 bn, or 80% of that year’s ceiling. The gap between this mechanical rule and the ceiling is the room in which forecasting and smarter policies operate; Scheme 3 shows how much of it the forecaster captures.

Scheme 3: The forecaster

No one knows tomorrow’s prices. But a trading desk is not guessing blindly either: electricity prices have a strong, stable daily rhythm. The figure below draws the operator’s raw material: every one of the 1,096 days the forecaster is allowed to learn from (2022 to 2024), each as one faint line, with the median training day in bold. Individual days are wild, and some blow straight through the top of the chart, yet the shape barely moves: cheap in the middle of the day, expensive after sunset. The bold line is what a disciplined price expectation looks like: not a guess about any particular tomorrow, but the pattern three years of days keep returning to. That repeating shape, together with the odds of moving from one price level to another as the day unfolds, is knowledge an operator genuinely has.

# The forecaster's raw material: all training days (2022 to 2024) as faint lines, one per
# day, with the median across days at each half-hour slot drawn in bold. The y axis is
# deliberately clipped: spike days genuinely run off the top of this chart (the worst reaches
# about $17,500/MWh), which is the point: levels are wild, the shape is stable.
train_days = pd.concat([prices_by_year[y] for y in TRAIN_YEARS])
day_matrix = (pd.DataFrame({"date": train_days.index.date,
                            "tod": train_days.index.hour + train_days.index.minute / 60,
                            "price": train_days.to_numpy()})
              .pivot(index="date", columns="tod", values="price"))
x = day_matrix.columns.to_numpy()
median_day = day_matrix.median(axis=0).to_numpy()             # the pattern the operator distils

fig, ax = plt.subplots()
ax.plot(x, day_matrix.T.to_numpy(), color=MUTED, alpha=0.06, lw=0.7)
ax.plot(x, median_day, color=INK, lw=2.8)
ymax = 800.0
spike_days = int((day_matrix.max(axis=1) > ymax).sum())       # days whose peak exceeds the frame
ax.set_ylim(-130, ymax)
ax.set(title="Daily Prices (2022 to 2024)", xlabel="Time of day", ylabel="Price ($/MWh)")
clock_axis(ax)
ax.yaxis.set_major_formatter(dollars)
ax.annotate(f"{spike_days} spike days run off the top of this chart", (0.3, ymax * 0.94),
            color=INK2, fontsize=9)
ax.legend(handles=[plt.Line2D([], [], color=INK, lw=2.8, label="Median price")],
          frameon=False, loc="upper left", bbox_to_anchor=(0.01, 0.90))
finish(fig, "08-daily-rhythm.png")

Summarised into bands, the same 1,096 days become the expectation the forecaster disciplines itself with. The wild individual days of the previous figure smooth into a stable daily shape: half of all days stay inside the narrow inner band, eight in ten inside the outer one, and the bands themselves tell the operator when the cheap and expensive parts of a typical day arrive.

# The same day-by-slot matrix as above, summarised: at each half-hour of the day, the median
# across the 1,096 training days with the middle 50% and middle 80% of days as bands. This is
# the smoothed view of the spaghetti above: the day-shape the forecaster's beliefs encode.
q10, q25, q75, q90 = (day_matrix.quantile(q, axis=0).to_numpy() for q in (0.10, 0.25, 0.75, 0.90))

fig, ax = plt.subplots()
ax.fill_between(x, q10, q90, color=MUTED, alpha=0.25, label="Middle 80% of days")
ax.fill_between(x, q25, q75, color=MUTED, alpha=0.5, label="Middle 50% of days")
ax.plot(x, median_day, color=INK, lw=2.8, label="Median price")
ax.set(title="Daily Prices, Summarised (2022 to 2024)", xlabel="Time of day", ylabel="Price ($/MWh)")
clock_axis(ax)
ax.yaxis.set_major_formatter(dollars)
ax.legend(frameon=False, loc="upper left")
finish(fig, "08a-daily-rhythm-bands.png")

We give the operator exactly the knowledge in that water-value example, learned from real data rather than assumed. From that cloud of days to a decision, the forecaster is built in four steps:

1. Compress prices into 16 buckets, from “strongly negative” up to “extreme spike”, so that “where is the price right now?” always has one of sixteen answers.
2. Estimate the odds of the next move: for every half-hour of the day, a 16-by-16 table of probabilities for where the price sits next, given where it sits now (the transition matrix, drawn below once fitted).
3. Attach a dollar value to stored energy in every situation (time of day, price bucket, store level) by backward induction: the water value.
4. Act by comparison: generate when the price beats the water value, pump when it sits below the water value less the round-trip loss, and otherwise wait.

Steps 1 and 2 are the operator’s beliefs; steps 3 and 4 turn those beliefs into a rule. The paragraphs below expand each in turn.

What the operator believes. Prices are grouped into the 16 buckets above; from history we estimate, for every half hour of the day, the odds of moving from each bucket to each other bucket in the next half hour. This is a Markov model of the price: it knows the daily rhythm and how sticky prices are, but it never knows what will actually happen.

NoteHow the beliefs are estimated, step by step

1. Pool the training prices (2022 to 2024: 52,608 half-hours) and cut them into 16 buckets with equal head-count, so each bucket holds exactly 6.25% of the observations. Equal head-count rather than equal dollar-width spends the model’s resolution where the data actually lives: the everyday $38 to $150 range gets many fine buckets of its own, while a single wide bucket absorbs everything above $300.
2. Count the moves. For each of the 48 half-hour slots of the day and each pair of buckets (i, j), count how often the price sat in bucket i at that slot and in bucket j one half-hour later. Only genuinely consecutive half-hours count (a pair straddling a data gap does not), and no move is counted across a year boundary.
3. Keep rare moves possible. Add one phantom observation to every cell (Laplace smoothing), so a jump never seen in three years keeps a small positive probability rather than being declared impossible. Three years is long enough to learn the rhythm, not long enough to have seen every extreme.
4. Turn counts into odds. Divide each row by its total. The result is 48 tables of 16 × 16 probabilities, and that is the operator’s entire view of the future: no gas prices, no weather forecasts, no order book. Just: given where the price is now and the time of day, where does it tend to go next?

What the operator decides. At every half hour the operator sees the time of day, the current price bucket, and how full the store is, and picks one of three actions: pump at full power, do nothing, or generate at full power (Assumption A13).

How the best decisions are found. The plant is a going concern: there is no last day, and the best rule for a Tuesday in March is the best rule for any day with the same price and store level. So we solve for a stationary rule: one map from time of day, price bucket and store level to an action, used every day, forever. The method is backward induction applied to a single day, repeatedly: sweep backwards through the day’s 48 half-hours to value every situation, feed that valuation back in as the value of reaching tomorrow, and sweep again, until the shape of the valuation stops changing. (Its level keeps growing by one day’s average profit each sweep, which is why only the shape can settle; a reference constant is subtracted each sweep to handle this, a standard method called relative value iteration.) Every earlier decision is then a comparison of “cash now” against the expected value of the energy later. The result attaches a dollar value to every stored megawatt-hour in every situation, which hydro traders call the water value. The optimal behaviour follows automatically: generate when the price is above the water value, pump when it is below the water value minus the round-trip loss, and otherwise wait.

NoteThe Bellman equation

Let \(V_t(s, k)\) be the expected future profit at slot \(t\) of the day with store level \(s\) and price bucket \(k\). Backward induction repeatedly evaluates \[ V_t(s,k) \;=\; \max_{a \,\in\, \{\text{pump},\, \text{idle},\, \text{generate}\}} \Big\{ r(a, s, k) \;+\; \mathbb{E}\big[\, V_{t+1}(s', k') \,\big|\, k \,\big] \Big\}, \] where \(r\) is the immediate cash flow of the action, \(s'\) is the store level it leads to, and the expectation runs over the estimated odds of the next price bucket \(k'\). The horizon is infinite: the plant has no final day, so the equation is solved to a fixed point rather than swept over a fixed number of days. The one-day update above is applied repeatedly, a reference constant is subtracted after each sweep (relative value iteration), and the iteration stops when the one-sweep change in \(V\) is flat across all states to within $100. The retained rule is then checked directly: solving on past the stopping point no longer moves its scored profit (the convergence evidence under Performance). The implementation is written out in full just below, folded by default.

Specification choices

Choice	Value	Why this value, and the tradeoff
Price buckets	16, equal head-count	Enough to separate “negative” from “cheap” from “spike”; more buckets thin out the history each estimate rests on.
Store grid	61 levels	The value of stored energy is smooth, so a modest grid with interpolation is accurate and fast.
Horizon	Infinite (stationary rule)	The plant has no final day, so the one-day update is repeated to a fixed point: a $100 span tolerance on the value shape, with the scored profit’s stability shown under Performance.
Annual boundary	Repeated-year (cyclic)	Each scored year starts at the store level it ends at, mirroring Scheme 1’s cyclic year; the level is found ex post, as disclosed under Performance.
Training years	2022, 2023, 2024	The pattern the operator is allowed to learn from (A14).
Test year	2025, held out	The honest check: a year the model has never seen. (The dispatch rule is held out; the annual boundary normalisation is ex post, as disclosed under Performance.)

The training years include the 2022 fuel crisis and two calmer years, so the learned pattern is neither all-crisis nor all-calm.

The functions that do this work are the heart of the report’s central estimate, so they are written out in full below (folded by default): a small PriceModel type that holds the operator’s beliefs, then four functions. In order: fit_price_model learns those beliefs, solve_known_distribution turns them into the stationary decision rule by relative value iteration, run_policy scores the rule against real prices, and cyclic_score applies the repeated-year boundary convention explained under Performance. The third is the honesty guarantee of the whole exercise: the buckets only ever choose the action; every dollar is paid and earned at the realised market price.

# ---- Part 1. What the operator believes: the price model and its fitting ----
@dataclass(frozen=True)
class PriceModel:
    '''A time-of-day Markov model of the spot price: what the forecaster believes about
    tomorrow. Prices are grouped into `n_bins` quantile bins; `transitions[t, i, j]` is the
    estimated probability that the price moves from bin i at half-hour slot t of the day to
    bin j at slot t+1. The operator knows this DISTRIBUTION, never the realised path.'''
    edges: np.ndarray        # (K+1,) bin edges, open at both ends so every price lands somewhere
    rep_price: np.ndarray    # (K,) representative price of each bin (its training-data mean)
    transitions: np.ndarray  # (slots, K, K) row-stochastic transition matrices, one per slot

    @property
    def n_bins(self) -> int:
        return len(self.rep_price)

    @property
    def slots(self) -> int:
        return self.transitions.shape[0]

    def to_bin(self, x: np.ndarray) -> np.ndarray:
        '''Map prices to bin indices 0..K-1.'''
        return np.clip(np.digitize(x, self.edges[1:-1]), 0, self.n_bins - 1)


def fit_price_model(yearly_series: list, n_bins: int = 16, slots: int = 48) -> PriceModel:
    '''Fit the time-of-day Markov model to training data. `yearly_series` is a list of price
    series, one per training year, so no transition is ever counted across a year boundary.
    Bin edges are quantiles of the pooled data (equal head-count bins); transition counts get
    add-one (Laplace) smoothing so unseen moves keep a small positive probability.

    Each series may be a pandas Series with a DatetimeIndex (preferred: the half-hour-of-day
    slot is read off each timestamp, so a missing interval cannot silently shift later
    observations into the wrong slot) or a plain 1-D array (positional slot numbering,
    valid only for gap-free data). NaN prices are rejected outright rather than silently
    binned.'''
    # Prepare each series: values, each observation's slot of day, and which adjacent pairs
    # are genuinely consecutive half-hours (a pair straddling a data gap is not a transition).
    prepared = []
    for s in yearly_series:
        if isinstance(s, pd.Series) and isinstance(s.index, pd.DatetimeIndex):
            values = s.to_numpy(dtype=float)
            slot = (s.index.hour * 2 + s.index.minute // 30).to_numpy()
            consecutive = np.diff(s.index.to_numpy()) == np.timedelta64(30, "m")
        else:
            values = np.asarray(s, dtype=float)
            slot = np.arange(len(values)) % slots           # positional fallback: assumes no gaps
            consecutive = np.ones(max(len(values) - 1, 0), dtype=bool)
        if np.isnan(values).any():
            raise ValueError("price series contains NaNs; fill or drop them before fitting")
        prepared.append((values, slot, consecutive))

    pooled = np.concatenate([values for values, _, _ in prepared])
    edges = np.quantile(pooled, np.linspace(0, 1, n_bins + 1))
    edges[0], edges[-1] = -np.inf, np.inf
    pooled_bins = np.clip(np.digitize(pooled, edges[1:-1]), 0, n_bins - 1)
    rep_price = np.array([pooled[pooled_bins == k].mean() for k in range(n_bins)])

    counts = np.ones((slots, n_bins, n_bins))               # the Laplace prior
    for values, slot, consecutive in prepared:
        b = np.clip(np.digitize(values, edges[1:-1]), 0, n_bins - 1)
        # tally observed slot-t -> slot-t+1 moves, only across truly consecutive half-hours
        np.add.at(counts, (slot[:-1][consecutive], b[:-1][consecutive], b[1:][consecutive]), 1)
    transitions = counts / counts.sum(axis=2, keepdims=True)
    return PriceModel(edges=edges, rep_price=rep_price, transitions=transitions)


# ---- Part 2. Turn beliefs into a decision rule: relative value iteration ----
def solve_known_distribution(model: PriceModel, plant: StoragePlant, n_soc: int = 61,
                             span_tol: float = 100.0, max_sweeps: int = 2000,
                             checkpoint_sweeps: tuple = ()) -> dict:
    '''The optimal STATIONARY dispatch policy when the operator knows the price distribution
    (the fitted Markov model) but not the realisation. The plant is a going concern with no
    final day, so the right formulation is infinite-horizon: the policy that maximises average
    profit per day, the same rule for every day of every year.

    Solved by relative value iteration on the one-day operator: one sweep is a backward
    induction through the `slots` half-hours of a day, using the current start-of-day value
    V as the continuation; after each sweep a reference constant (V at the first grid point)
    is subtracted so V stays bounded while its SHAPE (the only part that drives decisions)
    converges. The stopping rule is the span (max minus min) of the change in V across one
    sweep, in dollars: the span bounds how far the retained policy's average daily profit can
    sit below the optimum, and it is immune to the constant drift that plain value iteration
    would show. Returns action[slot, soc_index, price_bin] with 0 = idle, 1 = pump at full
    power, 2 = generate at full power, plus the convergence record; `checkpoint_sweeps` asks
    for snapshots of the policy at given sweep counts (for a published convergence table).'''
    dt, P = plant.interval_hr, plant.power_mw
    emax, ep, eg, vom = plant.max_stored_energy, plant.eta_pump, plant.eta_gen, plant.var_om_per_mwh
    K, slots = model.n_bins, model.slots
    soc_grid = np.linspace(0.0, emax, n_soc)
    rep = model.rep_price

    # What each action does, per state-of-charge grid point (identical at every slot):
    pump_mw = np.minimum(P, (emax - soc_grid) / (ep * dt))          # pumping, clipped to the room left
    gen_mw  = np.minimum(P, soc_grid * eg / dt)                     # generating, clipped to the energy held
    next_soc = [soc_grid,                                           # idle: the store is unchanged
                soc_grid + ep * pump_mw * dt,                       # pump: only eta_pump of the buy arrives
                soc_grid - (gen_mw / eg) * dt]                      # generate: draws 1/eta_gen per unit sold
    reward = [np.zeros((n_soc, K)),                                 # idle earns nothing
              -np.outer(pump_mw * dt, rep),                         # pumping buys at the bin price
              np.outer(gen_mw * dt, rep) - (vom * gen_mw * dt)[:, None]]  # selling, net of variable O&M

    V = np.zeros((n_soc, K))                                        # start-of-day relative value
    policy = np.zeros((slots, n_soc, K), dtype=np.int8)
    span_history, checkpoints = [], {}
    sweep, span = 0, float("inf")                                   # in case max_sweeps is 0
    for sweep in range(1, max_sweeps + 1):
        V_prev, W = V, V                                            # W walks backward through the day
        for t in reversed(range(slots)):
            # Value of each action = its immediate reward + the expected value of where it
            # leaves us: interpolate W at the landing state of charge, then average over where
            # the price bin can go next (the slot-t transition row).
            action_value = np.empty((3, n_soc, K))
            for a in range(3):
                landed = np.empty((n_soc, K))
                for j in range(K):
                    landed[:, j] = np.interp(next_soc[a], soc_grid, W[:, j])
                action_value[a] = reward[a] + landed @ model.transitions[t].T
            policy[t] = action_value.argmax(axis=0)                 # this sweep's rule for slot t
            W = action_value.max(axis=0)
        # Span of the one-sweep change: constant growth (the average daily profit) cancels out,
        # so the span measures only how much the SHAPE of V is still moving.
        diff = W - V_prev
        span = float(diff.max() - diff.min())
        span_history.append(span)
        V = W - W[0, 0]                                             # subtract the reference constant
        if sweep in checkpoint_sweeps:
            checkpoints[sweep] = policy.copy()
        if span < span_tol:
            break
    return {"policy": policy, "soc_grid": soc_grid, "sweeps": sweep,
            "span_history": np.array(span_history), "converged": span < span_tol,
            "checkpoints": checkpoints}


# ---- Part 3. Score the rule against real prices, with the start = end boundary ----
def run_policy(prices: np.ndarray, plant: StoragePlant, policy: np.ndarray,
               soc_grid: np.ndarray, model: PriceModel, start_soc: float = 0.0) -> dict:
    '''Score a (slot, soc, bin) policy against a REAL price path. The action is looked up from
    the policy, but every dollar uses the realised price, never the bin's representative price.
    `start_soc` is the reservoir's starting potential energy in MWh (default empty); see
    `cyclic_score` for the repeated-year convention that chooses it.'''
    dt, P = plant.interval_hr, plant.power_mw
    emax, ep, eg, vom = plant.max_stored_energy, plant.eta_pump, plant.eta_gen, plant.var_om_per_mwh
    slots = model.slots
    prices = np.asarray(prices, dtype=float)
    bins = model.to_bin(prices)
    n = len(prices)
    pump, gen, stored = np.zeros(n), np.zeros(n), np.zeros(n + 1)
    soc, profit = float(start_soc), 0.0
    stored[0] = soc
    for i in range(n):
        # look up the action for (this half-hour of the day, nearest store level, price bin)
        action = policy[i % slots, np.abs(soc_grid - soc).argmin(), bins[i]]
        if action == 1 and soc < emax:                              # pump: buy what fits
            mw = min(P, (emax - soc) / (ep * dt))                   # full power, or just enough to fill up
            soc += ep * mw * dt                                     # only eta_pump of the buy arrives
            profit -= prices[i] * mw * dt                           # pay the REAL price, not the bin price
            pump[i] = mw
        elif action == 2 and soc > 0:                               # generate: sell what we hold
            mw = min(P, soc * eg / dt)                              # full power, or just what is left
            soc -= (mw / eg) * dt                                   # each MWh sold draws 1/eta_gen from store
            profit += prices[i] * mw * dt - vom * mw * dt           # earn the REAL price, net of O&M
            gen[i] = mw
        stored[i + 1] = soc
    return {"profit": profit, "pump": pump, "gen": gen, "stored": stored,
            "generation_twh": gen.sum() * dt / 1e6,
            "capacity_factor": gen.sum() * dt / (P * n * dt)}


def cyclic_score(prices: np.ndarray, plant: StoragePlant, policy: np.ndarray,
                 soc_grid: np.ndarray, model: PriceModel) -> dict:
    '''Score a (slot, soc, bin) policy under the repeated-year accounting convention: choose a
    starting reservoir level whose end-of-year level equals it, so no un-cashed inventory
    build-up or draw-down distorts the annual profit. This mirrors the cyclic boundary the
    perfect-foresight LP already uses. The search itself is `find_cyclic_boundary`, defined
    with Scheme 2, so both trained schemes share exactly the same boundary machinery and its
    convergence and uniqueness checks.

    One honesty note, mirrored in the report: the fixed point belongs to the policy COMBINED
    WITH the realised annual price path, not to the policy alone. The whole year is used, ex
    post, to normalise the boundary, though every individual dispatch decision stays a causal
    function of the current time, price bucket and store level.'''
    return find_cyclic_boundary(
        lambda start: run_policy(prices, plant, policy, soc_grid, model, start_soc=start),
        plant)

import time

start = time.time()
# Fit on the full pandas Series (not bare arrays) so each observation's time-of-day slot comes
# from its own timestamp; a gap in the data can then never shift observations into wrong slots.
price_model = fit_price_model([prices_by_year[y] for y in TRAIN_YEARS])
# Solve for the stationary decision rule by relative value iteration, keeping snapshots of
# the rule part-way through the solve for the convergence evidence under Performance.
learned = solve_known_distribution(price_model, snowy, checkpoint_sweeps=(14, 28, 56, 112))
assert learned["converged"], "the value iteration hit its sweep cap before meeting tolerance"
elapsed = time.time() - start
took = f"{elapsed:.1f} seconds" if elapsed >= 1 else "under a second"
print(f"Fitted the price pattern and solved for the stationary decision rule in {took}: "
      f"{learned['sweeps']} one-day sweeps to reach the $100 span tolerance.")

Fitted the price pattern and solved for the stationary decision rule in 1.0 seconds: 178 one-day sweeps to reach the $100 span tolerance.

The beliefs, fitted

Before the rule, the beliefs themselves. First, the buckets: sixteen price levels, cut so that each holds an equal 6.25% share of the training half-hours. The everyday range gets fine resolution, the extremes get one wide bucket each way.

NoteThe 16 buckets, in dollars

# The bucket boundaries are quantiles of the pooled training prices, so every bucket holds the
# same number of half-hours; the averages below are the prices the solver plans with.
finite_edges = np.quantile(np.concatenate([prices_by_year[y].to_numpy() for y in TRAIN_YEARS]),
                           np.linspace(0, 1, price_model.n_bins + 1))
bucket_rows = []
for k in range(price_model.n_bins):
    if k == 0:
        span = f"below ${finite_edges[1]:,.0f}"
    elif k == price_model.n_bins - 1:
        span = f"${finite_edges[k]:,.0f} and above"
    else:
        span = f"${finite_edges[k]:,.0f} to ${finite_edges[k+1]:,.0f}"
    bucket_rows.append({"Bucket": k + 1, "Range ($/MWh)": span,
                        "Average ($/MWh)": f"${price_model.rep_price[k]:,.0f}"})
pd.DataFrame(bucket_rows).set_index("Bucket")

	Range ($/MWh)	Average ($/MWh)
Bucket
1	below $6	$-26
2	$6 to $38	$24
3	$38 to $55	$47
4	$55 to $62	$58
5	$62 to $70	$66
6	$70 to $78	$74
7	$78 to $86	$82
8	$86 to $93	$89
9	$93 to $103	$98
10	$103 to $115	$109
11	$115 to $129	$122
12	$129 to $150	$139
13	$150 to $181	$164
14	$181 to $240	$209
15	$240 to $300	$277
16	$300 and above	$647

Second, the counting. Everything in the fitted model is estimated by sorting and counting, and the three figures below walk one slice of that estimation end to end: the step from 7:00 am to 7:30 am. Start with a single morning. The model’s raw material is a pair of prices, where the price sat at 7:00 am and where it sat one half-hour later; the morning below contributes exactly one such pair, and everything else about the day is ignored.

# One real training morning, zoomed to the window around 7:00 am. The model's raw material
# is the PAIR (price at 7:00 am, price at 7:30 am); the rest of the morning is context it
# ignores. The morning is chosen deterministically: among mornings whose 7:00 am price sits
# in the busiest 7:00 am bucket and whose 7:30 am price lands one bucket lower, take the
# median morning by 7:00 am price (a typical example of the commonest starting point).
BELIEF = "#7b52ab"                   # purple, reserved for the forecaster's beliefs and results

slot_of_obs = train_days.index.hour * 2 + train_days.index.minute // 30
at_700, at_730 = train_days[slot_of_obs == 14], train_days[slot_of_obs == 15]  # 7:00 am / 7:30 am
assert (at_700.index.normalize() == at_730.index.normalize()).all()   # every day has both, in order
pairs = pd.DataFrame({"now": at_700.to_numpy(), "next": at_730.to_numpy()},
                     index=at_700.index.normalize())
assert len(pairs) == 1096, "expected one (7:00, 7:30) pair per training day"
bin_now  = price_model.to_bin(pairs["now"].to_numpy())
bin_next = price_model.to_bin(pairs["next"].to_numpy())
FOCUS = int(np.bincount(bin_now, minlength=price_model.n_bins).argmax())  # busiest 7:00 am bucket

falls = np.flatnonzero((bin_now == FOCUS) & (bin_next == FOCUS - 1))  # mornings stepping one bucket down
assert falls.size > 0
example = pairs.index[falls[np.argsort(pairs["now"].to_numpy()[falls])[len(falls) // 2]]]
morning = train_days.loc[example + pd.Timedelta(hours=5.5): example + pd.Timedelta(hours=9)]
hour_x = morning.index.hour + morning.index.minute / 60
p_now, p_next = pairs.loc[example, "now"], pairs.loc[example, "next"]

fig, ax = plt.subplots(figsize=(8, 4.2))
for x0, shade in ((7.0, 0.24), (7.5, 0.12)):                # the two half-hours of the pair
    ax.fill_between([x0, x0 + 0.5], 0, 1, transform=ax.get_xaxis_transform(),
                    color=BELIEF, alpha=shade)
ax.step(hour_x, morning.to_numpy(), where="post", color=INK, lw=2.2)
ax.plot([7.25, 7.75], [p_now, p_next], "o", color=BELIEF, ms=7, zorder=3)
def signed_dollars(v):
    '''$17, or -$10 for a negative price (never the awkward $-10).'''
    return f"-${abs(v):,.0f}" if v < 0 else f"${v:,.0f}"

ax.annotate(f"7:00 am: {signed_dollars(p_now)}", xy=(7.25, p_now), xytext=(6.0, p_now + 28),
            color=BELIEF, fontsize=10, fontweight="bold",
            arrowprops=dict(arrowstyle="->", color=BELIEF))
ax.annotate(f"7:30 am: {signed_dollars(p_next)}", xy=(7.75, p_next), xytext=(8.15, p_next + 34),
            color=BELIEF, fontsize=10, fontweight="bold",
            arrowprops=dict(arrowstyle="->", color=BELIEF))
ax.set(title="One Price Transition (7:00 to 7:30 am)", xlabel="Time of day", ylabel="Price ($/MWh)")
ax.set_xticks([6, 6.5, 7, 7.5, 8, 8.5, 9])
ax.set_xticklabels(["6am", "6:30", "7am", "7:30", "8am", "8:30", "9am"])
ax.set_xlim(5.5, 9)
ax.yaxis.set_major_formatter(dollars)
fig.text(0.005, -0.01, f"The morning of {example:%-d %B %Y}. Each half-hour has one settled price; "
         "the pair the model records is the shaded one.",
         ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "08a2-one-transition.png")

Now all of them at once. Each of the 1,096 training days contributes its own pair, drawn below as one line from its 7:00 am price to its 7:30 am price. On the future side the sixteen fitted buckets take over from exact dollars: all the model asks of 7:30 am is which bucket the price landed in. On the current side one bucket is singled out, the busiest at 7:00 am ($6 to $38): the purple bundle holds the 128 mornings that began there, with the morning above as the bold line among them.

# Every training morning as one faint line from its 7:00 am price to its 7:30 am price.
# The sixteen bucket edges divide the right (future) side; the busiest 7:00 am bucket is
# shaded on the left, and the mornings that start inside it (the bundle the next figure
# counts) are drawn in purple, with the example morning above in bold.
in_focus = bin_now == FOCUS
assert FOCUS == 1 and in_focus.sum() == 128, "the busiest 7:00 am bucket moved; update the prose"
lo, hi = price_model.edges[FOCUS], price_model.edges[FOCUS + 1]     # the $6 to $38 bucket

fig, ax = plt.subplots(figsize=(8, 5))
for sel, colour, alpha, lw in ((~in_focus, MUTED, 0.13, 0.6), (in_focus, BELIEF, 0.3, 0.8)):
    ax.plot([0, 1], np.vstack([pairs["now"][sel], pairs["next"][sel]]),
            color=colour, alpha=alpha, lw=lw)
ax.plot([0, 1], [p_now, p_next], color=BELIEF, lw=2.6, zorder=4)    # the morning above, bold
for e in price_model.edges[1:-1]:                                   # the bucket edges, future side
    ax.plot([1.01, 1.05], [e, e], color=INK2, lw=1, clip_on=False)
    ax.axhline(e, color=INK2, lw=0.5, ls=":", alpha=0.18)
for k in (1, 5, 12, 15):                                            # a few edges labelled in dollars
    ax.text(1.07, price_model.edges[k], f"${price_model.edges[k]:,.0f}",
            fontsize=7.5, color=INK2, va="center", clip_on=False)
ax.fill_betweenx([lo, hi], -0.06, 0.01, color=BELIEF, alpha=0.3, clip_on=False)
ax.text(-0.09, (lo + hi) / 2, "the \\$6 to \\$38 bucket:\n128 mornings", ha="right", va="center",
        color=BELIEF, fontsize=8.5)
n_clipped = int(((pairs["now"] > 340) | (pairs["next"] > 340)).sum())
ax.set(title="All 1,096 Training Mornings (7:00 to 7:30 am)", ylabel="Price ($/MWh)")
ax.set_xticks([0, 1])
ax.set_xticklabels(["7:00 am", "7:30 am"])
ax.set_xlim(-0.5, 1.22)
ax.set_ylim(-80, 340)
ax.yaxis.set_major_formatter(dollars)
ax.grid(False)
fig.text(0.005, -0.01, f"Bucket edges marked on the right. {n_clipped} spike mornings run off the top "
         "of the frame; the outermost buckets extend to the data's extremes.",
         ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "08a3-all-mornings.png")

Counting is all that remains. The third figure follows the purple bundle. Of the 128 mornings that began in the $6 to $38 bucket, 75 were still there at 7:30 am, 33 had slipped into the bucket below, 14 had climbed one bucket, and a handful moved further. Divide each count by 128 and the flows become probabilities: one full column of the 7:00 am transition table.

# Left: the 128 mornings fan out from their starting bucket to the bucket they landed in
# at 7:30 am (ribbon width = count). Right: each count divided by 128, the estimated odds;
# the black ticks show the same column AFTER the add-one smoothing the fitted table applies,
# which is why the heatmaps below differ slightly from the raw shares.
flow = np.bincount(bin_next[in_focus], minlength=price_model.n_bins)
assert (flow[FOCUS - 1], flow[FOCUS], flow[FOCUS + 1]) == (33, 75, 14), "counts moved; update the prose"
share = flow / flow.sum()
fitted_col = price_model.transitions[14, FOCUS, :]                  # the smoothed column the solver uses

def bucket_label(k):
    '''Y-axis label for bucket k: its number and dollar range. Dollar signs are escaped so
    matplotlib never reads a $...$ pair as mathtext.'''
    e = price_model.edges
    if k == 0:
        return f"1: below \\${e[1]:,.0f}"
    if k == price_model.n_bins - 1:
        return f"{k + 1}: \\${e[k]:,.0f} and above"
    return f"{k + 1}: \\${e[k]:,.0f} to \\${e[k + 1]:,.0f}"

fig, (flow_ax, odds_ax) = plt.subplots(1, 2, figsize=(9, 5.2), sharey=True,
                                       gridspec_kw={"width_ratios": [1.5, 1]})
# ribbons in bucket-index space: a smooth cosine ramp from the source bucket's node (its
# segments stacked bottom-to-top in destination order) to each destination bucket
xs_ramp = np.linspace(0.06, 0.9, 60)
ramp = (1 - np.cos(np.pi * (xs_ramp - xs_ramp[0]) / (xs_ramp[-1] - xs_ramp[0]))) / 2
seg_bottom = FOCUS - 0.42                                           # the source node spans +-0.42 buckets
for k in np.flatnonzero(flow):
    h = share[k] * 0.84                                             # ribbon width, in bucket units
    y_src = seg_bottom + np.array([0.0, h])
    seg_bottom += h
    y_dst = np.array([k - h / 2, k + h / 2])
    flow_ax.fill_between(xs_ramp, y_src[0] + (y_dst[0] - y_src[0]) * ramp,
                         y_src[1] + (y_dst[1] - y_src[1]) * ramp, color=BELIEF, alpha=0.42, lw=0)
    flow_ax.fill_between([0.9, 0.96], k - h / 2, k + h / 2, color=BELIEF, alpha=0.85)  # arrival stub
    if flow[k] >= 3:
        flow_ax.text(0.88, k, f"{flow[k]}", ha="right", va="center", fontsize=8.5, color=INK)
flow_ax.fill_between([0.0, 0.06], FOCUS - 0.42, FOCUS + 0.42, color=BELIEF, alpha=0.95)  # source node
flow_ax.text(-0.03, FOCUS, "128\nmornings", ha="right", va="center",
             color=BELIEF, fontsize=8.5, fontweight="bold")
flow_ax.set_yticks(range(price_model.n_bins))
flow_ax.set_yticklabels([bucket_label(k) for k in range(price_model.n_bins)], fontsize=7.4)
flow_ax.set_xticks([0.03, 0.93])
flow_ax.set_xticklabels(["7:00 am", "7:30 am"])
flow_ax.set_xlim(-0.28, 1.0)
flow_ax.set_ylim(-0.6, price_model.n_bins - 0.4)
flow_ax.grid(False)
flow_ax.set_title("Where the 128 mornings went", fontsize=10.5)
odds_ax.barh(range(price_model.n_bins), share * 100, color=BELIEF, alpha=0.85, height=0.72)
odds_ax.plot(fitted_col * 100, range(price_model.n_bins), marker="|", ms=9, ls="none", color=INK)
for k in np.flatnonzero(share >= 0.02):                             # label the visible odds
    odds_ax.text(share[k] * 100 + 1.2, k, f"{share[k] * 100:.0f}%", va="center", fontsize=8.5, color=INK)
odds_ax.set_xlabel("Probability (%)")
odds_ax.set_xlim(0, 68)
odds_ax.grid(axis="x")
odds_ax.set_title("Divided by 128: the odds", fontsize=10.5)
odds_ax.legend(handles=[plt.Line2D([], [], marker="|", ms=9, ls="none", color=INK,
                                   label="After add-one smoothing")],
               frameon=False, loc="upper right", fontsize=7.5)
fig.suptitle("From Counts to Odds (7:00 am)", fontsize=13, fontweight="bold")
fig.text(0.005, -0.01, "The fitted table adds one phantom observation to every cell (Laplace smoothing) "
         "before dividing, so moves never seen in three years keep small positive odds.",
         ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "08a4-counts-to-odds.png")

Assembled for every starting bucket, those columns of odds become one 16-by-16 transition table for the 7:00 am slot; repeated for each half-hour of the day, they become the model’s 48 tables. Two are drawn below; to read them, pick the current price on the horizontal axis, and the column above it gives the odds of each price one half-hour later (the column just built is the second from the left in the first panel). The dark diagonal in both says prices are sticky: the single most likely next price is almost always the current one. The lean around the diagonal is the forecast: at 7:00 am the mass sits below the diagonal (cheaper next), because solar is about to flood the market; at 4:00 pm it sits above the diagonal (more expensive next), because the evening peak is building. That lean, estimated separately for every half-hour of the day, is all the “forecasting” the forecaster ever does.

# Each panel is one fitted transition table, drawn with the CURRENT price on the x axis and
# the NEXT half-hour's price on the y axis, shade = fitted probability (so each vertical
# column of cells is one distribution over what happens next). A single purple ramp (light =
# rare, dark = likely); purple is reserved throughout for the forecaster's beliefs and results.
from matplotlib.colors import LinearSegmentedColormap
belief_ramp = LinearSegmentedColormap.from_list("belief", [SURFACE, BELIEF])

show_slots = ((14, "About to fall", "7:00 am", "7:30 am"),
              (32, "About to climb", "4:00 pm", "4:30 pm"))
vmax = max(price_model.transitions[t].max() for t, *_ in show_slots)
tick_buckets = [0, 4, 8, 12, 15]
fig, axes = plt.subplots(1, 2, figsize=(9, 4.7))
fig.subplots_adjust(wspace=0.34)              # room for the right panel's own y label
for ax, (t, label, now_time, next_time) in zip(axes, show_slots):
    mesh = ax.pcolormesh(np.arange(17), np.arange(17), price_model.transitions[t].T,
                         cmap=belief_ramp, vmin=0, vmax=vmax)     # transposed: columns = now
    ax.plot([0, 16], [0, 16], color=INK2, ls=":", lw=1)          # the stay-put diagonal
    ax.set_xticks([b + 0.5 for b in tick_buckets])
    ax.set_xticklabels([f"${price_model.rep_price[b]:,.0f}" for b in tick_buckets], fontsize=8)
    ax.set_yticks([b + 0.5 for b in tick_buckets])
    ax.set_yticklabels([f"${price_model.rep_price[b]:,.0f}" for b in tick_buckets], fontsize=8)
    ax.grid(False)
    ax.set_title(label, fontsize=10.5)
    ax.set_xlabel(f"Price at {now_time} ($/MWh)")
    ax.set_ylabel(f"Price at {next_time} ($/MWh)")
fig.suptitle("Transition Probability for Prices", fontsize=13, fontweight="bold")
cbar = fig.colorbar(mesh, ax=axes, shrink=0.85, pad=0.02)
cbar.set_label("Fitted probability")
finish(fig, "08b-transition-odds.png")

_slot, _price = 35, 100.0                                 # 5:30 pm, with the price near $100
_k = int(price_model.to_bin(np.array([_price]))[0])
_row = price_model.transitions[_slot, _k]
display(Markdown(
    f"One slice of the model, read aloud: at **5:30 pm**, with the price in the **\\${price_model.rep_price[_k]:,.0f} "
    f"bucket**, the fitted odds are **{_row[_k]*100:.0f}% stay put**, **{_row[_k+1:].sum()*100:.0f}% move up** "
    f"(including {_row[_k+2:].sum()*100:.0f}% by two buckets or more: the evening ramp), and "
    f"**{_row[:_k].sum()*100:.0f}% move down**. The operator never knows which will happen. It only ever "
    f"holds these odds, and the decision rule below is the best possible standing response to them."
))

One slice of the model, read aloud: at 5:30 pm, with the price in the $98 bucket, the fitted odds are 23% stay put, 46% move up (including 33% by two buckets or more: the evening ramp), and 31% move down. The operator never knows which will happen. It only ever holds these odds, and the decision rule below is the best possible standing response to them.

The estimated decision rule

The whole strategy compresses into one picture. The map below fixes the clock at 6 pm and shows the chosen action at every store level and price bucket: pump (blue) when the price is low, generate (orange) when it is high, and wait (grey) in the band between, where the price gap cannot cover the round-trip loss. The two boundaries are the water value at work, and they slide with the store level exactly as intuition says they should. Nearly empty, stored energy is precious: the rule buys at prices up to about $140/MWh and refuses to sell below the $200 price rows. Nearly full, stored energy is cheap to it: the rule only buys below about $70 and starts selling from about $100. Nothing here was told to do that; it fell out of the backward induction on its own.

One note for reading the vertical axis: the 16 price buckets each hold the same number of historical half-hours, not the same dollar width, so the top rows span far wider price ranges (the top row, labelled $647, covers every price from $300 up) than the tightly packed everyday rows in the middle. The labels mark each bucket’s average price.

from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

slot_6pm = 36                                                     # 6 pm; the map barely changes with the clock
action_map = learned["policy"][slot_6pm, :, :].T                  # rows = price buckets, cols = store levels
action_colours = ListedColormap(["#d8d6cf", BUY, SELL])           # 0 idle, 1 pump, 2 generate
store_gwh = learned["soc_grid"] * snowy.eta_gen / 1e3             # store axis in generatable GWh (Snowy's units)
step = store_gwh[1] - store_gwh[0]
store_edges = np.concatenate([[store_gwh[0] - step/2], store_gwh + step/2])   # cell edges for the mesh

fig, ax = plt.subplots(figsize=(9, 4.6))
ax.pcolormesh(store_edges, np.arange(17), action_map, cmap=action_colours, vmin=-0.5, vmax=2.5)
label_buckets = [0, 4, 8, 12, 15]                                 # a few representative price levels
ax.set_yticks([b + 0.5 for b in label_buckets])
ax.set_yticklabels([f"${price_model.rep_price[b]:,.0f}" for b in label_buckets])
ax.grid(False)
ax.set(title="Forecaster Decision Rule, 6 pm",
       xlabel="State of charge (GWh)", ylabel="Price level ($/MWh)")
ax.legend(handles=[Patch(color=BUY, label="Pump (buy)"), Patch(color="#d8d6cf", label="Wait"),
                   Patch(color=SELL, label="Generate (sell)")],
          frameon=False, loc="upper left", bbox_to_anchor=(0, -0.18), ncol=3)
finish(fig, "09-learned-rule.png")

One result we did not impose: re-drawing this map for any other hour of the day barely moves it (about 2% of its cells change). Snowy’s reservoir is so deep, roughly 160 hours at full power, that no single day can fill or drain much of it, so the clock hardly matters and the store level is everything. In effect the learned strategy is a pair of price triggers that adapt to how full the store is. That is also exactly where Scheme 2 leaves money on the table: its two triggers are frozen for the whole year, while this rule slides them with the store level and places them where the water value says they belong. The next figure makes the clock’s near-irrelevance visible directly, and shows it is a property of Snowy, not of the method.

Does the clock change the rule?

A natural intuition about storage says it should: you might pump at a higher price in the early afternoon than in the evening, because the evening is expected to be more expensive still. The figure below tests that intuition, one panel per plant. At every half-hour of the day, holding each store at half of its own capacity, it reads two prices off the fitted rule: the highest price at which the rule still pumps (blue) and the lowest at which it generates (orange). Between them sits the wait zone (grey), the prices at which the rule does nothing; prices below the blue line trigger pumping, prices above the orange line trigger generating. For Snowy (left) the two triggers barely move all day: with about 160 hours of storage behind every decision, no single half-hour’s outlook can change what a stored megawatt-hour is worth, so the time of day all but drops out of the rule. To show that this is the reservoir’s depth and not a limitation of the method, the right panel re-solves the identical problem for a store holding just four hours at full power: there the clock is everything, exactly as the intuition says.

def triggers_by_slot(rule, model, soc_frac=0.5):
    '''For each half-hour slot of the day, a rule's two working triggers at a given store
    fullness: the highest bucket average at which it still pumps, and the cheapest at which
    it generates. NaN where the rule takes no such action at any price in that slot.'''
    soc_ix = int(round(soc_frac * (len(rule["soc_grid"]) - 1)))
    buy_t, sell_t = np.full(model.slots, np.nan), np.full(model.slots, np.nan)
    for t in range(model.slots):
        actions = rule["policy"][t, soc_ix, :]
        pumps = np.flatnonzero(actions == 1)
        gens  = np.flatnonzero(actions == 2)
        if pumps.size:
            buy_t[t] = model.rep_price[pumps.max()]
        if gens.size:
            sell_t[t] = model.rep_price[gens.min()]
    return buy_t, sell_t

# The same beliefs and solver, but a store with four hours at full power instead of ~160.
battery_depth = replace(snowy, dischargeable_mwh=4 * snowy.power_mw)
learned_battery = solve_known_distribution(price_model, battery_depth)
assert learned_battery["converged"], "the battery-depth value iteration hit its sweep cap"

snowy_buy, snowy_sell = triggers_by_slot(learned, price_model)
batt_buy,  batt_sell  = triggers_by_slot(learned_battery, price_model)
slot_hours = np.arange(price_model.slots) / 2

fig, axes = plt.subplots(1, 2, figsize=(9, 4.4), sharey=True)
panels = ((axes[0], snowy_buy, snowy_sell, "Snowy 2.0: ~160 hours of storage"),
          (axes[1], batt_buy, batt_sell, "Same power, 4 hours of storage"))
for ax, buy_t, sell_t, label in panels:
    ax.fill_between(slot_hours, buy_t, sell_t, step="post", color="#d8d6cf", alpha=0.55)
    ax.step(slot_hours, sell_t, where="post", color=SELL, lw=2.4)
    ax.step(slot_hours, buy_t, where="post", color=BUY, lw=2.4)
    ax.set_title(label, fontsize=10.5)
    ax.set_xlabel("Time of day")
    clock_axis(ax, step=6)
axes[0].set_ylim(0, 300)                          # shared; the battery's sell trigger exits on purpose
axes[0].yaxis.set_major_formatter(dollars)
axes[0].set_ylabel("Price ($/MWh)")
# direct labels on the left panel carry the reading for both
axes[0].text(12, 150, "generates at or above", color=SELL, fontsize=9, ha="center")
axes[0].text(12, 100, "waits", color=INK2, fontsize=9, ha="center")
axes[0].text(12, 55, "pumps at or below", color=BUY, fontsize=9, ha="center")
# the battery's story, annotated where it happens (the peak pump trigger, read live); the
# two notes sit in the panel's empty corners so their arrows stay clear of each other
_pk = int(np.nanargmax(batt_buy))
axes[1].annotate(f"pays up to ${np.nanmax(batt_buy):,.0f} to fill\nahead of the peak",
                 xy=(slot_hours[_pk], batt_buy[_pk]), xytext=(0.6, 255), va="top",
                 color=BUY, fontsize=8.5, arrowprops=dict(arrowstyle="->", color=BUY))
axes[1].annotate("sell trigger leaves the chart:\nnothing is for sale before the peak",
                 xy=(15.1, 293), xytext=(8.8, 22), color=SELL, fontsize=8.5,
                 arrowprops=dict(arrowstyle="->", color=SELL))
fig.suptitle("Trigger Prices Through the Day", fontsize=13, fontweight="bold")
fig.text(0.005, -0.01, "Each store is read at half of its own capacity: about 80 hours of energy for "
         "Snowy, two for the battery. Each trigger is the fitted rule's marginal price "
         "bucket: the highest at which it still pumps, the lowest at which it generates.",
         ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "09b-triggers-by-clock.png")

The contrast carries the answer. A four-hour store lives inside a single day, so its rule is all clock: by mid-afternoon it will pay about $164/MWh to top itself up ahead of the peak, roughly twice what it pays overnight, and through the afternoon it refuses to sell at any everyday price at all (its sell trigger jumps to the extreme-spike bucket, off the top of the chart), hoarding its four hours for the evening. Once the peak passes, both triggers collapse back. Snowy, holding nearly a week of energy at full power, pumps below about $90 and generates above about $120 at every hour of the day; its wait zone is a flat, narrow band. Depth substitutes for timing: Snowy’s rule needs to know how full the store is, not what time it is. This is also why Scheme 2’s frozen dollar-triggers get as close as they do for Snowy, and why the same two-trigger idea would fail badly for a battery.

Performance in and out of sample

Now comes the real test. We let the rule trade each of the four years, paying and earning the actual prices, not its buckets. Three of those years it trained on; 2025 it has never seen.

Scheme 3 is scored under a repeated-year accounting convention, mirroring the cyclic year already imposed on Scheme 1: for each year we choose a starting reservoir level that equals the level the year ends at (within a stated tolerance of 0.1 GWh of dischargeable energy) when the rule trades that year’s realised prices. This removes uncompensated inventory build-up from the annual profit comparison; without it, a year that quietly banks a nearly full reservoir gets no credit for the stored value. The convention uses the full year only to normalise the boundary: every dispatch decision within the year remains a causal function of the current time, price bucket and store level.

One disclosure follows. The 2025 decision rule is fully held out, trained only on 2022 to 2024, but the boundary level is found ex post, from the 2025 path itself: it is a property of the rule combined with the year it trades, not of the rule alone. The diagnostics after the table show what the convention is worth, by also scoring 2025 from an arbitrary empty store.

(Figures in the table are rounded to two decimals, so the gap column can differ slightly from subtracting the rounded columns by hand.)

# Score the stationary rule on each year under the repeated-year convention. cyclic_score
# searches for the start = end reservoir level from five different starting levels, so a
# non-unique fixed point or a failure to settle cannot pass silently (both are reported).
scheme3_by_year, comparison_rows = {}, []
for year in YEARS:
    year_prices = prices_by_year[year].to_numpy()
    outcome = cyclic_score(year_prices, snowy, learned["policy"], learned["soc_grid"], price_model)
    scheme3_by_year[year] = outcome
    ceiling = ceiling_by_year[year]["profit"]
    comparison_rows.append({
        "Year": f"{year}{'*' if year not in TRAIN_YEARS else ''}",
        "Rule of thumb (S2), $bn": round(s2_by_year[year] / 1e9, 2),
        "Forecaster (S3), $bn": round(outcome["profit"] / 1e9, 2),
        "Perfect foresight (S1), $bn": round(ceiling / 1e9, 2),
        "Forecaster share of S1": f"{outcome['profit'] / ceiling * 100:.0f}%",
        "Foresight gap ($bn)": round((ceiling - outcome["profit"]) / 1e9, 2),
        "Start = end store (GWh)": round(outcome["start_soc"] * snowy.eta_gen / 1e3, 1)})
comparison = pd.DataFrame(comparison_rows).set_index("Year")
comparison

	Rule of thumb (S2), $bn	Forecaster (S3), $bn	Perfect foresight (S1), $bn	Forecaster share of S1	Foresight gap ($bn)	Start = end store (GWh)
Year
2022	0.27	0.47	0.89	53%	0.42	258.5
2023	0.44	0.60	0.71	84%	0.11	347.8
2024	1.02	1.32	1.43	92%	0.12	322.7
2025*	0.90	1.04	1.13	93%	0.08	347.5

* Out of sample: the two trained schemes commit to their rules on 2022 to 2024 before trading 2025.

held_out = scheme3_by_year[2025]["profit"] / ceiling_by_year[2025]["profit"] * 100
gap_2025 = (ceiling_by_year[2025]["profit"] - scheme3_by_year[2025]["profit"]) / 1e9
gap_2022 = (ceiling_by_year[2022]["profit"] - scheme3_by_year[2022]["profit"]) / 1e9
# Diagnostics for the accounting convention: the same rule scored from an arbitrary empty
# store (the naive boundary), and the boundary-search record for the held-out year.
empty_2025 = run_policy(prices_by_year[2025].to_numpy(), snowy,
                        learned["policy"], learned["soc_grid"], price_model)   # start_soc = 0
b = scheme3_by_year[2025]
assert b["converged"], "the 2025 boundary search did not meet its tolerance"
agreement = ("reaching the same level from all five starting levels tried"
             if b["unique_fixed_point"] else
             f"with the five starting levels tried settling within a "
             f"{b['basin_dis_mwh'] / 1e3:.1f} GWh band whose profits agree to "
             f"{b['profit_spread'] * 100:.2f}%")
display(Markdown(
    f"The headline is the held-out year: trading 2025 with a rule learned entirely from 2022 to "
    f"2024, the operator captures **{held_out:.0f}% of the perfect-foresight ceiling**. A perfect "
    f"forecast of every 2025 price would have added only **\\${gap_2025:.2f} bn** (the realised "
    f"perfect-foresight gap). The exception is 2022: the fuel crisis pushed prices far outside "
    f"any learned pattern, and the gap blows out to **\\${gap_2022:.2f} bn**. Foresight is worth "
    f"little in a normal year and a fortune in a crisis, because what it really buys is advance "
    f"warning of extreme years.\n\n"
    f"Two diagnostics for the accounting convention. Scored from an arbitrary empty store "
    f"instead, the 2025 profit is \\${empty_2025['profit'] / 1e9:.2f} bn: the difference is "
    f"purely the boundary, never the trading. And the 2025 boundary search settled at "
    f"{b['start_soc'] * snowy.eta_gen / 1e3:.1f} GWh in {b['iterations']} passes with a "
    f"start-to-end residual of {b['residual_dis_mwh']:.0f} MWh (tolerance 100 MWh), {agreement}."
))

The headline is the held-out year: trading 2025 with a rule learned entirely from 2022 to 2024, the operator captures 93% of the perfect-foresight ceiling. A perfect forecast of every 2025 price would have added only $0.08 bn (the realised perfect-foresight gap). The exception is 2022: the fuel crisis pushed prices far outside any learned pattern, and the gap blows out to $0.42 bn. Foresight is worth little in a normal year and a fortune in a crisis, because what it really buys is advance warning of extreme years.

Two diagnostics for the accounting convention. Scored from an arbitrary empty store instead, the 2025 profit is $1.01 bn: the difference is purely the boundary, never the trading. And the 2025 boundary search settled at 347.5 GWh in 4 passes with a start-to-end residual of 88 MWh (tolerance 100 MWh), reaching the same level from all five starting levels tried.

The same result as one picture: for each year, the grey square is the mechanical rule of thumb, the purple dot is the forecaster, and the open circle is perfect foresight (the ceiling). The distance from the purple dot up to the circle is the realised perfect-foresight gap.

RESULT = "#7b52ab"    # purple: the forecaster's results, here and in the figures below, the
                      # same hue as its beliefs. Green stays reserved for stored energy, so no
                      # colour changes meaning between figures.
# One dumbbell per year: rule of thumb (grey square) up to the ceiling (open circle), with
# the forecaster's dot between them; the annotated arrow is the perfect-foresight gap.
fig, ax = plt.subplots(figsize=(9, 4.8))
for i, year in enumerate(YEARS):
    bench_bn   = s2_by_year[year] / 1e9
    real_bn    = scheme3_by_year[year]["profit"] / 1e9
    ceiling_bn = ceiling_by_year[year]["profit"] / 1e9
    ax.plot([i, i], [bench_bn, ceiling_bn], color=GRID, lw=2.5, zorder=1)          # the year's range
    ax.plot(i, bench_bn,   marker="s", color=MUTED, ms=9,  zorder=2, ls="none")
    ax.plot(i, ceiling_bn, marker="o", mfc="white", mec=INK, mew=2, ms=11, zorder=2, ls="none")
    ax.plot(i, real_bn,    marker="o", color=RESULT, ms=12, zorder=3, ls="none")
    ax.annotate(f"${real_bn:.2f}bn", (i, real_bn), xytext=(14, -4), textcoords="offset points",
                color=RESULT, fontsize=9.5, fontweight="bold")
# label the perfect-foresight gap on the year where it is largest (the 2022 crisis)
gap_top = ceiling_by_year[2022]["profit"] / 1e9
gap_bottom = scheme3_by_year[2022]["profit"] / 1e9
ax.annotate("", xy=(0.18, gap_top), xytext=(0.18, gap_bottom),
            arrowprops=dict(arrowstyle="<->", color=SPIKE, lw=1.6))
ax.text(0.26, (gap_top + gap_bottom) / 2, "Perfect-foresight gap:\nwhat a perfect forecast\nwould have added",
        color=SPIKE, fontsize=9, va="center")
ax.set_xticks(range(len(YEARS)))
ax.set_xticklabels([f"{y}*" if y not in TRAIN_YEARS else str(y) for y in YEARS])
ax.set(title="Rule of Thumb, Forecaster, and Perfect Foresight by Year", ylabel="Annual arbitrage profit ($bn)")
ax.yaxis.set_major_formatter(lambda x, _: f"${x:,.1f}")
ax.legend(handles=[plt.Line2D([], [], marker="s", color=MUTED, ls="none", ms=9, label="Rule of thumb (S2)"),
                   plt.Line2D([], [], marker="o", color=RESULT, ls="none", ms=11, label="Forecaster (S3)"),
                   plt.Line2D([], [], marker="o", mfc="white", mec=INK, mew=2, ls="none", ms=10, label="Perfect foresight (S1)")],
          frameon=False, loc="upper left")
fig.text(0.985, 0.02, "* 2025 is out of sample: the trained schemes commit to their rules on 2022 to 2024",
         ha="right", fontsize=8, color=INK2, style="italic")
finish(fig, "10-scheme-comparison.png")

NoteConvergence evidence

Trusting a rule produced by an iterative solve needs two dials to have stopped moving: the solver’s internal one (the value shape, which is its stopping rule) and the one that matters to a reader (the scored result). The table shows the held-out 2025 score for snapshots of the rule taken part-way through the solve, and for a fresh solve pushed to twice as many sweeps as the stopping rule asked for. The profit settles to well inside 0.1% by the retained rule, and solving on past the stopping point buys nothing.

# Convergence of the answer, not just of the solver: score each mid-solve snapshot of the
# rule on held-out 2025 under the repeated-year convention, then push a fresh solve to twice
# the sweeps and confirm the retained rule's number no longer moves.
p25 = prices_by_year[2025].to_numpy()
snapshots = [*sorted(learned["checkpoints"].items()), (learned["sweeps"], learned["policy"])]
overrun = solve_known_distribution(price_model, snowy, span_tol=0.0,
                                   max_sweeps=2 * learned["sweeps"])
snapshots.append((2 * learned["sweeps"], overrun["policy"]))
conv_profits = [cyclic_score(p25, snowy, rule, learned["soc_grid"], price_model)["profit"]
                for _, rule in snapshots]
# The second half of the dual stopping criterion, asserted: the retained rule (second-last
# row) and the double-length solve (last row) must agree to better than 0.1%.
assert abs(conv_profits[-1] - conv_profits[-2]) / conv_profits[-2] < 0.001, \
    "the scored profit is still moving at the stopping point; raise the sweep budget"
convergence = pd.DataFrame({
    "One-day sweeps": [s for s, _ in snapshots],
    "2025 profit ($bn)": [round(v / 1e9, 3) for v in conv_profits],
    "Move vs previous row": ["", *(f"{abs(b - a) / a * 100:.2f}%"
                                   for a, b in zip(conv_profits, conv_profits[1:]))],
}).set_index("One-day sweeps")
convergence

	2025 profit ($bn)	Move vs previous row
One-day sweeps
14	1.098
28	1.069	2.59%
56	1.047	2.05%
112	1.044	0.30%
178	1.043	0.07%
356	1.043	0.01%

Decomposing the forecaster’s profit

The comparison above says how much the forecaster earned; the two views below open those earnings up. First by year: each bar is the year’s generation revenue, split into what was paid for pumping energy, variable operating cost, and what was left as profit. Two things stand out. The pumping bill is nearly constant, about $0.3 bn in every year, while the revenue side swings by a factor of two; years differ on what the plant sells for, not on what it pays. And the margin is fat: between 61% and 79% of every revenue dollar survives as profit, because the plant buys its energy at a small fraction of its selling price.

# Decompose each forecaster year from its own scored dispatch: generation revenue = pumping
# cost + variable O&M + profit. The identity is asserted against the year's reported profit,
# so the segments can never drift from the headline numbers.
dt = snowy.interval_hr
decomp_rows = []
for year in YEARS:
    out = scheme3_by_year[year]
    px = prices_by_year[year].to_numpy()
    revenue = float((out["gen"] * px).sum() * dt)                # dollars earned generating
    pump_cost = float((out["pump"] * px).sum() * dt)             # dollars paid to pump
    var_om = float(snowy.var_om_per_mwh * out["gen"].sum() * dt) # variable O&M on generation
    profit = revenue - pump_cost - var_om
    assert abs(profit - out["profit"]) < 1e-3 * abs(out["profit"]), f"{year} decomposition drifted"
    decomp_rows.append({"year": year, "revenue": revenue, "pump": pump_cost,
                        "om": var_om, "profit": profit})
decomp = pd.DataFrame(decomp_rows).set_index("year") / 1e9       # everything in $bn

fig, ax = plt.subplots(figsize=(9, 4.6))
xs = np.arange(len(YEARS))
ax.bar(xs, decomp["pump"], width=0.6, color=BUY, alpha=0.85,
       edgecolor=INK, linewidth=1.1, label="Pumping cost")
ax.bar(xs, decomp["om"], width=0.6, bottom=decomp["pump"], color=MUTED, alpha=0.9,
       edgecolor=INK, linewidth=1.1, label="Variable O&M")
ax.bar(xs, decomp["profit"], width=0.6, bottom=decomp["pump"] + decomp["om"], color=RESULT,
       alpha=0.9, edgecolor=INK, linewidth=1.1, label="Profit")
for i, year in enumerate(YEARS):                                 # label the takeaway numbers
    r = decomp.loc[year]
    ax.text(i, r["pump"] + r["om"] + r["profit"] / 2, f"${r['profit']:.2f}bn",
            ha="center", va="center", color="white", fontsize=9.5, fontweight="bold")
    ax.text(i, r["pump"] / 2, f"${r['pump']:.2f}bn",
            ha="center", va="center", color="white", fontsize=9.5, fontweight="bold")
    ax.text(i, r["revenue"] + 0.03, f"revenue ${r['revenue']:.2f}bn",
            ha="center", fontsize=8.5, color=INK2)
ax.set_xticks(xs)
ax.set_xticklabels([f"{y}*" if y not in TRAIN_YEARS else str(y) for y in YEARS])
ax.set(title="Forecaster Revenue: Decomposition", ylabel="$ billion per year",
       ylim=(0, decomp["revenue"].max() * 1.14))
ax.yaxis.set_major_formatter(lambda x, _: f"${x:,.1f}")
ax.legend(frameon=False, loc="upper left")
fig.text(0.985, 0.02, "* 2025 is out of sample: the rule commits on 2022 to 2024",
         ha="right", fontsize=8, color=INK2, style="italic")
finish(fig, "10b-forecaster-decomposition.png")

Second, by month inside the held-out year. Winter carries 2025: June alone contributes about $304m of the year’s $1,043m, essentially the entire perfect-foresight June ($305m), because in the month of the spike sequence the right move is obvious to foresight and forecaster alike: sell everything the store holds. The collapse to $22m in July is the other half of the same event, not a failure of the rule, and it has two parts. July 2025 was the mildest month of the year (its worst half-hour reached $395; every other month reached about $600 or more), so there was little worth selling into. And the plant entered July nearly empty, the store drained to 11.5 GWh by the end of June, so it spent the month rebuilding: 396 pumping hours, the year’s most, restored about 269 GWh, and the monthly ledger books that bill ($49m, the year’s largest) in July while the energy it bought earns its revenue only in the months that follow. In all three training years July out-earned June, so this is the June 2025 spikes at work, not a seasonal pattern. The hours columns also make the round-trip loss visible on the clock: the plant pumps roughly four hours for every three it generates, because almost a quarter of every stored megawatt-hour evaporates on the way through.

# The held-out year month by month, from the scored dispatch: how many hours the plant pumped
# and generated, what pumping cost, what generating earned, and the profit left over.
out25 = scheme3_by_year[2025]
px25 = prices_by_year[2025]
by_month = pd.DataFrame({"price": px25.to_numpy(), "pump": out25["pump"], "gen": out25["gen"]},
                        index=px25.index).groupby(px25.index.month)
monthly_2025 = pd.DataFrame({
    "Pumping hours": by_month.apply(lambda d: (d["pump"] > 0).sum() * dt),
    "Generating hours": by_month.apply(lambda d: (d["gen"] > 0).sum() * dt),
    "Pumping cost ($m)": by_month.apply(lambda d: (d["pump"] * d["price"]).sum() * dt / 1e6),
    "Variable O&M ($m)": by_month.apply(lambda d: snowy.var_om_per_mwh * d["gen"].sum() * dt / 1e6),
    "Generation revenue ($m)": by_month.apply(lambda d: (d["gen"] * d["price"]).sum() * dt / 1e6),
})
monthly_2025["Profit ($m)"] = (monthly_2025["Generation revenue ($m)"]
                               - monthly_2025["Pumping cost ($m)"] - monthly_2025["Variable O&M ($m)"])
monthly_2025.index = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
                      "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
monthly_2025.loc["Total"] = monthly_2025.sum()
monthly_2025.rename_axis("Month (2025)").style.format("{:,.0f}")

	Pumping hours	Generating hours	Pumping cost ($m)	Variable O&M ($m)	Generation revenue ($m)	Profit ($m)
Month (2025)
Jan	260	204	13	2	81	67
Feb	227	194	14	2	78	63
Mar	271	186	21	2	78	55
Apr	311	275	24	2	112	86
May	319	262	33	2	136	101
Jun	282	318	35	3	343	304
Jul	396	183	49	2	73	22
Aug	352	252	33	2	103	67
Sep	278	210	11	2	74	61
Oct	273	250	9	2	90	79
Nov	324	202	7	2	88	79
Dec	208	160	7	1	67	58
Total	3,502	2,696	256	24	1,323	1,043

The bars below repeat the table’s profit column; the green line adds the store itself, whose drain into June and refill through July is the mechanism just described. The lower panel prices each month’s trades and makes the same point in dollars: the plant sells at $160 to $240 a megawatt-hour in every ordinary month and at $490 in June, while July pairs an ordinary sell price ($181) with buying at $56, the dearest refill of the year bar June itself.

# Top panel: the profit column above as one bar per month, with the store level at each
# month's end on a second axis (the June drain and July refill that explain the swing).
# Bottom panel: the volume-weighted average price of each month's actual trades, one line
# per leg, which shows WHY: June sells at $490/MWh, July sells at an ordinary price while
# buying its refill at nearly the year's dearest.
fig, (ax, trade_ax) = plt.subplots(2, 1, figsize=(9, 6.4), sharex=True,
                                   gridspec_kw={"height_ratios": [2.0, 1.0], "hspace": 0.08})
months = monthly_2025.index[:-1]                                 # drop the Total row
profit_m = monthly_2025["Profit ($m)"].iloc[:-1]
ax.bar(range(12), profit_m, color=RESULT, alpha=0.9, edgecolor=INK, linewidth=1.1)
for i, v in enumerate(profit_m):                                 # the value on each bar
    ax.text(i, v + 6, f"${v:,.0f}m", ha="center", fontsize=8.5, color=INK)
ax.set(title="Forecaster Profit by Month (2025)", ylabel="Profit ($m)",
       ylim=(0, profit_m.max() * 1.14))
soc_end = (pd.Series(out25["stored"][:-1] * snowy.eta_gen / 1e3, index=px25.index)
           .groupby(px25.index.month).last())                    # store at each month's end, GWh
store_ax = ax.twinx()
store_ax.plot(range(12), soc_end.to_numpy(), color=STORE, lw=2.2, marker="o", ms=4.5)
store_ax.set_ylabel("State of charge,\nend of month (GWh)", color=STORE)
store_ax.set_ylim(bottom=0)
store_ax.tick_params(axis="y", colors=STORE)
store_ax.spines["right"].set_color(STORE)
store_ax.grid(False)
store_ax.annotate(f"drained to {soc_end[6]:,.1f} GWh", xy=(5.15, soc_end[6]), xytext=(6.7, 130),
                  color=STORE, fontsize=8.5, ha="left",
                  arrowprops=dict(arrowstyle="->", color=STORE))
# the trades that produced those profits: volume-weighted average price per month and leg
buy_px  = by_month.apply(lambda d: (d["pump"] * d["price"]).sum() / d["pump"].sum())
sell_px = by_month.apply(lambda d: (d["gen"]  * d["price"]).sum() / d["gen"].sum())
trade_ax.plot(range(12), sell_px.to_numpy(), color=SELL, lw=2.2, marker="o", ms=4.5,
              label="Average sell price")
trade_ax.plot(range(12), buy_px.to_numpy(), color=BUY, lw=2.2, marker="o", ms=4.5,
              label="Average buy price")
trade_ax.annotate(f"${sell_px[6]:,.0f}", xy=(5, sell_px[6]), xytext=(5.55, sell_px[6] - 28),
                  color=SELL, fontsize=9, fontweight="bold", va="top")
trade_ax.annotate(f"${buy_px[7]:,.0f}", xy=(6, buy_px[7]), xytext=(6, buy_px[7] + 30),
                  color=BUY, fontsize=9, fontweight="bold", ha="center")
trade_ax.set_xticks(range(12))
trade_ax.set_xticklabels(months)
trade_ax.set_ylim(0, 560)
trade_ax.set_ylabel("Trade price ($/MWh)")
trade_ax.yaxis.set_major_formatter(dollars)
trade_ax.legend(frameon=False, loc="upper left", fontsize=8.5, ncol=2)
fig.text(0.005, -0.01, "Trade prices are the volume-weighted averages of the month's actual trades.",
         ha="left", va="top", fontsize=8.5, color=INK2, style="italic")
finish(fig, "10c-monthly-profit-2025.png")

The forecaster in action

The two views that closed the perfect-foresight section, repeated for the forecaster’s own dispatch on the held-out year. First, its profit ranked hour by hour. The shape is the ceiling’s shape, slightly exaggerated: half the profit sits in the best ~250 hours (3% of the year, against the ceiling’s 4%), and the curve peaks at 128% before giving back the pumping bill, a give-back two points deeper than the ceiling’s. Both follow from the same fact: the forecaster earns less in total from the same spike hours, so those hours weigh more heavily in what it does earn, and its imperfect timing pays more for the energy it stores.

# The ceiling's concentration figure, repeated for the forecaster: rank held-out 2025's
# half-hours by the profit each contributed and accumulate. The per-half-hour profits are
# rebuilt from the scored dispatch and asserted against the year's headline number.
s3_close = scheme3_by_year[CLOSE_YEAR]
px_close = price_close.to_numpy()
s3_net = ((s3_close["gen"] - s3_close["pump"]) * px_close
          - snowy.var_om_per_mwh * s3_close["gen"]) * snowy.interval_hr
assert abs(s3_net.sum() - s3_close["profit"]) < 1e-3 * abs(s3_close["profit"]), "net dollars drifted"
s3_sorted = np.sort(s3_net)[::-1]
s3_cum = np.cumsum(s3_sorted) / s3_sorted.sum() * 100
s3_hours = np.arange(1, len(s3_sorted) + 1) * 0.5
s3_half = s3_hours[np.searchsorted(s3_cum, 50)]                  # hours that earn half the profit
s3_peak = s3_cum.argmax()                                        # the last hour that still adds money

fig, ax = plt.subplots()
ax.plot(s3_hours, s3_cum, color=RESULT, lw=2.5)                  # purple: the forecaster's result
ax.axhline(100, color=MUTED, lw=1, ls=":")
ax.axhline(50, color=MUTED, lw=1, ls=":")
ax.axvline(s3_half, color=SPIKE, lw=1, ls=":")
ax.fill_between(s3_hours, s3_cum, where=(s3_hours <= s3_half), color=SPIKE, alpha=0.2)
ax.annotate(f"The best ~{s3_half:,.0f} hours\n(~{s3_half / (len(s3_sorted) * 0.5) * 100:.0f}% of the year)\nearn half the profit",
            xy=(s3_half, 50), xytext=(s3_half + 380, 22), color=SPIKE, fontsize=10,
            arrowprops=dict(arrowstyle="->", color=SPIKE))
ax.annotate(f"every money-making hour counted: {s3_cum[s3_peak]:.0f}%",
            xy=(s3_hours[s3_peak], s3_cum[s3_peak]), xytext=(2500, 137),
            color=INK2, fontsize=9, arrowprops=dict(arrowstyle="->", color=MUTED))
ax.annotate("the descent is the pumping bill",
            xy=(7900, float(s3_cum[int(7900 / 0.5) - 1])), xytext=(5500, 62),
            color=BUY, fontsize=9, arrowprops=dict(arrowstyle="->", color=BUY))
ax.set(title=f"Profit Distribution, Forecaster ({CLOSE_YEAR})",
       ylabel="Cumulative share of the year's profit (%)",
       xlim=(0, len(s3_sorted) * 0.5), ylim=(0, 144))
# the x label as one line with reading arrows at each end, in place of a plain set_xlabel
ax.text(0.5, -0.13, "Hours of the year, ranked by profit", transform=ax.transAxes, ha="center", color=INK)
ax.text(0, -0.13, "← most profitable", transform=ax.transAxes, ha="left", fontsize=9, color=INK2)
ax.text(1, -0.13, "least profitable →", transform=ax.transAxes, ha="right", fontsize=9, color=INK2)
finish(fig, "10d-forecaster-profit-concentration.png")

Second, the same winter week the optimal schedule was drawn on. The price line is identical; only the behaviour differs. The forecaster spends 30 June selling its last stored energy into the tail of June’s spikes, touches near-empty (under 5 GWh) on 1 July, and passes the mild days that follow rebuilding: 152 pumping half-hours against 113 generating across the week, with the store climbing to about 70 GWh by the weekend. Perfect foresight, knowing the quiet week ahead, swung the same store between empty and 129 GWh; the forecaster, seeing only odds, rebuilds steadily and keeps selling the evening peaks as they come.

# The perfect-foresight week figure, repeated for the forecaster's dispatch over the same
# seven days (the June to July handover), so the two panels compare like for like. `week`,
# `week_index` and `pos` carry over from that figure's cell.
s3_store_gwh = s3_close["stored"][pos] * snowy.eta_gen / 1e3     # start-of-half-hour level, generatable GWh

fig, price_ax = plt.subplots(figsize=(9, 4.8))
price_ax.fill_between(week_index, 0, 1, where=s3_close["pump"][pos] > 0, step="mid",
                      transform=price_ax.get_xaxis_transform(), color=BUY, alpha=0.15)
price_ax.fill_between(week_index, 0, 1, where=s3_close["gen"][pos] > 0, step="mid",
                      transform=price_ax.get_xaxis_transform(), color=SELL, alpha=0.15)
price_ax.plot(week_index, price_close.loc[week].values, color=INK, lw=1.6, zorder=3)
price_ax.set(title="The Forecaster's Week (30 June to 6 July 2025)", ylabel="Price ($/MWh)")
price_ax.yaxis.set_major_formatter(dollars)
store_ax = price_ax.twinx()                                      # state of charge as a bold green line
store_ax.plot(week_index, s3_store_gwh, color=STORE, lw=2.4)
store_ax.set_ylabel("State of charge (GWh)", color=STORE)
store_ax.set_ylim(bottom=0)
store_ax.tick_params(axis="y", colors=STORE)
store_ax.spines["right"].set_color(STORE)
# in-figure legend, matching the perfect-foresight panel, placed under the plot
price_ax.legend(handles=[plt.Line2D([], [], color=INK, lw=1.6, label="Price"),
                         plt.Line2D([], [], color=STORE, lw=2.4, label="State of charge"),
                         Patch(color=BUY, alpha=0.3, label="Pumping"),
                         Patch(color=SELL, alpha=0.3, label="Generating")],
                frameon=False, loc="upper center", ncol=4, fontsize=8.5, bbox_to_anchor=(0.5, -0.22))
fig.autofmt_xdate()
finish(fig, "10e-forecaster-week.png")

Implied prices and utilisation

Three statistics summarise a whole year of dispatch: the volume-weighted average buy price, the volume-weighted average sell price, and the capacity factor. They compress the backtest into the units any project cash-flow model consumes. The table reports them for the ceiling and for the forecaster, both measured on the held-out year (2025), so the forecaster’s column is not flattered by in-sample knowledge; each scheme’s implied prices come from its own realised trades, and Scheme 3’s column uses the repeated-year boundary disclosed above. The pattern to notice: the plant earns its money on the spread, buying near the bottom of the market and selling at several times its buy price, while running well below full utilisation.

# The three summary statistics for 2025, the held-out test year. Implied prices come from each
# scheme's own realised trades: total dollars paid or received divided by the energy traded.
def implied_prices(outcome, year_prices):
    '''Average realised buy and sell price of a dispatch outcome.'''
    buy  = (outcome["pump"] * year_prices).sum() / outcome["pump"].sum()
    sell = (outcome["gen"]  * year_prices).sum() / outcome["gen"].sum()
    return buy, sell

s3_buy, s3_sell = implied_prices(scheme3_by_year[2025], prices_by_year[2025].to_numpy())
s3_cf = scheme3_by_year[2025]["capacity_factor"]
s1 = ceiling_by_year[2025]
summary_stats = pd.DataFrame({
    "Perfect foresight (S1)": [f"{s1['mean_buy']:.0f}", f"{s1['mean_sell']:.0f}", f"{s1['capacity_factor']:.0%}"],
    "Forecaster (S3)*": [f"{s3_buy:.0f}", f"{s3_sell:.0f}", f"{s3_cf:.0%}"],
}, index=["Mean buy (pump) price ($/MWh)", "Mean sell (generate) price ($/MWh)", "Capacity factor"])
display(summary_stats)

	Perfect foresight (S1)	Forecaster (S3)*
Mean buy (pump) price ($/MWh)	29	33
Mean sell (generate) price ($/MWh)	210	223
Capacity factor	35%	31%

* Both columns are measured on 2025; the forecaster’s rule is out of sample there, trained only on 2022 to 2024.

Interpretation

What we can say. On four years of real prices, an ideally informed Snowy-sized plant could not have made more than about $0.7 to $1.4 bn a year in arbitrage profit (the ceiling), and a realistic operator who knows the market’s pattern but not its future captures most of that ceiling in a normal year, including 84 to 93% in the years after the fuel crisis, with 2025 itself fully held out of the rule’s training. The revenue opportunity is real, large, and mostly capturable without foresight.

What foresight is worth. The realised perfect-foresight gap between the realistic operator and the ceiling is modest in normal years and explodes in 2022. Perfect forecasting is not where the project’s value hides; anticipating rare crises is, and no one can promise that.

What we cannot say. Every scheme still over-states revenue, for one dominant reason and several small ones. The dominant one is Assumption A1: Snowy is not a price taker. At 2,200 MW it pushes evening prices down when it sells and afternoon prices up when it buys, shrinking its own margin, most of all in the extreme hours where the profit concentrates (the concentration chart in Scheme 1). The small ones all lean the same way: no outages (A11), no ramp or start-up limits (A8), no network constraints (A12), no drought (A9). A separate limit is the evidence base itself: four years is a short record, and the central scheme’s edge assumes the future keeps resembling the 2022 to 2024 prices it learned from. A structural shift, from new storage, fresh transmission, or a changed generation mix, could move that pattern out from under it. A first-pass, reduced-form price-impact band in the appendix quantifies the direction: every plausible impact slope takes a material bite out of the ceiling. A serious bid-stack correction remains the natural next step of this work.

From profits to project value. A simulation-grounded risk pass in the appendix pushes the backtest through a standard project cash-flow model. On the held-out 2025 inputs and the official $12 billion capital cost the project returns 4.85% real after tax: comfortably positive at a 3% real discount rate, deeply negative at 7%, with a 7% break-even sell price ($324/MWh) above anything a full-scale operator achieved at comparable utilisation in four years of real prices. Drawing whole years from the backtest and letting capital cost and efficiency vary within defensible bounds does not soften the shape of that verdict: at 7% essentially no draw of ten thousand comes out positive, while at 3% the project is roughly a coin flip, decided by capital-cost discipline and by whether future years look like 2025 or like the 2022 crisis. The discount rate and the capital cost decide this project; the operating strategy defends it.

Interactive companion. The page snowy-arbitrage-demo.html lets the reader vary the two trigger prices of Scheme 2 and watch the profit respond, with the ceiling drawn in for comparison.

Appendix

Glossary

Term	Plain meaning
Arbitrage	Buy low, sell high. Here: store cheap electricity, sell it when it is expensive.
Pumped hydro	A “water battery”. Cheap power pumps water uphill; later it runs back down through turbines to make power.
Spot price	The market-clearing price for electricity, set every five minutes for each NEM region.
Dispatch	The operating schedule: when, and how hard, the plant pumps or generates.
Round-trip efficiency	How much comes back out. Snowy returns about 77% of what it stores; 23% is lost.
State of charge	How full the store is right now, in megawatt-hours (MWh).
Price taker	A player too small to move the price. Snowy is not one.
Perfect foresight	Pretending you knew every future price. Impossible in reality; it gives the ceiling.
MW vs MWh	MW is a rate (how fast); MWh is an amount (how much). 2,200 MW for one hour is 2,200 MWh.
Capacity factor	Average output as a share of maximum. 20% means it generates, on average, a fifth of flat out.
Linear program	An optimisation problem made entirely of straight-line relationships; solvers crack these exactly and fast.
Markov model	A model where the odds of what happens next depend only on the current situation, not the full history.
Backward induction	Solve the last decision first, then step backwards; each earlier decision becomes easy once the future’s value is known.
Water value	The opportunity-cost value of energy held in the reservoir: what a stored MWh should earn if kept instead of sold now.
Perfect-foresight gap	The dollars a perfect forecast would have added in a given year: the ceiling minus the realistic scheme. The realised, one-path cousin of the textbook EVPI (Expected Value of Perfect Information).
Stationary rule	A decision rule with no calendar: the same map from time of day, price and store level to an action, every day. The natural rule for a plant with no final day.
In sample / out of sample	Tested on the data a model learned from (flattering) versus on data it never saw (honest).

Full assumption audit

The Model states up front only the assumptions that move the answer most; this is the complete list. For each: the reason we make it, the direction it pushes the estimate, and whether we test it.

#	Assumption	Why we make it	Effect on the estimate	Tested?
A1	Snowy is a price taker	Standard first pass; keeps every scheme fast and clean	Over-states profit, probably a lot	First pass: the reduced-form band in the appendix
A2	Perfect foresight of prices (S1 only)	Defines the ceiling we measure everything against	Over-states profit	Yes: S3 removes it
A3	Round trip is 77%, split evenly across the two legs	We only know the product; profit depends on the product alone	Neutral	Not needed
A4	Storage is 350 GWh, taken at face value	Snowy’s headline figure, though the cyclable capacity is disputed	Possibly over-states	Yes: shown non-binding in the Analysis
A5	Constant head (the height the water falls, which sets the turbine power)	Head changes little relative to the deep reservoirs	Roughly neutral	Could add head-dependence later
A6	Pump and generate ratings both 2,200 MW	Reported figures; kept symmetric	Minor	Could split
A7	Continuous part-load operation	Snowy uses variable-speed units, so it can modulate	Realistic	Not needed
A8	No ramp, minimum-load or start-up limits	Keeps the ceiling clean	Over-states profit	Pending (an operational model)
A9	Closed water battery: no river inflows or environmental limits	Snowy 2.0 recirculates water; we value only the extra cycling it adds	Slight over-state (drought could cap pumping)	Noted, not modelled
A10	The year repeats (each scheme starts the year at the level it ends at)	Avoids draining the reservoir for free on 31 December, and stops S2 or S3 banking uncredited inventory	Small for S1 (0.01% to 1.9% depending on the year); about 3% for S3, and comparable for S2, where in each case it replaces an arbitrary empty start	Yes: the S3 boundary is reported in the Analysis, and S2 uses the same start = end search
A11	Always available: no planned or forced outages	The model lets the plant run every half hour	Over-states profit	Could apply an availability derate
A12	Clears at the NSW regional price, with no network limits	Keeps the price signal simple	Over-states profit: transmission limits (including HumeLink, the new line being built to carry Snowy’s output), loss factors, and times the plant is constrained off	Not modelled
A13	S3 acts at full power or not at all	The ceiling’s optimal plan already sits at the power limits almost everywhere, so part-load adds little	Negligible	Not needed
A14	The two trained schemes assume 2025 still resembles 2022 to 2024 (S2’s frozen triggers, S3’s learned pattern)	Any trained rule must learn from history; the honest check is a year it never saw	Neutral if the pattern holds; a structural shift would erode both	Yes: 2025 held out entirely for both

Hindsight benchmark for the fixed rule

Scheme 2 freezes its two triggers before 2025. To see how much of its shortfall comes from that frozen pair rather than from the shape of a two-trigger rule, we re-select the triggers with 2025 hindsight over a finer grid built on the 2025 prices themselves, scoring each candidate under the same start = end boundary. This is an ex-post grid-search benchmark over that candidate set, not a formal ceiling over every possible fixed rule.

# The best a fixed two-trigger rule could have done WITH hindsight: search a fine grid of dollar
# triggers built from 2025's OWN price percentiles, plus the frozen deployable pair so the
# benchmark can never fall below Scheme 2. Every candidate is scored under the same start = end
# boundary, and every candidate is asserted to settle to a unique cyclic fixed point. This is a
# benchmark over this candidate set, not a proof-level upper bound over every possible fixed rule.
def best_fixed_rule(year_prices):
    buys  = np.percentile(year_prices, np.arange(5, 46, 5))       # candidate buy triggers (cheap end)
    sells = np.percentile(year_prices, np.arange(55, 96, 5))      # candidate sell triggers (dear end)
    grid  = [(b, s) for b in buys for s in sells if b < s]
    grid.append((buy_trigger, sell_trigger))                      # include the frozen pair itself
    scored = [cyclic_threshold_score(year_prices, snowy, b, s) for b, s in grid]
    assert all(r["converged"] and r["unique_fixed_point"] for r in scored), \
        "some hindsight-grid candidate did not settle to a unique cyclic fixed point"
    return max(r["profit"] for r in scored)

s2_hindsight_2025 = best_fixed_rule(prices_by_year[2025].to_numpy())
ceiling_2025    = ceiling_by_year[2025]["profit"]
forecaster_2025 = scheme3_by_year[2025]["profit"]
depl_gap    = (ceiling_2025 - s2_by_year[2025]) / 1e9             # deployable rule -> ceiling
retune_gain = (s2_hindsight_2025 - s2_by_year[2025]) / 1e9        # gain from reselecting over the finer grid
fc_recovers = (forecaster_2025 - s2_hindsight_2025) / 1e9         # what the forecaster adds on top
resid_gap   = (ceiling_2025 - forecaster_2025) / 1e9             # realised forecaster -> ceiling gap
display(Markdown(
    f"Re-selecting the triggers with 2025 hindsight over this finer grid raises the frozen rule's "
    f"2025 profit from **\\${s2_by_year[2025]/1e9:.2f} bn** to **\\${s2_hindsight_2025/1e9:.2f} bn**, "
    f"a gain of about **\\${retune_gain:.2f} bn**. Most of its **\\${depl_gap:.2f} bn** gap to the "
    f"perfect-foresight ceiling therefore remains even after re-selecting the pair on this grid, which suggests the "
    f"policy *form* matters more than the particular frozen pair. The fitted forecaster adds a "
    f"further **\\${fc_recovers:.2f} bn**, and the remaining **\\${resid_gap:.2f} bn** is the "
    f"realised gap between this forecaster and the ceiling."
))

Re-selecting the triggers with 2025 hindsight over this finer grid raises the frozen rule’s 2025 profit from $0.90 bn to $0.95 bn, a gain of about $0.05 bn. Most of its $0.22 bn gap to the perfect-foresight ceiling therefore remains even after re-selecting the pair on this grid, which suggests the policy form matters more than the particular frozen pair. The fitted forecaster adds a further $0.09 bn, and the remaining $0.08 bn is the realised gap between this forecaster and the ceiling.

A first pass at price impact (relaxing A1)

Everything above assumes the plant can trade without moving the price (A1). At 2,200 MW that is false: selling into the evening peak pushes that peak down, pumping through the cheap midday pulls the trough up, and both squeeze exactly the price gap the plant lives on. This section asks the question that follows: how much of the 2025 profit survives once the plant pays for its own footprint?

The model is deliberately simple. The price the plant actually clears at moves in proportion to its own net injection: adjusted price = spot − β × (generation − pumping). The slope β is the only new parameter, and no public estimate of it exists for NSW (a serious one needs residual-demand or bid-stack data), so rather than pretend to know it we sweep four values and label each by the most tangible number it implies: how far the price moves when the plant trades at full power. The four scenarios move the full-power price by $5.50, $11, $22 and $44/MWh, a range we believe brackets any defensible estimate.

With the quadratic revenue term the perfect-foresight problem stops being a linear program, but it stays concave, so it can still be solved exactly up to a stated tolerance: we replace the quadratic with a fine piecewise-linear outer approximation (33 tangent lines), solve the resulting linear program with the same solver as Scheme 1, and then recompute the schedule’s profit under the exact quadratic prices. The linear program’s objective is an upper bound on the true optimum and the recomputed profit is a lower bound, so the pair brackets the truth; the code asserts the bracket is tighter than half a percent. The re-optimised operator adapts to its own footprint, trading fewer megawatt-hours and spreading them across more moderate hours; for comparison we also price the unadapted schedule (the price-taker’s plan executed under impact pricing), which shows what re-optimisation recovers.

def solve_with_price_impact(prices: np.ndarray, plant: StoragePlant, beta: float,
                            knots: int = 33) -> dict:
    '''Perfect-foresight dispatch when the plant's own net injection moves the price by
    beta $/MWh per MW. The concave quadratic profit is maximised through a piecewise-linear
    outer approximation (lower tangents of net^2), so HiGHS still solves it exactly as an LP;
    the schedule's profit is then recomputed under the exact quadratic prices.'''
    n, dt, P = len(prices), plant.interval_hr, plant.power_mw
    pump   = cp.Variable(n, nonneg=True)
    gen    = cp.Variable(n, nonneg=True)
    stored = cp.Variable(n + 1, nonneg=True)
    w      = cp.Variable(n, nonneg=True)          # stands in for net^2 via tangent cuts
    net = gen - pump
    constraints = [pump <= P, gen <= P, stored <= plant.max_stored_energy,
                   stored[1:] == stored[:-1] + plant.eta_pump * pump * dt - (gen / plant.eta_gen) * dt,
                   stored[0] == stored[n]]
    for qk in np.linspace(-P, P, knots):          # w >= every lower tangent of q^2 at knot qk
        constraints.append(w >= 2 * qk * net - qk**2)
    # A vanishing $0.0001/MWh charge on both legs breaks the pump-and-generate-at-once tie
    # that impact pricing can otherwise create at negative prices; economically it is nil.
    objective = (cp.sum(cp.multiply(prices, net) * dt) - beta * cp.sum(w) * dt
                 - plant.var_om_per_mwh * cp.sum(gen * dt) - 1e-4 * cp.sum(pump + gen) * dt)
    problem = cp.Problem(cp.Maximize(objective), constraints)
    problem.solve(solver=cp.HIGHS)
    g, u = gen.value, pump.value
    q = g - u
    realised = float((prices * q * dt - beta * q**2 * dt - plant.var_om_per_mwh * g * dt).sum())
    upper = float(problem.value + 1e-4 * (u + g).sum() * dt)     # bound without the tie-break charge
    assert problem.status == "optimal", f"price-impact LP ended {problem.status}"
    assert float(np.minimum(u, g).max()) < 0.5, "simultaneous pump and generate in the impact solve"
    assert (upper - realised) / realised < 0.005, "PWL bracket looser than 0.5%; raise the knot count"
    return {"profit": realised, "gen": g, "pump": u, "sold_gwh": g.sum() * dt / 1e3}

# Sweep the slope band on the held-out year, and price the UNADAPTED price-taker schedule
# under the same impact for comparison (what an operator who ignored its footprint would earn).
p25 = prices_by_year[2025].to_numpy()
base25 = ceiling_by_year[2025]
betas = [0.0025, 0.005, 0.01, 0.02]              # $/MWh per MW of net injection
impact_rows = []
for beta in betas:
    adapted = solve_with_price_impact(p25, snowy, beta)
    q0 = base25["gen"] - base25["pump"]          # the price-taker's schedule, unchanged
    naive = float((p25 * q0 - beta * q0**2 - snowy.var_om_per_mwh * base25["gen"]).sum() * snowy.interval_hr)
    assert adapted["profit"] >= naive - 1e-6, "re-optimising under impact cannot do worse"
    impact_rows.append({
        "Slope ($/MWh per 1,000 MW)": f"{beta * 1e3:g}",
        "Price move at full power ($/MWh)": f"{beta * snowy.power_mw:.1f}",
        "Re-optimised profit ($bn)": round(adapted["profit"] / 1e9, 2),
        "Erosion vs ceiling": f"{(1 - adapted['profit'] / base25['profit']) * 100:.0f}%",
        "Unadapted schedule ($bn)": round(naive / 1e9, 2),
        "Energy sold (GWh)": round(adapted["sold_gwh"], 0)})
impact_2025 = pd.DataFrame(impact_rows).set_index("Slope ($/MWh per 1,000 MW)")
impact_profits = [r["Re-optimised profit ($bn)"] for r in impact_rows]
assert all(a > b for a, b in zip(impact_profits, impact_profits[1:])), \
    "profit should fall monotonically in the impact slope"
display(impact_2025)

	Price move at full power ($/MWh)	Re-optimised profit ($bn)	Erosion vs ceiling	Unadapted schedule ($bn)	Energy sold (GWh)
Slope ($/MWh per 1,000 MW)
2.5	5.5	1.05	7%	1.04	6018.0
5	11.0	0.98	13%	0.96	5482.0
10	22.0	0.87	23%	0.79	4635.0
20	44.0	0.70	38%	0.45	3317.0

Read the figure left to right: each bar assumes the plant’s trading moves the price a little more. The grey bar is what the plant still earns after re-planning its whole year around that self-inflicted price move; the faint red block above it, up to the dashed line, is what the price-taker headline over-stated; the short dark dash inside each bar is a naive operator who executes the price-taker plan regardless.

# One bar per impact scenario. Grey = profit after re-optimising around the plant's own
# footprint; faint red block = the bite out of the price-taker headline (dashed line); the
# dark dash = the naive operator (price-taker schedule executed unchanged), so the gap from
# dash to bar top is what re-planning recovers.
ceiling_bn = base25["profit"] / 1e9
labels = ["No impact\n(price taker)"] + [f"${b * snowy.power_mw:g}/MWh" for b in betas]
adapted_bn = [ceiling_bn] + [r["Re-optimised profit ($bn)"] for r in impact_rows]
naive_bn = [r["Unadapted schedule ($bn)"] for r in impact_rows]

fig, ax = plt.subplots(figsize=(9, 4.5))
xs = np.arange(len(labels))
ax.bar(xs, adapted_bn, width=0.62, color=MUTED, alpha=0.9, edgecolor=INK, linewidth=1.1)
ax.bar(xs[1:], [ceiling_bn - a for a in adapted_bn[1:]], bottom=adapted_bn[1:], width=0.62,
       color=SPIKE, alpha=0.15, edgecolor=SPIKE, linewidth=0.8, linestyle=":")
ax.axhline(ceiling_bn, color=INK, lw=1.2, ls="--")
ax.plot(xs[1:], naive_bn, ls="none", marker="_", ms=24, mew=2.6, color=INK)
for x, a in zip(xs[1:], adapted_bn[1:]):                       # the erosion, above each bar
    ax.annotate(f"-{(1 - a / ceiling_bn) * 100:.0f}%", (x, ceiling_bn), xytext=(0, 7),
                textcoords="offset points", ha="center", color=SPIKE, fontsize=10.5, fontweight="bold")
ax.set_xticks(xs)
ax.set_xticklabels(labels, fontsize=9)
ax.set(title="How Much of the 2025 Profit Survives the Plant's Own Price Impact?",
       xlabel="Assumed price move when the plant trades at full power (2,200 MW)",
       ylabel="Annual profit ($bn)", ylim=(0, ceiling_bn * 1.16))
ax.yaxis.set_major_formatter(lambda x, _: f"${x:,.1f}")
ax.legend(handles=[Patch(color=MUTED, alpha=0.9, label="Earned after re-planning"),
                   Patch(facecolor=SPIKE, alpha=0.15, edgecolor=SPIKE, ls=":", label="Lost to own price impact"),
                   plt.Line2D([], [], ls="--", lw=1.2, color=INK, label="Price-taker headline"),
                   plt.Line2D([], [], ls="none", marker="_", ms=14, mew=2.6, color=INK,
                              label="Naive: not re-planned")],
          frameon=False, loc="upper center", bbox_to_anchor=(0.5, -0.22), ncol=2, fontsize=8.5)
finish(fig, "10f-price-impact.png")

_ero_lo = (1 - impact_rows[0]["Re-optimised profit ($bn)"] * 1e9 / base25["profit"]) * 100
_ero_hi = (1 - impact_rows[-1]["Re-optimised profit ($bn)"] * 1e9 / base25["profit"]) * 100
display(Markdown(
    f"Across the band, the ceiling's 2025 profit falls by **{_ero_lo:.0f}% to {_ero_hi:.0f}%**, "
    f"and even the smallest slope costs several times more than perfect foresight was worth that "
    f"year. Three readings. First, the direction and rough size of the A1 bias is now quantified: "
    f"headline profits in this report are over-statements by an amount that plausibly starts near "
    f"a tenth and could exceed a third. Second, re-optimisation matters: an operator who re-plans "
    f"around its own footprint trades less, trades wider, and keeps materially more than one who "
    f"executes the price-taker plan regardless. Third, the same haircut applies in spirit to the "
    f"forecaster, whose dispatch is also predominantly full-power: a cash-flow model consuming "
    f"this report's implied prices should stress the margin per megawatt-hour sold by the same "
    f"fractions. An exact forecaster correction needs partial-power actions in its rule and is "
    f"future work, as is estimating the slope itself from bid-stack data."))

Across the band, the ceiling’s 2025 profit falls by 7% to 38%, and even the smallest slope costs several times more than perfect foresight was worth that year. Three readings. First, the direction and rough size of the A1 bias is now quantified: headline profits in this report are over-statements by an amount that plausibly starts near a tenth and could exceed a third. Second, re-optimisation matters: an operator who re-plans around its own footprint trades less, trades wider, and keeps materially more than one who executes the price-taker plan regardless. Third, the same haircut applies in spirit to the forecaster, whose dispatch is also predominantly full-power: a cash-flow model consuming this report’s implied prices should stress the margin per megawatt-hour sold by the same fractions. An exact forecaster correction needs partial-power actions in its rule and is future work, as is estimating the slope itself from bid-stack data.

Results store

Every render writes the backtest’s outputs to python/results/, so simple questions can be answered without re-running the model: the half-hourly dispatch of all three schemes as one tidy table (dispatch.parquet), the annual summary (summary.csv), the financial-modelling table below (fin_inputs.csv), the frozen parameters (parameters.csv), and an Excel workbook of the small tables (results.xlsx). The companion script workbench.py loads the store for quick queries; the summary numbers are recomputed from the dispatch table itself and checked against the schemes’ reported profits before anything is written.

# Persist the backtest's outputs so they can be queried later (by workbench.py, or straight
# from Excel) without re-running the solvers. Everything below is derived from objects the
# report has already computed; nothing new is decided here.
RESULTS_DIR = DATA_DIR.parent / "results"
RESULTS_DIR.mkdir(parents=True, exist_ok=True)
dt_hr = snowy.interval_hr

# Step 1. One long, tidy table of every scheme's half-hourly dispatch. Scheme 2's simulator
# returns only the reservoir path, so its two legs are reconstructed from the level changes:
# a rise in the store came from pumping, a fall came from generating (exact by construction,
# because the simulator never does both in the same half hour).
def _legs_from_store(stored_path):
    ds = np.diff(stored_path)
    pump = np.where(ds > 0, ds / (snowy.eta_pump * dt_hr), 0.0)
    gen  = np.where(ds < 0, -ds * snowy.eta_gen / dt_hr, 0.0)
    return pump, gen

_frames = []
for _year in YEARS:
    _px = prices_by_year[_year]
    for _scheme, _out in (("S1 perfect foresight", ceiling_by_year[_year]),
                          ("S2 rule of thumb",     s2_scored[_year]),
                          ("S3 forecaster",        scheme3_by_year[_year])):
        _pump, _gen = (_out["pump"], _out["gen"]) if "pump" in _out else _legs_from_store(_out["stored"])
        _frames.append(pd.DataFrame({
            "year": _year, "scheme": _scheme,
            "interval_start": _px.index, "price_per_mwh": _px.to_numpy(),
            "pump_mw": _pump, "gen_mw": _gen,
            "stored_mwh": np.asarray(_out["stored"])[:-1],   # level at the START of each half hour
        }))
dispatch = pd.concat(_frames, ignore_index=True)

# Step 2. The annual summary, recomputed FROM the dispatch table (profit = sales - purchases
# - variable O&M, the report's definition throughout) so the store is self-consistent...
dispatch["profit_dollars"] = ((dispatch["gen_mw"] - dispatch["pump_mw"]) * dispatch["price_per_mwh"]
                              - snowy.var_om_per_mwh * dispatch["gen_mw"]) * dt_hr
summary = (dispatch.groupby(["year", "scheme"])
           .agg(profit_bn=("profit_dollars", lambda s: s.sum() / 1e9),
                energy_sold_gwh=("gen_mw", lambda s: s.sum() * dt_hr / 1e3),
                energy_bought_gwh=("pump_mw", lambda s: s.sum() * dt_hr / 1e3))
           .reset_index())
_g = dispatch.groupby(["year", "scheme"])
summary["capacity_factor"] = _g.apply(
    lambda d: d["gen_mw"].sum() / (snowy.power_mw * len(d)), include_groups=False).to_numpy()
summary["mean_buy_per_mwh"] = _g.apply(
    lambda d: (d["pump_mw"] * d["price_per_mwh"]).sum() / d["pump_mw"].sum(), include_groups=False).to_numpy()
summary["mean_sell_per_mwh"] = _g.apply(
    lambda d: (d["gen_mw"] * d["price_per_mwh"]).sum() / d["gen_mw"].sum(), include_groups=False).to_numpy()
_ceiling_bn = summary.loc[summary["scheme"] == "S1 perfect foresight"].set_index("year")["profit_bn"]
summary["share_of_ceiling"] = summary["profit_bn"] / summary["year"].map(_ceiling_bn)

# ... and checked against the profits the schemes themselves reported, so a store that
# disagrees with the report can never be written.
_reported = {("S1 perfect foresight", y): ceiling_by_year[y]["profit"] for y in YEARS}
_reported |= {("S2 rule of thumb", y): s2_by_year[y] for y in YEARS}
_reported |= {("S3 forecaster", y): scheme3_by_year[y]["profit"] for y in YEARS}
for _row in summary.itertuples():
    _want = _reported[(_row.scheme, _row.year)] / 1e9
    assert abs(_row.profit_bn - _want) <= max(1e-4, 1e-6 * abs(_want)), \
        f"store profit for {_row.scheme} {_row.year} disagrees with the report"

# Step 3. The frozen parameters a downstream user needs to interpret the store.
parameters = pd.DataFrame([
    ("power_mw", snowy.power_mw, "plant power rating, either direction"),
    ("max_stored_energy_mwh", snowy.max_stored_energy, "reservoir capacity in stored potential-energy units; 350 GWh dischargeable (A4)"),
    ("round_trip_efficiency", snowy.eta_pump * snowy.eta_gen, "eta, split evenly across the legs"),
    ("var_om_per_mwh", snowy.var_om_per_mwh, "variable O&M charged on generation"),
    ("interval_hr", dt_hr, "length of one dispatch step, hours"),
    ("s2_buy_trigger_per_mwh", buy_trigger, "frozen on 2022 to 2024 pooled prices"),
    ("s2_sell_trigger_per_mwh", sell_trigger, "frozen on 2022 to 2024 pooled prices"),
    ("train_years", "2022 2023 2024", "S2 and S3 training window; 2025 held out"),
], columns=["parameter", "value", "note"])

# Step 3b. The financial-modelling table: the summary statistics plus the volume-weighted
# spread of the prices each scheme actually traded at. The percentiles weight each half-hour's
# price by the energy moved in it, so they describe the dollars, not the clock.
def _wq(values, weights, qs):
    '''Volume-weighted quantiles: the price below which the given share of energy traded.'''
    order = np.argsort(values)
    v, w = values[order], weights[order]
    cum = np.cumsum(w) - 0.5 * w                    # midpoint convention
    return np.interp(np.asarray(qs) * w.sum(), cum, v)

_fin_rows = []
for (_year, _scheme), _d in dispatch.groupby(["year", "scheme"]):
    _px, _pump, _gen = (_d[c].to_numpy() for c in ("price_per_mwh", "pump_mw", "gen_mw"))
    _b10, _b50, _b90 = _wq(_px, _pump, [0.1, 0.5, 0.9])
    _s10, _s50, _s90 = _wq(_px, _gen,  [0.1, 0.5, 0.9])
    _srow = summary.set_index(["year", "scheme"]).loc[(_year, _scheme)]
    _fin_rows.append({
        "year": _year, "scheme": _scheme,
        "energy_sold_gwh": round(_srow["energy_sold_gwh"], 1),
        "capacity_factor": round(_srow["capacity_factor"], 4),
        "buy_p10": round(_b10, 2), "buy_mean": round(_srow["mean_buy_per_mwh"], 2),
        "buy_p50": round(_b50, 2), "buy_p90": round(_b90, 2),
        "sell_p10": round(_s10, 2), "sell_mean": round(_srow["mean_sell_per_mwh"], 2),
        "sell_p50": round(_s50, 2), "sell_p90": round(_s90, 2),
        "margin_per_mwh_sold": round(_srow["mean_sell_per_mwh"]
                                     - _srow["mean_buy_per_mwh"] / (snowy.eta_pump * snowy.eta_gen)
                                     - snowy.var_om_per_mwh, 2),
        "gross_profit_bn": round(_srow["profit_bn"], 3)})
fin_inputs = pd.DataFrame(_fin_rows)

# Step 4. Write everything: parquet for the big table, CSV for the small ones, and one Excel
# workbook (summary, monthly profit, financial-model inputs, parameters) for spreadsheet users.
dispatch.to_parquet(RESULTS_DIR / "dispatch.parquet", index=False)
summary.to_csv(RESULTS_DIR / "summary.csv", index=False)
parameters.to_csv(RESULTS_DIR / "parameters.csv", index=False)
fin_inputs.to_csv(RESULTS_DIR / "fin_inputs.csv", index=False)
monthly = (dispatch.assign(month=dispatch["interval_start"].dt.to_period("M").astype(str))
           .groupby(["scheme", "month"], as_index=False)["profit_dollars"].sum()
           .assign(profit_m=lambda d: d["profit_dollars"] / 1e6)
           .pivot(index="month", columns="scheme", values="profit_m").round(2))
with pd.ExcelWriter(RESULTS_DIR / "results.xlsx", engine="openpyxl") as _xl:
    summary.round(4).to_excel(_xl, sheet_name="Summary", index=False)
    monthly.to_excel(_xl, sheet_name="Monthly profit ($m)")
    fin_inputs.to_excel(_xl, sheet_name="Financial model inputs", index=False)
    parameters.to_excel(_xl, sheet_name="Parameters", index=False)
print(f"Results store written to {RESULTS_DIR}/: dispatch.parquet "
      f"({len(dispatch):,} rows), summary.csv, fin_inputs.csv, parameters.csv, results.xlsx")

Results store written to results/: dispatch.parquet (210,384 rows), summary.csv, fin_inputs.csv, parameters.csv, results.xlsx

Annual results for financial modelling

Any project cash-flow model consumes a year of dispatch through a handful of numbers: how much energy was sold, at what average prices each leg traded (pumping = buying, generating = selling), and how hard the plant ran. The table gives those numbers for every scheme and year, plus the volume-weighted spread of the realised prices (the price below which 10%, 50% and 90% of each leg’s energy traded). The averages are the point inputs; the percentile columns are what a point input hides, and they are the honest way to parameterise a sensitivity or risk analysis rather than guessing symmetric bands. The margin column grosses the buy leg up by the round-trip loss and subtracts variable operating cost, so margin times energy sold reproduces each year’s gross profit.

_display_names = {"S1 perfect foresight": "Perfect foresight (S1)",
                  "S2 rule of thumb": "Rule of thumb (S2)",
                  "S3 forecaster": "Forecaster (S3)"}
fin_view = (fin_inputs
            .assign(scheme=[_display_names[s] + ("*" if y == 2025 and s != "S1 perfect foresight" else "")
                            for s, y in zip(fin_inputs["scheme"], fin_inputs["year"])])
            .set_index(["year", "scheme"]).rename_axis(["Year", "Scheme"]))
fin_view.columns = pd.MultiIndex.from_tuples([
    ("Energy", "sold GWh"), ("Energy", "capacity factor"),
    ("Buy $/MWh", "P10"), ("Buy $/MWh", "mean"), ("Buy $/MWh", "P50"), ("Buy $/MWh", "P90"),
    ("Sell $/MWh", "P10"), ("Sell $/MWh", "mean"), ("Sell $/MWh", "P50"), ("Sell $/MWh", "P90"),
    ("Result", "margin $/MWh sold"), ("Result", "gross profit $bn")])
(fin_view.style
 .format({("Energy", "sold GWh"): "{:,.0f}", ("Energy", "capacity factor"): "{:.0%}",
          ("Buy $/MWh", "P10"): "{:,.0f}", ("Buy $/MWh", "mean"): "{:,.0f}",
          ("Buy $/MWh", "P50"): "{:,.0f}", ("Buy $/MWh", "P90"): "{:,.0f}",
          ("Sell $/MWh", "P10"): "{:,.0f}", ("Sell $/MWh", "mean"): "{:,.0f}",
          ("Sell $/MWh", "P50"): "{:,.0f}", ("Sell $/MWh", "P90"): "{:,.0f}",
          ("Result", "margin $/MWh sold"): "{:,.0f}", ("Result", "gross profit $bn"): "{:.2f}"})
 # centre the merged Year cells vertically and horizontally so each year reads as one block
 .set_table_styles([
     {"selector": "th.row_heading.level0", "props": [("vertical-align", "middle"), ("text-align", "center")]},
     {"selector": "th.row_heading.level1", "props": [("text-align", "left")]},
     {"selector": "th.col_heading", "props": [("text-align", "center")]},
     {"selector": "td", "props": [("text-align", "right")]}]))

		Energy		Buy $/MWh				Sell $/MWh				Result
		sold GWh	capacity factor	P10	mean	P50	P90	P10	mean	P50	P90	margin $/MWh sold	gross profit $bn
Year	Scheme
2022	Perfect foresight (S1)	5,147	27%	9	101	83	247	119	308	241	508	172	0.89
	Rule of thumb (S2)	1,548	8%	-38	23	31	65	138	207	173	299	173	0.27
	Forecaster (S3)	3,410	18%	-7	64	70	127	118	224	190	300	137	0.47
2023	Perfect foresight (S1)	5,745	30%	-22	39	49	81	92	179	133	265	124	0.71
	Rule of thumb (S2)	2,375	12%	-3	38	49	65	138	238	171	297	184	0.44
	Forecaster (S3)	4,839	25%	-1	49	57	86	106	192	139	277	125	0.60
2024	Perfect foresight (S1)	6,529	34%	-23	40	51	84	93	276	158	293	220	1.43
	Rule of thumb (S2)	3,017	16%	-30	32	46	64	140	382	185	300	337	1.02
	Forecaster (S3)	5,516	29%	-16	45	54	85	106	301	164	294	239	1.32
2025	Perfect foresight (S1)	6,680	35%	-19	29	35	73	102	210	140	242	169	1.13
	Rule of thumb (S2)*	3,587	19%	-18	25	31	63	138	289	167	279	252	0.90
	Forecaster (S3)*	5,932	31%	-18	33	40	78	107	223	146	250	176	1.04

* Out of sample: rule fixed on 2022 to 2024 before trading 2025.

The same information as a picture, for the deployable forecaster: the bar is the P10 to P90 band of the prices its energy actually traded at, the tick is the volume-weighted mean. The two legs live in different parts of the price distribution, which is the entire business, and the generate band is wide because a large share of the revenue rides on a small number of expensive hours. Note 2024: a few extreme half-hours carried so much of the sold volume that the mean sits above the P90 of the band. That is not an error; it is what spike-driven revenue looks like, and it is why the mean and the median tell different stories in a spike year.

# P10-P90 volume-weighted bands per year for the forecaster's two legs, mean marked.
# BUY and SELL keep their fixed meanings (blue = pumping/buying, orange = generating/selling).
s3fin = fin_inputs[fin_inputs["scheme"] == "S3 forecaster"].set_index("year")
fig, ax = plt.subplots(figsize=(9, 4.4))
for i, year in enumerate(YEARS):
    r = s3fin.loc[year]
    for leg, colour, off in (("buy", BUY, -0.16), ("sell", SELL, +0.16)):
        ax.plot([r[f"{leg}_p10"], r[f"{leg}_p90"]], [i + off, i + off],
                color=colour, lw=7, solid_capstyle="butt", alpha=0.45, zorder=1)
        ax.plot(r[f"{leg}_mean"], i + off, marker="|", color=colour, ms=16, mew=3, zorder=2)
ax.set_yticks(range(len(YEARS)))
ax.set_yticklabels([f"{y}*" if y not in TRAIN_YEARS else str(y) for y in YEARS])
ax.invert_yaxis()
ax.set(title="Trading Prices", xlabel="Realised price ($/MWh)")
ax.xaxis.set_major_formatter(dollars)
ax.set_xlim(right=395)                                        # room for the 2024 annotation
# 2024's mean sits above its own P90: a handful of extreme half-hours carry that much volume
r24 = s3fin.loc[2024]
ax.annotate("2024: spike half-hours pull the\nmean above the band's P90",
            xy=(r24["sell_mean"], 2 + 0.16), xytext=(r24["sell_mean"] + 18, 1.05),
            color=INK2, fontsize=8.5, arrowprops=dict(arrowstyle="->", color=MUTED))
ax.legend(handles=[plt.Line2D([], [], color=BUY, lw=7, alpha=0.45, label="Pump price, P10 to P90 of volume"),
                   plt.Line2D([], [], color=SELL, lw=7, alpha=0.45, label="Generate price, P10 to P90 of volume"),
                   plt.Line2D([], [], marker="|", color=INK2, ls="none", ms=12, mew=3, label="Volume-weighted mean")],
          frameon=False, loc="upper center", bbox_to_anchor=(0.5, -0.16), ncol=3, fontsize=8.5)
fig.text(0.985, 0.9, "* out of sample", ha="right", fontsize=8, color=INK2, style="italic")
finish(fig, "11-financial-price-bands.png")

Finally, every scheme-year on one plane. A cash-flow model that takes a buy price, a sell price and a capacity factor can land anywhere on this chart; the twelve points are where four years of real prices actually put the three schemes. The dashed lines are constant gross margin per MWh sold (sell price minus the loss-adjusted buy price minus variable cost), so a point’s line tells you what each sold megawatt-hour cleared before fixed costs.

RESULT = "#7b52ab"                          # the forecaster keeps its purple from the comparison chart
eta_rt = snowy.eta_pump * snowy.eta_gen
fig, ax = plt.subplots(figsize=(9, 5.2))
bx = np.linspace(0, 115, 2)
for m in (50, 100, 150, 200, 250):          # iso-margin lines: sell = margin + buy/eta + var O&M
    ax.plot(bx, m + bx / eta_rt + snowy.var_om_per_mwh, ls="--", color=GRID, lw=1.2, zorder=1)
    ax.annotate(f"${m}/MWh margin", (1.5, m + 1.5 / eta_rt + snowy.var_om_per_mwh),
                color=MUTED, fontsize=8.5, ha="left", va="bottom")
_marks = {"S1 perfect foresight": dict(marker="o", mfc="white", mec=INK, mew=1.8, ms=9, color=INK),
          "S2 rule of thumb":     dict(marker="s", color=MUTED, ms=8),
          "S3 forecaster":        dict(marker="o", color=RESULT, ms=9)}
# default label offset, with hand nudges for the crowded lower-left cluster
_nudges = {("S1 perfect foresight", 2025): (-36, -3), ("S2 rule of thumb", 2022): (-10, -14)}
for scheme, mk in _marks.items():
    d = fin_inputs[fin_inputs["scheme"] == scheme]
    ax.plot(d["buy_mean"], d["sell_mean"], ls="none", zorder=3, **mk)
    for _, r in d.iterrows():
        held_out = r["year"] == 2025 and scheme != "S1 perfect foresight"
        ax.annotate(f"{int(r['year'])}{'*' if held_out else ''}", (r["buy_mean"], r["sell_mean"]),
                    xytext=_nudges.get((scheme, r["year"]), (7, 5)),
                    textcoords="offset points", fontsize=8, color=INK2)
ax.set(title="Every Scheme-Year as a Buy Price and a Sell Price",
       xlabel="Volume-weighted mean buy (pump) price ($/MWh)",
       ylabel="Volume-weighted mean sell (generate) price ($/MWh)", xlim=(0, 115),
       ylim=(40, fin_inputs["sell_mean"].max() + 35))
ax.xaxis.set_major_formatter(dollars)
ax.yaxis.set_major_formatter(dollars)
ax.legend(handles=[plt.Line2D([], [], ls="none", label=_display_names[s], **mk)
                   for s, mk in _marks.items()], frameon=False, loc="lower right")
fig.text(0.985, 0.02, "* out of sample: rule fixed on 2022 to 2024", ha="right",
         fontsize=8, color=INK2, style="italic")
finish(fig, "12-buy-sell-plane.png")

A simulation-grounded risk pass

The backtest compresses each year of dispatch into three numbers. A project decision needs those numbers pushed through a cash-flow model, against capital cost, operating cost, tax, depreciation and a discount rate, and it needs the risk question asked properly: not “what is the NPV of the base case?” but “what does the NPV distribution look like across everything the backtest saw?” This section does both, in the open, with every parameter stated.

The cash-flow model is a standard single-asset DCF, in real dollars throughout. The plant builds for nine years (20% of the capital in year 1, 10% in each of the next eight) and then operates for the remainder of a 150-year design life: generation equals capacity factor × 2,200 MW × 8,760 hours; revenue is the sell price on that generation; costs are variable O&M at $4/MWh, fixed O&M at 1% of the capital cost per year, and the pumping energy bought at the buy price, grossed up by the 77% round trip. Depreciation is straight-line over 40 years, company tax is 30%, and free cash flow is discounted at a real, after-tax WACC. The capital cost is Snowy 2.0’s official $12 billion. The long horizon flatters the project less than it seems: cutting the operating life to 60 years lowers the 3% NPV by about a third and leaves the 7% verdict unchanged.

# A standard single-asset DCF, small enough to state completely. All dollars are real; the
# discount rate is a real, after-tax WACC. Everything the model uses is set right here.
POWER_MW, HOURS_YR = snowy.power_mw, 8760.0
VAR_OM, TAX, FIXED_FRAC = snowy.var_om_per_mwh, 0.30, 0.01
BUILD_SHARE = np.array([0.20] + [0.10] * 8)                  # nine construction years
LIFE_YEARS, DEPN_YEARS = 150, 40
_years = np.arange(1, LIFE_YEARS + 1)
_operating = _years >= 10                                    # operations start in year 10
_depreciating = _operating & (_years - 9 <= DEPN_YEARS)      # 40 years of straight-line depreciation

def project_npv(buy, sell, cf, capex, wacc, rte=0.77, fcas=0.0):
    '''NPV of the plant as a project: build for nine years, then operate every year at the
    given buy and sell prices and capacity factor. Returns the NPV in dollars.'''
    gen_mwh = cf * POWER_MW * HOURS_YR                       # energy sold per operating year
    ebitda = sell * gen_mwh + fcas - (VAR_OM * gen_mwh       # revenue less variable O&M,
                                      + FIXED_FRAC * capex   # fixed O&M (1% of capex a year),
                                      + buy * gen_mwh / rte) # and pumping energy, grossed up for losses
    depn = np.where(_depreciating, capex / DEPN_YEARS, 0.0)
    fcf = np.where(_operating, (ebitda - depn) * (1 - TAX) + depn, 0.0)
    fcf[:9] = -capex * BUILD_SHARE                           # the build years pay out the capital
    return (fcf / (1 + wacc) ** (_years - 1)).sum()

BASE = dict(buy=s3fin.loc[2025, "buy_mean"], sell=s3fin.loc[2025, "sell_mean"],
            cf=s3fin.loc[2025, "capacity_factor"], capex=12e9)
npv3, npv7 = project_npv(**BASE, wacc=0.03), project_npv(**BASE, wacc=0.07)
lo, hi = 0.0, 0.20                                           # IRR by bisection (NPV falls as the rate rises)
for _ in range(60):
    mid = (lo + hi) / 2
    lo, hi = (mid, hi) if project_npv(**BASE, wacc=mid) > 0 else (lo, mid)
irr = lo

def breakeven_sell(wacc):
    '''The sell price at which NPV = 0. NPV is linear in the sell price, so one step suffices.'''
    annuity = (1 / (1 + wacc) ** (_years[_operating] - 1)).sum()
    gen_mwh = BASE["cf"] * POWER_MW * HOURS_YR
    return BASE["sell"] - project_npv(**BASE, wacc=wacc) / (gen_mwh * (1 - TAX) * annuity)

be3, be7 = breakeven_sell(0.03), breakeven_sell(0.07)
# the highest sell price achieved at a comparable utilisation, across every scheme and year
best_at_volume = fin_inputs.loc[fin_inputs["capacity_factor"] >= 0.25, "sell_mean"].max()
assert abs(npv3 - 7.56e9) < 5e7 and abs(npv7 + 3.49e9) < 5e7, "the cash-flow model lost its anchors"
assert 0.048 < irr < 0.049, "the IRR moved off its anchor"
display(Markdown(
    f"Fed the forecaster's held-out 2025 numbers (buy **\\${BASE['buy']:.2f}/MWh**, sell "
    f"**\\${BASE['sell']:.2f}/MWh**, capacity factor **{BASE['cf']:.1%}**) at the official \\$12 bn "
    f"capital cost, the project returns **{irr*100:.2f}% real after tax**: NPV "
    f"**{npv3/1e9:+.2f} bn at a 3% real WACC**, **{npv7/1e9:+.2f} bn at 7%**. The break-even sell "
    f"price is **\\${be3:.0f}/MWh at 3%** and **\\${be7:.0f}/MWh at 7%**. For scale: across every "
    f"scheme and year in the backtest, the highest volume-weighted sell price achieved at a "
    f"comparable capacity factor is **\\${best_at_volume:.0f}/MWh**, so the 7% hurdle sits above "
    f"anything four years of real prices ever paid a full-scale operator at that utilisation."
))

Fed the forecaster’s held-out 2025 numbers (buy $33.21/MWh, sell $223.02/MWh, capacity factor 30.8%) at the official $12 bn capital cost, the project returns 4.85% real after tax: NPV +7.56 bn at a 3% real WACC, -3.49 bn at 7%. The break-even sell price is $153/MWh at 3% and $324/MWh at 7%. For scale: across every scheme and year in the backtest, the highest volume-weighted sell price achieved at a comparable capacity factor is $308/MWh, so the 7% hurdle sits above anything four years of real prices ever paid a full-scale operator at that utilisation.

The distributions come from the backtest, not from guesses. Three rules govern the Monte Carlo:

1. Prices and utilisation are drawn as whole years. One draw picks one of the forecaster’s four backtested years and takes its buy price, sell price and capacity factor together. This preserves the co-movement the backtest actually found: cheap-pumping years are high-cycling years (the two are perfectly rank-correlated across our four years), while the sell price moves roughly independently of both. Drawing each input separately would invent physically incoherent years, such as crisis-year prices at boom-year utilisation. It also refuses to invent a distribution shape: with four observed years, the honest distribution is those four years, equally weighted.
2. Engineering and cost inputs get Pert distributions. Round-trip efficiency: Pert(0.67, 0.77, 0.80), the published figure as the mode, tunnel-friction pessimism below, best-case operation above. Capital cost: Pert($12 bn, $15 bn, $22 bn), with the official estimate as the floor (it is already contracted spend), a mode above it because this project has already moved from $2 bn announced through $5.9 bn at approval to $12 bn, and a tail acknowledging that large hydro overruns of this scale are historically common.
3. The WACC is deliberately not drawn. The discount rate is the price of risk, not a project risk: mixing it into the same Monte Carlo muddles “how uncertain are the cash flows?” with “what should uncertain cash flows cost?”, and produces a distribution no one can interpret. The whole simulation is priced twice instead, at 3% and at 7% real.

# The recipe above, executed: a year-draw for (buy, sell, capacity factor), independent Perts
# for round-trip efficiency and capital cost, and the whole distribution priced at two WACCs.
rng = np.random.default_rng(20260716)
N_DRAWS = 10_000
s3_triplets = (fin_inputs.query("scheme == 'S3 forecaster'")
               [["buy_mean", "sell_mean", "capacity_factor"]].to_numpy())
year_pick = rng.integers(0, len(s3_triplets), N_DRAWS)       # rows are 2022, 2023, 2024, 2025
buy_d, sell_d, cf_d = s3_triplets[year_pick].T

def pert(a, m, b, size):
    '''A Pert draw: smooth, bounded on [a, b], weighted toward the mode m.'''
    alpha, beta = 1 + 4 * (m - a) / (b - a), 1 + 4 * (b - m) / (b - a)
    return a + rng.beta(alpha, beta, size) * (b - a)

rte_d = pert(0.67, 0.77, 0.80, N_DRAWS)
capex_d = pert(12e9, 15e9, 22e9, N_DRAWS)

gen_d = cf_d * POWER_MW * HOURS_YR
ebitda_d = sell_d * gen_d - (VAR_OM * gen_d + FIXED_FRAC * capex_d + buy_d * gen_d / rte_d)
depn_d = capex_d / DEPN_YEARS
npv_by_wacc = {}
for wacc in (0.03, 0.07):
    disc = 1 / (1 + wacc) ** (_years - 1)
    pv_build = (BUILD_SHARE * disc[:9]).sum()
    a_depn = disc[_depreciating].sum()                       # operating years still depreciating
    a_post = disc[_operating & ~_depreciating].sum()         # operating years after depreciation ends
    npv_by_wacc[wacc] = (-capex_d * pv_build
                         + ((ebitda_d - depn_d) * (1 - TAX) + depn_d) * a_depn
                         + ebitda_d * (1 - TAX) * a_post)

RESULT = "#7b52ab"                                           # purple: the forecaster-based result
fig, ax = plt.subplots(figsize=(9, 4.4))
bins = np.linspace(min(v.min() for v in npv_by_wacc.values()) / 1e9,
                   max(v.max() for v in npv_by_wacc.values()) / 1e9, 70)
for wacc, colour, label in ((0.03, RESULT, "Discounted at 3% real"),
                            (0.07, MUTED, "Discounted at 7% real")):
    ax.hist(npv_by_wacc[wacc] / 1e9, bins=bins, color=colour, alpha=0.55, label=label)
ax.axvline(0, color=INK, lw=1.4, ls="--")
ax.annotate("NPV = 0", (0, ax.get_ylim()[1] * 0.97), xytext=(5, -6), textcoords="offset points",
            color=INK2, fontsize=8.5)
ax.set(title="Ten Thousand Futures, Priced Twice",
       xlabel="Project NPV ($bn)", ylabel=f"Draws (of {N_DRAWS:,})")
ax.xaxis.set_major_formatter(lambda x, _: f"${x:,.0f}")
ax.legend(frameon=False, loc="upper left")
finish(fig, "12b-npv-distribution.png")

q = {w: np.percentile(v, [5, 50, 95]) / 1e9 for w, v in npv_by_wacc.items()}
pos = {w: (v > 0).mean() for w, v in npv_by_wacc.items()}
assert pos[0.07] < 0.01, "the 7% verdict is supposed to be decisive"
# One robustness stat: the same draws with the 2022 crisis year excluded from the pool.
calm = year_pick != 0                                        # triplet row 0 is 2022
med_calm = np.percentile(npv_by_wacc[0.03][calm], 50) / 1e9
pos_calm = (npv_by_wacc[0.03][calm] > 0).mean()
annuity7 = (1 / 1.07 ** (_years[_operating] - 1)).sum()
fcas_needed = -npv7 / ((1 - TAX) * annuity7)                 # flat extra revenue closing the 7% gap
seven_word = "not one" if pos[0.07] == 0 else f"only {pos[0.07]:.1%}"
display(Markdown(
    f"At **7% real**, the verdict is not just negative at the base case: **{seven_word} of "
    f"{N_DRAWS:,} draws comes out positive** (P5 {q[0.07][0]:+.1f}, median {q[0.07][1]:+.1f}, "
    f"P95 {q[0.07][2]:+.1f} \\$bn). No plausible side revenue changes that: closing the base-case "
    f"gap would take about **\\${fcas_needed/1e6:,.0f}M of extra revenue every operating year**, "
    f"roughly twice the size of the entire NEM frequency-control (FCAS) market in a typical "
    f"recent year. At **3% real** the project is close to a coin flip: **{pos[0.03]:.0%} of draws "
    f"are positive** (P5 {q[0.03][0]:+.1f}, median {q[0.03][1]:+.1f}, P95 {q[0.03][2]:+.1f} \\$bn), "
    f"and the coin is weighted by the two inputs the operator does not control: whether the "
    f"capital cost holds near \\$12 bn, and whether future years look like 2025 or like the 2022 "
    f"crisis (excluding 2022 from the draw pool lifts the 3% median to {med_calm:+.1f} \\$bn and "
    f"the positive share to {pos_calm:.0%}). The conclusion to carry forward: **the discount rate "
    f"and the capital cost decide this project; the operating strategy defends it.** Any "
    f"Monte Carlo built on this same recipe should reproduce these percentiles to within sampling "
    f"noise; treat any material disagreement as a bug, not a difference of opinion."
))

At 7% real, the verdict is not just negative at the base case: not one of 10,000 draws comes out positive (P5 -11.5, median -7.6, P95 -3.5 $bn). No plausible side revenue changes that: closing the base-case gap would take about $600M of extra revenue every operating year, roughly twice the size of the entire NEM frequency-control (FCAS) market in a typical recent year. At 3% real the project is close to a coin flip: 49% of draws are positive (P5 -8.0, median -0.8, P95 +10.6 $bn), and the coin is weighted by the two inputs the operator does not control: whether the capital cost holds near $12 bn, and whether future years look like 2025 or like the 2022 crisis (excluding 2022 from the draw pool lifts the 3% median to +4.0 $bn and the positive share to 65%). The conclusion to carry forward: the discount rate and the capital cost decide this project; the operating strategy defends it. Any Monte Carlo built on this same recipe should reproduce these percentiles to within sampling noise; treat any material disagreement as a bug, not a difference of opinion.

Data notes

Prices are the AEMO dispatch price for the NSW1 region (the regional reference price), downloaded from AEMO’s public NEMWeb archive with the open-source NEMOSIS package. Rows flagged as intervention pricing re-runs are dropped (we keep INTERVENTION = 0), each 5-minute price is labelled by the start of its interval, and the series is averaged to 30-minute steps. The first run of the loader saves one small snapshot file per year under python/data/; those snapshots are committed to the repository, so every later run (and every reader) reproduces identical numbers without touching the network. Delete a snapshot to force a fresh download.

Every figure in this report is also saved as a PNG under python/figures/ on each run, ready for reuse elsewhere.

References

• Snowy Hydro, About Snowy 2.0 (accessed July 2026): capacity, storage, and project description.
• Wikipedia, Snowy 2.0 Pumped Storage Power Station (accessed July 2026): cost and timeline history with sources.
• The Conversation (2019), Snowy 2.0 will not produce nearly as much electricity as claimed: the case that cyclable storage is far below 350 GWh.
• SBS News (2023), “Bloody hard to build”: the Australian mega project at the centre of an energy fight: cost debate and responses to the storage critique.
• AEMO NEMWeb public data archive, accessed through NEMOSIS (UNSW Collaboration on Energy and Environmental Markets).