ARIMA vs. XGBoost vs. LightGBM: A Full Walk-Forward Comparison

A classical ARIMA baseline was tested, honestly, against XGBoost and LightGBM on real quarterly loan-delinquency data: same walk-forward evaluation for every model, real significance testing, and no result trusted until it survived scrutiny. A leakage bug in the first ML implementation produced an implausibly good early result; once fixed, ARIMA won every comparison, and SHAP explains why.

1. ARIMA beat XGBoost and LightGBM on all five FRED loan-delinquency series, walk-forward backtested with an identical evaluation protocol for every model.
2. A real data-leakage bug was caught before any ML result was trusted: a boundary condition in the direct-horizon target construction let a handful of training rows see their own test-window answer. Fixed and covered by a regression test.
3. ARIMA’s advantage is statistically significant in 27 of 40 series/model/horizon comparisons (Diebold-Mariano, Harvey-Leybourne-Newbold corrected), universal at the 1-quarter horizon and fading at longer horizons.
4. SHAP interpretability shows why: in all 10 of 10 series/model combinations, the single most recent lag carries 47–59% of total feature importance: both approaches lean on essentially the same signal, and ARIMA simply expresses it more efficiently on ~140 rows than a tree ensemble can.

Five accounts, one evaluation rule for every model

Data: FRED, quarterly, 1991-01-01 through 2026-04-01 (142 observations per series), truncated to that range at pull time.

Series

DRALACBS (All Loans), DRCCLACBS (Credit Card), DRBLACBS (Business), DRSFRMACBS (Mortgage), DRCRELEXFACBS (CRE)

Backtest window

min_train_size = 80, h = 4 quarters ahead, expanding window, refit every fold: 59 graded folds per series

Benchmarks

Seasonal-naive (lag-4) and simple-naive (lag-1)

Scoring metric

MASE: mean absolute error of the model, divided by mean absolute error of the naive benchmark, averaged per fold

Below 1.0 means the model beats the naive benchmark. Every model in this study (ARIMA, XGBoost, and LightGBM) is scored on the exact same 59 folds per series, so every number that follows is directly comparable across models.

Stationarity, autocorrelation, decomposition

ADF (H₀: unit root) and KPSS (H₀: stationary) were run on each raw series, at the conventional α = 0.05:

A conflicting ADF/KPSS verdict typically indicates trend-stationary behavior rather than test failure: real short-run mean reversion riding on top of longer, crisis-driven trend moves. ACF/PACF showed the classic AR signature (sharp PACF cutoff after lag 1–2, gradually decaying ACF) in every series, motivating AR(2) as the initial order guess for All Loans/Business/Credit Card/Mortgage and AR(1) for CRE. STL decomposition (period=4) confirmed trend dominates seasonality by an order of magnitude in every series, with both seasonal amplitude and residual variance expanding specifically during crisis windows, evidence of heteroscedasticity that motivated using conformal, error-history-based prediction intervals for the ML models rather than a fixed-variance assumption (see Prediction Intervals, below).

Full EDA write-up, all five ACF/PACF and STL plots: docs/phase1_eda.md in the project repository.

Per-series order search, diagnostics, backtest

Per series: a 5-candidate order grid built around the AR guess: (p,0,0), (p,1,0), (p,0,1), (p,1,1), (p+1,1,0), fit by maximum likelihood, lowest-AIC order selected, then walk-forward backtested.

All five pass Ljung-Box comfortably at every lag, no significant residual autocorrelation left unmodeled. One caveat worth stating plainly: by the conventional Burnham & Anderson rule of thumb (ΔAIC < 2 implies no decisive separation), four of the five winning orders are not decisively separated from their runner-up. CRE is the exception, and not on the top-two orders: it’s d=1 vs. d=0. ΔAIC of over 220, an unambiguous signal despite CRE’s formally ambiguous ADF/KPSS verdict above.

Feature construction and direct-horizon models

Univariate only: lags and rolling statistics of each target series itself, no exogenous variables, for the cleanest possible apples-to-apples comparison against ARIMA:

Lags

[1, 2, 3, 4, 8] quarters: 1–4 covers the AR order Phase 2 already selected as best; 8 gives a 2-year lookback

Rolling features

4-quarter trailing mean and standard deviation

Forecast strategy

Direct, per-horizon: one model trained per forecast step (h=1..4), not recursive; no error compounding across the horizon

Model scope

Per-series (5 independent models) and one global model pooling all 5 series with a one-hot series_id feature

Hyperparameters

max_depth=3, learning_rate=0.05, subsampling, L2 regularization: shallow and regularized by design, given ~80–140 rows per fold

The bug: direct-horizon target leakage

target_h{h} = value.shift(-h) only produces NaN at the tail of the entire series, never at a fold boundary. The first implementation trained on every row in a fold’s training window, including the last h rows, whose targets land at positions i, i+1, …, i+h-1, inside that same fold’s own test window.

Caught because an early Optuna tuning run produced an implausible result: one series’ outer-fold MASE dropped from a plausible 0.18 to 0.025 after only 5 trials. Verified directly (a fold’s last training row’s target position exceeded the fold’s own cutoff), fixed with the guard above, and pinned down with a dedicated regression test so it can’t silently reappear.

MASE vs. seasonal-naive, all five series

Pooling all five series into one global model (one-hot series_id) helped both ML models on 4 of 5 accounts over their own per-series version (more effective training rows), but did not close the gap to ARIMA on any account:

Diebold-Mariano test: is ARIMA’s edge real?

The Diebold-Mariano test compares two models’ forecast-loss series and asks whether the difference is real or within sampling noise. With only ~59 folds per series (not thousands), the Harvey-Leybourne-Newbold small-sample correction was applied, using a Newey-West-style long-run variance with h-1 lags to account for the serial correlation h-step-ahead forecast errors carry.

Run per series × model × horizon (40 comparisons, 59 folds each). ARIMA is significantly more accurate (p<0.05) in 27 of 40 comparisons, and directionally favored in all 40 of 40: the non-significant cells are “not proven at 5%,” never reversed.

Series-level pattern: Business (LightGBM) and Mortgage (XGBoost) go 4-for-4 significant; Credit Card loses significance cleanly at h=3/h=4 on both models; All Loans is the weakest of the five, with LightGBM losing significance starting at h=2.

Native ARIMA intervals vs. walk-forward conformal

ARIMA gets calibrated intervals almost for free via get_forecast().conf_int(). Tree models don’t: a walk-forward-safe split-conformal method was used instead: the interval half-width for horizon h is the (1−α)-quantile of that horizon’s past absolute errors only (folds strictly before the current one), so there’s no leakage into the interval itself, on top of the point-forecast leakage guard above.

SHAP: what are the models actually looking at?

One final model per horizon, fit on all available history (there’s no future fold left to hold out for a model meant for actual deployment), explained with shap.TreeExplainer.

This is the mechanism behind the significance results above: ARIMA and the ML models converge on the same signal, but ARIMA writes it down directly in its coefficients while the tree ensembles spend a few hundred splits and ~80–140 training rows rediscovering roughly the same relationship.

Nested Optuna search: does tuning actually help here?

The fold sequence for each series was split by time order (not randomly) into the earliest ~70% (inner, 41 folds: Optuna searched here, 20 trials, minimizing mean backtest MAE) and the latest ~30% (outer, 18 folds, never touched by the search). Both fixed-default and tuned numbers below are evaluated on the same untouched outer folds.

Tuning helped in only 4 of 10 combinations, and made the outer-fold result worse in 6 of 10, consistent with the small-N overfitting risk that motivated trying fixed defaults first.

Model complexity has to match data size

For this problem (univariate, ~140-quarter macro series, direct-horizon forecasting), classical ARIMA is the stronger model, significantly so on the majority of series/horizon combinations. Global pooling and hyperparameter tuning both help the ML side somewhat but don’t close the gap.

With only ~140 observations per account, a 3-parameter model can be estimated reliably; a model built from hundreds of tree splits needs more data than that to find structure ARIMA doesn’t already capture. SHAP confirms it isn’t a training or tuning problem: both approaches are extracting the same signal, one just does it more efficiently at this sample size. The one lever that could plausibly change this result is genuinely new information (an exogenous macro regressor ARIMA never had access to either), not further feature engineering on the same short series.