Demand Forecasting for FMCG: Time Series, Promotions, Seasonality, and Supply Chain Integration

Demand forecasting in Fast-Moving Consumer Goods (FMCG) is one of those problems where the theoretical elegance of time series modeling collides with the messy reality of supply chains. After building forecasting systems for multiple FMCG companies — from regional food manufacturers to multinational consumer goods corporations — I’ve learned that the forecast model itself accounts for maybe 30% of the solution. The other 70% is data pipeline engineering, business process integration, and managing stakeholder expectations.

The challenge is straightforward to state: predict how much of each product will be sold at each location over the next 1-12 weeks. The difficulty lies in the scale (millions of SKU-location combinations), the noise (retail demand is inherently volatile), the external factors (promotions, holidays, weather, competitor actions), and the hierarchical consistency requirements (your forecasts need to add up correctly from store to region to national level).

This post covers the technical approach that actually works in production, the pitfalls that catch first-time forecasters, and the integration patterns that make forecasts useful to supply chain teams.

The Data Landscape

FMCG demand forecasting draws on multiple data sources, each with its own quirks:

Sales Data

The foundation. But “sales” is not the same as “demand”:

import pandas as pd
import numpy as np

def clean_sales_data(sales: pd.DataFrame) -> pd.DataFrame:
    """Clean sales data to approximate true demand."""
    # Remove stockout periods — zero sales during stockout don't mean zero demand
    sales["is_stockout"] = sales["inventory_on_hand"] == 0
    sales["adjusted_sales"] = sales["units_sold"].where(
        ~sales["is_stockout"],
        sales["units_sold"].rolling(4, min_periods=1).mean(),  # Impute with recent average
    )

    # Remove outliers (likely data entry errors)
    q99 = sales["adjusted_sales"].quantile(0.99)
    sales["adjusted_sales"] = sales["adjusted_sales"].clip(upper=q99 * 3)

    # Account for returns
    sales["net_sales"] = sales["adjusted_sales"] - sales["units_returned"].fillna(0)

    return sales

Stockout adjustment is the single most impactful preprocessing step in FMCG forecasting. If you forecast on raw sales data, you’re predicting sell-through, not demand. During stockouts, demand is censored — you don’t observe the sales that would have happened. Ignoring this creates biased forecasts that systematically underestimate demand for popular products.

Promotion Data

Promotions account for 20-40% of FMCG sales and are the largest source of demand variability:

def engineer_promotion_features(
    sales: pd.DataFrame,
    promotions: pd.DataFrame,
    lookback_weeks: int = 52,
) -> pd.DataFrame:
    """Engineer promotion-related features."""
    features = sales[["sku_id", "location_id", "week"]].copy()

    # Merge promotion data
    features = features.merge(
        promotions[["sku_id", "location_id", "week", "promo_type", "discount_pct", "display_type"]],
        on=["sku_id", "location_id", "week"],
        how="left",
    )

    # Current promotion features
    features["is_on_promo"] = features["promo_type"].notna().astype(int)
    features["discount_depth"] = features["discount_pct"].fillna(0)

    # Promotion history
    for lag in [1, 2, 3, 4, 8, 12]:
        features[f"was_on_promo_{lag}w_ago"] = (
            features.groupby(["sku_id", "location_id"])["is_on_promo"].shift(lag)
        )

    # Post-promotion dip indicator (demand often drops below baseline after promotion)
    features["post_promo"] = (
        (features["is_on_promo"] == 0)
        & (features["was_on_promo_1w_ago"] == 1)
    ).astype(int)

    # Promotion frequency (how often is this SKU promoted?)
    promo_frequency = features.groupby(["sku_id", "location_id"])["is_on_promo"].mean()
    features = features.merge(
        promo_frequency.rename("promo_frequency"),
        on=["sku_id", "location_id"],
        how="left",
    )

    return features

Post-promotion dips are real and forecastable. When a promotion ends, demand typically drops below the normal baseline for 1-2 weeks as consumers work through their pantry stock. A model that doesn’t account for this will overforecast the weeks following a promotion, leading to excess inventory.

Modeling Approaches

The Model Zoo Strategy

No single model dominates across all SKU-location combinations. The best-performing approach I’ve found is a model zoo — train multiple models and select the best one per series:

from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
import lightgbm as lgb
from prophet import Prophet

class ModelZoo:
    """Ensemble of forecasting models with automatic selection."""

    def __init__(self):
        self.models = {
            "ets": self._fit_ets,
            "sarimax": self._fit_sarimax,
            "prophet": self._fit_prophet,
            "lightgbm": self._fit_lightgbm,
        }

    def forecast(self, series: pd.Series, horizon: int, features: pd.DataFrame = None) -> dict:
        """Generate forecasts from all models and select best."""
        results = {}

        for name, fit_fn in self.models.items():
            try:
                result = fit_fn(series, horizon, features)
                results[name] = result
            except Exception as e:
                logger.warning(f"Model {name} failed: {e}")
                continue

        # Select best model based on cross-validated MAPE
        best_model = min(results, key=lambda k: results[k]["cv_mape"])
        return {
            "forecast": results[best_model]["forecast"],
            "model_used": best_model,
            "all_results": results,
        }

    def _fit_ets(self, series: pd.Series, horizon: int, features=None) -> dict:
        """Exponential Smoothing with trend and seasonality."""
        model = ExponentialSmoothing(
            series,
            trend="add",
            seasonal="mul",
            seasonal_periods=52,  # Weekly data, annual seasonality
        ).fit(optimized=True)

        forecast = model.forecast(horizon)
        cv_mape = self._cross_validate(series, "ets", horizon)

        return {"forecast": forecast, "cv_mape": cv_mape}

    def _fit_prophet(self, series: pd.Series, horizon: int, features=None) -> dict:
        """Facebook Prophet with custom seasonality."""
        df = pd.DataFrame({"ds": series.index, "y": series.values})

        model = Prophet(
            yearly_seasonality=True,
            weekly_seasonality=True,
            changepoint_prior_scale=0.05,
        )

        # Add promotion regressors if available
        if features is not None and "is_on_promo" in features.columns:
            model.add_regressor("is_on_promo")

        model.fit(df)
        future = model.make_future_dataframe(periods=horizon, freq="W")
        forecast = model.predict(future)

        cv_mape = self._cross_validate(series, "prophet", horizon)

        return {"forecast": forecast["yhat"].tail(horizon), "cv_mape": cv_mape}

    def _fit_lightgbm(self, series: pd.Series, horizon: int, features: pd.DataFrame) -> dict:
        """Gradient boosting with lag features."""
        # Create lag features
        df = pd.DataFrame({"y": series})
        for lag in [1, 2, 3, 4, 8, 12, 52]:
            df[f"lag_{lag}"] = df["y"].shift(lag)
        df["rolling_mean_4"] = df["y"].rolling(4).mean()
        df["rolling_std_4"] = df["y"].rolling(4).std()
        df = df.dropna()

        X = df.drop(columns=["y"])
        y = df["y"]

        model = lgb.LGBMRegressor(
            n_estimators=500,
            max_depth=6,
            learning_rate=0.05,
            subsample=0.8,
        )
        model.fit(X, y)

        # Recursive forecasting
        forecast = []
        last_values = series.tail(52).tolist()
        for _ in range(horizon):
            lag_features = self._create_lag_features(last_values)
            pred = model.predict([lag_features])[0]
            forecast.append(pred)
            last_values.append(pred)
            last_values.pop(0)

        cv_mape = self._cross_validate(series, "lightgbm", horizon)

        return {"forecast": pd.Series(forecast), "cv_mape": cv_mape}

ETS (Exponential Smoothing) is surprisingly competitive for FMCG data. It handles trend and seasonality naturally, it’s fast to fit, and it doesn’t require feature engineering. For many SKU-location combinations, especially those with limited history, ETS outperforms more complex models.

LightGBM for Forecasting

Gradient-boosted trees have become the dominant approach for demand forecasting in practice, not because they’re theoretically optimal for time series, but because they handle the messy reality of FMCG data:

def create_forecast_features(
    sales: pd.DataFrame,
    promotions: pd.DataFrame,
    calendar: pd.DataFrame,
) -> pd.DataFrame:
    """Create feature matrix for gradient boosting forecast model."""
    features = sales[["sku_id", "location_id", "week", "units_sold"]].copy()

    # Lag features
    for lag in [1, 2, 3, 4, 8, 12, 26, 52]:
        features[f"lag_{lag}"] = (
            features.groupby(["sku_id", "location_id"])["units_sold"].shift(lag)
        )

    # Rolling statistics
    for window in [4, 8, 12, 26]:
        features[f"rolling_mean_{window}"] = (
            features.groupby(["sku_id", "location_id"])["units_sold"]
            .transform(lambda x: x.shift(1).rolling(window).mean())
        )
        features[f"rolling_std_{window}"] = (
            features.groupby(["sku_id", "location_id"])["units_sold"]
            .transform(lambda x: x.shift(1).rolling(window).std())
        )

    # Calendar features
    features = features.merge(calendar, on="week", how="left")

    # Promotion features
    features = features.merge(
        promotions[["sku_id", "location_id", "week", "is_on_promo", "discount_pct"]],
        on=["sku_id", "location_id", "week"],
        how="left",
    )

    # Price features
    features["price"] = sales["price"]
    features["price_change_pct"] = (
        features.groupby(["sku_id", "location_id"])["price"].pct_change()
    )

    return features

The key insight for LightGBM forecasting: create features that encode temporal patterns explicitly. Lags capture recent trends, rolling statistics capture local patterns, and calendar features capture seasonality. The model doesn’t know it’s looking at a time series — you have to give it temporal context through features.

Hierarchical Forecasting

FMCG companies need forecasts at multiple levels: SKU-store, SKU-region, SKU-national, category-region, total-company. These forecasts must be coherent — the sum of store-level forecasts must equal the regional forecast.

Bottom-Up vs. Top-Down vs. Optimal Reconciliation

import numpy as np

class HierarchicalForecaster:
    """Hierarchical forecast reconciliation."""

    def __init__(self, hierarchy: dict):
        """
        hierarchy: mapping of aggregation structure
        e.g., {"national": ["region_1", "region_2"],
               "region_1": ["store_1", "store_2"],
               "region_2": ["store_3", "store_4"]}
        """
        self.hierarchy = hierarchy

    def bottom_up(self, base_forecasts: dict) -> dict:
        """Sum leaf-level forecasts to get higher-level forecasts."""
        result = {}

        # Start from leaves
        for parent, children in self.hierarchy.items():
            result[parent] = sum(base_forecasts.get(child, 0) for child in children)

        # Add leaf forecasts
        for key, value in base_forecasts.items():
            if key not in self.hierarchy:
                result[key] = value

        return result

    def top_down(self, total_forecast: float, proportions: dict) -> dict:
        """Distribute top-level forecast using historical proportions."""
        result = {}
        for key, proportion in proportions.items():
            result[key] = total_forecast * proportion
        return result

    def optimal_reconciliation(
        self,
        base_forecasts: dict,
        forecast_errors: dict,
    ) -> dict:
        """
        MinT (Minimum Trace) reconciliation.
        Uses forecast error covariance to find optimal combination.
        """
        # This is the theoretically optimal approach
        # Implementation depends on the hierarchy structure
        # For large hierarchies, use the trace minimization approach
        from scipy.optimize import minimize

        n = len(base_forecasts)
        forecasts_array = np.array(list(base_forecasts.values()))

        # Error covariance matrix (simplified: diagonal)
        error_vars = np.array([forecast_errors.get(k, 1.0) for k in base_forecasts.keys()])
        W = np.diag(error_vars)

        # Summing matrix S (maps bottom-level to all levels)
        S = self._build_summing_matrix()

        # Optimal weights: G = (S'W^{-1}S)^{-1}S'W^{-1}
        W_inv = np.diag(1.0 / error_vars)
        G = np.linalg.solve(S.T @ W_inv @ S, S.T @ W_inv)

        # Reconciled forecasts
        reconciled = S @ G @ forecasts_array

        return dict(zip(self._get_all_labels(), reconciled))

Optimal reconciliation (MinT) consistently outperforms bottom-up and top-down approaches. Bottom-up forecasts are noisy at the leaf level but unbiased. Top-down forecasts are smooth but biased. MinT finds the compromise that minimizes overall forecast variance across the hierarchy.

Evaluation Metrics

Beyond MAPE

MAPE (Mean Absolute Percentage Error) is the most commonly used forecast metric, but it has serious problems:

def evaluate_forecast(actual: np.ndarray, forecast: np.ndarray) -> dict:
    """Comprehensive forecast evaluation metrics."""
    # MAPE — problematic when actuals are near zero
    non_zero_mask = actual > 0
    mape = np.mean(np.abs((actual[non_zero_mask] - forecast[non_zero_mask]) / actual[non_zero_mask])) * 100

    # WAPE — Weighted Absolute Percentage Error (handles zeros)
    wape = np.sum(np.abs(actual - forecast)) / np.sum(actual) * 100

    # MASE — Mean Absolute Scaled Error (scale-independent)
    naive_errors = np.abs(np.diff(actual))
    mase = np.mean(np.abs(actual - forecast)) / np.mean(naive_errors)

    # Bias — systematic over or under-forecasting
    bias = np.mean(forecast - actual) / np.mean(actual) * 100

    # Service level metrics — what supply chain actually cares about
    # Forecast accuracy at the quantile used for inventory decisions
    forecast_ratio = forecast / actual.clip(lower=1)
    fill_rate = np.mean(forecast >= actual) * 100  # % of time forecast covers actual demand

    return {
        "mape": mape,
        "wape": wape,
        "mase": mase,
        "bias_pct": bias,
        "fill_rate_pct": fill_rate,
    }

WAPE is the metric I recommend for FMCG forecasting. It handles zero demand periods, it’s interpretable (weighted average of percentage errors), and it correlates well with supply chain costs. MASE is better for comparing across series with different scales, but it’s harder to explain to business stakeholders.

Supply Chain Integration

From Forecast to Order

The forecast is only useful if it drives inventory and ordering decisions:

def forecast_to_order(
    forecast: dict,
    current_inventory: float,
    lead_time_days: int,
    service_level: float = 0.95,
) -> dict:
    """Convert demand forecast to order recommendation."""
    from scipy.stats import norm

    # Forecast for lead time + review period
    lt_demand = sum(f["mean"] for f in forecast[:lead_time_days // 7 + 1])
    lt_std = np.sqrt(sum(f["std"] ** 2 for f in forecast[:lead_time_days // 7 + 1]))

    # Safety stock for target service level
    z_score = norm.ppf(service_level)
    safety_stock = z_score * lt_std

    # Reorder point
    reorder_point = lt_demand + safety_stock

    # Order quantity (if below reorder point)
    if current_inventory < reorder_point:
        order_quantity = reorder_point - current_inventory + lt_demand
    else:
        order_quantity = 0

    return {
        "order_quantity": max(0, round(order_quantity)),
        "reorder_point": round(reorder_point),
        "safety_stock": round(safety_stock),
        "expected_demand": round(lt_demand),
    }

Forecast Value Added (FVA) Analysis

Track where forecast accuracy is gained or lost in the process:

def forecast_value_added(forecasts: pd.DataFrame) -> pd.DataFrame:
    """Analyze value added at each stage of forecasting."""
    stages = ["statistical_forecast", "analyst_adjustment", "consensus_forecast", "final_forecast"]

    fva_results = []
    for stage in stages:
        mape = evaluate_forecast(forecasts["actual"], forecasts[stage])["wape"]
        fva_results.append({"stage": stage, "wape": mape})

    # Calculate incremental value
    for i in range(1, len(fva_results)):
        fva_results[i]["incremental_value"] = (
            fva_results[i - 1]["wape"] - fva_results[i]["wape"]
        )

    return pd.DataFrame(fva_results)

In my experience, manual analyst adjustments improve forecasts only about 40% of the time. The other 60%, analysts add noise by overreacting to recent demand spikes or incorporating anecdotes that don’t represent systematic patterns. FVA analysis makes this visible and helps organizations calibrate the right level of human intervention.

Production Considerations

Scale Management

A typical FMCG company might have 10,000 SKUs × 500 locations = 5 million SKU-location series to forecast weekly. This demands efficient infrastructure:

import dask.dataframe as dd
from dask.distributed import Client

def parallel_forecast(
    sales_data: dd.DataFrame,
    forecast_fn,
    n_workers: int = 16,
) -> dd.DataFrame:
    """Parallel forecast across SKU-location combinations."""
    client = Client(n_workers=n_workers, threads_per_worker=2)

    # Group by SKU-location and forecast in parallel
    results = sales_data.groupby(["sku_id", "location_id"]).apply(
        forecast_fn,
        meta=pd.DataFrame({
            "sku_id": pd.Series(dtype="str"),
            "location_id": pd.Series(dtype="str"),
            "week": pd.Series(dtype="datetime64[ns]"),
            "forecast": pd.Series(dtype="float"),
        }),
    ).compute()

    return results

Forecast Refresh Cadence

Not all forecasts need the same refresh frequency:

• Fast-moving SKUs (top 20%): Weekly refresh
• Medium-moving SKUs (middle 60%): Biweekly refresh
• Slow-moving SKUs (bottom 20%): Monthly refresh with intermittent demand models

This tiered approach reduces compute costs by 40-60% while maintaining forecast accuracy where it matters most.

Lessons from the Field

1. Data quality trumps model complexity. Spend 80% of your effort on clean, stockout-adjusted, consistently coded data. The fanciest model can’t overcome garbage inputs.

2. Promotions are the hardest part. Promotion effects are non-linear, vary by product category, interact with each other, and create post-promotion dips. Budget significant time for promotion modeling.

3. Forecast accuracy is a lagging indicator. You won’t know if your forecast was good until the actual sales happen. Build fast feedback loops.

4. Bias is worse than noise. A forecast that’s consistently 10% too high is more damaging than one that fluctuates ±20% around the actual. Systematic bias leads to systematic overstock or stockouts.

5. The supply chain team is your customer. They don’t care about MAPE. They care about fill rates, inventory turns, and stockout frequency. Translate your accuracy metrics into supply chain outcomes.

Demand forecasting is one of those problems where the 80/20 rule applies aggressively. A simple model with clean data and good business process integration will outperform a sophisticated model that’s poorly integrated. Start simple, prove value, then add complexity where it demonstrably helps.