CAISO Outage Model Specification Search

Authors: Brian Nguyen · Claire Kim · Created: 2026-08-25 · Tags: #Experiment #Outages #Model-Selection

Notebook 11 fits one county-day NB2 specification. This notebook tests whether it is the right one: which engineered features earn a slot, which variance family fits the counts, whether the weights hold as history accrues, and how the extremes behave.

from pathlib import Path
import sys
import warnings

import numpy as np
import pandas as pd
import polars as pl
from bokeh.io import output_notebook, show
from bokeh.layouts import row
from bokeh.models import ColumnDataSource, HoverTool, Range1d, Span
from bokeh.plotting import figure
from bokeh.resources import INLINE
from bokeh.transform import dodge
from scipy.stats import nbinom, poisson
from statsmodels.tools.sm_exceptions import (
    ConvergenceWarning,
    HessianInversionWarning,
)

PROJECT_ROOT = Path.cwd()
if PROJECT_ROOT.name == "notebooks":
    PROJECT_ROOT = PROJECT_ROOT.parent
sys.path.insert(0, str(PROJECT_ROOT))

SRC_DIR = PROJECT_ROOT / "src"
if str(SRC_DIR) not in sys.path:
    sys.path.insert(0, str(SRC_DIR))

from src.caiso_outage_nb_pipeline import (
    COUNTY_DAILY_FEATURES,
    build_caiso_county_nb_panel,
    build_caiso_market_zone_nb_panel,
    fit_negative_binomial_model,
    likelihood_ratio_block_tests,
    model_performance_by_group,
    overdispersion_summary,
    poisson_negative_binomial_comparison,
    resource_mapping_coverage_summary,
    rollup_county_predictions_to_market_zone,
    rolling_backtest_negative_binomial_model,
)
from src.caiso_outage_specification_search import (
    add_candidate_interaction_columns,
    add_extreme_weather_columns,
    candidate_specifications,
    coefficient_stability_across_expanding_fits,
    compare_specifications,
    extreme_specifications,
    fit_nb1_model,
    l1_regularization_path,
    purged_walk_forward_family_comparison,
    variance_family_comparison,
    walk_forward_specification_metrics,
)
from chart_theme import SERIES, corr_heatmap, html_table, style

DATA_DIR = PROJECT_ROOT / "data"
PROCESSED_DIR = DATA_DIR / "processed"
PROCESSED_DIR.mkdir(parents=True, exist_ok=True)

pd.options.display.max_columns = 80
pd.options.display.float_format = "{:.4f}".format

# The county fits sit on the alpha = 0 dispersion boundary, so statsmodels
# reports a failed convergence check and an uninvertible Hessian while still
# returning the Poisson-limit coefficients. The tables below carry the fitted
# alpha and the held-out scores, which is where a failed fit shows.
warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", category=HessianInversionWarning)

output_notebook(INLINE, hide_banner=True)

Setup: The Panel and the Incumbent Fit

The candidates are scored against notebook 11’s fitted model on the same rows and the same chronological split.

county_panel, mapped_events, resource_crosswalk, county_zone_weights = build_caiso_county_nb_panel(DATA_DIR)
county_fit = fit_negative_binomial_model(
    county_panel,
    feature_cols=COUNTY_DAILY_FEATURES,
    high_count_quantile=0.90,
    maxiter=300,
)
county_coef = county_fit.coefficients.copy()

An NB2 fit keeps every feature it is given. Whether each feature deserves its slot, and whether NB2 itself is the right distribution, are separate exercises, both run here through src/caiso_outage_specification_search.py and judged on held-out performance first:

• (i) named candidate feature sets, compared on AIC/BIC and held-out surprise;

• (ii) an L1 penalty path that drives weak coefficients exactly to zero, showing the order in which features die;

• (iii) a variance-family comparison - Poisson, NB1, NB2, and a zero-inflated NB2 - on the same features;

• (iv) coefficient stability across expanding training sets.

1. Named Candidate Specifications

Notebook 11 section 4’s block tests picked the candidates, and each is fitted on the same rows and the same chronological split, so the scores are directly comparable. aic and bic are the training fit penalized for parameter count (lower is better; BIC charges more per parameter), and the rule throughout is to prefer the simplest specification within noise of the best held-out NLL.

county_panel_extended = add_candidate_interaction_columns(county_panel)
SPECIFICATIONS = candidate_specifications()

specification_table = pd.DataFrame(
    {
        "specification": list(SPECIFICATIONS),
        "features": [len(c.features) for c in SPECIFICATIONS.values()],
        "what it tests": [c.tests for c in SPECIFICATIONS.values()],
    }
)
html_table(pl.from_pandas(specification_table), kind="reference")

reference · 7 rows × 3 columns

specification	features	what it tests
full (incumbent)	15	the section 3 model, unchanged - the reference point
no load	14	drops the one block the section 4 tests found insignificant
no weather	11	whether the weather block earns its slots out of sample
persistence + size	5	how far the two dominant blocks go alone
persistence only	3	the floor - yesterday's outage state and nothing else
full + CDD squared	16	whether outage risk accelerates in extreme heat
full + CDD x capacity	16	whether heat matters more where more MW is installed

specification_comparison = compare_specifications(county_panel_extended, SPECIFICATIONS)

comparison_sorted = specification_comparison.sort_values("test_nll", ascending=False)
fig = figure(
    title="Held-Out Surprise by Candidate Specification",
    y_range=comparison_sorted["specification"].tolist(),
    x_range=Range1d(
        comparison_sorted["test_nll"].min() - 0.005,
        comparison_sorted["test_nll"].max() + 0.005,
    ),
    x_axis_label="Held-Out NLL (lower is better)",
    height=340,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.hbar(
    y=comparison_sorted["specification"].tolist(),
    right=comparison_sorted["test_nll"].tolist(),
    height=0.6,
    color=[
        SERIES[1] if name == "full (incumbent)" else SERIES[0]
        for name in comparison_sorted["specification"]
    ],
)
fig.add_tools(HoverTool(tooltips=[("specification", "@y"), ("test NLL", "@right{0.0000}")]))
show(style(fig))

html_table(pl.from_pandas(specification_comparison.round(4)))

7 rows × 10 columns

specification	features	train_llf	aic	bic	alpha	test_nll	test_mae	test_auc_outage	test_auc_high_outage
full (incumbent)	15	−79,124.85	158,283.70	158,441.23	0.013	1.015	0.745	0.885	0.979
no load	14	−79,134.22	158,300.43	158,448.69	0.013	1.015	0.744	0.885	0.979
no weather	11	−79,148.55	158,323.10	158,443.56	0.012	1.016	0.745	0.885	0.979
persistence + size	5	−80,236.70	160,487.39	160,552.26	0.027	1.022	0.770	0.893	0.978
persistence only	3	−82,439.68	164,889.36	164,935.69	0.056	1.062	0.798	0.846	0.970
full + CDD squared	16	−79,122.82	158,281.63	158,448.43	0.013	1.015	0.745	0.885	0.979
full + CDD x capacity	16	−79,049.37	158,134.74	158,301.53	0.012	1.015	0.743	0.884	0.979

The incumbent is the orange bar. Held-out surprise falls from 1.062 (persistence only) to 1.022 (persistence + size) and then flattens: every full-size variant lands at 1.015 or 1.016, CDD x capacity included. The interaction separates itself on the training criteria instead, about 149 AIC points and 140 BIC points below the incumbent.

engineered = specification_comparison[
    specification_comparison["specification"].str.startswith("full + ")
]
best_engineered = engineered.sort_values("test_nll").iloc[0]["specification"]
finalists = ["full (incumbent)", "no load", "no weather", best_engineered]

walk_forward_finalists = walk_forward_specification_metrics(
    county_panel_extended, SPECIFICATIONS, finalists
)
html_table(pl.from_pandas(walk_forward_finalists.round(4)))

4 rows × 8 columns

specification	features	nll	mae	auc_outage	auc_high_outage	brier_outage	brier_high_outage
full (incumbent)	15	1.029	0.753	0.891	0.979	0.129	0.027
no load	14	1.029	0.752	0.891	0.979	0.129	0.027
no weather	11	1.029	0.752	0.891	0.979	0.129	0.027
full + CDD x capacity	16	1.028	0.751	0.891	0.979	0.128	0.027

Under notebook 11’s purged walk-forward, CDD x capacity takes AIC, BIC, the walk-forward surprise (1.028 vs 1.029) and count error (0.751 vs 0.753); held-out surprise on the single split ties at 1.015 and the AUCs sit within 0.001. At that margin the choice is not forced: notebook 11 keeps the incumbent as its model of record, and the interaction is the recommended addition when the rollup feeds pricing. Dropping the weather block entirely costs a thousandth of a nat out of sample.

2. The L1 Path: Which Coefficients Die First

The penalty is charged on the standardized weights, never on the intercept or \(\alpha\), so sweeping it from nothing to punishing traces the model’s priority order.

l1_bundle = l1_regularization_path(
    county_panel, list(COUNTY_DAILY_FEATURES), [0, 5, 15, 50, 150, 500, 1500, 5000]
)

penalties = sorted(l1_bundle.path_metrics["penalty"])
spine_penalty = penalties[-2]
spine_features = l1_bundle.coefficient_paths.query(
    "penalty == @spine_penalty and abs(coefficient) > 1e-6"
)["feature"].tolist()
spine_colors = dict(zip(spine_features, SERIES))

positive_penalty_paths = l1_bundle.coefficient_paths[
    l1_bundle.coefficient_paths["penalty"] > 0
]
fig = figure(
    title="L1 Path: Coefficients as the Penalty Grows",
    x_axis_type="log",
    x_axis_label="Penalty per Unit of Coefficient (log scale)",
    y_axis_label="Standardized Coefficient",
    height=400,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
renderers = []
for feature, group in positive_penalty_paths.groupby("feature"):
    in_spine = feature in spine_colors
    line_kwargs = {
        "color": spine_colors.get(feature, "#c3c2b7"),
        "line_width": 2.4 if in_spine else 1.3,
        "name": feature,
    }
    if in_spine:
        line_kwargs["legend_label"] = feature
    renderers.append(
        fig.line(group["penalty"].to_numpy(), group["coefficient"].to_numpy(), **line_kwargs)
    )
fig.add_layout(Span(location=0, dimension="width", line_color="black", line_width=1))
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("feature", "$name"), ("penalty", "@x{0}"),
              ("coefficient", "@y{0.000}")])
)
fig.legend.location = "top_right"
show(style(fig))

html_table(pl.from_pandas(l1_bundle.path_metrics.round(4)))

8 rows × 5 columns

penalty	surviving features	alpha	test_nll	converged
0.00	15	0.013	1.015	True
5.00	15	0.013	1.015	True
15.00	15	0.013	1.015	True
50.00	15	0.013	1.015	True
150.00	13	0.012	1.015	True
500.00	12	0.012	1.016	True
1,500.00	9	0.011	1.017	True
5,000.00	7	0.018	1.028	True

Colored lines are the last survivors; gray is everything else. The spine is persistence plus county size, and the degree-day and seasonal weights are the first pulled toward the axis. surviving features counts coefficients above 1e-6 in size, so it holds at fifteen until penalty 150 - where the shrunken weights finally cross - then falls to thirteen, twelve at 500, nine at 1,500 and seven at 5,000. Held-out surprise is flat at 1.015 through penalty 150 and reaches 1.028 only at 5,000, where the spine itself is hit. No pruning is needed for performance, but a lean operational version costs almost nothing: nine features at penalty 1,500 give up 0.002 nats.

3. Is the Variance Family Itself Right?

Same features, same split, four distributions scored on their own held-out surprise: Poisson, NB1 (extra variance linear in the mean, \(\mathrm{Var} = (1+\alpha)\,\mu\)), NB2 (quadratic), and a zero-inflated NB2.

family_comparison = variance_family_comparison(county_panel, list(COUNTY_DAILY_FEATURES))
html_table(pl.from_pandas(family_comparison.round(4)))

4 rows × 8 columns

model	parameters	train_llf	aic	bic	alpha	test_nll	converged
Poisson	16	−79,152.40	158,336.80	158,485.06	—	1.019	True
NB1 (variance linear in the mean)	17	−78,298.78	156,631.56	156,789.08	0.242	1.001	True
NB2 (variance quadratic in the mean)	17	−79,124.85	158,283.70	158,441.23	0.013	1.015	True
zero-inflated NB2 (constant inflation)	18	—	—	—	—	—	False

• NB1 wins. Training likelihood is about 830 log points above NB2’s, AIC and BIC both prefer it, and held-out surprise drops from 1.015 to 1.001. The extra scatter is a constant premium of roughly 24% over Poisson variance (\(\alpha \approx 0.24\)), not dispersion growing with the mean.

• Zero inflation offers nothing usable. The fit fails to converge on this panel, so it returns no scores at all; the 60% zero share is already carried by the features and the count model’s own zero mass.

• The pipeline’s probability formulas are NB2-specific, so src/caiso_outage_specification_search.py implements the NB1 pair (\(\Pr[Y = 0] = (1 + \alpha)^{-\mu/\alpha}\) and the matching tail):

county_fit_nb1 = fit_nb1_model(county_panel, list(COUNTY_DAILY_FEATURES))

scorecard_keys = ["alpha", "test_nll", "test_mae", "test_rmse", "test_auc_outage",
                  "test_auc_high_outage", "test_brier_outage", "test_brier_high_outage"]
family_scorecard = pd.DataFrame(
    [
        {"model": "NB2 (incumbent)",
         **{key: county_fit.metrics[key] for key in scorecard_keys}},
        {"model": "NB1",
         **{key: county_fit_nb1.metrics[key] for key in scorecard_keys}},
    ]
)
html_table(pl.from_pandas(family_scorecard.round(4)))

2 rows × 9 columns

model	alpha	test_nll	test_mae	test_rmse	test_auc_outage	test_auc_high_outage	test_brier_outage	test_brier_high_outage
NB2 (incumbent)	0.013	1.015	0.745	1.426	0.885	0.979	0.131	0.025
NB1	0.242	1.001	0.744	1.426	0.885	0.979	0.130	0.025

Both families through the same purged expanding folds on notebook 11’s section 5 settings (the embargo differs slightly, so the NB2 row is close to but not exactly that table). Every row and fold is shared, so any difference is the variance form:

family_walk_forward = purged_walk_forward_family_comparison(
    county_panel, list(COUNTY_DAILY_FEATURES)
)
html_table(pl.from_pandas(family_walk_forward.round(4)))

2 rows × 9 columns

model	rows	mae	rmse	nll	auc_outage	auc_high_outage	brier_outage	brier_high_outage
NB2 (incumbent)	39,096	0.753	1.406	1.029	0.891	0.979	0.129	0.027
NB1	39,096	0.752	1.407	1.017	0.891	0.979	0.128	0.027

NB1’s surprise is lower on both the single split (1.001 vs 1.015) and the identical folds (1.017 vs 1.029), with Brier scores equal or a touch better. Point error is a wash (walk-forward MAE 0.752 vs 0.753) and AUCs sit within 0.001.

fit_nb1_model returns the same fitted-frame contract as the NB2 bundle (mu_hat, p_outage_hat, p_high_outage_hat, sample), so notebook 11’s section 6 rollup consumes it unchanged.

3.1 The Two Changes Together

Section 1 chose a feature and section 3 chose a variance family, each holding the other fixed. Whether they compose is a separate question, so all four combinations are fitted on the same rows and the same split:

interaction_features = list(COUNTY_DAILY_FEATURES) + ["cdd65_x_log_capacity"]
family_feature_grid = {
    ("NB2", "full"): (fit_negative_binomial_model, list(COUNTY_DAILY_FEATURES)),
    ("NB2", "full + CDD x capacity"): (fit_negative_binomial_model, interaction_features),
    ("NB1", "full"): (fit_nb1_model, list(COUNTY_DAILY_FEATURES)),
    ("NB1", "full + CDD x capacity"): (fit_nb1_model, interaction_features),
}
pair_rows = []
for (family, features_label), (fitter, features) in family_feature_grid.items():
    bundle = fitter(county_panel_extended, features)
    pair_rows.append(
        {"variance family": family, "features": features_label,
         **{key: bundle.metrics[key] for key in scorecard_keys}}
    )
family_feature_scorecard = pd.DataFrame(pair_rows)
# The whole question here lives in the fourth decimal of test_nll, which the table
# renders to three, so the gain is shown in thousandths of a nat as well.
baseline_nll = family_feature_scorecard["test_nll"].iloc[0]
family_feature_scorecard.insert(
    3,
    "surprise gain vs NB2 full (millinats)",
    (baseline_nll - family_feature_scorecard["test_nll"]) * 1000,
)
html_table(pl.from_pandas(family_feature_scorecard.round(4)))

4 rows × 11 columns

variance family	features	alpha	surprise gain vs NB2 full (millinats)	test_nll	test_mae	test_rmse	test_auc_outage	test_auc_high_outage	test_brier_outage	test_brier_high_outage
NB2	full	0.013	0.00	1.015	0.745	1.426	0.885	0.979	0.131	0.025
NB2	full + CDD x capacity	0.012	0.44	1.015	0.743	1.428	0.884	0.979	0.131	0.025
NB1	full	0.242	13.95	1.001	0.744	1.426	0.885	0.979	0.130	0.025
NB1	full + CDD x capacity	0.240	14.42	1.001	0.742	1.429	0.884	0.979	0.130	0.025

The two changes do not interact: the gain column is additive. Moving to NB1 is worth about 14 thousandths of a nat with or without the interaction, and adding the interaction is worth about 0.4 with or without NB1. The pair is the best of the four cells - 14.4 against 13.9 for NB1 alone - and it takes count error too (MAE 0.742 against 0.745 for the incumbent), which is what the recommendation rests on. But the interaction is worth about a thirtieth of the family change and gives a little back on RMSE (1.429 against 1.426) and any-outage AUC (0.884 against 0.885). NB1 is the change that matters; the interaction is a refinement on top of it, not a second lever.

4. Do the Weights Hold Up Over Time?

Refitting on the first 60, 70, 80 and 90 percent of dates - expanding and strictly out-of-time - shows whether the weights are stable properties of the process or artifacts of one sample.

stability = coefficient_stability_across_expanding_fits(
    county_panel, list(COUNTY_DAILY_FEATURES)
)

lead_features = (
    county_coef[county_coef["feature"].ne("const")]
    .sort_values("abs_coefficient", ascending=False)["feature"]
    .head(4)
    .tolist()
)
watch_features = list(dict.fromkeys([*lead_features, "cdd65", "hdd65"]))

fig = figure(
    title="Coefficient Stability Across Expanding Training Sets",
    x_axis_type="datetime",
    x_axis_label="Training Set End Date",
    y_axis_label="Standardized Coefficient",
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
renderers = []
for color, feature in zip(SERIES, watch_features):
    feature_history = stability[stability["feature"] == feature].sort_values("train_end")
    x = pd.to_datetime(feature_history["train_end"])
    fig.segment(x, feature_history["conf_low"].to_numpy(),
                x, feature_history["conf_high"].to_numpy(),
                color=color, line_width=2, alpha=0.35)
    fig.line(x, feature_history["coefficient"].to_numpy(), color=color, line_width=2,
             legend_label=feature)
    renderers.append(
        fig.scatter(x, feature_history["coefficient"].to_numpy(), color=color, size=7,
                    name=feature)
    )
fig.add_layout(Span(location=0, dimension="width", line_color="black", line_width=1))
fig.add_tools(
    HoverTool(renderers=renderers, tooltips=[("feature", "$name"),
              ("train end", "@x{%F}"), ("coefficient", "@y{0.000}")],
              formatters={"@x": "datetime"})
)
fig.legend.location = "top_right"
show(style(fig))

Whiskers are each fit’s 95% range. The structural weights are flat: yesterday’s outage count holds between 0.78 and 0.80 at every horizon, county capacity climbs from 0.32 to 0.36, and nothing flips sign. The degree-day weights are indistinguishable from zero throughout (cdd65 drifts between 0.004 and 0.009, hdd65 between 0.007 and 0.009, whiskers spanning zero at all four cutoffs), so weather’s significance in notebook 11 is a property of the block - jointly with temperature_f_std and the multi-proxy flag - not of any single degree-day coefficient.

5. The Extremes: Hottest and Coldest Days

Two questions: does the model miss the most extreme county-days, and do tail-only features earn a slot? First the audit: mean actual against mean predicted by degree-day bin.

extreme_bin_edges = [-0.01, 0.0, 5.0, 10.0, 15.0, 100.0]
extreme_bin_labels = ["0", "(0, 5]", "(5, 10]", "(10, 15]", "15+"]
calibration_pieces = []
for weather_col, weather_label in (("cdd65", "CDD65"), ("hdd65", "HDD65")):
    binned = county_fit.fitted_df.copy()
    binned["bin"] = pd.cut(
        binned[weather_col], bins=extreme_bin_edges, labels=extreme_bin_labels
    )
    grouped = (
        binned.groupby(["bin", "sample"], observed=True)
        .agg(
            days=("Y", "size"),
            actual=("Y", "mean"),
            predicted=("mu_hat", "mean"),
        )
        .reset_index()
    )
    grouped.insert(0, "degree days", weather_label)
    calibration_pieces.append(grouped)
extreme_calibration = pd.concat(calibration_pieces, ignore_index=True).round(3)

held_out_heat = extreme_calibration[
    extreme_calibration["degree days"].eq("CDD65")
    & extreme_calibration["sample"].eq("test")
]
heat_source = ColumnDataSource(
    {
        "bin": held_out_heat["bin"].astype(str).tolist(),
        "actual": held_out_heat["actual"].tolist(),
        "predicted": held_out_heat["predicted"].tolist(),
        "days": held_out_heat["days"].tolist(),
    }
)
fig = figure(
    title="Held-Out Calibration on the Hottest Days",
    x_range=extreme_bin_labels,
    x_axis_label="County-Day CDD65 Bin",
    y_axis_label="Mean Outage Starts",
    height=340,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.vbar(x=dodge("bin", -0.16, range=fig.x_range), top="actual", width=0.28,
         color=SERIES[0], legend_label="actual", source=heat_source)
fig.vbar(x=dodge("bin", 0.16, range=fig.x_range), top="predicted", width=0.28,
         color=SERIES[1], legend_label="predicted", source=heat_source)
fig.add_tools(
    HoverTool(tooltips=[("bin", "@bin"), ("actual", "@actual{0.00}"),
              ("predicted", "@predicted{0.00}"), ("held-out days", "@days{0,0}")])
)
fig.legend.location = "top_left"
show(style(fig))

html_table(pl.from_pandas(extreme_calibration))

20 rows × 6 columns

degree days	bin	sample	days	actual	predicted
CDD65	0	test	9,964	1.144	1.174
CDD65	0	train	42,267	1.387	1.392
CDD65	(0, 5]	test	3,015	1.604	1.679
CDD65	(0, 5]	train	10,005	1.828	1.829
CDD65	(5, 10]	test	2,758	1.314	1.352
CDD65	(5, 10]	train	9,134	1.818	1.849
CDD65	(10, 15]	test	1,909	1.164	1.195
CDD65	(10, 15]	train	7,088	1.435	1.425
CDD65	15+	test	1,902	1.478	1.468
CDD65	15+	train	9,644	1.778	1.799
HDD65	0	test	9,298	1.372	1.407
HDD65	0	train	34,977	1.690	1.700
HDD65	(0, 5]	test	3,503	1.484	1.536
HDD65	(0, 5]	train	8,718	1.927	1.934
HDD65	(5, 10]	test	2,436	1.214	1.230
HDD65	(5, 10]	train	9,290	1.730	1.723
HDD65	(10, 15]	test	1,742	1.051	1.060
HDD65	(10, 15]	train	11,422	1.316	1.334
HDD65	15+	test	2,569	0.838	0.876
HDD65	15+	train	13,731	1.006	1.011

Tail candidates, all built from constants fixed up front so nothing leaks from later dates: hinges at 8 and 15 cooling degrees (the deep one covers roughly the hottest tenth of county-days), at 10 and 18 heating degrees, both moderate hinges together, and a backward three-day county mean of CDD65. Single split first, then the purged walk-forward:

county_panel_extremes = add_extreme_weather_columns(county_panel)
EXTREME_SPECIFICATIONS = extreme_specifications()

extreme_comparison = compare_specifications(county_panel_extremes, EXTREME_SPECIFICATIONS)
html_table(
    pl.from_pandas(
        extreme_comparison[
            ["specification", "features", "aic", "bic", "test_nll", "test_mae",
             "test_auc_high_outage"]
        ].round(4)
    )
)

7 rows × 7 columns

specification	features	aic	bic	test_nll	test_mae	test_auc_high_outage
full (incumbent)	15	158,283.70	158,441.23	1.015	0.745	0.979
full + CDD hinge (8)	16	158,285.13	158,451.92	1.015	0.745	0.979
full + deep CDD hinge (15)	16	158,282.13	158,448.92	1.015	0.745	0.979
full + HDD hinge (10)	16	158,279.99	158,446.78	1.016	0.746	0.979
full + deep HDD hinge (18)	16	158,261.37	158,428.16	1.016	0.746	0.979
full + both hinges (8, 10)	17	158,281.97	158,458.03	1.016	0.746	0.979
full + 3-day CDD	16	158,273.64	158,440.43	1.015	0.745	0.979

extreme_walk_forward = walk_forward_specification_metrics(
    county_panel_extremes, EXTREME_SPECIFICATIONS, list(EXTREME_SPECIFICATIONS)
)
html_table(pl.from_pandas(extreme_walk_forward.round(4)))

7 rows × 8 columns

specification	features	nll	mae	auc_outage	auc_high_outage	brier_outage	brier_high_outage
full (incumbent)	15	1.029	0.753	0.891	0.979	0.129	0.027
full + CDD hinge (8)	16	1.029	0.753	0.891	0.979	0.129	0.027
full + deep CDD hinge (15)	16	1.029	0.753	0.891	0.979	0.129	0.027
full + HDD hinge (10)	16	1.029	0.753	0.891	0.979	0.129	0.027
full + deep HDD hinge (18)	16	1.029	0.754	0.891	0.979	0.129	0.027
full + both hinges (8, 10)	17	1.029	0.754	0.891	0.979	0.129	0.027
full + 3-day CDD	16	1.029	0.753	0.891	0.979	0.129	0.027

• The model tracks the extremes. Five summers put 1,902 held-out county-days in the hottest bin (more than 15 cooling degrees) and 2,569 in the deepest cold bin. Prediction runs mildly high in every bin except the hottest, where it sits 0.7% under (1.468 against 1.478 actual).

• The tail candidates stay marginal. The training criteria mildly prefer the deep heating-degree hinge (22 AIC points, 13 BIC) and the sustained-heat term (10 and 1), but none of the seven separates out of sample: walk-forward surprise is 1.029 for every one of them and MAE spans 0.753 to 0.754. One parameter spent on the CDD x capacity interaction buys about seven times the AIC.

• The sample was not the limit. The panel starts 2021-06-17; the hottest bin holds about 11,500 county-days against about 1,700 on a single-summer sample, and the tail answer barely moved.

Findings

• CDD65 interacted with installed capacity is the one engineered feature that improves both in-sample criteria and walk-forward performance.

• NB1’s constant dispersion premium fits the counts better than NB2’s quadratic over five summers.

• The persistence-plus-size spine is stable across expanding fits, and the extremes track closely.

• Recommended production pair when the rollup feeds pricing: NB1 with the CDD x capacity interaction, fitted together in section 3.1 and the best of the four cells on held-out surprise and count error - though nearly all of that gain is NB1’s. Notebook 11 keeps NB2 for continuity.