03. Weather - Station Degree-Day Indices

03. Weather - Station Degree-Day Indices#

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

2. Objective & Input Framework#

Objective:

Build daily and monthly HDD/CDD indices at the contract-relevant California stations, on the settlement convention.
Form a CAISO composite index.
Measure the covariance triple - CDD against price, CDD against temperature-derate outages, and derate against price.

Data & Feature Catalog:

Dataset	Granularity	Role
`openmeteo_temperature_daily`	Pacific day × 4 stations	temperature leg (Sacramento & Burbank are CME contract cities; San Diego & Fresno complete load-center coverage)
`caiso_lmp_real_time_15_min` (SP15 hub)	15-min	price leg
`caiso_curtailed_non_operational_generator_report`	outage segment × edition	temp-derate MW (notebook 02’s validated recipe)

3. Mathematical Model & Analytical Formulation#

Settlement mean and degree days per station $s$ and Pacific day $d$ (base 65°F):

\[T_{s,d} = \frac{T^{max}_{s,d} + T^{min}_{s,d}}{2}, \qquad \text{HDD}_{s,d} = \max(0,\ 65 - T_{s,d}), \qquad \text{CDD}_{s,d} = \max(0,\ T_{s,d} - 65)\]

The base is DEGREE_DAY_BASE_FAHRENHEIT in markets.definitions.caiso; the floor is the library’s pricing.core.numeric.positive_part.

Contracts settle on the cumulative monthly sum; the composite is an equal-weighted basket (an index of indices):

\[M_{s,m} = \sum_{d \in m} \text{CDD}_{s,d}, \qquad \bar C_d = \frac{1}{N}\sum_{s} \text{CDD}_{s,d}\]

The floor sits inside the basket: each station index is floored before averaging (TRANSFORM_THEN_BLEND), while a regional index floors the blended temperature (BLEND_THEN_TRANSFORM, the library default). Since $\max(0, \sum w T - 65) \ne \sum w \max(0, T - 65)$, Jensen puts the basket at or above the blended alternative; section 5.3 measures the gap, and notebook 04’s regional index takes the other convention.

Variable Definition Table

Symbol	Meaning	Type / Units	Bounds	Code name
$T^{max}_{s,d}, T^{min}_{s,d}$	daily max / min 2 m temperature	Float64, °F	plausible CA range	`temp_max_f`, `temp_min_f`
$T_{s,d}$	settlement mean	Float64, °F	-	`temp_mean_f`
$65$	degree-day base temperature	Float64, °F	market constant	`DEGREE_DAY_BASE_FAHRENHEIT`
$\text{HDD}_{s,d}, \text{CDD}_{s,d}$	daily degree days	Float64, °F-days	$\ge 0$	`hdd`, `cdd`
$M_{s,m}$	monthly cumulative CDD	Float64, °F-days	$\ge 0$	`monthly` `cdd`
$\bar C_d$	composite daily CDD (equal weights)	Float64, °F-days	$\ge 0$	`composite` `cdd`
$P_d$	daily max SP15 RT 15-min LMP	Float64, $/MWh	can be $< 0$	`rt_max`
$W_d$	daily mean temp-derate forced MW	Float64, MW	$\ge 0$	`temp_derate_mw`
$K_w$	tail threshold on $\bar C_d$	Float64, °F-days	90th pctile	`q90`

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, html_table, style
from markets.definitions.caiso import CAISO_CLOCK
from markets.definitions.caiso import DEGREE_DAY_BASE_FAHRENHEIT
from pricing.core.numeric import positive_part
from pricing.indices.weather import DegreeDayConvention
from settings import config

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

# A basket of separately settled station indices floors each station
# before averaging, so it takes the convention the library does not
# default to. Section 6.3 measures what the ordering costs.
COMPOSITE_CONVENTION = DegreeDayConvention.TRANSFORM_THEN_BLEND

STATION_LABELS = {
    "SAC": "Sacramento (CME)",
    "BUR": "Burbank (CME)",
    "SAN": "San Diego",
    "FAT": "Fresno",
}

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

# The cache also carries the notebook-04 mapping stations (SJC, ONT);
# this notebook's settlement indices use the original four only.
temp = (
    pl.read_parquet(OS_DIR / "openmeteo_temperature_daily.parquet")
    .filter(pl.col("station").is_in(list(STATION_LABELS)))
    .with_columns(PACIFIC.dt.date().alias("date"))
    .with_columns(
        hdd=positive_part(DEGREE_DAY_BASE_FAHRENHEIT - pl.col("temp_mean_f")),
        cdd=positive_part(pl.col("temp_mean_f") - DEGREE_DAY_BASE_FAHRENHEIT),
    )
)

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 The Station Temperature Data#

Open-Meteo ERA5 reanalysis at station coordinates, lagging a few days behind present, not the official station feed contracts settle on. The cache spans about 8.6 years per station; across the four stations total HDD is 63,160 °F-days against 58,292 CDD. The market legs cover 2023-06-01 through 2026-07-31, all on Pacific calendar days.

summary = (
    temp.group_by("station")
    .agg(
        pl.len().alias("days"),
        pl.col("temp_mean_f").mean().round(1).alias("mean temp (°F)"),
        pl.col("temp_max_f").max().alias("hottest day max (°F)"),
        pl.col("temp_min_f").min().alias("coldest day min (°F)"),
        pl.col("hdd").sum().round(0).alias("total HDD"),
        pl.col("cdd").sum().round(0).alias("total CDD"),
    )
    .sort("station")
)
html_table(summary)

4 rows × 7 columns

station	days	mean temp (°F)	hottest day max (°F)	coldest day min (°F)	total HDD	total CDD
BUR	3,160	66.00	115.20	32.90	11,600.00	14,831.00
FAT	3,160	66.10	117.40	26.60	19,418.00	23,018.00
SAC	3,160	63.40	117.40	27.00	19,774.00	14,687.00
SAN	3,160	63.00	97.80	35.30	12,368.00	6,095.00

quality = temp.null_count().transpose(
    include_header=True, header_name="column", column_names=["null rows"]
)
html_table(quality)

9 rows × 2 columns

column	null rows
interval_start_utc	0
interval_end_utc	0
station	0
temp_max_f	0
temp_min_f	0
temp_mean_f	0
date	0
hdd	0
cdd	0

4.2 Distribution Audit#

Station CDD is censored at zero; the price leg is the heavy-tail check.

rt_daily = (
    pl.scan_parquet(OS_DIR / "caiso_lmp_real_time_15_min.parquet")
    .filter(pl.col("location") == "TH_SP15_GEN-APND")
    .group_by(PACIFIC.dt.date().alias("date"))
    .agg(pl.col("lmp").max().alias("rt_max"))
    .collect()
)

dist_stats = pl.concat(
    [
        temp.group_by(pl.col("station").alias("series"))
        .agg(
            pl.col("cdd").mean().round(2).alias("mean"),
            pl.col("cdd").std().round(2).alias("std"),
            pl.col("cdd").skew().round(2).alias("skew"),
            pl.col("cdd").kurtosis().round(2).alias("kurtosis"),
            pl.col("cdd").max().round(1).alias("max"),
        )
        .sort("series"),
        rt_daily.select(
            pl.lit("SP15 RT daily max").alias("series"),
            pl.col("rt_max").mean().round(2).alias("mean"),
            pl.col("rt_max").std().round(2).alias("std"),
            pl.col("rt_max").skew().round(2).alias("skew"),
            pl.col("rt_max").kurtosis().round(2).alias("kurtosis"),
            pl.col("rt_max").max().round(1).alias("max"),
        ),
    ]
)
html_table(dist_stats)

5 rows × 6 columns

series	mean	std	skew	kurtosis	max
BUR	4.69	6.30	1.260	0.76	33.30
FAT	7.28	9.28	0.900	−0.64	35.10
SAC	4.65	6.36	1.230	0.65	31.80
SAN	1.93	3.16	1.720	2.43	19.20
SP15 RT daily max	96.62	151.25	6.670	55.57	2,000.00

4.3 The Temperature Series#

fig = ts_figure("Daily Settlement-Mean Temperature by Station", "Temperature (°F)")
renderers = []
for color, (station, label) in zip(SERIES, STATION_LABELS.items()):
    d = temp.filter(pl.col("station") == station).sort("date")
    renderers.append(
        fig.line(d["date"].to_numpy(), d["temp_mean_f"].to_numpy(), color=color,
                 line_width=1.5, legend_label=label, name=label)
    )
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("station", "$name"), ("day", "$x{%F}"),
              ("T", "$y{0.0} °F")], formatters={"$x": "datetime"})
)
fig.legend.location = "bottom_left"
show(style(fig))

5. Implementation & Simulation Logic#

5.1 Daily Degree Days#

Each daily CDD value is already an option-style payoff on temperature.

fig = ts_figure("Daily CDD by Station", "CDD (°F-days)")
renderers = []
for color, (station, label) in zip(SERIES, STATION_LABELS.items()):
    d = temp.filter(pl.col("station") == station).sort("date")
    renderers.append(
        fig.line(d["date"].to_numpy(), d["cdd"].to_numpy(), color=color,
                 line_width=1.5, legend_label=label, name=label)
    )
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("station", "$name"), ("day", "$x{%F}"),
              ("CDD", "$y{0.0}")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

5.2 Monthly Cumulative Indices - The Contract Settlement Object#

monthly = (
    temp.group_by(pl.col("date").dt.truncate("1mo").alias("month"), "station")
    .agg(pl.col("cdd").sum().round(1).alias("cdd"), pl.col("hdd").sum().round(1).alias("hdd"))
    .sort("month")
)

fig = ts_figure("Monthly CDD by Station (Contract Settlement View)", "CDD (°F-days)")
renderers = []
for color, (station, label) in zip(SERIES, STATION_LABELS.items()):
    d = monthly.filter(pl.col("station") == station)
    x = d["month"].to_numpy()
    r = fig.line(x, d["cdd"].to_numpy(), color=color, line_width=2,
                 legend_label=label, name=label)
    fig.scatter(x, d["cdd"].to_numpy(), color=color, size=5)
    renderers.append(r)
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("station", "$name"), ("month", "$x{%Y-%m}"),
              ("CDD", "$y{0.0}")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

5.3 The CAISO Composite Index#

Equal-weighted across the four stations - a placeholder until TAC-area load weighting - with Sacramento kept as the CME-settled NorCal proxy although SMUD sits outside CAISO. The dashed series is the BLEND_THEN_TRANSFORM alternative on the same stations and days.

composite = (
    temp.group_by("date")
    .agg(
        pl.col("cdd").mean().alias("cdd"),
        pl.col("hdd").mean().alias("hdd"),
        positive_part(
            pl.col("temp_mean_f").mean() - DEGREE_DAY_BASE_FAHRENHEIT
        ).alias("cdd_blended_temperature"),
    )
    .sort("date")
)
COMPOSITE_SERIES = {
    "cdd": "CDD (transform then blend)",
    "hdd": "HDD (transform then blend)",
    "cdd_blended_temperature": "CDD (blend then transform)",
}

fig = ts_figure("Composite Daily Degree Days (equal-weighted, 4 stations)", "°F-days")
x = composite["date"].to_numpy()
renderers = [
    fig.line(x, composite[name].to_numpy(), color=color, line_width=1.8,
             legend_label=label, name=label,
             line_dash="dashed" if name == "cdd_blended_temperature" else "solid")
    for color, (name, label) in zip(SERIES, COMPOSITE_SERIES.items())
]
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("index", "$name"), ("day", "$x{%F}"),
              ("value", "$y{0.0}")], formatters={"$x": "datetime"})
)
fig.legend.location = "top_left"
show(style(fig))

Inputs: the four station indices and the equal-weighted temperature over the same stations and days. Operation: average the floored station values against flooring the averaged temperature. Validation: the basket sits at or above the blended alternative, with the gap coming entirely from shoulder days on which some stations are above 65°F while the blend is not.

convention_receipt = pl.DataFrame(
    {
        "convention": [
            str(COMPOSITE_CONVENTION),
            str(DegreeDayConvention.BLEND_THEN_TRANSFORM),
        ],
        "mean composite CDD": [
            composite["cdd"].mean(),
            composite["cdd_blended_temperature"].mean(),
        ],
        "max composite CDD": [
            composite["cdd"].max(),
            composite["cdd_blended_temperature"].max(),
        ],
        "days above zero": [
            (composite["cdd"] > 0).sum(),
            (composite["cdd_blended_temperature"] > 0).sum(),
        ],
    }
).with_columns(pl.col(pl.Float64).round(3), pl.lit(composite.height).alias("days"))
html_table(convention_receipt)

2 rows × 5 columns

convention	mean composite CDD	max composite CDD	days above zero	days
transform_then_blend	4.639	26.25	1,831	3,160
blend_then_transform	4.485	26.25	1,539	3,160

5.4 Temp-Derate MW and the Daily Market Data#

$W_d$ uses notebook 02’s edition/interval-overlap recipe, restricted to FORCED AMBIENT_DUE_TO_TEMP rows:

report = (
    pl.read_parquet(OS_DIR / "caiso_curtailed_non_operational_generator_report.parquet")
    .filter(
        (pl.col("outage_type") == "FORCED")
        & (pl.col("nature_of_work") == "AMBIENT_DUE_TO_TEMP")
    )
    .with_columns(
        # The market's own clock; the day back is the reporting convention.
        trade_date=CAISO_CLOCK.trade_day(pl.col("publish_time_utc")).dt.offset_by("-1d"),
    )
)
grid = (
    report.select("trade_date").unique().sort("trade_date")
    .with_columns(
        day_start=pl.col("trade_date").cast(pl.Datetime).dt.replace_time_zone("US/Pacific")
    )
    .with_columns(day_end=pl.col("day_start").dt.offset_by("1d"))
    .with_columns(
        hour=pl.datetime_ranges("day_start", "day_end", interval="1h", closed="left")
    )
    .explode("hour", empty_as_null=False)
    .with_columns(
        interval_start=pl.col("hour").dt.convert_time_zone("UTC").dt.replace_time_zone(None)
    )
    .with_columns(interval_end=pl.col("interval_start").dt.offset_by("1h"))
    .select("trade_date", "interval_start", "interval_end")
)
derate_daily = (
    grid.join_where(
        report.select(
            "trade_date", "curtailment_mw",
            pl.col("curtailment_start_time_utc").alias("out_start"),
            pl.col("curtailment_end_time_utc").alias("out_end"),
        ),
        pl.col("trade_date") == pl.col("trade_date_right"),
        pl.col("out_start") <= pl.col("interval_end"),
        pl.col("out_end").is_null() | (pl.col("out_end") >= pl.col("interval_start")),
    )
    .group_by("trade_date", "interval_start")
    .agg(pl.col("curtailment_mw").sum())
    .group_by(pl.col("trade_date").alias("date"))
    .agg(pl.col("curtailment_mw").mean().alias("temp_derate_mw"))
)

market_data = composite.join(rt_daily, on="date").join(derate_daily, on="date").sort("date")

6. Quantitative Analysis & Validation#

6.1 The Covariance Triple#

Heat raises A/C load and derates thermal capacity, so a high-CDD day tightens the system from both sides. That joint move is the covariance term, and it is what the quanto prices.

correlations = market_data.select(
    pl.corr("cdd", "rt_max").alias("corr(CDD, RT max)"),
    pl.corr("cdd", "temp_derate_mw").alias("corr(CDD, temp-derate MW)"),
    pl.corr("temp_derate_mw", "rt_max").alias("corr(temp-derate MW, RT max)"),
)
html_table(correlations)

1 rows × 3 columns

corr(CDD, RT max)	corr(CDD, temp-derate MW)	corr(temp-derate MW, RT max)
0.334	0.885	0.308

6.2 Tail Conditioning - The Quanto View#

The price leg’s fat right tail makes a full-sample correlation understate joint-tail dependence, so condition on the hottest days: top-decile CDD days against all others.

q90 = market_data["cdd"].quantile(0.9)
tails = (
    market_data.with_columns(
        pl.when(pl.col("cdd") >= q90)
        .then(pl.lit("Top-decile CDD days"))
        .otherwise(pl.lit("All other days"))
        .alias("bucket")
    )
    .group_by("bucket")
    .agg(
        pl.len().alias("days"),
        pl.col("cdd").mean().round(1).alias("mean CDD"),
        pl.col("rt_max").mean().round(1).alias("mean RT max ($/MWh)"),
        pl.col("rt_max").quantile(0.95).round(1).alias("p95 RT max ($/MWh)"),
        pl.col("temp_derate_mw").mean().round(0).alias("mean temp-derate MW"),
    )
    .sort("bucket", descending=True)
)
html_table(tails)

2 rows × 6 columns

bucket	days	mean CDD	mean RT max ($/MWh)	p95 RT max ($/MWh)	mean temp-derate MW
Top-decile CDD days	121	16.70	266.10	996.80	2,960.00
All other days	1,081	3.60	76.60	140.50	1,377.00

fig = figure(
    title="Daily Max SP15 RT LMP vs Composite CDD",
    x_axis_label="Composite CDD (°F-days)",
    y_axis_label="SP15 RT 15-min LMP daily max ($/MWh)",
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.scatter(market_data["cdd"].to_numpy(), market_data["rt_max"].to_numpy(),
            size=6, color=SERIES[0], alpha=0.45)
show(style(fig))

Top-decile CDD days carry a materially higher mean and p95 RT max, and the scatter shows the same tail-driven shape. The station panel has not graduated to the CAISO weather adapter, whose build starts in 2021 and carries no Burbank, so the settlement indices stay on openmeteo_temperature_daily and its 2018 start; the transform itself is pricing.indices.weather.DegreeDayIndex. The quanto payoff $\max(0, \bar C - K_w) \times \max(0, P - K_p)$ is priced in notebooks 13 and 14 as a heat rate call conditioned on a degree-day trigger; a load-weighted composite is notebook 04’s subject.

Categories

Tags