CAISO Negative Binomial Outage Model#
Authors: Claire Kim · Brian Nguyen ·
Created: 2026-08-21 ·
Tags: #Experiment #Outages #Count-Model
The response is the daily count of new forced outage starts, modelled on a market-zone by day panel for NP15, SP15, and ZP26 with a negative-binomial count model. On the extended five-year panel the fitted dispersion comes out at zero, so the estimated model is effectively Poisson; the raw zone-day counts still show variance-to-mean of 1.4 to 2.8, which is why the NB2 family is fitted rather than assumed.
The model outputs a conditional expected outage count, mu, and a high-outage probability by market zone and day.
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.models import Range1d, Span
from bokeh.plotting import figure
from bokeh.resources import INLINE
from statsmodels.tools.sm_exceptions import (
ConvergenceWarning,
HessianInversionWarning,
)
pd.set_option("display.max_columns", 140)
pd.set_option("display.width", 180)
warnings.filterwarnings("ignore", category=FutureWarning)
# The daily and 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)
ROOT = Path.cwd()
if ROOT.name == "notebooks":
ROOT = ROOT.parent
if str(ROOT / "src") not in sys.path:
sys.path.insert(0, str(ROOT / "src"))
from caiso_outage_nb_pipeline import (
DEFAULT_FEATURES,
MARKET_ZONE_FEATURES,
WEATHER_SENSITIVE_TARGET,
build_caiso_market_zone_nb_panel,
build_caiso_nb_panel,
poisson_negative_binomial_comparison,
likelihood_ratio_block_tests,
fit_logistic_probability_model,
binary_model_performance_by_group,
feature_correlation_table,
fit_negative_binomial_model,
model_performance_by_group,
overdispersion_summary,
resource_mapping_coverage_summary,
rolling_backtest_negative_binomial_model,
top_nature_of_work,
)
from chart_theme import SECONDARY, SERIES, html_table, style
DATA_DIR = ROOT / "data"
PROCESSED_DIR = DATA_DIR / "processed"
PROCESSED_DIR.mkdir(parents=True, exist_ok=True)
output_notebook(INLINE, hide_banner=True)
1. Source Tables#
Inputs come from CAISO/GridStatus market caches, Open-Meteo weather aggregates, active outage MW, load, and the plant-level generator-weather inventory.
source_paths = {
"CAISO curtailed generator report": DATA_DIR / "opensource" / "caiso_curtailed_non_operational_generator_report.parquet",
"CAISO hourly TAC load": DATA_DIR / "opensource" / "caiso_load_hourly.parquet",
"CAISO active outage MW, hourly": DATA_DIR / "opensource" / "caiso_outages_aggregated_hourly_derived.parquet",
"Open-Meteo market-zone weather": DATA_DIR / "processed" / "caiso_weather_daily_region.parquet",
"Generator-weather inventory": DATA_DIR / "processed" / "caiso_generator_weather_imputed.parquet",
"Weather-zone crosswalk": DATA_DIR / "processed" / "weather_zone_crosswalk.parquet",
}
source_summary = []
for name, path in source_paths.items():
if path.exists():
df = pd.read_parquet(path)
source_summary.append({
"source": name,
"path": str(path.relative_to(ROOT)),
"rows": len(df),
"columns": len(df.columns),
"available": True,
})
else:
source_summary.append({
"source": name,
"path": str(path.relative_to(ROOT)),
"rows": np.nan,
"columns": np.nan,
"available": False,
})
display(html_table(pl.from_pandas(pd.DataFrame(source_summary)), kind="reference"))
reference · 6 rows × 5 columns
| source | path | rows | columns | available |
|---|---|---|---|---|
| CAISO curtailed generator report | data\opensource\caiso_curtailed_non_operational_generator_report.parquet | 3,655,547 | 11 | True |
| CAISO hourly TAC load | data\opensource\caiso_load_hourly.parquet | 1,944,693 | 4 | True |
| CAISO active outage MW, hourly | data\opensource\caiso_outages_aggregated_hourly_derived.parquet | 45,479 | 4 | True |
| Open-Meteo market-zone weather | data\processed\caiso_weather_daily_region.parquet | 8,032 | 29 | True |
| Generator-weather inventory | data\processed\caiso_generator_weather_imputed.parquet | 1,661 | 53 | True |
| Weather-zone crosswalk | data\processed\weather_zone_crosswalk.parquet | 7 | 11 | True |
2. Resource-to-Market-Zone Mapping#
CAISO outage records use resource labels, while the plant map uses EIA plant names. The crosswalk matches on name with a capacity sanity check, assigning a zone only where the plant match is strong.
market_panel, mapped_events, resource_crosswalk = build_caiso_market_zone_nb_panel(DATA_DIR)
coverage = resource_mapping_coverage_summary(
mapped_events,
mapped_events,
resource_crosswalk,
start_date=market_panel["date"].min(),
end_date=market_panel["date"].max(),
)
coverage_display = coverage[
[
"outage_type",
"events",
"matched_events",
"event_coverage",
"unique_resources",
"matched_unique_resources",
"resource_coverage",
"crosswalk_rows",
]
].copy()
display(html_table(pl.from_pandas(coverage_display.round(3)), kind="reference"))
region_breakout = (
mapped_events[mapped_events["region"].isin(["NP15", "SP15", "ZP26"])]
.query("outage_type == 'FORCED'")
.groupby("region")
.agg(
matched_forced_events=("outage_mrid_str", "count"),
unique_resources=("resource_name", "nunique"),
total_start_mw=("max_curtailment_mw", "sum"),
)
.reset_index()
)
display(html_table(pl.from_pandas(region_breakout.round(2))))
crosswalk_preview_cols = [
"resource_id",
"resource_name",
"matched_plant_name",
"region",
"nameplate_capacity_mw",
"resource_max_curtailment_mw",
"match_score",
"match_overlap_tokens",
]
display(html_table(pl.from_pandas(resource_crosswalk[crosswalk_preview_cols].head(12)), kind="reference"))
reference · 3 rows × 8 columns
| outage_type | events | matched_events | event_coverage | unique_resources | matched_unique_resources | resource_coverage | crosswalk_rows |
|---|---|---|---|---|---|---|---|
| FORCED | 243,709 | 147,616 | 0.606 | 1,779 | 740 | 0.416 | 754 |
| PLANNED | 16,736 | 8,812 | 0.527 | 1,245 | 630 | 0.506 | 754 |
| ALL | 260,445 | 156,428 | 0.601 | 1,838 | 747 | 0.406 | 754 |
3 rows × 4 columns
| region | matched_forced_events | unique_resources | total_start_mw |
|---|---|---|---|
| NP15 | 46,637 | 166 | 1,717,829.09 |
| SP15 | 56,563 | 294 | 3,407,634.54 |
| ZP26 | 48,880 | 288 | 1,493,131.96 |
reference · 12 rows × 8 columns
| resource_id | resource_name | matched_plant_name | region | nameplate_capacity_mw | resource_max_curtailment_mw | match_score | match_overlap_tokens |
|---|---|---|---|---|---|---|---|
| LMEC_1_PL1X3 | Los Medanos Energy Center AGGREGATE | Los Medanos Energy Center | NP15 | 678.30 | 580.00 | 100.00 | LOS, MEDANOS |
| OTMESA_2_PL1X3 | OTAY MESA ENERGY CENTER | Otay Mesa Generating Project | SP15 | 710.70 | 603.68 | 100.00 | MESA, OTAY |
| PIOPIC_2_CTG1 | Pio Pico Unit 1 | Pio Pico Energy Center | SP15 | 395.40 | 111.30 | 95.00 | PICO, PIO |
| PIOPIC_2_CTG2 | Pio Pico Unit 2 | Pio Pico Energy Center | SP15 | 395.40 | 112.70 | 95.00 | PICO, PIO |
| METEC_2_PL1X3 | Metcalf Energy Center | Metcalf Energy Center | NP15 | 635.00 | 597.05 | 100.00 | METCALF |
| PIOPIC_2_CTG3 | Pio Pico Unit 3 | Pio Pico Energy Center | SP15 | 395.40 | 112.00 | 95.00 | PICO, PIO |
| GATWAY_2_PL1X3 | GATEWAY GENERATING STATION | Gateway Generating Station | NP15 | 619.70 | 585.00 | 100.00 | GATEWAY |
| COCOPP_2_CTG3 | Marsh Landing 3 | Marsh Landing Generating Station | NP15 | 828.00 | 208.96 | 95.00 | LANDING, MARSH |
| SENTNL_2_CTG6 | Sentinel Unit 6 | Sentinel Energy Center, LLC | SP15 | 800.00 | 105.00 | 95.00 | SENTINEL |
| SENTNL_2_CTG8 | Sentinel Unit 8 | Sentinel Energy Center, LLC | SP15 | 800.00 | 106.85 | 95.00 | SENTINEL |
| ELKHIL_2_PL1X3 | ELK HILLS COMBINED CYCLE (AGGREGATE) | Elk Hills Power LLC | ZP26 | 623.00 | 551.70 | 95.00 | ELK, HILLS |
| SENTNL_2_CTG4 | Sentinel Unit 4 | Sentinel Energy Center, LLC | SP15 | 800.00 | 106.55 | 95.00 | SENTINEL |
The zonal panel is a matched-resource sample, not a forced assignment of every CAISO outage, so the unmapped share stays visible as sampling error.
unmatched_forced = (
mapped_events
.query("outage_type == 'FORCED' and region.isna()")
.groupby("resource_name")
.agg(
forced_events=("outage_mrid_str", "count"),
max_curtailment_mw=("max_curtailment_mw", "max"),
common_reason=("nature_of_work", lambda s: s.mode().iloc[0] if not s.mode().empty else np.nan),
)
.sort_values("forced_events", ascending=False)
.head(12)
.reset_index()
)
display(html_table(pl.from_pandas(unmatched_forced)))
12 rows × 4 columns
| resource_name | forced_events | max_curtailment_mw | common_reason |
|---|---|---|---|
| KERKHOFF PH 2 UNIT #1 | 1,951 | 153.90 | AMBIENT_NOT_DUE_TO_TEMP |
| BCTSYS_5_PWXDYN | 1,855 | 428.00 | TRANSITIONAL_LIMITATION |
| LODI GAS TURBINE | 1,830 | 25.00 | AMBIENT_DUE_TO_TEMP |
| Colusa Generating Station | 1,808 | 641.00 | AMBIENT_DUE_TO_TEMP |
| PIT PH 6 UNIT 1 | 1,788 | 39.00 | AMBIENT_NOT_DUE_TO_TEMP |
| PIT PH 6 UNIT 2 | 1,785 | 40.00 | AMBIENT_NOT_DUE_TO_TEMP |
| DELTA ENERGY CENTER AGGREGATE | 1,757 | 880.00 | AMBIENT_DUE_TO_TEMP |
| PIT PH 7 UNIT 1 | 1,749 | 55.70 | AMBIENT_NOT_DUE_TO_TEMP |
| Lodi Energy Center | 1,606 | 302.58 | AMBIENT_DUE_TO_TEMP |
| Barre Peaker | 1,565 | 49.00 | AMBIENT_DUE_TO_TEMP |
| PIT PH 7 UNIT 2 | 1,511 | 54.60 | AMBIENT_NOT_DUE_TO_TEMP |
| GILROY COGEN AGGREGATE | 1,489 | 120.00 | AMBIENT_DUE_TO_TEMP |
The largest unmatched resources are mostly hydro aggregates, pseudo-CISO shares, and resources without a clean plant-name counterpart.
3. Market-Zone Daily Panel#
Each row is one market zone on one Pacific calendar day. The target Y is the matched count of new forced outage starts in that zone. Weather is joined at the same region x date level through the Open-Meteo market-zone aggregates.
panel_preview_cols = [
"date",
"region",
"Y",
"planned_start_count",
"weather_sensitive_start_count",
"fuel_insufficiency_start_count",
"system_load_mw_mean",
"temperature_f",
"hdd65",
"cdd65",
"temperature_f_std",
"market_weather_basis_risk_flag",
]
display(html_table(pl.from_pandas(market_panel[panel_preview_cols].head(9))))
panel_coverage = pd.Series({
"first_model_date": market_panel["date"].min(),
"last_model_date": market_panel["date"].max(),
"region_day_rows": len(market_panel),
"regions": ", ".join(sorted(market_panel["region"].unique())),
"target": "matched new forced outage starts per market-zone day",
})
display(html_table(pl.from_pandas(panel_coverage.to_frame("value").reset_index().astype(str)), kind="reference"))
9 rows × 12 columns
| date | region | Y | planned_start_count | weather_sensitive_start_count | fuel_insufficiency_start_count | system_load_mw_mean | temperature_f | hdd65 | cdd65 | temperature_f_std | market_weather_basis_risk_flag |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2021-06-17 | NP15 | 38 | 2 | 19 | 0 | 32,513.46 | 86.14 | 0.000 | 21.14 | 5.481 | 1 |
| 2021-06-18 | NP15 | 37 | 0 | 17 | 0 | 32,603.62 | 84.35 | 0.000 | 19.35 | 3.929 | 1 |
| 2021-06-19 | NP15 | 31 | 3 | 16 | 0 | 29,559.54 | 76.84 | 0.000 | 11.84 | 3.833 | 1 |
| 2021-06-20 | NP15 | 25 | 0 | 14 | 0 | 27,179.83 | 73.91 | 0.000 | 8.91 | 3.135 | 1 |
| 2021-06-21 | NP15 | 34 | 4 | 14 | 0 | 28,705.00 | 71.81 | 0.000 | 6.81 | 1.810 | 1 |
| 2021-06-22 | NP15 | 35 | 4 | 14 | 0 | 28,442.33 | 71.96 | 0.000 | 6.96 | 1.873 | 1 |
| 2021-06-23 | NP15 | 28 | 2 | 15 | 0 | 28,560.62 | 70.43 | 0.000 | 5.43 | 3.685 | 1 |
| 2021-06-24 | NP15 | 31 | 1 | 15 | 0 | 27,672.58 | 68.41 | 0.000 | 3.41 | 5.169 | 1 |
| 2021-06-25 | NP15 | 32 | 2 | 17 | 0 | 27,887.83 | 69.78 | 0.000 | 4.78 | 5.163 | 1 |
reference · 5 rows × 2 columns
| index | value |
|---|---|
| first_model_date | 2021-06-17 |
| last_model_date | 2026-06-30 |
| region_day_rows | 5520 |
| regions | NP15, SP15, ZP26 |
| target | matched new forced outage starts per market-zone day |
market_weather_basis_risk_flag mainly identifies that NP15 spans multiple weather-zone proxies, so the specification uses zone fixed effects and temperature dispersion instead of the flag.
4. Overdispersion by Market Zone#
Outage intensity is heterogeneous across zones, so the table checks whether each zone-day count has variance above the mean.
zone_dispersion = (
market_panel.groupby("region")["Y"]
.agg(
rows="size",
mean="mean",
variance="var",
zero_share=lambda s: (s == 0).mean(),
p90=lambda s: s.quantile(0.90),
max="max",
)
)
zone_dispersion["variance_to_mean"] = zone_dispersion["variance"] / zone_dispersion["mean"]
display(html_table(pl.from_pandas(zone_dispersion.round(3).reset_index())))
regions = sorted(market_panel["region"].unique())
bin_edges = np.histogram_bin_edges(market_panel["Y"].to_numpy(), bins=25)
fig = figure(
title="Forced Outage Starts by Market Zone",
x_axis_label="zone-day count",
y_axis_label="zone-days",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for region, color in zip(regions, SERIES):
counts, _ = np.histogram(
market_panel.loc[market_panel["region"] == region, "Y"].to_numpy(), bins=bin_edges
)
fig.quad(
top=counts,
bottom=0,
left=bin_edges[:-1],
right=bin_edges[1:],
fill_alpha=0,
line_color=color,
line_width=1.8,
legend_label=region,
)
show(style(fig))
fig = figure(
title="Matched Forced Starts Over Time",
x_axis_type="datetime",
y_axis_label="daily count",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for region, color in zip(regions, SERIES):
zone_rows = market_panel[market_panel["region"] == region].sort_values("date")
fig.line(
pd.to_datetime(zone_rows["date"]), zone_rows["Y"],
color=color, line_width=1.2, legend_label=region,
)
fig.legend.location = "top_left"
show(style(fig))
3 rows × 8 columns
| region | rows | mean | variance | zero_share | p90 | max | variance_to_mean |
|---|---|---|---|---|---|---|---|
| NP15 | 1,840 | 24.62 | 33.94 | 0.000 | 32.00 | 52 | 1.379 |
| SP15 | 1,840 | 29.75 | 84.13 | 0.000 | 42.00 | 65 | 2.828 |
| ZP26 | 1,840 | 25.86 | 48.24 | 0.000 | 34.00 | 63 | 1.866 |
The zone-day counts are dense, but the scale is more useful than the all-CAISO count. SP15 and ZP26 show clearer overdispersion than NP15 in this matched sample.
5. Conditioning Set#
The feature vector follows the outage-count workflow: weather stress, load, lagged outage counts, active curtailment MW, calendar effects, and zone fixed effects.
feature_roles = {
"cdd65": "Cooling stress from Open-Meteo market-zone weather",
"hdd65": "Heating stress from Open-Meteo market-zone weather",
"log_system_load_mean": "CAISO system load level; zone load is not yet mapped to NP15/SP15/ZP26",
"Y_lag1": "Prior-day matched forced outage starts in the same market zone",
"planned_start_count_lag1": "Prior-day planned outage starts in the same market zone",
"active_forced_mw_mean_lag1": "Prior-day active forced curtailment MW at CAISO system level",
"temperature_f_std": "Within-zone temperature dispersion across representative weather points",
"weekend": "Weekend indicator",
"month_sin": "Seasonal sine term",
"month_cos": "Seasonal cosine term",
"region_NP15": "NP15 fixed effect; ZP26 is the omitted zone",
"region_SP15": "SP15 fixed effect; ZP26 is the omitted zone",
}
feature_table = pd.DataFrame({"feature": MARKET_ZONE_FEATURES})
feature_table["role"] = feature_table["feature"].map(feature_roles)
display(html_table(pl.from_pandas(feature_table), kind="reference"))
display(html_table(pl.from_pandas(feature_correlation_table(market_panel, target_col="Y", feature_cols=MARKET_ZONE_FEATURES).round(3))))
reference · 12 rows × 2 columns
| feature | role |
|---|---|
| cdd65 | Cooling stress from Open-Meteo market-zone weather |
| hdd65 | Heating stress from Open-Meteo market-zone weather |
| log_system_load_mean | CAISO system load level; zone load is not yet mapped to NP15/SP15/ZP26 |
| Y_lag1 | Prior-day matched forced outage starts in the same market zone |
| planned_start_count_lag1 | Prior-day planned outage starts in the same market zone |
| active_forced_mw_mean_lag1 | Prior-day active forced curtailment MW at CAISO system level |
| temperature_f_std | Within-zone temperature dispersion across representative weather points |
| weekend | Weekend indicator |
| month_sin | Seasonal sine term |
| month_cos | Seasonal cosine term |
| region_NP15 | NP15 fixed effect; ZP26 is the omitted zone |
| region_SP15 | SP15 fixed effect; ZP26 is the omitted zone |
12 rows × 2 columns
| feature | pearson_corr_with_target |
|---|---|
| Y_lag1 | 0.639 |
| log_system_load_mean | 0.464 |
| month_sin | −0.334 |
| hdd65 | −0.279 |
| cdd65 | 0.277 |
| region_SP15 | 0.274 |
| weekend | −0.258 |
| month_cos | −0.212 |
| region_NP15 | −0.193 |
| active_forced_mw_mean_lag1 | −0.169 |
| temperature_f_std | 0.088 |
| planned_start_count_lag1 | −0.013 |
Weather enters alongside load, lagged outage state, and geography, not as a standalone explanation for forced outages.
6. Market-Zone NB2 Fit#
The fitted specification is:
\[
Y_{z,t} \mid X_{z,t} \sim \mathrm{NB2}(\mu_{z,t},\ \alpha), \qquad
\log(\mu_{z,t}) = X_{z,t}\beta, \qquad
\mathrm{Var}[Y_{z,t} \mid X_{z,t}] = \mu_{z,t} + \alpha\,\mu_{z,t}^{2}
\]
The split is chronological: the final 20% of observations are held out.
zone_fit = fit_negative_binomial_model(
market_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
)
metrics_to_show = [
"alpha",
"n_train",
"n_test",
"train_mean_actual",
"test_mean_actual",
"train_outage_share",
"test_outage_share",
"train_mae",
"test_mae",
"train_median_ae",
"test_median_ae",
"train_p95_ae",
"test_p95_ae",
"train_rmse",
"test_rmse",
"train_nll",
"test_nll",
"train_max_actual",
"test_max_actual",
"train_max_predicted",
"test_max_predicted",
"train_high_count_threshold",
"test_high_count_threshold",
"train_auc_outage",
"test_auc_outage",
"train_auc_high_outage",
"test_auc_high_outage",
"train_brier_high_outage",
"test_brier_high_outage",
]
zone_metrics_df = pd.Series(zone_fit.metrics).loc[metrics_to_show].rename("value").to_frame()
display(html_table(pl.from_pandas(zone_metrics_df.round(4).reset_index())))
zone_oos_by_region = model_performance_by_group(
zone_fit.fitted_df,
group_cols=["sample", "region"],
target_col="Y",
high_count_threshold=zone_fit.high_count_threshold,
)
zone_oos_display_cols = [
"region",
"rows",
"mean_actual",
"mean_predicted",
"zero_share",
"mae",
"median_ae",
"p95_ae",
"rmse",
"count_r2",
"max_actual",
"max_predicted",
"high_outage_share",
"auc_high_outage",
"brier_high_outage",
]
display(html_table(pl.from_pandas(zone_oos_by_region.query("sample == 'test'")[zone_oos_display_cols].round(4))))
zone_coef = (
zone_fit.coefficients
.query("feature != 'const'")
.sort_values("abs_coefficient", ascending=False)
.copy()
)
coefficient_display_cols = [
"feature",
"coefficient",
"std_error",
"z_stat",
"p_value",
"incidence_rate_ratio",
"irr_low",
"irr_high",
]
display(html_table(pl.from_pandas(zone_coef[coefficient_display_cols].round(4))))
zone_significant = zone_coef.query("p_value < 0.05").copy()
display(html_table(pl.from_pandas(zone_significant[coefficient_display_cols].round(4))))
29 rows × 2 columns
| index | value |
|---|---|
| alpha | 0.00 |
| n_train | 4,341.00 |
| n_test | 1,086.00 |
| train_mean_actual | 27.84 |
| test_mean_actual | 22.92 |
| train_outage_share | 1.00 |
| test_outage_share | 1.00 |
| train_mae | 4.05 |
| test_mae | 4.61 |
| train_median_ae | 3.20 |
| test_median_ae | 3.83 |
| train_p95_ae | 10.34 |
| test_p95_ae | 11.26 |
| train_rmse | 5.28 |
| test_rmse | 5.89 |
| train_nll | 3.07 |
| test_nll | 3.19 |
| train_max_actual | 65.00 |
| test_max_actual | 63.00 |
| train_max_predicted | 58.60 |
| test_max_predicted | 45.92 |
| train_high_count_threshold | 38.00 |
| test_high_count_threshold | 38.00 |
| train_auc_outage | — |
| test_auc_outage | — |
| train_auc_high_outage | 0.88 |
| test_auc_high_outage | 0.82 |
| train_brier_high_outage | 0.06 |
| test_brier_high_outage | 0.03 |
3 rows × 15 columns
| region | rows | mean_actual | mean_predicted | zero_share | mae | median_ae | p95_ae | rmse | count_r2 | max_actual | max_predicted | high_outage_share | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NP15 | 362 | 23.05 | 22.91 | 0.000 | 3.724 | 3.006 | 9.22 | 4.761 | 0.123 | 43.00 | 33.82 | 0.008 | 0.538 | 0.009 |
| SP15 | 362 | 23.96 | 26.39 | 0.000 | 5.045 | 4.130 | 12.07 | 6.261 | 0.327 | 50.00 | 45.92 | 0.052 | 0.863 | 0.044 |
| ZP26 | 362 | 21.76 | 23.65 | 0.000 | 5.058 | 4.369 | 12.04 | 6.488 | 0.152 | 63.00 | 45.61 | 0.028 | 0.771 | 0.026 |
12 rows × 8 columns
| feature | coefficient | std_error | z_stat | p_value | incidence_rate_ratio | irr_low | irr_high |
|---|---|---|---|---|---|---|---|
| Y_lag1 | 0.111 | 0.004 | 31.24 | 0.000 | 1.117 | 1.109 | 1.125 |
| log_system_load_mean | 0.048 | 0.006 | 7.46 | 0.000 | 1.049 | 1.036 | 1.062 |
| weekend | −0.046 | 0.004 | −13.07 | 0.000 | 0.955 | 0.948 | 0.961 |
| month_sin | −0.035 | 0.005 | −6.60 | 0.000 | 0.966 | 0.956 | 0.976 |
| region_SP15 | 0.035 | 0.004 | 8.08 | 0.000 | 1.035 | 1.027 | 1.044 |
| month_cos | −0.031 | 0.005 | −5.86 | 0.000 | 0.970 | 0.960 | 0.980 |
| region_NP15 | −0.025 | 0.004 | −6.25 | 0.000 | 0.975 | 0.967 | 0.983 |
| active_forced_mw_mean_lag1 | −0.022 | 0.003 | −6.95 | 0.000 | 0.979 | 0.973 | 0.985 |
| planned_start_count_lag1 | −0.010 | 0.003 | −3.27 | 0.001 | 0.990 | 0.984 | 0.996 |
| cdd65 | −0.010 | 0.005 | −1.79 | 0.073 | 0.990 | 0.979 | 1.001 |
| hdd65 | 0.003 | 0.006 | 0.51 | 0.608 | 1.003 | 0.992 | 1.014 |
| temperature_f_std | 0.003 | 0.004 | 0.80 | 0.424 | 1.003 | 0.996 | 1.010 |
9 rows × 8 columns
| feature | coefficient | std_error | z_stat | p_value | incidence_rate_ratio | irr_low | irr_high |
|---|---|---|---|---|---|---|---|
| Y_lag1 | 0.111 | 0.004 | 31.24 | 0.000 | 1.117 | 1.109 | 1.125 |
| log_system_load_mean | 0.048 | 0.006 | 7.46 | 0.000 | 1.049 | 1.036 | 1.062 |
| weekend | −0.046 | 0.004 | −13.07 | 0.000 | 0.955 | 0.948 | 0.961 |
| month_sin | −0.035 | 0.005 | −6.60 | 0.000 | 0.966 | 0.956 | 0.976 |
| region_SP15 | 0.035 | 0.004 | 8.08 | 0.000 | 1.035 | 1.027 | 1.044 |
| month_cos | −0.031 | 0.005 | −5.86 | 0.000 | 0.970 | 0.960 | 0.980 |
| region_NP15 | −0.025 | 0.004 | −6.25 | 0.000 | 0.975 | 0.967 | 0.983 |
| active_forced_mw_mean_lag1 | −0.022 | 0.003 | −6.95 | 0.000 | 0.979 | 0.973 | 0.985 |
| planned_start_count_lag1 | −0.010 | 0.003 | −3.27 | 0.001 | 0.990 | 0.984 | 0.996 |
On the extended five-year panel the estimated alpha is 0.00 and the NB2-vs-Poisson likelihood-ratio test is 0.000: the conditional features absorb the overdispersion the raw counts show, and fitted variance equals fitted mean row for row. The probability layer is stated as P(high outage day) because every held-out zone-day has at least one matched forced start. Lagged outage starts, load, weekend effects, the seasonal terms and zone effects are the clearest conditional signals. None of the three weather terms survives the five-year panel: cdd65 at p 0.07, hdd65 at 0.61, and temperature dispersion - significant on the one-year panel - at 0.42.
coef_plot = zone_coef.iloc[::-1].copy()
coef_features = coef_plot["feature"].tolist()
fig = figure(
title="Market-Zone NB2 Coefficients, Standardized Features",
y_range=coef_features,
x_axis_label="log-count coefficient",
height=420,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.add_layout(Span(location=0, dimension="height", line_color="#0b0b0b", line_width=1))
fig.segment(
coef_plot["conf_low"], coef_features, coef_plot["conf_high"], coef_features,
line_color=SECONDARY, line_width=1.5,
)
fig.scatter(coef_plot["coefficient"], coef_features, size=8, color=SERIES[0])
show(style(fig))
7. Statistical Diagnostics#
The formal checks below compare the NB2 count specification against a Poisson benchmark, test feature blocks with likelihood-ratio chi-square tests, and compare NB-implied probabilities with a direct logistic probability model.
daily_feature_blocks = {
"outage state": ["Y_lag1", "planned_start_count_lag1", "active_forced_mw_mean_lag1"],
"weather": ["cdd65", "hdd65", "temperature_f_std"],
"load": ["log_system_load_mean"],
"calendar": ["weekend", "month_sin", "month_cos"],
"zone effects": ["region_NP15", "region_SP15"],
}
daily_nb_vs_poisson = poisson_negative_binomial_comparison(
market_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
maxiter=250,
)
display(html_table(pl.from_pandas(daily_nb_vs_poisson.round(4))))
daily_block_tests = likelihood_ratio_block_tests(
market_panel,
feature_cols=MARKET_ZONE_FEATURES,
feature_blocks=daily_feature_blocks,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
maxiter=250,
)
display(html_table(pl.from_pandas(daily_block_tests.round(4))))
daily_logit = fit_logistic_probability_model(
market_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
threshold=zone_fit.high_count_threshold,
test_size=0.20,
maxiter=200,
positive_label="high_outage_day",
)
daily_logit_oos = binary_model_performance_by_group(
daily_logit.fitted_df,
group_cols=["sample", "region"],
actual_col="binary_actual",
probability_col="binary_probability_hat",
)
daily_probability_model_compare = pd.DataFrame([
{
"model": "NB2-implied P(high outage day)",
"threshold": zone_fit.high_count_threshold,
"test_positive_share": zone_fit.metrics["test_high_outage_share"],
"test_auc": zone_fit.metrics["test_auc_high_outage"],
"test_brier": zone_fit.metrics["test_brier_high_outage"],
},
{
"model": "Logistic P(high outage day)",
"threshold": daily_logit.metrics["threshold"],
"test_positive_share": daily_logit.metrics["test_positive_share"],
"test_auc": daily_logit.metrics["test_auc"],
"test_brier": daily_logit.metrics["test_brier"],
},
])
display(html_table(pl.from_pandas(daily_probability_model_compare.round(4))))
display(
html_table(pl.from_pandas(
daily_logit_oos.query("sample == 'test'")[
["region", "rows", "positive_share", "mean_probability", "auc", "brier", "log_loss", "accuracy_at_50"]
].round(4)
))
)
display(
html_table(pl.from_pandas(
daily_logit.coefficients
.query("feature != 'const'")
.sort_values("abs_coefficient", ascending=False)
[["feature", "coefficient", "std_error", "z_stat", "p_value", "odds_ratio", "odds_ratio_low", "odds_ratio_high"]]
.head(12)
.round(4)
))
)
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 | −13,323.84 | 26,673.69 | 4.609 | 5.887 | 3.187 | — | — | — |
| Negative binomial NB2 | 0.000 | −13,323.84 | 26,675.68 | 4.609 | 5.887 | 3.187 | 0.002 | 1.000 | 0.482 |
5 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 | Y_lag1, planned_start_count_lag1, active_forced_mw_mean_lag1 | 3 | −13,323.84 | −13,790.77 | 933.85 | 0.000 | 3.187 | 3.628 | 0.441 |
| zone effects | region_NP15, region_SP15 | 2 | −13,323.84 | −13,427.77 | 207.85 | 0.000 | 3.187 | 3.147 | −0.040 |
| calendar | weekend, month_sin, month_cos | 3 | −13,323.84 | −13,417.47 | 187.26 | 0.000 | 3.187 | 3.224 | 0.037 |
| load | log_system_load_mean | 1 | −13,323.84 | −13,351.45 | 55.22 | 0.000 | 3.187 | 3.183 | −0.004 |
| weather | cdd65, hdd65, temperature_f_std | 3 | −13,323.84 | −13,326.12 | 4.55 | 0.208 | 3.187 | 3.183 | −0.004 |
2 rows × 5 columns
| model | threshold | test_positive_share | test_auc | test_brier |
|---|---|---|---|---|
| NB2-implied P(high outage day) | 38.00 | 0.029 | 0.818 | 0.026 |
| Logistic P(high outage day) | 38.00 | 0.029 | 0.810 | 0.026 |
3 rows × 8 columns
| region | rows | positive_share | mean_probability | auc | brier | log_loss | accuracy_at_50 |
|---|---|---|---|---|---|---|---|
| NP15 | 362 | 0.008 | 0.014 | 0.496 | 0.009 | 0.054 | 0.992 |
| SP15 | 362 | 0.052 | 0.104 | 0.862 | 0.046 | 0.171 | 0.945 |
| ZP26 | 362 | 0.028 | 0.024 | 0.769 | 0.025 | 0.113 | 0.975 |
12 rows × 8 columns
| feature | coefficient | std_error | z_stat | p_value | odds_ratio | odds_ratio_low | odds_ratio_high |
|---|---|---|---|---|---|---|---|
| Y_lag1 | 0.833 | 0.076 | 11.02 | 0.000 | 2.300 | 1.984 | 2.668 |
| log_system_load_mean | 0.772 | 0.134 | 5.76 | 0.000 | 2.164 | 1.664 | 2.814 |
| region_SP15 | 0.568 | 0.105 | 5.40 | 0.000 | 1.765 | 1.437 | 2.169 |
| region_NP15 | −0.313 | 0.125 | −2.51 | 0.012 | 0.732 | 0.573 | 0.934 |
| weekend | −0.256 | 0.079 | −3.23 | 0.001 | 0.774 | 0.663 | 0.904 |
| month_cos | −0.253 | 0.115 | −2.20 | 0.028 | 0.776 | 0.620 | 0.973 |
| cdd65 | −0.218 | 0.119 | −1.84 | 0.066 | 0.804 | 0.637 | 1.015 |
| active_forced_mw_mean_lag1 | −0.201 | 0.066 | −3.03 | 0.003 | 0.818 | 0.719 | 0.932 |
| planned_start_count_lag1 | −0.198 | 0.070 | −2.85 | 0.004 | 0.820 | 0.716 | 0.940 |
| month_sin | −0.064 | 0.127 | −0.50 | 0.617 | 0.938 | 0.731 | 1.204 |
| temperature_f_std | −0.062 | 0.080 | −0.77 | 0.439 | 0.940 | 0.804 | 1.099 |
| hdd65 | −0.008 | 0.159 | −0.05 | 0.962 | 0.992 | 0.727 | 1.355 |
8. Purged Expanding-Window Backtest#
Four expanding windows: each fold trains on history, drops the day immediately before the test window, and tests on the next unseen block, with a one-day embargo between scored windows.
zone_backtest = rolling_backtest_negative_binomial_model(
market_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
time_col="date",
n_splits=4,
initial_train_fraction=0.60,
purge_gap="1D",
embargo_gap="1D",
high_count_quantile=0.90,
maxiter=250,
)
backtest_fold_cols = [
"fold",
"train_start",
"train_end",
"test_start",
"test_end",
"n_train",
"n_test",
"alpha",
"mae",
"median_ae",
"p95_ae",
"rmse",
"nll",
"auc_high_outage",
"brier_high_outage",
]
display(html_table(pl.from_pandas(zone_backtest.fold_metrics[backtest_fold_cols].round(4))))
backtest_oos_cols = [
"rows",
"mean_actual",
"mean_predicted",
"mae",
"median_ae",
"p95_ae",
"rmse",
"nll",
"count_r2",
"high_outage_share",
"auc_high_outage",
"brier_high_outage",
]
display(html_table(pl.from_pandas(zone_backtest.overall_metrics[backtest_oos_cols].round(4))))
display(
html_table(pl.from_pandas(
zone_backtest.region_metrics[
[
"region",
"rows",
"mean_actual",
"mean_predicted",
"mae",
"median_ae",
"p95_ae",
"rmse",
"count_r2",
"auc_high_outage",
"brier_high_outage",
]
].round(4)
))
)
fig = figure(
title="Daily Backtest Count Error",
x_axis_label="fold",
y_axis_label="MAE",
height=320,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.line(zone_backtest.fold_metrics["fold"], zone_backtest.fold_metrics["mae"],
color=SERIES[0], line_width=2)
fig.scatter(zone_backtest.fold_metrics["fold"], zone_backtest.fold_metrics["mae"],
size=8, color=SERIES[0])
show(style(fig))
fig = figure(
title="Daily Backtest High-Outage Ranking",
x_axis_label="fold",
y_axis_label="AUC",
y_range=Range1d(0.45, 1.0),
height=320,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.add_layout(Span(location=0.50, dimension="width", line_color=SECONDARY,
line_width=1, line_dash="dotted"))
fig.line(zone_backtest.fold_metrics["fold"], zone_backtest.fold_metrics["auc_high_outage"],
color=SERIES[1], line_width=2)
fig.scatter(zone_backtest.fold_metrics["fold"], zone_backtest.fold_metrics["auc_high_outage"],
size=8, color=SERIES[1])
show(style(fig))
4 rows × 15 columns
| fold | train_start | train_end | test_start | test_end | n_train | n_test | alpha | mae | median_ae | p95_ae | rmse | nll | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2021-06-17 00:00 | 2024-06-04 00:00 | 2024-06-06 00:00 | 2024-12-03 00:00 | 3,252 | 543 | 0.000 | 4.360 | 3.605 | 10.81 | 5.718 | 3.117 | 0.842 | 0.099 |
| 2 | 2021-06-17 00:00 | 2024-12-03 00:00 | 2024-12-05 00:00 | 2025-07-04 00:00 | 3,798 | 543 | 0.000 | 5.640 | 5.062 | 12.81 | 6.849 | 3.472 | 0.813 | 0.023 |
| 3 | 2021-06-17 00:00 | 2025-07-04 00:00 | 2025-07-06 00:00 | 2026-01-02 00:00 | 4,344 | 543 | 0.000 | 4.836 | 4.065 | 11.48 | 6.133 | 3.238 | 0.874 | 0.035 |
| 4 | 2021-06-17 00:00 | 2026-01-02 00:00 | 2026-01-04 00:00 | 2026-06-30 00:00 | 4,890 | 534 | 0.002 | 4.231 | 3.388 | 10.95 | 5.515 | 3.103 | 0.671 | 0.018 |
1 rows × 12 columns
| rows | mean_actual | mean_predicted | mae | median_ae | p95_ae | rmse | nll | count_r2 | high_outage_share | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2,163 | 24.59 | 26.09 | 4.769 | 3.928 | 11.60 | 6.077 | 3.233 | 0.372 | 0.056 | 0.870 | 0.044 |
3 rows × 11 columns
| region | rows | mean_actual | mean_predicted | mae | median_ae | p95_ae | rmse | count_r2 | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|
| NP15 | 721 | 24.36 | 24.20 | 3.865 | 3.008 | 9.65 | 4.993 | 0.276 | 0.754 | 0.024 |
| SP15 | 721 | 26.31 | 28.78 | 5.350 | 4.593 | 12.19 | 6.660 | 0.447 | 0.894 | 0.072 |
| ZP26 | 721 | 23.10 | 25.28 | 5.092 | 4.196 | 12.05 | 6.443 | 0.267 | 0.808 | 0.035 |
Count error is stable across folds (MAE 4.2 to 5.6). High-outage ranking is not: AUC runs 0.84, 0.81, 0.87 and 0.67, so the pooled 0.87 rests on the three older folds and the most recent half-year is the weak one. The daily model holds as a count-intensity baseline, and notebook 12 scores this probability layer against an hourly alternative on the same folds.
9. Fitted Counts and High-Outage Probability#
P(Y > 0) is uninformative at this density, so the probability layer uses P(Y >= training p90) as the outage-stress indicator.
fitted = zone_fit.fitted_df.copy()
fitted["date_ts"] = pd.to_datetime(fitted["date"])
fitted_regions = sorted(fitted["region"].unique())
fig = figure(
title="Actual Counts and NB Expected Counts by Market Zone",
x_axis_type="datetime",
y_axis_label="matched forced starts",
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.add_layout(Span(location=zone_fit.high_count_threshold, dimension="width",
line_color=SECONDARY, line_width=1.2, line_dash="dotted"))
for region, color in zip(fitted_regions, SERIES):
zone_rows = fitted[fitted["region"] == region].sort_values("date_ts")
fig.line(zone_rows["date_ts"], zone_rows["Y"],
color=color, line_width=1.2, legend_label=region)
fig.line(zone_rows["date_ts"], zone_rows["mu_hat"],
color=color, line_width=1.4, line_dash="dashed")
fig.legend.location = "top_left"
show(style(fig))
fig = figure(
title="NB-Implied High-Outage Probability by Market Zone",
x_axis_type="datetime",
y_axis_label="P(zone-day count >= training p90)",
height=380,
sizing_mode="stretch_width",
toolbar_location="above",
)
high_days = fitted[fitted["high_outage_day_actual"].eq(1)]
for region, color in zip(fitted_regions, SERIES):
zone_rows = fitted[fitted["region"] == region].sort_values("date_ts")
fig.line(zone_rows["date_ts"], zone_rows["p_high_outage_hat"],
color=color, line_width=1.4, legend_label=region)
zone_high_days = high_days[high_days["region"] == region]
fig.scatter(zone_high_days["date_ts"], zone_high_days["p_high_outage_hat"],
size=5, color=color)
fig.legend.location = "top_left"
show(style(fig))
display(
html_table(pl.from_pandas(
fitted[["date", "region", "sample", "Y", "mu_hat", "var_hat", "p_high_outage_hat", "high_outage_day_actual"]]
.tail(12)
.round(4)
))
)
12 rows × 8 columns
| date | region | sample | Y | mu_hat | var_hat | p_high_outage_hat | high_outage_day_actual |
|---|---|---|---|---|---|---|---|
| 2026-06-27 | NP15 | test | 19 | 22.64 | 22.66 | 0.002 | 0 |
| 2026-06-27 | SP15 | test | 17 | 24.16 | 24.19 | 0.006 | 0 |
| 2026-06-27 | ZP26 | test | 23 | 22.98 | 23.00 | 0.003 | 0 |
| 2026-06-28 | NP15 | test | 17 | 20.45 | 20.46 | 0.000 | 0 |
| 2026-06-28 | SP15 | test | 19 | 22.67 | 22.69 | 0.002 | 0 |
| 2026-06-28 | ZP26 | test | 20 | 22.77 | 22.79 | 0.002 | 0 |
| 2026-06-29 | NP15 | test | 19 | 22.88 | 22.89 | 0.002 | 0 |
| 2026-06-29 | SP15 | test | 24 | 27.05 | 27.07 | 0.027 | 0 |
| 2026-06-29 | ZP26 | test | 17 | 24.97 | 24.99 | 0.009 | 0 |
| 2026-06-30 | NP15 | test | 18 | 23.22 | 23.24 | 0.003 | 0 |
| 2026-06-30 | SP15 | test | 23 | 28.91 | 28.94 | 0.060 | 0 |
| 2026-06-30 | ZP26 | test | 30 | 23.84 | 23.86 | 0.004 | 0 |
10. Weather-Sensitive Diagnostic Target#
The main target mixes outage causes, so the same NB2 setup is fitted to AMBIENT_DUE_TO_FUEL_INSUFFICIENCY starts as a weather-specific check.
weather_sensitive_features = [
"cdd65",
"hdd65",
"log_system_load_mean",
"fuel_insufficiency_start_count_lag1",
"active_forced_mw_mean_lag1",
"temperature_f_std",
"weekend",
"month_sin",
"month_cos",
"region_NP15",
"region_SP15",
]
fuel_fit = fit_negative_binomial_model(
market_panel,
feature_cols=weather_sensitive_features,
target_col=WEATHER_SENSITIVE_TARGET,
test_size=0.20,
high_count_quantile=0.90,
)
fuel_dispersion = (
market_panel.groupby("region")[WEATHER_SENSITIVE_TARGET]
.agg(
mean="mean",
variance="var",
zero_share=lambda s: (s == 0).mean(),
max="max",
)
)
fuel_dispersion["variance_to_mean"] = fuel_dispersion["variance"] / fuel_dispersion["mean"].replace(0, np.nan)
display(html_table(pl.from_pandas(fuel_dispersion.round(3).reset_index())))
display(
html_table(pl.from_pandas(
pd.Series(fuel_fit.metrics)
.loc[["alpha", "test_mae", "test_median_ae", "test_p95_ae", "test_rmse", "test_high_count_threshold", "test_auc_high_outage", "test_brier_high_outage"]]
.rename("value")
.to_frame()
.round(4)
.reset_index()
))
)
fuel_oos_by_region = model_performance_by_group(
fuel_fit.fitted_df,
group_cols=["sample", "region"],
target_col=WEATHER_SENSITIVE_TARGET,
high_count_threshold=fuel_fit.high_count_threshold,
)
display(
html_table(pl.from_pandas(
fuel_oos_by_region.query("sample == 'test'")[zone_oos_display_cols]
.round(4)
))
)
display(
html_table(pl.from_pandas(
fuel_fit.coefficients
.query("feature != 'const'")
.sort_values("abs_coefficient", ascending=False)
[["feature", "coefficient", "conf_low", "conf_high", "p_value"]]
.round(4)
))
)
3 rows × 6 columns
| region | mean | variance | zero_share | max | variance_to_mean |
|---|---|---|---|---|---|
| NP15 | 0.254 | 0.365 | 0.818 | 5 | 1.434 |
| SP15 | 0.241 | 0.539 | 0.860 | 8 | 2.233 |
| ZP26 | 0.193 | 0.378 | 0.856 | 12 | 1.957 |
8 rows × 2 columns
| index | value |
|---|---|
| alpha | 2.357 |
| test_mae | 0.299 |
| test_median_ae | 0.182 |
| test_p95_ae | 0.839 |
| test_rmse | 0.452 |
| test_high_count_threshold | 1.000 |
| test_auc_high_outage | 0.573 |
| test_brier_high_outage | 0.100 |
3 rows × 15 columns
| region | rows | mean_actual | mean_predicted | zero_share | mae | median_ae | p95_ae | rmse | count_r2 | max_actual | max_predicted | high_outage_share | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NP15 | 362 | 0.185 | 0.230 | 0.856 | 0.344 | 0.189 | 0.869 | 0.528 | −0.089 | 3.000 | 2.001 | 0.144 | 0.510 | 0.125 |
| SP15 | 362 | 0.121 | 0.221 | 0.906 | 0.283 | 0.191 | 0.817 | 0.422 | 0.004 | 4.000 | 3.694 | 0.094 | 0.596 | 0.085 |
| ZP26 | 362 | 0.119 | 0.193 | 0.903 | 0.269 | 0.162 | 0.851 | 0.396 | −0.051 | 2.000 | 1.084 | 0.097 | 0.612 | 0.089 |
11 rows × 5 columns
| feature | coefficient | conf_low | conf_high | p_value |
|---|---|---|---|---|
| fuel_insufficiency_start_count_lag1 | 0.527 | 0.454 | 0.601 | 0.000 |
| hdd65 | 0.126 | −0.024 | 0.276 | 0.099 |
| weekend | −0.118 | −0.215 | −0.022 | 0.016 |
| region_SP15 | 0.083 | −0.032 | 0.198 | 0.159 |
| month_cos | 0.081 | −0.059 | 0.222 | 0.256 |
| log_system_load_mean | 0.051 | −0.120 | 0.222 | 0.560 |
| region_NP15 | 0.035 | −0.070 | 0.141 | 0.509 |
| month_sin | 0.017 | −0.124 | 0.159 | 0.809 |
| temperature_f_std | 0.009 | −0.090 | 0.107 | 0.865 |
| cdd65 | −0.006 | −0.167 | 0.155 | 0.942 |
| active_forced_mw_mean_lag1 | 0.003 | −0.080 | 0.086 | 0.945 |
This target is sparse, so the estimates are less stable. It still helps separate a weather-aligned reason-code channel from the broader forced-outage process.
11. System-Level Benchmark#
The all-CAISO model is kept as a benchmark: full event coverage, but heterogeneous resources and zones collapse into one daily count.
system_panel, system_events = build_caiso_nb_panel(DATA_DIR)
system_fit = fit_negative_binomial_model(
system_panel,
feature_cols=DEFAULT_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
)
system_dispersion = overdispersion_summary(system_panel, target_col="Y")
zone_overall = overdispersion_summary(market_panel, target_col="Y")
comparison = pd.DataFrame([
{
"sampling_unit": "market zone × day, matched resources",
"rows": len(market_panel),
"mean_count": market_panel["Y"].mean(),
"variance_to_mean": zone_overall["variance_to_mean"],
"test_mae": zone_fit.metrics["test_mae"],
"test_rmse": zone_fit.metrics["test_rmse"],
"test_high_outage_auc": zone_fit.metrics["test_auc_high_outage"],
"coverage_note": "partial resource coverage, spatial signal retained",
},
{
"sampling_unit": "CAISO system × day, all resources",
"rows": len(system_panel),
"mean_count": system_panel["Y"].mean(),
"variance_to_mean": system_dispersion["variance_to_mean"],
"test_mae": system_fit.metrics["test_mae"],
"test_rmse": system_fit.metrics["test_rmse"],
"test_high_outage_auc": system_fit.metrics["test_auc_high_outage"],
"coverage_note": "full event coverage, spatial signal aggregated away",
},
])
display(html_table(pl.from_pandas(comparison.round(3))))
2 rows × 8 columns
| sampling_unit | rows | mean_count | variance_to_mean | test_mae | test_rmse | test_high_outage_auc | coverage_note |
|---|---|---|---|---|---|---|---|
| market zone × day, matched resources | 5,520 | 26.74 | 2.251 | 4.61 | 5.89 | 0.818 | partial resource coverage, spatial signal retained |
| CAISO system × day, all resources | 1,840 | 132.45 | 6.211 | 16.13 | 20.33 | 0.892 | full event coverage, spatial signal aggregated away |
12. Modeling Notes#
Primary panel: market zone by day for
NP15,SP15, andZP26.Outage-zone assignment: conservative matched-resource sample, with unmatched resources retained as a coverage diagnostic.
Weather-zone mismatch: kept as a basis-risk diagnostic;
NP15is the main zone where weather-zone proxies are mixed.Current probability object:
P(high outage day)from the NB2 conditional distribution.Robustness check: four-fold purged expanding-window backtest on future date blocks.
Statistical checks: coefficient z-statistics/p-values, NB2-vs-Poisson LR test, feature-block LR chi-square tests, and a logistic high-outage benchmark.
Main data need: a vetted CAISO resource-to-zone mapping table to raise coverage and reduce sampling error.