Time Series

A time series is a sequence of observations ordered in time: {x_1, x_2, ..., x_T}. Unlike tabular data where rows are i.i.d. samples, time series observations are temporally dependent — the order matters, past values inform future ones. Key properties: trend (long-term direction), seasonality (periodic patterns), noise (random fluctuations), and autocorrelation (correlation of a series with its own past). Standard ML algorithms assume i.i.d. data; applying them naively to time series ignores temporal structure and leads to data leakage and poor generalization.

A time series is (weakly) stationary if its mean, variance, and autocovariance are constant over time: E[X_t] = μ, Var(X_t) = σ², Cov(X_t, X_{t+k}) = γ(k). Stationarity is required by classical models (ARIMA, ACF analysis) because non-stationary series have time-varying statistics that make parameter estimation unreliable. Most real-world series are non-stationary (stock prices, temperature trends). Tests: Augmented Dickey-Fuller (ADF) — H0: unit root (non-stationary). Transformation to stationarity: differencing, log transformation, detrending, seasonal decomposition.

from statsmodels.tsa.stattools import adfuller
result = adfuller(series)
print('p-value:', result[1])  # p < 0.05 → stationary

Differencing computes the change between consecutive observations: Δx_t = x_t − x_{t-1}. It removes trends and makes non-stationary series stationary. First differencing removes linear trends; second differencing (Δ²x_t) removes quadratic trends. Seasonal differencing: Δ_s x_t = x_t − x_{t-s} removes periodic seasonality of period s. If a series requires d differences to become stationary, it's called I(d) (integrated of order d). ARIMA's 'I' stands for this integration order. Differencing is a simple, effective preprocessing step before applying classical time series models.

# First differencing
diff_series = series.diff().dropna()
# Seasonal differencing (period=12)
seasonal_diff = series.diff(12).dropna()

Autocorrelation (ACF): ACF(k) = Corr(X_t, X_{t-k}) — correlation between the series and its lag-k version. Measures total linear dependence between observations k steps apart. Partial Autocorrelation (PACF): PACF(k) — correlation between X_t and X_{t-k} after removing the effects of intermediate lags. ACF and PACF plots guide ARIMA model selection: MA(q) has ACF cutting off at lag q; AR(p) has PACF cutting off at lag p; ARMA has both tailing off. They reveal the temporal dependence structure and suggest appropriate model orders.

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plot_acf(series, lags=40)
plot_pacf(series, lags=40)

ARIMA(p,d,q) combines: AR(p) — autoregressive part (weighted sum of past p values), I(d) — differencing d times to achieve stationarity, MA(q) — moving average part (weighted sum of past q errors). The model: Δ^d x_t = c + Σ_{i=1}^p φ_i Δ^d x_{t-i} + Σ_{j=1}^q θ_j ε_{t-j} + ε_t. Box-Jenkins methodology: identify p,d,q from ACF/PACF plots → estimate parameters → diagnostic checking (Ljung-Box test on residuals). ARIMA is a classical workhorse for univariate forecasting, interpretable and well-understood but limited to linear relationships.

from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(series, order=(2,1,2))
result = model.fit()
forecast = result.forecast(steps=10)

SARIMA(p,d,q)(P,D,Q)_s extends ARIMA to handle seasonal patterns with period s. It adds seasonal AR (Φ_i), seasonal MA (Θ_j), and seasonal differencing (D) terms. Example: monthly data with yearly seasonality uses s=12. The model captures both within-season and cross-season patterns. Notation: (p,d,q) for non-seasonal; (P,D,Q)_s for seasonal. SARIMA is appropriate when ACF shows spikes at multiples of s (12, 24, 36 for monthly). It's widely used for retail sales, electricity consumption, and weather forecasting with seasonal data.

from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(series, order=(1,1,1), seasonal_order=(1,1,1,12))
result = model.fit()

Exponential smoothing forecasts by taking weighted averages of past observations, with weights decaying exponentially: ŷ_{t+1} = α x_t + α(1−α) x_{t-1} + α(1−α)² x_{t-2} + ... = α x_t + (1−α) ŷ_t. α ∈ (0,1) is the smoothing parameter. Small α = slow adaptation (smoother, less responsive); large α = fast adaptation (tracks recent values closely). Extensions: Holt's (double) adds trend component; Holt-Winters (triple) adds seasonality. ETS (Error-Trend-Seasonality) framework encompasses all variants and selects the best via AIC.

from statsmodels.tsa.holtwinters import ExponentialSmoothing
model = ExponentialSmoothing(series, trend='add', seasonal='add', seasonal_periods=12)
result = model.fit()

Standard k-fold cross-validation randomly shuffles data, violating temporal order and causing data leakage (future predicting past). Time series CV uses expanding window (walk-forward): train on 1...t, test on t+1...t+h, expand, repeat. Or rolling window: train on t-W...t, test on t+1. This respects temporal order and mimics real forecasting. scikit-learn's TimeSeriesSplit implements this. Never shuffle time series. Always test on the most recent period. The gap between train and test sets should match the real forecasting horizon to avoid leakage.

from sklearn.model_selection import TimeSeriesSplit
tss = TimeSeriesSplit(n_splits=5)
for train_idx, test_idx in tss.split(X):
    # train on X[train_idx], test on X[test_idx]

Decomposition separates a time series into components: Trend T_t (long-term direction), Seasonal S_t (periodic patterns), and Residual R_t (noise). Two models: additive (x_t = T_t + S_t + R_t, use when seasonal amplitude is constant) and multiplicative (x_t = T_t × S_t × R_t, use when seasonal amplitude grows with trend). Classical decomposition uses moving averages for trend; STL (Seasonal-Trend decomposition using LOESS) is more robust to outliers. Decomposition aids interpretation, anomaly detection, and preprocessing before forecasting models.

from statsmodels.tsa.seasonal import seasonal_decompose
result = seasonal_decompose(series, model='additive', period=12)
result.plot()

Rolling (sliding) window features extract statistics over a fixed past window of size w: rolling mean, std, min, max, lag values. Example: for row at time t, feature = mean(x_{t-w},...,x_{t-1}). This converts a time series into a tabular regression problem where each row is a snapshot with engineered temporal features. Must avoid look-ahead bias: window must use strictly past observations. Common features: rolling mean (trend), rolling std (volatility), lag-k (direct past value), diff (rate of change). This enables standard ML algorithms (gradient boosting, RF) on time series.

import pandas as pd
df['roll_mean_7'] = df['value'].rolling(7).mean().shift(1)  # shift to avoid leakage
df['lag_1'] = df['value'].shift(1)

The ADF test tests for a unit root (non-stationarity). H0: the series has a unit root (non-stationary). Ha: the series is stationary. The test regresses Δx_t on x_{t-1} and lagged differences. Test statistic is compared to critical values (or p-value). If p < 0.05, reject H0 (series is stationary). Augmentation: adding lagged differences to control for serial correlation. A small p-value means stationarity. Complement with KPSS test (H0: stationary) for confirmation — if ADF rejects H0 AND KPSS fails to reject, strong evidence of stationarity.

from statsmodels.tsa.stattools import adfuller
adf_stat, p_val, _, _, crit_vals, _ = adfuller(series)
print(f'ADF Stat: {adf_stat:.4f}, p-value: {p_val:.4f}')

Prophet (Facebook/Meta, 2017) is an additive model for time series with trend, seasonality, and holiday effects: y(t) = g(t) + s(t) + h(t) + ε(t). g(t): piecewise linear or logistic growth trend with automatic changepoint detection. s(t): Fourier series seasonality (daily, weekly, yearly). h(t): user-specified holiday effects. Advantages: handles missing data, outliers, and multiple seasonalities automatically; interpretable components; works well with only 100+ data points. Designed for business analysts without deep time series expertise. Parameters are fit via Stan (MAP estimation).

from prophet import Prophet
df = df.rename(columns={'date':'ds','value':'y'})
model = Prophet(yearly_seasonality=True)
model.fit(df)
future = model.make_future_dataframe(periods=30)
forecast = model.predict(future)

Long Short-Term Memory networks (LSTMs) process sequences using cell state and gating: forget gate discards irrelevant past, input gate adds new information, output gate determines what to expose. For time series: input sequences of length W, output forecast H steps ahead. LSTMs can capture complex nonlinear temporal dependencies that ARIMA misses. Training: teacher forcing (feed true past), then inference with rolling prediction. Key choice: univariate (one feature, one output) vs multivariate (multiple features, e.g., exogenous variables). LSTMs outperform ARIMA on nonlinear and volatile series but require more data.

import torch.nn as nn
class LSTM_TS(nn.Module):
    def __init__(self, input=1, hidden=64, out=1):
        super().__init__()
        self.lstm = nn.LSTM(input, hidden, batch_first=True)
        self.fc = nn.Linear(hidden, out)
    def forward(self, x):
        out, _ = self.lstm(x)
        return self.fc(out[:,-1,:])

TCN (Bai et al. 2018) uses causal dilated 1D convolutions for sequence modeling. Causal: conv at step t uses only x_{≤t} (no future leakage). Dilated: exponentially growing dilation rates (1,2,4,8,...) give receptive field exponential in depth without many layers — TCN with 10 layers and dilation up to 512 has receptive field 1024. Advantages over LSTM: parallel training (no sequential dependency), stable gradients (no vanishing). TCN matches or beats LSTM on many time series benchmarks. Residual connections handle deep networks. WaveNet (audio generation) uses the same principle.

Granger causality tests whether time series X helps predict time series Y beyond Y's own past. Procedure: (1) Fit VAR(p) model predicting Y from its own lags only → RSS_1. (2) Fit VAR(p) predicting Y from its own lags + X's lags → RSS_2. (3) F-test whether RSS_2 < RSS_1 significantly. If yes, X 'Granger-causes' Y. Important caveat: Granger causality ≠ true causality — it measures predictive ability, not structural causation. Confounding variables can appear to Granger-cause. Used in econometrics, neuroscience (brain connectivity), and climate science.

from statsmodels.tsa.stattools import grangercausalitytests
grangercausalitytests(data[['Y','X']], maxlag=4)

VAR(p) models multiple time series simultaneously, where each variable is a linear function of its own past p lags AND the past p lags of all other variables: y_t = c + Σ_{i=1}^p A_i y_{t-i} + ε_t. This captures cross-series dependencies that univariate models miss. Example: GDP, inflation, interest rate interact — VAR captures all these relationships. Model selection: use AIC/BIC to choose p. For n series and p lags, VAR has n² × p parameters — can overfit with many series. VARMA adds MA terms; VEC (Vector Error Correction) handles cointegrated non-stationary series.

from statsmodels.tsa.vector_ar.var_model import VAR
model = VAR(data)  # data: T x n array
result = model.fit(maxlags=4, ic='aic')

Rolling mean (simple moving average, SMA): unweighted average of the past w observations: SMA_t = (x_t + x_{t-1} + ... + x_{t-w+1}) / w. Equal weight to all past w values. MA in ARIMA: a model component where x_t = μ + ε_t + θ_1 ε_{t-1} + ... (weighted past errors, not values). They sound similar but are different: SMA smooths the series for visualization/preprocessing; MA(q) is a stochastic process model for autocorrelated noise. Exponential Moving Average (EMA) gives exponentially decaying weights — more responsive to recent values than SMA.

# SMA
sma = series.rolling(window=7).mean()
# EMA (alpha=0.3 ↔ span=5.67)
ema = series.ewm(alpha=0.3).mean()

Two non-stationary series X_t and Y_t are cointegrated if a linear combination aX_t + bY_t is stationary — they share a long-run equilibrium relationship even though individually they're non-stationary. Example: stock price of Coca-Cola and Pepsi may both trend randomly but stay close together. Test: Engle-Granger (regress Y on X, test residuals for stationarity) or Johansen test. Cointegrated series should be modeled with Vector Error Correction Models (VECM), not standard VAR on levels. Pairs trading exploits cointegration: short the overvalued, long the undervalued, profit from mean reversion.

Time series anomaly detection identifies points/subsequences that deviate from expected patterns. Approaches: (1) Statistical: z-score thresholding (mean ± 3σ), Isolation Forest on rolling features. (2) Forecast-based: train ARIMA/LSTM, flag high residuals (|x_t − x̂_t| > k·σ). (3) Unsupervised: autoencoder with high reconstruction error. (4) Distance-based: matrix profile (STUMPY) identifies similar subsequence patterns. Applications: server metrics (CPU spike), fraud detection, industrial equipment failure (vibration sensors). Anomaly definition is context-dependent — seasonal peaks may be normal.

# Forecast-based anomaly detection
residuals = series - model.predict()
threshold = residuals.mean() + 3 * residuals.std()
anomalies = series[abs(residuals) > threshold]

The Discrete Fourier Transform (DFT) decomposes a time series into a sum of sinusoids at different frequencies: X_k = Σ_{t=0}^{N-1} x_t e^{-2πi kt/N}. The power spectrum |X_k|² shows which frequencies dominate — peaks indicate seasonality periods. A peak at k=1/12 (monthly data) signals yearly seasonality. Fourier features (sin/cos at known frequencies) can be added as features for ML models. Prophet uses Fourier series to model seasonality: s(t) = Σ_{n=1}^N [a_n cos(2πnt/P) + b_n sin(2πnt/P)].

import numpy as np
fft = np.fft.fft(series)
freqs = np.fft.fftfreq(len(series), d=1)  # d=sampling interval
power = np.abs(fft)**2

N-BEATS (Neural Basis Expansion Analysis for Interpretable Time Series Forecasting, Oreshkin et al. 2020) uses a fully connected network with backward (backcast) and forward (forecast) outputs per block. Stack of blocks: each block produces a backcast (explains past) and forecast; residuals propagate forward. Trend and seasonality stacks use polynomial and Fourier basis expansions, giving interpretability. N-BEATS won the M4 forecasting competition without any domain-specific feature engineering. It's purely data-driven, interpretable, and achieves strong results on diverse univariate forecasting benchmarks.

Look-ahead bias occurs when features or preprocessing use future information unavailable at prediction time, artificially inflating model performance. Common mistakes: scaling with global mean/std (includes test data), rolling stats without shifting (rolling mean of row t uses x_t itself), computing target-correlated features using future data. Prevention: (1) Shift rolling features by one step (shift(1)). (2) Fit scalers on training data only. (3) Use TimeSeriesSplit instead of random split. (4) Audit every feature for temporal validity. Look-ahead bias is the most common error in time series ML and leads to models that fail completely in production.

# WRONG: includes current value
df['roll'] = df['value'].rolling(7).mean()
# CORRECT: shift by 1
df['roll'] = df['value'].rolling(7).mean().shift(1)

Vanilla Transformers struggle with long time series due to O(n²) attention. Informer (2021) uses ProbSparse attention (O(n log n)) — select only the top-u queries with highest attention score sparsity. Autoformer (2022) replaces attention with auto-correlation (FFT-based dependency discovery) + decomposition blocks (trend + seasonality). PatchTST (2023) patches time series into subseries-level tokens (like ViT) + channel independence, achieving SOTA. TimesNet treats 1D time series as 2D by reshaping (FFT-detected period as one axis). These architectures handle hundreds to thousands of time steps.

MAPE = (1/T) Σ |y_t − ŷ_t| / |y_t| × 100%. Advantages: scale-independent, easily interpretable as percentage error. Limitations: undefined when y_t = 0 (division by zero); asymmetric (underforecasts penalized more than overforecasts); biased toward low values. Alternatives: sMAPE (symmetric MAPE): 2|y−ŷ|/(|y|+|ŷ|); MASE (Mean Absolute Scaled Error): MAE normalized by in-sample naïve forecast MAE — scale-independent and well-defined. MASE < 1 means better than naïve; used in M4 competition. Choose metric based on business cost of over/under-forecasting.

import numpy as np
def mape(y, y_hat):
    return np.mean(np.abs((y - y_hat) / y)) * 100
def mase(y, y_hat, y_train):
    naive_mae = np.mean(np.abs(np.diff(y_train)))
    return np.mean(np.abs(y - y_hat)) / naive_mae

Naïve forecast: predict next value = last observed value: ŷ_{t+h} = x_t. Simple baseline; works well for random walks (stock prices). Seasonal naïve: predict next period = same period last season: ŷ_{t+h} = x_{t+h-s} (e.g., predict January using last January's value). Strong baseline for seasonal data. Any model should beat naïve forecasts to be useful. Random walk with drift: ŷ_{t+h} = x_t + h·(x_t − x_{t-n})/n. These baselines are surprisingly competitive — many complex models fail to significantly outperform them on certain datasets (especially stock data).

TFT (Lim et al. 2020) is a Transformer-based architecture for multi-horizon probabilistic time series forecasting. Key components: gating mechanisms select relevant variables; variable selection networks identify important inputs; static covariate encoders condition on time-invariant features; LSTM encoder for local processing; self-attention for long-range dependencies; quantile outputs for prediction intervals. TFT handles: multiple related series, exogenous variables, static metadata, different frequencies. It won M5 competition variants and provides interpretable attention patterns showing which time steps and variables matter most for each forecast.

Changepoint detection identifies times when the statistical properties of a time series change: mean shift, variance change, or trend break. Methods: CUSUM (cumulative sum of deviations — triggers when sum exceeds threshold), PELT (Pruned Exact Linear Time — exact minimum across all possible changepoints), Prophet's changepoint detection (MAP estimation of piecewise linear trends). Online vs offline: offline uses full series; online detects changes in real-time. Applications: detecting business events (marketing campaign effect), equipment degradation, policy change impacts, or A/B test effects.

import ruptures as rpt
algo = rpt.Pelt(model='rbf').fit(signal)
breakpoints = algo.predict(pen=10)

GARCH(p,q) (Generalized Autoregressive Conditional Heteroskedasticity) models time-varying variance (volatility): σ²_t = ω + Σ_{i=1}^p α_i ε²_{t-i} + Σ_{j=1}^q β_j σ²_{t-j}. This captures 'volatility clustering' — periods of high volatility tend to cluster (stock markets). ARCH(q) is the special case with β_j=0. GARCH is used for financial risk (VaR), option pricing, and portfolio management. The residual ε_t = σ_t z_t where z_t ~ N(0,1). GARCH is typically combined with an ARIMA mean equation for joint mean+variance modeling.

from arch import arch_model
gm = arch_model(returns, vol='GARCH', p=1, q=1)
res = gm.fit()
forecasts = res.forecast(horizon=5)

Conformal prediction provides distribution-free prediction intervals with guaranteed coverage. For time series: (1) Use an exchange exchangeability assumption (approximate for TS via recent calibration window). (2) Fit a point forecaster and compute residuals on calibration set. (3) Prediction interval: ŷ_{t+h} ± quantile(|residuals|, 1−α). Guarantees: P(y_{t+h} ∈ interval) ≥ 1−α without distributional assumptions. EnCQR (Ensemble Cross Quantile Regression with Conformal) handles non-exchangeable time series by weighting calibration samples by recency. Important for risk-aware forecasting in medical, financial, and safety-critical applications.

Multi-step forecasting predicts h steps ahead. Strategies: (1) Recursive (MIMO 1-step): forecast 1 step, feed it as input, repeat h times — errors accumulate. (2) Direct (DIRMO): train separate model per horizon h — no error accumulation but ignores dependencies between horizons. (3) MIMO (Multi-Input Multi-Output): one model maps input window to all h outputs simultaneously — efficient and consistent. MIMO is typically best. Neural nets naturally support MIMO by outputting vector of length h. For ARIMA: analytic multi-step forecasts exist but CI grows with h. Long-horizon forecasting (h >> 1) is fundamentally harder than 1-step.

# MIMO forecasting with PyTorch
# Model outputs h=12 steps from window of w=48
class MIMOForecast(nn.Module):
    def __init__(self, w, h):
        super().__init__()
        self.fc = nn.Linear(w, h)
    def forward(self, x):  # x: (batch, w)
        return self.fc(x)  # (batch, h)

The matrix profile (MP) is a vector where MP[i] = Euclidean distance between subsequence starting at i and its nearest neighbor in the same time series (excluding trivial self-match). Computing all pairwise subsequence distances is O(n²) naïvely; STAMP/STOMP/STUMPY compute MP in O(n² log n) via FFT. MP enables: anomaly detection (high MP = distant from all neighbors), motif discovery (low MP = frequently recurring pattern), discord detection, and segmentation. STUMPY is the Python library. MP is exact, non-parametric, and works without domain knowledge.

import stumpy
mp = stumpy.stump(series.values, m=50)  # m=subsequence length
# mp[:,0] = distances, mp[:,1] = nearest-neighbor indices

State space models (SSMs) represent time series via a latent state x_t: Transition: x_t = F x_{t-1} + Q w_t (how state evolves). Observation: y_t = H x_t + R v_t (how state relates to observations). Kalman Filter computes optimal state estimates recursively under Gaussian assumptions. SSMs unify: ARIMA (specific F,H,Q,R), exponential smoothing, local level model, and structural time series. They handle missing data naturally and allow uncertainty quantification. Modern deep SSMs (Mamba, S4) use structured SSMs for long-sequence modeling, achieving competitive performance with Transformers.

from filterpy.kalman import KalmanFilter
kf = KalmanFilter(dim_x=2, dim_z=1)
kf.F = np.array([[1,1],[0,1]])  # constant velocity
kf.H = np.array([[1,0]])

A point forecast gives a single predicted value ŷ_{t+h}. Probabilistic forecasts quantify uncertainty: (1) Prediction intervals: [lower, upper] with coverage probability 1−α (e.g., 95% PI means 95% of future values should fall within). (2) Quantile forecasts: predict multiple quantiles (10th, 50th, 90th percentile). (3) Full distribution: output parameters of a parametric distribution (Gaussian mean+std, or student-t). Probabilistic forecasts are more useful for decision-making (inventory management needs upper-bound estimate; risk management needs tails). CRPS (Continuous Ranked Probability Score) evaluates probabilistic forecasts.

Seasonal adjustment removes predictable seasonal patterns from a time series to reveal the underlying trend and irregular components. Example: retail sales in December are always high — adjusting removes this effect to see whether underlying demand is growing. Methods: (1) X-13ARIMA-SEATS (US Census Bureau standard). (2) STL decomposition then remove seasonal component. Seasonally adjusted series enable fair comparison across months. Used by government agencies (GDP, unemployment figures) and businesses. The caveat: seasonal adjustment requires the seasonal pattern to be stable — changing seasons (COVID disruption) can break models.

from statsmodels.tsa.seasonal import STL
stl = STL(series, period=12)
result = stl.fit()
deseasonalized = series - result.seasonal

AR(p) (Autoregressive): x_t = φ_1 x_{t-1} + ... + φ_p x_{t-p} + ε_t. Current value is a linear function of p past values + white noise. ACF tails off; PACF cuts off at lag p. MA(q) (Moving Average): x_t = ε_t + θ_1 ε_{t-1} + ... + θ_q ε_{t-q}. Current value is a linear function of current + q past error terms. ACF cuts off at lag q; PACF tails off. Key difference: AR has infinite memory (current value affects all future); MA has finite memory (shock disappears after q steps). ARMA(p,q) combines both for parsimony.

Exogenous variables are external inputs that help predict the target series but are not modeled themselves. ARIMAX adds exogenous variables to ARIMA: y_t = ARIMA + β x_t. Example: predicting electricity demand (y) using temperature (x) — temperature helps but is not predicted. SARIMAX is the seasonal variant. In ML models, exogenous features are simply additional columns. Prophet allows regressors: model.add_regressor('temperature'). Challenge: to forecast multiple steps ahead, you also need future exogenous values — often themselves uncertain. This requires forecasting the exogenous variables too or using known future values (e.g., planned promotions).

from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(y, exog=X_train, order=(1,1,1))
result = model.fit()
forecast = result.forecast(steps=10, exog=X_test)

Local models (ARIMA, ETS, Prophet) fit a separate model per time series — capturing each series' unique patterns but cannot share information. Global models (LSTM, Transformer, N-BEATS) train one model on all series simultaneously — sharing statistical strength, especially useful for short series with limited data. Global models can generalize to new, unseen series (zero-shot forecasting). Challenge: series can have very different scales (normalize per-series), patterns, or frequencies. LightGBM (gradient boosted tree) used as a global model with lag features is a strong baseline. The winning M5 competition solution used LightGBM globally.

The M competitions (M1–M5, Makridakis) are the benchmark forecasting competitions that shaped the field. M1 (1982): 1001 series, showed simple methods beat complex ones. M3 (2000): 3003 series, Theta method won, hybrid methods strong. M4 (2018): 100,000 series, ES-RNN hybrid won — showed ML must be combined with statistical methods. M5 (2020): Walmart sales data (42,840 series), probabilistic forecasting, LightGBM-based methods dominated. Key lesson from M competitions: ensemble of methods consistently outperforms single models; feature engineering and domain knowledge matter enormously; pure deep learning often loses to well-tuned statistical+ML hybrids.

For gradient boosted tree models on time series (LightGBM, XGBoost), feature importance reveals which temporal features drive predictions: gain importance (average improvement in split criterion), permutation importance (drop in performance when feature is shuffled). Common important features: recent lags (lag_1, lag_7), rolling statistics (mean/std over different windows), date features (day_of_week, month, is_holiday), target-encoded categoricals. SHAP values give per-sample attribution — important for time series debugging (e.g., 'this high forecast is driven by last week's spike'). Temporal SHAP plots reveal how feature importance evolves over time.

import shap
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)
shap.summary_plot(shap_values, X_test)

Hierarchical time series are organized in a tree: national total → state → city → store. Forecasts at each level must be coherent: national = sum of states. Bottom-up: forecast lowest level, aggregate up (coherent but noisy). Top-down: disaggregate top-level forecast down (smooth but loses local patterns). Optimal reconciliation (MinT, Wickramasuriya 2019): forecast all levels independently then project to the coherent subspace minimizing covariance of revision. Used in retail (SKU to store to region), energy (substations to grid), and web traffic (pages to sites to total). MinT is available in the `hts` R package and `sktime` Python library.

The ACF plot shows Corr(X_t, X_{t-k}) for k=1,2,...,max_lag on the y-axis vs. lag k on x-axis. Blue shaded area is the 95% confidence interval: ±1.96/√T. Spikes outside this band are statistically significant autocorrelations. AR processes: ACF tails off geometrically; PACF cuts off at p. MA processes: ACF cuts off at q; PACF tails off. White noise: ACF ≈ 0 for all lags (all inside CI). Seasonal patterns: ACF shows spikes at multiples of season period s (12, 24, 52). ACF/PACF plots are the primary diagnostic tool before fitting ARIMA.

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(series, lags=40, alpha=0.05)

Rolling origin (walk-forward) evaluation simulates real forecasting by iterating: train on data up to time t, forecast h steps ahead, advance origin by one step. Collects forecast errors at each origin: e_{t+h} = y_{t+h} − ŷ_{t+h|t}. Aggregate errors across origins for RMSE/MAE/MAPE. Key parameter: gap (between train end and test start) should match the operational forecast horizon. A one-step-ahead gap for h=12-step forecasting overestimates performance. Rolling origin evaluation is the correct protocol for model selection in time series — it measures out-of-sample performance honestly.

Box-Jenkins is a systematic approach for ARIMA model identification: (1) Identification — check stationarity (ADF test, unit root); difference if needed; plot ACF/PACF to guess p and q. (2) Estimation — estimate parameters via MLE or conditional least squares. (3) Diagnostic checking — test residuals for white noise (Ljung-Box test: H0 that residuals are i.i.d.); check residual ACF/PACF; assess normality. If residuals have structure, revise model. Repeat until residuals are white noise. Auto-ARIMA (statsmodels, pmdarima) automates this via information criteria (AIC/BIC minimization).

from pmdarima import auto_arima
model = auto_arima(series, seasonal=True, m=12, information_criterion='aic', trace=True)
print(model.summary())

In-sample fit measures how well the model fits the data it was trained on — always optimistic because the model has already seen this data. Out-of-sample forecast measures performance on future unseen data — the true test. For time series: reserve the last h observations as the test set (holdout period); train on the rest. Never use the test set for model selection. A model can have near-perfect in-sample fit (overfit) but poor out-of-sample performance. Information criteria (AIC/BIC) penalize complexity and serve as in-sample proxies for out-of-sample performance.

# Train-test split preserving temporal order
n = len(series)
train, test = series[:int(0.8*n)], series[int(0.8*n):]

DTW is a similarity measure for time series that allows elastic alignment — matching sequences even if they are shifted or stretched in time. Standard Euclidean distance requires same-length sequences aligned point-to-point. DTW finds the optimal alignment (warping path) using dynamic programming: DTW(i,j) = ||x_i − y_j|| + min(DTW(i-1,j), DTW(i,j-1), DTW(i-1,j-1)). Lower DTW = more similar. Used for: time series classification (k-NN with DTW), clustering, and anomaly detection. Limitation: O(n²) computation per pair; FastDTW approximates in O(n). DTW is the gold-standard distance for irregular-speed or phase-shifted time series.

from dtaidistance import dtw
distance = dtw.distance(series_a, series_b)

The KPSS (Kwiatkowski-Phillips-Schmidt-Shin) test has H0: series is stationary, H_a: series has a unit root. Opposite of ADF (where H0 is non-stationarity). KPSS statistic measures the variance of partial sums; high variance → non-stationary → reject H0. Use both ADF and KPSS together: if ADF rejects H0 AND KPSS fails to reject → clearly stationary; if both reject H0 → might be fractionally integrated; if ADF fails to reject AND KPSS fails to reject → non-stationary; both failing to reject → conflicting (might need more data). Two KPSS variants: level stationarity (no trend) and trend stationarity (stationary around trend).

from statsmodels.tsa.stattools import kpss
stat, p_val, lags, crit = kpss(series, regression='c')
print(f'KPSS stat: {stat:.4f}, p-value: {p_val:.4f}')

Classical decomposition uses a centered moving average for trend and additive/multiplicative decomposition for seasonal component — inflexible and not robust to outliers. STL (Seasonal and Trend decomposition using LOESS) uses locally weighted regression (LOESS) for both trend and seasonal components: more flexible, handles non-linear trends, and is robust to outliers (they go into the remainder). STL also allows the seasonal component to change over time. STL parameters: seasonal window (s-window), trend window (t-window). Used as preprocessing before forecasting (remove seasonality, model trend+remainder, reseasonalize). Stats libraries: statsmodels.tsa.seasonal.STL.

from statsmodels.tsa.seasonal import STL
stl_result = STL(series, period=12, seasonal=7).fit()
stl_result.plot()

Multivariate time series has multiple co-evolving variables {x_1t, x_2t, ..., x_Kt} observed at each time step. Classical model: VAR(p) — vector autoregression. Each variable is regressed on its own and all others' past p lags. Neural models: multivariate LSTM (all variables in the input vector), Temporal Fusion Transformer (handles exogenous + static features), N-BEATS-M (multivariate extension). Key challenge: cross-variable dependencies (one variable leads another). Feature: channel-independent vs channel-mixing — some Transformer models (PatchTST) treat each channel independently for scalability; others (iTransformer) model cross-channel interactions explicitly.

# Multivariate LSTM: all variables as features
# X: (samples, seq_len, n_variables)
# y: (samples, horizon, n_variables or 1)

The Ljung-Box test checks whether autocorrelations of residuals at multiple lags are collectively different from zero. H0: residuals are white noise (no autocorrelation up to lag h). Test statistic: Q = T(T+2) Σ_{k=1}^h r_k²/(T−k), where r_k is residual autocorrelation at lag k. Under H0: Q ~ χ²(h−p−q). If p-value < 0.05, reject H0 — residuals have significant autocorrelation, meaning the model hasn't captured all temporal structure and needs revision. Run after fitting ARIMA or any time series model. A well-fit model should have Ljung-Box p-value > 0.05.

from statsmodels.stats.diagnostic import acorr_ljungbox
result = acorr_ljungbox(residuals, lags=10, return_df=True)
print(result)  # check p-values > 0.05

The spectral density S(f) of a stationary time series decomposes variance by frequency f. A peak at frequency f_0 means the series has a cycle of period 1/f_0. For monthly data with yearly seasonality: peak at f=1/12 cycles per month. For weekly data (daily observations): peak at f=1/7. The periodogram (sample spectral density via FFT) estimates S(f). Smoothed periodogram (Bartlett, Welch) reduces estimation variance. Spectral analysis can detect multiple hidden periodicities (daily + weekly + annual) simultaneously. Used in detecting seasonality before model selection, and in frequency-domain feature engineering for ML.

import numpy as np
from scipy.signal import periodogram
f, psd = periodogram(series, fs=1.0)  # fs=1 for unit time steps
peak_period = 1/f[np.argmax(psd)]