04. Regional Temperature Mapping

04. Regional Temperature Mapping#

Authors: Brian Nguyen · Claire Kim · Created: 2026-07-23 · Tags: #Experiment #Quanto #Index-Construction #Weather

2. Objective & Input Framework#

Objective:

Build mapped regional temperatures $T^{NP15}_d$ and $T^{SP15}_d$ by blending each hub footprint’s stations on utility-load weights.
Measure the region-matched CDD-to-price dependence per trading hub.
Compare against notebook 03’s statewide-blend baseline and against a deliberately mismatched pairing.

Data & Feature Catalog:

Dataset	Granularity	Role
`openmeteo_temperature_daily`	Pacific day × 6 stations	city temperatures (SAC/SJC/FAT north; BUR/ONT/SAN south)
`caiso_load_hourly` (TAC areas)	hourly × TAC	utility-load weights (PGE-TAC, SCE-TAC, SDGE-TAC)
`caiso_lmp_real_time_15_min` (NP15 & SP15 hubs)	15-min	region-matched price legs

3. Mathematical Model & Analytical Formulation#

Mapped regional temperature for hub region $R$ with station set $S_R$:

\[T^{R}_d = \sum_{s \in S_R} w_s \, T_{s,d}, \qquad \sum_{s \in S_R} w_s = 1\]

Weights come from mean utility TAC load: each utility’s share of its region’s load, split equally across that utility’s stations (no city-level load exists):

\[w_s = \frac{\bar L_{u(s)}}{\sum_{u \in R} \bar L_u} \cdot \frac{1}{|S_{u(s)}|}\]

Regional degree days are computed on the mapped temperature, blending before the nonlinear payoff:

\[\text{CDD}^{R}_d = \max(0,\ T^{R}_d - 65)\]

That ordering is REGIONAL_CONVENTION (DegreeDayConvention.BLEND_THEN_TRANSFORM), with the floor applied by pricing.core.numeric.positive_part. Notebook 03’s statewide basket floors each station first (BASELINE_CONVENTION, rebuilt in section 6.2), a different object since $\max(0, \sum w T - 65) \ne \sum w \max(0, T - 65)$.

Variable Definition Table

Symbol	Meaning	Type / Units	Bounds	Code name
$T_{s,d}$	station settlement-mean temperature	Float64, °F	-	`temp_mean_f`
$65$	degree-day base temperature	Float64, °F	market constant	`DEGREE_DAY_BASE_FAHRENHEIT`
$\bar L_u$	mean hourly load of utility TAC $u$	Float64, MW	$> 0$	`tac_means`
$w_s$	station weight within its region	Float64	$[0,1]$, sums to 1	`WEIGHTS`
$T^{R}_d$	mapped regional temperature	Float64, °F	-	`t_np15`, `t_sp15`
$\text{CDD}^{R}_d$	regional cooling degree days	Float64, °F-days	$\ge 0$	`cdd_np15`, `cdd_sp15`
$P^{R}_d$	daily max RT 15-min LMP at hub $R$	Float64, $/MWh	can be $< 0$	`rt_np15`, `rt_sp15`

import polars as pl
from bokeh.io import output_notebook, show
from bokeh.models import HoverTool
from bokeh.plotting import figure
from bokeh.resources import INLINE

from chart_theme import SERIES, corr_heatmap, html_table, style
from markets.definitions.caiso import CAISO_CLOCK
from markets.definitions.caiso import CAISO, DEGREE_DAY_BASE_FAHRENHEIT
from pricing.core.numeric import guarded_ratio, positive_part
from pricing.indices.weather import DegreeDayConvention
from settings import config

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

# The regional index floors the blended temperature, which is the library
# default; notebook 03's basket of station indices takes the other
# convention, and section 6.2 rebuilds it here as the baseline.
REGIONAL_CONVENTION = DegreeDayConvention.BLEND_THEN_TRANSFORM
BASELINE_CONVENTION = DegreeDayConvention.TRANSFORM_THEN_BLEND

# Region -> station -> utility TAC. FAT (Fresno) sits at the Path 15
# boundary; it is assigned north with its PG&E service territory.
REGION_STATIONS = {
    "NP15": {"SAC": "PGE-TAC", "SJC": "PGE-TAC", "FAT": "PGE-TAC"},
    "SP15": {"BUR": "SCE-TAC", "ONT": "SCE-TAC", "SAN": "SDGE-TAC"},
}
# The APNode ids come from the market definition, not a retyped literal.
HUB_APNODES = {label: CAISO.hub(label).identifier for label in ("NP15", "SP15")}

# 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"))

temp = (
    pl.read_parquet(OS_DIR / "openmeteo_temperature_daily.parquet")
    .with_columns(PACIFIC.dt.date().alias("date"))
)

output_notebook(INLINE, hide_banner=True)


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 Microclimate Evidence#

“Region” means three things in CAISO - the trading hub (NP15/SP15), the utility load area (PGE/SCE/SDGE TACs), and the microclimate. A 90°F day in coastal San Jose stresses the north very differently than one in inland Sacramento. Do the added stations carry information the originals don’t?

micro = (
    temp.group_by("station")
    .agg(
        pl.len().alias("days"),
        pl.col("temp_mean_f").mean().round(1).alias("mean T (°F)"),
        pl.col("temp_mean_f").filter(PACIFIC.dt.month().is_in([6, 7, 8, 9]))
        .mean().round(1).alias("mean summer T (°F)"),
        positive_part(pl.col("temp_mean_f") - DEGREE_DAY_BASE_FAHRENHEIT)
        .mean().round(2).alias("mean CDD"),
    )
    .sort("station")
)
html_table(micro)

6 rows × 5 columns

station	days	mean T (°F)	mean summer T (°F)	mean CDD
BUR	3,160	66.00	75.80	4.690
FAT	3,160	66.10	83.40	7.280
ONT	3,152	67.40	78.70	5.940
SAC	3,160	63.40	76.60	4.650
SAN	3,160	63.00	69.80	1.930
SJC	3,152	61.80	71.00	2.490

wide_t = temp.pivot(on="station", index="date", values="temp_mean_f").sort("date")
stations = [s for r in REGION_STATIONS.values() for s in r]
corr_matrix = pl.DataFrame(
    {
        "station": stations,
        **{
            other: [
                wide_t.select(pl.corr(a, other)).item() for a in stations
            ]
            for other in stations
        },
    }
)
corr_heatmap(corr_matrix, "station",
             "Cross-Station Correlation of Daily Settlement-Mean Temperature")

Cross-Station Correlation of Daily Settlement-Mean Temperature

4.2 The Utility Load Weights#

tac_means = (
    pl.scan_parquet(OS_DIR / "caiso_load_hourly.parquet")
    .filter(pl.col("tac_area_name").is_in(["PGE-TAC", "SCE-TAC", "SDGE-TAC"]))
    .group_by("tac_area_name")
    .agg(pl.col("load").mean().round(0).alias("mean load (MW)"))
    .sort("mean load (MW)", descending=True)
    .collect()
)
html_table(tac_means)

3 rows × 2 columns

tac_area_name	mean load (MW)
SCE-TAC	11,723.00
PGE-TAC	11,342.00
SDGE-TAC	2,218.00

5. Theoretical Assumptions & Limitations#

Hub footprints and utility TAC areas do not coincide exactly (the NP15/SP15 split is at Path 15); Fresno is assigned north with PG&E, the within-utility split is equal across each utility’s stations, and SMUD load sits outside the PGE-TAC total. Joins are inner, so the panel ends at the shortest station.
Inherited from notebook 03: ERA5 reanalysis rather than settlement feeds, DST day lengths, and a price leg with skew ≈ 6.7 and excess kurtosis ≈ 56, so tail conditioning is mandatory. The weather panel spans about 8.6 years with multiple full winters.

6. Implementation & Simulation Logic#

6.1 Weights and Mapped Temperatures#

tac_lookup = dict(tac_means.iter_rows())
WEIGHTS = {}
for region, station_map in REGION_STATIONS.items():
    region_total = sum(tac_lookup[tac] for tac in set(station_map.values()))
    utility_counts = {tac: list(station_map.values()).count(tac) for tac in station_map.values()}
    WEIGHTS[region] = {
        s: tac_lookup[tac] / region_total / utility_counts[tac]
        for s, tac in station_map.items()
    }
weights_table = pl.DataFrame(
    [
        {"region": region, "station": s, "weight": round(w, 4)}
        for region, ws in WEIGHTS.items()
        for s, w in ws.items()
    ]
)
html_table(weights_table)

6 rows × 3 columns

region	station	weight
NP15	SAC	0.333
NP15	SJC	0.333
NP15	FAT	0.333
SP15	BUR	0.420
SP15	ONT	0.420
SP15	SAN	0.159

mapped = (
    wide_t.drop_nulls()
    .with_columns(
        t_np15=sum(pl.col(s) * w for s, w in WEIGHTS["NP15"].items()),
        t_sp15=sum(pl.col(s) * w for s, w in WEIGHTS["SP15"].items()),
    )
    .with_columns(
        cdd_np15=positive_part(pl.col("t_np15") - DEGREE_DAY_BASE_FAHRENHEIT),
        cdd_sp15=positive_part(pl.col("t_sp15") - DEGREE_DAY_BASE_FAHRENHEIT),
    )
    .select("date", "t_np15", "t_sp15", "cdd_np15", "cdd_sp15")
)

fig = ts_figure("Mapped Regional Temperature (load-weighted station blend)", "Temperature (°F)")
x = mapped["date"].to_numpy()
renderers = [
    fig.line(x, mapped[colname].to_numpy(), color=color, line_width=1.8,
             legend_label=label, name=label)
    for color, (colname, label) in zip(
        SERIES, {"t_np15": "NP15 mapped T", "t_sp15": "SP15 mapped T"}.items()
    )
]
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("region", "$name"), ("day", "$x{%F}"),
              ("T", "$y{0.0} °F")], formatters={"$x": "datetime"})
)
fig.legend.location = "bottom_left"
show(style(fig))

6.2 Region-Matched Price Legs and the Notebook-03 Baseline#

rt_hub = (
    pl.scan_parquet(OS_DIR / "caiso_lmp_real_time_15_min.parquet")
    .filter(pl.col("location").is_in(list(HUB_APNODES.values())))
    .group_by(PACIFIC.dt.date().alias("date"), "location")
    .agg(pl.col("lmp").max())
    .collect()
    .pivot(on="location", index="date", values="lmp")
    .rename({HUB_APNODES["NP15"]: "rt_np15", HUB_APNODES["SP15"]: "rt_sp15"})
)

# Notebook 03's statewide baseline: equal-weighted mean of the four
# original stations' CDDs, which is the BASELINE_CONVENTION ordering
baseline = (
    temp.filter(pl.col("station").is_in(["SAC", "BUR", "SAN", "FAT"]))
    .with_columns(
        cdd=positive_part(pl.col("temp_mean_f") - DEGREE_DAY_BASE_FAHRENHEIT)
    )
    .group_by("date")
    .agg(pl.col("cdd").mean().alias("cdd_statewide"))
)

market_data = mapped.join(rt_hub, on="date").join(baseline, on="date").sort("date")

7. Quantitative Analysis & Validation#

7.1 Performance Evaluation Matrix#

Each pairing is scored on full-sample correlation and on tail uplift: mean daily RT max on top-decile CDD days over the mean on all other days.

PAIRINGS = {
    "NP15 mapped CDD → NP15 price (matched)":
        ("cdd_np15", "rt_np15", REGIONAL_CONVENTION),
    "SP15 mapped CDD → SP15 price (matched)":
        ("cdd_sp15", "rt_sp15", REGIONAL_CONVENTION),
    "Statewide CDD → SP15 price (nb03 baseline)":
        ("cdd_statewide", "rt_sp15", BASELINE_CONVENTION),
    "NP15 mapped CDD → SP15 price (mismatched)":
        ("cdd_np15", "rt_sp15", REGIONAL_CONVENTION),
}


def score_pairing(
    label: str, cdd_col: str, price_col: str, convention: DegreeDayConvention
) -> dict:
    """Correlation and tail uplift of one CDD index against one price leg."""
    q90 = market_data[cdd_col].quantile(0.9)
    tail = market_data.filter(pl.col(cdd_col) >= q90)[price_col].mean()
    rest = market_data.filter(pl.col(cdd_col) < q90)[price_col].mean()
    return {
        "pairing": label,
        "convention": str(convention),
        "corr": round(market_data.select(pl.corr(cdd_col, price_col)).item(), 3),
        "mean RT max, top-decile CDD": round(tail, 1),
        "mean RT max, other days": round(rest, 1),
        "tail uplift": round(guarded_ratio(tail, rest), 2),
    }


evaluation = pl.DataFrame(
    [score_pairing(label, *terms) for label, terms in PAIRINGS.items()]
)
html_table(evaluation)

4 rows × 6 columns

pairing	convention	corr	mean RT max, top-decile CDD	mean RT max, other days	tail uplift
NP15 mapped CDD → NP15 price (matched)	blend_then_transform	0.320	247.40	79.50	3.110
SP15 mapped CDD → SP15 price (matched)	blend_then_transform	0.361	271.10	77.40	3.500
Statewide CDD → SP15 price (nb03 baseline)	transform_then_blend	0.345	278.40	76.60	3.640
NP15 mapped CDD → SP15 price (mismatched)	blend_then_transform	0.312	226.00	82.40	2.740

fig = figure(
    title="Daily Max RT LMP vs Mapped Regional CDD (region-matched pairs)",
    x_axis_label="Mapped regional CDD (°F-days)",
    y_axis_label="RT 15-min LMP daily max ($/MWh)",
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.scatter(market_data["cdd_np15"].to_numpy(), market_data["rt_np15"].to_numpy(),
            size=6, color=SERIES[0], alpha=0.45, legend_label="NP15")
fig.scatter(market_data["cdd_sp15"].to_numpy(), market_data["rt_sp15"].to_numpy(),
            size=6, color=SERIES[1], alpha=0.45, legend_label="SP15")
fig.legend.location = "top_left"
show(style(fig))

The matched pairings beat the mismatched pairing on both metrics, but the notebook-03 statewide baseline posts the highest tail uplift of the four, and only SP15’s matched correlation is above it. The mapping’s case rests on interpretability - each hub priced off its own region’s weather - not on measured uplift.

The construction has graduated to markets.adapters.caiso.weather.CaisoWeatherSource and pricing.indices.weather.WeatherBlendIndex under REGIONAL_CONVENTION, on a seven-city profile with market-wide load weights; this notebook stays the six-station derivation. The quanto payoff $\max(0, \text{CDD}^R - K_w) \times \max(0, P^R - K_p)$ prices on the matched pairs in notebooks 13 and 14. Queued: city-level load weights, a ZP26/valley treatment, and scoring the library’s blend on this matrix.

Categories

Tags