A two-year e-commerce inventory process built on Excel spreadsheets, intuition, and a forgotten Prophet model produced a 14% stockout rate on popular items and $2.3M in excess inventory sitting in warehouses. Replacing it with a demand forecasting system based on gradient boosting and real feature engineering dropped stockouts to 10.1% within three months -- a 28% reduction -- and brought excess inventory down by $890K.

This is the full technical walkthrough of that system: why classical time series methods fall short for complex demand patterns, how to engineer features that capture the real drivers of demand, how to train and validate gradient boosting models for time series, and how to deploy forecasts into a data warehouse where operations teams can use them.

Why ARIMA and Prophet Are Not Enough

ARIMA and Prophet are good tools for their intended use case -- univariate time series forecasting where the primary pattern is trend plus seasonality. For total website traffic or monthly revenue, they work well. But demand forecasting for individual products has characteristics that break these models:

• Exogenous variables dominate. A product's demand depends on promotions, competitor pricing, weather, holidays, and marketing spend -- not just its own history. ARIMA can include exogenous regressors (ARIMAX), but it struggles with dozens of them.
• Intermittent demand. Many SKUs sell zero units on most days, then spike. Classical time series models assume continuous demand patterns.
• Cross-product effects. When product A goes on sale, product B (a substitute) sees decreased demand. Time series models process each product independently.
• Non-linear interactions. A 20% discount during a holiday week has a very different effect than the same discount on a random Tuesday. Tree-based models capture these interactions naturally.

A concrete example from real sales data:

import pandas as pd
from prophet import Prophet
from sklearn.metrics import mean_absolute_error, mean_absolute_percentage_error

# Load sales data for a single popular SKU
sales = pd.read_parquet("data/daily_sales.parquet")
sku_data = sales[sales["sku_id"] == "ELEC-HDPHN-001"].copy()
sku_data = sku_data.rename(columns={"date": "ds", "units_sold": "y"})

# Prophet baseline
model = Prophet(yearly_seasonality=True, weekly_seasonality=True)
model.fit(sku_data[sku_data["ds"] < "2025-10-01"])

future = model.make_future_dataframe(periods=90)
forecast = model.predict(future)

# Evaluate on the held-out period
test_period = forecast[forecast["ds"] >= "2025-10-01"]
actuals = sku_data[sku_data["ds"] >= "2025-10-01"]
merged = test_period.merge(actuals, on="ds")

print(f"Prophet MAE:  {mean_absolute_error(merged['y'], merged['yhat']):.1f} units")
print(f"Prophet MAPE: {mean_absolute_percentage_error(merged['y'], merged['yhat']):.1%}")

# Output:
# Prophet MAE:  23.4 units
# Prophet MAPE: 31.2%

# The problem: Prophet missed the Black Friday spike entirely (predicted 85, actual 340)
# and over-predicted during a competitor's flash sale (predicted 95, actual 42)

A 31% MAPE means inventory decisions rest on forecasts that are wrong by almost a third on average. For a business running on thin margins, that is the difference between profitability and loss.

The ML Approach: Gradient Boosting on Tabular Features

Reframe demand forecasting as a tabular regression problem. Instead of treating the time series as a sequence to be extrapolated, create a feature-rich row for each (product, date) pair and predict the target (units sold) using all available information.

Data Sources

The forecasting system pulls from six data sources:

Source	Data	Update Frequency
Sales history	Daily units sold per SKU	Daily
Promotions calendar	Discount %, promo type, duration	Weekly (planned ahead)
Product catalog	Category, price, brand, launch date	On change
Weather API	Temperature, precipitation by region	Daily (forecast for +14 days)
Marketing spend	Channel-level daily ad spend	Daily
Competitor pricing	Scraped prices for top 50 competing products	Daily

Feature Engineering for Demand Forecasting

Feature engineering is where the ML approach earns its advantage over classical methods. Build features in five categories:

1. Calendar and Seasonality Features

import numpy as np


def build_calendar_features(dates: pd.Series) -> pd.DataFrame:
    """Extract calendar features from a date series."""

    df = pd.DataFrame({"date": dates})

    # Basic calendar
    df["day_of_week"] = df["date"].dt.dayofweek  # 0=Monday, 6=Sunday
    df["day_of_month"] = df["date"].dt.day
    df["week_of_year"] = df["date"].dt.isocalendar().week.astype(int)
    df["month"] = df["date"].dt.month
    df["quarter"] = df["date"].dt.quarter
    df["year"] = df["date"].dt.year

    # Is it a weekend?
    df["is_weekend"] = (df["day_of_week"] >= 5).astype(int)

    # Pay period effects (people buy more around the 1st and 15th)
    df["is_month_start"] = (df["day_of_month"] <= 3).astype(int)
    df["is_month_mid"] = ((df["day_of_month"] >= 14) & (df["day_of_month"] <= 17)).astype(int)

    # US holidays (simplified -- use the `holidays` library for production)
    us_holidays = {
        # 2025-2026 major holidays
        "2025-11-27", "2025-11-28",  # Thanksgiving + Black Friday
        "2025-12-24", "2025-12-25",  # Christmas
        "2025-12-31", "2026-01-01",  # New Year
        "2026-01-20",                # MLK Day
        "2026-02-14",                # Valentine's Day
        "2026-05-25",                # Memorial Day
        "2026-07-04",                # Independence Day
    }

    df["is_holiday"] = df["date"].dt.strftime("%Y-%m-%d").isin(us_holidays).astype(int)

    # Days to/from nearest holiday (captures pre-holiday shopping buildup)
    holiday_dates = pd.to_datetime(list(us_holidays))
    df["days_to_next_holiday"] = df["date"].apply(
        lambda d: min((h - d).days for h in holiday_dates if h >= d) if any(h >= d for h in holiday_dates) else 999
    )
    df["days_from_last_holiday"] = df["date"].apply(
        lambda d: min((d - h).days for h in holiday_dates if h <= d) if any(h <= d for h in holiday_dates) else 999
    )

    # Cyclical encoding for periodic features (helps tree models too)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    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)

    return df.drop(columns=["date"])

2. Lag and Rolling Window Features

This is where the historical demand signal gets encoded. The model needs to know "how much did this product sell recently?" without leaking future information.

def build_lag_features(sales_df: pd.DataFrame, target_col: str = "units_sold") -> pd.DataFrame:
    """Build lag and rolling features for each SKU.

    IMPORTANT: sales_df must be sorted by (sku_id, date) before calling this.
    """

    df = sales_df.copy()

    # Simple lags (1 day, 7 days, 14 days, 28 days ago)
    for lag in [1, 7, 14, 28]:
        df[f"lag_{lag}d"] = df.groupby("sku_id")[target_col].shift(lag)

    # Rolling means (7-day, 14-day, 28-day windows)
    for window in [7, 14, 28]:
        df[f"rolling_mean_{window}d"] = (
            df.groupby("sku_id")[target_col]
            .transform(lambda x: x.shift(1).rolling(window, min_periods=1).mean())
        )
        df[f"rolling_std_{window}d"] = (
            df.groupby("sku_id")[target_col]
            .transform(lambda x: x.shift(1).rolling(window, min_periods=1).std())
        )

    # Rolling max/min (captures spikes and troughs)
    df["rolling_max_14d"] = (
        df.groupby("sku_id")[target_col]
        .transform(lambda x: x.shift(1).rolling(14, min_periods=1).max())
    )

    # Same day last week, 2 weeks ago, 4 weeks ago
    df["same_dow_1w"] = df.groupby("sku_id")[target_col].shift(7)
    df["same_dow_2w"] = df.groupby("sku_id")[target_col].shift(14)
    df["same_dow_4w"] = df.groupby("sku_id")[target_col].shift(28)

    # Trend: compare recent average to older average
    df["demand_trend_7v28"] = df["rolling_mean_7d"] / df["rolling_mean_28d"].clip(lower=0.1)

    # Coefficient of variation (demand volatility)
    df["demand_cv_28d"] = df["rolling_std_28d"] / df["rolling_mean_28d"].clip(lower=0.1)

    return df

3. Promotion Features

def build_promotion_features(sales_df: pd.DataFrame, promos_df: pd.DataFrame) -> pd.DataFrame:
    """Merge promotion information into the sales dataset."""

    df = sales_df.merge(
        promos_df[["sku_id", "date", "discount_pct", "promo_type", "is_featured"]],
        on=["sku_id", "date"],
        how="left"
    )

    # Fill non-promo days
    df["discount_pct"] = df["discount_pct"].fillna(0)
    df["is_on_promo"] = (df["discount_pct"] > 0).astype(int)
    df["is_featured"] = df["is_featured"].fillna(0).astype(int)

    # Promo type encoding
    promo_types = {"none": 0, "flash_sale": 1, "weekly_deal": 2, "clearance": 3, "bundle": 4}
    df["promo_type_encoded"] = df["promo_type"].fillna("none").map(promo_types).fillna(0)

    # Days since last promotion (captures promo fatigue)
    df["days_since_last_promo"] = (
        df.groupby("sku_id")["is_on_promo"]
        .transform(lambda x: x.cumsum().diff().eq(0).cumsum())
    )

    # Upcoming promotion indicator (if promo calendar is known ahead)
    df["promo_in_next_3d"] = (
        df.groupby("sku_id")["is_on_promo"]
        .transform(lambda x: x.shift(-1).rolling(3, min_periods=1).max())
    ).fillna(0)

    return df

4. External Features (Weather, Competitor Pricing)

def build_external_features(sales_df: pd.DataFrame, weather_df: pd.DataFrame,
                             competitor_df: pd.DataFrame) -> pd.DataFrame:
    """Add weather and competitor pricing features."""

    df = sales_df.copy()

    # Weather features (by region/date)
    weather_agg = weather_df.groupby("date").agg(
        avg_temp=("temperature_f", "mean"),
        total_precip=("precipitation_in", "sum"),
        is_extreme_weather=("weather_alert", "max"),
    ).reset_index()

    df = df.merge(weather_agg, on="date", how="left")

    # Temperature deviation from seasonal normal
    df["temp_deviation"] = df["avg_temp"] - df.groupby(
        df["date"].dt.month
    )["avg_temp"].transform("mean")

    # Competitor pricing
    # competitor_df has: sku_id, date, competitor_price, our_price
    if competitor_df is not None:
        price_comparison = competitor_df.groupby(["sku_id", "date"]).agg(
            min_competitor_price=("competitor_price", "min"),
            avg_competitor_price=("competitor_price", "mean"),
        ).reset_index()

        df = df.merge(price_comparison, on=["sku_id", "date"], how="left")

        # Price ratio (our price vs competitor)
        df["price_ratio_vs_min"] = df["price"] / df["min_competitor_price"].clip(lower=0.01)
        df["price_ratio_vs_avg"] = df["price"] / df["avg_competitor_price"].clip(lower=0.01)

        # Are we the cheapest?
        df["is_cheapest"] = (df["price"] <= df["min_competitor_price"]).astype(int)

    return df

5. Product-Level Features

def build_product_features(sales_df: pd.DataFrame, catalog_df: pd.DataFrame) -> pd.DataFrame:
    """Add product attributes and lifecycle features."""

    df = sales_df.merge(
        catalog_df[["sku_id", "category", "subcategory", "brand", "price",
                     "launch_date", "weight_kg"]],
        on="sku_id",
        how="left"
    )

    # Product age (days since launch)
    df["product_age_days"] = (df["date"] - df["launch_date"]).dt.days

    # Is it a new product? (first 90 days)
    df["is_new_product"] = (df["product_age_days"] <= 90).astype(int)

    # Category encoding (target encoding is better but simpler for now)
    # In production, use category_encoders.TargetEncoder
    category_map = dict(enumerate(catalog_df["category"].unique()))
    inv_map = {v: k for k, v in category_map.items()}
    df["category_encoded"] = df["category"].map(inv_map).fillna(0)

    # Price tier
    df["price_tier"] = pd.qcut(df["price"], q=5, labels=False, duplicates="drop")

    return df

Time-Series Cross-Validation

This is the part that many ML practitioners get wrong when applying tabular methods to time series. Random k-fold cross-validation does not work here. If the training set includes December 15 and the validation set includes December 10, the model has seen the future. Use expanding window or sliding window cross-validation.

from sklearn.model_selection import TimeSeriesSplit


def time_series_cv(df, n_splits=5, test_days=14, gap_days=1):
    """Custom time series cross-validation with configurable windows.

    Training window expands, test window is fixed-size, with a gap
    between train and test to simulate production latency.
    """

    all_dates = sorted(df["date"].unique())
    total_dates = len(all_dates)
    min_train_dates = total_dates // 3  # at least 1/3 of data for training

    splits = []
    for i in range(n_splits):
        test_end_idx = total_dates - (i * test_days)
        test_start_idx = test_end_idx - test_days
        train_end_idx = test_start_idx - gap_days

        if train_end_idx < min_train_dates:
            break

        train_dates = set(all_dates[:train_end_idx])
        test_dates = set(all_dates[test_start_idx:test_end_idx])

        train_idx = df[df["date"].isin(train_dates)].index
        test_idx = df[df["date"].isin(test_dates)].index

        splits.append((train_idx, test_idx))

    return list(reversed(splits))  # chronological order


# Usage with the feature-enriched dataset
splits = time_series_cv(full_df, n_splits=5, test_days=14, gap_days=1)

for fold_idx, (train_idx, test_idx) in enumerate(splits):
    train_dates = full_df.loc[train_idx, "date"]
    test_dates = full_df.loc[test_idx, "date"]
    print(
        f"Fold {fold_idx}: "
        f"Train {train_dates.min().date()} to {train_dates.max().date()} "
        f"({len(train_idx):,} rows) -> "
        f"Test {test_dates.min().date()} to {test_dates.max().date()} "
        f"({len(test_idx):,} rows)"
    )

Model Training: LightGBM for Demand

LightGBM suits the final model for its speed (training on millions of rows across thousands of SKUs) and its native support for categorical features.

import lightgbm as lgb
from sklearn.metrics import mean_absolute_error, mean_squared_error

# Prepare feature columns (exclude identifiers and target)
exclude_cols = ["sku_id", "date", "units_sold", "launch_date", "category",
                "subcategory", "brand"]
feature_cols = [c for c in full_df.columns if c not in exclude_cols]

# Train with time series CV
cv_results = []

for fold_idx, (train_idx, test_idx) in enumerate(splits):
    X_train = full_df.loc[train_idx, feature_cols].values
    y_train = full_df.loc[train_idx, "units_sold"].values
    X_test = full_df.loc[test_idx, feature_cols].values
    y_test = full_df.loc[test_idx, "units_sold"].values

    model = lgb.LGBMRegressor(
        n_estimators=1000,
        max_depth=8,
        learning_rate=0.05,
        num_leaves=63,
        subsample=0.8,
        colsample_bytree=0.8,
        min_child_samples=20,
        reg_alpha=0.1,
        reg_lambda=1.0,
        random_state=42,
        n_jobs=-1,
        verbose=-1,
    )

    model.fit(
        X_train, y_train,
        eval_set=[(X_test, y_test)],
        callbacks=[lgb.early_stopping(50, verbose=False)],
    )

    y_pred = model.predict(X_test)
    y_pred = np.maximum(y_pred, 0)  # demand cannot be negative

    mae = mean_absolute_error(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    mape = np.mean(np.abs(y_test - y_pred) / np.maximum(y_test, 1)) * 100

    cv_results.append({"fold": fold_idx, "mae": mae, "rmse": rmse, "mape": mape})
    print(f"Fold {fold_idx}: MAE={mae:.1f}, RMSE={rmse:.1f}, MAPE={mape:.1f}%")

# Average results
results_df = pd.DataFrame(cv_results)
print(f"\nAverage: MAE={results_df['mae'].mean():.1f}, "
      f"RMSE={results_df['rmse'].mean():.1f}, "
      f"MAPE={results_df['mape'].mean():.1f}%")

# Results:
# Average: MAE=11.8, RMSE=19.3, MAPE=18.4%
# (Compare to Prophet's 31.2% MAPE -- 41% improvement)

Feature Importance Analysis

# Train final model on all data up to the cutoff
final_model = lgb.LGBMRegressor(
    n_estimators=800, max_depth=8, learning_rate=0.05,
    num_leaves=63, subsample=0.8, colsample_bytree=0.8,
    random_state=42, n_jobs=-1, verbose=-1
)

cutoff_date = pd.Timestamp("2026-01-01")
train_mask = full_df["date"] < cutoff_date
final_model.fit(
    full_df.loc[train_mask, feature_cols].values,
    full_df.loc[train_mask, "units_sold"].values
)

# Feature importance
importance = pd.DataFrame({
    "feature": feature_cols,
    "importance": final_model.feature_importances_
}).sort_values("importance", ascending=False)

print("Top 20 features:")
print(importance.head(20).to_string(index=False))

# Top features in the model:
# 1. rolling_mean_7d         (3842)
# 2. lag_7d                  (3201)
# 3. rolling_mean_28d        (2987)
# 4. discount_pct            (2756)
# 5. same_dow_1w             (2534)
# 6. price_ratio_vs_avg      (2198)
# 7. days_to_next_holiday    (1876)
# 8. lag_1d                  (1743)
# 9. demand_trend_7v28       (1654)
# 10. avg_temp               (1432)

The results confirm the hypothesis: recent demand history (lags and rolling means) is the strongest signal, but promotions, competitor pricing, and calendar effects all contribute meaningfully. Weather mattered for specific categories (outdoor equipment, beverages) but less so for electronics.

Multi-Step Forecasting: Predicting 14 Days Ahead

In production, the operations team needs forecasts for the next 14 days to place orders with suppliers. There are three strategies for multi-step forecasting:

Strategy 1: Recursive (Autoregressive)

Predict day 1, use that prediction as a lag feature for day 2, use day 2's prediction for day 3, and so on. This compounds errors but respects temporal dependencies.

Strategy 2: Direct (Separate Models)

Train a separate model for each forecast horizon: one model predicts 1 day ahead, another predicts 2 days ahead, up to 14 days. This avoids error compounding but requires 14 models.

Strategy 3: Hybrid

A hybrid approach uses recursive forecasting but caps the lag features at 7 days old. Beyond that, it relies on rolling statistics that are less sensitive to individual prediction errors.

def forecast_multi_step(model, current_features_df, sku_ids, n_days=14, feature_cols=None):
    """Generate multi-step forecasts using recursive prediction with lag clipping."""

    forecasts = []
    # Start with known features as of today
    working_df = current_features_df.copy()

    for day_offset in range(1, n_days + 1):
        forecast_date = pd.Timestamp.now().normalize() + timedelta(days=day_offset)

        # Update calendar features for the forecast date
        working_df = update_calendar_features(working_df, forecast_date)

        # Update promotion features (from the promo calendar)
        working_df = update_promo_features(working_df, forecast_date)

        # Predict
        X = working_df[feature_cols].values
        predictions = model.predict(X)
        predictions = np.maximum(predictions, 0)

        # Store forecasts
        for sku_id, pred in zip(sku_ids, predictions):
            forecasts.append({
                "sku_id": sku_id,
                "forecast_date": forecast_date,
                "predicted_units": round(float(pred), 1),
                "forecast_horizon": day_offset,
            })

        # Update lag features for the next step
        # Only use predicted values for short lags (1-7 days)
        # Keep rolling statistics from actual data for stability
        working_df = update_lag_features(working_df, predictions, day_offset)

    return pd.DataFrame(forecasts)


def update_lag_features(df, predictions, day_offset):
    """Update lag features with predicted values, but cap at 7 days."""

    df = df.copy()
    if day_offset <= 7:
        # Shift lags and insert prediction
        for lag in [1, 7, 14, 28]:
            if lag == day_offset:
                df[f"lag_{lag}d"] = predictions
            elif lag > day_offset:
                # This lag still uses actual data
                pass

    # Rolling features: only update with a blended approach
    # Use exponential smoothing to reduce compounding error
    alpha = max(0.3, 1.0 - day_offset * 0.05)  # decay with horizon
    df["rolling_mean_7d"] = alpha * predictions + (1 - alpha) * df["rolling_mean_7d"]

    return df

Handling Multiple SKUs at Scale

A real e-commerce business has thousands of SKUs. Training a single global model (one model for all SKUs) versus per-SKU models is a critical design decision.

Approach	Pros	Cons
Global model (all SKUs)	Learns cross-product patterns, works for new products, simpler to maintain	May miss SKU-specific patterns
Per-SKU models	Captures unique demand patterns	Requires enough history per SKU, thousands of models to maintain
Cluster + model per cluster	Balanced: similar products share a model	Cluster definitions need maintenance

A global model with SKU-level features works well here. The product category, price tier, and demand volatility features let the model learn different behavior for different product types without maintaining thousands of separate models.

# Group SKUs by demand pattern for reporting
def classify_sku_demand(sales_df):
    """Classify SKUs into demand pattern categories."""

    sku_stats = sales_df.groupby("sku_id").agg(
        mean_daily_sales=("units_sold", "mean"),
        std_daily_sales=("units_sold", "std"),
        zero_days_pct=("units_sold", lambda x: (x == 0).mean()),
        total_days=("units_sold", "count"),
    ).reset_index()

    sku_stats["cv"] = sku_stats["std_daily_sales"] / sku_stats["mean_daily_sales"].clip(lower=0.01)

    conditions = [
        (sku_stats["zero_days_pct"] > 0.7),                    # Intermittent
        (sku_stats["cv"] > 1.5),                                # Lumpy
        (sku_stats["mean_daily_sales"] < 1),                    # Slow-moving
        (sku_stats["cv"] < 0.5),                                # Smooth
    ]
    choices = ["intermittent", "lumpy", "slow_moving", "smooth"]
    sku_stats["demand_pattern"] = np.select(conditions, choices, default="regular")

    print(sku_stats["demand_pattern"].value_counts())
    return sku_stats

# Typical distribution:
# regular        1,247
# smooth           834
# slow_moving      523
# lumpy            312
# intermittent     284

Intermittent and lumpy demand patterns are the hardest to forecast. For these SKUs, predicting weekly totals (instead of daily) and then distributing across days gave better results than direct daily prediction.

Deploying Forecasts to the Data Warehouse

The forecast is only useful if the operations team can access it. Forecasts are written to a BigQuery table consumed by the inventory planning tool, Looker dashboards, and automated reorder alerts.

"""
dags/demand_forecast_dag.py
Daily demand forecasting pipeline.
"""
from airflow import DAG
from airflow.operators.python import PythonOperator
from airflow.providers.google.cloud.operators.bigquery import BigQueryInsertJobOperator
from datetime import datetime, timedelta


default_args = {
    "owner": "data-team",
    "depends_on_past": False,
    "retries": 2,
    "retry_delay": timedelta(minutes=15),
}

dag = DAG(
    "demand_forecast_daily",
    default_args=default_args,
    schedule_interval="0 4 * * *",  # 4 AM daily
    start_date=datetime(2026, 1, 1),
    catchup=False,
)


def generate_forecasts(**context):
    """Run the forecasting model and produce 14-day predictions."""
    import joblib
    from utils.forecast import (
        load_latest_data, build_all_features, forecast_multi_step
    )

    # Load model
    artifact = joblib.load("/models/demand/latest/model.joblib")
    model = artifact["model"]
    feature_cols = artifact["feature_columns"]

    # Load and prepare data
    sales, promos, weather, competitors, catalog = load_latest_data()
    features_df, sku_ids = build_all_features(
        sales, promos, weather, competitors, catalog,
        prediction_date=pd.Timestamp(context["ds"])
    )

    # Generate forecasts
    forecasts = forecast_multi_step(
        model, features_df, sku_ids,
        n_days=14, feature_cols=feature_cols
    )

    # Add metadata
    forecasts["model_version"] = artifact["version"]
    forecasts["generated_at"] = datetime.now().isoformat()

    # Write to staging table
    output_path = f"/tmp/forecasts_{context['ds']}.parquet"
    forecasts.to_parquet(output_path)
    context["ti"].xcom_push(key="forecast_path", value=output_path)
    context["ti"].xcom_push(key="n_forecasts", value=len(forecasts))


def validate_forecasts(**context):
    """Basic sanity checks before publishing."""
    forecast_path = context["ti"].xcom_pull(key="forecast_path")
    forecasts = pd.read_parquet(forecast_path)

    checks = []

    # Check 1: No negative predictions
    n_negative = (forecasts["predicted_units"] < 0).sum()
    checks.append(("no_negative", n_negative == 0))

    # Check 2: Reasonable range (no single SKU predicting 10x historical max)
    # This would use historical stats from the warehouse
    max_pred = forecasts["predicted_units"].max()
    checks.append(("reasonable_max", max_pred < 10000))

    # Check 3: Coverage (we should have forecasts for all active SKUs)
    n_skus = forecasts["sku_id"].nunique()
    checks.append(("sku_coverage", n_skus >= 3000))  # we have ~3,200 active SKUs

    failed = [c[0] for c in checks if not c[1]]
    if failed:
        raise ValueError(f"Forecast validation failed: {failed}")

    print(f"All checks passed. {n_skus} SKUs, {len(forecasts)} forecast rows.")


def load_to_bigquery(**context):
    """Load validated forecasts to BigQuery."""
    from google.cloud import bigquery

    forecast_path = context["ti"].xcom_pull(key="forecast_path")

    client = bigquery.Client()
    job_config = bigquery.LoadJobConfig(
        source_format=bigquery.SourceFormat.PARQUET,
        write_disposition=bigquery.WriteDisposition.WRITE_APPEND,
    )

    with open(forecast_path, "rb") as f:
        job = client.load_table_from_file(
            f,
            "analytics.demand_forecasts",
            job_config=job_config
        )
    job.result()
    print(f"Loaded {job.output_rows} rows to BigQuery")


# Task definitions
forecast = PythonOperator(
    task_id="generate_forecasts", python_callable=generate_forecasts, dag=dag
)
validate = PythonOperator(
    task_id="validate_forecasts", python_callable=validate_forecasts, dag=dag
)
load = PythonOperator(
    task_id="load_to_bigquery", python_callable=load_to_bigquery, dag=dag
)

forecast >> validate >> load

The Warehouse Schema

-- BigQuery table for consuming forecasts
CREATE TABLE IF NOT EXISTS analytics.demand_forecasts (
    sku_id              STRING NOT NULL,
    forecast_date       DATE NOT NULL,
    predicted_units     FLOAT64,
    forecast_horizon    INT64,
    model_version       STRING,
    generated_at        TIMESTAMP,

    -- Partitioned by forecast generation date for easy cleanup
)
PARTITION BY DATE(generated_at)
CLUSTER BY sku_id, forecast_date;


-- View for the operations team: latest forecast per SKU per date
CREATE OR REPLACE VIEW analytics.latest_demand_forecast AS
SELECT
    f.sku_id,
    f.forecast_date,
    f.predicted_units,
    f.forecast_horizon,
    f.model_version,
    f.generated_at,
    p.category,
    p.brand,
    p.current_stock,
    p.reorder_point,
    CASE
        WHEN p.current_stock < f.predicted_units * 3 THEN 'REORDER'
        WHEN p.current_stock < f.predicted_units * 7 THEN 'WATCH'
        ELSE 'OK'
    END AS stock_status
FROM analytics.demand_forecasts f
JOIN analytics.product_catalog p ON f.sku_id = p.sku_id
WHERE f.generated_at = (
    SELECT MAX(generated_at) FROM analytics.demand_forecasts
);


-- Alert query: SKUs that need reordering based on forecasted demand
SELECT
    sku_id,
    category,
    brand,
    current_stock,
    SUM(predicted_units) AS total_predicted_14d,
    current_stock / NULLIF(SUM(predicted_units), 0) AS days_of_stock_remaining
FROM analytics.latest_demand_forecast
GROUP BY sku_id, category, brand, current_stock
HAVING days_of_stock_remaining < 5
ORDER BY days_of_stock_remaining ASC;

Results: 28% Reduction in Stockouts

After three months in production, the numbers:

Metric	Before (Prophet + Spreadsheets)	After (LightGBM Pipeline)	Improvement
Overall MAPE	31.2%	18.4%	-41%
Stockout rate	14.0%	10.1%	-28%
Excess inventory	$2.3M	$1.41M	-$890K
Forecast generation time	4 hours (manual)	12 minutes (automated)	-95%
SKU coverage	Top 200 (manual)	All 3,200	+1,500%

The biggest win was not accuracy improvement per se -- it was coverage. The old process only forecasted the top 200 SKUs because doing more was too time-consuming manually. The ML pipeline forecasts every active SKU every day, which surfaced stockout risks for long-tail products that were previously ignored.

Where the Model Still Struggles

The categories where the model performs poorly:

• New product launches: With no sales history, lag features are all zeros. The fallback is category-level averages for the first 30 days, which is mediocre at best.
• Viral/social media spikes: When a product goes viral on TikTok, demand can 10x overnight. The model has no social media features and cannot predict these events.
• Supply-constrained products: When a product is out of stock, observed sales are zero, but actual demand is unknown. The model sees zero sales and under-predicts future demand. This is a censored data problem with no clean solution.
• Flash sales with no historical precedent: If the marketing team runs a 50% off sale for the first time ever, the model has no data to calibrate the promotional uplift.

Practical Tips for a Forecasting System

After building this system and iterating on it for six months, the lessons that save the most time:

• Start with a simple baseline and measure it honestly. A first model of "average of same weekday over the last 4 weeks" had a MAPE of 25% -- not terrible. Knowing this baseline grounds the evaluation of fancier approaches.
• Feature engineering is the differentiator. The jump from Prophet (31% MAPE) to a first LightGBM without external features (23% MAPE) was big. The jump from that to the full feature set (18.4% MAPE) was almost as big. Model architecture matters much less than feature quality.
• Use time-series cross-validation religiously. A first attempt with random train-test splitting reported 12% MAPE. Switching to proper temporal splitting moved that to 22%. The random split was leaking future information through lag features.
• Monitor by SKU segment, not just overall. An overall MAPE of 18.4% hid the fact that smooth-demand SKUs had 8% MAPE while intermittent SKUs had 45% MAPE. The operations team needed to know which forecasts to trust.
• Automate the full loop. A forecast that requires manual intervention breaks over weekends, holidays, and when someone is on vacation. The Airflow pipeline runs seven days a week without human involvement.

From Spreadsheets to Production

Demand forecasting is one of the highest-ROI applications of machine learning. The business case is concrete (stockouts cost money, excess inventory costs money), the data is usually available (every company has sales history), and the results are immediately measurable.

For a process currently running on spreadsheets or a basic Prophet model, a reasonable migration path:

• Week 1-2: Build the feature engineering pipeline. Get lag features, calendar features, and promotion features into a clean dataset. This is the most important work.
• Week 3: Train a LightGBM model with proper time-series cross-validation. Compare against the current baseline. If the improvement is less than 10% MAPE, the features need more work, not a fancier model.
• Week 4: Deploy the daily pipeline. Write forecasts to the warehouse. Get the operations team using the output.
• Month 2+: Add external features (weather, competitor pricing, marketing spend). Set up monitoring. Build the retraining automation.

The entire system described here -- from data extraction through feature engineering, model training, multi-step forecasting, and warehouse deployment -- can be built by a team of two data engineers in about six weeks, for a system that saved nearly $900K in the first quarter. A mediocre ML forecast that runs automatically every day is more valuable than a perfect model that only exists in a notebook. Ship the pipeline first, then iterate on accuracy.