Forecasting Methods¶

ForecastAI supports 10 forecasting algorithms across three families. This chapter provides full mathematical descriptions, implementation details, parameter tables, and supply chain guidance for each.

Method Families at a Glance¶

mindmap
  root((Forecasting\nMethods))
    Statistical
      AutoARIMA
      AutoETS
      AutoTheta
      MSTL
      CrostonOptimized
    Neural
      NHITS
      NBEATS
      PatchTST
      TFT
      DeepAR
    Foundation
      TimesFM

Method	Library	Probabilistic	Training	Best for
AutoARIMA	StatsForecast	Yes (analytical)	Per-series	Stationary / trended, no strong seasonality
AutoETS	StatsForecast	Yes (analytical)	Per-series	Trended series with additive/multiplicative seasonality
AutoTheta	StatsForecast	Yes (conformal)	Per-series	General demand with trend
MSTL	StatsForecast	Yes (via AutoARIMA)	Per-series	Multiple overlapping seasonal periods
CrostonOptimized	StatsForecast	Yes (conformal)	Per-series	Intermittent demand (ADI > 1.32)
NHITS	NeuralForecast	Yes (MQLoss)	Global	Long horizons, multiple seasonal patterns
NBEATS	NeuralForecast	Yes (MQLoss)	Global	Large portfolios, no domain knowledge needed
PatchTST	NeuralForecast	Yes (MQLoss)	Global	Long context windows, transformer
TFT	NeuralForecast	Yes (MQLoss)	Global	Rich covariate datasets
DeepAR	NeuralForecast	Yes (parametric)	Global	Large portfolios, similar patterns
TimesFM	timesfm	Yes (scaling)	Pre-trained	Cold-start, new products

Statistical Models (StatsForecast)¶

Statistical models are fast, interpretable, and trained per-series. They produce native prediction intervals from analytical theory or conformal calibration. Implementation: files/forecasting/statistical_models.py, class StatisticalForecaster.

Common Pre-Flight Guards¶

Before any statistical model is fit, minimum data requirements are checked:

Model	Minimum observations
SeasonalNaive, SeasonalWindowAverage	\(\geq\) `season_length`
HistoricAverage	\(\geq 2\)
MSTL	\(\geq 2 \times \text{season\_length} + 1\)
AutoETS, AutoARIMA, AutoTheta, AutoCES	\(\geq \max(2 \times \text{season\_length},\; 10)\)
CrostonOptimized, ADIDA, IMAPA	\(\geq 2\)
`Naive`	\(\geq 1\) — last-resort forecasting-time fallback when every recommended method is ineligible

If the check fails, the method is skipped for that series with a logged warning.

Season length resolution: The configured season_length is derived from the first detected seasonal period in the series characteristics. If no seasonal period was detected, it falls back to a frequency map: Weekly → 52, Monthly → 12, Quarterly → 4, etc.

Frequency alias mapping: Pandas 2.x deprecated some aliases ('M' → 'ME', 'Q' → 'QE', 'Y' → 'YE'). The StatisticalForecaster applies this mapping before passing to StatsForecast.

AutoARIMA¶

What it is: ARIMA (AutoRegressive Integrated Moving Average) is the foundational model for univariate time series. The "Auto" variant uses the Hyndman-Khandakar algorithm to automatically search for the best ARIMA orders using information criteria.

Model: The seasonal ARIMA model is written \(\text{ARIMA}(p,d,q)(P,D,Q)_s\):

\[\phi_p(B)\Phi_P(B^s) \nabla^d \nabla_s^D y_t = c + \theta_q(B)\Theta_Q(B^s) \epsilon_t\]

where: - \(B\) is the backshift operator: \(B y_t = y_{t-1}\) - \(\nabla^d = (1 - B)^d\) is the \(d\)-th order difference operator - \(\nabla_s^D = (1 - B^s)^D\) is the \(D\)-th order seasonal difference - \(\phi_p(B) = 1 - \phi_1 B - \cdots - \phi_p B^p\) — non-seasonal AR polynomial - \(\Phi_P(B^s) = 1 - \Phi_1 B^s - \cdots - \Phi_P B^{Ps}\) — seasonal AR polynomial - \(\theta_q(B) = 1 + \theta_1 B + \cdots + \theta_q B^q\) — non-seasonal MA polynomial - \(\Theta_Q(B^s)\) — seasonal MA polynomial - \(\epsilon_t\) — white noise with variance \(\sigma^2\) - \(s\) — seasonal period (52 for weekly data)

Components explained:

Component	Symbol	Effect
Autoregression	\(p\)	Current value depends on \(p\) past values
Differencing	\(d\)	Removes \(d\)-order polynomial trend (d=1 removes linear trend)
Moving average	\(q\)	Current value depends on \(q\) past forecast errors
Seasonal AR	\(P\)	Seasonal AR at multiples of \(s\)
Seasonal differencing	\(D\)	Removes seasonal trend (D=1 removes one seasonal difference)
Seasonal MA	\(Q\)	Seasonal MA at multiples of \(s\)

Auto-selection: The Hyndman-Khandakar (2008) algorithm:

Determine \(d\) and \(D\) using KPSS and Canova-Hansen tests
Run a stepwise search over \((p, q, P, Q)\) starting from \((2, 2, 1, 1)\)
Select the model minimizing AICc (corrected Akaike Information Criterion):

\[\text{AICc} = -2\ln(\hat{L}) + 2k + \frac{2k(k+1)}{n - k - 1}\]

where \(k\) is the number of parameters and \(\hat{L}\) is the maximized likelihood.

Prediction intervals: Analytically derived from the MA representation of the model. The \(h\)-step forecast variance is:

\[V(e_h) = \sigma^2 \sum_{j=0}^{h-1} \psi_j^2\]

where \(\psi_j\) are the infinite MA coefficients obtained from the model.

Parameters in ForecastAI:

Parameter	Default / Source	Description
`season_length`	From characteristics	Seasonal period \(s\)
`approximation`	`True` if n > 150	Speed approximation for long series
`p, d, q, P, D, Q`	Auto-selected by AICc	ARIMA orders

When manual orders are provided via the hyperparameter override system, the model switches from AutoARIMA() to a fixed ARIMA() with exact specified orders.

Overridable parameters: p, d, q, P, D, Q, max_p, max_q, max_P, max_Q, max_order, max_d, max_D, start_p, start_q, start_P, start_Q, stationary, seasonal, stepwise, allowdrift, allowmean.

Strengths: - Analytically rigorous prediction intervals - Handles trended and integrated series naturally via differencing - Interpretable model orders - Very fast fitting

Weaknesses: - Assumes linear relationships - A single seasonal period only (use MSTL for multiple) - Can over-fit on short series with high orders

When selected: standard group (default path when no seasonality, high complexity, or intermittency is detected). Also in seasonal group as a secondary candidate.

AutoETS (Exponential Smoothing State Space Model)¶

What it is: ETS stands for Error-Trend-Seasonal. The state space framework unifies all exponential smoothing methods into a single probabilistic model. The "Auto" variant selects the best combination of components by minimizing AICc.

Model: The general ETS model is specified by three components, each of which can be additive (A), multiplicative (M), none (N), or for trend: damped (D):

Measurement equation (additive error, additive trend, additive seasonal — ETS(A,A,A)):

\[y_t = l_{t-1} + b_{t-1} + s_{t-m} + \epsilon_t\]

State equations:

\[l_t = l_{t-1} + b_{t-1} + \alpha \epsilon_t \quad \text{(level)}\]

\[b_t = b_{t-1} + \beta \epsilon_t \quad \text{(trend)}\]

\[s_t = s_{t-m} + \gamma \epsilon_t \quad \text{(seasonal)}\]

\(h\)-step forecast (point):

\[\hat{y}_{t+h|t} = l_t + h \cdot b_t + s_{t+h-m(k+1)}\]

where \(k = \lfloor (h-1)/m \rfloor\) and \(m\) is the season length.

Smoothing parameters (\(\alpha\), \(\beta\), \(\gamma\)) control how quickly each component adapts: - \(\alpha \in (0,1)\): Level smoothing. High = reactive, low = smooth. - \(\beta \in (0,\alpha)\): Trend smoothing. - \(\gamma \in (0,1-\alpha)\): Seasonal smoothing.

Damped trend (model type 'A' for damped): A damping parameter \(\phi \in (0,1)\) prevents explosive trend extrapolation:

\[\hat{y}_{t+h|t} = l_t + (\phi + \phi^2 + \cdots + \phi^h) b_t + s_{t+h-m(k+1)}\]

As \(h \to \infty\), the trend contribution converges to a constant.

Auto-selection: Tries up to 30 model combinations (Error \(\times\) Trend \(\times\) Seasonal from the allowed set) and selects by AICc. Some combinations are excluded (multiplicative error with damped trend, etc.) to avoid non-identifiable models.

Parameters in ForecastAI:

Parameter	Default	Description
`model`	`'ZZZ'`	Auto-select all three components
`season_length`	From characteristics	Seasonal period \(m\)

Extracted fitted parameters stored for review: alpha, beta, gamma, phi, sigma2, errortype, trendtype, seasontype, damped, aic, aicc, bic.

Overridable: model (e.g., 'AAA', 'MNM'), damped, phi.

Strengths: - Handles additive and multiplicative seasonality - Damped trend prevents over-extrapolation on short horizons - Native prediction intervals from state space theory - Very robust across diverse demand patterns

Weaknesses: - Single seasonal period only - Multiplicative models can produce negative forecasts for very small series — guarded by corrected_qty > 0 checks

When selected: All groups except intermittent and sparse. The most commonly selected method in the Bicycle dataset (222 out of 587 series).

AutoTheta¶

What it is: The Theta method (Assimakopoulos & Nikolopoulos, 2000) decomposes a series into two "theta lines" by modifying the local curvature. The optimized version (AutoTheta) automatically selects the best variant.

Theta decomposition: Starting from a series \(y_t\), a theta line \(Z_t(\theta)\) is defined by the second-order differential equation modification:

\[\frac{d^2 Z_t}{dt^2} = \theta \frac{d^2 y_t}{dt^2}\]

Two special theta lines:

\(\theta_1 = 0\): Extrapolates only the long-term linear trend (amplifies \(d^2y/dt^2 = 0\))
\(\theta_2 = 2\): Doubles the local curvature, capturing short-term dynamics

Combination: The final forecast is the mean of the two theta lines:

\[\hat{y}_t = \frac{1}{2}Z_t(\theta_1 = 0) + \frac{1}{2}Z_t(\theta_2)\]

In practice: \(Z_t(\theta_1)\) is a simple linear regression fitted to the data (extrapolated trend), and \(Z_t(\theta_2)\) is a Simple Exponential Smoothing (SES) forecast (captures level).

AutoTheta variants: The library implements DynamicTheta, DynamicOptimizedTheta, and OptimizedTheta. These extend the basic method by allowing \(\theta\) to be optimized over time rather than fixed.

Decomposition type — controlled by the characteristics:

Complexity level	`decomposition_type`
High	`'multiplicative'`
Low / Medium	`'additive'`

Parameters in ForecastAI:

Parameter	Default	Description
`season_length`	From characteristics	Seasonal period
`decomposition_type`	Based on complexity	How seasonal component is modelled

Extracted fitted parameters: theta, alpha, drift.

Overridable: model ('OptimizedTheta', 'DynamicTheta', 'DynamicOptimizedTheta').

Strengths: - Simple and interpretable - Robust to over-fitting - Good for medium-length series - Runner-up to AutoETS on the Bicycle dataset (174 series)

Weaknesses: - Less expressive than full ARIMA or ETS for complex patterns - Does not model multiplicative error natively

When selected: standard group alongside AutoETS and AutoARIMA.

MSTL (Multiple Seasonal-Trend Decomposition using LOESS)¶

What it is: MSTL extends classic STL decomposition to handle multiple seasonal periods simultaneously. This is essential for weekly data that exhibits both a weekly pattern (period=52 for annual cycle) and, say, a quarterly pattern (period=13).

Figure: MSTL seasonal-decomposition flow — the trend is forecast by non-seasonal AutoARIMA; seasonal components are projected by repeating the last cycle.

flowchart LR
    Y["y_t raw series"] --> STL["MSTL decomposition<br/>LOESS, robust"]
    STL --> T["T_t trend"]
    STL --> SK["S_t^(k) seasonal<br/>k = 1..K"]
    STL --> R["R_t remainder"]
    T --> FA["AutoARIMA(season_length=1)<br/>trend forecast"]
    SK --> REP["Repeat last complete cycle"]
    FA --> COMB["Sum components"]
    REP --> COMB
    COMB --> OUT["y_hat_(t+h)"]

Decomposition model:

\[y_t = T_t + \sum_{k=1}^{K} S_t^{(k)} + R_t\]

where: - \(T_t\) — trend component (LOESS smoother over the de-seasonalized series) - \(S_t^{(k)}\) — \(k\)-th seasonal component with period \(m_k\) - \(R_t\) — remainder (residual) - \(K\) — number of seasonal periods

LOESS (Locally Estimated Scatterplot Smoothing): Each seasonal component is estimated by a locally weighted regression. For each point \(t\), a polynomial is fit to a neighbourhood of nearby points, weighted by their distance from \(t\):

\[\hat{S}_t^{(k)} = \text{LOESS}(y_t - T_t - \sum_{j \neq k} S_t^{(j)})\]

The algorithm iterates between estimating \(T_t\) and each \(S_t^{(k)}\) until convergence.

Trend forecasting: After decomposition, the trend component \(T_t\) is extrapolated using an AutoARIMA(season_length=1) (non-seasonal ARIMA). The seasonal components are projected by repeating the last complete cycle. The remainder is assumed zero.

Final forecast:

\[\hat{y}_{t+h} = \hat{T}_{t+h} + \sum_{k=1}^{K} S_{t+h \bmod m_k}^{(k)}\]

Minimum data: \(\geq 2 \times \text{season\_length} + 1\) observations.

Parameters in ForecastAI:

Parameter	Default	Description
`season_length`	From characteristics (first detected)	Primary seasonal period
`trend_forecaster`	`AutoARIMA(season_length=1)`	Fixed — non-seasonal ARIMA for trend

Dependency: Requires the supersmoother Python package (pip install supersmoother). The pipeline checks for this at import time and logs a warning if unavailable.

Strengths: - Only method that can model multiple overlapping seasonal periods - Robust decomposition (LOESS is robust to outliers) - Interpretable: trend and seasonal components are separately visualizable

Weaknesses: - Computationally more expensive than ETS/ARIMA - Requires more data (2 full cycles minimum) - Seasonal components assumed stable — does not handle evolving seasonality well

When selected: seasonal group as the primary candidate. Also explicitly checked for compatibility in check_method_compatibility (warning if has_seasonality = False).

CrostonOptimized (and ADIDA / IMAPA)¶

What it is: Croston's method (1972) is the foundational technique for intermittent demand forecasting. It separately smooths the demand sizes and the inter-demand intervals using independent exponential smoothing, then combines them.

Figure: Croston / TSB intermittent-demand flow — separate smoothing of sizes and intervals, SBA bias correction, optional conformal intervals.

flowchart TD
    D["Demand series<br/>with zero periods"] --> SEP["Separate into<br/>demand sizes z_t & intervals p_t"]
    SEP --> SZ["SES on sizes<br/>z_hat = alpha*z + (1-alpha)*z_hat"]
    SEP --> IT["SES on intervals<br/>p_hat = alpha*p + (1-alpha)*p_hat"]
    SZ --> YC["y_hat = z_hat / p_hat"]
    IT --> YC
    YC --> SBA["SBA bias correction<br/>y_SBA = (1 - alpha/2) * y_hat"]
    SBA --> OUT["CrostonOptimized forecast<br/>(point)"]
    YC --> CONF{"n >= 2h+1?"}
    CONF -- yes --> CI["Conformal intervals<br/>n_windows=2, h=horizon"]
    CONF -- no --> PO["Point forecast only"]
    CI --> OUT
    PO --> OUT

The Croston model:

Let \(z_t\) be the demand size at occurrence \(t\) (non-zero only), and \(p_t\) be the inter-demand interval (number of periods since last demand).

Exponential smoothing for sizes:

\[\hat{z}_{t+1} = \alpha z_t + (1-\alpha)\hat{z}_t\]

Exponential smoothing for intervals:

\[\hat{p}_{t+1} = \alpha p_t + (1-\alpha)\hat{p}_t\]

Forecast demand per period:

\[\hat{y} = \frac{\hat{z}}{\hat{p}}\]

Syntetos-Boylan bias correction (the "Optimized" in CrostonOptimized): The original Croston method is positively biased. Syntetos & Boylan (2001) showed the corrected estimate is:

\[\hat{y}_{\text{SBA}} = \left(1 - \frac{\alpha}{2}\right) \frac{\hat{z}}{\hat{p}}\]

The CrostonOptimized model automatically applies this bias correction and optimizes \(\alpha\) by minimizing MSE.

Prediction intervals: When \(n_{\text{obs}} \geq 2 \times \text{horizon} + 1\), conformal prediction intervals are generated using ConformalIntervals(n_windows=2, h=horizon). Without sufficient data, only point forecasts are produced.

Related intermittent methods:

ADIDA (Aggregate-Disaggregate Intermittent Demand Approach): Aggregates the series to a lower frequency to reduce intermittency, forecasts at the aggregate level, then dis-aggregates back. Handles both demand occurrence and size uncertainty.

IMAPA (Intermittent Multiple Aggregation Prediction Algorithm): Forecasts at multiple aggregation levels simultaneously, then combines the predictions. Generally more robust than single-level approaches.

TSB (Teunter-Syntetos-Babai): An alternative to Croston that models demand probability directly, with separate smoothing for demand occurrence probability (\(\alpha_d\), default 0.1) and demand size (\(\alpha_p\), default 0.1). Better tracks obsolescence risk (declining demand probability).

Parameters in ForecastAI:

Parameter	Default	Description
`alpha`	Auto-optimized	Smoothing parameter
Conformal intervals	Enabled if n ≥ 2h+1	n_windows=2

Strengths: - Purpose-built for intermittent demand — no other method handles it as gracefully - Separate smoothing for sizes and intervals allows each to have its own memory - Bias correction prevents systematic over-stocking

Weaknesses: - Does not model seasonality (assumes stationary intermittent demand) - Flat forecast — does not capture trend or seasonality in demand - Conformal intervals may be wide on very sparse series

When selected: intermittent group. Triggered when is_intermittent = True (any of: ADI > 1.32 AND CoV > 0.49, zero_ratio > 0.5, or fewer than 5 non-zero periods).

Neural Models (NeuralForecast)¶

Neural models are trained globally across all series in a pipeline, sharing learned representations. They produce probabilistic forecasts via multi-quantile regression (MQLoss). Implementation: files/forecasting/neural_models.py, class NeuralForecaster.

Common Neural Model Parameters¶

All neural models share these base settings:

Parameter	Low/Medium complexity	High complexity
`h`	`horizon` (52 weeks)	Same
`loss`	`MQLoss(level=quantiles)`	Same
`max_steps`	500 (n≤200) / 1000 (n>200)	Same
`val_check_steps`	100	Same
`early_stop_patience_steps`	3	Same

MQLoss: Multi-Quantile Loss trains the model to simultaneously produce forecasts at all configured quantile levels \([0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]\):

\[\mathcal{L}_{MQ} = \frac{1}{|Q||H|} \sum_{q \in Q} \sum_{h=1}^{H} \text{QL}(y_{t+h}, \hat{y}_{t+h}^{(q)}, q)\]

where the pinball (quantile) loss is:

\[\text{QL}(y, \hat{y}, q) = \max\left(q(y - \hat{y}),\; (q-1)(y - \hat{y})\right)\]

Minimum data: 50 observations (hardcoded); additionally filtered by sufficient_for_deep_learning.

NHITS (Neural Hierarchical Interpolation for Time Series)¶

What it is: NHITS (Challu et al., 2022) is a purely MLP-based architecture that uses multi-rate input sampling and hierarchical signal interpolation to efficiently capture patterns at multiple time scales without attention mechanisms.

Architecture:

flowchart LR
    X[Input window\n5×horizon] --> S1[Stack 1\nPooling k=1\nLarge scale]
    S1 --> I1[Interpolation\nLong-term]
    X --> S2[Stack 2\nPooling k=2\nMedium scale]
    S2 --> I2[Interpolation\nMedium-term]
    X --> S3[Stack 3\nPooling k=1\nShort scale]
    S3 --> I3[Interpolation\nShort-term]
    I1 & I2 & I3 --> OUT[Forecast\n+ Quantiles]

Each stack consists of: 1. MaxPooling with kernel size \(k\) — downsamples the input, capturing patterns at scale \(k\) 2. MLP blocks — dense layers with ReLU activations 3. Basis expansion — projects to forecast horizon via linear combination of basis functions 4. Hierarchical interpolation — upsamples the pooled output back to the full horizon

Multi-rate sampling: Different stacks receive the same full input but pool at different rates, allowing each stack to specialize at different temporal frequencies.

Double residual learning: Each stack receives both the original input and the residual from previous stacks, allowing deep stacks to focus on what earlier stacks missed.

Parameters in ForecastAI:

Parameter	Low/Medium complexity	High complexity
`input_size`	`min(5 × horizon, n // 2)`	Same
`n_blocks`	`[1, 1, 1]`	`[3, 3, 3]`
`mlp_units`	`[[256, 256]] × 3`	`[[512, 512]] × 3`
`n_pool_kernel_size`	`[1, 1, 1]` (non-seasonal)	`[2, 2, 1]` (seasonal)
`learning_rate`	1e-3	1e-3
`batch_size`	32	32

For seasonal series, the pool kernel sizes [2, 2, 1] focus two stacks on coarser scales where seasonal patterns manifest.

Strengths: - No attention mechanism → much faster than Transformers at long horizons - Multi-resolution: naturally handles weekly + annual patterns simultaneously - Strong empirical performance on M4 and M5 benchmarks - Interpretable: each stack's contribution can be visualized

Weaknesses: - No covariate support (unlike TFT) - Requires moderate amounts of data (50+ observations) - Global training: performance depends on portfolio having similar patterns

When selected: complex, seasonal, and standard groups. One of the most broadly applicable neural methods.

NBEATS (Neural Basis Expansion Analysis for Interpretable Time Series Forecasting)¶

What it is: N-BEATS (Oreshkin et al., 2019) uses a deep network with basis expansion to decompose forecasts into interpretable trend and seasonality components. It was the winning method of the M4 competition.

Architecture: A stack of blocks, each performing double residual learning:

\[\hat{y}_{t+1:t+H} = \sum_{k=1}^{K} \hat{y}_{t+1:t+H}^{(k)}\]

Each block \(k\):

Fully-connected network: \([x_{t-L:t}; r_{t-L:t}^{(k-1)}] \to \theta_{f}^{(k)},\; \theta_{b}^{(k)}\)
Basis expansion — forecast: \(\hat{y}^{(k)} = V_f \theta_f^{(k)}\), backcast: \(\hat{x}^{(k)} = V_b \theta_b^{(k)}\)
Residual: \(r^{(k)} = x^{(k-1)} - \hat{x}^{(k)}\) (the part of the input not explained by this block)

Basis functions \(V\):

Generic stack (identity): \(V = I\) — unconstrained basis, learns any shape
Trend stack: \(V_{ij} = i^j / H^j\) — polynomial basis of degree \(p\) (forces smooth trends)
Seasonality stack: \(V_{ij} = \cos(2\pi i j / m)\) or \(\sin(2\pi i j / m)\) — Fourier basis

Parameters in ForecastAI:

Parameter	Low/Medium complexity	High complexity
`input_size`	`min(5 × horizon, n // 2)`	Same
`n_blocks`	`[2, 2, 2]`	`[3, 3, 3]`
`stack_types`	`['identity', 'trend']` (non-seasonal)	`['identity', 'trend', 'seasonality']` (seasonal)
`learning_rate`	1e-3	1e-3
`batch_size`	32	32

For seasonal series, a seasonality stack using Fourier basis is added.

Strengths: - Interpretable: trend and seasonal contributions separated in the stack architecture - Proven on large benchmark datasets (won M4) - No external covariates needed - Double residual learning prevents gradient vanishing

Weaknesses: - No covariate support - Seasonality stack requires sufficient data to learn Fourier basis - More computationally intensive than NHITS

When selected: complex group. Used when the series is flagged as high complexity.

PatchTST (Patch Time Series Transformer)¶

What it is: PatchTST (Nie et al., 2023) applies the Transformer architecture to time series by segmenting the input into patches (sub-sequences) and treating each patch as a token. This dramatically reduces the attention complexity and allows longer context windows.

Architecture:

Patching: Divide the input window of length \(L\) into patches of length \(P\) with stride \(S\):

Number of patches: \(N = \lfloor (L - P) / S \rfloor + 2\)
Linear embedding: Each patch is projected to a \(d_{\text{model}}\)-dimensional embedding.
Positional encoding: Learnable or sinusoidal positions added.
Multi-head self-attention across patches (not across time steps):

\[\text{Attention}(Q, K, V) = \text{softmax}\left(\frac{QK^T}{\sqrt{d_k}}\right)V\]

Complexity: \(O(N^2 \cdot d) = O((L/P)^2 \cdot d)\) vs. \(O(L^2 \cdot d)\) for full Transformers. For \(L=520\), \(P=16\): ~33× reduction.
Head: Linear projection from patch tokens to the forecast horizon.

Parameters in ForecastAI:

Parameter	Low/Medium complexity	High complexity
`input_size`	`min(10 × horizon, n // 2)`	Same (longer context)
`patch_len`	16	16
`stride`	8 (50% overlap)	8
`n_layers`	2	3
`d_model`	64	128
`n_heads`	4	8
`learning_rate`	1e-4	1e-4

Strengths: - Longer context windows than NHITS/NBEATS: input_size = 10 × horizon (vs. 5×) - Efficient attention: captures long-range dependencies without quadratic scaling - Strong on series where distant history matters (long lead-time products)

Weaknesses: - Transformer training instability (requires careful tuning) - No covariate support in the NeuralForecast PatchTST implementation - Lower learning rate (1e-4) means slower convergence

When selected: complex and seasonal groups.

TFT (Temporal Fusion Transformer)¶

What it is: TFT (Lim et al., 2020) is a purpose-built Transformer for time series with explicit support for static covariates (item-level attributes), known future inputs (e.g., planned promotions, calendar features), and observed historical inputs. It produces interpretable attention weights.

Architecture components:

Variable Selection Networks (VSN): Learn to select which input features are most important at each time step. Uses gated residual networks.
Gated Residual Networks (GRN):

\[\text{GRN}(a, c) = \text{LayerNorm}(a + \text{GLU}(\phi_1 + \text{ELU}(\phi_1 \cdot a + \phi_2 \cdot c)))\]
LSTM encoder-decoder: Processes the selected features through sequence-to-sequence LSTM.
Gated skip connections: Allow the model to skip over unnecessary components.
Multi-head attention (interpretable):

\[\hat{A}(Q, K) = \text{softmax}\left(\frac{QK^T}{\sqrt{d_{\text{attn}}}}\right)\]

The attention weights \(A_{ij}\) show how much weight the model places on time step \(j\) when forecasting step \(i\). These are aggregated across heads to produce interpretable attention patterns.
Quantile output layer: Separate outputs for each quantile level.

Parameters in ForecastAI:

Parameter	Low/Medium complexity	High complexity
`input_size`	`min(5 × horizon, n // 2)`	Same
`hidden_size`	64	128
`dropout`	0.1	0.1
`num_attention_heads`	4	4
`learning_rate`	1e-3	1e-3

Strengths: - Handles static covariates (item attributes), known futures (promotions), and observed covariates simultaneously - Interpretable attention: identify which past periods drive the forecast - Designed from the ground up for supply chain time series

Weaknesses: - Most complex architecture: requires more data and training time - Benefits are only realized with rich covariate inputs — without covariates, simpler models often perform comparably - Sensitive to feature engineering quality

When selected: complex group. Best exploited when item-level attributes and calendar features are passed as covariates.

DeepAR¶

What it is: DeepAR (Salinas et al., 2020) is an autoregressive RNN that models the distribution of future values rather than point predictions. It trains a global model across all series, learning shared representations.

Architecture:

\[h_t = \text{LSTM}(h_{t-1},\; z_{t-1},\; x_t, \theta)\]

where \(z_{t-1}\) is the previous value (autoregressive conditioning) and \(x_t\) are covariates.

Probabilistic output: The LSTM hidden state is mapped to the parameters of a parametric distribution:

Gaussian: \(p(z_t | \mu_t, \sigma_t)\) — for continuous demand
Negative Binomial: \(p(z_t | \mu_t, \alpha_t)\) — for count data (integer demand)

Forecast sampling: Monte Carlo sampling from the predicted distribution produces trajectories:

\[\hat{y}_t^{(s)} \sim \mathcal{N}(\mu_t, \sigma_t^2) \quad \text{for } s = 1, \ldots, S\]

Quantiles are computed empirically from the sampled trajectories.

Global training: A single model is trained across all series in the portfolio. Each series is identified by a learned embedding (item-level features). This allows the model to transfer patterns learned from data-rich series to data-sparse series.

Parameters in ForecastAI:

Parameter	Low/Medium complexity	High complexity
`input_size`	`min(5 × horizon, n // 2)`	Same
`encoder_hidden_size`	64	128
`decoder_hidden_size`	64	128
`encoder_n_layers`	2	3
`decoder_n_layers`	2	3
`learning_rate`	1e-3	1e-3

Strengths: - Natively probabilistic: full predictive distributions, not just quantile approximations - Global model: good for large portfolios (hundreds to thousands of series) - Handles cold-start via shared embeddings - Proven at scale (Amazon production forecasting)

Weaknesses: - Requires many series for effective global training (few series → poor embeddings) - Slower than MLP-based methods - Distribution assumptions may not match all demand patterns (e.g., Gaussian for intermittent)

When selected: complex group. Most effective in large portfolios.

Foundation Model (TimesFM)¶

TimesFM (Google)¶

What it is: TimesFM (Time Series Foundation Model, Google Research, 2024) is a pre-trained transformer trained on 100 billion time points from diverse real-world datasets. It can forecast new series zero-shot (without any fine-tuning on the target series), making it uniquely suited for cold-start scenarios.

Figure: TimesFM zero-shot inference flow — native quantile output with a scaling fallback when native quantiles fail.

flowchart LR
    C["Context window<br/>512 points"] --> M["TimesFM<br/>200M-param transformer<br/>(pre-trained, zero-shot)"]
    M --> Q["Native quantile forecast<br/>Q10..Q99"]
    Q --> OK{"quantile output ok?"}
    OK -- yes --> OUT["Fitted quantiles"]
    OK -- no --> FB["Scaling fallback<br/>y * (1 + |q - 0.5| * 2)"]
    FB --> OUT

Architecture: A decoder-only transformer using patching (similar to PatchTST but at the foundation model scale):

Architecture component	Value
Model ID	`google/timesfm-1.0-200m-pytorch`
Total parameters	200 million
Context length	512
Horizon length	128
Input patch length	32
Output patch length	128
Transformer layers	20
Model dimensions	1280

Zero-shot operation: The model was pre-trained on a massive corpus of time series. At inference time, it takes a raw context window and directly outputs a forecast — no training or fine-tuning on the target series is needed.

Frequency hint: A frequency hint integer helps the model adapt to different data frequencies:

Data frequency	Frequency hint
High (H/D/B)	0
Medium (W/M)	1
Low (Q/Y/A)	2

Quantile outputs: TimesFM produces quantile forecasts natively. If the native quantile output fails, a scaling fallback is applied:

\[\hat{y}^{(q)} = \hat{y}^{(0.5)} \times \left(1 + |q - 0.5| \times 2\right)\]

Integration in ForecastAI: files/forecasting/foundation_models.py, class FoundationForecaster. The model is loaded once and cached in memory. GPU inference (CUDA) is used if available, otherwise CPU.

Backend detection:

backend = 'gpu' if torch.cuda.is_available() else 'cpu'

Parameters in ForecastAI:

Parameter	Value	Source
`context_length`	512	Config: `forecasting.timesfm.context_length`
`horizon_length`	128	Config: `forecasting.timesfm.horizon_length`
`quantiles`	[0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]	Config: `forecasting.timesfm.quantiles`

Strengths: - Zero-shot: no training required → works immediately on new series with no history - Trained on enormous diverse corpus: robust to unusual patterns - Very fast inference (no per-series training loop) - Excellent for New Product Introductions (NPI) with limited history

Weaknesses: - Cannot be fine-tuned per-series (black box) - May not outperform a well-tuned statistical model on series with stable, long history - Large memory footprint (200M parameter model loaded in GPU memory) - Quantile fallback is heuristic — less reliable than native quantile training

When selected: sparse_data, intermittent, seasonal, and standard groups. Always available as a fallback for any series where other methods are excluded.

Configuration:

forecasting:
  timesfm:
    model_name: timesfm-1.0-200m
    context_length: 512
    horizon_length: 128
    quantiles: [0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]

Benchmark / Fallback Models¶

These simple models serve as baselines and fallbacks when data is insufficient for statistical or neural models.

Model	Forecast	Use case
`HistoricAverage`	Mean of all historical values	Ultimate fallback; no pattern
`Naive`	Last observed value	Random walk baseline
`SeasonalNaive`	Last complete seasonal cycle	Sparse seasonal data
`SeasonalWindowAverage`	Mean of last \(w\) seasonal cycles	Smoother seasonal baseline

HistoricAverage is the characterization-layer recommended_methods fallback when all candidates are filtered out (it needs ≥ 2 observations).

At forecasting time, StatisticalForecaster.forecast_multiple_series (files/forecasting/statistical_models.py) applies a stronger last-resort: after eligibility filtering, if a series still has no eligible statistical method (e.g. a 2-observation series routed to AutoETS/AutoARIMA/AutoTheta/MSTL, all of which need ≥ 10), it is reassigned to Naive (needs only ≥ 1 observation) rather than aborting the whole pipeline with "0 forecasts produced". This keeps very short series forecastable where HistoricAverage alone would still fail (it requires ≥ 2).

Parameter Override System¶

All neural model parameters are individually overridable via the hyperparameter_overrides table (unique_id, method, overrides JSONB):

Overridable neural parameters: input_size, max_steps, learning_rate, batch_size, dropout, hidden_size, encoder_hidden_size, encoder_n_layers, decoder_hidden_size, decoder_n_layers, n_layers, d_model, n_heads, patch_len, stride, num_attention_heads, val_check_steps, early_stop_patience_steps.

Overridable statistical parameters: model, damped, phi (ETS); p, d, q, P, D, Q, max_p, max_q, etc. (ARIMA); decomposition_type (Theta).

When manual orders are provided for AutoARIMA, the model switches from AutoARIMA() to a fixed ARIMA(), bypassing the stepwise search entirely.