CAISO County Outage Model

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

County-day is the intermediate spatial layer between generator-level outage records and NP15/SP15/ZP26 pricing zones. The panel keeps market-zone labels for spark-spread integration, but lets weather enter closer to the plant location.

from pathlib import Path
import sys

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

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

output_notebook(INLINE, hide_banner=True)

1. County Panel Construction

The mapping path is CAISO outage resource → matched EIA plant → county and CAISO market zone. County weather comes from the plant inventory as capacity-weighted Open-Meteo proxy locations.

county_panel, mapped_events, resource_crosswalk, county_zone_weights = build_caiso_county_nb_panel(DATA_DIR)
market_panel, _, _ = build_caiso_market_zone_nb_panel(DATA_DIR)

coverage = resource_mapping_coverage_summary(
    mapped_events.drop(columns=[col for col in ["region", "county"] if col in mapped_events.columns], errors="ignore"),
    mapped_events,
    resource_crosswalk,
    start_date=county_panel["date"].min(),
    end_date=county_panel["date"].max(),
)

panel_shape = pd.DataFrame([
    {
        "start_date": county_panel["date"].min(),
        "end_date": county_panel["date"].max(),
        "county_days": len(county_panel),
        "county_zone_pairs": county_panel[["region", "county"]].drop_duplicates().shape[0],
        "market_zone_days": len(market_panel),
        "matched_forced_event_coverage": coverage.loc[coverage["outage_type"].eq("FORCED"), "event_coverage"].iloc[0],
    }
])

html_table(pl.from_pandas(panel_shape))

1 rows × 6 columns

start_date	end_date	county_days	county_zone_pairs	market_zone_days	matched_forced_event_coverage
2021-06-17	2026-06-30	99,360	54	5,520	0.606

The county panel keeps the market-zone model’s matched-resource sample; the gain is spatial resolution, not broader coverage.

preview_cols = [
    "date",
    "region",
    "county",
    "Y",
    "forced_start_mw_sum",
    "planned_start_count",
    "temperature_f_mean",
    "cdd65",
    "hdd65",
    "county_nameplate_capacity_mw",
    "county_weather_basis_risk_flag",
]
html_table(pl.from_pandas(county_panel.loc[:, preview_cols].head(10)))

10 rows × 11 columns

date	region	county	Y	forced_start_mw_sum	planned_start_count	temperature_f_mean	cdd65	hdd65	county_nameplate_capacity_mw	county_weather_basis_risk_flag
2021-06-17	NP15	Alameda	3	10.15	0	80.66	15.66	0.000	1,347.50	0
2021-06-18	NP15	Alameda	2	6.00	0	80.42	15.42	0.000	1,347.50	0
2021-06-19	NP15	Alameda	2	8.00	0	73.00	8.00	0.000	1,347.50	0
2021-06-20	NP15	Alameda	2	8.00	0	70.77	5.77	0.000	1,347.50	0
2021-06-21	NP15	Alameda	3	100.04	0	70.00	5.00	0.000	1,347.50	0
2021-06-22	NP15	Alameda	6	92.00	0	70.09	5.09	0.000	1,347.50	0
2021-06-23	NP15	Alameda	3	100.04	0	66.74	1.74	0.000	1,347.50	0
2021-06-24	NP15	Alameda	2	8.00	0	63.25	0.00	1.754	1,347.50	0
2021-06-25	NP15	Alameda	2	7.00	1	64.62	0.00	0.379	1,347.50	0
2021-06-26	NP15	Alameda	2	8.00	0	65.85	0.85	0.000	1,347.50	0

The weather fields are daily county exposures, while region remains the pricing anchor.

2. Target Sparsity

County-day counts are sparser than market-zone-day counts but not as thin as county-hour, which makes them a workable first robustness layer.

county_target_summary = (
    county_panel.groupby("region")
    .agg(
        counties=("county", "nunique"),
        rows=("Y", "size"),
        mean_count=("Y", "mean"),
        variance=("Y", "var"),
        zero_share=("Y", lambda s: float((s == 0).mean())),
        outage_share=("Y", lambda s: float((s > 0).mean())),
        p90=("Y", lambda s: float(s.quantile(0.90))),
        max_count=("Y", "max"),
    )
    .reset_index()
)
county_target_summary["variance_to_mean"] = county_target_summary["variance"] / county_target_summary["mean_count"]
html_table(pl.from_pandas(county_target_summary))

3 rows × 10 columns

region	counties	rows	mean_count	variance	zero_share	outage_share	p90	max_count	variance_to_mean
NP15	30	55,200	0.821	2.42	0.668	0.332	3.00	14	2.95
SP15	8	14,720	3.719	15.55	0.277	0.723	10.00	24	4.18
ZP26	16	29,440	1.616	17.54	0.648	0.352	4.00	46	10.85

target_regions = county_target_summary["region"].tolist()

fig = figure(
    title="County-Day Zero Share",
    x_range=target_regions,
    y_range=Range1d(0, 1),
    y_axis_label="Share of County-Days With Y = 0",
    height=340,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.vbar(
    x=target_regions,
    top=county_target_summary["zero_share"].tolist(),
    width=0.6,
    color=SERIES[0],
)
fig.add_tools(HoverTool(tooltips=[("region", "@x"), ("zero share", "@top{0.000}")]))
show(style(fig))

fig = figure(
    title="County-Day Overdispersion",
    x_range=target_regions,
    y_axis_label="Variance / Mean",
    height=340,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.vbar(
    x=target_regions,
    top=county_target_summary["variance_to_mean"].tolist(),
    width=0.6,
    color=SERIES[1],
)
fig.add_layout(
    Span(location=1, dimension="width", line_color="black", line_width=1, line_dash="dashed")
)
fig.add_tools(HoverTool(tooltips=[("region", "@x"), ("variance / mean", "@top{0.00}")]))
show(style(fig))

SP15 has the densest county-day outage process, while NP15 and ZP26 are more zero-heavy. Variance remains above the mean, so the negative binomial variance structure is still warranted at this finer spatial layer.

3. County NB2 Specification

The target is \(Y_{c,d}\), forced outage starts in county \(c\) on day \(d\) (Y), negative binomial with mean \(\mu_{c,d}\) (mu_hat), dispersion \(\alpha\) (alpha) and variance var_hat:

\[\mu_{c,d} = \exp\!\big(\beta_{0} + x_{c,d}^{\top}\beta\big), \qquad \mathrm{Var}[Y] = \mu + \alpha\,\mu^{2}\]

\[\Pr[Y = y] = \frac{\Gamma(y + 1/\alpha)}{\Gamma(1/\alpha)\,y!}\left(\frac{1}{1 + \alpha\mu}\right)^{1/\alpha}\left(\frac{\alpha\mu}{1 + \alpha\mu}\right)^{y}\]

Section 2’s variance-to-mean ratios of 3 to 11 rule out the Poisson restriction \(\alpha = 0\); section 4 tests the two head-to-head. Features are standardized, lag terms enter as \(\log(1 + x)\), size and load as logs, and \(\text{CDD65} = \max(0,\ T - 65)\) / \(\text{HDD65} = \max(0,\ 65 - T)\) are computed on the county’s capacity-weighted proxy temperature.

The scorecards grade \(\Pr[Y \geq 1] = 1 - (1 + \alpha\,\mu_{c,d})^{-1/\alpha}\) (p_outage_hat) and \(\Pr[Y \geq q]\) (p_high_outage_hat), \(q\) the smallest whole number at or above the 0.90 quantile (test_high_count_threshold).

def count_probabilities(mu, alpha, counts):
    """P[Y = k] under NB2(mu, alpha); alpha = 0 is exactly Poisson."""
    if alpha == 0:
        return poisson.pmf(counts, mu)
    return nbinom.pmf(counts, 1.0 / alpha, 1.0 / (1.0 + alpha * mu))


# 0.013 mirrors the fitted alpha in the section 3 scorecard below
ALPHA_SERIES = {
    "Poisson (α = 0)": 0.0,
    "NB2, fitted α = 0.013": 0.013,
    "NB2, α = 1 (exaggerated)": 1.0,
}
shape_figures = []
for mu, count_max, day_label in [(1, 8, "a Quiet County-Day"), (10, 30, "a Busy County-Day")]:
    counts = np.arange(count_max + 1)
    fig = figure(
        title=f"Distribution Shapes on {day_label} (μ = {mu})",
        x_axis_label="Outages Started in the Day",
        y_axis_label="Probability",
        height=360,
        sizing_mode="stretch_width",
        toolbar_location="above",
    )
    renderers = []
    for color, (label, alpha) in zip(SERIES, ALPHA_SERIES.items()):
        pmf = count_probabilities(mu, alpha, counts)
        fig.line(counts, pmf, color=color, line_width=1.2, line_dash="dotted")
        renderers.append(
            fig.scatter(counts, pmf, color=color, size=7, legend_label=label, name=label)
        )
    fig.add_tools(
        HoverTool(renderers=renderers, tooltips=[("model", "$name"), ("count", "@x"),
                  ("probability", "@y{0.0000}")])
    )
    fig.legend.location = "top_right"
    shape_figures.append(style(fig))
show(row(*shape_figures, sizing_mode="stretch_width"))

The fitted \(\alpha = 0.013\) is visually Poisson - on a busy day the variance is 11.3 against Poisson’s 10 - so most of section 2’s raw overdispersion was differences in \(\mu\) across counties and seasons that the features absorb.

The Features

The feature vector \(x_{c,d}\): one row per input, one fitted weight each in the coefficient chart.

feature_table = pd.DataFrame({
    "feature": COUNTY_DAILY_FEATURES,
    "role": [
        "cooling-degree weather exposure",
        "heating-degree weather exposure",
        "system demand level",
        "county outage persistence",
        "planned-maintenance persistence",
        "lagged forced MW persistence",
        "within-county weather-proxy dispersion",
        "county generation footprint",
        "county plant-count footprint",
        "multi-proxy county flag",
        "calendar effect",
        "seasonal cycle",
        "seasonal cycle",
        "NP15 market-zone control",
        "SP15 market-zone control",
    ],
    "why it should matter": [
        "hot days stress cooling equipment and force units below their rated output",
        "cold snaps freeze instrumentation and squeeze gas supply",
        "high demand puts more units in service, so more can fail",
        "trouble clusters in time - the stressors behind yesterday's outage starts usually persist into today",
        "active maintenance campaigns keep crews on site doing equipment work",
        "the MW size of yesterday's outages adds information beyond the count",
        "uneven weather inside a county means its average understates local stress",
        "more installed MW is more equipment that can break",
        "more distinct plants is more independent chances of a failure",
        "plants spanning multiple weather proxies or climate zones blur the county's single weather series",
        "staffing and maintenance scheduling differ on weekends",
        "smooth annual seasonality, first coordinate",
        "smooth annual seasonality, second coordinate",
        "level difference of NP15 counties the other features miss (ZP26 is the baseline)",
        "level difference of SP15 counties the other features miss",
    ],
})
html_table(pl.from_pandas(feature_table), kind="reference")

reference · 15 rows × 3 columns

feature	role	why it should matter
cdd65	cooling-degree weather exposure	hot days stress cooling equipment and force units below their rated output
hdd65	heating-degree weather exposure	cold snaps freeze instrumentation and squeeze gas supply
log_system_load_mean	system demand level	high demand puts more units in service, so more can fail
log1p_Y_lag1	county outage persistence	trouble clusters in time - the stressors behind yesterday's outage starts usually persist into today
log1p_planned_start_count_lag1	planned-maintenance persistence	active maintenance campaigns keep crews on site doing equipment work
log1p_forced_start_mw_sum_lag1	lagged forced MW persistence	the MW size of yesterday's outages adds information beyond the count
temperature_f_std	within-county weather-proxy dispersion	uneven weather inside a county means its average understates local stress
log_county_nameplate_capacity_mw	county generation footprint	more installed MW is more equipment that can break
log_county_plant_count	county plant-count footprint	more distinct plants is more independent chances of a failure
county_weather_basis_risk_flag	multi-proxy county flag	plants spanning multiple weather proxies or climate zones blur the county's single weather series
weekend	calendar effect	staffing and maintenance scheduling differ on weekends
month_sin	seasonal cycle	smooth annual seasonality, first coordinate
month_cos	seasonal cycle	smooth annual seasonality, second coordinate
region_NP15	NP15 market-zone control	level difference of NP15 counties the other features miss (ZP26 is the baseline)
region_SP15	SP15 market-zone control	level difference of SP15 counties the other features miss

How to Read the Scorecards

Lower nll, mae and rmse is better, higher auc_... and count_r2; brier_... is benchmarked against the base rate \(r\), which scores \(r(1-r)\). All are computed on held-out county-days - the final stretch of dates, withheld before fitting.

county_fit = fit_negative_binomial_model(
    county_panel,
    feature_cols=COUNTY_DAILY_FEATURES,
    high_count_quantile=0.90,
    maxiter=300,
)

county_metrics = pd.Series(county_fit.metrics).to_frame("value")
html_table(
    pl.from_pandas(
        county_metrics.loc[
            [
                "alpha",
                "n_train",
                "n_test",
                "test_mae",
                "test_rmse",
                "test_nll",
                "test_auc_outage",
                "test_auc_high_outage",
                "test_brier_high_outage",
                "test_high_count_threshold",
                "test_high_outage_share",
                "test_max_predicted",
            ]
        ].reset_index()
    )
)

12 rows × 2 columns

index	value
alpha	0.01
n_train	78,138.00
n_test	19,548.00
test_mae	0.74
test_rmse	1.43
test_nll	1.02
test_auc_outage	0.88
test_auc_high_outage	0.98
test_brier_high_outage	0.03
test_high_count_threshold	5.00
test_high_outage_share	0.08
test_max_predicted	51.35

The same held-out county-days, scored against two do-nothing forecasts:

test_fitted = county_fit.fitted_df[county_fit.fitted_df["sample"].eq("test")]
train_fitted = county_fit.fitted_df[county_fit.fitted_df["sample"].eq("train")]
train_mean = float(train_fitted["Y"].mean())
train_high_rate = float((train_fitted["Y"] >= county_fit.high_count_threshold).mean())
test_high_actual = (test_fitted["Y"] >= county_fit.high_count_threshold).astype(float)

baseline_table = pd.DataFrame(
    [
        {
            "forecast": "NB2 model",
            "test MAE (outages)": (test_fitted["Y"] - test_fitted["mu_hat"]).abs().mean(),
            "test Brier (high-outage)": ((test_high_actual - test_fitted["p_high_outage_hat"]) ** 2).mean(),
        },
        {
            "forecast": "always predict the training mean",
            "test MAE (outages)": (test_fitted["Y"] - train_mean).abs().mean(),
            "test Brier (high-outage)": ((test_high_actual - train_high_rate) ** 2).mean(),
        },
        {
            "forecast": "always predict zero outages",
            "test MAE (outages)": test_fitted["Y"].abs().mean(),
            "test Brier (high-outage)": (test_high_actual**2).mean(),
        },
    ]
).round(4)
html_table(pl.from_pandas(baseline_table))

3 rows × 3 columns

forecast	test MAE (outages)	test Brier (high-outage)
NB2 model	0.745	0.025
always predict the training mean	1.790	0.078
always predict zero outages	1.274	0.084

Ranking is the strength: AUC 0.88 any-outage, 0.98 high-outage, and a high-outage Brier under a third of the training-mean forecast’s. Point forecasts are modest - an average miss of 0.74 outages against 1.27 for always-zero, the harder benchmark because six in ten county-days are zero - so use the ranking, not the exact county count.

Reading the Coefficient Chart and Table

coefficient is the fitted weight \(\beta_k\) per standardized step, reported multiplicatively as incidence_rate_ratio \(= e^{\beta}\) with its 95% range; conf_low / conf_high are the chart’s whiskers, and a range containing zero means the direction is unsettled. With roughly 78,000 training rows, read p_value as direction, not size.

county_coef = county_fit.coefficients.copy()
coef_view = county_coef[county_coef["feature"].ne("const")].sort_values("coefficient")

coef_src = ColumnDataSource(
    {
        "feature": coef_view["feature"].tolist(),
        "coefficient": coef_view["coefficient"].to_numpy(),
        "conf_low": coef_view["conf_low"].to_numpy(),
        "conf_high": coef_view["conf_high"].to_numpy(),
    }
)
fig = figure(
    title="County-Day NB2 Coefficients, Standardized Features",
    y_range=coef_view["feature"].tolist(),
    x_axis_label="Log-Count Coefficient",
    height=460,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.segment(
    "conf_low", "feature", "conf_high", "feature",
    source=coef_src, color=SERIES[0], line_width=2, alpha=0.45,
)
coef_points = fig.scatter("coefficient", "feature", source=coef_src, size=7, color=SERIES[0])
fig.add_layout(Span(location=0, dimension="height", line_color="black", line_width=1))
fig.add_tools(
    HoverTool(
        renderers=[coef_points],
        tooltips=[
            ("feature", "@feature"),
            ("coefficient", "@coefficient{0.000}"),
            ("conf low", "@conf_low{0.000}"),
            ("conf high", "@conf_high{0.000}"),
        ],
    )
)
show(style(fig))

html_table(pl.from_pandas(county_coef.sort_values("p_value").head(12)))

12 rows × 11 columns

feature	coefficient	std_error	z_stat	conf_low	conf_high	p_value	abs_coefficient	incidence_rate_ratio	irr_low	irr_high
const	−0.590	0.006	−102.54	−0.602	−0.579	0.000	0.590	0.554	0.548	0.560
log1p_Y_lag1	0.779	0.007	114.41	0.766	0.793	0.000	0.779	2.180	2.151	2.209
region_NP15	0.227	0.006	40.16	0.216	0.238	0.000	0.227	1.255	1.241	1.269
log_county_nameplate_capacity_mw	0.345	0.011	30.66	0.323	0.367	0.000	0.345	1.413	1.382	1.444
log_county_plant_count	0.215	0.010	21.93	0.196	0.234	0.000	0.215	1.240	1.216	1.264
region_SP15	0.058	0.004	16.54	0.051	0.065	0.000	0.058	1.060	1.053	1.067
log1p_forced_start_mw_sum_lag1	0.094	0.007	13.66	0.081	0.108	0.000	0.094	1.099	1.084	1.114
weekend	−0.045	0.004	−12.44	−0.052	−0.038	0.000	0.045	0.956	0.949	0.963
temperature_f_std	−0.018	0.003	−6.04	−0.024	−0.012	0.000	0.018	0.982	0.976	0.988
county_weather_basis_risk_flag	0.019	0.003	5.57	0.012	0.026	0.000	0.019	1.019	1.013	1.026
log_system_load_mean	0.026	0.006	4.34	0.014	0.037	0.000	0.026	1.026	1.014	1.038
log1p_planned_start_count_lag1	−0.004	0.002	−2.00	−0.007	−0.000	0.045	0.004	0.996	0.993	1.000

The leading signals are persistence, the county generation footprint, and the zone effects. Degree-day weights are indistinguishable from zero on five summers (cdd65 at +0.005 and hdd65 at +0.007, both whiskers spanning zero); what survives of weather is the temperature-spread term and the multi-proxy flag, both small. Notebook 15.1 shows the form weather actually takes - heat interacted with installed capacity, not a uniform level shift.

4. Specification Tests

Poisson is the benchmark; NB2 adds the dispersion parameter the heavier tails call for.

county_nb_vs_poisson = poisson_negative_binomial_comparison(
    county_panel,
    feature_cols=COUNTY_DAILY_FEATURES,
    high_count_quantile=0.90,
    maxiter=300,
)
html_table(pl.from_pandas(county_nb_vs_poisson))

2 rows × 10 columns

model	alpha	train_llf	aic	test_mae	test_rmse	test_nll	lr_chi2_vs_poisson	lr_df	lr_p_value_boundary_mixture
Poisson	0.000	−79,152.40	158,336.80	0.745	1.416	1.019	—	—	—
Negative binomial NB2	0.013	−79,124.85	158,283.70	0.745	1.426	1.015	55.10	1.000	0.000

lr_chi2_vs_poisson scores 55 against a benchmark near 4, and NB2 wins on likelihood, AIC and held-out NLL. But the fitted dispersion is only 0.013 and the two models nearly coincide (MAE identical to three decimals, RMSE within 0.010, NLL within 0.005 nats) because the features absorb most of the dispersion. Notebook 15.3 tests the variance shape directly and settles on the linear NB1 form.

county_feature_blocks = {
    "outage state": [
        "log1p_Y_lag1",
        "log1p_planned_start_count_lag1",
        "log1p_forced_start_mw_sum_lag1",
    ],
    "size": ["log_county_nameplate_capacity_mw", "log_county_plant_count"],
    "weather": ["cdd65", "hdd65", "temperature_f_std", "county_weather_basis_risk_flag"],
    "load": ["log_system_load_mean"],
    "calendar": ["weekend", "month_sin", "month_cos"],
    "zone effects": ["region_NP15", "region_SP15"],
}

county_block_tests = likelihood_ratio_block_tests(
    county_panel,
    feature_cols=COUNTY_DAILY_FEATURES,
    feature_blocks=county_feature_blocks,
    high_count_quantile=0.90,
    maxiter=300,
)
html_table(pl.from_pandas(county_block_tests))

6 rows × 10 columns

block	dropped_features	df	full_train_llf	reduced_train_llf	lr_chi2	p_value	full_test_nll	reduced_test_nll	delta_test_nll_reduced_minus_full
outage state	log1p_Y_lag1, log1p_planned_start_count_lag1, log1p_forced_start_mw_sum_lag1	3	−79,124.85	−94,932.81	31,615.92	0.000	1.015	1.183	0.167
size	log_county_nameplate_capacity_mw, log_county_plant_count	2	−79,124.85	−81,910.15	5,570.59	0.000	1.015	1.055	0.040
zone effects	region_NP15, region_SP15	2	−79,124.85	−79,957.56	1,665.41	0.000	1.015	1.019	0.004
calendar	weekend, month_sin, month_cos	3	−79,124.85	−79,209.06	168.42	0.000	1.015	1.017	0.002
weather	cdd65, hdd65, temperature_f_std, county_weather_basis_risk_flag	4	−79,124.85	−79,148.55	47.40	0.000	1.015	1.016	0.000
load	log_system_load_mean	1	−79,124.85	−79,134.22	18.73	0.000	1.015	1.015	−0.000

lr_chi2 is how much fit collapses when a block is deleted; delta_test_nll_reduced_minus_full is the held-out version. Outage persistence dominates (LR 31,616, 0.167 held-out nats), county size is second (5,571, 0.040 nats), and zone, calendar and weather register modestly (1,665, 168, 47) with weather’s held-out contribution rounding to zero. Load is significant yet worthless out of sample.

5. Purged Walk-Forward Backtest

The model is refit four times, each fold trained only on dates before the days it scores, with the day before each fold purged and one day after embargoed so the lagged variables cannot leak.

county_backtest = rolling_backtest_negative_binomial_model(
    county_panel,
    feature_cols=COUNTY_DAILY_FEATURES,
    n_splits=4,
    initial_train_fraction=0.60,
    purge_gap="1D",
    embargo_gap="1D",
    high_count_quantile=0.90,
    maxiter=300,
)

html_table(pl.from_pandas(county_backtest.overall_metrics))

1 rows × 22 columns

backtest	rows	mean_actual	mean_predicted	zero_share	outage_share	high_count_threshold	high_outage_share	mae	median_ae	p95_ae	rmse	nll	count_r2	max_actual	max_predicted	mean_p_outage	auc_outage	brier_outage	mean_p_high_outage	auc_high_outage	brier_high_outage
all folds	38,934	1.366	1.393	0.603	0.397	5.000	0.091	0.753	0.335	2.690	1.407	1.029	0.755	46.00	51.32	0.434	0.891	0.129	0.087	0.979	0.027

Walk-forward barely moves off the section 3 single split: AUC 0.89 any-outage and 0.98 high-outage, Brier 0.13 against a 0.24 base-rate benchmark (event rate 0.40) and 0.027 against 0.083 (event rate 0.09), mean predicted 1.39 against 1.36 actual. Counts remain the weak spot - count_r2 0.76 and a 95th-percentile miss near 2.7 outages.

html_table(pl.from_pandas(county_backtest.region_metrics))

3 rows × 22 columns

region	rows	mean_actual	mean_predicted	zero_share	outage_share	high_count_threshold	high_outage_share	mae	median_ae	p95_ae	rmse	nll	count_r2	max_actual	max_predicted	mean_p_outage	auc_outage	brier_outage	mean_p_high_outage	auc_high_outage	brier_high_outage
NP15	21,630	0.812	0.849	0.681	0.319	5.000	0.043	0.550	0.291	1.775	0.865	0.834	0.708	14.00	13.89	0.394	0.894	0.120	0.040	0.975	0.016
SP15	5,768	3.289	3.335	0.286	0.714	5.000	0.287	1.409	0.708	4.806	2.253	1.736	0.653	24.00	26.46	0.702	0.928	0.100	0.283	0.942	0.086
ZP26	11,536	1.444	1.443	0.615	0.385	5.000	0.085	0.807	0.349	2.922	1.656	1.042	0.784	46.00	51.32	0.374	0.840	0.159	0.078	0.986	0.019

backtest_regions = county_backtest.region_metrics["region"].tolist()

fig = figure(
    title="Walk-Forward AUC: Any County Outage",
    x_range=backtest_regions,
    y_range=Range1d(0.5, 1.0),
    y_axis_label="AUC",
    height=340,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.vbar(
    x=backtest_regions,
    top=county_backtest.region_metrics["auc_outage"].tolist(),
    width=0.6,
    color=SERIES[0],
)
fig.add_tools(HoverTool(tooltips=[("region", "@x"), ("AUC", "@top{0.000}")]))
show(style(fig))

fig = figure(
    title="Walk-Forward AUC: High-Outage County-Day",
    x_range=backtest_regions,
    y_range=Range1d(0.5, 1.0),
    y_axis_label="AUC",
    height=340,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
fig.vbar(
    x=backtest_regions,
    top=county_backtest.region_metrics["auc_high_outage"].tolist(),
    width=0.6,
    color=SERIES[1],
)
fig.add_tools(HoverTool(tooltips=[("region", "@x"), ("AUC", "@top{0.000}")]))
show(style(fig))

ZP26 ranks any-outage days worst but high-outage county-days well. Calibration by zone: ZP26 1.44 predicted against 1.44 actual, SP15 1% high (3.33 vs 3.28), NP15 5% high (0.85 vs 0.81) - over-prediction would inflate the conditioning variable the rollup feeds to pricing.

6. Market-Zone Rollup for Spark-Spread Inputs

County predictions roll up into market-zone features: a summed county intensity, and a weighted version that is a size-adjusted signal joinable to NP15/SP15/ZP26 spreads.

For market zone \(z\) on day \(d\) with counties \(c \in z\), the rollup builds

\[\mu^{(z)}_{d} = \sum_{c \in z} \mu_{c,d}, \qquad \tilde{\mu}^{(z)}_{d} = \sum_{c \in z} w_{c}\,\mu_{c,d}, \qquad \Pr\big[\text{any high-outage county in } z\big] = 1 - \prod_{c \in z}\big(1 - p^{\text{high}}_{c,d}\big)\]

where \(w_{c}\) is the county’s share of the zone, by capacity or by resource count. The any-high-outage probability assumes conditional independence, so prob_any_county_high_outage is a ranking feature, not a calibrated probability.

county_zone_rollup = rollup_county_predictions_to_market_zone(
    county_fit.fitted_df,
    county_zone_weights,
)

rollup_preview_cols = [
    "date",
    "region",
    "sample",
    "county_count",
    "actual_forced_count",
    "county_mu_sum",
    "prob_any_county_high_outage",
    "capacity_weighted_county_mu",
    "resource_weighted_county_mu",
]
html_table(pl.from_pandas(county_zone_rollup.loc[:, rollup_preview_cols].tail(12)))

12 rows × 9 columns

date	region	sample	county_count	actual_forced_count	county_mu_sum	prob_any_county_high_outage	capacity_weighted_county_mu	resource_weighted_county_mu
2026-06-27	NP15	test	30	19.00	22.53	0.762	2.527	1.310
2026-06-27	SP15	test	8	17.00	19.51	0.977	3.589	3.572
2026-06-27	ZP26	test	16	23.00	18.94	0.979	5.904	5.511
2026-06-28	NP15	test	30	17.00	20.13	0.668	2.357	1.171
2026-06-28	SP15	test	8	19.00	16.70	0.917	3.032	3.024
2026-06-28	ZP26	test	16	20.00	19.70	0.977	5.986	5.549
2026-06-29	NP15	test	30	19.00	22.16	0.765	2.662	1.330
2026-06-29	SP15	test	8	24.00	21.42	0.993	4.137	4.126
2026-06-29	ZP26	test	16	17.00	18.92	0.949	5.339	4.970
2026-06-30	NP15	test	30	18.00	23.87	0.833	2.810	1.474
2026-06-30	SP15	test	8	23.00	25.16	0.997	4.581	4.532
2026-06-30	ZP26	test	16	30.00	18.50	0.974	5.862	5.463

plot_df = county_zone_rollup[county_zone_rollup["sample"].eq("test")].copy()

fig = figure(
    title="County Model Rolled Up to Market Zones",
    x_axis_type="datetime",
    y_axis_label="Forced Outage Starts",
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
for color, (region, group) in zip(SERIES, plot_df.groupby("region")):
    fig.line(
        pd.to_datetime(group["date"]), group["actual_forced_count"],
        color=color, alpha=0.75, legend_label=f"{region} actual",
    )
    fig.line(
        pd.to_datetime(group["date"]), group["county_mu_sum"],
        color=color, line_dash="dashed", legend_label=f"{region} fitted",
    )
fig.legend.location = "top_left"
show(style(fig))

fig = figure(
    title="Probability of Any High-Outage County Within Zone",
    x_axis_type="datetime",
    y_range=Range1d(0, 1),
    y_axis_label="Probability",
    height=380,
    sizing_mode="stretch_width",
    toolbar_location="above",
)
for color, (region, group) in zip(SERIES, plot_df.groupby("region")):
    fig.line(
        pd.to_datetime(group["date"]), group["prob_any_county_high_outage"],
        color=color, legend_label=region,
    )
fig.legend.location = "top_left"
show(style(fig))

The rollup produces market-zone daily variables joinable to spark spreads, DART spreads, or quanto conditioning terms.

Findings

• County-day keeps the market-zone anchor while giving weather a more local exposure channel.

• Lagged outage activity and county generation footprint dominate; plain degree-day levels carry no weight once they are in.

• The rollup produces daily NP15/SP15/ZP26 outage-intensity features. The model ranks and assigns probabilities; it does not forecast exact counts.

• Notebook 15 searches the specification (NB1, CDD x capacity); notebook 16 locates the heat signal in derate depth by reason code.