Step 4 - The Prophet + XGB hybrid pipeline¶

This is one of three step4_model_* notebooks, one per model architecture. Read them after step 1 (data splitting), step 2 (cross-validation), and step 3 (hyperparameter tuning -- demoed on this very model). Step 3 shows where a final set of hyperparameters comes from; here we fit with the package defaults from params.yaml. To use the tuned parameters from step 3, pass config_overrides={"models": {"prophet_xgb": tuning_result.best_params}} into run_single_its().

Goal: walk through run_single_its() using the ProphetXGBHybridModel.

Sections:

4a. Load the pre-built dummy data.
4b. Fit ProphetXGBHybridModel manually and inspect the FitResult.
4c. Understand the two-stage fit: Prophet for trend+seasonality, XGB for residuals.
4d. Run the full pipeline via run_single_its().
4e. Inspect PipelineResult: metrics, excess table, ATE.
4f. Reproduce the counterfactual plot with annotations.

In [1]:

Copied!





%matplotlib inline

from IPython.display import display
import logging
import warnings
from pathlib import Path

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

warnings.filterwarnings("ignore")  # suppress Prophet / XGB verbosity
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s  %(levelname)s  %(name)s - %(message)s",
    datefmt="%H:%M:%S",
)
logging.getLogger("cmdstanpy").setLevel(logging.WARNING)
logging.getLogger("its2s").setLevel(logging.WARNING)

OUT_DIR = Path.cwd() / "figures"
OUT_DIR.mkdir(exist_ok=True)
INTERVENTION = "2022-03-15"
TEST_DAYS    = 365
HOLDOUT_DAYS = 42
%matplotlib inline

from IPython.display import display
import logging
import warnings
from pathlib import Path

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

warnings.filterwarnings("ignore")  # suppress Prophet / XGB verbosity
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s  %(levelname)s  %(name)s - %(message)s",
    datefmt="%H:%M:%S",
)
logging.getLogger("cmdstanpy").setLevel(logging.WARNING)
logging.getLogger("its2s").setLevel(logging.WARNING)

OUT_DIR = Path.cwd() / "figures"
OUT_DIR.mkdir(exist_ok=True)
INTERVENTION = "2022-03-15"
TEST_DAYS    = 365
HOLDOUT_DAYS = 42

4a. Load the pre-built dummy data¶

The series has a +8/day intervention effect baked in for 42 days after 2022-03-15.

In [2]:

Copied!





df = pd.read_csv("data/dummy_data.csv", parse_dates=["ds"])

print("=" * 60)
print("Dummy dataset (with +8/day intervention effect)")
print("=" * 60)
print(df.tail())
df = pd.read_csv("data/dummy_data.csv", parse_dates=["ds"])

print("=" * 60)
print("Dummy dataset (with +8/day intervention effect)")
print("=" * 60)
print(df.tail())

============================================================
Dummy dataset (with +8/day intervention effect)
============================================================
             ds          y  covar_linear  covar_dow  covar_noise
1571 2022-04-21  71.975041      0.995182        3.0    -0.356611
1572 2022-04-22  73.042307      1.033630        4.0     0.247103
1573 2022-04-23  72.275013      1.003513        5.0     1.129482
1574 2022-04-24  69.332534      1.004775        6.0    -0.321536
1575 2022-04-25  69.651058      0.981820        0.0    -1.057655

4b. Manual fit¶

This replicates what run_single_its does internally, so we can poke at the FitResult.

In [3]:

Copied!





from its2s.data_prep import prepare_splits
from its2s.models.prophet_xgb import ProphetXGBHybridModel
from its2s.settings import get_model_config, load_config

config = load_config()
splits = prepare_splits(df, INTERVENTION, split_method="days", test_days=TEST_DAYS, holdout_days=HOLDOUT_DAYS)

model_params = get_model_config(config, "prophet_xgb")
model = ProphetXGBHybridModel(params=model_params)

print("Fitting ProphetXGBHybridModel on training data ...")
print(f"  Training rows : {len(splits.train_df)}")
print(f"  Training range: {splits.train_df['ds'].min().date()} -> {splits.train_df['ds'].max().date()}")

fit_result = model.fit(splits.train_df, target_col="y", date_col="ds")

print("\nFitResult fields:")
print(f"  fitted_values  shape = {fit_result.fitted_values.shape}")
print(f"  residuals      shape = {fit_result.residuals.shape}")
print(f"  residuals  mean={fit_result.residuals.mean():.4f}  std={fit_result.residuals.std():.4f}")
print(f"  model_object keys: {list(fit_result.model_object.keys())}")
from its2s.data_prep import prepare_splits
from its2s.models.prophet_xgb import ProphetXGBHybridModel
from its2s.settings import get_model_config, load_config

config = load_config()
splits = prepare_splits(df, INTERVENTION, split_method="days", test_days=TEST_DAYS, holdout_days=HOLDOUT_DAYS)

model_params = get_model_config(config, "prophet_xgb")
model = ProphetXGBHybridModel(params=model_params)

print("Fitting ProphetXGBHybridModel on training data ...")
print(f"  Training rows : {len(splits.train_df)}")
print(f"  Training range: {splits.train_df['ds'].min().date()} -> {splits.train_df['ds'].max().date()}")

fit_result = model.fit(splits.train_df, target_col="y", date_col="ds")

print("\nFitResult fields:")
print(f"  fitted_values  shape = {fit_result.fitted_values.shape}")
print(f"  residuals      shape = {fit_result.residuals.shape}")
print(f"  residuals  mean={fit_result.residuals.mean():.4f}  std={fit_result.residuals.std():.4f}")
print(f"  model_object keys: {list(fit_result.model_object.keys())}")

Fitting ProphetXGBHybridModel on training data ...
  Training rows : 1169
  Training range: 2018-01-01 -> 2021-03-14

FitResult fields:
  fitted_values  shape = (1169,)
  residuals      shape = (1169,)
  residuals  mean=0.0018  std=1.5274
  model_object keys: ['prophet', 'xgb']

4c. Inside the hybrid - Prophet components vs XGB residuals¶

Stage 1 (Prophet) captures trend + seasonality. Stage 2 (XGB) fits the residuals using time features (day_of_week, day_of_year, month, week_of_year). The stage-2 residuals are what drives the Moving Block Bootstrap in step 3.

In [4]:

Copied!





prophet_model = fit_result.model_object["prophet"]
xgb_model     = fit_result.model_object["xgb"]

prophet_input = splits.train_df[["ds"]].copy()
prophet_input["ds"] = pd.to_datetime(prophet_input["ds"])
prophet_pred = prophet_model.predict(prophet_input)

prophet_yhat     = prophet_pred["yhat"].values
raw_y            = splits.train_df["y"].values
stage1_residuals = raw_y - prophet_yhat     # what XGB was trained on
final_residuals  = fit_result.residuals     # what remained after XGB too

print(f"  Stage-1 residuals (y - Prophet)          std = {stage1_residuals.std():.4f}")
print(f"  Stage-2 residuals (y - Prophet - XGB)    std = {final_residuals.std():.4f}")
print("  -> XGB reduces residual variance from stage 1 to stage 2")

xgb_feat_names = ["day_of_week", "day_of_year", "month", "week_of_year"]
importances = pd.Series(xgb_model.feature_importances_, index=xgb_feat_names).sort_values(ascending=False)
print("\n  XGB feature importances:")
print(importances.to_string())
prophet_model = fit_result.model_object["prophet"]
xgb_model     = fit_result.model_object["xgb"]

prophet_input = splits.train_df[["ds"]].copy()
prophet_input["ds"] = pd.to_datetime(prophet_input["ds"])
prophet_pred = prophet_model.predict(prophet_input)

prophet_yhat     = prophet_pred["yhat"].values
raw_y            = splits.train_df["y"].values
stage1_residuals = raw_y - prophet_yhat     # what XGB was trained on
final_residuals  = fit_result.residuals     # what remained after XGB too

print(f"  Stage-1 residuals (y - Prophet)          std = {stage1_residuals.std():.4f}")
print(f"  Stage-2 residuals (y - Prophet - XGB)    std = {final_residuals.std():.4f}")
print("  -> XGB reduces residual variance from stage 1 to stage 2")

xgb_feat_names = ["day_of_week", "day_of_year", "month", "week_of_year"]
importances = pd.Series(xgb_model.feature_importances_, index=xgb_feat_names).sort_values(ascending=False)
print("\n  XGB feature importances:")
print(importances.to_string())

  Stage-1 residuals (y - Prophet)          std = 1.9444
  Stage-2 residuals (y - Prophet - XGB)    std = 1.5274
  -> XGB reduces residual variance from stage 1 to stage 2

  XGB feature importances:
day_of_year     0.267140
day_of_week     0.248502
week_of_year    0.246037
month           0.238320

In [5]:

Copied!





fig, axes = plt.subplots(1, 2, figsize=(13, 3.5))
for ax, resid, title, color in [
    (axes[0], stage1_residuals, "Stage-1 residuals (y - Prophet)", "#4C72B0"),
    (axes[1], final_residuals,  "Stage-2 residuals (y - Prophet - XGB)", "#DD8452"),
]:
    ax.plot(splits.train_df["ds"], resid, linewidth=0.6, color=color, alpha=0.7)
    ax.axhline(0, color="black", linewidth=0.8, linestyle="--")
    ax.set_title(title, fontsize=10)
    ax.set_ylabel("Residual")
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
fig.suptitle("Prophet+XGB: residuals before and after XGB stage", fontsize=11)
plt.tight_layout()
plt.savefig(OUT_DIR / "step4a_residuals.png", dpi=150)
display(fig)
fig, axes = plt.subplots(1, 2, figsize=(13, 3.5))
for ax, resid, title, color in [
    (axes[0], stage1_residuals, "Stage-1 residuals (y - Prophet)", "#4C72B0"),
    (axes[1], final_residuals,  "Stage-2 residuals (y - Prophet - XGB)", "#DD8452"),
]:
    ax.plot(splits.train_df["ds"], resid, linewidth=0.6, color=color, alpha=0.7)
    ax.axhline(0, color="black", linewidth=0.8, linestyle="--")
    ax.set_title(title, fontsize=10)
    ax.set_ylabel("Residual")
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
fig.suptitle("Prophet+XGB: residuals before and after XGB stage", fontsize=11)
plt.tight_layout()
plt.savefig(OUT_DIR / "step4a_residuals.png", dpi=150)
display(fig)

No description has been provided for this image

4d. Full pipeline run¶

MBB bootstrap runs with n_sim=100 for speed; production use should set this to 1000+.

In [6]:

Copied!





from its2s import run_single_its

result = run_single_its(
    df=df,
    intervention_date=INTERVENTION,
    model_name="prophet_xgb",
    config_overrides={
        "bootstrap": {"n_sim": 100},
        "periods":   {"split_method": "days", "test_days": TEST_DAYS, "holdout_days": HOLDOUT_DAYS},
    },
    output_dir=OUT_DIR,
    seed=42,
)

print("PipelineResult fields:")
print(f"  model_name       : {result.model_name}")
print(f"  fit_result       : FitResult with {len(result.fit_result.fitted_values)} fitted values")
print(f"  bootstrap_result : BootstrapCIResult  pred_matrix shape = {result.bootstrap_result.pred_matrix.shape}")
print(f"  metrics_train    : {result.metrics_train}")
print(f"  metrics_test     : {result.metrics_test}")
from its2s import run_single_its

result = run_single_its(
    df=df,
    intervention_date=INTERVENTION,
    model_name="prophet_xgb",
    config_overrides={
        "bootstrap": {"n_sim": 100},
        "periods":   {"split_method": "days", "test_days": TEST_DAYS, "holdout_days": HOLDOUT_DAYS},
    },
    output_dir=OUT_DIR,
    seed=42,
)

print("PipelineResult fields:")
print(f"  model_name       : {result.model_name}")
print(f"  fit_result       : FitResult with {len(result.fit_result.fitted_values)} fitted values")
print(f"  bootstrap_result : BootstrapCIResult  pred_matrix shape = {result.bootstrap_result.pred_matrix.shape}")
print(f"  metrics_train    : {result.metrics_train}")
print(f"  metrics_test     : {result.metrics_test}")

PipelineResult fields:
  model_name       : prophet_xgb
  fit_result       : FitResult with 1169 fitted values
  bootstrap_result : BootstrapCIResult  pred_matrix shape = (407, 100)
  metrics_train    : MetricsResult(rmse=1.527386354794742, mae=1.2009044881186752, mape=2.360806292904574, mase=None, mase_m=7, mase_denominator=None)
  metrics_test     : MetricsResult(rmse=2.318135012278161, mae=1.8528514238121523, mape=3.4921501819923497, mase=0.8060230118919026, mase_m=7, mase_denominator=2.298757475252632)

4e. Metrics and excess table¶

MASE is reported for the held-out test window only (the Train cell is NaN by design): it is the ratio of the model's MAE to the in-sample MAE of the seasonal-naive forecast at the resolved period m (mase_m and mase_denominator on the result). A value below 1 means the model beats the seasonal-naive benchmark. The test window serves as a single-use adequacy check of the fitted model, not a retuning target.

In [7]:

Copied!





metrics_df = pd.DataFrame({
    "RMSE":  [result.metrics_train.rmse,  result.metrics_test.rmse],
    "MAE":   [result.metrics_train.mae,   result.metrics_test.mae],
    "MAPE":  [result.metrics_train.mape,  result.metrics_test.mape],
    "MASE":  [result.metrics_train.mase,  result.metrics_test.mase],
}, index=["Train", "Test"])
print(metrics_df.round(3).to_string())

mt = result.metrics_test
print(f"\nMASE benchmark: in-sample seasonal-naive MAE = "
      f"{mt.mase_denominator:.3f} at m = {mt.mase_m}")
metrics_df = pd.DataFrame({
    "RMSE":  [result.metrics_train.rmse,  result.metrics_test.rmse],
    "MAE":   [result.metrics_train.mae,   result.metrics_test.mae],
    "MAPE":  [result.metrics_train.mape,  result.metrics_test.mape],
    "MASE":  [result.metrics_train.mase,  result.metrics_test.mase],
}, index=["Train", "Test"])
print(metrics_df.round(3).to_string())

mt = result.metrics_test
print(f"\nMASE benchmark: in-sample seasonal-naive MAE = "
      f"{mt.mase_denominator:.3f} at m = {mt.mase_m}")

        RMSE    MAE   MAPE   MASE
Train  1.527  1.201  2.361    NaN
Test   2.318  1.853  3.492  0.806

MASE benchmark: in-sample seasonal-naive MAE = 2.299 at m = 7

In [8]:

Copied!





print("Period-level excess:")
print(result.excess_table.period_excess.to_string(index=False))
print("\nDaily excess - first 10 holdout days:")
print(result.excess_table.obs_excess.head(10).to_string(index=False))
print("Period-level excess:")
print(result.excess_table.period_excess.to_string(index=False))
print("\nDaily excess - first 10 holdout days:")
print(result.excess_table.obs_excess.head(10).to_string(index=False))

Period-level excess:
      period start_date   end_date  n_obs  total_observed  total_expected  total_excess  excess_ci_lo  excess_ci_hi  excess_pct
Full holdout 2022-03-15 2022-04-25     42     3029.800778     2727.650623    302.150155      282.0523    337.650058   11.077304

Daily excess - first 10 holdout days:
      date  observed  expected  expected_ci_lo  expected_ci_hi    excess  excess_ci_lo  excess_ci_hi  excess_pct  excess_pct_ci_lo  excess_pct_ci_hi
2022-03-15 76.449700 63.998524       63.053132       65.508686 12.451176     10.941014     13.396569   19.455411         17.095729         20.932621
2022-03-16 70.565458 64.569162       63.231448       65.523850  5.996296      5.041608      7.334010    9.286625          7.808074         11.358378
2022-03-17 74.032036 64.508181       63.478247       65.659495  9.523855      8.372541     10.553789   14.763794         12.979037         16.360388
2022-03-18 69.205309 64.976148       63.787714       65.997934  4.229161      3.207375      5.417595    6.508790          4.936234          8.337821
2022-03-19 70.763015 65.717400       64.349104       66.509370  5.045615      4.253645      6.413911    7.677746          6.472631          9.759837
2022-03-20 75.696935 64.939114       63.869905       66.179721 10.757821      9.517214     11.827030   16.566012         14.655596         18.212491
2022-03-21 72.180154 65.828221       64.440184       66.706735  6.351933      5.473419      7.739970    9.649256          8.314699         11.757829
2022-03-22 74.416647 65.021015       63.558993       66.067205  9.395632      8.349442     10.857654   14.450147         12.841144         16.698685
2022-03-23 73.242643 64.778175       63.509081       65.689452  8.464468      7.553191      9.733563   13.066852         11.660086         15.025991
2022-03-24 72.587821 64.554763       63.225188       65.546803  8.033058      7.041018      9.362634   12.443789         10.907047         14.503397

In [9]:

Copied!





from its2s.metrics.excess import calc_ate_summary

ate = calc_ate_summary(result.excess_table.obs_excess)
print("Average Treatment Effect (ATE) summary:")
print(ate.to_string(index=False))
print("\n  Total ATE      = sum of daily excess over full holdout")
print("  Mean ATE per obs = average excess per observation")
print(f"  Simulated effect was +8/day for {HOLDOUT_DAYS} days -> expected total excess ~{8 * HOLDOUT_DAYS}")
from its2s.metrics.excess import calc_ate_summary

ate = calc_ate_summary(result.excess_table.obs_excess)
print("Average Treatment Effect (ATE) summary:")
print(ate.to_string(index=False))
print("\n  Total ATE      = sum of daily excess over full holdout")
print("  Mean ATE per obs = average excess per observation")
print(f"  Simulated effect was +8/day for {HOLDOUT_DAYS} days -> expected total excess ~{8 * HOLDOUT_DAYS}")

Average Treatment Effect (ATE) summary:
          metric   estimate      ci_lo      ci_hi  n_obs
       Total ATE 302.150155 257.378400 356.434864     42
Mean ATE per obs   7.194051   6.128057   8.486544     42

  Total ATE      = sum of daily excess over full holdout
  Mean ATE per obs = average excess per observation
  Simulated effect was +8/day for 42 days -> expected total excess ~336

4f. Counterfactual plot (annotated)¶

In [10]:

Copied!





br = result.bootstrap_result
pred_dates = pd.to_datetime(br.dates)
intervention_ts = pd.Timestamp(INTERVENTION)

fig, ax = plt.subplots(figsize=(14, 5))

for part in [splits.train_df, splits.test_df, splits.holdout_df]:
    ax.plot(part["ds"], part["y"], color="#333333", linewidth=0.6, alpha=0.7)
ax.plot([], [], color="#333333", linewidth=0.6, alpha=0.7, label="Observed")

ax.plot(pred_dates, br.predicted, color="#B2182B", linewidth=1.4,
        label="Counterfactual (no-intervention)")
ax.fill_between(pred_dates, br.conf_lo, br.conf_hi,
                color="#B2182B", alpha=0.15, label="95% CI (MBB)")

ax.axvspan(intervention_ts, splits.holdout_df["ds"].max(),
           color="#FEE08B", alpha=0.25, label="Holdout (post-intervention)")
ax.axvline(intervention_ts, color="#4DAF4A", linestyle="--", linewidth=1.3,
           label="Intervention date")

last_date = pred_dates[pred_dates >= intervention_ts][-1]
last_obs  = splits.holdout_df.loc[splits.holdout_df["ds"] == last_date, "y"].values
last_pred = br.predicted[pred_dates == last_date]
if len(last_obs) and len(last_pred):
    ax.annotate(
        f"Excess ~ {float(last_obs[0] - last_pred[0]):.1f}",
        xy=(last_date, float(last_pred[0])),
        xytext=(last_date - pd.Timedelta(days=90), float(last_pred[0]) + 6),
        arrowprops=dict(arrowstyle="->", color="black"),
        fontsize=9,
    )

ax.set_xlabel("Date")
ax.set_ylabel("y (daily outcome)")
ax.set_title(
    f"Prophet+XGB counterfactual  |  Test RMSE: {result.metrics_test.rmse:.2f}"
    f"  |  Test MAPE: {result.metrics_test.mape:.1f}%",
    fontsize=10,
)
ax.legend(loc="upper left", fontsize=8)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
plt.tight_layout()
plt.savefig(OUT_DIR / "step4a_counterfactual.png", dpi=150)
display(fig)
br = result.bootstrap_result
pred_dates = pd.to_datetime(br.dates)
intervention_ts = pd.Timestamp(INTERVENTION)

fig, ax = plt.subplots(figsize=(14, 5))

for part in [splits.train_df, splits.test_df, splits.holdout_df]:
    ax.plot(part["ds"], part["y"], color="#333333", linewidth=0.6, alpha=0.7)
ax.plot([], [], color="#333333", linewidth=0.6, alpha=0.7, label="Observed")

ax.plot(pred_dates, br.predicted, color="#B2182B", linewidth=1.4,
        label="Counterfactual (no-intervention)")
ax.fill_between(pred_dates, br.conf_lo, br.conf_hi,
                color="#B2182B", alpha=0.15, label="95% CI (MBB)")

ax.axvspan(intervention_ts, splits.holdout_df["ds"].max(),
           color="#FEE08B", alpha=0.25, label="Holdout (post-intervention)")
ax.axvline(intervention_ts, color="#4DAF4A", linestyle="--", linewidth=1.3,
           label="Intervention date")

last_date = pred_dates[pred_dates >= intervention_ts][-1]
last_obs  = splits.holdout_df.loc[splits.holdout_df["ds"] == last_date, "y"].values
last_pred = br.predicted[pred_dates == last_date]
if len(last_obs) and len(last_pred):
    ax.annotate(
        f"Excess ~ {float(last_obs[0] - last_pred[0]):.1f}",
        xy=(last_date, float(last_pred[0])),
        xytext=(last_date - pd.Timedelta(days=90), float(last_pred[0]) + 6),
        arrowprops=dict(arrowstyle="->", color="black"),
        fontsize=9,
    )

ax.set_xlabel("Date")
ax.set_ylabel("y (daily outcome)")
ax.set_title(
    f"Prophet+XGB counterfactual  |  Test RMSE: {result.metrics_test.rmse:.2f}"
    f"  |  Test MAPE: {result.metrics_test.mape:.1f}%",
    fontsize=10,
)
ax.legend(loc="upper left", fontsize=8)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
plt.tight_layout()
plt.savefig(OUT_DIR / "step4a_counterfactual.png", dpi=150)
display(fig)

Key takeaways¶

ProphetXGBHybridModel.fit() runs two models in sequence: Prophet for trend+seasonality, then XGB on y - Prophet_yhat.
The FitResult.residuals are the raw material for the Moving Block Bootstrap in step 5.
run_single_its() orchestrates: load_config -> prepare_splits -> fit -> bootstrap -> metrics -> excess -> save.
Excess = observed - counterfactual_predicted. With a true +8/day effect over 42 days, total excess should land near 336 (noise aside).