01. Data Exploration - CAISO

Authors: Brian Nguyen · Claire Kim · Created: 2026-07-23 · Tags: #Exploration #Survey #Tail-Risk #Weather

2. Objective & Input Framework

Objective: cover what the chartbook skips - it already carries the DART spread by trading hub, SP15 DAM vs RT, the SP15 market-implied heat rate and the fuel-mix views. The path here:

1. Audit the coverage of every cached dataset.

2. Audit both tails of the hourly DAM price distribution against a matched normal benchmark, and check where the 1 $/MWh log-return floor of the Margrabe volatility estimator binds.

3. Map the zonal weather layer against the market geography.

4. Survey location spreads, the TB4 storage spread, net load, forced vs planned outages, and the citygate gas regions.

5. Measure weather’s coincidence with load, outages and hub prices - the covariance term a quanto pays for.

Data & Feature Catalog: all datasets are cached as typed parquet (coverage in section 4.1); timestamps on disk are naive UTC. The weather tables come from Claire Kim’s Open-Meteo pipeline (the collab_weather_pipeline doit task), not from src/pulls/.

Dataset	Granularity	Role
every data/opensource/*.parquet	varies	the 4.1 coverage audit reads all of them
caiso_lmp_day_ahead_hourly	hourly (hub APNode)	DAM LMPs - the tail and floor audits, location spreads, TB4
caiso_lmp_real_time_15_min	15-min (hub APNode)	RT leg of the DAM-RT spread in 6.8
caiso_load	sub-hourly (system)	system load, averaged to hourly for net load and daily for 6.7
caiso_load_hourly	hourly (TAC area)	the TAC layer’s unit count in the 4.4 geography audit
caiso_renewables_hourly	hourly	solar and wind generation for net load
caiso_fuel_mix	5-min (system)	gas, solar and wind generation for the 6.7 coincidence checks
caiso_outages_aggregated_hourly_derived	hourly (system)	forced vs planned outage MW
caiso_fuel_prices	hourly (fuel region)	citygate gas indices and the gas-leg floor audit
caiso_weather_daily_region / caiso_weather_hourly_region	daily / hourly (market zone)	weighted zonal weather indices - CDD/HDD and CDH/HDH
caiso_weather_daily_weather_zone	daily (NOAA proxy zone)	the weather-side geography
weather_zone_crosswalk / weather_market_zone_basis_risk	reference	Open-Meteo point assignments and the zone-mapping audit

Setup: Inputs and Source Data

Declare the hub universe and the Pacific clock, then load the shared DAM power leg.

from datetime import date
from pathlib import Path

import numpy as np
import polars as pl
from bokeh.io import output_notebook, show
from bokeh.layouts import row
from bokeh.models import FactorRange, HoverTool, Span
from bokeh.plotting import figure
from bokeh.resources import INLINE
from scipy.stats import norm

from chart_theme import SERIES, corr_heatmap, html_table, style
from markets.definitions.caiso import CAISO_CLOCK
from markets.definitions.caiso import CAISO
from markets import MarketContext
from pricing.math.regression import least_squares_weights, partial_correlation
from settings import config

DATA_DIR = config("DATA_DIR")
OS_DIR = DATA_DIR / "opensource"
PROCESSED_DIR = DATA_DIR / "processed"

output_notebook(INLINE, hide_banner=True)

Declare: Trading Hubs and the Pacific Clock

• HUBS - the three CAISO trading-hub APNodes with full history.

• ALL_HUBS - adds the two PacifiCorp WEIM hubs, live since 2026-05-01.

• PACIFIC - the US/Pacific view of the naive-UTC timestamps; trade days are Pacific.

• DATE_RANGE - 2023-06-01, the first clean month above the OASIS LMP retention floor after the 2026-08 backfill, through 2026-07-31. Fixing it stops the daily cache top-up from drifting the numbers cited in prose.

# The market defines its own hubs, so the APNode ids are read off the definition
# rather than retyped. Definition order is preserved, which the colour zips below
# rely on. HUBS is the three with full history; the WEIM pair went live 2026-05-01.
ALL_HUBS = {CAISO.hub(label).identifier: label for label in CAISO.hub_labels}
HUBS = {
    CAISO.hub(label).identifier: label
    for label in CAISO.hub_labels
    if not CAISO.hub(label).is_external
}

# The market's own clock rather than a retyped conversion, so this notebook and
# the library's trade_day key cannot drift apart.
PACIFIC = CAISO_CLOCK.local(pl.col("interval_start_utc"))
DATE_RANGE = (date(2023, 6, 1), date(2026, 7, 31))

Source: The DAM Power Leg

• dam - hourly day-ahead LMPs for the five hub APNodes on the Pacific trade-day clock; core narrows to the three full-history hubs. Every other dataset is scanned inside the one section that uses it.

dam = (
    pl.read_parquet(OS_DIR / "caiso_lmp_day_ahead_hourly.parquet")
    .filter(pl.col("location").is_in(list(ALL_HUBS)))
    .with_columns(
        pl.col("location").replace(ALL_HUBS).alias("hub"),
        PACIFIC.dt.date().alias("date"),
    )
    .filter(pl.col("date").is_between(*DATE_RANGE))
)
core = dam.filter(pl.col("hub").is_in(list(HUBS.values())))

Source: The Zonal Weather Layer

Shared by sections 4.4 and 6.6-6.8; timestamps arrive tz-aware UTC and are normalized to the market data’s naive-UTC convention.

• weather_daily - one row per day per market zone (CAISO, NP15, SP15, ZP26): weighted temperature, HDD65/CDD65, and temperature dispersion across the zone’s Open-Meteo points.

• weather_hourly - the hourly version, with CDH65/HDH65 degree hours.

• weather_zone_daily - the same daily features on the NOAA proxy zones.

• weather_crosswalk - the Open-Meteo points with their market-zone and NOAA-zone assignments and weights.

• basis_risk - the market-zone to weather-zone mapping audit.

weather_daily = (
    pl.read_parquet(PROCESSED_DIR / "caiso_weather_daily_region.parquet")
    .filter(pl.col("date").is_between(*DATE_RANGE))
)
weather_hourly = (
    pl.read_parquet(PROCESSED_DIR / "caiso_weather_hourly_region.parquet")
    .filter(pl.col("date").is_between(*DATE_RANGE))
    .with_columns(
        pl.col("timestamp_utc")
        .dt.replace_time_zone(None)
        .dt.cast_time_unit("us")
        .alias("interval_start_utc")
    )
)
weather_zone_daily = (
    pl.read_parquet(PROCESSED_DIR / "caiso_weather_daily_weather_zone.parquet")
    .filter(pl.col("date").is_between(*DATE_RANGE))
)
weather_crosswalk = pl.read_parquet(PROCESSED_DIR / "weather_zone_crosswalk.parquet")
basis_risk = pl.read_parquet(PROCESSED_DIR / "weather_market_zone_basis_risk.parquet")

Derive: The Market Context

Resolving DATE_RANGE against the CAISO market definition gives the object that settles the spread instruments in sections 6.1 and 6.2.

• market - CAISO over DATE_RANGE, carrying the trading hubs, the Pacific trade-day clock and the settlement conventions. The frames above stay for the distribution audits, which need raw hourly prints.

market = MarketContext.for_market("CAISO", span=DATE_RANGE)

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 Coverage Audit

Row counts and time span of every cached dataset, read as the cache stands. Everything below fixes the analysis range to DATE_RANGE.

def coverage_row(path: Path) -> pl.DataFrame:
    lf = pl.scan_parquet(path)
    schema = lf.collect_schema()
    time_col = next(
        (c for c in ("interval_start_utc", "publish_time_utc") if c in schema), None
    )
    time_expr = (
        pl.col(time_col) if time_col else pl.lit(None, dtype=pl.Datetime("us"))
    )
    return (
        lf.select(
            pl.lit(path.stem).alias("dataset"),
            pl.len().alias("rows"),
            # datasets mix nanosecond/microsecond datetime units; cast to one
            # unit so the rows concat
            time_expr.min().cast(pl.Datetime("us")).alias("start_utc"),
            time_expr.max().cast(pl.Datetime("us")).alias("end_utc"),
        )
        .collect()
    )


coverage = pl.concat([coverage_row(p) for p in sorted(OS_DIR.glob("*.parquet"))])
html_table(coverage)

15 rows × 4 columns

dataset	rows	start_utc	end_utc
caiso_curtailed_non_operational_generator_report	3,655,547	2021-06-18 07:00	2026-08-26 07:00
caiso_fuel_mix	880,737	2018-04-10 07:00	2026-08-27 06:55
caiso_fuel_prices	10,876,916	2018-01-01 08:00	2026-08-27 06:00
caiso_fuel_regions	184	—	—
caiso_lmp_day_ahead_hourly	92,349	2023-05-11 07:00	2026-08-27 06:00
caiso_lmp_real_time_15_min	369,394	2023-05-11 07:00	2026-08-27 06:45
caiso_lmp_real_time_5_min	1,108,184	2023-05-11 07:00	2026-08-27 06:55
caiso_load	880,751	2018-04-10 07:00	2026-08-27 06:55
caiso_load_hourly	1,944,693	2018-01-01 08:00	2026-08-27 06:00
caiso_outages_aggregated_hourly_derived	45,479	2021-06-17 07:00	2026-08-26 06:00
caiso_renewables_hourly	116,279	2023-05-01 07:00	2026-08-27 06:00
cec_gas_service_areas	13	—	—
eia860m_generators	28,318	—	—
openmeteo_temperature_daily	18,944	2018-01-01 08:00	2026-08-26 07:00
us_ba_boundaries	85	—	—

4.2 Price Tails Against a Normal Benchmark

Tail behavior decides what the later option notebooks are worth: HRCO and quanto payoffs monetize extremes, and Margrabe assumes roughly lognormal legs. This section measures how the hourly DAM price distribution departs from a normal benchmark with the same mean and standard deviation, per trading hub.

tail_stats = (
    core.group_by("hub")
    .agg(
        pl.len().alias("hours"),
        pl.col("lmp").mean().alias("mean"),
        pl.col("lmp").std().alias("sd"),
        pl.col("lmp").skew().alias("skewness"),
        pl.col("lmp").kurtosis().alias("excess kurtosis"),
        pl.col("lmp").min().alias("min"),
        pl.col("lmp").quantile(0.999).alias("p99.9"),
        pl.col("lmp").max().alias("max"),
        (100 * (pl.col("lmp") < 0).mean()).alias("pct hours negative"),
    )
    .sort("hub")
)
tail_curves = (
    core.select("hub", "lmp")
    .sort("hub", "lmp")
    .with_columns(
        (100 * (pl.int_range(pl.len()).over("hub") + 1) / pl.len().over("hub"))
        .alias("pct_below"),
        (100 * (pl.len().over("hub") - pl.int_range(pl.len()).over("hub"))
         / pl.len().over("hub")).alias("pct_above"),
    )
)

The full hourly price distribution per hub on a log count scale, with the matched normal dashed; the axis is truncated at 400 $/MWh for readability.

histogram_fig = figure(
    title="Hourly DAM Price Distribution Against the Matched Normal",
    x_axis_label="hourly DAM LMP ($/MWh)",
    y_axis_label="share of hours per $5 bin (%)",
    y_axis_type="log",
    y_range=(0.0008, 40),
    x_range=(-85, 400),
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
bin_edges = np.arange(-85, 405, 5.0)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
histogram_renderers = []
for color, hub in zip(SERIES, HUBS.values()):
    values = core.filter(pl.col("hub") == hub)["lmp"].to_numpy()
    counts, _ = np.histogram(values, bins=bin_edges)
    share = 100 * counts / len(values)
    share_plot = np.where(share > 0, share, np.nan)
    histogram_renderers.append(
        histogram_fig.line(bin_centers, share_plot, color=color, line_width=1.8,
                           legend_label=hub, name=hub)
    )
    mu, sd = tail_stats.filter(pl.col("hub") == hub).select("mean", "sd").row(0)
    histogram_fig.line(bin_centers, 100 * 5.0 * norm.pdf(bin_centers, mu, sd),
                       color=color, line_width=1.2, line_dash="dashed")
histogram_fig.add_tools(
    HoverTool(renderers=histogram_renderers,
              tooltips=[("hub", "$name"), ("price bin", "$x{0} $/MWh"),
                        ("share of hours", "$y{0.0000}%")])
)
histogram_fig.legend.location = "top_right"
show(style(histogram_fig))

left_fig = figure(
    title="Left Tail: Share of Hours at or Below Price",
    x_axis_label="hourly DAM LMP ($/MWh)",
    y_axis_label="share of hours (%)",
    y_axis_type="log",
    y_range=(0.008, 130),
    x_range=(float(tail_stats["min"].min()) - 2, 30),
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
right_fig = figure(
    title="Right Tail: Share of Hours at or Above Price",
    x_axis_label="hourly DAM LMP ($/MWh)",
    y_axis_label="share of hours (%)",
    y_axis_type="log",
    y_range=(0.008, 130),
    x_range=(25, float(tail_stats["max"].max()) + 5),
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
left_fig.add_layout(Span(location=0, dimension="height", line_color="#c3c2b7",
                         line_width=1))
renderers_left, renderers_right = [], []
for color, hub in zip(SERIES, HUBS.values()):
    sub = tail_curves.filter(pl.col("hub") == hub)
    mu, sd = tail_stats.filter(pl.col("hub") == hub).select("mean", "sd").row(0)
    grid = np.linspace(float(sub["lmp"].min()), float(sub["lmp"].max()), 400)
    renderers_left.append(
        left_fig.line(sub["lmp"].to_numpy(), sub["pct_below"].to_numpy(),
                      color=color, line_width=2, legend_label=hub, name=hub)
    )
    left_fig.line(grid, 100 * norm.cdf(grid, mu, sd), color=color, line_width=1.5,
                  line_dash="dashed")
    renderers_right.append(
        right_fig.line(sub["lmp"].to_numpy(), sub["pct_above"].to_numpy(),
                       color=color, line_width=2, legend_label=hub, name=hub)
    )
    right_fig.line(grid, 100 * norm.sf(grid, mu, sd), color=color, line_width=1.5,
                   line_dash="dashed")
left_fig.add_tools(
    HoverTool(renderers=renderers_left, tooltips=[("hub", "$name"),
              ("price", "$x{0.00} $/MWh"), ("share below", "$y{0.000}%")])
)
right_fig.add_tools(
    HoverTool(renderers=renderers_right, tooltips=[("hub", "$name"),
              ("price", "$x{0.00} $/MWh"), ("share above", "$y{0.000}%")])
)
left_fig.legend.location = "top_left"
right_fig.legend.location = "bottom_left"
show(row(style(left_fig), style(right_fig), sizing_mode="stretch_width"))

tail_audit = (
    tail_stats.with_columns(
        pl.Series(
            "normal pct negative",
            100 * norm.cdf(0.0, tail_stats["mean"].to_numpy(),
                           tail_stats["sd"].to_numpy()),
        ),
        pl.Series(
            "normal p99.9 ($/MWh)",
            tail_stats["mean"].to_numpy()
            + norm.ppf(0.999) * tail_stats["sd"].to_numpy(),
        ),
    )
    .select(
        "hub",
        "hours",
        pl.col("mean").round(2).alias("mean ($/MWh)"),
        pl.col("skewness").round(2),
        pl.col("excess kurtosis").round(2),
        pl.col("min").alias("min ($/MWh)"),
        pl.col("pct hours negative").round(1),
        pl.col("normal pct negative").round(1),
        pl.col("p99.9").round(2).alias("p99.9 ($/MWh)"),
        pl.col("normal p99.9 ($/MWh)").round(2),
    )
)
html_table(tail_audit)

3 rows × 10 columns

hub	hours	mean ($/MWh)	skewness	excess kurtosis	min ($/MWh)	pct hours negative	normal pct negative	p99.9 ($/MWh)	normal p99.9 ($/MWh)
NP15	27,767	39.34	8.760	243.82	−38.47	3.40	6.90	244.00	121.24
SP15	27,767	33.95	7.420	189.66	−67.08	11.70	15.00	311.97	135.26
ZP26	27,767	33.14	6.040	164.14	−68.01	11.50	12.70	224.66	122.78

What the Tails Actually Look Like

• The right tail is fat. Observed p99.9 runs 244 $/MWh at NP15, 312 at SP15 and 225 at ZP26, roughly double the matched-normal benchmarks of 121 to 135, with excess kurtosis of 164 to 244: the day-ahead auction reprices scarcity rather than smoothing it.

• The below-zero shoulder is real but not the anomaly. Hours below zero are 3.4% at NP15, 11.7% at SP15 and 11.5% at ZP26, below the matched-normal shares of 6.9%, 15.0% and 12.7% because the spike-inflated standard deviation expects more negative mass than the market delivers.

• The left tail is wide. The worst hours are -38 $/MWh at NP15, -67 at SP15 and -68 at ZP26, near the benchmark’s 0.1st percentile.

4.3 Negative Prices and the Log-Return Floor

pricing.math.volatility.RealizedVolatilityEstimator measures the Margrabe spread volatility from daily log returns, so the power leg is floored at the CAISO definition’s negative_price_floor before the log. How often does the hourly below-zero shoulder survive daily averaging and make that floor bind? Daily mean DAM price per trading hub (US/Pacific days) against the floor, then the same check on gas.

POWER_PRICE_FLOOR = market.definition.negative_price_floor

hub_daily_mean = (
    dam.group_by("date", "hub")
    .agg(pl.col("lmp").mean().alias("daily_mean"))
    .sort("date", "hub")
)

fig = ts_figure(
    "Daily Mean DAM Price by Trading Hub (US/Pacific days)",
    "Daily mean DAM LMP ($/MWh)",
)
renderers = []
for color, hub in zip(SERIES, ALL_HUBS.values()):
    d = hub_daily_mean.filter(pl.col("hub") == hub)
    renderers.append(
        fig.line(d["date"].cast(pl.Datetime("ms")).to_numpy(),
                 d["daily_mean"].to_numpy(), color=color,
                 line_width=1.6, legend_label=hub, name=hub)
    )
# the floor the volatility estimator applies before taking the log
fig.add_layout(Span(location=POWER_PRICE_FLOOR, dimension="width",
                    line_color="#e34948", line_dash="dashed", line_width=1.5))
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("hub", "$name"), ("day", "$x{%F}"),
              ("mean", "$y{0.00} $/MWh")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

Where the Floor Binds

The dashed line is the floor. Exact per-hub counts, then the gas legs:

floor_audit = (
    hub_daily_mean.group_by("hub")
    .agg(
        pl.len().alias("days"),
        (pl.col("daily_mean") < POWER_PRICE_FLOOR).sum().alias("days below floor"),
        (pl.col("daily_mean") < 0.0).sum().alias("days negative"),
        pl.col("daily_mean").min().round(2).alias("min daily mean ($/MWh)"),
    )
    .sort("hub")
)
html_table(floor_audit)

5 rows × 5 columns

hub	days	days below floor	days negative	min daily mean ($/MWh)
NP15	1,157	4	2	−2.23
PACE	92	8	8	−10.67
PACW	92	0	0	4.97
SP15	1,157	11	11	−10.23
ZP26	1,157	9	8	−7.66

# The pairing the market defines, which is also what leaves the GHG variants out.
GAS_LEGS = {
    CAISO.gas_hub_for(hub).identifier: CAISO.gas_hub_for(hub).label
    for hub in CAISO.hub_labels
}
gas_floor_audit = (
    pl.scan_parquet(OS_DIR / "caiso_fuel_prices.parquet")
    .filter(pl.col("fuel_region_id").is_in(list(GAS_LEGS)))
    .with_columns(PACIFIC.dt.date().alias("date"))
    .filter(pl.col("date").is_between(*DATE_RANGE))
    .group_by("fuel_region_id")
    .agg(
        pl.len().alias("rows"),
        (pl.col("price") <= 0.0).sum().alias("rows at or below zero"),
        pl.col("price").min().round(3).alias("min price ($/MMBtu)"),
    )
    .sort("fuel_region_id")
    .collect()
    .with_columns(pl.col("fuel_region_id").replace(GAS_LEGS).alias("gas leg"))
    .select("gas leg", "fuel_region_id", "rows", "rows at or below zero",
            "min price ($/MMBtu)")
)
html_table(gas_floor_audit)

4 rows × 5 columns

gas leg	fuel_region_id	rows	rows at or below zero	min price ($/MMBtu)
PacifiCorp East	FRPACE	27,765	0	0.760
PacifiCorp West	FRPACW	27,765	120	−0.490
PG&E Citygate	FRPGE1	27,765	0	2.480
SoCal Citygate	FRSCE1	27,765	0	5.050

The floor is close to inert for the three CAISO trading hubs - it binds on 4 of 1,157 Pacific days at NP15, 11 at SP15 and 9 at ZP26 - so the SP15 estimate feeding the Margrabe pricing is clean. Daily averaging erases the below-zero shoulder from 4.2: the evening ramp pays for the midday surplus within the same day.

PACE is different: eight of its 92 live days (about 9%) settle at a negative daily mean, so a log-return volatility there is largely an artifact of the floor, and PacifiCorp West gas prints at or below zero on about 0.4% of hours. PACE/PACW pricing needs an arithmetic (Bachelier-style) spread volatility rather than a bigger floor.

4.4 The Weather Layer: Geography and Coverage

The market data settles on trading hubs; the weather layer lives on two geographies of its own - market-zone indices built from weighted Open-Meteo city points, and NOAA climate-division proxies. (Notebooks 03 and 04 carry the station-level degree-day build.) The point assignments and weights below are analyst choices:

html_table(
    weather_crosswalk.select(
        "location",
        "region",
        pl.col("noaa_climate_division_proxy").alias("NOAA proxy zone"),
        pl.col("prototype_caiso_weight").alias("weight"),
        pl.col("latitude").round(3),
        pl.col("longitude").round(3),
    ),
    kind="reference",
)

reference · 7 rows × 6 columns

location	region	NOAA proxy zone	weight	latitude	longitude
sacramento	NP15	Sacramento Drainage	0.180	38.58	−121.49
san_jose	NP15	Central Coast Drainage	0.180	37.34	−121.89
los_angeles	SP15	South Coast Drainage	0.280	34.05	−118.24
san_diego	SP15	South Coast Drainage	0.130	32.72	−117.16
bakersfield	ZP26	San Joaquin Drainage	0.080	35.37	−119.02
fresno	ZP26	San Joaquin Drainage	0.100	36.74	−119.79
merced	ZP26	San Joaquin Drainage	0.050	37.30	−120.48

tac_areas = (
    pl.scan_parquet(OS_DIR / "caiso_load_hourly.parquet")
    .select(pl.col("tac_area_name").unique().sort())
    .collect()["tac_area_name"]
)
geography_layers = pl.DataFrame(
    {
        "layer": [
            "CAISO trading hubs",
            "TAC load areas",
            "market-zone weather index",
            "NOAA climate-division proxy",
        ],
        "distinct units": [
            dam["hub"].n_unique(),
            len(tac_areas),
            weather_daily["region"].n_unique(),
            weather_zone_daily["weather_zone"].n_unique(),
        ],
        "examples": [
            ", ".join(sorted(dam["hub"].unique())),
            ", ".join(tac_areas.head(4).to_list()) + ", ...",
            ", ".join(sorted(weather_daily["region"].unique())),
            ", ".join(sorted(weather_zone_daily["weather_zone"].unique())),
        ],
        "role": [
            "price and spread settlement",
            "lower-level load geography; needs a crosswalk before weather joins",
            "primary weather join to hub-level market data",
            "weather-exposure geography and basis-risk diagnostics",
        ],
    }
)
html_table(geography_layers, kind="reference")

reference · 4 rows × 4 columns

layer	distinct units	examples	role
CAISO trading hubs	5	NP15, PACE, PACW, SP15, ZP26	price and spread settlement
TAC load areas	38	AVA, AVRN, AZPS, BANC, ...	lower-level load geography; needs a crosswalk before weather joins
market-zone weather index	4	CAISO, NP15, SP15, ZP26	primary weather join to hub-level market data
NOAA climate-division proxy	4	Central Coast Drainage, Sacramento Drainage, San Joaquin Drainage, South Coast Drainage	weather-exposure geography and basis-risk diagnostics

basis_flags = basis_risk.filter(
    pl.col("basis_direction") == "market_zone_to_weather_zones"
).select(
    pl.col("market_region").alias("market zone"),
    pl.col("weather_zones").alias("weather zones inside"),
    pl.col("basis_risk_flag").alias("basis risk"),
    pl.col("basis_risk_note").alias("note"),
)
html_table(basis_flags, kind="reference")

reference · 4 rows × 4 columns

market zone	weather zones inside	basis risk	note
CAISO	Central Coast Drainage, Sacramento Drainage, San Joaquin Drainage, South Coast Drainage	True	All-CAISO index combines multiple market and weather-zone proxies; treat as a broad system exposure.
NP15	Central Coast Drainage, Sacramento Drainage	True	Market zone mixes multiple weather-zone proxies; include dispersion or interaction terms.
SP15	South Coast Drainage	False	Market zone maps to one weather-zone proxy in this prototype.
ZP26	San Joaquin Drainage	False	Market zone maps to one weather-zone proxy in this prototype.

NP15 is the flagged mismatch: one market zone blending two weather regimes, which notebook 09 (CAISO Basis Risk Exploration) treats in full. Coverage of the weather tables inside DATE_RANGE - the joins below stop where the weather build ends:

WEATHER_TABLES = {
    "caiso_weather_daily_region": weather_daily,
    "caiso_weather_hourly_region": weather_hourly,
    "caiso_weather_daily_weather_zone": weather_zone_daily,
}
weather_coverage = pl.concat(
    [
        frame.select(
            pl.lit(name).alias("table"),
            pl.len().alias("rows"),
            pl.col("date").min().alias("start"),
            pl.col("date").max().alias("end"),
        )
        for name, frame in WEATHER_TABLES.items()
    ]
)
html_table(weather_coverage)

3 rows × 4 columns

table	rows	start	end
caiso_weather_daily_region	4,504	2023-06-01	2026-06-30
caiso_weather_hourly_region	107,976	2023-06-01	2026-06-30
caiso_weather_daily_weather_zone	4,504	2023-06-01	2026-06-30

6. Quantitative Analysis & Validation

6.1 Location Spreads Between the Trading Hubs

Price at Hub A minus price at Hub B - pure congestion and losses. Each pair is a settled instrument off market on the trade-day reporting key.

HUB_PAIRS = (("NP15", "SP15"), ("ZP26", "SP15"), ("NP15", "ZP26"))
location_spreads = pl.concat(
    [
        market.location_spread(long_hub, short_hub)
        .by("trade_day")
        .results.frame.select(
            pl.lit(f"{long_hub}-{short_hub}").alias("pair"),
            pl.col("trade_day").alias("date"),
            pl.col("value").alias("spread"),
        )
        for long_hub, short_hub in HUB_PAIRS
    ]
).sort("pair", "date")

fig = ts_figure("Daily Mean DAM Location Spreads Between CAISO Trading Hubs",
                "Spread ($/MWh)")
fig.add_layout(Span(location=0, dimension="width", line_color="#c3c2b7", line_width=1))
renderers = []
for color, (long_hub, short_hub) in zip(SERIES, HUB_PAIRS):
    pair = f"{long_hub}-{short_hub}"
    d = location_spreads.filter(pl.col("pair") == pair)
    renderers.append(
        fig.line(d["date"].cast(pl.Datetime("ms")).to_numpy(), d["spread"].to_numpy(),
                 color=color, line_width=2, legend_label=pair, name=pair)
    )
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("pair", "$name"), ("day", "$x{%F}"),
              ("spread", "$y{0.00} $/MWh")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

What Do the Spreads Mean?

Positive means the first hub trades at a premium; near zero means transmission is unconstrained and prices equalize.

6.2 TB4 - The Daily Storage Spread for Batteries

Average of the 4 most expensive hours minus the average of the 4 cheapest hours of each US/Pacific trading day, per hub.

\[\text{TB4}_{h,d} = \frac{1}{4}\sum_{i \,\in\, \text{Top-4}(h,d)} P_{h,i} \;-\; \frac{1}{4}\sum_{i \,\in\, \text{Bottom-4}(h,d)} P_{h,i}\]

where \(P_{h,i}\) is the DAM hourly LMP at hub \(h\) in hour \(i\) of Pacific trade day \(d\). market.top_bottom settles this on a perfect-foresight engine.

tb4 = pl.concat(
    [
        market.top_bottom(hub, duration_hours=4)
        .by("trade_day")
        .results.frame.select(
            pl.col("point").alias("hub"),
            pl.col("trade_day").alias("date"),
            pl.col("value").alias("tb4"),
        )
        for hub in HUBS.values()
    ]
).sort("date")

fig = ts_figure("TB4 by Trading Hub (DAM, Top-4 minus Bottom-4 hours per Pacific day)",
                "TB4 ($/MWh)")
renderers = []
for color, hub in zip(SERIES, HUBS.values()):
    d = tb4.filter(pl.col("hub") == hub)
    renderers.append(
        fig.line(d["date"].cast(pl.Datetime("ms")).to_numpy(), d["tb4"].to_numpy(),
                 color=color, line_width=2, legend_label=hub, name=hub)
    )
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("hub", "$name"), ("day", "$x{%F}"),
              ("TB4", "$y{0.00} $/MWh")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

What the Chart Tells You

Baseline TB4 sits around 20-55 $/MWh at all three hubs (medians: 28 NP15, 41 SP15, 41 ZP26). The largest single-day spreads all land on 2023-08-16 - about 892 at SP15, 769 at NP15, 758 at ZP26 - a statewide scarcity day, so the extreme TB4 days are system events, not localized congestion. SP15 and ZP26 track each other closely (correlation 0.94) and share a seasonal shape; NP15 correlates with each at only 0.62.

6.3 Net Load - The Duck Curve

System load minus solar and wind: what the dispatchable fleet has to serve. The solar build-out is what shapes the other series here - midday surplus pushes prices toward and below zero, and the evening ramp concentrates scarcity into a few hours. Average profile by Pacific hour of day, four sample months.

\[L^{\text{net}}_h = L_h - S_h - W_h\]

where \(L_h\) is system load, \(S_h\) solar and \(W_h\) wind generation in hour \(h\).

load_hourly = (
    pl.scan_parquet(OS_DIR / "caiso_load.parquet")
    .group_by(pl.col("interval_start_utc").dt.truncate("1h"))
    .agg(pl.col("load").mean())
)
renewables = (
    pl.scan_parquet(OS_DIR / "caiso_renewables_hourly.parquet")
    .filter(pl.col("location") == "CAISO")
    .select("interval_start_utc", "solar", "wind")
)
SAMPLE_MONTHS = ["2025-08", "2025-12", "2026-04", "2026-07"]
net_profile = (
    load_hourly.join(renewables, on="interval_start_utc")
    .with_columns(
        (pl.col("load") - pl.col("solar") - pl.col("wind")).alias("net_load"),
        PACIFIC.dt.strftime("%Y-%m").alias("month"),
        PACIFIC.dt.hour().alias("hour"),
    )
    .filter(pl.col("month").is_in(SAMPLE_MONTHS))
    .group_by("month", "hour")
    .agg(pl.col("net_load").mean())
    .sort("month", "hour")
    .collect()
)

fig = figure(
    title="Average Net Load by Pacific Hour of Day",
    x_axis_label="Hour of day (US/Pacific)",
    y_axis_label="Net load (MW)",
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
for color, month in zip(SERIES, SAMPLE_MONTHS):
    d = net_profile.filter(pl.col("month") == month)
    fig.line(d["hour"].to_numpy(), d["net_load"].to_numpy(), color=color,
             line_width=2, legend_label=month)
fig.legend.location = "top_left"
show(style(fig))

Anatomy of the Plot

Solar through hours 9-15 absorbs almost all system load: April 2026 dips to about 1,400 MW at 1 PM, the hours that print negative RT prices and curtailment. The evening ramp is the neck of the duck - August 2025 and July 2026 climb roughly 24,000-25,000 MW between hour 10 and hour 19. December 2025 is flatter and double-peaked: about 22,100 MW at 6 AM and 26,100 MW at 5 PM.

6.4 Outages - Forced vs Planned

Daily mean of the derived hourly outage aggregate. Forced (unscheduled) outages are the quanto conditioning variable candidate: weather-driven generator failures coinciding with price spikes.

outages_daily = (
    pl.scan_parquet(OS_DIR / "caiso_outages_aggregated_hourly_derived.parquet")
    .with_columns(PACIFIC.dt.date().alias("date"))
    .filter(pl.col("date").is_between(*DATE_RANGE))
    .group_by("date")
    .agg(
        pl.col("forced_curtailment_mw").mean().alias("forced"),
        pl.col("planned_curtailment_mw").mean().alias("planned"),
    )
    .sort("date")
    .collect()
)

fig = ts_figure("Daily Mean Outage MW, Forced vs Planned", "Outage (MW)")
x = outages_daily["date"].cast(pl.Datetime("ms")).to_numpy()
renderers = [
    fig.line(x, outages_daily[name].to_numpy(), color=color, line_width=2,
             legend_label=name, name=name)
    for color, name in zip(SERIES, ["forced", "planned"])
]
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("type", "$name"), ("day", "$x{%F}"),
              ("MW", "$y{0,0}")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

CAISO’s Curtailed and Non-Operational Generator (CNOG) report files ambient temperature derates, hydro head-loss and air density limits as “FORCED”, so the series measures available nameplate capacity reduction, not unit reliability failures. Against roughly 107 GW of nameplate across about 1,500 generators, 15-25 GW of forced derates is 15-20% of the fleet.

6.5 The Gas Leg - Citygate Fuel Regions

Daily mean gas price in the two main index regions

GAS_REGIONS = {
    CAISO.gas_hub_for(hub).identifier: CAISO.gas_hub_for(hub).label
    for hub in HUBS.values()
}
gas_daily = (
    pl.scan_parquet(OS_DIR / "caiso_fuel_prices.parquet")
    .filter(pl.col("fuel_region_id").is_in(list(GAS_REGIONS)))
    .with_columns(PACIFIC.dt.date().alias("date"))
    .filter(pl.col("date").is_between(*DATE_RANGE))
    .group_by("date", "fuel_region_id")
    .agg(pl.col("price").mean())
    .sort("date")
    .collect()
)

fig = ts_figure("Daily Mean Citygate Gas Prices by Fuel Region", "Gas price ($/MMBtu)")
renderers = []
for color, (region, label) in zip(SERIES, GAS_REGIONS.items()):
    d = gas_daily.filter(pl.col("fuel_region_id") == region)
    renderers.append(
        fig.line(d["date"].cast(pl.Datetime("ms")).to_numpy(), d["price"].to_numpy(),
                 color=color, line_width=2, legend_label=label, name=label)
    )
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("region", "$name"), ("day", "$x{%F}"),
              ("price", "$y{0.000} $/MMBtu")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

SoCal Citygate trades at a premium over PG&E Citygate on 99.4% of days in the sample, usually 1.80 to 3.20 $/MMBtu higher (median 2.51), with a handful of inversion days reaching -1.80 at the deepest. Southern California produces almost no gas locally and depends on long-distance interstate pipelines.

6.6 Zonal Weather - Cooling Stress on Two Geographies

The same cooling-degree signal on both geographies: 14-day rolling mean of daily CDD65 by market zone (left) and by NOAA climate-division proxy (right). If the market zones were weather-homogeneous, each left-panel curve would match one on the right.

zone_cdd = (
    weather_daily.filter(pl.col("region").is_in(list(HUBS.values())))
    .sort("region", "date")
    .with_columns(
        pl.col("cdd65_from_weighted_temp")
        .rolling_mean(14, min_samples=1)
        .over("region")
        .alias("cdd65_14d")
    )
)
proxy_cdd = weather_zone_daily.sort("weather_zone", "date").with_columns(
    pl.col("cdd65_from_weighted_temp")
    .rolling_mean(14, min_samples=1)
    .over("weather_zone")
    .alias("cdd65_14d")
)

zone_fig = ts_figure("14-Day Rolling CDD65 by Market Zone", "CDD65 (14-day mean)")
proxy_fig = ts_figure("14-Day Rolling CDD65 by NOAA Proxy Zone", "CDD65 (14-day mean)")
for fig, frame, key in ((zone_fig, zone_cdd, "region"),
                        (proxy_fig, proxy_cdd, "weather_zone")):
    renderers = []
    for color, name in zip(SERIES, frame[key].unique().sort()):
        d = frame.filter(pl.col(key) == name)
        renderers.append(
            fig.line(d["date"].cast(pl.Datetime("ms")).to_numpy(),
                     d["cdd65_14d"].to_numpy(), color=color, line_width=2,
                     legend_label=name, name=name)
        )
    fig.add_tools(
        HoverTool(renderers=renderers, tooltips=[("zone", "$name"), ("day", "$x{%F}"),
                  ("CDD65", "$y{0.00}")], formatters={"$x": "datetime"})
    )
    fig.legend.location = "top_left"
show(row(style(zone_fig), style(proxy_fig), sizing_mode="stretch_width"))

def zone_weather_summary(frame: pl.DataFrame, key: str, layer: str) -> pl.DataFrame:
    return (
        frame.group_by(pl.col(key).alias("zone"))
        .agg(
            pl.len().alias("days"),
            pl.col("temperature_f_mean_weighted").mean().round(1).alias("mean temp (F)"),
            pl.col("temperature_f_std").mean().round(1).alias("mean dispersion (F)"),
            pl.col("cdd65_from_weighted_temp").mean().round(2).alias("mean CDD65"),
            pl.col("cdd65_from_weighted_temp").max().round(1).alias("max CDD65"),
            pl.col("hdd65_from_weighted_temp").mean().round(2).alias("mean HDD65"),
            pl.col("hdd65_from_weighted_temp").max().round(1).alias("max HDD65"),
        )
        .with_columns(pl.lit(layer).alias("layer"))
        .sort("zone")
    )


html_table(
    pl.concat(
        [
            zone_weather_summary(weather_daily, "region", "market zone"),
            zone_weather_summary(weather_zone_daily, "weather_zone", "NOAA proxy"),
        ]
    ).select("layer", "zone", "days", "mean temp (F)", "mean dispersion (F)",
             "mean CDD65", "max CDD65", "mean HDD65", "max HDD65")
)

8 rows × 9 columns

layer	zone	days	mean temp (F)	mean dispersion (F)	mean CDD65	max CDD65	mean HDD65	max HDD65
market zone	CAISO	1,126	64.30	4.600	3.890	20.50	4.600	21.30
market zone	NP15	1,126	62.50	2.200	3.250	23.80	5.730	22.90
market zone	SP15	1,126	64.10	1.300	2.490	23.90	3.420	18.30
market zone	ZP26	1,126	67.50	2.000	7.890	34.20	5.420	25.40
NOAA proxy	Central Coast Drainage	1,126	61.50	0.000	1.940	21.80	5.450	22.70
NOAA proxy	Sacramento Drainage	1,126	63.60	0.000	4.720	27.40	6.160	25.40
NOAA proxy	San Joaquin Drainage	1,126	67.20	2.000	7.670	33.90	5.480	25.20
NOAA proxy	South Coast Drainage	1,126	63.70	1.300	2.180	21.10	3.440	17.90

Where the Weather Geographies Disagree

• ZP26 and its San Joaquin Drainage counterpart carry the strongest cooling signal (mean CDD65 7.86 and 7.64) - the Central Valley heat the generator map in notebook 07 shows plant by plant.

• NP15 is the mismatch flagged in 4.4: its 3.21 mean CDD65 averages Sacramento Drainage heat (4.69) against Central Coast marine cooling (1.89), two right-panel curves that separate all summer, and it carries the widest dispersion of the three hub zones.

• The market-zone index therefore understates the exposure of anything sited inland in NP15; notebook 09 prices that gap.

6.7 Weather Against Load and Outages

The quanto conditioning story needs weather stress to coincide with market stress. The CAISO-wide daily weather index against system load, generation and the outage aggregates, on Pacific days. The joins are inner, so the sample ends where the weather build ends.

load_daily = (
    pl.scan_parquet(OS_DIR / "caiso_load.parquet")
    .with_columns(PACIFIC.dt.date().alias("date"))
    .group_by("date")
    .agg(pl.col("load").mean().alias("system_load_mw"))
)
generation_daily = (
    pl.scan_parquet(OS_DIR / "caiso_fuel_mix.parquet")
    .with_columns(PACIFIC.dt.date().alias("date"))
    .group_by("date")
    .agg(
        pl.col("natural_gas").mean().alias("gas_generation_mw"),
        pl.col("solar").mean().alias("solar_mw"),
        pl.col("wind").mean().alias("wind_mw"),
    )
)
outages_daily_mean = (
    pl.scan_parquet(OS_DIR / "caiso_outages_aggregated_hourly_derived.parquet")
    .with_columns(PACIFIC.dt.date().alias("date"))
    .group_by("date")
    .agg(
        pl.col("forced_curtailment_mw").mean().alias("forced_mw"),
        pl.col("planned_curtailment_mw").mean().alias("planned_mw"),
    )
)
coincidence_daily = (
    weather_daily.filter(pl.col("region") == "CAISO")
    .select("date", "temperature_f_mean_weighted", "temperature_f_std",
            "cdd65_from_weighted_temp", "hdd65_from_weighted_temp")
    .lazy()
    .join(load_daily, on="date")
    .join(generation_daily, on="date")
    .join(outages_daily_mean, on="date")
    .collect()
    .drop_nulls()
    .sort("date")
)


def coincidence_scatter(y_col: str, title: str, y_label: str):
    """Scatter of one daily aggregate against CDD65, with its OLS line."""
    x = coincidence_daily["cdd65_from_weighted_temp"].to_numpy()
    y = coincidence_daily[y_col].to_numpy()
    fit = least_squares_weights({"cdd65": x}, y)
    fig = figure(title=title, x_axis_label="daily CDD65 (CAISO weighted index)",
                 y_axis_label=y_label, height=380, sizing_mode="stretch_width",
                 toolbar_location="above")
    fig.scatter(x, y, size=5, color=SERIES[0], alpha=0.45)
    grid = np.array([x.min(), x.max()])
    fig.line(grid, fit.weights["cdd65"] * grid + fit.weights["intercept"],
             color="#52514e", line_width=2, line_dash="dashed")
    return fig


show(row(
    style(coincidence_scatter("system_load_mw", "Daily CDD65 vs System Load",
                              "daily mean load (MW)")),
    style(coincidence_scatter("forced_mw", "Daily CDD65 vs Forced Outage MW",
                              "daily mean forced curtailment (MW)")),
    sizing_mode="stretch_width",
))

The dashed lines are OLS fits. The full pairwise picture:

COINCIDENCE_LABELS = {
    "CDD65": "cdd65_from_weighted_temp",
    "HDD65": "hdd65_from_weighted_temp",
    "dispersion": "temperature_f_std",
    "system load": "system_load_mw",
    "solar": "solar_mw",
    "wind": "wind_mw",
    "gas generation": "gas_generation_mw",
    "forced MW": "forced_mw",
    "planned MW": "planned_mw",
}
coincidence_matrix = np.corrcoef(
    coincidence_daily.select(list(COINCIDENCE_LABELS.values())).to_numpy(),
    rowvar=False,
)
corr_square = pl.DataFrame(
    {"series": list(COINCIDENCE_LABELS)}
    | {
        label: coincidence_matrix[:, i]
        for i, label in enumerate(COINCIDENCE_LABELS)
    }
)
corr_heatmap(corr_square, "series",
             "Daily Correlations: Weather Against the Market Aggregates")

Exact correlations against the two degree-day series, then the same outage relationship residualized on load, solar, wind, gas generation and temperature dispersion, which strips out the demand channel:

CORR_TARGETS = {"system load": "system_load_mw",
                "gas generation": "gas_generation_mw",
                "forced outage": "forced_mw",
                "planned outage": "planned_mw"}
weather_corr = (
    coincidence_daily.select(
        [pl.corr("cdd65_from_weighted_temp", col).round(3).alias(label)
         for label, col in CORR_TARGETS.items()]
    )
    .unpivot(variable_name="daily series", value_name="corr with CDD65")
    .join(
        coincidence_daily.select(
            [pl.corr("hdd65_from_weighted_temp", col).round(3).alias(label)
             for label, col in CORR_TARGETS.items()]
        ).unpivot(variable_name="daily series", value_name="corr with HDD65"),
        on="daily series",
    )
    .with_columns(pl.lit(coincidence_daily.height).alias("days"))
)
html_table(weather_corr)

4 rows × 4 columns

daily series	corr with CDD65	corr with HDD65	days
system load	0.867	−0.435	1,125
gas generation	0.459	−0.048	1,125
forced outage	0.046	−0.044	1,125
planned outage	−0.592	0.340	1,125

CONTROLS = ["system_load_mw", "solar_mw", "wind_mw", "gas_generation_mw",
            "temperature_f_std"]
controls = {name: coincidence_daily[name].to_numpy() for name in CONTROLS}
forced_mw = coincidence_daily["forced_mw"].to_numpy()
resid_audit = pl.DataFrame(
    {
        "relationship": ["CDD65 vs forced outage MW", "HDD65 vs forced outage MW"],
        "raw corr": [
            coincidence_daily.select(
                pl.corr("cdd65_from_weighted_temp", "forced_mw")
            ).item(),
            coincidence_daily.select(
                pl.corr("hdd65_from_weighted_temp", "forced_mw")
            ).item(),
        ],
        "residualized corr": [
            partial_correlation(
                coincidence_daily["cdd65_from_weighted_temp"].to_numpy(),
                forced_mw,
                controls,
            ),
            partial_correlation(
                coincidence_daily["hdd65_from_weighted_temp"].to_numpy(),
                forced_mw,
                controls,
            ),
        ],
    }
).with_columns(pl.col(pl.Float64).round(3), pl.lit(coincidence_daily.height).alias("days"))
html_table(resid_audit)

2 rows × 4 columns

relationship	raw corr	residualized corr	days
CDD65 vs forced outage MW	0.046	−0.094	1,125
HDD65 vs forced outage MW	−0.044	0.011	1,125

Cooling stress is a demand story in this aggregate, not an outage story: CDD65 moves with load (0.89) and the gas generation serving it (0.48), planned outages correlate negatively (-0.60) as maintenance concentrates in the mild shoulder seasons, and forced outages sit flat (+0.05). Residualizing leaves -0.09 for CDD65 and 0.00 for HDD65. The system file is also the wrong resolution for the outage question: reason codes go to notebook 02, outage-count models to notebooks 10-12, plant-level heat exposure to notebook 07.

6.8 The Quanto Covariance Term - Weather Against Hub Prices

A quanto’s value over pricing the legs independently is exactly the covariance term:

\[E[XY] = E[X]E[Y] + \text{Cov}(X,Y)\]

Measured hourly per trading hub: cooling degree hours (CDH65) from the market-zone weather index against the DAM and RT power legs and the DAM-RT spread (hourly DAM minus the hourly mean of RT 15-min prices).

rt_hourly = (
    pl.scan_parquet(OS_DIR / "caiso_lmp_real_time_15_min.parquet")
    .filter(pl.col("location").is_in(list(HUBS)))
    .with_columns(
        pl.col("location").replace(HUBS).alias("hub"),
        # the market parquets mix nanosecond/microsecond datetime units;
        # align every join key to the weather layer's microseconds
        pl.col("interval_start_utc").dt.truncate("1h").dt.cast_time_unit("us"),
    )
    .group_by("interval_start_utc", "hub")
    .agg(pl.col("lmp").mean().alias("rt_lmp"))
)
price_weather_hourly = (
    core.select(
        pl.col("interval_start_utc").dt.cast_time_unit("us"),
        "hub",
        pl.col("lmp").alias("da_lmp"),
    )
    .lazy()
    .join(rt_hourly, on=["interval_start_utc", "hub"])
    .join(
        weather_hourly.filter(pl.col("region").is_in(list(HUBS.values())))
        .select("interval_start_utc", pl.col("region").alias("hub"),
                "cdh65_from_weighted_temp", "temperature_f_std")
        .lazy(),
        on=["interval_start_utc", "hub"],
    )
    .with_columns((pl.col("da_lmp") - pl.col("rt_lmp")).alias("dam_rt_spread"))
    .collect()
)

hub_corr = (
    price_weather_hourly.group_by("hub")
    .agg(
        pl.corr("cdh65_from_weighted_temp", "rt_lmp").alias("CDH65 vs RT LMP"),
        pl.corr("cdh65_from_weighted_temp", "da_lmp").alias("CDH65 vs DAM LMP"),
        pl.corr("cdh65_from_weighted_temp", "dam_rt_spread").alias("CDH65 vs DAM-RT"),
        pl.corr("temperature_f_std", "dam_rt_spread").alias("dispersion vs DAM-RT"),
    )
    .sort("hub")
)
corr_long = hub_corr.unpivot(index="hub", variable_name="pair", value_name="correlation")
factors = list(corr_long.select("hub", "pair").iter_rows())
pair_color = dict(zip(hub_corr.columns[1:], SERIES))

fig = figure(
    y_range=FactorRange(*factors[::-1]),
    x_axis_label="hourly correlation",
    title="Weather-Price Coincidence by Trading Hub",
    height=420,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.hbar(y=factors, right=corr_long["correlation"].to_list(), height=0.7,
         color=[pair_color[p] for p in corr_long["pair"]])
fig.add_layout(Span(location=0, dimension="height", line_color="#c3c2b7", line_width=1))
show(style(fig))

The covariance decomposition for the stylized quanto pair - \(X\) the hourly CDH65, \(Y\) the hourly RT price. Inputs: the joined hourly prices and weather above. Operation: the mean product \(E[XY]\) less the product of means \(E[X]E[Y]\). Validation: the covariance share is the fraction of the naive product an independence assumption misses.

quanto_cov = (
    price_weather_hourly.group_by("hub")
    .agg(
        pl.len().alias("hours"),
        (pl.col("cdh65_from_weighted_temp") * pl.col("rt_lmp")).mean()
        .alias("mean product"),
        (pl.col("cdh65_from_weighted_temp").mean() * pl.col("rt_lmp").mean())
        .alias("product of means"),
    )
    .with_columns(
        (pl.col("mean product") - pl.col("product of means")).alias("covariance")
    )
    .with_columns(
        (100 * pl.col("covariance") / pl.col("mean product")).alias("covariance share (pct)")
    )
    .sort("hub")
    .with_columns(pl.col(pl.Float64).round(2))
)
html_table(quanto_cov)

3 rows × 6 columns

hub	hours	mean product	product of means	covariance	covariance share (pct)
NP15	26,993	172.86	165.38	7.48	4.33
SP15	26,993	155.12	110.36	44.75	28.85
ZP26	26,993	319.98	265.90	54.08	16.90

The covariance term is real and, on three years, uniformly positive: the coincidence adds about 15% to the naive product at NP15, 7% at SP15 and 11% at ZP26. One mild year told a duck-curve story in which solar-flooded afternoons made hourly heat and hourly price coincide negatively at SP15; the extended sample, which contains the 2023-24 scarcity era, overturns that sign - hot hours coincide with expensive hours at every hub once real spike days are in the average. Two handoffs follow. The quanto pricing in notebooks 13 and 14 still picks its settlement granularity deliberately: SettlementFrequency on the priced instrument names that choice, and the hourly and daily covariance terms remain different objects even when they agree in sign. The interaction terms this section suggests - CDD with load, forced outages, and the DAM-RT spread - are candidate regressors for notebooks 10-12.