Time Series Forecasting: From ARIMA to Production-Grade Deep Learning

Most tutorials stop at ARIMA on toy datasets. Meanwhile, real-world temporal data from IoT sensors and real-time pipelines drives operational decisions across finance, supply chain, energy, and monitoring. This article bridges the gap between textbook models and production reality: how to engineer features, validate without leaking the future, choose between classical and deep learning approaches, and deploy a system that doesn't break when the data distribution shifts. We'll cover the full stack, from statistical foundations to MLOps for time series.

The Nature of Time Series: Stationarity, Trend, Seasonality, and Noise

A time series is a sequence of observations ordered in time. Unlike cross-sectional data, time series have a natural temporal ordering, which means observations close together in time are more related than those far apart. This dependency is the core property that differentiates time series analysis from standard regression. Every series can be decomposed into four components: trend (long-term direction), seasonality (repeating cycles of fixed period), cyclical patterns (longer, non-fixed cycles), and residual noise (irreducible random variation). For example, daily temperature data shows a clear annual seasonality and a possible warming trend, while stock prices are dominated by noise and trend with no fixed seasonality.

Stationarity is a foundational concept: a stationary series has constant mean, constant variance, and autocovariance that depends only on lag, not on time. Most forecasting models require stationarity or explicit differencing to achieve it. The Augmented Dickey-Fuller (ADF) test is the standard statistical test for stationarity; a p-value below 0.05 suggests the series is stationary. Non-stationary series often exhibit trends or changing variance, which must be removed via differencing (e.g., first difference y_t - y_{t-1}) or transformations like log or Box-Cox. For instance, monthly airline passenger counts are non-stationary due to increasing trend and multiplicative seasonality; log-transforming and differencing at lag 12 yields a stationary series.

Seasonality is a special case of periodic behavior with known, fixed period (e.g., 12 for monthly, 4 for quarterly, 7 for daily). It must be modeled explicitly, either by seasonal differencing or by including seasonal dummies. Noise, or irregular component, is the residual after removing trend and seasonality; it should ideally resemble white noise (zero mean, constant variance, no autocorrelation). The Ljung-Box test checks whether residuals are independently distributed. In practice, a good model leaves residuals that are indistinguishable from white noise. Understanding these components is not academic—it directly dictates which forecasting model class (ARIMA, ETS, Prophet, or deep learning) is appropriate.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller

# Generate synthetic monthly data with trend + seasonality + noise
np.random.seed(42)
dates = pd.date_range('2015-01-01', periods=60, freq='M')
trend = np.linspace(100, 200, 60)
seasonality = 10 * np.sin(2 * np.pi * np.arange(60) / 12)
noise = np.random.normal(0, 5, 60)
series = trend + seasonality + noise

# Decompose
df = pd.DataFrame({'value': series}, index=dates)
result = seasonal_decompose(df['value'], model='additive', period=12)

# ADF test
adf_stat, p_value, _, _, critical_values, _ = adfuller(series)
print(f"ADF Statistic: {adf_stat:.3f}")
print(f"p-value: {p_value:.5f}")
print(f"Critical values: {critical_values}")
if p_value < 0.05:
    print("Series is stationary (reject H0)")
else:
    print("Series is non-stationary (fail to reject H0)")

Classical Forecasting Models: ARIMA, ETS, and Their Assumptions

ARIMA (AutoRegressive Integrated Moving Average) is the standard tool of univariate time series forecasting. It combines three components: autoregressive (AR) terms that model the series as a linear combination of its own lags, differencing (I) to achieve stationarity, and moving average (MA) terms that model the error as a linear combination of past forecast errors. The notation ARIMA(p,d,q) specifies the order of each component. For example, ARIMA(2,1,2) means two AR lags, one differencing, and two MA lags. The model assumes linearity, stationarity after differencing, and that residuals are white noise. It is estimated via maximum likelihood or conditional sum-of-squares. In practice, the Box-Jenkins methodology guides model selection: identify p,d,q via ACF/PACF plots, estimate parameters, diagnose residuals, and iterate.

ETS (Error, Trend, Seasonal) models are exponential smoothing state-space models that handle trend and seasonality directly without differencing. The notation ETS(A,A,M) means additive error, additive trend, multiplicative seasonality. ETS models are robust to non-stationarity because they model the level, trend, and seasonal components explicitly. They are particularly effective for series with clear trend and seasonality, and they automatically handle missing values. Unlike ARIMA, ETS does not require stationarity, but it assumes the components evolve smoothly. The model is estimated via maximum likelihood, and forecast intervals are derived from the state-space structure.

Both models have strong assumptions: ARIMA assumes linear relationships and constant variance; ETS assumes additive or multiplicative components that evolve slowly. Neither handles complex nonlinear patterns, multiple seasonalities, or external regressors natively (though ARIMAX and regressor variants exist). In practice, ARIMA often outperforms on economic and financial data with complex autocorrelation, while ETS excels on retail demand and inventory data with clear seasonal patterns. The M3 and M4 forecasting competitions showed that simple methods like ETS and ARIMA often beat more complex machine learning models on univariate data. Always compare AIC or BIC for model selection, and never trust a model that fails residual diagnostics.

import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt

# Generate synthetic data with trend and seasonality
np.random.seed(42)
dates = pd.date_range('2020-01-01', periods=100, freq='M')
trend = np.linspace(50, 150, 100)
seasonality = 10 * np.sin(2 * np.pi * np.arange(100) / 12)
noise = np.random.normal(0, 3, 100)
series = trend + seasonality + noise

# Train/test split
train = series[:80]
test = series[80:]

# ARIMA(2,1,2) model
arima_model = ARIMA(train, order=(2,1,2))
arima_fit = arima_model.fit()
arima_forecast = arima_fit.forecast(steps=20)

# ETS(A,A,A) model
ets_model = ExponentialSmoothing(train, trend='add', seasonal='add', seasonal_periods=12)
ets_fit = ets_model.fit()
ets_forecast = ets_fit.forecast(steps=20)

print("ARIMA AIC:", arima_fit.aic)
print("ETS AIC:", ets_fit.aic)
print("ARIMA RMSE:", np.sqrt(np.mean((arima_forecast - test)**2)))
print("ETS RMSE:", np.sqrt(np.mean((ets_forecast - test)**2)))

Feature Engineering for Time Series: Lags, Rolling Windows, Calendar Features

Feature engineering transforms raw time series into a supervised learning problem. The most fundamental features are lagged values: y_{t-1}, y_{t-2}, ..., y_{t-k}. These capture autocorrelation and allow models like XGBoost or linear regression to learn temporal dependencies. The optimal number of lags can be determined from the partial autocorrelation function (PACF): significant PACF values indicate which lags to include. For monthly data, lag 12 is almost always significant due to seasonality. Rolling window statistics (mean, std, min, max over a window of size w) capture local trends and volatility. For example, a 7-day rolling mean smooths daily noise, while a 30-day rolling standard deviation captures volatility regimes. Exponentially weighted moving averages (EWMA) give more weight to recent observations and are often more robust than simple rolling windows.

Calendar features encode time-based patterns: hour of day, day of week, month, quarter, day of year, week of year, and boolean flags for holidays or special events. These are critical for data with strong human-driven patterns (retail, web traffic, energy consumption). For instance, retail sales spike on weekends and during holidays; ignoring day-of-week leads to large errors. Use cyclical encoding (sin/cos) for circular features like hour or month to preserve periodicity: sin(2π hour/24), cos(2π hour/24). This prevents models from treating 23 and 0 as far apart. Additionally, time since an event (e.g., days since last promotion) can capture carryover effects.

More advanced features include Fourier terms for multiple seasonalities (e.g., daily and weekly cycles), lagged differences (y_{t-1} - y_{t-2}) to capture momentum, and interaction features between lags and calendar variables. For example, the effect of a holiday may depend on the day of week. In practice, feature engineering for time series is more art than science: start with domain knowledge, add lags from PACF, add calendar features, then use feature importance from a tree model to prune. Beware of lookahead bias: never use future information to create features. Rolling windows must be computed only on past data at each time step. In production, implement feature computation as a pipeline that updates incrementally to avoid recomputing from scratch.

import pandas as pd
import numpy as np

# Generate daily data with weekly seasonality
dates = pd.date_range('2023-01-01', periods=365, freq='D')
values = 10 + 5 * np.sin(2 * np.pi * dates.dayofyear / 7) + np.random.normal(0, 2, 365)
df = pd.DataFrame({'date': dates, 'value': values})

# Lag features (from PACF: lag 1 and 7 significant)
df['lag_1'] = df['value'].shift(1)
df['lag_7'] = df['value'].shift(7)

# Rolling window features
df['rolling_mean_7'] = df['value'].rolling(window=7).mean()
df['rolling_std_7'] = df['value'].rolling(window=7).std()
df['ewma_alpha_0.3'] = df['value'].ewm(alpha=0.3).mean()

# Calendar features
df['day_of_week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['is_weekend'] = (df['day_of_week'] >= 5).astype(int)

# Cyclical encoding for day_of_week
df['dow_sin'] = np.sin(2 * np.pi * df['day_of_week'] / 7)
df['dow_cos'] = np.cos(2 * np.pi * df['day_of_week'] / 7)

# Drop NaN rows from shifting
df_clean = df.dropna()
print(df_clean.head())
print(df_clean.shape)

Evaluation Without Leakage: Walk-Forward Validation and Metrics

Standard k-fold cross-validation leaks temporal information: future data contaminates training sets, leading to overly optimistic performance. Time series requires walk-forward validation (also called time series cross-validation or rolling origin evaluation). The procedure: train on an expanding window of past data, predict the next h steps, then expand the training window to include those h steps, and repeat. The number of folds is typically 5-10, with a fixed forecast horizon. For example, with 100 data points and horizon 10, fold 1 trains on points 1-70, forecasts 71-80; fold 2 trains on 1-80, forecasts 81-90; etc. This mimics real-world forecasting where you only know the past. The final error is the average across all forecast steps.

Choosing the right metric is critical. For point forecasts, RMSE (Root Mean Squared Error) penalizes large errors quadratically, making it sensitive to outliers. MAE (Mean Absolute Error) is robust to outliers but ignores magnitude. MAPE (Mean Absolute Percentage Error) is scale-independent but undefined when actuals are zero and asymmetric (over-penalizes negative errors). SMAPE (symmetric MAPE) addresses asymmetry but still has issues with near-zero values. MASE (Mean Absolute Scaled Error) compares against a naive forecast (e.g., seasonal naive) and is scale-independent, making it the recommended metric for comparing across series. For probabilistic forecasts, use pinball loss or Continuous Ranked Probability Score (CRPS).

A common pitfall is using the same metric for model selection and final evaluation, which leads to overfitting to the metric. Instead, use a validation set (e.g., last 20% of data) for hyperparameter tuning, then evaluate on a holdout test set (e.g., last 10%) only once. Never look at the test set until the final model is chosen. In production, implement a backtesting pipeline that runs walk-forward validation daily or weekly to detect model degradation. Monitor forecast error distributions over time; a sudden increase in RMSE or a shift in bias (mean error) signals concept drift. Always compute confidence intervals for your metrics via bootstrapping over the forecast steps.

import numpy as np
import pandas as pd
from sklearn.metrics import mean_squared_error, mean_absolute_error
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# Generate synthetic data
np.random.seed(42)
dates = pd.date_range('2020-01-01', periods=200, freq='M')
series = np.linspace(100, 200, 200) + 10 * np.sin(2 * np.pi * np.arange(200) / 12) + np.random.normal(0, 5, 200)

# Walk-forward validation
n_splits = 5
horizon = 12  # forecast 12 months ahead
train_size = len(series) - n_splits * horizon

rmse_list = []
mae_list = []

for i in range(n_splits):
    train_end = train_size + i * horizon
    train = series[:train_end]
    test = series[train_end:train_end + horizon]
    
    model = ExponentialSmoothing(train, trend='add', seasonal='add', seasonal_periods=12)
    fit = model.fit()
    forecast = fit.forecast(steps=horizon)
    
    rmse = np.sqrt(mean_squared_error(test, forecast))
    mae = mean_absolute_error(test, forecast)
    rmse_list.append(rmse)
    mae_list.append(mae)
    print(f"Fold {i+1}: RMSE={rmse:.2f}, MAE={mae:.2f}")

print(f"\nAverage RMSE: {np.mean(rmse_list):.2f} +/- {np.std(rmse_list):.2f}")
print(f"Average MAE: {np.mean(mae_list):.2f} +/- {np.std(mae_list):.2f}")

Deep Learning for Time Series: LSTM, CNN, and Transformer Architectures

Deep learning has fundamentally changed time series forecasting by eliminating manual feature engineering and capturing complex, non-linear dependencies. Long Short-Term Memory (LSTM) networks, with their gated cell state and forget/input/output gates, directly address the vanishing gradient problem in vanilla RNNs. An LSTM cell maintains a memory vector c_t and hidden state h_t: f_t = σ(W_f·[h_{t-1}, x_t] + b_f), i_t = σ(W_i·[h_{t-1}, x_t] + b_i), o_t = σ(W_o·[h_{t-1}, x_t] + b_o), c_t = f_t ⊙ c_{t-1} + i_t ⊙ tanh(W_c·[h_{t-1}, x_t] + b_c), h_t = o_t ⊙ tanh(c_t). For univariate forecasting with lookback window L, an LSTM with 64–128 hidden units typically outperforms ARIMA on series with >5000 points and non-stationary variance. However, LSTMs are sequential by nature—they cannot parallelize over the time dimension during training, making them slower than CNNs or Transformers for long sequences.

Convolutional Neural Networks (CNNs) for time series use 1D dilated convolutions to expand receptive fields exponentially without stacking many layers. A WaveNet-style architecture stacks residual blocks with dilation factors 1, 2, 4, ..., 2^k, achieving receptive field size 2^{k+1} - 1. Causal padding ensures no leakage from future time steps. CNNs train faster than LSTMs (often 2–5x on GPU) and are more stable, but they lack an explicit memory mechanism—they rely on stacking many layers to capture long-range dependencies. For series with strong local patterns (e.g., hourly electricity load), a CNN with 8–16 filters and kernel size 3–5 often matches LSTM accuracy at lower latency.

Transformer architectures, originally from NLP, treat time steps as tokens and use self-attention to model all pairwise interactions. The scaled dot-product attention: Attention(Q,K,V) = softmax(QK^T/√d_k)V. For a sequence of length T, self-attention has O(T^2) complexity, which is prohibitive for long series (e.g., 10k+ steps). Variants like Informer (ProbSparse attention) and Autoformer (series decomposition + auto-correlation) reduce complexity to O(T log T) or O(T). In practice, Transformers excel when there are complex, long-range dependencies (e.g., multivariate retail demand with promotions, holidays, weather). But they require large datasets (100k+ time steps) and careful hyperparameter tuning—learning rate warmup, dropout (0.1–0.3), and label smoothing. For most production forecasting with <10k points, a well-tuned LSTM or TCN (temporal convolutional network) is simpler and more robust.

import numpy as np
import torch
import torch.nn as nn

class LSTMForecaster(nn.Module):
    def __init__(self, input_size=1, hidden_size=64, num_layers=2, output_size=1):
        super().__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, output_size)

    def forward(self, x):
        # x: (batch, seq_len, input_size)
        out, (h_n, c_n) = self.lstm(x)
        # take last hidden state
        return self.fc(out[:, -1, :])

# Generate synthetic sine wave with noise
np.random.seed(42)
t = np.linspace(0, 100, 1000)
data = np.sin(0.1 * t) + 0.1 * np.random.randn(1000)

# Create sequences: lookback=20, predict 1 step ahead
def create_sequences(data, seq_len=20):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len])
    return np.array(X).reshape(-1, seq_len, 1), np.array(y)

X, y = create_sequences(data)
split = int(0.8 * len(X))
X_train, y_train = torch.FloatTensor(X[:split]), torch.FloatTensor(y[:split])
X_test, y_test = torch.FloatTensor(X[split:]), torch.FloatTensor(y[split:])

model = LSTMForecaster()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
criterion = nn.MSELoss()

for epoch in range(50):
    model.train()
    optimizer.zero_grad()
    output = model(X_train)
    loss = criterion(output.squeeze(), y_train)
    loss.backward()
    torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
    optimizer.step()
    if epoch % 10 == 0:
        print(f'Epoch {epoch}, Loss: {loss.item():.6f}')

model.eval()
with torch.no_grad():
    preds = model(X_test).squeeze().numpy()
mae = np.mean(np.abs(preds - y_test.numpy()))
print(f'Test MAE: {mae:.4f}')

Handling Multiple Series: Global vs Local Models and Hierarchical Forecasting

When forecasting hundreds or thousands of related time series (e.g., SKU-level demand across stores), you face a choice: train one model per series (local) or a single model across all series (global). Local models (e.g., individual ARIMA or Prophet per SKU) are simple to implement and isolate failures, but they ignore cross-series patterns and scale poorly—training 10,000 models sequentially can take hours. Global models (e.g., a single LSTM or LightGBM trained on all series with a series ID feature) share statistical strength across series, especially beneficial for cold-start items with sparse history. A 2020 study by Salinas et al. (DeepAR) showed global RNNs reduce RMSE by 15–30% compared to local models on retail data with >1000 series. The trade-off: global models can suffer from negative transfer if series have fundamentally different dynamics (e.g., fast-moving vs slow-moving SKUs). A practical middle ground is clustering series into groups (e.g., by category, velocity, or seasonality pattern) and training a global model per cluster.

Hierarchical forecasting addresses the need for coherent predictions across aggregation levels—e.g., total company demand, regional, store, SKU. The hierarchy forms a tree: bottom-level series (SKU-store) sum to higher levels. Naively forecasting each level independently leads to incoherence: the sum of bottom forecasts ≠ top forecast. Reconciliation methods adjust forecasts to enforce consistency. The optimal combination approach (Wickramasuriya et al., 2019) uses a weighted least squares solution: ỹ = S(S' W^{-1} S)^{-1} S' W^{-1} ŷ, where S is the summing matrix, ŷ is the vector of base forecasts, and W is the covariance matrix of forecast errors. In practice, using W = diag(ŷ) (variance scaling) works well and avoids estimating full covariance. MinT (Minimum Trace) estimator further improves by shrinking the covariance. For production, implement reconciliation as a post-processing step—it's O(n^3) in the number of bottom series, so for >10k bottom series, use iterative or sparse methods (e.g., hts R package or PySpark-based reconciliation).

A critical nuance: global models can be combined with hierarchical reconciliation. Train a global model at the bottom level (SKU-store), then reconcile up. This captures granular patterns while ensuring top-level forecasts are coherent. Alternatively, train separate global models at each level and reconcile—this is more flexible but increases maintenance. In practice, the bottom-up approach (forecast bottom, sum up) is simplest and often performs within 1–2% of optimal reconciliation, especially when bottom-level models are accurate. For sparse hierarchies (e.g., many zero-demand days), use a two-stage model: first predict demand occurrence (binary classifier), then predict magnitude conditional on occurrence.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error

# Simulate 50 time series with shared trend + individual noise
np.random.seed(42)
n_series = 50
n_timesteps = 200
t = np.arange(n_timesteps)

series_data = []
for i in range(n_series):
    trend = 0.01 * t
    seasonality = 5 * np.sin(2 * np.pi * t / 12)
    noise = np.random.randn(n_timesteps) * 2
    series_data.append(trend + seasonality + noise + i * 0.1)  # series-specific offset

# Local model: fit LinearRegression per series (using last 12 lags)
def create_lags(series, lags=12):
    X, y = [], []
    for i in range(lags, len(series)):
        X.append(series[i-lags:i])
        y.append(series[i])
    return np.array(X), np.array(y)

local_maes = []
for s in series_data:
    X, y = create_lags(s)
    split = int(0.8 * len(X))
    model = LinearRegression()
    model.fit(X[:split], y[:split])
    preds = model.predict(X[split:])
    local_maes.append(mean_absolute_error(y[split:], preds))
print(f'Local model avg MAE: {np.mean(local_maes):.3f}')

# Global model: stack all series, add series_id as feature
X_global, y_global = [], []
for idx, s in enumerate(series_data):
    X_lags, y_lags = create_lags(s)
    # Add one-hot series ID (simplified: use scalar ID)
    id_feat = np.full((X_lags.shape[0], 1), idx)
    X_global.append(np.hstack([X_lags, id_feat]))
    y_global.append(y_lags)
X_global = np.vstack(X_global)
y_global = np.concatenate(y_global)

# Train/test split per series (maintain temporal order)
split_idx = int(0.8 * n_timesteps) - 12  # account for lags
X_train = np.vstack([X_global[i*split_idx:(i+1)*split_idx] for i in range(n_series)])
y_train = np.concatenate([y_global[i*split_idx:(i+1)*split_idx] for i in range(n_series)])
X_test = np.vstack([X_global[(i+1)*split_idx:] for i in range(n_series)])
y_test = np.concatenate([y_global[(i+1)*split_idx:] for i in range(n_series)])

global_model = LinearRegression()
global_model.fit(X_train, y_train)
global_preds = global_model.predict(X_test)
global_mae = mean_absolute_error(y_test, global_preds)
print(f'Global model MAE: {global_mae:.3f}')

Production Deployment: Retraining, Monitoring Drift, and Fallback Strategies

Deploying a forecasting model is not a one-time event—it's a continuous cycle. The first decision is retraining frequency: fixed schedule (daily, weekly) vs trigger-based (when drift is detected). For retail demand, daily retraining with a rolling window (e.g., last 365 days) is common. The retraining pipeline must be idempotent and versioned: use a DAG orchestrator (Airflow, Prefect) to fetch fresh data, validate schema, train model, evaluate on holdout, and push to registry. A typical pipeline takes 10–30 minutes for 1000 SKUs with LightGBM. For deep learning models, consider incremental training (warm-start from previous weights) to reduce time by 50–70%. But beware of catastrophic forgetting—always validate on a fixed test set from the last N days.

Monitoring drift is critical. Two types: data drift (input distribution changes) and concept drift (relationship between inputs and target changes). For time series, use a sliding window of prediction residuals: compute mean absolute error (MAE) over the last 7 days and compare to a baseline (e.g., 95th percentile of historical MAE). If MAE exceeds threshold, trigger alert and retrain. More sophisticated: use a two-sample Kolmogorov-Smirnov test on feature distributions (e.g., lag-1 values) between training and recent windows. For concept drift, monitor the cumulative sum (CUSUM) of signed errors: S_t = max(0, S_{t-1} + e_t - k), where e_t is the error and k is a slack parameter (typically 0.5 std of errors). If S_t exceeds threshold h (e.g., 5 std), flag drift. In practice, a simple rule like "retrain if MAE > 1.5x baseline for 3 consecutive days" catches 90% of drifts with low false positives.

Fallback strategies are essential. When the primary model fails (e.g., NaN predictions, API timeout, drift alert), you need a degraded mode. Common fallbacks: (1) Last known good forecast—cache the previous forecast and use it for up to 7 days. (2) Simple statistical model—fit a naive seasonal model (e.g., same day last week) as a lightweight backup. (3) Ensemble fallback—average of last 3 model versions from registry. Implement a circuit breaker pattern: if primary model errors >5% in a 10-minute window, switch to fallback for 1 hour, then retry. Log all fallback activations for post-mortem. For critical applications (e.g., energy grid balancing), use a redundant deployment across two cloud regions with automatic failover.

import numpy as np
from scipy.stats import ks_2samp

class DriftMonitor:
    def __init__(self, baseline_mae, threshold_factor=1.5, consecutive_days=3):
        self.baseline_mae = baseline_mae
        self.threshold = baseline_mae * threshold_factor
        self.consecutive_days = consecutive_days
        self.recent_errors = []
        self.drift_count = 0

    def update(self, y_true, y_pred):
        mae = np.mean(np.abs(y_true - y_pred))
        self.recent_errors.append(mae)
        if len(self.recent_errors) > 7:
            self.recent_errors.pop(0)
        # Check if last 3 days exceed threshold
        if len(self.recent_errors) >= self.consecutive_days:
            if all(e > self.threshold for e in self.recent_errors[-3:]):
                self.drift_count += 1
                return True
        return False

    def ks_drift(self, train_features, recent_features):
        # Compare distribution of first lag feature
        stat, p_value = ks_2samp(train_features, recent_features)
        return p_value < 0.05  # reject null if significant drift

# Simulate
np.random.seed(42)
baseline_mae = 0.5
monitor = DriftMonitor(baseline_mae)

# Normal period
for _ in range(10):
    y_true = np.random.randn(100)
    y_pred = y_true + np.random.randn(100) * 0.4  # MAE ~0.5
    drift = monitor.update(y_true, y_pred)
    print(f'Normal: drift={drift}')

# Drift period (sudden shift)
for _ in range(5):
    y_true = np.random.randn(100) + 2.0  # shift
    y_pred = y_true + np.random.randn(100) * 0.4  # model still predicts old distribution
    drift = monitor.update(y_true, y_pred)
    print(f'Drift: drift={drift}')

Case Study: Building a Real-Time Demand Forecasting Pipeline

Consider a grocery chain with 500 stores, each selling 10,000 SKUs. They need hourly demand forecasts for the next 48 hours to optimize inventory replenishment and reduce waste. The pipeline must handle 5 million time series (store-SKU-hour) with updates every hour. The architecture: data ingestion via Kafka streams (point-of-sale transactions, weather, promotions), feature engineering in Spark Structured Streaming (compute rolling 7-day averages, holiday flags, price elasticity), model inference with a global LightGBM model (trained on 6 months of data, 200 features), and post-processing for hierarchical reconciliation (store → region → total).

The model: LightGBM with 500 trees, max_depth=8, learning_rate=0.05, trained on 100M rows (sampled 20% of store-SKU combinations). Features include: hour-of-day, day-of-week, month, lagged demand (1, 2, 3, 7, 14, 28 days), rolling mean (7, 28 days), rolling std (7 days), price discount ratio, temperature, and holiday proximity. Training takes 4 hours on a 32-core machine with 256GB RAM. Inference is sub-100ms per batch of 10k series using LightGBM's predict() in C API via Python bindings. To handle 5M series per hour, we batch them into 500 chunks of 10k, run on 50 workers (10 chunks each) in a Spark cluster, achieving 3-minute total inference time.

Reconciliation: bottom-level forecasts (store-SKU-hour) are summed to store, region, and total. We use bottom-up reconciliation because it's simple and the bottom-level model is accurate (MAPE 12%). For the top-level (total company), we also run a separate global model as a sanity check—if the sum of bottom forecasts deviates >5% from the top model, we alert. The pipeline runs every hour on the hour, with a 10-minute SLA. Fallback: if LightGBM inference fails, we serve the previous hour's forecast with a decay factor (0.95 per hour) for up to 6 hours. If Kafka stream is down, we use the last 7 days' average for that hour. All predictions are written to a PostgreSQL database and served via a REST API (FastAPI) with Redis caching (TTL=1 hour).

Results: After deployment, the chain reduced out-of-stock incidents by 18% and perishable waste by 22%. The key lessons: (1) Global model with LightGBM was 3x faster to train and 10x faster to infer than an LSTM alternative, with comparable accuracy (MAPE 12% vs 11.5%). (2) Feature engineering mattered more than model architecture—adding weather and promotion features reduced MAPE by 4 percentage points. (3) Monitoring drift on top-100 SKUs (by revenue) caught 80% of issues with 1% of compute. (4) The fallback strategy was activated 12 times in 6 months, each time preventing a forecast outage. The pipeline cost $500/month in cloud compute, saving $200k/month in waste reduction.

import lightgbm as lgb
import pandas as pd
import numpy as np
from datetime import datetime, timedelta

# Simulate a batch of 10k store-SKU-hour forecasts
np.random.seed(42)
n_series = 10000

# Generate features (simplified)
df = pd.DataFrame({
    'hour': np.random.randint(0, 24, n_series),
    'day_of_week': np.random.randint(0, 7, n_series),
    'lag_1': np.random.randn(n_series) * 10 + 50,
    'lag_7': np.random.randn(n_series) * 10 + 50,
    'rolling_mean_7': np.random.randn(n_series) * 5 + 50,
    'price_discount': np.random.uniform(0, 0.3, n_series),
    'temperature': np.random.uniform(-5, 35, n_series),
    'is_holiday': np.random.randint(0, 2, n_series),
})

# Load pre-trained model (simulate with random forest)
# In production, load from MLflow: model = mlflow.lightgbm.load_model('models:/demand_model/Production')
# For demo, create a dummy model
params = {
    'objective': 'regression',
    'metric': 'mae',
    'num_leaves': 31,
    'learning_rate': 0.05,
    'n_estimators': 100,
    'verbose': -1
}
# Train on fake data for demo
X_train = df.values
y_train = np.random.randn(n_series) * 10 + 50
model = lgb.LGBMRegressor(**params)
model.fit(X_train, y_train)

# Predict
preds = model.predict(df)
df['forecast'] = preds

# Bottom-up reconciliation: sum to store level (simulate 100 stores)
df['store_id'] = np.random.randint(0, 100, n_series)
store_forecast = df.groupby('store_id')['forecast'].sum().reset_index()
store_forecast.columns = ['store_id', 'forecast_store']

print(f'Total forecast (sum of stores): {store_forecast["forecast_store"].sum():.2f}')
print(f'Average store forecast: {store_forecast["forecast_store"].mean():.2f}')
print(f'Inference time for 10k series: ~100ms (simulated)')