CAISO Outage Model Granularity Comparison#
Authors: Claire Kim · Brian Nguyen ·
Created: 2026-08-21 ·
Tags: #Experiment #Outages #Validation
Daily, hourly, and 15-minute sampling of the CAISO outage model are compared. The daily market-zone NB2 is the baseline; the hourly version holds the mapping and market zones fixed and asks whether outage probability improves when timing is modeled directly.
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 ColumnDataSource, FactorRange, HoverTool, Span
from bokeh.plotting import figure
from bokeh.resources import INLINE
from statsmodels.tools.sm_exceptions import (
ConvergenceWarning,
HessianInversionWarning,
)
from bokeh.transform import dodge
from sklearn.metrics import roc_auc_score, r2_score
pd.set_option("display.max_columns", 160)
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 (
MARKET_ZONE_FEATURES,
MARKET_ZONE_HOURLY_FEATURES,
build_caiso_market_zone_nb_panel,
build_caiso_market_zone_hourly_nb_panel,
fit_negative_binomial_model,
binary_model_performance_by_group,
likelihood_ratio_block_tests,
poisson_negative_binomial_comparison,
rolling_backtest_logistic_probability_model,
fit_logistic_probability_model,
rolling_backtest_negative_binomial_model,
model_performance_by_group,
interval_market_zone_outage_counts,
interval_market_zone_target_summary,
resource_mapping_coverage_summary,
)
from chart_theme import 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. Baseline Mapping Layer#
The same conservative resource-to-plant-to-market-zone crosswalk is used across granularities, so the comparison isolates time sampling.
daily_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=daily_panel["date"].min(),
end_date=daily_panel["date"].max(),
)
coverage_display = coverage[
[
"outage_type",
"events",
"matched_events",
"event_coverage",
"unique_resources",
"matched_unique_resources",
"resource_coverage",
"crosswalk_rows",
]
]
display(html_table(pl.from_pandas(coverage_display.round(3))))
print(f"Daily baseline rows: {len(daily_panel):,}")
print(f"Date window: {daily_panel['date'].min()} to {daily_panel['date'].max()}")
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 |
Daily baseline rows: 5,520
Date window: 2021-06-17 to 2026-06-30
Coverage is unchanged from the daily model. Unmatched resources remain sampling error and basis risk.
2. Target Sparsity Across Granularities#
Daily counts are dense, hourly counts sparse but workable, and 15-minute counts sparse enough to push the problem toward event ranking or a hurdle model rather than count prediction.
hourly_panel, _, _ = build_caiso_market_zone_hourly_nb_panel(DATA_DIR)
summary_rows = []
for label, panel in [("daily", daily_panel), ("hourly", hourly_panel)]:
grouped = panel.groupby("region")["Y"].agg(
rows="size",
mean="mean",
variance="var",
zero_share=lambda s: (s == 0).mean(),
outage_share=lambda s: (s > 0).mean(),
p90=lambda s: s.quantile(0.90),
p95=lambda s: s.quantile(0.95),
p99=lambda s: s.quantile(0.99),
max="max",
).reset_index()
grouped["frequency"] = label
grouped["variance_to_mean"] = grouped["variance"] / grouped["mean"].replace(0, np.nan)
summary_rows.append(grouped)
fifteen_summary = interval_market_zone_target_summary(
mapped_events,
start_date=daily_panel["date"].min(),
end_date=daily_panel["date"].max(),
freq="15min",
)
fifteen_summary["frequency"] = "15min"
granularity_summary = pd.concat(summary_rows + [fifteen_summary], ignore_index=True)
cols = ["frequency", "region", "rows", "mean", "variance", "variance_to_mean", "zero_share", "outage_share", "p90", "p95", "p99", "max"]
display(html_table(pl.from_pandas(granularity_summary[cols].round(4))))
hourly_start_hour_summary = (
hourly_panel.groupby("hour")["Y"]
.agg(
region_hour_rows="size",
forced_starts="sum",
mean_starts="mean",
zero_share=lambda s: (s == 0).mean(),
max_starts="max",
)
.reset_index()
)
hourly_start_hour_summary["start_share"] = (
hourly_start_hour_summary["forced_starts"] / hourly_start_hour_summary["forced_starts"].sum()
)
display(html_table(pl.from_pandas(hourly_start_hour_summary.round(4))))
9 rows × 12 columns
| frequency | region | rows | mean | variance | variance_to_mean | zero_share | outage_share | p90 | p95 | p99 | max |
|---|---|---|---|---|---|---|---|---|---|---|---|
| daily | NP15 | 1,840 | 24.62 | 33.94 | 1.38 | 0.000 | 1.000 | 32.00 | 35.00 | 41.00 | 52 |
| daily | SP15 | 1,840 | 29.75 | 84.13 | 2.83 | 0.000 | 1.000 | 42.00 | 46.00 | 52.00 | 65 |
| daily | ZP26 | 1,840 | 25.85 | 48.24 | 1.87 | 0.000 | 1.000 | 34.00 | 37.00 | 42.00 | 63 |
| hourly | NP15 | 44,160 | 1.03 | 10.49 | 10.22 | 0.680 | 0.320 | 2.00 | 3.00 | 18.00 | 35 |
| hourly | SP15 | 44,160 | 1.24 | 8.59 | 6.93 | 0.567 | 0.433 | 3.00 | 5.00 | 16.00 | 38 |
| hourly | ZP26 | 44,160 | 1.08 | 9.74 | 9.04 | 0.633 | 0.367 | 2.00 | 4.00 | 18.00 | 33 |
| 15min | NP15 | 176,640 | 0.26 | 2.75 | 10.72 | 0.908 | 0.092 | 0.00 | 1.00 | 8.00 | 34 |
| 15min | SP15 | 176,640 | 0.31 | 2.31 | 7.45 | 0.868 | 0.132 | 1.00 | 1.00 | 7.00 | 38 |
| 15min | ZP26 | 176,640 | 0.27 | 2.57 | 9.55 | 0.892 | 0.108 | 1.00 | 1.00 | 6.00 | 31 |
24 rows × 7 columns
| hour | region_hour_rows | forced_starts | mean_starts | zero_share | max_starts | start_share |
|---|---|---|---|---|---|---|
| 0 | 5,520 | 79,445 | 14.39 | 0.001 | 38 | 0.538 |
| 1 | 5,535 | 5,071 | 0.92 | 0.699 | 10 | 0.034 |
| 2 | 5,505 | 844 | 0.15 | 0.872 | 8 | 0.006 |
| 3 | 5,520 | 959 | 0.17 | 0.858 | 7 | 0.006 |
| 4 | 5,520 | 1,239 | 0.22 | 0.820 | 9 | 0.008 |
| 5 | 5,520 | 2,354 | 0.43 | 0.700 | 14 | 0.016 |
| 6 | 5,520 | 5,423 | 0.98 | 0.449 | 10 | 0.037 |
| 7 | 5,520 | 4,801 | 0.87 | 0.483 | 9 | 0.033 |
| 8 | 5,520 | 4,874 | 0.88 | 0.461 | 10 | 0.033 |
| 9 | 5,520 | 4,634 | 0.84 | 0.488 | 10 | 0.031 |
| 10 | 5,520 | 4,463 | 0.81 | 0.498 | 13 | 0.030 |
| 11 | 5,520 | 3,897 | 0.71 | 0.541 | 9 | 0.026 |
| 12 | 5,520 | 3,797 | 0.69 | 0.539 | 17 | 0.026 |
| 13 | 5,520 | 2,902 | 0.53 | 0.622 | 8 | 0.020 |
| 14 | 5,520 | 2,682 | 0.49 | 0.646 | 14 | 0.018 |
| 15 | 5,520 | 2,474 | 0.45 | 0.672 | 8 | 0.017 |
| 16 | 5,520 | 2,637 | 0.48 | 0.662 | 10 | 0.018 |
| 17 | 5,520 | 2,737 | 0.50 | 0.651 | 9 | 0.018 |
| 18 | 5,520 | 2,906 | 0.53 | 0.641 | 9 | 0.020 |
| 19 | 5,520 | 2,253 | 0.41 | 0.707 | 7 | 0.015 |
| 20 | 5,520 | 2,213 | 0.40 | 0.710 | 11 | 0.015 |
| 21 | 5,520 | 1,768 | 0.32 | 0.759 | 7 | 0.012 |
| 22 | 5,520 | 1,595 | 0.29 | 0.787 | 8 | 0.011 |
| 23 | 5,520 | 1,648 | 0.30 | 0.772 | 10 | 0.011 |
Zero share is the transition: a zone almost always has at least one matched forced outage start per day, while roughly 57-68% of region-hours and 87-91% of region-intervals at 15-minute granularity are zero.
plot_data = granularity_summary.copy()
plot_data["frequency"] = pd.Categorical(plot_data["frequency"], ["daily", "hourly", "15min"], ordered=True)
plot_data = plot_data.sort_values(["frequency", "region"])
zero_share_factors = list(zip(plot_data["frequency"].astype(str), plot_data["region"]))
zero_share_src = ColumnDataSource(
{"x": zero_share_factors, "zero_share": plot_data["zero_share"].to_list()}
)
fig = figure(
title="Zero-Count Share by Granularity and Market Zone",
x_range=FactorRange(*zero_share_factors),
y_range=(0, 1),
y_axis_label="Zero Share",
height=360,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.vbar(x="x", top="zero_share", width=0.8, color=SERIES[0], source=zero_share_src)
fig.add_tools(
HoverTool(tooltips=[("granularity, zone", "@x"), ("zero share", "@zero_share{0.000}")])
)
show(style(fig))
start_hour_src = ColumnDataSource(
{
"hour": hourly_start_hour_summary["hour"].to_list(),
"start_share": hourly_start_hour_summary["start_share"].to_list(),
}
)
fig = figure(
title="Matched Forced Starts by Pacific Start Hour",
x_axis_label="Pacific Hour",
y_axis_label="Share of Starts",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.vbar(x="hour", top="start_share", width=0.8, color=SERIES[0], source=start_hour_src)
fig.add_tools(
HoverTool(tooltips=[("hour", "@hour"), ("share of starts", "@start_share{0.000}")])
)
show(style(fig))
A large share of recorded starts lands at midnight Pacific, so the hourly model is a probability-ranking prototype pending a timestamp-quality review.
3. Daily NB2 Baseline#
The daily model estimates zone-day outage intensity.
daily_fit = fit_negative_binomial_model(
daily_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
)
daily_metrics = pd.Series(daily_fit.metrics).loc[
[
"alpha",
"n_train",
"n_test",
"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_auc_high_outage",
"test_auc_high_outage",
"train_brier_high_outage",
"test_brier_high_outage",
"train_high_count_threshold",
"test_high_count_threshold",
]
].rename("daily_value").to_frame()
display(html_table(pl.from_pandas(daily_metrics.round(4).reset_index())))
daily_oos_by_region = model_performance_by_group(
daily_fit.fitted_df,
group_cols=["sample", "region"],
target_col="Y",
high_count_threshold=daily_fit.high_count_threshold,
)
oos_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(daily_oos_by_region.query("sample == 'test'")[oos_cols].round(4))))
display(
html_table(
pl.from_pandas(
daily_fit.coefficients
.query("feature != 'const'")
.sort_values("abs_coefficient", ascending=False)
[["feature", "coefficient", "std_error", "z_stat", "p_value", "incidence_rate_ratio", "irr_low", "irr_high"]]
.head(10)
.round(4)
)
)
)
19 rows × 2 columns
| index | daily_value |
|---|---|
| alpha | 0.00 |
| n_train | 4,341.00 |
| n_test | 1,086.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_auc_high_outage | 0.88 |
| test_auc_high_outage | 0.82 |
| train_brier_high_outage | 0.06 |
| test_brier_high_outage | 0.03 |
| train_high_count_threshold | 38.00 |
| test_high_count_threshold | 38.00 |
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 |
10 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 |
Daily tracks the level (test MAE 4.61 on a mean of 27) and ranks high-count days well (test AUC 0.82). What it cannot express is P(Y > 0): every held-out region-day has a positive matched count, so the daily probability readout is P(high outage day) and nothing else. A high outage period is one whose count reaches the 90th percentile of the training counts, rounded up to a whole number: high_count_threshold = ceil(quantile(Y_train, 0.90)), reported per fit as train_high_count_threshold / test_high_count_threshold. The threshold is fixed on training data, so the held-out share sitting above it is a result rather than a definition, and it is a count of outage starts - not megawatts and not a severity measure. Sections 4 to 7 ask whether hourly sampling buys a sharper probability layer than this one.
4. Hourly NB2 Extension#
The target is new forced outage starts per market-zone hour, with hourly load, active outage MW, and Open-Meteo weather joined at the same timestamp. Lagged outage and active-MW terms are log-scaled so the NB log link is not dominated by a few spike hours.
hourly_fit = fit_negative_binomial_model(
hourly_panel,
feature_cols=MARKET_ZONE_HOURLY_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
maxiter=250,
)
hourly_metrics = pd.Series(hourly_fit.metrics).loc[
[
"alpha",
"n_train",
"n_test",
"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_auc_outage",
"test_auc_outage",
"train_auc_high_outage",
"test_auc_high_outage",
"train_brier_outage",
"test_brier_outage",
"train_brier_high_outage",
"test_brier_high_outage",
"train_high_count_threshold",
"test_high_count_threshold",
]
].rename("hourly_value").to_frame()
display(html_table(pl.from_pandas(hourly_metrics.round(4).reset_index())))
hourly_oos_by_region = model_performance_by_group(
hourly_fit.fitted_df,
group_cols=["sample", "region"],
target_col="Y",
high_count_threshold=hourly_fit.high_count_threshold,
)
display(html_table(pl.from_pandas(hourly_oos_by_region.query("sample == 'test'")[oos_cols].round(4))))
display(
html_table(
pl.from_pandas(
hourly_fit.coefficients
.query("feature != 'const'")
.sort_values("abs_coefficient", ascending=False)
[["feature", "coefficient", "std_error", "z_stat", "p_value", "incidence_rate_ratio", "irr_low", "irr_high"]]
.head(12)
.round(4)
)
)
)
23 rows × 2 columns
| index | hourly_value |
|---|---|
| alpha | 0.45 |
| n_train | 104,100.00 |
| n_test | 26,028.00 |
| train_mae | 0.82 |
| test_mae | 0.80 |
| train_median_ae | 0.40 |
| test_median_ae | 0.38 |
| train_p95_ae | 2.65 |
| test_p95_ae | 2.65 |
| train_rmse | 1.44 |
| test_rmse | 1.40 |
| train_nll | 1.08 |
| test_nll | 1.07 |
| train_auc_outage | 0.69 |
| test_auc_outage | 0.68 |
| train_auc_high_outage | 0.78 |
| test_auc_high_outage | 0.78 |
| train_brier_outage | 0.21 |
| test_brier_outage | 0.21 |
| train_brier_high_outage | 0.10 |
| test_brier_high_outage | 0.10 |
| train_high_count_threshold | 2.00 |
| test_high_count_threshold | 2.00 |
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 | 8,676 | 0.961 | 0.789 | 0.667 | 0.766 | 0.331 | 2.667 | 1.438 | 0.737 | 27.00 | 21.25 | 0.117 | 0.771 | 0.075 |
| SP15 | 8,676 | 0.999 | 0.975 | 0.594 | 0.828 | 0.438 | 2.622 | 1.374 | 0.608 | 27.00 | 24.43 | 0.180 | 0.797 | 0.111 |
| ZP26 | 8,676 | 0.907 | 0.792 | 0.622 | 0.804 | 0.384 | 2.672 | 1.378 | 0.580 | 24.00 | 18.78 | 0.149 | 0.754 | 0.099 |
12 rows × 8 columns
| feature | coefficient | std_error | z_stat | p_value | incidence_rate_ratio | irr_low | irr_high |
|---|---|---|---|---|---|---|---|
| log1p_Y_lag24h | 0.871 | 0.003 | 257.39 | 0.000 | 2.388 | 2.373 | 2.404 |
| hour_sin | 0.101 | 0.006 | 17.76 | 0.000 | 1.107 | 1.094 | 1.119 |
| log_system_load_mean | 0.101 | 0.007 | 14.83 | 0.000 | 1.106 | 1.091 | 1.121 |
| weekend | −0.095 | 0.005 | −19.71 | 0.000 | 0.909 | 0.900 | 0.918 |
| region_SP15 | 0.063 | 0.005 | 11.62 | 0.000 | 1.065 | 1.054 | 1.076 |
| log1p_active_forced_mw_lag1h | −0.053 | 0.004 | −14.91 | 0.000 | 0.948 | 0.942 | 0.955 |
| month_sin | 0.044 | 0.007 | 6.68 | 0.000 | 1.045 | 1.032 | 1.059 |
| region_NP15 | −0.039 | 0.006 | −6.88 | 0.000 | 0.962 | 0.951 | 0.973 |
| hour_cos | −0.029 | 0.006 | −4.64 | 0.000 | 0.971 | 0.959 | 0.983 |
| log1p_Y_lag1h | 0.027 | 0.004 | 6.08 | 0.000 | 1.027 | 1.018 | 1.036 |
| log1p_planned_start_count_lag1h | 0.024 | 0.004 | 5.53 | 0.000 | 1.024 | 1.015 | 1.032 |
| temperature_f_std | −0.012 | 0.005 | −2.43 | 0.015 | 0.988 | 0.979 | 0.998 |
Both probability objects are meaningful at hourly granularity: P(Y > 0) for any outage-start hour and P(Y >= 2) for a higher-count hour. The strongest signals are the 24-hour lag, load, and hour-of-day structure.
5. Specification Tests#
The count models are checked against Poisson benchmarks and nested feature-block likelihood-ratio tests.
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"],
}
hourly_feature_blocks = {
"outage state": ["log1p_Y_lag1h", "log1p_Y_lag24h", "log1p_planned_start_count_lag1h", "log1p_active_forced_mw_lag1h"],
"weather": ["cdh65", "hdh65", "temperature_f_std"],
"load": ["log_system_load_mean"],
"calendar": ["hour_sin", "hour_cos", "weekend", "month_sin", "month_cos"],
"zone effects": ["region_NP15", "region_SP15"],
}
daily_nb_vs_poisson = poisson_negative_binomial_comparison(
daily_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
maxiter=250,
)
hourly_nb_vs_poisson = poisson_negative_binomial_comparison(
hourly_panel,
feature_cols=MARKET_ZONE_HOURLY_FEATURES,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
maxiter=250,
)
nb_vs_poisson = pd.concat(
[
daily_nb_vs_poisson.assign(frequency="daily"),
hourly_nb_vs_poisson.assign(frequency="hourly"),
],
ignore_index=True,
)
display(html_table(pl.from_pandas(nb_vs_poisson[["frequency", "model", "alpha", "train_llf", "aic", "test_mae", "test_rmse", "test_nll", "lr_chi2_vs_poisson", "lr_p_value_boundary_mixture"]].round(4))))
daily_block_tests = likelihood_ratio_block_tests(
daily_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,
).assign(frequency="daily")
hourly_block_tests = likelihood_ratio_block_tests(
hourly_panel,
feature_cols=MARKET_ZONE_HOURLY_FEATURES,
feature_blocks=hourly_feature_blocks,
target_col="Y",
test_size=0.20,
high_count_quantile=0.90,
maxiter=250,
).assign(frequency="hourly")
block_tests = pd.concat([daily_block_tests, hourly_block_tests], ignore_index=True)
display(
html_table(
pl.from_pandas(
block_tests[
["frequency", "block", "df", "lr_chi2", "p_value", "full_test_nll", "reduced_test_nll", "delta_test_nll_reduced_minus_full"]
].round(4)
)
)
)
plot_blocks = block_tests.copy()
plot_blocks["minus_log10_p"] = -np.log10(plot_blocks["p_value"].clip(lower=1e-300))
block_order = list(dict.fromkeys(plot_blocks["block"]))
block_wide = (
plot_blocks.pivot(index="block", columns="frequency", values="minus_log10_p")
.reindex(block_order)
.reset_index()
)
block_src = ColumnDataSource(
{
"block": block_wide["block"].to_list(),
"daily": block_wide["daily"].to_list(),
"hourly": block_wide["hourly"].to_list(),
}
)
fig = figure(
title="NB2 Feature-Block Likelihood-Ratio Tests",
x_range=FactorRange(*block_order),
y_axis_label="-log10 p-value",
height=360,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.vbar(x=dodge("block", -0.17, range=fig.x_range), top="daily", width=0.3,
color=SERIES[0], legend_label="daily", source=block_src)
fig.vbar(x=dodge("block", 0.17, range=fig.x_range), top="hourly", width=0.3,
color=SERIES[1], legend_label="hourly", source=block_src)
fig.add_layout(Span(location=-np.log10(0.05), dimension="width",
line_color="#555555", line_width=1, line_dash="dotted"))
fig.add_tools(
HoverTool(tooltips=[("block", "@block"), ("daily", "@daily{0.00}"), ("hourly", "@hourly{0.00}")])
)
fig.legend.location = "top_right"
show(style(fig))
4 rows × 10 columns
| frequency | model | alpha | train_llf | aic | test_mae | test_rmse | test_nll | lr_chi2_vs_poisson | lr_p_value_boundary_mixture |
|---|---|---|---|---|---|---|---|---|---|
| daily | Poisson | 0.000 | −13,323.84 | 26,673.69 | 4.609 | 5.887 | 3.187 | — | — |
| daily | Negative binomial NB2 | 0.000 | −13,323.84 | 26,675.68 | 4.609 | 5.887 | 3.187 | 0.00 | 0.482 |
| hourly | Poisson | 0.000 | −113,745.99 | 227,523.99 | 0.756 | 1.261 | 1.098 | — | — |
| hourly | Negative binomial NB2 | 0.453 | −112,148.28 | 224,330.55 | 0.799 | 1.397 | 1.066 | 3,195.44 | 0.000 |
10 rows × 8 columns
| frequency | block | df | lr_chi2 | p_value | full_test_nll | reduced_test_nll | delta_test_nll_reduced_minus_full |
|---|---|---|---|---|---|---|---|
| daily | outage state | 3 | 933.85 | 0.000 | 3.187 | 3.628 | 0.441 |
| daily | zone effects | 2 | 207.85 | 0.000 | 3.187 | 3.147 | −0.040 |
| daily | calendar | 3 | 187.26 | 0.000 | 3.187 | 3.224 | 0.037 |
| daily | load | 1 | 55.22 | 0.000 | 3.187 | 3.183 | −0.004 |
| daily | weather | 3 | 4.55 | 0.208 | 3.187 | 3.183 | −0.004 |
| hourly | outage state | 4 | 41,972.55 | 0.000 | 1.066 | 1.236 | 0.169 |
| hourly | calendar | 5 | 935.33 | 0.000 | 1.066 | 1.072 | 0.006 |
| hourly | zone effects | 2 | 372.05 | 0.000 | 1.066 | 1.067 | 0.000 |
| hourly | load | 1 | 221.92 | 0.000 | 1.066 | 1.068 | 0.002 |
| hourly | weather | 3 | 7.28 | 0.063 | 1.066 | 1.066 | 0.000 |
The daily fit puts alpha at 0.000 and the likelihood-ratio test against Poisson at 0.00, so the daily NB2 is a Poisson in everything but name - the conditional features absorb the overdispersion the raw zone-day counts show. Only the hourly target keeps genuine residual dispersion, at alpha 0.45. The block tests separate in-sample significance from holdout usefulness: a block can have a small p-value while adding only a small improvement in test NLL.
6. Daily vs Hourly Performance View#
Daily estimates zone-day outage intensity; hourly estimates outage-start timing and short-window event risk.
comparison = pd.DataFrame([
{
"model": "daily market-zone NB2",
"rows": len(daily_fit.fitted_df),
"target_mean": daily_fit.fitted_df["Y"].mean(),
"zero_share": (daily_fit.fitted_df["Y"] == 0).mean(),
"alpha": daily_fit.metrics["alpha"],
"test_mae": daily_fit.metrics["test_mae"],
"test_median_ae": daily_fit.metrics["test_median_ae"],
"test_p95_ae": daily_fit.metrics["test_p95_ae"],
"test_rmse": daily_fit.metrics["test_rmse"],
"test_max_predicted": daily_fit.metrics["test_max_predicted"],
"test_nll": daily_fit.metrics["test_nll"],
"test_auc_outage": daily_fit.metrics["test_auc_outage"],
"test_auc_high_outage": daily_fit.metrics["test_auc_high_outage"],
"test_brier_high_outage": daily_fit.metrics["test_brier_high_outage"],
"probability_use": "P(high outage day)",
},
{
"model": "hourly market-zone NB2",
"rows": len(hourly_fit.fitted_df),
"target_mean": hourly_fit.fitted_df["Y"].mean(),
"zero_share": (hourly_fit.fitted_df["Y"] == 0).mean(),
"alpha": hourly_fit.metrics["alpha"],
"test_mae": hourly_fit.metrics["test_mae"],
"test_median_ae": hourly_fit.metrics["test_median_ae"],
"test_p95_ae": hourly_fit.metrics["test_p95_ae"],
"test_rmse": hourly_fit.metrics["test_rmse"],
"test_max_predicted": hourly_fit.metrics["test_max_predicted"],
"test_nll": hourly_fit.metrics["test_nll"],
"test_auc_outage": hourly_fit.metrics["test_auc_outage"],
"test_auc_high_outage": hourly_fit.metrics["test_auc_high_outage"],
"test_brier_high_outage": hourly_fit.metrics["test_brier_high_outage"],
"probability_use": "P(any outage hour), P(high outage hour)",
},
])
display(html_table(pl.from_pandas(comparison.round(4))))
auc_plot = pd.concat(
[
daily_oos_by_region.query("sample == 'test'").assign(model="daily: high-outage day", auc=lambda df: df["auc_high_outage"]),
hourly_oos_by_region.query("sample == 'test'").assign(model="hourly: any-outage hour", auc=lambda df: df["auc_outage"]),
hourly_oos_by_region.query("sample == 'test'").assign(model="hourly: high-outage hour", auc=lambda df: df["auc_high_outage"]),
],
ignore_index=True,
)
auc_plot = auc_plot[["model", "region", "auc"]].dropna()
auc_models = list(dict.fromkeys(auc_plot["model"]))
auc_regions = sorted(auc_plot["region"].unique())
auc_wide = (
auc_plot.pivot(index="region", columns="model", values="auc")
.reindex(auc_regions)
.reset_index()
)
auc_src = ColumnDataSource(
{
"region": auc_wide["region"].to_list(),
**{model: auc_wide[model].to_list() for model in auc_models},
}
)
fig = figure(
title="Out-of-Sample Probability Ranking by Market Zone",
x_range=FactorRange(*auc_regions),
y_range=(0.45, 1.0),
y_axis_label="Test AUC",
height=360,
sizing_mode="stretch_width",
toolbar_location="above",
)
for offset, color, model in zip([-0.25, 0.0, 0.25], SERIES, auc_models):
fig.vbar(x=dodge("region", offset, range=fig.x_range), top=model, width=0.22,
color=color, legend_label=model, source=auc_src)
fig.add_layout(Span(location=0.50, dimension="width", line_color="#555555",
line_width=1, line_dash="dotted"))
fig.legend.location = "bottom_right"
show(style(fig))
2 rows × 15 columns
| model | rows | target_mean | zero_share | alpha | test_mae | test_median_ae | test_p95_ae | test_rmse | test_max_predicted | test_nll | test_auc_outage | test_auc_high_outage | test_brier_high_outage | probability_use |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| daily market-zone NB2 | 5,427 | 26.85 | 0.000 | 0.000 | 4.609 | 3.830 | 11.26 | 5.887 | 45.92 | 3.187 | — | 0.818 | 0.026 | P(high outage day) |
| hourly market-zone NB2 | 130,128 | 1.12 | 0.626 | 0.453 | 0.799 | 0.384 | 2.65 | 1.397 | 24.43 | 1.066 | 0.680 | 0.777 | 0.095 | P(any outage hour), P(high outage hour) |
Three of the four granularity-event combinations are rankable. Daily zone-days are never zero (zero_share 0.000, every zone starts at least one outage each day), so P(any outage day) has a single class and no AUC is defined - the dash in the table above. Hourly is 62.6% zeros, which is what makes the any-outage event rankable at all, and P(any outage hour) is the one question only hourly can answer. Read each AUC by its own event: daily ranks high-count days, hourly ranks any-outage hours and high-outage hours. The two high-outage columns are the only near-comparison, and even there the events differ - a day above 38 starts against an hour above 2 - so section 7 scores them on repeated folds rather than on this single split.
hourly_fitted = hourly_fit.fitted_df.copy()
hourly_fitted["interval_start_utc"] = pd.to_datetime(hourly_fitted["interval_start_utc"])
recent_start = hourly_fitted["interval_start_utc"].max() - pd.Timedelta(21, "D")
recent = hourly_fitted[hourly_fitted["interval_start_utc"] >= recent_start]
recent_regions = sorted(recent["region"].unique())
fig = figure(
title="Hourly Actual Counts and NB Expected Counts, Last 21 Days",
x_axis_type="datetime",
y_axis_label="Forced Starts/Hour",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for color, region in zip(SERIES, recent_regions):
sub = recent[recent["region"].eq(region)].sort_values("interval_start_utc")
fig.line(sub["interval_start_utc"], sub["Y"], color=color, line_width=1.0,
legend_label=f"{region} actual")
fig.line(sub["interval_start_utc"], sub["mu_hat"], color=color, line_width=1.2,
line_dash="dashed", legend_label=f"{region} expected")
fig.legend.location = "top_left"
show(style(fig))
fig = figure(
title="Hourly NB-Implied P(any outage start)",
x_axis_type="datetime",
y_axis_label="P(Y > 0)",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for color, region in zip(SERIES, recent_regions):
sub = recent[recent["region"].eq(region)].sort_values("interval_start_utc")
fig.line(sub["interval_start_utc"], sub["p_outage_hat"], color=color,
line_width=1.2, legend_label=region)
fig.legend.location = "top_left"
show(style(fig))
The hourly fitted line is a conditional intensity, not a spike reconstruction: it ranks outage-prone hours rather than matching counts one-for-one.
7. Purged Walk-Forward Backtests#
The robustness check uses four future test windows, each trained on earlier observations with the immediate lookback window purged before the test block and an embargo gap between scored blocks.
daily_backtest = rolling_backtest_negative_binomial_model(
daily_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,
)
hourly_backtest = rolling_backtest_negative_binomial_model(
hourly_panel,
feature_cols=MARKET_ZONE_HOURLY_FEATURES,
target_col="Y",
time_col="interval_start_utc",
n_splits=4,
initial_train_fraction=0.60,
purge_gap="24h",
embargo_gap="24h",
high_count_quantile=0.90,
maxiter=250,
)
daily_logit_backtest = rolling_backtest_logistic_probability_model(
daily_panel,
feature_cols=MARKET_ZONE_FEATURES,
target_col="Y",
threshold=daily_fit.high_count_threshold,
time_col="date",
n_splits=4,
initial_train_fraction=0.60,
purge_gap="1D",
embargo_gap="1D",
maxiter=200,
positive_label="high_outage_day",
)
hourly_logit_backtest = rolling_backtest_logistic_probability_model(
hourly_panel,
feature_cols=MARKET_ZONE_HOURLY_FEATURES,
target_col="Y",
threshold=hourly_fit.high_count_threshold,
time_col="interval_start_utc",
n_splits=4,
initial_train_fraction=0.60,
purge_gap="24h",
embargo_gap="24h",
maxiter=200,
positive_label="high_outage_hour",
)
backtest_comparison = pd.DataFrame([
{
"model": "daily NB2: P(high outage day)",
"folds": len(daily_backtest.fold_metrics),
"rows": daily_backtest.overall_metrics.loc[0, "rows"],
"mae": daily_backtest.overall_metrics.loc[0, "mae"],
"median_ae": daily_backtest.overall_metrics.loc[0, "median_ae"],
"p95_ae": daily_backtest.overall_metrics.loc[0, "p95_ae"],
"rmse": daily_backtest.overall_metrics.loc[0, "rmse"],
"nll": daily_backtest.overall_metrics.loc[0, "nll"],
"auc_outage": daily_backtest.overall_metrics.loc[0, "auc_outage"],
"auc_high_outage": daily_backtest.overall_metrics.loc[0, "auc_high_outage"],
"brier_high_outage": daily_backtest.overall_metrics.loc[0, "brier_high_outage"],
},
{
"model": "hourly NB2: P(high outage hour)",
"folds": len(hourly_backtest.fold_metrics),
"rows": hourly_backtest.overall_metrics.loc[0, "rows"],
"mae": hourly_backtest.overall_metrics.loc[0, "mae"],
"median_ae": hourly_backtest.overall_metrics.loc[0, "median_ae"],
"p95_ae": hourly_backtest.overall_metrics.loc[0, "p95_ae"],
"rmse": hourly_backtest.overall_metrics.loc[0, "rmse"],
"nll": hourly_backtest.overall_metrics.loc[0, "nll"],
"auc_outage": hourly_backtest.overall_metrics.loc[0, "auc_outage"],
"auc_high_outage": hourly_backtest.overall_metrics.loc[0, "auc_high_outage"],
"brier_high_outage": hourly_backtest.overall_metrics.loc[0, "brier_high_outage"],
},
])
display(html_table(pl.from_pandas(backtest_comparison.round(4))))
probability_backtest_comparison = pd.DataFrame([
{
"model": "daily NB2-implied",
"target": "high outage day",
"rows": daily_backtest.overall_metrics.loc[0, "rows"],
"positive_share": daily_backtest.overall_metrics.loc[0, "high_outage_share"],
"auc": daily_backtest.overall_metrics.loc[0, "auc_high_outage"],
"brier": daily_backtest.overall_metrics.loc[0, "brier_high_outage"],
},
{
"model": "daily logistic",
"target": "high outage day",
"rows": daily_logit_backtest.overall_metrics.loc[0, "rows"],
"positive_share": daily_logit_backtest.overall_metrics.loc[0, "positive_share"],
"auc": daily_logit_backtest.overall_metrics.loc[0, "auc"],
"brier": daily_logit_backtest.overall_metrics.loc[0, "brier"],
},
{
"model": "hourly NB2-implied",
"target": "high outage hour",
"rows": hourly_backtest.overall_metrics.loc[0, "rows"],
"positive_share": hourly_backtest.overall_metrics.loc[0, "high_outage_share"],
"auc": hourly_backtest.overall_metrics.loc[0, "auc_high_outage"],
"brier": hourly_backtest.overall_metrics.loc[0, "brier_high_outage"],
},
{
"model": "hourly logistic",
"target": "high outage hour",
"rows": hourly_logit_backtest.overall_metrics.loc[0, "rows"],
"positive_share": hourly_logit_backtest.overall_metrics.loc[0, "positive_share"],
"auc": hourly_logit_backtest.overall_metrics.loc[0, "auc"],
"brier": hourly_logit_backtest.overall_metrics.loc[0, "brier"],
},
])
display(html_table(pl.from_pandas(probability_backtest_comparison.round(4))))
fold_plot = pd.concat(
[
daily_backtest.fold_metrics.assign(model="daily NB2", auc=lambda df: df["auc_high_outage"]),
hourly_backtest.fold_metrics.assign(model="hourly NB2", auc=lambda df: df["auc_high_outage"]),
daily_logit_backtest.fold_metrics.assign(model="daily logistic", auc=lambda df: df["auc"]),
hourly_logit_backtest.fold_metrics.assign(model="hourly logistic", auc=lambda df: df["auc"]),
],
ignore_index=True,
)
fold_display_cols = ["model", "fold", "test_start", "test_end", "purge_gap", "embargo_gap", "n_test", "mae", "median_ae", "p95_ae", "rmse", "auc_outage", "auc_high_outage", "auc", "brier_high_outage", "brier"]
display(html_table(pl.from_pandas(fold_plot[[col for col in fold_display_cols if col in fold_plot.columns]].round(4))))
count_plot = fold_plot[fold_plot["model"].isin(["daily NB2", "hourly NB2"])]
fig = figure(
title="Purged WFO Count Error",
x_axis_label="Fold",
y_axis_label="MAE",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for color, model in zip(SERIES, ["daily NB2", "hourly NB2"]):
sub = count_plot[count_plot["model"].eq(model)].sort_values("fold")
fig.line(sub["fold"], sub["mae"], color=color, line_width=2, legend_label=model)
fig.scatter(sub["fold"], sub["mae"], color=color, size=7)
fig.legend.location = "top_right"
show(style(fig))
fold_models = list(dict.fromkeys(fold_plot["model"]))
fig = figure(
title="Purged WFO High-Outage Ranking",
x_axis_label="Fold",
y_axis_label="AUC",
y_range=(0.45, 0.90),
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for color, model in zip(SERIES, fold_models):
sub = fold_plot[fold_plot["model"].eq(model)].sort_values("fold")
fig.line(sub["fold"], sub["auc"], color=color, line_width=2, legend_label=model)
fig.scatter(sub["fold"], sub["auc"], color=color, size=7)
fig.add_layout(Span(location=0.50, dimension="width", line_color="#555555",
line_width=1, line_dash="dotted"))
fig.legend.location = "bottom_right"
show(style(fig))
hourly_region_backtest = hourly_backtest.region_metrics[
[
"region",
"rows",
"mean_actual",
"mean_predicted",
"mae",
"median_ae",
"p95_ae",
"rmse",
"auc_outage",
"auc_high_outage",
"brier_high_outage",
]
]
display(html_table(pl.from_pandas(hourly_region_backtest.round(4))))
region_src = ColumnDataSource(
{
"region": hourly_region_backtest["region"].to_list(),
"auc_high_outage": hourly_region_backtest["auc_high_outage"].to_list(),
}
)
fig = figure(
title="Hourly Backtest High-Outage AUC by Market Zone",
x_range=FactorRange(*hourly_region_backtest["region"].to_list()),
y_range=(0.45, 0.90),
y_axis_label="AUC",
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
fig.vbar(x="region", top="auc_high_outage", width=0.6, color=SERIES[0], source=region_src)
fig.add_layout(Span(location=0.50, dimension="width", line_color="#555555",
line_width=1, line_dash="dotted"))
fig.add_tools(
HoverTool(tooltips=[("zone", "@region"), ("AUC", "@auc_high_outage{0.000}")])
)
show(style(fig))
2 rows × 11 columns
| model | folds | rows | mae | median_ae | p95_ae | rmse | nll | auc_outage | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|
| daily NB2: P(high outage day) | 4 | 2,163 | 4.769 | 3.928 | 11.60 | 6.077 | 3.233 | — | 0.870 | 0.044 |
| hourly NB2: P(high outage hour) | 4 | 51,837 | 0.819 | 0.399 | 2.67 | 1.428 | 1.090 | 0.683 | 0.774 | 0.100 |
4 rows × 6 columns
| model | target | rows | positive_share | auc | brier |
|---|---|---|---|---|---|
| daily NB2-implied | high outage day | 2,163 | 0.056 | 0.870 | 0.044 |
| daily logistic | high outage day | 2,163 | 0.056 | 0.864 | 0.043 |
| hourly NB2-implied | high outage hour | 51,837 | 0.155 | 0.774 | 0.100 |
| hourly logistic | high outage hour | 51,837 | 0.155 | 0.778 | 0.099 |
16 rows × 16 columns
| model | fold | test_start | test_end | purge_gap | embargo_gap | n_test | mae | median_ae | p95_ae | rmse | auc_outage | auc_high_outage | auc | brier_high_outage | brier |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| daily NB2 | 1 | 2024-06-06 00:00 | 2024-12-03 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 543 | 4.360 | 3.605 | 10.81 | 5.718 | — | 0.842 | 0.842 | 0.099 | — |
| daily NB2 | 2 | 2024-12-05 00:00 | 2025-07-04 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 543 | 5.640 | 5.062 | 12.81 | 6.849 | — | 0.813 | 0.813 | 0.023 | — |
| daily NB2 | 3 | 2025-07-06 00:00 | 2026-01-02 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 543 | 4.836 | 4.065 | 11.48 | 6.133 | — | 0.874 | 0.874 | 0.035 | — |
| daily NB2 | 4 | 2026-01-04 00:00 | 2026-06-30 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 534 | 4.231 | 3.388 | 10.95 | 5.515 | — | 0.671 | 0.671 | 0.018 | — |
| hourly NB2 | 1 | 2024-06-06 02:00 | 2024-12-03 21:00 | 1 days 00:00:00 | 1 days 00:00:00 | 13,011 | 0.880 | 0.493 | 2.81 | 1.552 | 0.692 | 0.777 | 0.777 | 0.111 | — |
| hourly NB2 | 2 | 2024-12-04 22:00 | 2025-07-04 16:00 | 1 days 00:00:00 | 1 days 00:00:00 | 13,011 | 0.795 | 0.394 | 2.59 | 1.356 | 0.679 | 0.768 | 0.768 | 0.100 | — |
| hourly NB2 | 3 | 2025-07-05 17:00 | 2026-01-02 11:00 | 1 days 00:00:00 | 1 days 00:00:00 | 13,011 | 0.814 | 0.381 | 2.68 | 1.434 | 0.693 | 0.775 | 0.775 | 0.100 | — |
| hourly NB2 | 4 | 2026-01-03 12:00 | 2026-06-30 09:00 | 1 days 00:00:00 | 1 days 00:00:00 | 12,804 | 0.788 | 0.394 | 2.59 | 1.360 | 0.665 | 0.778 | 0.778 | 0.091 | — |
| daily logistic | 1 | 2024-06-06 00:00 | 2024-12-03 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 543 | — | — | — | — | — | — | 0.836 | — | 0.098 |
| daily logistic | 2 | 2024-12-05 00:00 | 2025-07-04 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 543 | — | — | — | — | — | — | 0.802 | — | 0.023 |
| daily logistic | 3 | 2025-07-06 00:00 | 2026-01-02 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 543 | — | — | — | — | — | — | 0.862 | — | 0.034 |
| daily logistic | 4 | 2026-01-04 00:00 | 2026-06-30 00:00 | 1 days 00:00:00 | 1 days 00:00:00 | 534 | — | — | — | — | — | — | 0.677 | — | 0.018 |
| hourly logistic | 1 | 2024-06-06 02:00 | 2024-12-03 21:00 | 1 days 00:00:00 | 1 days 00:00:00 | 13,011 | — | — | — | — | — | — | 0.782 | — | 0.107 |
| hourly logistic | 2 | 2024-12-04 22:00 | 2025-07-04 16:00 | 1 days 00:00:00 | 1 days 00:00:00 | 13,011 | — | — | — | — | — | — | 0.772 | — | 0.099 |
| hourly logistic | 3 | 2025-07-05 17:00 | 2026-01-02 11:00 | 1 days 00:00:00 | 1 days 00:00:00 | 13,011 | — | — | — | — | — | — | 0.778 | — | 0.099 |
| hourly logistic | 4 | 2026-01-03 12:00 | 2026-06-30 09:00 | 1 days 00:00:00 | 1 days 00:00:00 | 12,804 | — | — | — | — | — | — | 0.779 | — | 0.090 |
3 rows × 11 columns
| region | rows | mean_actual | mean_predicted | mae | median_ae | p95_ae | rmse | auc_outage | auc_high_outage | brier_high_outage |
|---|---|---|---|---|---|---|---|---|---|---|
| NP15 | 17,279 | 1.015 | 0.852 | 0.774 | 0.331 | 2.642 | 1.484 | 0.661 | 0.774 | 0.074 |
| SP15 | 17,279 | 1.097 | 1.063 | 0.877 | 0.572 | 2.680 | 1.417 | 0.700 | 0.771 | 0.127 |
| ZP26 | 17,279 | 0.962 | 0.858 | 0.806 | 0.378 | 2.668 | 1.381 | 0.675 | 0.758 | 0.100 |
On the five-year panel the daily high-outage layer ranks better than the hourly one and is steadier fold to fold than the earlier one-year reading suggested: backtest AUC 0.870 against hourly’s 0.774, ahead in three of the four folds. The two are not scored on the same event, though - a high outage day is a day above 38 starts, met on 5.6% of held-out zone-days, while a high outage hour is an hour above 2 starts, met on 15.5% of held-out zone-hours - so the Brier gap (0.044 daily against 0.100 hourly in the walk-forward, 0.026 against 0.095 on the single split) mostly reflects the rarer event, not a better model. What the comparison does support: modelling the timing directly did not buy a sharper high-outage signal, and hourly’s own claim rests on P(any outage hour), which daily cannot express at all. Daily’s weak fold is the most recent one (AUC 0.671 against hourly’s 0.778 in fold 4). The logistic benchmark does not beat the NB-implied probability at either granularity (0.864 against 0.870 daily; 0.778 against 0.774 hourly, inside noise).
8. Fifteen-Minute Feasibility Check#
A 15-minute target is constructible from curtailment start timestamps, and CAISO market data exists at 15-minute and 5-minute frequencies; the constraint is sparsity.
fifteen_panel = interval_market_zone_outage_counts(
mapped_events,
start_date=daily_panel["date"].min(),
end_date=daily_panel["date"].max(),
freq="15min",
)
fifteen_summary_display = fifteen_summary[cols].round(4)
display(html_table(pl.from_pandas(fifteen_summary_display)))
fifteen_times = pd.Index(sorted(pd.to_datetime(fifteen_panel["interval_start_utc"]).unique()))
initial_count = int(np.floor(len(fifteen_times) * 0.60))
fifteen_chunks = np.array_split(fifteen_times[initial_count:].to_numpy(), 4)
fifteen_window_rows = []
for fold, chunk in enumerate(fifteen_chunks, start=1):
if len(chunk) == 0:
continue
start = pd.Timestamp(chunk[0])
end = pd.Timestamp(chunk[-1])
sub = fifteen_panel[pd.to_datetime(fifteen_panel["interval_start_utc"]).between(start, end, inclusive="both")]
grouped = (
sub.groupby("region")["Y"]
.agg(
rows="size",
mean="mean",
variance="var",
zero_share=lambda s: (s == 0).mean(),
outage_share=lambda s: (s > 0).mean(),
p95=lambda s: s.quantile(0.95),
p99=lambda s: s.quantile(0.99),
max="max",
)
.reset_index()
)
grouped["fold"] = fold
grouped["test_start"] = start
grouped["test_end"] = end
grouped["variance_to_mean"] = grouped["variance"] / grouped["mean"].replace(0, np.nan)
fifteen_window_rows.append(grouped)
fifteen_window_summary = pd.concat(fifteen_window_rows, ignore_index=True)
display(
html_table(
pl.from_pandas(
fifteen_window_summary[
["fold", "region", "rows", "mean", "variance_to_mean", "zero_share", "outage_share", "p95", "p99", "max"]
].round(4)
)
)
)
fig = figure(
title="15-Minute Target Sparsity Across Future Windows",
x_axis_label="Future Window",
y_axis_label="Zero-Count Share",
y_range=(0.80, 0.96),
height=340,
sizing_mode="stretch_width",
toolbar_location="above",
)
for color, region in zip(SERIES, sorted(fifteen_window_summary["region"].unique())):
sub = fifteen_window_summary[fifteen_window_summary["region"].eq(region)].sort_values("fold")
fig.line(sub["fold"], sub["zero_share"], color=color, line_width=2, legend_label=region)
fig.scatter(sub["fold"], sub["zero_share"], color=color, size=7)
fig.legend.location = "bottom_right"
show(style(fig))
rt15_path = DATA_DIR / "opensource" / "caiso_lmp_real_time_15_min.parquet"
if rt15_path.exists():
rt15 = pd.read_parquet(rt15_path)
rt15_summary = pd.Series({
"rt15_rows": len(rt15),
"rt15_first_interval_utc": rt15["interval_start_utc"].min(),
"rt15_last_interval_utc": rt15["interval_start_utc"].max(),
"rt15_locations": rt15["location"].nunique(),
})
display(html_table(pl.from_pandas(rt15_summary.astype(str).to_frame("value").reset_index())))
3 rows × 12 columns
| frequency | region | rows | mean | variance | variance_to_mean | zero_share | outage_share | p90 | p95 | p99 | max |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 15min | NP15 | 176,640 | 0.256 | 2.749 | 10.72 | 0.908 | 0.092 | 0.000 | 1.000 | 8.000 | 34 |
| 15min | SP15 | 176,640 | 0.310 | 2.310 | 7.45 | 0.868 | 0.132 | 1.000 | 1.000 | 7.000 | 38 |
| 15min | ZP26 | 176,640 | 0.269 | 2.571 | 9.55 | 0.892 | 0.108 | 1.000 | 1.000 | 6.000 | 31 |
12 rows × 10 columns
| fold | region | rows | mean | variance_to_mean | zero_share | outage_share | p95 | p99 | max |
|---|---|---|---|---|---|---|---|---|---|
| 1 | NP15 | 17,664 | 0.294 | 13.31 | 0.906 | 0.094 | 1.000 | 11.37 | 33 |
| 1 | SP15 | 17,664 | 0.337 | 6.12 | 0.837 | 0.163 | 1.000 | 7.00 | 26 |
| 1 | ZP26 | 17,664 | 0.300 | 10.08 | 0.888 | 0.112 | 1.000 | 9.37 | 30 |
| 2 | NP15 | 17,664 | 0.230 | 7.89 | 0.904 | 0.096 | 1.000 | 6.00 | 30 |
| 2 | SP15 | 17,664 | 0.238 | 4.40 | 0.873 | 0.127 | 1.000 | 6.00 | 24 |
| 2 | ZP26 | 17,664 | 0.191 | 3.83 | 0.892 | 0.108 | 1.000 | 4.00 | 25 |
| 3 | NP15 | 17,664 | 0.243 | 9.06 | 0.906 | 0.094 | 1.000 | 8.00 | 27 |
| 3 | SP15 | 17,664 | 0.269 | 5.37 | 0.866 | 0.134 | 1.000 | 6.00 | 28 |
| 3 | ZP26 | 17,664 | 0.239 | 5.63 | 0.888 | 0.112 | 1.000 | 6.00 | 21 |
| 4 | NP15 | 17,664 | 0.237 | 8.24 | 0.902 | 0.098 | 1.000 | 7.00 | 22 |
| 4 | SP15 | 17,664 | 0.233 | 5.05 | 0.888 | 0.112 | 1.000 | 7.00 | 26 |
| 4 | ZP26 | 17,664 | 0.214 | 4.97 | 0.891 | 0.109 | 1.000 | 6.00 | 23 |
4 rows × 2 columns
| index | value |
|---|---|
| rt15_rows | 369394 |
| rt15_first_interval_utc | 2023-05-11 07:00:00 |
| rt15_last_interval_utc | 2026-08-27 06:45:00 |
| rt15_locations | 5 |
About nine out of ten zone-intervals are zero, which points to a hurdle, zero-inflated, or ranking model rather than plain NB2, especially while outage timestamp conventions are unverified.
9. Findings#
Daily is the baseline for count intensity and, on the five-year panel, the better high-outage ranker: walk-forward AUC 0.870 against hourly’s 0.774, ahead in three of four folds. The two score different events (a day above 38 starts, an hour above 2), so the Brier gap between them is not a like-for-like comparison.
Modelling the timing directly did not buy a sharper high-outage signal. What hourly adds is P(any outage hour), an event daily cannot express because no held-out zone-day is ever zero.
Daily’s weakness is recency, not level: fold 4 falls to 0.671 while hourly holds 0.778.
Logistic benchmarks are a direct binary alternative to the NB-implied probabilities and do not beat them at either granularity.
The daily fit’s dispersion is zero, so the daily NB2 is a Poisson in practice; only the hourly target carries residual overdispersion (alpha 0.45).
Fifteen-minute counts are constructible from outage timestamps but sparse enough that the model class should change.
Mapping remains the main sampling risk; higher-confidence resource-to-zone mapping should raise coverage.