Forecasting Models: From SARIMA to Foundation Models

TL;DR — There is no single best forecasting model. The landscape has four families: classical (SARIMA, ETS, Prophet — great when the signal is clean seasonality and trend), ML on lag/calendar features (Ridge, LightGBM/XGBoost — the workhorses when you have rich features like weather and holidays), deep learning (LSTM, N-BEATS — shine with lots of data and long/complex patterns), and foundation time-series models (Chronos, TimesFM — pretrained, forecast zero-shot with no training). Because each family fails differently, practitioners run a heterogeneous panel and combine them (a simple average is a shockingly strong baseline; weight by validation error to do better). Wrap every model behind one fit/predict interface and a small registry so you can swap and compare them fairly.

---

1. Simple explanation

A forecasting model takes what you know up to now (the forecast origin) and produces the next H values. But models "think" in very different ways.

Analogy — a panel of weather forecasters. Imagine four forecasters predicting tomorrow's temperature:

• The almanac reader (classical) knows the seasons cold: "it's July, it's usually warm, here's the typical curve." Superb when the pattern is regular; lost when something unusual happens.
• The data analyst (ML on features) ignores theory and reads the dials: yesterday's temp, humidity, the calendar, the satellite. Give them good instruments (features) and they win.
• The pattern-memorizer (deep learning) has seen millions of days and remembers long, subtle sequences. Brilliant with huge data; useless with a week of it.
• The well-traveled generalist (foundation model) has forecasted thousands of different places and can give a decent guess for a city they've never seen — zero-shot, no local training.

No single forecaster is best on every day. So you ask all of them and average their answers. The average is calmer and usually more accurate than any one voice — because their mistakes point in different directions and partly cancel. That is a heterogeneous panel.

We keep the running domain: hourly electricity load (MW), forecasting a 24-hour horizon.

---

2. Diagram

                     THE FORECASTING MODEL LANDSCAPE
   known history to t0 ─────────────┬───────────────► forecast t0+1 … t0+H

   ┌───────────────┐ ┌──────────────────┐ ┌────────────────┐ ┌──────────────────┐
   │  CLASSICAL    │ │  ML on FEATURES  │ │  DEEP LEARNING │ │   FOUNDATION TS  │
   │  SARIMA       │ │  Ridge (linear)  │ │  LSTM (RNN)    │ │  Chronos         │
   │  ETS          │ │  LightGBM        │ │  N-BEATS       │ │  TimesFM         │
   │  Prophet      │ │  XGBoost         │ │                │ │  (pretrained)    │
   └──────┬────────┘ └────────┬─────────┘ └───────┬────────┘ └────────┬─────────┘
          │ needs clean       │ needs good        │ needs LOTS       │ needs NOTHING
          │ seasonality       │ features          │ of data          │ (zero-shot)
          ▼                   ▼                    ▼                  ▼
        ŷ_A                 ŷ_B                  ŷ_C                ŷ_D
          └─────────────────────┴───────────┬────────┴──────────────┘
                                            ▼
                          ┌────────────────────────────────────┐
                          │   COMBINE (heterogeneous panel)    │
                          │   simple mean:  ŷ = (A+B+C+D)/4    │
                          │   weighted:     ŷ = Σ w_i ŷ_i      │
                          │   weights from validation error    │
                          └────────────────────────────────────┘
                                            ▼
                              final forecast (steadier than any one model)

   COMMON INTERFACE so they're interchangeable:
        model.fit(history) ──► model.predict(H) ──► array of length H

---

3. How it works

3.1 The landscape at a glance

3.2 Classical models

• SARIMA (Seasonal ARIMA). Models the series as autoregression (AR: depends on its own past values), differencing (I: to remove trend/seasonality and reach stationarity), and moving average (MA: depends on past errors), each with a seasonal counterpart. Notation (p,d,q)(P,D,Q)_s. For hourly load with a daily season, s = 24. Strong when seasonality is stable; painful to fit with multiple seasonalities (daily and weekly).
• ETS (Error, Trend, Seasonal — exponential smoothing). A weighted average of past values where recent points count more, with explicit components for level, trend, and seasonal. Fast, robust, few parameters. Holt-Winters is the classic seasonal form.
• Prophet (from Meta). A curve-fitting model: trend + seasonality + holidays + noise, using flexible seasonal terms (Fourier series) and easy holiday handling. Great for business series with strong calendar effects and analysts who want knobs they understand.

3.3 ML on lag/calendar features

Turn the forecasting problem into a tabular regression problem. Each row is a timestamp; the target is load_t; the features are lags (lag_1, lag_24, lag_168), rolling stats, calendar (hour, day-of-week, holiday), and weather (temperature).

• Ridge — linear regression with L2 regularization. A fast, strong, interpretable baseline. If Ridge on good features is close to your fancy model, use Ridge.
• LightGBM / XGBoost — gradient-boosted decision trees. They build many small trees, each correcting the previous one's errors. They capture non-linear effects (temperature's U-shaped effect on load) and feature interactions automatically, handle mixed feature types, and train fast. This family wins a large share of real-world tabular forecasting.

Because you control the feature matrix, this family absorbs external drivers (weather, holidays) most naturally — its biggest advantage over classical models.

3.4 Deep learning

• LSTM (Long Short-Term Memory) — a recurrent network that reads the sequence step by step and carries a "memory" cell, letting it remember patterns over long spans. Good for complex temporal dynamics and multi-output horizons, but data-hungry and fiddly.
• N-BEATS — a deep model built from stacks of fully-connected blocks that expand the forecast onto learned basis functions (trend and seasonality bases), producing both a forecast and a "backcast." Strong pure-forecasting architecture; competitive on benchmarks without hand-built features.

Deep models pay off when you have many series (thousands of meters) or long, subtle patterns and enough data to train.

3.5 Foundation time-series models

• Chronos (Amazon) — tokenizes the numeric series and runs it through a pretrained language-model-style architecture; forecasts zero-shot.
• TimesFM (Google) — a decoder-only transformer pretrained on a huge, diverse time-series corpus; forecasts new series with no training.

The pitch: like a large language model, they were pretrained on enormous, varied data, so they give a reasonable forecast for a series they have never seen — perfect for cold start or a fast, strong baseline. The catch: less local control, and they can miss idiosyncrasies like your holiday calendar or a plant-specific load spike.

3.6 Why run a heterogeneous panel — and how to combine

No family dominates across all conditions. Classical nails clean seasonality; trees exploit weather; foundation models cover cold start; deep nets capture long patterns. Running a panel and combining them is a well-established way to cut error and variance, because different models make different, partly-cancelling errors.

A practical weighting: set w_i ∝ 1 / error_i on a validation window, then normalize. An ensemble also gives you a spread across members (how much they disagree), which is a rough, textbook signal of forecast uncertainty — see Ensemble Methods for the general theory of why averaging diverse models helps. Keep it simple: the mean is your baseline; earn any added complexity.

---

4. The math

SARIMA (compact form). With backshift operator B (where B·y_t = y_{t−1}) and seasonal period s:

φ(B) · Φ(B^s) · (1−B)^d · (1−B^s)^D · y_t  =  θ(B) · Θ(B^s) · ε_t

where φ, θ are the non-seasonal AR/MA polynomials of orders p, q; Φ, Θ the seasonal ones of orders P, Q; d, D the non-seasonal/seasonal differencing; ε_t white noise. You don't compute this by hand — you fit it — but the shape shows it is AR + differencing + MA, times a seasonal copy.

ETS (Holt-Winters additive). Level ℓ, trend b, seasonal s, smoothing constants α, β, γ ∈ (0,1):

ℓ_t = α (y_t − s_{t−m}) + (1−α)(ℓ_{t−1} + b_{t−1})
b_t = β (ℓ_t − ℓ_{t−1})  + (1−β) b_{t−1}
s_t = γ (y_t − ℓ_{t−1} − b_{t−1}) + (1−γ) s_{t−m}
ŷ_{t+h} = ℓ_t + h·b_t + s_{t+h−m}      (m = season length)

Ridge. Minimize squared error plus an L2 penalty (λ controls shrinkage):

min_w  Σ_t (y_t − w·x_t)²  +  λ ‖w‖²

Ensemble combination. With M model forecasts ŷ_i and weights w_i ≥ 0, Σ w_i = 1:

ŷ_ensemble = Σ_{i=1}^{M} w_i · ŷ_i          (simple mean when all w_i = 1/M)

Worked numeric example — a 3-model panel for one hour.

Three models forecast the load (GW) at hour t0+1. The actual turns out to be 4.00.

model            forecast ŷ_i     abs error |ŷ_i − 4.00|
SARIMA  (A)        3.70               0.30
LightGBM (B)       4.20               0.20
Chronos  (C)       4.10               0.10

Simple mean:

ŷ = (3.70 + 4.20 + 4.10) / 3 = 12.00 / 3 = 4.00   → error 0.00

The mean lands exactly on target here because the members err in opposite directions and cancel — the whole point of a panel.

Weighted mean using inverse validation error. Suppose on a validation window the models had MAEs A=0.40, B=0.25, C=0.20. Weights ∝ 1/MAE:

raw:   1/0.40 = 2.50   1/0.25 = 4.00   1/0.20 = 5.00     sum = 11.50
w_A = 2.50/11.50 = 0.217   w_B = 4.00/11.50 = 0.348   w_C = 5.00/11.50 = 0.435

ŷ = 0.217·3.70 + 0.348·4.20 + 0.435·4.10
  = 0.803       + 1.462      + 1.784
  = 4.049       → error 0.049

Here the simple mean happened to win, but the weighted mean is usually steadier over many hours because it trusts the historically better models. Always check on a validation window which combiner wins — do not assume.

---

5. Real code

Every model hides behind one fit(history) / predict(H) interface, registered in a small dict so you can add, swap, and ensemble them uniformly. The LightGBM path is a full lag-feature forecaster; the LSTM and foundation-model paths are honest sketches you can flesh out.

"""A tiny forecasting REGISTRY: many models behind ONE fit/predict interface,
then a heterogeneous-panel ensemble. Domain: hourly electricity load."""
import numpy as np
import pandas as pd

# ---------- data: a realistic hourly load series ----------
rng = np.random.default_rng(0)
n = 24 * 90
idx = pd.date_range("2024-01-01", periods=n, freq="h")
load = (5.0
        + 2.0 * np.sin((idx.hour - 6) / 24 * 2 * np.pi)     # daily cycle
        + np.where(idx.dayofweek >= 5, -0.8, 0.0)           # weekend dip
        + np.linspace(0, 0.5, n)                            # trend
        + rng.normal(0, 0.15, n))
series = pd.Series(load, index=idx, name="load")
H = 24


def make_features(s: pd.Series) -> pd.DataFrame:
    """Time-safe lag + calendar features (shift before rolling — no peeking)."""
    df = pd.DataFrame({"load": s})
    df["lag_1"]   = df["load"].shift(1)
    df["lag_24"]  = df["load"].shift(24)
    df["lag_168"] = df["load"].shift(168)
    df["roll24"]  = df["load"].shift(1).rolling(24).mean()
    df["hour"]    = df.index.hour
    df["dow"]     = df.index.dayofweek
    return df


# ---------- common interface: every model implements fit(history), predict(H) ----------
class SeasonalNaive:
    """Baseline: tomorrow's hour = value one season (24h) ago."""
    def fit(self, s):     self.s = s; return self
    def predict(self, H): return self.s.iloc[-24:].values[:H] if H <= 24 else \
                                 np.resize(self.s.iloc[-24:].values, H)


class LGBMForecaster:
    """ML-on-features model: gradient-boosted trees on lag/calendar features.
    Recursive multi-step: predict one hour, append it, re-featurize, repeat."""
    def __init__(self):
        from lightgbm import LGBMRegressor
        self.model = LGBMRegressor(n_estimators=300, learning_rate=0.05,
                                   num_leaves=31, min_child_samples=20, verbose=-1)

    def fit(self, s):
        self.s = s.copy()
        df = make_features(s).dropna()
        self.cols = ["lag_1", "lag_24", "lag_168", "roll24", "hour", "dow"]
        self.model.fit(df[self.cols], df["load"])
        return self

    def predict(self, H):
        hist = self.s.copy()
        preds = []
        for _ in range(H):
            feats = make_features(hist).iloc[[-1]][self.cols]   # newest row
            yhat = float(self.model.predict(feats)[0])
            preds.append(yhat)
            nxt = hist.index[-1] + pd.Timedelta(hours=1)
            hist.loc[nxt] = yhat                                # feed prediction back
        return np.array(preds)


class LSTMForecaster:
    """Deep-learning SKETCH (PyTorch). Windows of past load -> next value, multi-output.
    Shown structurally; train longer with more data for real use."""
    def fit(self, s, lookback=48, epochs=3):
        import torch, torch.nn as nn
        self.torch, self.nn, self.lookback = torch, nn, lookback
        v = ((s - s.mean()) / s.std()).values.astype("float32")
        self.mu, self.sd = float(s.mean()), float(s.std())
        X = np.stack([v[i:i+lookback] for i in range(len(v) - lookback - 1)])
        y = v[lookback+1: len(v)]
        Xt = torch.tensor(X).unsqueeze(-1); yt = torch.tensor(y).unsqueeze(-1)
        self.net = nn.Sequential(nn.LSTM(1, 32, batch_first=True))   # (see forward below)
        self.head = nn.Linear(32, 1)
        opt = torch.optim.Adam(list(self.net.parameters()) + list(self.head.parameters()), 1e-2)
        loss_fn = nn.MSELoss()
        for _ in range(epochs):
            out, _ = self.net(Xt); pred = self.head(out[:, -1, :])
            loss = loss_fn(pred, yt); opt.zero_grad(); loss.backward(); opt.step()
        self.tail = v[-lookback:]
        return self

    def predict(self, H):
        torch = self.torch
        window = list(self.tail); preds = []
        for _ in range(H):
            x = torch.tensor(np.array(window[-self.lookback:], "float32")).view(1, -1, 1)
            out, _ = self.net(x); yhat = float(self.head(out[:, -1, :]).item())
            preds.append(yhat * self.sd + self.mu); window.append(yhat)   # recursive
        return np.array(preds)


class FoundationForecaster:
    """Foundation TS model (e.g. Chronos/TimesFM): ZERO-SHOT, no training on your data.
    Pattern shown; requires the respective package + weights installed."""
    def fit(self, s):
        self.context = s.values.astype("float32"); return self   # no training!

    def predict(self, H):
        # --- with Chronos (illustrative) ---
        # from chronos import ChronosPipeline; import torch
        # pipe = ChronosPipeline.from_pretrained("amazon/chronos-t5-small")
        # fc = pipe.predict(torch.tensor(self.context), prediction_length=H)  # samples
        # return fc.median(dim=1).values.numpy().ravel()
        # fallback so this file runs without the weights: seasonal-naive echo
        return np.resize(self.context[-24:], H)


# ---------- the REGISTRY ----------
REGISTRY = {
    "seasonal_naive": SeasonalNaive,
    "lightgbm":       LGBMForecaster,
    "lstm":           LSTMForecaster,
    "foundation":     FoundationForecaster,
}


def run_panel(series, H, names):
    """Fit each named model on the same history, collect H-step forecasts."""
    out = {}
    for name in names:
        try:
            out[name] = REGISTRY[name]().fit(series).predict(H)
        except Exception as e:          # a missing lib shouldn't kill the panel
            print(f"skip {name}: {e}")
    return out


def ensemble(forecasts: dict, weights: dict | None = None) -> np.ndarray:
    names = list(forecasts)
    M = np.vstack([forecasts[k] for k in names])
    if weights is None:
        return M.mean(axis=0)                       # simple mean
    w = np.array([weights[k] for k in names]); w = w / w.sum()
    return (w[:, None] * M).sum(axis=0)             # weighted mean


if __name__ == "__main__":
    train = series.iloc[:-H]
    panel = run_panel(train, H, ["seasonal_naive", "lightgbm", "foundation"])
    combo = ensemble(panel)                          # simple average
    actual = series.iloc[-H:].values
    for name, fc in panel.items():
        print(f"{name:14s} MAE = {np.mean(np.abs(fc - actual)):.3f}")
    print(f"{'ENSEMBLE':14s} MAE = {np.mean(np.abs(combo - actual)):.3f}")

Typical output (exact numbers depend on installed libraries and seed):

seasonal_naive MAE = 0.2xx
lightgbm       MAE = 0.1xx
foundation     MAE = 0.2xx
ENSEMBLE       MAE = 0.1xx      # usually at or below the best single model

---

6. Real-world example

Building a load-forecasting service for a utility, model by model.

• Cold start (week 1). A new region just came online — you have days of data, not years. You cannot train an LSTM. You deploy a foundation model (Chronos/TimesFM) zero-shot plus seasonal-naïve, and ship a usable day-ahead forecast immediately.
• A few months in. Enough history to engineer features. You add a LightGBM model on lag_24, lag_168, rolling means, hour, dow, holiday, and temperature. It beats naïve most on hot days because it learns temperature's non-linear effect. You add SARIMA(2,0,1)(1,1,1)_24 as a clean-seasonality reference and a Ridge on the same features as an interpretable floor.
• Year one, lots of meters. With thousands of related series you train an N-BEATS/LSTM across all of them (a "global" model), sharing patterns between similar substations.
• Production. You run all of them as a heterogeneous panel behind the same fit/predict interface and combine with a validation-error-weighted mean. When LightGBM stumbles on an unusual holiday, SARIMA and the foundation model keep the ensemble sane. The ensemble MAPE is lower and, crucially, less volatile week to week than any single model.
• The payoff of the interface. Because every model implements fit/predict, swapping a model, adding a new one, or dropping one for an ablation is a one-line registry change — not a rewrite. That is what makes fair comparison and Forecast Evaluation & Backtesting tractable.

Concretely: on the evening peak (hour 19), seasonal-naïve MAPE ≈ 4.0%, LightGBM ≈ 2.8%, the panel ≈ 2.5% — and on the two summer heat-wave days that broke naïve entirely, the panel held near 3% while naïve blew past 9%.

---

7. Interview questions companies actually ask

Q [Amazon Forecast / SageMaker] "When would you choose SARIMA over gradient-boosted trees?"
  A SARIMA when the series is driven mainly by its own clean, stable seasonality and trend, you
    have limited history, and you want an interpretable, statistically-grounded model with
    natural intervals. Trees (LightGBM/XGBoost) when you have rich EXTERNAL drivers — weather,
    holidays, promotions — because you control the feature matrix and trees capture their
    non-linear effects and interactions automatically.

Q [Google / TimesFM team] "What is a foundation time-series model and when is zero-shot useful?"
  A A model pretrained on a huge, diverse corpus of series (Chronos, TimesFM) that forecasts a
    NEW series with no task-specific training. Zero-shot shines at cold start (little/no local
    history), for a fast strong baseline, or across many heterogeneous series. Trade-off: less
    local control — it may miss your specific holidays or plant-specific spikes.

Q [Meta / Prophet] "Why do practitioners ensemble forecasting models instead of picking one?"
  A Different families make different, partly-cancelling errors: classical nails clean
    seasonality, trees exploit weather, foundation models cover cold start, deep nets capture
    long patterns. Averaging diverse models reduces both error and variance. A simple mean is a
    very strong baseline; weighting by validation error usually improves it. The ensemble is
    also steadier over time than any single model.

Q [Uber / DoorDash] "How do you produce a multi-step forecast from a one-step tree model?"
  A Recursively: predict t+1, append the prediction to the history, rebuild the lag features,
    predict t+2, and so on — but beware error compounding. Alternatively train direct models per
    horizon (one 'predict t+3' model), which avoids feedback at the cost of maintaining H models.

Q [Netflix] "What's the simplest model that could work, and why start there?"
  A Seasonal-naïve (copy the value one season back) plus Ridge on lag/calendar features. They're
    fast, interpretable, and on strongly seasonal load they're surprisingly hard to beat. They
    set the bar every fancier model must clear — if LightGBM barely beats Ridge, ship Ridge.

Q [Two Sigma / quant] "Explain the components of SARIMA."
  A AR (depends on its own past values, order p), I (differencing d times to reach
    stationarity), MA (depends on past forecast errors, order q), each with a seasonal
    counterpart (P, D, Q) at period s. Written (p,d,q)(P,D,Q)_s. For hourly load with a daily
    cycle, s = 24; the seasonal difference removes the daily pattern.

Q [Microsoft / energy] "LSTM vs N-BEATS — when and why?"
  A Both are deep and data-hungry. LSTM reads the sequence step-by-step with a memory cell —
    flexible for complex dynamics and multi-output horizons. N-BEATS uses stacked
    fully-connected blocks projecting onto learned trend/seasonality bases, giving strong
    pure-forecasting accuracy without hand-built features. Use them when you have lots of data
    (ideally many related series); otherwise trees usually win with less fuss.

Q [Anany DS role] "How would you weight models in an ensemble?"
  A Score each model on a held-out VALIDATION window (walk-forward, no leakage), then set weights
    inversely proportional to its error and normalize so they sum to 1. Compare the weighted mean
    against the simple mean on validation — sometimes the simple mean wins, so don't add
    complexity unless it pays off. Consider a trimmed mean/median if one model occasionally
    blows up.

Q [Stripe] "Why hide every model behind the same fit/predict interface?"
  A So models are interchangeable: you can add, swap, drop, or ensemble them without rewriting
    the pipeline, and you can backtest them under identical conditions for a fair comparison. A
    small registry (name -> class) makes ablations a one-line change and keeps evaluation honest.

---

8. When to use / tradeoffs

   PICK BY SITUATION:
     ✓ clean seasonality, short history, need intervals + interpretability → SARIMA / ETS
     ✓ strong calendar/holiday effects, business series, analyst knobs       → Prophet
     ✓ rich features (weather, holidays), moderate data, tabular tooling     → LightGBM/XGBoost
     ✓ fast interpretable floor / sanity baseline                            → Ridge
     ✓ lots of data, many related series, long patterns                      → LSTM / N-BEATS
     ✓ cold start, no time to train, strong instant baseline                 → Chronos / TimesFM
   HONEST LIMITS:
     ✗ classical struggles with many external drivers and multiple seasonalities
     ✗ trees have no native 'sequence' notion — you MUST engineer lags/rolling features
     ✗ deep models are data-hungry, slow, and easy to overfit on one series
     ✗ foundation models are less controllable and may miss local quirks (your holidays)
     ✗ ensembles cost more compute and can hide WHY a forecast moved
   THE PANEL REMINDER:
     no family wins everywhere. Run a heterogeneous panel behind one interface, combine with a
     simple mean first, and only add weighting/complexity when validation says it pays.

---

9. Summary + related articles

• Four model families: classical (SARIMA, ETS, Prophet), ML on features (Ridge, LightGBM, XGBoost), deep learning (LSTM, N-BEATS), foundation TS (Chronos, TimesFM — zero-shot).
• Classical loves clean seasonality; trees love engineered features (weather, holidays); deep loves big data; foundation loves cold start.
• No family wins everywhere, so run a heterogeneous panel and combine — the simple mean is a strong default; weight by validation error to improve.
• An ensemble also gives a spread (member disagreement) as a rough uncertainty signal — kept textbook-general here.
• Wrap everything in one fit/predict interface and a small registry so models are interchangeable and fairly comparable.
• Always keep seasonal-naïve in the panel as the baseline to beat.

Related: Time-Series Forecasting Fundamentals · Forecast Evaluation & Backtesting · Ensemble Methods

Resources

• Hyndman & Athanasopoulos, Forecasting: Principles and Practice (ARIMA, ETS, Prophet chapters) — https://otexts.com/fpp3/
• statsmodels SARIMAX documentation — https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html
• Meta Prophet — https://facebook.github.io/prophet/
• LightGBM documentation — https://lightgbm.readthedocs.io/
• Oreshkin et al., "N-BEATS" — https://arxiv.org/abs/1905.10437
• Ansari et al., "Chronos: Learning the Language of Time Series" — https://arxiv.org/abs/2403.07815
• Das et al., "A decoder-only foundation model for time-series forecasting" (TimesFM) — https://arxiv.org/abs/2310.10688