Time Series Forecasting in Production: Feature Engineering, Model Selection, and Deployment Patterns for Real-Time Predictions

Most time series forecasting projects fail in the same way. The model scores well on a held-out test set, ships to production, and quietly degrades over the next three months until a business analyst notices the predictions are useless. The test set was constructed incorrectly. The features leak future information. The model was never retrained after the initial deploy. Monitoring checks accuracy on a sample that is already stale.

This article covers the parts that cause those failures: how to engineer features that do not leak, how to choose between model families based on actual tradeoffs rather than benchmarks, how to validate without data contamination, and how to structure the serving layer so the model stays maintainable after the first deploy.

Feature Engineering

The most common source of silent failure in time series models is feature leakage. When you create lag features or rolling statistics incorrectly, you inadvertently give the model access to future information during training. The model learns a spurious correlation, tests perfectly, and breaks the moment you remove the future signal in production.

The rule is simple: every feature computed at time t must use only values from times strictly before t.

import pandas as pd
import numpy as np

def build_lag_features(df: pd.DataFrame, target_col: str, lags: list[int]) -> pd.DataFrame:
    """
    Builds lag features with correct temporal alignment.
    df must be sorted ascending by timestamp before calling this.
    """
    df = df.sort_values("timestamp").copy()
    for lag in lags:
        df[f"{target_col}_lag_{lag}"] = df[target_col].shift(lag)
    return df

def build_rolling_features(
    df: pd.DataFrame,
    target_col: str,
    windows: list[int],
    min_periods: int = 1
) -> pd.DataFrame:
    """
    Rolling statistics. The `shift(1)` is critical: it prevents the current
    observation from being included in its own rolling window.
    """
    df = df.sort_values("timestamp").copy()
    for window in windows:
        rolled = df[target_col].shift(1).rolling(window=window, min_periods=min_periods)
        df[f"{target_col}_rolling_mean_{window}"] = rolled.mean()
        df[f"{target_col}_rolling_std_{window}"] = rolled.std()
        df[f"{target_col}_rolling_min_{window}"] = rolled.min()
        df[f"{target_col}_rolling_max_{window}"] = rolled.max()
    return df

def build_calendar_features(df: pd.DataFrame, timestamp_col: str) -> pd.DataFrame:
    ts = df[timestamp_col]
    df["hour"] = ts.dt.hour
    df["day_of_week"] = ts.dt.dayofweek
    df["day_of_month"] = ts.dt.day
    df["month"] = ts.dt.month
    df["is_weekend"] = (ts.dt.dayofweek >= 5).astype(int)
    df["quarter"] = ts.dt.quarter
    return df

def build_fourier_features(
    df: pd.DataFrame,
    timestamp_col: str,
    period: float,
    n_terms: int,
    feature_prefix: str
) -> pd.DataFrame:
    """
    Fourier terms capture smooth seasonality better than one-hot calendar features
    for periodic signals (weekly: period=7, annual: period=365.25).
    """
    t = (df[timestamp_col] - df[timestamp_col].min()).dt.total_seconds()
    period_seconds = period * 86400
    for k in range(1, n_terms + 1):
        df[f"{feature_prefix}_sin_{k}"] = np.sin(2 * np.pi * k * t / period_seconds)
        df[f"{feature_prefix}_cos_{k}"] = np.cos(2 * np.pi * k * t / period_seconds)
    return df

Fourier terms deserve a specific mention. Calendar features like day_of_week force the model to learn a step function for seasonality. If you have daily data with weekly seasonality, Fourier terms with period=7 and two to four terms will model the smooth seasonal curve more accurately, with fewer parameters. Use both: calendar features for sharp discontinuities (holidays, shift changes), Fourier terms for smooth periodic patterns.

For lag selection: start with the natural periodicity of your series. Weekly data with annual seasonality needs lags at 1, 2, 4, 8, 26, and 52. Adding lags beyond the forecast horizon is not useful and increases overfitting risk. If you are forecasting 7 days ahead, your shortest lag should be at least 7.

Walk-Forward Validation

Random train-test splits do not work for time series. They allow the model to train on data from the future relative to the test period, which inflates every accuracy metric. Walk-forward (expanding window) validation is the correct approach.

from sklearn.metrics import mean_absolute_error, mean_squared_error
from typing import Iterator
import numpy as np

def walk_forward_splits(
    n_samples: int,
    min_train_size: int,
    test_size: int,
    step_size: int
) -> Iterator[tuple[np.ndarray, np.ndarray]]:
    """
    Yields (train_indices, test_indices) for each fold.
    train grows with each fold; test window slides forward.
    """
    start = min_train_size
    while start + test_size <= n_samples:
        train_idx = np.arange(0, start)
        test_idx = np.arange(start, start + test_size)
        yield train_idx, test_idx
        start += step_size

def evaluate_walk_forward(model, X: np.ndarray, y: np.ndarray, config: dict) -> dict:
    maes, rmses = [], []
    for train_idx, test_idx in walk_forward_splits(
        n_samples=len(X),
        min_train_size=config["min_train_size"],
        test_size=config["test_size"],
        step_size=config["step_size"]
    ):
        model.fit(X[train_idx], y[train_idx])
        preds = model.predict(X[test_idx])
        maes.append(mean_absolute_error(y[test_idx], preds))
        rmses.append(np.sqrt(mean_squared_error(y[test_idx], preds)))
    return {
        "mae_mean": float(np.mean(maes)),
        "mae_std": float(np.std(maes)),
        "rmse_mean": float(np.mean(rmses)),
        "rmse_std": float(np.std(rmses)),
        "n_folds": len(maes)
    }

The mae_std is as important as mae_mean. High variance across folds means the model’s performance depends heavily on which time window it sees during training. That instability will show up in production as inconsistent forecast quality across different periods of the year.

A common mistake is using too small a min_train_size. The model needs enough history to learn seasonal patterns. If your series has annual seasonality, min_train_size should cover at least two to three full cycles.

Model Selection Tradeoffs

No model family dominates across all forecasting problems. The right choice depends on your data volume, update frequency, interpretability requirements, and how much engineering effort you can spend on the serving layer.

Model Family	Best For	Weaknesses	Serving Complexity
ARIMA / SARIMA	Single series, short history, interpretability required	Manual order selection, does not scale to many series	Low (stateless, fast inference)
Prophet	Business time series with holidays, non-technical users	Slow on many series, weak on non-seasonal data	Low
XGBoost / LightGBM	Many series, exogenous features, tabular patterns	Requires careful feature engineering, no native uncertainty	Medium
N-BEATS	Pure time series, no exogenous features needed	Training cost, harder to debug	Medium
Temporal Fusion Transformer (TFT)	Mixed covariates, multiple related series, uncertainty needed	Training cost, complex architecture	High
Exponential Smoothing (ETS)	Baseline, simple seasonal patterns	Limited expressivity	Very Low

For most production forecasting problems at mid-scale (hundreds to thousands of series, hourly or daily cadence), LightGBM with well-engineered features outperforms all the others on the time-cost tradeoff. It trains fast, serves fast, is debuggable, handles missing values, and integrates naturally with a tabular feature pipeline.

Reserve deep learning models (N-BEATS, TFT) for cases where you have dense training data (years of sub-hourly observations), the relationship between series is important, or you need calibrated prediction intervals. The operational overhead is not worth it for a demand forecasting problem with daily data going back eighteen months.

import lightgbm as lgb
from dataclasses import dataclass

@dataclass
class ForecastingConfig:
    lags: list[int]
    rolling_windows: list[int]
    fourier_period_days: float
    fourier_n_terms: int
    forecast_horizon: int
    lgb_params: dict

DEFAULT_CONFIG = ForecastingConfig(
    lags=[1, 2, 3, 7, 14, 21, 28],
    rolling_windows=[7, 14, 28],
    fourier_period_days=7.0,
    fourier_n_terms=3,
    forecast_horizon=7,
    lgb_params={
        "objective": "regression",
        "metric": "mae",
        "num_leaves": 63,
        "learning_rate": 0.05,
        "feature_fraction": 0.8,
        "bagging_fraction": 0.8,
        "bagging_freq": 5,
        "min_child_samples": 20,
        "n_estimators": 500,
        "early_stopping_rounds": 50,
        "verbosity": -1
    }
)

def train_lgbm_forecaster(
    X_train: np.ndarray,
    y_train: np.ndarray,
    X_val: np.ndarray,
    y_val: np.ndarray,
    config: ForecastingConfig
) -> lgb.LGBMRegressor:
    model = lgb.LGBMRegressor(**config.lgb_params)
    model.fit(
        X_train, y_train,
        eval_set=[(X_val, y_val)],
        callbacks=[lgb.early_stopping(config.lgb_params["early_stopping_rounds"], verbose=False)]
    )
    return model

One pattern that works well at scale is recursive multi-step forecasting: train a single one-step-ahead model, then feed its predictions as lagged inputs for subsequent steps. It is simpler than training a separate model per horizon and degrades gracefully as the horizon grows. The alternative (direct multi-output forecasting) requires a separate model per step and increases operational surface area.

Deployment Patterns

The serving layer is where most of the operational complexity lives. Two decisions drive everything: batch versus online, and how to handle prediction caching.

Batch forecasting runs on a schedule (hourly, daily), writes predictions to a store, and serves reads from that store. It is the right pattern when forecast consumers can tolerate some staleness and when the feature pipeline runs on historical data. Computational cost is predictable. Latency at read time is near zero.

Online forecasting computes predictions on demand, using a feature store that holds precomputed lag and rolling features. It is necessary when consumers need predictions at arbitrary future times or when the model must reflect the most recent observations. Latency is higher and the feature computation path must be fast.

Here is a TypeScript serving layer that handles both patterns, reading from a prediction cache for batch forecasts and falling through to the model service for online requests:

interface PredictionRequest {
  seriesId: string;
  forecastHorizon: number;
  asOf?: Date;
}

interface PredictionResult {
  seriesId: string;
  predictions: Array<{ timestamp: Date; value: number; lower?: number; upper?: number }>;
  modelVersion: string;
  generatedAt: Date;
  source: "cache" | "online";
}

interface PredictionStore {
  get(seriesId: string, horizon: number): Promise<PredictionResult | null>;
}

interface ModelService {
  predict(request: PredictionRequest): Promise<PredictionResult>;
}

class ForecastingService {
  private readonly cacheTtlMs: number;

  constructor(
    private readonly store: PredictionStore,
    private readonly model: ModelService,
    cacheTtlSeconds: number = 3600
  ) {
    this.cacheTtlMs = cacheTtlSeconds * 1000;
  }

  async forecast(request: PredictionRequest): Promise<PredictionResult> {
    const cached = await this.store.get(request.seriesId, request.forecastHorizon);

    if (cached !== null) {
      const age = Date.now() - cached.generatedAt.getTime();
      if (age < this.cacheTtlMs) {
        return { ...cached, source: "cache" };
      }
    }

    const result = await this.model.predict(request);
    return { ...result, source: "online" };
  }
}

// Health check endpoint pattern for the model service
interface ForecastHealthStatus {
  modelVersion: string;
  lastTrainedAt: Date;
  featurePipelineLastRunAt: Date;
  stalePredictionCount: number;
  status: "healthy" | "degraded" | "unhealthy";
}

async function checkForecastHealth(
  modelVersion: string,
  lastTrainedAt: Date,
  featurePipelineLastRunAt: Date,
  stalePredictionCount: number,
  maxStalenessHours: number = 25
): Promise<ForecastHealthStatus> {
  const pipelineAge = (Date.now() - featurePipelineLastRunAt.getTime()) / (1000 * 60 * 60);
  const modelAge = (Date.now() - lastTrainedAt.getTime()) / (1000 * 60 * 60 * 24);

  let status: ForecastHealthStatus["status"] = "healthy";
  if (pipelineAge > maxStalenessHours || stalePredictionCount > 0) {
    status = "degraded";
  }
  if (pipelineAge > maxStalenessHours * 2 || modelAge > 30) {
    status = "unhealthy";
  }

  return { modelVersion, lastTrainedAt, featurePipelineLastRunAt, stalePredictionCount, status };
}

The source field on every prediction response is not optional. Without it, you cannot distinguish a stale cache hit from a freshly computed prediction in your observability layer.

Production Considerations

Monitoring forecast degradation. MAE or RMSE on a rolling window of actuals vs predictions is the core metric. Track it per series segment: a model that performs well on average may be badly wrong for your high-value series. Alert when the rolling error exceeds 1.5x the baseline you measured during walk-forward validation.

Feature pipeline reliability. The model is only as good as its inputs. Build explicit freshness checks into the feature pipeline: if lag features are more than N periods stale (because upstream data was delayed), emit a staleness flag and optionally fall back to the previous cached prediction rather than serving predictions built on stale features. Serving a known-stale prediction is preferable to serving a silently wrong one.

Data drift vs concept drift. Data drift: the distribution of input features shifts (new customer segment added, geographic expansion). Concept drift: the relationship between features and target changes (demand pattern changes post-COVID, pricing algorithm changed). Both degrade forecast quality, but the remediation differs. For data drift, retraining on recent data usually recovers performance. Concept drift may require re-engineering features or reconsidering the model family.

from scipy import stats

def detect_feature_drift(
    reference: np.ndarray,
    current: np.ndarray,
    alpha: float = 0.05
) -> dict:
    """
    KS test for distributional drift on a single feature.
    Returns whether drift is detected and the statistic.
    Run this per feature after each retraining cycle.
    """
    stat, p_value = stats.ks_2samp(reference, current)
    return {
        "drift_detected": bool(p_value < alpha),
        "ks_statistic": float(stat),
        "p_value": float(p_value)
    }

Retraining cadence. For most business forecasting problems, weekly retraining with an expanding window strikes the right balance. Daily retraining adds operational cost without much accuracy gain unless the series has very high non-stationarity. Monthly is often too slow to adapt to structural changes. The expanding window approach (keeping all historical data) beats rolling window for series with long-range seasonal patterns because the model keeps learning the full seasonality structure.

Cold start for new series. When a new series is added (new product, new location), you have no lag history to build features from. Three options in order of operational simplicity: (1) fall back to a statistical model like ETS that needs no lagged features, (2) find a similar series and use its history as a proxy, or (3) use a global model trained across all series that can transfer patterns. Option three (global models, particularly LightGBM trained on all series simultaneously with a series_id feature) is the most scalable approach if you have enough series.

Prediction intervals. Point forecasts are rarely sufficient for production use cases. Downstream systems that react to forecasts (inventory replenishment, capacity planning, staffing) need uncertainty bounds. For LightGBM, quantile regression (set objective="quantile" and train separate models for the 10th and 90th percentiles) is the practical approach. It is more honest than multiplying the RMSE by a fixed factor.

The Gap Between Offline and Online

The failure mode this article opened with comes down to a single gap: the offline evaluation environment does not match the online serving environment. The test set was built incorrectly, or the feature computation path in training differed from the path in serving.

The fix is to treat the feature pipeline as the first-class artifact, not the model. The same code that builds features during training should run during serving, with the only difference being whether it reads from a historical dataset or a live feature store. When that is true, a model that validates well offline will generally hold up in production, and degradation becomes diagnosable rather than mysterious.

Time series forecasting is not an ML problem that gets solved once. It is an ongoing engineering problem: the world changes, data pipelines break, feature distributions shift, and models need retraining. Build the monitoring and retraining infrastructure before you optimize the model, not after.