Step 4 - The ARIMA 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 prophet_xgb). Step 3 shows where a final set of hyperparameters comes from; here we fit the model with the package defaults from params.yaml and inspect its behavior. To run this notebook with tuned parameters instead, pass config_overrides={"models": {"arima": best_params}} into run_single_its().

Goal: walk through run_single_its() using the ARIMAModel.

Sections:

• 4a. Load the pre-built dummy data.
• 4b. Fit ARIMAModel manually and inspect the FitResult.
• 4c. Inside ARIMA -- automatic order selection and in-sample fit.
• 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]:

%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")
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") 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]:

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 inspect the FitResult.

In [3]:

from its2s.data_prep import prepare_splits
from its2s.models.arima import ARIMAModel
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, "arima")
model = ARIMAModel(params=model_params)

print("Fitting ARIMAModel 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"  metadata: {fit_result.metadata}")

from its2s.data_prep import prepare_splits from its2s.models.arima import ARIMAModel 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, "arima") model = ARIMAModel(params=model_params) print("Fitting ARIMAModel 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" metadata: {fit_result.metadata}")

Fitting ARIMAModel 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.0166  std=2.1331
  metadata: {'order': (2, 0, 1), 'seasonal_order': (1, 0, 1, 7)}

4c. Inside ARIMA -- automatic order selection and in-sample fit¶

auto_arima selects the best (p,d,q) and seasonal order via stepwise search during the initial fit. The discovered order is stored in fit_result.metadata and preserved by clone_fresh(), so Moving Block Bootstrap refits use the same model structure without repeating the expensive search on each simulation.

In [4]:

arima_model = fit_result.model_object

print("Discovered ARIMA order:")
print(f"  Non-seasonal (p,d,q)  : {fit_result.metadata['order']}")
print(f"  Seasonal (P,D,Q,m)    : {fit_result.metadata['seasonal_order']}")
print()
print(arima_model.summary())

arima_model = fit_result.model_object print("Discovered ARIMA order:") print(f" Non-seasonal (p,d,q) : {fit_result.metadata['order']}") print(f" Seasonal (P,D,Q,m) : {fit_result.metadata['seasonal_order']}") print() print(arima_model.summary())

Discovered ARIMA order:
  Non-seasonal (p,d,q)  : (2, 0, 1)
  Seasonal (P,D,Q,m)    : (1, 0, 1, 7)

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                                 y   No. Observations:                 1169
Model:             SARIMAX(2, 0, 1)x(1, 0, 1, 7)   Log Likelihood               -2544.364
Date:                           Sat, 01 Aug 2026   AIC                           5102.727
Time:                                   11:11:54   BIC                           5138.175
Sample:                                        0   HQIC                          5116.097
                                          - 1169                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.0273      0.018      1.482      0.138      -0.009       0.063
ar.L1          1.0515      0.037     28.564      0.000       0.979       1.124
ar.L2         -0.0552      0.036     -1.524      0.128      -0.126       0.016
ma.L1         -0.8038      0.024    -33.259      0.000      -0.851      -0.756
ar.S.L7        0.8623      0.078     11.028      0.000       0.709       1.016
ma.S.L7       -0.7903      0.092     -8.616      0.000      -0.970      -0.610
sigma2         4.5021      0.184     24.466      0.000       4.141       4.863
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):                 0.32
Prob(Q):                              0.96   Prob(JB):                         0.85
Heteroskedasticity (H):               1.03   Skew:                            -0.04
Prob(H) (two-sided):                  0.78   Kurtosis:                         3.02
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

In [5]:

# clone_fresh() preserves the discovered order -- verify
fresh_clone = model.clone_fresh()
print("clone_fresh() preserves order for MBB:")
print(f"  Original  _fixed_order         : {model._fixed_order}")
print(f"  Clone     _fixed_order         : {fresh_clone._fixed_order}")
print(f"  Original  _fixed_seasonal_order: {model._fixed_seasonal_order}")
print(f"  Clone     _fixed_seasonal_order: {fresh_clone._fixed_seasonal_order}")
print()
print("When MBB fits the clone, it skips auto_arima and goes directly to pm.ARIMA(order=...)")

# clone_fresh() preserves the discovered order -- verify fresh_clone = model.clone_fresh() print("clone_fresh() preserves order for MBB:") print(f" Original _fixed_order : {model._fixed_order}") print(f" Clone _fixed_order : {fresh_clone._fixed_order}") print(f" Original _fixed_seasonal_order: {model._fixed_seasonal_order}") print(f" Clone _fixed_seasonal_order: {fresh_clone._fixed_seasonal_order}") print() print("When MBB fits the clone, it skips auto_arima and goes directly to pm.ARIMA(order=...)")

clone_fresh() preserves order for MBB:
  Original  _fixed_order         : (2, 0, 1)
  Clone     _fixed_order         : (2, 0, 1)
  Original  _fixed_seasonal_order: (1, 0, 1, 7)
  Clone     _fixed_seasonal_order: (1, 0, 1, 7)

When MBB fits the clone, it skips auto_arima and goes directly to pm.ARIMA(order=...)

In [6]:

fig, ax = plt.subplots(figsize=(13, 3.5))
ax.plot(splits.train_df["ds"], fit_result.residuals,
        linewidth=0.6, color="#4C72B0", alpha=0.7)
ax.axhline(0, color="black", linewidth=0.8, linestyle="--")
ax.set_title("ARIMA residuals (y - fitted)", fontsize=10)
ax.set_ylabel("Residual")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
plt.tight_layout()
plt.savefig(OUT_DIR / "arima_residuals.png", dpi=150)
display(fig)

fig, ax = plt.subplots(figsize=(13, 3.5)) ax.plot(splits.train_df["ds"], fit_result.residuals, linewidth=0.6, color="#4C72B0", alpha=0.7) ax.axhline(0, color="black", linewidth=0.8, linestyle="--") ax.set_title("ARIMA residuals (y - fitted)", fontsize=10) ax.set_ylabel("Residual") ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y")) plt.tight_layout() plt.savefig(OUT_DIR / "arima_residuals.png", dpi=150) display(fig)

4d. Full pipeline run¶

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

In [7]:

from its2s import run_single_its

result = run_single_its(
    df=df,
    intervention_date=INTERVENTION,
    model_name="arima",
    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="arima", 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       : arima
  fit_result       : FitResult with 1169 fitted values
  bootstrap_result : BootstrapCIResult  pred_matrix shape = (407, 100)
  metrics_train    : MetricsResult(rmse=8.64647436271905, mae=7.262726163373645, mape=15.18148711267987, mase=None, mase_m=7, mase_denominator=None)
  metrics_test     : MetricsResult(rmse=9.071485233554489, mae=7.3625843280321535, mape=15.050509066174161, mase=3.202853892720027, 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 [8]:

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  8.646  7.263  15.181    NaN
Test   9.071  7.363  15.051  3.203

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

In [9]:

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     2379.594922    650.205856     517.42889    755.251793   27.324224

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 56.877922       54.305913       59.945263 19.571778     16.504437     22.143787   34.410150         29.017300         38.932131
2022-03-16 70.565458 56.866622       54.298199       59.924017 13.698836     10.641441     16.267259   24.089415         18.712982         28.605989
2022-03-17 74.032036 56.855252       54.290277       59.915793 17.176784     14.116242     19.741759   30.211429         24.828388         34.722841
2022-03-18 69.205309 56.844260       54.282442       59.907675 12.361049      9.297633     14.922867   21.745465         16.356328         26.252197
2022-03-19 70.763015 56.832941       54.274828       59.913454 13.930074     10.849561     16.488187   24.510563         19.090268         29.011672
2022-03-20 75.696935 56.821815       54.267083       59.917739 18.875120     15.779196     21.429852   33.218088         27.769610         37.714128
2022-03-21 72.180154 56.810573       54.259288       59.913658 15.369581     12.266496     17.920866   27.054086         21.591925         31.544949
2022-03-22 74.416647 56.799459       54.251537       59.877442 17.617188     14.539205     20.165109   31.016471         25.597436         35.502291
2022-03-23 73.242643 56.788456       54.244116       59.887226 16.454188     13.355418     18.998528   28.974529         23.517839         33.454912
2022-03-24 72.587821 56.777396       54.236519       59.884320 15.810425     12.703501     18.351302   27.846337         22.374223         32.321493

In [10]:

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 650.205856 517.428890 755.251793     42
Mean ATE per obs  15.481092  12.319735  17.982186     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 [11]:

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"ARIMA 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 / "arima_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"ARIMA 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 / "arima_counterfactual.png", dpi=150) display(fig)

Known issue -- smooth counterfactual curve. The counterfactual line above appears as a smooth, monotonically declining curve rather than tracking the annual seasonal pattern visible in the observed data. This happens because params.yaml sets m: 7, so auto_arima searches for weekly (7-day) seasonal structure only. The dominant annual cycle (period ~ 365 days) is invisible to the model, and multi-step ARIMA forecasts without that component converge toward the long-run trend. To capture annual seasonality, set m=365 via config_overrides={"models": {"arima": {"m": 365}}} — note that this increases fitting time substantially.

Key takeaways¶

1. ARIMAModel.fit() calls auto_arima on the first run to discover (p,d,q) and seasonal order; the selected order is stored in fit_result.metadata.
2. clone_fresh() preserves the discovered order so Moving Block Bootstrap refits do not repeat the expensive stepwise search on each simulation.
3. Unlike Prophet-based models, ARIMA has no decomposition stages -- residuals are simply y - fitted.
4. run_single_its() orchestrates: load_config -> prepare_splits -> fit -> bootstrap -> metrics -> excess -> save.
5. Excess = observed - counterfactual_predicted. With a true +8/day effect over 42 days, total excess should land near 336 (noise aside).