05. Spark Spreads - CAISO Gas Generator Economics#
Authors: Brian Nguyen · Claire Kim ·
Created: 2026-07-28 ·
Tags: #Experiment #Pricing #MIHR #Spark-Spread
2. Objective & Input Framework#
Objective: measure what the daily average misses. The path:
Build daily implied heat rates and spark spreads across the CAISO hub APNodes (aggregated pricing nodes), each paired with the fuel region supplying its gas.
Measure profitability for an efficient combined-cycle gas turbine (CCGT, HR 7.0) and a gas peaker (HR 10.5).
Quantify regional gas basis between the matched fuel regions.
Price the hourly HRCO strip across strikes, and compare hourly against daily settlement to measure the intraday optionality premium \(\Pi(H)\).
Examine how that premium varies by season.
Measure how much of the strip value concentrates in the extreme tail of hours.
Pricing layer: every index, payoff and settlement value below comes
from the pricing and markets packages.
Data & Feature Catalog:
Dataset |
Granularity |
Role |
|---|---|---|
|
hourly (hub APNode) |
day-ahead locational marginal prices (LMPs) - the power leg of the spark spread |
|
hourly (fuel region) |
natural gas price indices ($/MMBtu) from the OASIS commitment costs - the gas leg of the spark spread |
|
snapshot (fuel region) |
lookup table mapping each fuel region to its pricing hub and balancing authority - pairs the power and gas legs |
|
snapshot (generating unit) |
EIA-860M generator inventory locating gas-fired units on the 4.1 map |
|
snapshot (gas utility) |
CEC service-area polygons for California’s gas utilities - the 4.1 map’s territory layer |
3. Mathematical Model & Analytical Formulation#
The model compares electricity revenue with the cost of gas needed to generate one MWh. The power price is hourly, while the gas price is a daily index carried across that day’s hours (verified in section 4.2).
3.1 Spark Spread and Implied Heat Rate#
For a heat rate \(H\), the hourly spark spread is
where \(P_t\) is the hourly power price and \(G_d\) is the daily gas price. The market-implied heat rate (MIHR) is the power price divided by the gas price:
It is the heat rate at which generation breaks even in that hour.
3.2 Daily-Average Benchmark#
The hourly spread is compared with a daily-average version:
where \(\bar P_d\) is the unweighted mean of the day’s hourly power prices; only the electricity price is averaged.
3.3 Profitability and Dispatch Value#
The strike \(H\) is a fixed heat-rate assumption - 7.0 for the CCGT and 10.5 for the peaker, swept from 3.5 to 12.0 in 6.2. Three measures describe profitability:
Always-on margin - the average spark spread across all hours.
ITM share - the share of hours or days with a positive spread.
Dispatch value - the average payoff from operating only when the spread is positive:
\(V(H)\) is the expected payoff of an hourly HRCO strip struck at \(H\).
3.5 Gas Basis#
Gas basis measures a region’s gas price relative to PG&E Citygate:
Positive values are a premium to PG&E Citygate, negative a discount. The market registry defines PG&E Citygate as the reference.
3.6 Tail Concentration of the Strip#
Sort the hourly payoffs \(\max(S_t(H),\ 0)\) largest first, so that \(S^+_{(1)} \geq S^+_{(2)} \geq \cdots \geq S^+_{(N)}\) over the \(N\) hours of the sample. Then
is the share of total strip value delivered by the top \(q\%\) of hours; \(C(q)\) far above \(q\) marks a tail-risk instrument. Section 6.5 computes it per hub.
3.7 Variable Definitions#
Symbol |
Meaning |
Units |
Code |
|---|---|---|---|
\(P_t\), \(\bar P_d\) |
hourly and daily-average power price |
$/MWh |
|
\(G_d\) |
daily gas price |
$/MMBtu |
|
\(H\) |
strike heat rate |
MMBtu/MWh |
|
\(\text{MIHR}_t\) |
market-implied heat rate |
MMBtu/MWh |
|
\(S_t(H)\), \(S_d(H)\) |
spark spread at strike \(H\) |
$/MWh |
|
\(V(H)\) |
hourly dispatch value |
$/MWh |
|
\(V_{\text{hourly}}\), \(V_{\text{daily}}\) |
expected HRCO value under each settlement convention |
$/MWh |
|
\(\Pi(H)\) |
intraday optionality premium |
$/MWh |
|
\(\text{Basis}_r\) |
regional gas basis |
$/MMBtu |
|
\(C(q)\) |
share of strip value from the top \(q\%\) of hours |
% |
|
Setup: Inputs, Source Data, and Hub-to-Gas Pairing#
Set the strike heat rates, open the CAISO market, and resolve each power hub to its matched fuel region.
from dataclasses import asdict
from datetime import date
import polars as pl
from bokeh.io import output_notebook, show
from bokeh.models import (
BoxAnnotation,
ColumnDataSource,
FactorRange,
HoverTool,
Label,
Span,
)
from bokeh.plotting import figure
from bokeh.resources import INLINE
from bokeh.transform import dodge
from chart_theme import (
INK,
SEASON_COLORS,
SECONDARY,
SERIES,
corr_heatmap,
html_table,
season_bands,
style,
)
from geo_utils import merc_x, merc_y, to_multipolygon
from markets.definitions.caiso import CAISO_CLOCK
from markets import MarketContext
from markets.adapters.caiso.gas_leg_resolver import (
CAISO_GAS_PAIRING_RULES,
audit_market_pairing,
resolve_gas_leg_candidates,
resolve_gas_legs,
)
from markets.definitions.caiso import CAISO, WEIM_AVAILABLE_FROM
from pricing.core.parameters import scenario
from pricing.core.schedule import SettlementFrequency
from pricing.engines.historical import HistoricalSettlementEngine
from pricing.engines.strikes import AtTheMoneyStrikeRule
from pricing.evaluation.concentration import value_concentration
from pricing.evaluation.sweeps import ScenarioGrid
from pricing.indices.heat_rate import SparkSpreadIndex
from settings import config
DATA_DIR = config("DATA_DIR")
OS_DIR = DATA_DIR / "opensource"
output_notebook(INLINE, hide_banner=True)
Declare: The Two Reference Units#
The analysis uses two fixed reference heat rates:
CCGT_HR- CCGT strike at 7.0.PEAKER_HR- gas peaker strike at 10.5.PHYSICAL_FLOOR- theoretical efficiency limit at 3.412.ASSET_STRIKES- the two strikes as one scenario grid.
PHYSICAL_FLOOR = 3.412
CCGT_HR = 7.0
PEAKER_HR = 10.5
ASSET_HEAT_RATES = {
f"Efficient CCGT - HR {CCGT_HR} (hypothetical)": CCGT_HR,
f"Gas peaker - HR {PEAKER_HR} (hypothetical)": PEAKER_HR,
}
UNIT_LABELS = {CCGT_HR: "CCGT", PEAKER_HR: "Peaker"}
ASSET_STRIKES = ScenarioGrid(strike_heat_rate=list(ASSET_HEAT_RATES.values()))
Declare: Which Gas Market Each Trading Hub Buys In#
Each trading hub is assigned to one CAISO fuel region by
markets.adapters.caiso.gas_leg_resolver, using three fields from the
reference table:
balancing authority- the balancing authority of the trading hub:CISO,PACE, orPACW.serving gas utility- the gas utility serving the hub.region id prefix- restricts candidates to the utility’s service territory; null for the WEIM hubs, which take their authority’s market-wide region because PacifiCorp has no California citygate.
The expected pricing hub column is withheld from the rule; 4.1 checks
the derivation against it.
HUBS = CAISO.hub_labels
HUB_ENUM = pl.Enum(HUBS)
CORE_HUBS = tuple(label for label in HUBS if not CAISO.hub(label).is_external)
GAS_RULE = pl.DataFrame([asdict(rule) for rule in CAISO_GAS_PAIRING_RULES]).rename(
{
"trading_hub": "trading hub",
"balancing_authority": "balancing authority",
"serving_utility": "serving gas utility",
"fuel_region_prefix": "region id prefix",
"expected_gas_hub": "expected pricing hub",
}
)
html_table(GAS_RULE, kind="reference")
| trading hub | balancing authority | serving gas utility | region id prefix | expected pricing hub |
|---|---|---|---|---|
| NP15 | CISO | PG&E | FRPGE | PG&E Citygate |
| SP15 | CISO | SoCalGas via SCE | FRSCE | SoCal Citygate |
| ZP26 | CISO | PG&E | FRPGE | PG&E Citygate |
| PACE | PACE | — | — | — |
| PACW | PACW | — | — | — |
Source: The CAISO Fuel-Region Reference Table#
The published reference table provides the hub-to-gas mapping.
regions_all = pl.read_parquet(OS_DIR / "caiso_fuel_regions.parquet")
html_table(
regions_all.select(
pl.len().alias("fuel regions"),
pl.col("balancing_authority").n_unique().alias("balancing authorities"),
pl.col("pricing_hub").is_not_null().sum().alias("regions naming a pricing hub"),
pl.col("pricing_hub").is_null().sum().alias("BA-wide regions"),
pl.col("pricing_hub")
.is_null()
.sum()
.over("balancing_authority")
.max()
.alias("most BA-wide regions in any one BA"),
),
kind="reference",
)
| fuel regions | balancing authorities | regions naming a pricing hub | BA-wide regions | most BA-wide regions in any one BA |
|---|---|---|---|---|
| 184 | 24 | 160 | 24 | 1 |
Derive: One Fuel Region per Trading Hub#
resolve_gas_legs applies the pairing rules and returns one fuel
region per trading hub; where several qualify the lowest
fuel_region_id wins, and candidate count above one marks a hub
resolved by that tie-break. resolve_gas_leg_candidates keeps the
regions that lost, which 4.1 prices.
derived_legs = resolve_gas_legs()
candidates = resolve_gas_leg_candidates()
GAS_OF_HUB = dict(derived_legs.as_pairing())
GAS_LABELS = list(dict.fromkeys(GAS_OF_HUB.values()))
REGION_ENUM = pl.Enum(GAS_LABELS)
REFERENCE_GAS = CAISO.reference_gas.label
html_table(
derived_legs.frame().rename(
{
"trading_hub": "trading hub",
"fuel_region_id": "fuel region",
"gas_hub_label": "gas hub",
"balancing_authority": "balancing authority",
"candidate_count": "candidate count",
}
),
kind="reference",
)
| trading hub | fuel region | gas hub | balancing authority | candidate count |
|---|---|---|---|---|
| NP15 | FRPGE1 | PG&E Citygate | CISO | 5 |
| SP15 | FRSCE1 | SoCal Citygate | CISO | 12 |
| ZP26 | FRPGE1 | PG&E Citygate | CISO | 5 |
| PACE | FRPACE | PacifiCorp East | PACE | 1 |
| PACW | FRPACW | PacifiCorp West | PACW | 1 |
audit_market_pairing re-derives the registry’s baseline \(g^*(z)\) map
from the published reference table and compares hub by hub;
raise_for_failure stops the notebook on a stale pairing.
pairing_audit = audit_market_pairing()
pairing_audit.raise_for_failure()
print(pairing_audit.report())
gas leg audit: 5 hubs agree, none diverge
agrees: NP15, PACE, PACW, SP15, ZP26
Source: The Two Price Legs#
MarketContext.for_market builds one row per trading hub per
market-local interval, each power leg joined to its own gas leg.
DATE_RANGE- the trade-day boundary this notebook prices over.market- the CAISO market over that range.
# The pull cache tops up daily, so that boundary is enforced
# here - without it every re-execution would widen the market data and drift
# the numbers cited in prose. The start is 2023-06-01, the first clean
# month above the OASIS LMP retention floor after the 2026-08 backfill.
# This range runs one month longer than notebooks 13 and 14.
DATE_RANGE = (date(2023, 6, 1), date(2026, 7, 31))
market = MarketContext.for_market("CAISO", span=DATE_RANGE)
power_legs = {hub: market.power(hub) for hub in HUBS}
gas_legs = {hub: market.gas(for_hub=hub) for hub in HUBS}
heat_rates = {hub: market.market_implied_heat_rate(hub) for hub in HUBS}
market_data = market.market_data(*heat_rates.values())
Derive: The Hourly Market Data#
Each hub’s power and gas columns are null on every other hub’s rows, so
carried coalesces one index set into the single column the charts and
tables below read.
interval_grid- the complete market-local clock grid, nulls included, which 4.2 audits.hourly- the intervals that carry a published power print.
def carried(indices: dict[str, object], name: str) -> pl.Expr:
"""Coalesce a per-hub index set into one canonical market data column.
Parameters
----------
indices : dict of str to Index
One index per trading hub, keyed by hub label.
name : str
Column name the coalesced series takes.
Returns
-------
pl.Expr
The index value at each row's own delivery point.
"""
return pl.coalesce(
[
index.expression(market_data.schema, market_data.clock)
for index in indices.values()
]
).alias(name)
interval_grid = (
market_data.frame.with_columns(
pl.col("point").cast(HUB_ENUM).alias("hub"),
pl.col("point").replace_strict(GAS_OF_HUB).cast(REGION_ENUM).alias("region"),
market_data.clock.season(pl.col("interval_start_utc")).alias("season"),
carried(power_legs, "lmp_dam"),
carried(gas_legs, "price_gas"),
carried(heat_rates, "mihr"),
)
.rename({"trade_day": "date", "hour_of_day": "hour"})
.collect()
)
hourly = interval_grid.drop_nulls("lmp_dam")
Derive: The Spark Spread at Both Reference Strikes#
SparkSpreadIndex is the project’s definition of \(S(H)\), applied once
against the scenario grid holding both reference strikes.
spark_spread = SparkSpreadIndex.stylized(
strike_heat_rate=scenario("strike_heat_rate")
)
implied_heat_rate = spark_spread.market_implied_heat_rate()
def with_asset_strikes(frame: pl.DataFrame) -> pl.DataFrame:
"""Cross a price frame with both reference strikes and spread it.
Parameters
----------
frame : pl.DataFrame
Any frame carrying ``lmp_dam`` and ``price_gas``.
Returns
-------
pl.DataFrame
The same rows once per strike, with ``strike_heat_rate``, a
readable ``unit`` label and ``spark_spread``.
"""
return frame.join(ASSET_STRIKES.frame, how="cross").with_columns(
spark_spread.combine(pl.col("lmp_dam"), pl.col("price_gas")).alias(
"spark_spread"
),
pl.col("strike_heat_rate").replace_strict(UNIT_LABELS).alias("unit"),
)
hourly_spark = with_asset_strikes(hourly)
def ts_figure(title: str, y_label: str):
"""One-line factory for the standard time-series figure frame."""
return figure(
title=title,
x_axis_type="datetime",
y_axis_label=y_label,
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
4. Quantitative Data Exploration & Assumption Testing#
4.1 Resolving the Hub-to-Gas Pairing#
Gas costs vary by fuel region and power prices follow California’s solar-heavy load profile, so profitability turns on both where the fuel is bought and when the generation occurs.
The map compares California gas service areas with the operating fleet.
MAP_BAS = ["CISO", "PACE", "PACW"]
SIMPLIFY_TOL_M = 1500
gas_fleet = (
pl.read_parquet(OS_DIR / "eia860m_generators.parquet")
.filter(
pl.col("balancing_authority").is_in(MAP_BAS),
pl.col("energy_source_code") == "NG",
pl.col("status") == "(OP) Operating",
)
.with_columns(
pl.col("longitude").map_elements(merc_x, return_dtype=pl.Float64).alias("mx"),
pl.col("latitude").map_elements(merc_y, return_dtype=pl.Float64).alias("my"),
)
)
GAS_TERRITORIES = {
"PG&E": ("PG&E Citygate", SERIES[0]),
"SCG": ("SoCal Citygate", SERIES[1]),
"SDG&E": ("SDG&E (SoCal-adjacent)", SERIES[2]),
}
territories = pl.read_parquet(OS_DIR / "cec_gas_service_areas.parquet").filter(
pl.col("abbrev").is_in(list(GAS_TERRITORIES))
)
ca_fleet = gas_fleet.filter(pl.col("state") == "CA")
fig = figure(
title="Which Citygate a Unit Buys At: CEC Gas Service Areas and the CAISO Gas Fleet",
height=620,
sizing_mode="stretch_width",
match_aspect=True,
toolbar_location="above",
tools="pan,wheel_zoom,box_zoom,reset",
active_scroll="wheel_zoom",
)
fig.add_tile("CartoDB Positron", retina=True)
for abbrev, (label, color) in GAS_TERRITORIES.items():
rows = territories.filter(pl.col("abbrev") == abbrev)
xs, ys = [], []
for geo in rows["geometry_geojson"].to_list():
gx, gy = to_multipolygon(geo, SIMPLIFY_TOL_M)
xs.extend(gx)
ys.extend(gy)
fig.multi_polygons(xs=[xs], ys=[ys], fill_color=color, fill_alpha=0.28,
line_color=color, line_width=1.5, legend_label=label)
plants = ColumnDataSource({
"x": ca_fleet["mx"].to_numpy(), "y": ca_fleet["my"].to_numpy(),
"plant": ca_fleet["plant_name"].to_list(),
"mw": ca_fleet["nameplate_capacity_mw"].to_list(),
"county": ca_fleet["county"].to_list(),
})
r = fig.scatter("x", "y", source=plants, size=6, color=SECONDARY, alpha=0.85,
line_color="#ffffff", line_width=0.5, legend_label="CAISO gas unit")
fig.add_tools(HoverTool(renderers=[r], tooltips=[
("plant", "@plant"), ("county", "@county"), ("capacity", "@mw{0,0} MW"),
]))
style(fig)
# after style(): it re-enables the grid, and a tile map wants neither
# gridlines nor axes
fig.legend.location = "top_right"
fig.axis.visible = False
fig.grid.grid_line_color = None
show(fig)
Sources: CEC California Natural Gas Service Areas (ArcGIS layer), EIA-860M generator inventory; basemap © CARTO, © OpenStreetMap contributors.
The table shows how many regions each hub matched, which was selected, and whether its pricing hub agrees with the withheld expected hub.
receipt = (
candidates.group_by("trading_hub", "serving_utility", "expected_gas_hub")
.agg(
pl.len().alias("matched regions"),
pl.col("pricing_hub").fill_null("(BA-wide)").n_unique().alias("distinct gas markets"),
pl.col("fuel_region_id").filter(pl.col("selected")).first().alias("selected region"),
pl.col("pricing_hub").filter(pl.col("selected")).first().alias("selected prices at"),
pl.col("cap_and_trade_credit").filter(pl.col("selected")).first().alias("cap-and-trade"),
)
.with_columns(
(pl.col("selected prices at") == pl.col("expected_gas_hub")).alias(
"agrees with expected hub"
)
)
.rename({"trading_hub": "hub", "serving_utility": "serving gas utility",
"expected_gas_hub": "expected pricing hub"})
.sort("hub")
)
html_table(receipt, kind="reference")
| hub | serving gas utility | expected pricing hub | matched regions | distinct gas markets | selected region | selected prices at | cap-and-trade | agrees with expected hub |
|---|---|---|---|---|---|---|---|---|
| NP15 | PG&E | PG&E Citygate | 5 | 1 | FRPGE1 | PG&E Citygate | N | True |
| PACE | — | — | 1 | 1 | FRPACE | — | — | — |
| PACW | — | — | 1 | 1 | FRPACW | — | — | — |
| SP15 | SoCalGas via SCE | SoCal Citygate | 12 | 1 | FRSCE1 | SoCal Citygate | N | True |
| ZP26 | PG&E | PG&E Citygate | 5 | 1 | FRPGE1 | PG&E Citygate | N | True |
SP15 has the most pairing discretion, so the chart reprices its daily-mean CCGT spark spread against every SCE fuel region that satisfies the rule. Those alternatives sit outside the market registry, so their gas prices are read from the published fuel-price file while the SP15 power leg still comes from the market data.
# The market's own trade day, not a retyped conversion.
PACIFIC_DATE = CAISO_CLOCK.trade_day(pl.col("interval_start_utc"))
sp15_candidates = (
pl.read_parquet(OS_DIR / "caiso_fuel_prices.parquet")
.join(
candidates.filter(pl.col("trading_hub") == "SP15").select(
"fuel_region_id", "selected"
),
on="fuel_region_id",
)
.with_columns(PACIFIC_DATE.alias("date"))
.filter(pl.col("date").is_between(*DATE_RANGE))
.group_by("fuel_region_id", "selected", "date")
.agg(pl.col("price").first())
)
dam_daily_sp15 = (
hourly.filter(pl.col("hub") == "SP15")
.group_by("date")
.agg(pl.col("lmp_dam").mean())
)
sp15_alternatives = (
dam_daily_sp15.join(sp15_candidates, on="date")
.group_by("fuel_region_id", "selected")
.agg(
pl.col("price").mean().round(2).alias("mean_gas"),
(100 * ((pl.col("lmp_dam") - CCGT_HR * pl.col("price")) > 0).mean())
.round(1)
.alias("pct_itm"),
)
.sort("pct_itm")
)
fig = figure(
title=f"SP15 CCGT Profitability by SCE Fuel Region (HR {CCGT_HR}, daily-mean spark spread)",
y_range=sp15_alternatives["fuel_region_id"].to_list(),
x_axis_label="share of days in the money (%)",
height=380,
sizing_mode="stretch_width",
toolbar_location=None,
tools="",
)
bar_src = ColumnDataSource(
{
"region": sp15_alternatives["fuel_region_id"].to_list(),
"pct": sp15_alternatives["pct_itm"].to_list(),
"gas": sp15_alternatives["mean_gas"].to_list(),
"color": [
SERIES[1] if sel else SERIES[0]
for sel in sp15_alternatives["selected"].to_list()
],
"role": [
"selected region (FRSCE1)" if sel else "alternative SCE region"
for sel in sp15_alternatives["selected"].to_list()
],
}
)
bars = fig.hbar(y="region", right="pct", height=0.7, source=bar_src, color="color",
legend_group="role")
fig.legend.location = "bottom_right"
fig.add_tools(
HoverTool(
renderers=[bars],
tooltips=[("region", "@region"), ("mean gas", "@gas{0.00} $/MMBtu"),
("days ITM", "@pct{0.0}%")],
)
)
show(style(fig))
4.2 Trade-Day Grid and Gas-Leg Integrity#
Section 3 assumes the gas leg is a daily index carried across the day’s hours. The first check is which trade days carry other than 24 published power prints.
day_lengths = hourly.group_by("date", "hub").agg(pl.len().alias("n_hours"))
odd_days = (
day_lengths.filter(pl.col("n_hours") != 24)
.group_by("date", "n_hours")
.agg(pl.len().alias("hubs affected"))
.sort("date")
)
html_table(odd_days)
| date | n_hours | hubs affected |
|---|---|---|
| 2023-11-05 | 25 | 3 |
| 2024-03-10 | 23 | 3 |
| 2024-11-03 | 25 | 3 |
| 2025-03-09 | 23 | 3 |
| 2025-11-02 | 25 | 3 |
| 2025-11-14 | 23 | 3 |
| 2026-03-08 | 23 | 3 |
The audit counts the market data against the clock grid the specification implies, so a duplicate natural key, a dropped interval or a missing print cannot hide inside an average.
def audit_row(hub: str) -> dict[str, object]:
"""One trading hub's market data audit, flattened for the receipt table.
Parameters
----------
hub : str
Trading hub label.
Returns
-------
dict
Counts from :class:`pricing.data.audit.MarketDataAudit`, with the
per-column null and non-positive tallies summed.
"""
audit = market.audit(heat_rates[hub])
return {
"hub": hub,
"rows": audit.rows,
"expected rows": audit.expected_rows,
"trade days": audit.trade_days,
"duplicate natural keys": audit.duplicate_natural_keys,
"missing intervals": audit.missing_intervals,
"interval count anomalies": audit.interval_count_anomalies,
"missing price prints": sum(count for _, count in audit.null_counts),
"non-positive gas prints": sum(
count for _, count in audit.non_positive_counts
),
}
html_table(pl.DataFrame([audit_row(hub) for hub in CORE_HUBS]))
| hub | rows | expected rows | trade days | duplicate natural keys | missing intervals | interval count anomalies | missing price prints | non-positive gas prints |
|---|---|---|---|---|---|---|---|---|
| NP15 | 27,768 | 27,768 | 1,157 | 0 | 0 | 0 | 1 | 0 |
| SP15 | 27,768 | 27,768 | 1,157 | 0 | 0 | 0 | 1 | 0 |
| ZP26 | 27,768 | 27,768 | 1,157 | 0 | 0 | 0 | 1 | 0 |
gas_daily = (
interval_grid.select("region", "date", "price_gas")
.unique(("region", "date"))
.sort("region", "date")
)
gas_audit = (
gas_daily.group_by("region")
.agg(
pl.len().alias("days"),
pl.col("price_gas").min().round(2).alias("min ($/MMBtu)"),
pl.col("price_gas").mean().round(2).alias("mean ($/MMBtu)"),
pl.col("price_gas").max().round(2).alias("max ($/MMBtu)"),
(pl.col("price_gas") <= 0).sum().alias("non-positive days"),
((pl.col("price_gas") <= 0) & (pl.col("date") < WEIM_AVAILABLE_FROM))
.sum()
.alias("non-positive before hub go-live"),
)
.sort("region")
)
html_table(gas_audit)
| region | days | min ($/MMBtu) | mean ($/MMBtu) | max ($/MMBtu) | non-positive days | non-positive before hub go-live |
|---|---|---|---|---|---|---|
| PG&E Citygate | 1,157 | 2.480 | 4.930 | 18.45 | 0 | 0 |
| SoCal Citygate | 1,157 | 5.050 | 7.440 | 16.73 | 0 | 0 |
| PacifiCorp East | 1,157 | 0.760 | 2.490 | 22.62 | 0 | 0 |
| PacifiCorp West | 1,157 | −0.490 | 1.480 | 20.77 | 5 | 5 |
The market data has already collapsed the gas leg to one print per trade day, so the collapse can only be checked against the published file.
gas_prints = (
pl.read_parquet(OS_DIR / "caiso_fuel_prices.parquet")
.filter(pl.col("fuel_region_id").is_in(list(derived_legs.as_fuel_regions().values())))
.with_columns(PACIFIC_DATE.alias("date"))
.filter(pl.col("date").is_between(*DATE_RANGE))
.group_by("fuel_region_id", "date")
.agg(pl.col("price").n_unique().alias("n_prices"), pl.col("price").first())
.sort("fuel_region_id", "date")
.with_columns(
changed=(pl.col("price") != pl.col("price").shift(1)).over("fuel_region_id"),
# Polars weekday(): Monday = 1 through Sunday = 7, so Tuesday
# through Saturday is 2-6
weekday=pl.col("date").dt.weekday(),
)
.group_by("fuel_region_id")
.agg(
pl.len().alias("days"),
(100 * (pl.col("n_prices") > 1).mean())
.round(1)
.alias("pct days with >1 intraday price"),
(100 * pl.col("changed").filter(pl.col("weekday").is_between(2, 6)).mean())
.round(1)
.alias("pct Tue-Sat prints changed"),
(100 * pl.col("changed").filter(~pl.col("weekday").is_between(2, 6)).mean())
.round(1)
.alias("pct Sun-Mon prints changed"),
)
.sort("fuel_region_id")
)
html_table(gas_prints)
| fuel_region_id | days | pct days with >1 intraday price | pct Tue-Sat prints changed | pct Sun-Mon prints changed |
|---|---|---|---|---|
| FRPACE | 1,157 | 0.000 | 92.10 | 0.000 |
| FRPACW | 1,157 | 0.000 | 92.30 | 0.000 |
| FRPGE1 | 1,157 | 0.000 | 94.60 | 0.000 |
| FRSCE1 | 1,157 | 0.000 | 95.50 | 0.000 |
Trade-day length - 1,150 of 1,157 days carry 24 published hours; the rest are the six DST days plus 2025-11-14, which is missing one source hour. The unweighted daily mean handles all seven.
Grid vs print - all three CAISO hubs match their 27,768 expected rows with no duplicate natural keys, missing intervals or interval-count anomalies, and each reports exactly one missing price print: the 2025-11-14 interval exists on the clock, the number does not.
Non-positive gas prices - PacifiCorp West prints a non-positive daily gas index on five September 2025 days, before the PACE/PACW hubs went live, carried across 120 hours where MIHR is null and the spark spread stays defined.
Gas frequency - no region prints more than one distinct price on any day, so the daily gas index loses nothing. The price does not change on Sunday or Monday, so weekend and Monday power hours use Friday’s gas price.
4.3 The Two Legs, Side by Side#
fig = ts_figure("Fuel-Region Gas Price (one daily index print per region)", "$/MMBtu")
renderers = []
for color, region_label in zip(SERIES, GAS_LABELS):
sub = gas_daily.filter(pl.col("region") == region_label)
renderers.append(
fig.line(
sub["date"].cast(pl.Datetime("ms")).to_numpy(),
sub["price_gas"].to_numpy(),
color=color,
line_width=2,
legend_label=region_label,
name=region_label,
)
)
fig.add_layout(Span(location=0, dimension="width", line_color=SECONDARY,
line_dash="dashed", line_width=1))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("date", "@x{%F}"), ("region", "$name"), ("gas", "@y{0.000} $/MMBtu")],
formatters={"@x": "datetime"},
)
)
fig.legend.location = "top_left"
show(style(fig))
In-state premium - SoCal Citygate averages 7.44 $/MMBtu against 4.93 at PG&E Citygate, on pipeline and storage constraints.
Out-of-state discount - PACE and PACW average 2.49 and 1.48 $/MMBtu, and PACW dipped below zero in September 2025. Both sit outside CAISO but participate in the shared WEIM/EDAM markets.
Winter spike - a late-January 2026 cold snap pushed PACE to 17.90 and SoCal Citygate to 11.40 $/MMBtu, narrowing the discount; section 5.4 shows the same event in basis terms.
daily = (
hourly.group_by("hub", "region", "date")
.agg(
pl.col("lmp_dam").mean().alias("lmp_dam"),
pl.col("lmp_dam").min().alias("lmp_min"),
pl.col("lmp_dam").max().alias("lmp_max"),
(pl.col("lmp_dam") < 0).sum().alias("neg_hours"),
pl.len().alias("n_hours"),
pl.col("price_gas").first().alias("price_gas"),
)
.with_columns(
implied_heat_rate.combine(pl.col("lmp_dam"), pl.col("price_gas")).alias("mihr")
)
.sort("date", "hub")
)
fig = ts_figure("Daily Mean DAM LMP by Trading Hub", "$/MWh")
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = daily.filter(pl.col("hub") == hub).sort("date")
renderers.append(
fig.line(
sub["date"].cast(pl.Datetime("ms")).to_numpy(),
sub["lmp_dam"].to_numpy(),
color=color,
line_width=2,
legend_label=hub,
name=hub,
)
)
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("date", "@x{%F}"), ("hub", "$name"), ("LMP", "@y{0.00} $/MWh")],
formatters={"@x": "datetime"},
)
)
fig.legend.location = "top_left"
show(style(fig))
power_audit = (
daily.group_by("hub")
.agg(
pl.len().alias("days"),
pl.col("lmp_dam").mean().round(2).alias("mean daily LMP ($/MWh)"),
pl.col("lmp_min").min().round(2).alias("min hourly"),
pl.col("lmp_max").max().round(2).alias("max hourly"),
(pl.col("lmp_max") - pl.col("lmp_min")).mean().round(2).alias("avg daily range"),
(100 * pl.col("neg_hours").sum() / pl.col("n_hours").sum())
.round(1)
.alias("pct hours negative"),
(100 * (pl.col("lmp_dam") < 0).mean()).round(1).alias("pct days negative"),
)
.sort("hub")
)
html_table(power_audit)
| hub | days | mean daily LMP ($/MWh) | min hourly | max hourly | avg daily range | pct hours negative | pct days negative |
|---|---|---|---|---|---|---|---|
| NP15 | 1,157 | 39.33 | −38.47 | 1,090.90 | 41.74 | 3.40 | 0.200 |
| SP15 | 1,157 | 33.95 | −67.08 | 1,247.63 | 58.41 | 11.70 | 1.000 |
| ZP26 | 1,157 | 33.13 | −68.01 | 1,092.81 | 53.77 | 11.50 | 0.700 |
| PACE | 92 | 18.11 | −42.18 | 91.81 | 31.61 | 9.60 | 8.700 |
| PACW | 92 | 22.75 | −24.50 | 90.34 | 26.10 | 1.50 | 0.000 |
One system - the three CAISO hubs move together across seasons, so profitability turns on the power-to-gas ratio more than on power prices alone.
Levels - CAISO hubs average 33-39 $/MWh against 18.11 at PACE, which is negative on 8.7% of days; SP15 and ZP26 carry the widest average daily ranges, 54-58 dollars per MWh.
Daily means hide the hourly story - SP15 and ZP26 are negative in roughly 12% of hours but only 1.0% and 0.7% of days, midday negatives offset by the evening ramp; PACE is negative on 9.6% of hours.
The heatmap compares the five full-history daily series: three CAISO hub LMPs and two citygate gas indices.
CORR_VARS = [*CORE_HUBS, REFERENCE_GAS, "SoCal Citygate"]
corr_market_data = (
daily.filter(pl.col("hub").is_in(CORE_HUBS))
.with_columns(pl.col("hub").cast(pl.Utf8))
.pivot(on="hub", index="date", values="lmp_dam")
.join(
gas_daily.filter(pl.col("region").is_in([REFERENCE_GAS, "SoCal Citygate"]))
.with_columns(pl.col("region").cast(pl.Utf8))
.pivot(on="region", index="date", values="price_gas"),
on="date",
)
)
corr_square = pl.DataFrame(
{"variable": CORR_VARS}
| {b: [corr_market_data.select(pl.corr(a, b)).item() for a in CORR_VARS] for b in CORR_VARS}
)
corr_heatmap(corr_square, "variable",
"Daily Correlation: Trading-Hub LMPs and Citygate Gas Indices")
Pairwise daily-mean LMP correlations are 0.91-0.96, so the three hubs share one duck-curve pattern; PG&E and SoCal gas correlate at 0.72. The WEIM hubs are excluded, with only 92 days of data.
leg_corr = (
daily.group_by("hub")
.agg(
pl.col("region").first().cast(pl.Utf8).alias("matched gas region"),
pl.len().alias("days"),
pl.corr("lmp_dam", "price_gas").round(2).alias("corr, daily LMP vs gas"),
)
.sort("hub")
)
html_table(leg_corr)
| hub | matched gas region | days | corr, daily LMP vs gas |
|---|---|---|---|
| NP15 | PG&E Citygate | 1,157 | 0.840 |
| SP15 | SoCal Citygate | 1,157 | 0.690 |
| ZP26 | PG&E Citygate | 1,157 | 0.780 |
| PACE | PacifiCorp East | 92 | 0.750 |
| PACW | PacifiCorp West | 92 | 0.460 |
CAISO hubs (0.69-0.84) - each in-state hub tracks its matched citygate closely; PACE (0.75) sits inside that range, above SP15’s 0.69.
PACW is the outlier (0.46) - cheap Rockies/Northwest gas moves less with WEIM power, giving a wider, more volatile spark margin.
5. Implementation & Simulation Logic#
5.1 The Daily Market Data#
The daily market data averages the hourly market data onto the trade
day and crosses it with both reference strikes. It is long in unit,
so the CCGT and the peaker are two rows of one spread definition.
daily_spark = with_asset_strikes(daily)
html_table(
daily_spark.select(
"hub", "region", "date", "lmp_dam", "price_gas", "mihr",
"unit", "strike_heat_rate", "spark_spread",
).head(8)
)
| hub | region | date | lmp_dam | price_gas | mihr | unit | strike_heat_rate | spark_spread |
|---|---|---|---|---|---|---|---|---|
| NP15 | PG&E Citygate | 2023-06-01 | 23.78 | 4.330 | 5.493 | CCGT | 7.00 | −6.53 |
| NP15 | PG&E Citygate | 2023-06-01 | 23.78 | 4.330 | 5.493 | Peaker | 10.50 | −21.68 |
| SP15 | SoCal Citygate | 2023-06-01 | 24.01 | 5.360 | 4.479 | CCGT | 7.00 | −13.51 |
| SP15 | SoCal Citygate | 2023-06-01 | 24.01 | 5.360 | 4.479 | Peaker | 10.50 | −32.27 |
| ZP26 | PG&E Citygate | 2023-06-01 | 21.12 | 4.330 | 4.877 | CCGT | 7.00 | −9.19 |
| ZP26 | PG&E Citygate | 2023-06-01 | 21.12 | 4.330 | 4.877 | Peaker | 10.50 | −24.35 |
| NP15 | PG&E Citygate | 2023-06-02 | 21.27 | 4.070 | 5.227 | CCGT | 7.00 | −7.22 |
| NP15 | PG&E Citygate | 2023-06-02 | 21.27 | 4.070 | 5.227 | Peaker | 10.50 | −21.46 |
5.2 The Daily Headline: Profitability by Trading Hub#
How often each reference unit is profitable on the daily-mean spread.
itm_by_hub = (
daily_spark.group_by("hub", "unit")
.agg((100 * (pl.col("spark_spread") > 0).mean()).round(1).alias("pct_itm"))
.pivot(on="unit", index="hub", values="pct_itm")
.sort("hub")
)
fig = figure(
title=f"Share of Days in the Money by Trading Hub (HR {CCGT_HR} vs {PEAKER_HR}, on the daily-mean spark spread)",
x_range=itm_by_hub["hub"].cast(pl.Utf8).to_list(),
y_axis_label="share of days (%)",
height=360,
sizing_mode="stretch_width",
toolbar_location=None,
tools="",
)
itm_src = ColumnDataSource(
{
"hub": itm_by_hub["hub"].cast(pl.Utf8).to_list(),
"CCGT": itm_by_hub["CCGT"].to_list(),
"Peaker": itm_by_hub["Peaker"].to_list(),
}
)
for offset, col, color in ((-0.16, "CCGT", SERIES[0]), (0.16, "Peaker", SERIES[1])):
r = fig.vbar(x=dodge("hub", offset, range=fig.x_range), top=col, width=0.3,
source=itm_src, color=color, legend_label=col)
fig.add_tools(HoverTool(renderers=[r],
tooltips=[("hub", "@hub"), (col, f"@{col}{{0.0}}%")]))
fig.legend.location = "top_left"
show(style(fig))
summary = (
daily.group_by("hub")
.agg(
pl.len().alias("days"),
pl.col("price_gas").mean().round(2).alias("mean gas ($/MMBtu)"),
pl.col("lmp_dam").mean().round(2).alias("mean LMP ($/MWh)"),
pl.col("mihr").median().round(2).alias("median MIHR"),
)
.join(
daily_spark.group_by("hub", "unit")
.agg(
pl.col("spark_spread").mean().round(2).alias("mean spark"),
(100 * (pl.col("spark_spread") > 0).mean()).round(1).alias("pct days ITM"),
)
.pivot(on="unit", index="hub", values=["mean spark", "pct days ITM"]),
on="hub",
)
.select(
"hub", "days", "mean gas ($/MMBtu)", "mean LMP ($/MWh)", "median MIHR",
pl.col("mean spark_CCGT").alias("mean CCGT spark ($/MWh)"),
pl.col("pct days ITM_CCGT").alias("pct days CCGT ITM"),
pl.col("mean spark_Peaker").alias("mean peaker spark ($/MWh)"),
pl.col("pct days ITM_Peaker").alias("pct days peaker ITM"),
)
.sort("hub")
)
html_table(summary)
| hub | days | mean gas ($/MMBtu) | mean LMP ($/MWh) | median MIHR | mean CCGT spark ($/MWh) | pct days CCGT ITM | mean peaker spark ($/MWh) | pct days peaker ITM |
|---|---|---|---|---|---|---|---|---|
| NP15 | 1,157 | 4.930 | 39.33 | 7.82 | 4.82 | 70.00 | −12.43 | 7.30 |
| SP15 | 1,157 | 7.440 | 33.95 | 4.40 | −18.14 | 8.30 | −44.18 | 0.90 |
| ZP26 | 1,157 | 4.930 | 33.13 | 6.61 | −1.38 | 42.60 | −18.63 | 4.10 |
| PACE | 92 | 1.680 | 18.11 | 11.49 | 6.36 | 79.30 | 0.48 | 63.00 |
| PACW | 92 | 1.350 | 22.75 | 16.23 | 13.30 | 95.70 | 8.58 | 88.00 |
SP15 is the weakest - CCGT spark averages -18.14 $/MWh and is profitable on 8.3% of days, on high gas against low daily-average power.
NP15 is close to breakeven - median MIHR 7.82 against the 7.0 strike, profitable on 70.0% of days but averaging only +4.82 dollars per MWh; ZP26 shares its gas region at lower power prices and is profitable on 42.6% of days.
PACW is the cheap-gas counterexample - median MIHR 16.2, CCGT profitable on 95.7% of days; PACE is a weaker version, pulled down by negative-price days.
The peaker rarely works on daily averages - 0.9-7.3% of days across the CAISO hubs.
5.3 Implied Heat Rate and Spark Spreads Over Time#
core = daily.filter(pl.col("hub").is_in(CORE_HUBS))
DATE_SPAN = (core["date"].min(), core["date"].max())
fig = ts_figure(
"Daily Market-Implied Heat Rate (hourly DAM LMP averaged per day ÷ daily gas index)",
"MIHR (MMBtu/MWh)",
)
season_bands(fig, *DATE_SPAN, label_y=float(core["mihr"].max()))
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = daily.filter(pl.col("hub") == hub)
renderers.append(
fig.line(
sub["date"].cast(pl.Datetime("ms")).to_numpy(),
sub["mihr"].to_numpy(),
color=color,
line_width=2,
legend_label=hub,
name=hub,
)
)
for name, hr in ASSET_HEAT_RATES.items():
fig.add_layout(Span(location=hr, dimension="width", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_layout(Label(x=8, y=hr, x_units="screen", y_units="data", y_offset=5,
text=name, text_font_size="11px", text_color=SECONDARY))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("date", "@x{%F}"), ("hub", "$name"), ("MIHR", "@y{0.00} MMBtu/MWh")],
formatters={"@x": "datetime"},
)
)
fig.legend.location = "top_left"
show(style(fig))
Median MIHR is 7.82 at NP15, 6.61 at ZP26 and 4.40 at SP15: only NP15 sits above the 7.0 CCGT breakeven, and SP15’s power revenue does not cover CCGT gas. MIHR exceeds the 10.5 peaker strike only on rare price-spike days.
ccgt_daily = daily_spark.filter(pl.col("strike_heat_rate") == CCGT_HR)
fig = ts_figure(
f"Daily CCGT Spark Spread (daily-mean hourly DAM LMP − {CCGT_HR} × daily gas index)",
"$/MWh",
)
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = ccgt_daily.filter(pl.col("hub") == hub).sort("date")
renderers.append(
fig.line(
sub["date"].cast(pl.Datetime("ms")).to_numpy(),
sub["spark_spread"].to_numpy(),
color=color,
line_width=2,
legend_label=hub,
name=hub,
)
)
fig.add_layout(Span(location=0, dimension="width", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_layout(BoxAnnotation(top=0, fill_color="#d1495b", fill_alpha=0.05))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("date", "@x{%F}"), ("hub", "$name"), ("spark", "@y{0.00} $/MWh")],
formatters={"@x": "datetime"},
)
)
fig.legend.location = "top_left"
show(style(fig))
NP15 stays close to zero, positive but thin, with some winter losses; SP15 sits below zero for almost the whole period, so a SoCal CCGT is rarely profitable on the daily average. The peaker is an hourly opportunity rather than a daily-average one, which section 6 prices.
5.4 Gas Basis Between the Matched Regions#
basis_indices = {
label: market.gas_basis(label) for label in GAS_LABELS if label != REFERENCE_GAS
}
basis_market_data = market.market_data(*basis_indices.values(), points=(REFERENCE_GAS,))
basis_daily = (
basis_market_data.frame.group_by(pl.col("trade_day").alias("date"))
.agg(
*[
index.expression(basis_market_data.schema, basis_market_data.clock)
.first()
.alias(label)
for label, index in basis_indices.items()
]
)
.sort("date")
.collect()
)
BASIS_SERIES = {
f"{label} - PG&E": (label, color)
for label, color in zip(basis_indices, SERIES[1:])
}
basis_top = basis_daily.select(
pl.max_horizontal(pl.col(label) for label in basis_indices).max()
).item()
fig = ts_figure("Daily Gas Basis to PG&E Citygate (difference of daily gas indices)",
"$/MMBtu")
season_bands(fig, basis_daily["date"].min(), basis_daily["date"].max(),
label_y=float(basis_top))
renderers = []
for label, (column, color) in BASIS_SERIES.items():
renderers.append(
fig.line(
basis_daily["date"].cast(pl.Datetime("ms")).to_numpy(),
basis_daily[column].to_numpy(),
color=color,
line_width=2,
legend_label=label,
name=label,
)
)
fig.add_layout(Span(location=0, dimension="width", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("date", "@x{%F}"), ("pair", "$name"), ("basis", "@y{0.00} $/MMBtu")],
formatters={"@x": "datetime"},
)
)
fig.legend.location = "top_left"
show(style(fig))
basis_means = basis_daily.select(
pl.col(column).mean().round(2).alias(label)
for label, (column, _) in BASIS_SERIES.items()
)
html_table(basis_means)
| SoCal Citygate - PG&E | PacifiCorp East - PG&E | PacifiCorp West - PG&E |
|---|---|---|
| 2.510 | −2.440 | −3.450 |
SoCal Citygate averages +2.51 $/MMBtu over PG&E Citygate, widening sharply in winter, while PacifiCorp East and West sit 2.44 and 3.45 below it - just under six dollars of basis end to end.
6. Quantitative Analysis & Validation#
6.1 The Hourly Rerun: What Daily Averaging Hides#
mihr_dist = (
hourly.group_by("hub")
.agg(
pl.len().alias("hours"),
pl.col("mihr").quantile(0.1).round(2).alias("MIHR p10"),
pl.col("mihr").median().round(2).alias("MIHR median"),
pl.col("mihr").quantile(0.9).round(2).alias("MIHR p90"),
(100 * (pl.col("lmp_dam") < 0).mean()).round(1).alias("pct hours LMP < 0"),
)
.join(
hourly_spark.group_by("hub", "unit")
.agg((100 * (pl.col("spark_spread") > 0).mean()).round(1).alias("pct_itm"))
.pivot(on="unit", index="hub", values="pct_itm")
.rename({"CCGT": "pct hours CCGT ITM", "Peaker": "pct hours peaker ITM"}),
on="hub",
)
.sort("hub")
)
# The hourly MIHR spread against the two strikes: whichever part of a bar
# sits above a dashed line is the part of the year that unit is in the money.
fig = figure(
title="Hourly Market-Implied Heat Rate, p10 to p90 by Trading Hub (hourly DAM LMP ÷ daily gas index)",
x_range=mihr_dist["hub"].cast(pl.Utf8).to_list(),
y_axis_label="MIHR (MMBtu/MWh)",
height=380,
sizing_mode="stretch_width",
toolbar_location=None,
tools="",
)
dist_src = ColumnDataSource(
{
"hub": mihr_dist["hub"].cast(pl.Utf8).to_list(),
"p10": mihr_dist["MIHR p10"].to_list(),
"median": mihr_dist["MIHR median"].to_list(),
"p90": mihr_dist["MIHR p90"].to_list(),
}
)
fig.vbar(x="hub", bottom="p10", top="p90", width=0.5, source=dist_src,
color=SERIES[0], fill_alpha=0.35, line_color=SERIES[0],
legend_label="p10-p90 hourly range")
med = fig.scatter(x="hub", y="median", source=dist_src, size=12, marker="dash",
line_color=INK, line_width=3, legend_label="median hour")
for name, hr in ASSET_HEAT_RATES.items():
fig.add_layout(Span(location=hr, dimension="width", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_layout(Label(x=8, y=hr, x_units="screen", y_units="data", y_offset=5,
text=name, text_font_size="11px", text_color=SECONDARY))
fig.legend.location = "top_left"
fig.add_tools(HoverTool(
renderers=[med],
tooltips=[("hub", "@hub"), ("p10", "@p10{0.00}"), ("median", "@median{0.00}"),
("p90", "@p90{0.00}")],
))
show(style(fig))
The hourly p10-p90 ranges are wide enough to cover most of the strike range within a single hub. SP15 is still the lowest hub, but 14.5% of its hours are above breakeven, putting its p90 above the 7.0 CCGT strike; NP15 and ZP26 have medians above that strike with p90s near the peaker strike, and PACE and PACW run higher still on cheap gas.
html_table(
mihr_dist.select(
"hub", "hours", "pct hours LMP < 0", "pct hours CCGT ITM", "pct hours peaker ITM"
)
)
| hub | hours | pct hours LMP < 0 | pct hours CCGT ITM | pct hours peaker ITM |
|---|---|---|---|---|
| NP15 | 27,767 | 3.40 | 68.20 | 14.10 |
| SP15 | 27,767 | 11.70 | 14.50 | 1.60 |
| ZP26 | 27,767 | 11.50 | 59.80 | 11.60 |
| PACE | 2,208 | 9.60 | 84.10 | 58.50 |
| PACW | 2,208 | 1.50 | 92.20 | 81.80 |
optionality = (
hourly_spark.group_by("hub", "unit")
.agg(
pl.col("spark_spread").mean().round(2).alias("always-on"),
pl.col("spark_spread").clip(lower_bound=0.0).mean().round(2).alias("perfect dispatch"),
)
.pivot(on="unit", index="hub", values=["always-on", "perfect dispatch"])
.select(
"hub",
pl.col("always-on_CCGT").alias("CCGT always-on ($/MWh)"),
pl.col("perfect dispatch_CCGT").alias("CCGT perfect dispatch ($/MWh)"),
pl.col("always-on_Peaker").alias("peaker always-on ($/MWh)"),
pl.col("perfect dispatch_Peaker").alias("peaker perfect dispatch ($/MWh)"),
)
.sort("hub")
)
html_table(optionality)
| hub | CCGT always-on ($/MWh) | CCGT perfect dispatch ($/MWh) | peaker always-on ($/MWh) | peaker perfect dispatch ($/MWh) |
|---|---|---|---|---|
| NP15 | 4.82 | 8.59 | −12.43 | 1.800 |
| SP15 | −18.13 | 2.66 | −44.18 | 1.000 |
| ZP26 | −1.38 | 7.24 | −18.63 | 1.430 |
| PACE | 6.36 | 9.43 | 0.48 | 5.070 |
| PACW | 13.30 | 13.81 | 8.58 | 9.650 |
AtTheMoneyStrikeRule reads each hub’s median hourly MIHR at the
nearest rank, so the at-the-money strike is one level per delivery
point and always a level the market printed.
AT_THE_MONEY = AtTheMoneyStrikeRule()
def prevailing_level(hub: str) -> float:
"""The market-implied heat rate this hub prevails at.
Parameters
----------
hub : str
Trading hub label.
Returns
-------
float
The median hourly MIHR at the nearest rank, in MMBtu/MWh.
"""
market_data = market.market_data(heat_rates[hub])
return AT_THE_MONEY.levels(heat_rates[hub], market_data).level(hub)
html_table(
pl.DataFrame(
{
"hub": list(CORE_HUBS),
"prevailing MIHR (median hour, MMBtu/MWh)": [
round(prevailing_level(hub), 2) for hub in CORE_HUBS
],
}
)
)
| hub | prevailing MIHR (median hour, MMBtu/MWh) |
|---|---|
| NP15 | 8.010 |
| SP15 | 4.960 |
| ZP26 | 7.670 |
6.2 Pricing the HRCO: Strip Value and Settlement Frequency#
One HistoricalSettlementEngine prices the same contract at both
settlement frequencies, so the chart’s solid and dashed curves are
\(V_{\text{hourly}}(H)\) and \(V_{\text{daily}}(H)\) and the gap between
them is \(\Pi(H)\).
LADDER = HistoricalSettlementEngine(
schedules=(SettlementFrequency.HOURLY, SettlementFrequency.DAILY)
)
def settlement_ladder(hubs, strikes) -> pl.DataFrame:
"""Price the hourly and daily rungs across a strike grid, per hub.
Parameters
----------
hubs : iterable of str
Trading hub labels, one priced instrument each.
strikes : iterable of float
Strike heat rates the scenario grid sweeps.
Returns
-------
pl.DataFrame
One row per hub per strike, carrying ``value_hourly``,
``value_daily`` and ``optionality_premium``.
"""
grid = ScenarioGrid(strike_heat_rate=list(strikes))
return (
pl.concat(
[
market.heat_rate_call(hub, scenario("strike_heat_rate"))
.with_engine(LADDER)
.sweep(grid)
for hub in hubs
]
)
.rename({"point": "hub"})
.with_columns(pl.col("hub").cast(HUB_ENUM))
)
STRIKE_GRID = [round(3.5 + 0.5 * i, 1) for i in range(18)]
sweep = settlement_ladder(CORE_HUBS, STRIKE_GRID)
fig = figure(
title="What an HRCO Is Worth: Expected Strip Value by Strike (Solid: Hourly settlement, Dashed: Daily-Average)",
x_axis_label="strike heat rate H (MMBtu/MWh)",
y_axis_label="expected payoff ($/MWh)",
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = sweep.filter(pl.col("hub") == hub).sort("strike_heat_rate")
renderers.append(
fig.line(sub["strike_heat_rate"].to_numpy(), sub["value_hourly"].to_numpy(),
color=color, line_width=2.5, legend_label=hub,
name=f"{hub} Hourly HRCO strip")
)
renderers.append(
fig.line(sub["strike_heat_rate"].to_numpy(), sub["value_daily"].to_numpy(),
color=color, line_width=1.5, line_dash="dashed",
name=f"{hub} Daily")
)
for hr in ASSET_HEAT_RATES.values():
fig.add_layout(Span(location=hr, dimension="height", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("strike", "@x{0.0}"), ("series", "$name"),
("value", "@y{0.00} $/MWh")],
)
)
fig.legend.location = "top_right"
show(style(fig))
At the CCGT strike the hourly strip is worth $8.59/MWh at NP15 (about 75,000 dollars per MW-year), 7.24 at ZP26 and 2.66 at SP15 (about 23,300 per MW-year). Value falls with the strike at every hub and the ordering NP15 > ZP26 > SP15 never changes.
The annual figure is the strip value times the settlement periods the schedule delivers in a year, read off the market clock: 8,765.8 for a Pacific hourly schedule, not 8,760.
INTERVALS_PER_YEAR = market.heat_rate_call(
CORE_HUBS[0], CCGT_HR
).instrument.schedule.settled_intervals_per_year(market_data.clock)
settle = (
settlement_ladder(HUBS, ASSET_HEAT_RATES.values())
.with_columns(
(pl.col("value_hourly") * INTERVALS_PER_YEAR)
.round(0)
.alias("hourly strip ($/MW-year)")
)
.with_columns(
pl.col("value_hourly", "value_daily", "optionality_premium").round(2)
)
.select("hub", "strike_heat_rate", "value_hourly", "value_daily",
"optionality_premium", "hourly strip ($/MW-year)")
.sort("hub", "strike_heat_rate")
)
html_table(settle)
| hub | strike_heat_rate | value_hourly | value_daily | optionality_premium | hourly strip ($/MW-year) |
|---|---|---|---|---|---|
| NP15 | 7.00 | 8.59 | 6.80 | 1.790 | 75,342.00 |
| NP15 | 10.50 | 1.80 | 0.70 | 1.100 | 15,739.00 |
| SP15 | 7.00 | 2.66 | 1.27 | 1.390 | 23,286.00 |
| SP15 | 10.50 | 1.00 | 0.29 | 0.710 | 8,743.00 |
| ZP26 | 7.00 | 7.24 | 3.69 | 3.550 | 63,454.00 |
| ZP26 | 10.50 | 1.43 | 0.43 | 0.990 | 12,502.00 |
| PACE | 7.00 | 9.43 | 8.25 | 1.180 | 82,673.00 |
| PACE | 10.50 | 5.07 | 3.68 | 1.380 | 44,431.00 |
| PACW | 7.00 | 13.81 | 13.38 | 0.430 | 121,090.00 |
| PACW | 10.50 | 9.65 | 8.96 | 0.690 | 84,600.00 |
value_hourly matches the perfect-dispatch values in 6.1. At the CCGT
strike, daily settlement gives up 1.39 dollars per MWh of SP15’s 2.66,
about half the strip, against $1.79/MWh at NP15, 3.55 at ZP26, and
only 1.18 at PACE and 0.43 at PACW: the premium is largest where the
spark spread changes most across the day.
fig = figure(
title="Intraday Optionality Premium by Strike (hourly vs daily settlement)",
x_axis_label="strike heat rate H (MMBtu/MWh)",
y_axis_label="premium ($/MWh)",
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = sweep.filter(pl.col("hub") == hub).sort("strike_heat_rate")
renderers.append(
fig.line(sub["strike_heat_rate"].to_numpy(),
sub["optionality_premium"].to_numpy(),
color=color, line_width=2, legend_label=hub, name=hub)
)
for hr in ASSET_HEAT_RATES.values():
fig.add_layout(Span(location=hr, dimension="height", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("strike", "@x{0.0}"), ("hub", "$name"),
("premium", "@y{0.00} $/MWh")],
)
)
fig.legend.location = "top_right"
show(style(fig))
The premium is small at very low and very high strikes, where the spread rarely changes sign within a day, and largest in the middle. ZP26 peaks highest, on NP15 gas against lower power prices. SP15’s premium is still rising at the lowest tested strike of 3.5 ($3.50/MWh), below the 3.412 physical floor; within the buildable range it is smaller, 1.39 dollars per MWh at HR 7.0.
6.3 The Hour-of-Day MIHR Profile#
profile = (
hourly.filter(pl.col("hub").is_in(CORE_HUBS))
.group_by("hour", "hub")
.agg(pl.col("mihr").median().alias("median_mihr"))
.sort("hub", "hour")
)
fig = figure(
title="Median MIHR by Hour of Day (hourly DAM LMP ÷ daily gas index, Pacific hours)",
x_axis_label="hour of day (Pacific)",
y_axis_label="median MIHR (MMBtu/MWh)",
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = profile.filter(pl.col("hub") == hub)
renderers.append(
fig.line(sub["hour"].to_numpy(), sub["median_mihr"].to_numpy(),
color=color, line_width=2, legend_label=hub, name=hub)
)
for name, hr in ASSET_HEAT_RATES.items():
fig.add_layout(Span(location=hr, dimension="width", line_color=SECONDARY,
line_dash="dashed", line_width=1.5))
fig.add_layout(Label(x=8, y=hr, x_units="screen", y_units="data", y_offset=5,
text=name, text_font_size="11px", text_color=SECONDARY))
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("hour", "@x"), ("hub", "$name"), ("median MIHR", "@y{0.00}")],
)
)
fig.legend.location = "top_left"
show(style(fig))
midday_floor = (
profile.filter(pl.col("hour").is_between(9, 15))
.group_by("hub")
.agg(
pl.col("median_mihr")
.max()
.round(2)
.alias("highest midday median MIHR (hours 9-15)")
)
.with_columns(
(pl.col("highest midday median MIHR (hours 9-15)") < PHYSICAL_FLOOR).alias(
"below the 3.412 physical floor"
)
)
.sort("hub")
)
html_table(midday_floor)
| hub | highest midday median MIHR (hours 9-15) | below the 3.412 physical floor |
|---|---|---|
| NP15 | 6.720 | False |
| SP15 | 2.690 | True |
| ZP26 | 4.070 | False |
All three hubs share the duck-curve shape: MIHR falls around midday and recovers in the evening. SP15’s median falls as low as 1.48 MMBtu/MWh at 1 p.m., below the 3.412 physical floor; the highest midday medians are 6.72 at NP15 and 4.07 at ZP26, both still under the 7.0 CCGT reference.
6.5 Tail Concentration of the Strip Value#
value_concentration sorts each hub’s settled hourly payoffs largest
first and cumulates them into \(C(q)\) from section 3.6, against a dashed
equal-contribution benchmark. Core hubs only: the WEIM hubs’ 92-day
history is too short for tail shares.
CCGT_STRIPS = {hub: market.heat_rate_call(hub, CCGT_HR) for hub in CORE_HUBS}
strip_book = pl.concat(
[
strip.instrument.compile(strip.market_data).book().select("point", "settled_value")
for strip in CCGT_STRIPS.values()
]
).collect()
concentration = value_concentration(strip_book, by=("point",))
tail_curve = concentration.curve.rename({"point": "hub"}).with_columns(
pl.col("hub").cast(HUB_ENUM),
(100 * pl.col("period_share")).alias("pct_hours"),
(100 * pl.col("value_share")).alias("pct_value"),
)
fig = figure(
title=f"Share of Strip Value From the Top Hours (CCGT HR {CCGT_HR}, hourly settlement)",
x_axis_label="top share of hours (%)",
y_axis_label="share of total strip value (%)",
x_axis_type="log",
x_range=(0.01, 100),
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
equal_grid = [0.01 * 10 ** (i / 10) for i in range(41)]
fig.line(equal_grid, equal_grid, color=SECONDARY, line_width=1.5,
line_dash="dashed", legend_label="equal contribution")
renderers = []
for color, hub in zip(SERIES, CORE_HUBS):
sub = tail_curve.filter(pl.col("hub") == hub)
renderers.append(
fig.line(sub["pct_hours"].to_numpy(), sub["pct_value"].to_numpy(),
color=color, line_width=2, legend_label=hub, name=hub)
)
fig.add_tools(
HoverTool(
renderers=renderers,
tooltips=[("hub", "$name"), ("top hours", "$x{0.000}%"),
("share of value", "$y{0.0}%")],
)
)
fig.legend.location = "top_left"
show(style(fig))
SP15 is a tail-risk instrument - half its strip value arrives in the top 1.0% of hours, 278 of the 27,767 in the sample, the top 5% carries 81.6%, and the single best hour carries 1.6%. Selling an SP15 HRCO is being short a few hundred evenings, not a year of margin.
NP15 and ZP26 are broad-based - they need the top 12.7% and 11.4% of hours to reach half value and their best hour carries 0.4-0.5%, but both curves still run well above equal contribution.
TAIL_SHARES = [0.001, 0.01, 0.05, 0.10]
tail_summary = (
tail_curve.group_by("hub", maintain_order=True)
.agg(
pl.len().alias("hours"),
(100 * pl.col("value_share").first()).round(1).alias("best single hour (%)"),
)
.join(
concentration.shares_at(*TAIL_SHARES)
.rename({"point": "hub"})
.with_columns(
pl.col("hub").cast(HUB_ENUM),
(100 * pl.col("value_share")).round(1).alias("value_share"),
pl.format(
"top {}% hours (%)", (100 * pl.col("period_share")).round(1)
).alias("bucket"),
)
.pivot(on="bucket", index="hub", values="value_share"),
on="hub",
)
.join(
concentration.half_value_periods.rename({"point": "hub"}).select(
pl.col("hub").cast(HUB_ENUM),
(100 * pl.col("share_to_half_value"))
.round(2)
.alias("hours to half the value (%)"),
),
on="hub",
)
.sort("hub")
)
html_table(tail_summary)
| hub | hours | best single hour (%) | top 0.1% hours (%) | top 1.0% hours (%) | top 5.0% hours (%) | top 10.0% hours (%) | hours to half the value (%) |
|---|---|---|---|---|---|---|---|
| NP15 | 27,767 | 0.400 | 4.40 | 13.50 | 30.30 | 44.00 | 12.69 |
| SP15 | 27,767 | 1.600 | 17.10 | 50.00 | 81.60 | 96.40 | 1.00 |
| ZP26 | 27,767 | 0.500 | 5.00 | 14.40 | 32.10 | 46.60 | 11.40 |