Table of Contents

    Sales Forecasting

    Machine Learning Project

    Sales Forecasting

    Build an end-to-end time-series forecasting project that predicts future sales using historical observations, trend and seasonality analysis, naïve baselines, lag features, rolling statistics, walk-forward validation, Gradient Boosting, forecast intervals, reproducible artifacts, and production monitoring.

    Project Introduction

    Sales forecasting estimates future sales from historical observations and other information available before the forecast is produced. Forecasts can support inventory planning, purchasing, staffing, production, budgeting, logistics, and commercial decision-making.

    Unlike ordinary regression, time-series forecasting must preserve chronological order. A model trained using future information can appear highly accurate during testing while failing in actual deployment.

    A reliable forecasting project therefore requires a precisely defined target, forecast frequency, horizon, data-availability rules, chronological validation, strong baselines, suitable error metrics, uncertainty estimates, and continuous monitoring.

    Project goal: Build a reproducible forecasting pipeline that predicts future sales without using information unavailable at the forecast creation time.

    Learning Objectives

    By the End of This Project

    • Define the sales target, frequency, horizon, and forecast origin
    • Distinguish sales, demand, orders, shipments, and revenue
    • Load and validate an authorized sales dataset
    • Aggregate transactional data into a regular time series
    • Detect duplicate records and missing periods
    • Distinguish zero sales from missing observations
    • Analyze trend, seasonality, volatility, and anomalies
    • Create naïve, seasonal naïve, and moving-average baselines
    • Engineer lag and rolling-window features safely
    • Avoid temporal and target leakage
    • Create chronological training, validation, and test sets
    • Use expanding-window and walk-forward validation
    • Train a Gradient Boosting forecasting model
    • Evaluate MAE, RMSE, WAPE, bias, and MASE
    • Compare errors across forecast horizons
    • Generate recursive multi-step forecasts
    • Create empirical forecast intervals
    • Analyze high-error periods
    • Save model artifacts and forecasting metadata
    • Monitor data quality, accuracy, bias, and drift

    Prerequisites

    Recommended Knowledge

    • Python functions, classes, loops, and dictionaries
    • pandas DataFrames and date-time indexes
    • NumPy arrays
    • Matplotlib and Seaborn
    • Regression fundamentals
    • MAE, RMSE, and residual analysis
    • Time-series trend and seasonality
    • Lag variables and rolling statistics
    • Train-validation-test separation
    • Basic business understanding of sales and demand

    Required Packages

    • pandas
    • NumPy
    • Matplotlib
    • Seaborn
    • scikit-learn
    • joblib
    python -m pip install pandas numpy matplotlib seaborn scikit-learn joblib

    Define the Forecasting Problem

    Before choosing a forecasting model, define exactly what must be predicted, when the forecast is generated, and how far into the future it must extend.

    Forecast Element Question to Answer
    Target Are you forecasting units, orders, gross sales, net sales, or revenue?
    Granularity Is the forecast required by product, store, channel, territory, or total business?
    Frequency Are observations daily, weekly, monthly, or quarterly?
    Forecast Origin At what timestamp is the forecast created?
    Horizon How many future periods must be forecast?
    Refresh Frequency How often will the forecast be recalculated?
    Decision Which action will use the forecast?
    Success Measure Which forecasting errors matter operationally?
    Example project definition: Predict total weekly units sold for the next four weeks using information available at the end of the current week.

    Sales, Demand, Orders, and Revenue

    Measure Meaning Important Limitation
    Units Sold Quantity successfully sold May be limited by stock availability
    Customer Demand Quantity customers wanted to purchase May not be fully observed during stockouts
    Orders Quantity ordered by customers May include cancellations or unfulfilled orders
    Shipments Quantity dispatched to customers Can be delayed relative to demand
    Gross Sales Sales value before deductions May differ from realized revenue
    Net Sales Sales value after defined returns or discounts Adjustments may become known later
    Revenue Recognized financial value under applicable rules Recognition timing may differ from transaction timing
    Incorrect Assumption Observed sales always equal unconstrained customer demand.
    Correct Understanding Stockouts, capacity limits, fulfillment failures, and commercial policies can cause observed sales to be lower than true demand.

    Forecast Origin and Horizon

    Forecast at Time \(t\)
    \[ \hat{y}_{t+h\mid t} \]

    The notation represents a forecast created at time \(t\) for a future value \(h\) periods ahead.

    Symbol Meaning
    \(t\) Forecast origin
    \(h\) Forecast horizon
    \(y_{t+h}\) Future observed sales
    \(\hat{y}_{t+h\mid t}\) Forecast for period \(t+h\) created at time \(t\)

    Expected Dataset

    Save an authorized transactional dataset as sales_history.csv. The project uses the following logical schema.

    Column Type Description
    date Date Transaction or reporting date
    sales Numerical Sales target in consistent units
    product_id Categorical Optional product or product-family identifier
    store_id Categorical Optional store or location identifier
    promotion Binary or categorical Promotion information known before the forecasted period
    price Numerical Price known or planned for the relevant period
    stock_available Numerical Optional availability information recorded at the correct time
    Data notice: This project does not create synthetic sales history. Use an authorized and documented dataset. Update the target name, dimensions, units, and aggregation logic when your schema differs.

    Data Governance

    Review Before Forecasting

    • Confirm authority to use transaction and customer data
    • Document the source systems and extraction timestamp
    • Define the currency, quantity, and area units
    • Record whether values are gross or net of returns
    • Document cancellations and correction handling
    • Confirm transaction and reporting time zones
    • Identify delayed postings and backdated transactions
    • Document product and store hierarchy changes
    • Protect confidential commercial information
    • Preserve the dataset version used for training

    Load the Dataset

    from pathlib import Path
    
    import numpy as np
    import pandas as pd
    
    
    data_path = Path(
        "sales_history.csv"
    )
    
    
    if not data_path.exists():
    
        raise FileNotFoundError(
            "sales_history.csv was not found. "
            "Add an authorized sales dataset before running the project."
        )
    
    
    sales_data = pd.read_csv(
        data_path
    )
    
    
    required_columns = {
        "date",
        "sales"
    }
    
    
    missing_columns = (
        required_columns
        - set(sales_data.columns)
    )
    
    
    if missing_columns:
    
        raise ValueError(
            "Required columns are missing: "
            f"{sorted(missing_columns)}"
        )
    
    
    sales_data["date"] = pd.to_datetime(
        sales_data["date"],
        errors="coerce"
    )
    
    
    sales_data["sales"] = pd.to_numeric(
        sales_data["sales"],
        errors="coerce"
    )
    
    
    invalid_rows = sales_data[
        sales_data["date"].isna()
        |
        sales_data["sales"].isna()
    ]
    
    
    print(
        "Invalid rows:",
        len(invalid_rows)
    )
    
    
    sales_data = sales_data.dropna(
        subset=[
            "date",
            "sales"
        ]
    ).copy()
    
    
    sales_data = sales_data.sort_values(
        by="date"
    )
    
    
    print(
        sales_data.head()
    )
    
    print(
        sales_data.shape
    )

    Initial Data Audit

    print(
        sales_data.info()
    )
    
    
    print(
        sales_data.describe(
            include="all"
        ).transpose()
    )
    
    
    duplicate_rows = sales_data.duplicated().sum()
    
    
    negative_sales_count = (
        sales_data["sales"]
        <
        0
    ).sum()
    
    
    print(
        "Complete duplicate rows:",
        duplicate_rows
    )
    
    
    print(
        "Negative sales rows:",
        negative_sales_count
    )
    
    
    print(
        "Minimum date:",
        sales_data["date"].min()
    )
    
    
    print(
        "Maximum date:",
        sales_data["date"].max()
    )
    Negative-value caution: Negative sales may represent returns, reversals, corrections, or data errors. Investigate their meaning before removing or changing them.

    Aggregate to a Regular Frequency

    The following example aggregates transaction-level values into weekly totals. Change the frequency only after confirming the planning requirement.

    weekly_sales = (
        sales_data
        .set_index(
            "date"
        )["sales"]
        .resample(
            "W-SUN"
        )
        .sum(
            min_count=1
        )
        .rename(
            "sales"
        )
        .to_frame()
    )
    
    
    print(
        weekly_sales.head()
    )
    
    print(
        weekly_sales.tail()
    )
    Frequency Possible Use Important Consideration
    Daily Short-term staffing or replenishment May contain strong weekday effects and noise
    Weekly Operational planning Week definitions must remain consistent
    Monthly Budgeting and capacity planning Months have different lengths
    Quarterly Strategic planning Provides fewer observations for modeling

    Detect Missing Periods

    complete_weekly_index = pd.date_range(
        start=weekly_sales.index.min(),
        end=weekly_sales.index.max(),
        freq="W-SUN"
    )
    
    
    missing_periods = complete_weekly_index.difference(
        weekly_sales.index
    )
    
    
    print(
        "Missing weekly periods:",
        len(missing_periods)
    )
    
    
    print(
        missing_periods[:20]
    )
    Incorrect Assumption A missing date automatically means zero sales.
    Correct Review Determine whether the period had no sales, the business was closed, data was not loaded, or reporting failed.

    Zero Sales vs Missing Data

    Condition Meaning Possible Handling
    Confirmed Zero Sales The business operated but no sales occurred Retain zero
    Planned Closure The business was not operating Represent closure explicitly
    Missing Extraction Sales may exist but data was not loaded Recover data or mark it missing
    Stockout Observed sales may be constrained below demand Add availability context where appropriate
    New Product No historical sales exist before introduction Use a cold-start strategy

    Plot the Sales Series

    import matplotlib.pyplot as plt
    import seaborn as sns
    
    
    plt.figure(
        figsize=(12, 5)
    )
    
    
    plt.plot(
        weekly_sales.index,
        weekly_sales["sales"]
    )
    
    
    plt.title(
        "Weekly Sales History"
    )
    
    plt.xlabel(
        "Week"
    )
    
    plt.ylabel(
        "Sales"
    )
    
    plt.tight_layout()
    
    plt.show()

    Inspect the Series For

    • Long-term growth or decline
    • Weekly, monthly, quarterly, or annual seasonality
    • Promotion-driven spikes
    • Stockout-related drops
    • Structural breaks
    • Long periods of zero sales
    • Unexplained outliers
    • Changes in volatility

    Time-Series Components

    Additive Representation
    \[ y_t = T_t + S_t + R_t \]
    Component Meaning
    Level Typical magnitude of the series
    Trend Persistent upward or downward movement
    Seasonality Pattern repeating at a known interval
    Cycle Longer variation without a fixed seasonal period
    Residual Variation not explained by the modeled structure

    Rolling Mean for Exploration

    exploration_data = weekly_sales.copy()
    
    
    exploration_data[
        "rolling_mean_4"
    ] = (
        exploration_data["sales"]
        .rolling(
            window=4
        )
        .mean()
    )
    
    
    exploration_data[
        "rolling_mean_13"
    ] = (
        exploration_data["sales"]
        .rolling(
            window=13
        )
        .mean()
    )
    
    
    plt.figure(
        figsize=(12, 6)
    )
    
    
    plt.plot(
        exploration_data.index,
        exploration_data["sales"],
        label="Weekly Sales",
        alpha=0.55
    )
    
    
    plt.plot(
        exploration_data.index,
        exploration_data["rolling_mean_4"],
        label="4-Week Mean"
    )
    
    
    plt.plot(
        exploration_data.index,
        exploration_data["rolling_mean_13"],
        label="13-Week Mean"
    )
    
    
    plt.title(
        "Sales and Rolling Means"
    )
    
    plt.xlabel(
        "Week"
    )
    
    plt.ylabel(
        "Sales"
    )
    
    plt.legend()
    
    plt.tight_layout()
    
    plt.show()
    Leakage warning: Rolling statistics used as model features must be shifted so the current target is not included in its own predictor.

    Prevent Temporal Leakage

    Potential Leakage Why It Is Invalid
    Centered rolling mean Uses observations from both past and future
    Unshifted current-period total Uses the target being predicted
    Actual future promotion result Becomes known only after the forecast origin
    Final future price May not be known when the forecast is created
    Full-dataset normalization Uses statistics from validation or test periods
    Random train-test split Allows later observations to influence training
    Revised future actuals Uses corrections unavailable during historical forecast production
    Availability Rule Every feature must be reproducible using only information known at the historical forecast origin.

    Chronological Train, Validation, and Test Split

    series = weekly_sales[
        "sales"
    ].dropna()
    
    
    number_of_periods = len(
        series
    )
    
    
    if number_of_periods < 40:
    
        raise ValueError(
            "The example requires at least 40 complete periods. "
            "Longer history is recommended."
        )
    
    
    training_end = int(
        number_of_periods
        * 0.70
    )
    
    
    validation_end = int(
        number_of_periods
        * 0.85
    )
    
    
    training_series = series.iloc[
        :training_end
    ]
    
    
    validation_series = series.iloc[
        training_end:validation_end
    ]
    
    
    test_series = series.iloc[
        validation_end:
    ]
    
    
    print(
        "Training periods:",
        len(training_series)
    )
    
    print(
        "Validation periods:",
        len(validation_series)
    )
    
    print(
        "Test periods:",
        len(test_series)
    )
    Time-order rule: Training must occur before validation, and validation must occur before testing. Never shuffle the time axis for this project.

    Naïve Forecast Baseline

    A naïve forecast predicts that the next period will equal the most recent observed value.

    Naïve Forecast
    \[ \hat{y}_{t+1\mid t} = y_t \]
    naive_predictions = (
        series
        .shift(
            1
        )
        .loc[
            test_series.index
        ]
    )
    
    
    naive_evaluation = pd.DataFrame(
        {
            "actual": test_series,
            "prediction": naive_predictions
        }
    ).dropna()
    
    
    print(
        naive_evaluation.head()
    )

    Seasonal Naïve Baseline

    A seasonal naïve forecast predicts that the next value will equal the value observed in the corresponding previous seasonal period.

    Seasonal Naïve Forecast
    \[ \hat{y}_{t+h\mid t} = y_{t+h-m} \]

    Here, \(m\) represents the assumed seasonal period. For weekly data, the appropriate value depends on the business calendar and observed periodicity.

    seasonal_period = 52
    
    
    seasonal_naive_predictions = (
        series
        .shift(
            seasonal_period
        )
        .loc[
            test_series.index
        ]
    )
    
    
    seasonal_naive_evaluation = pd.DataFrame(
        {
            "actual": test_series,
            "prediction": (
                seasonal_naive_predictions
            )
        }
    ).dropna()
    
    
    print(
        seasonal_naive_evaluation.head()
    )
    Baseline requirement: Do not assume a seasonal period only because it is common. Confirm it using business calendars, plots, domain knowledge, and validation evidence.

    Moving-Average Baseline

    moving_average_window = 4
    
    
    moving_average_predictions = (
        series
        .shift(
            1
        )
        .rolling(
            window=moving_average_window
        )
        .mean()
        .loc[
            test_series.index
        ]
    )
    
    
    moving_average_evaluation = pd.DataFrame(
        {
            "actual": test_series,
            "prediction": (
                moving_average_predictions
            )
        }
    ).dropna()
    Fair comparison: Every model and baseline must predict the same dates using the same information cutoff.

    Create Lag and Rolling Features

    def create_time_series_features(
        sales_series
    ):
        feature_data = pd.DataFrame(
            {
                "sales": sales_series
            }
        )
    
        feature_data[
            "lag_1"
        ] = feature_data["sales"].shift(
            1
        )
    
        feature_data[
            "lag_2"
        ] = feature_data["sales"].shift(
            2
        )
    
        feature_data[
            "lag_4"
        ] = feature_data["sales"].shift(
            4
        )
    
        feature_data[
            "lag_13"
        ] = feature_data["sales"].shift(
            13
        )
    
        feature_data[
            "lag_52"
        ] = feature_data["sales"].shift(
            52
        )
    
        shifted_sales = feature_data[
            "sales"
        ].shift(
            1
        )
    
        feature_data[
            "rolling_mean_4"
        ] = shifted_sales.rolling(
            window=4
        ).mean()
    
        feature_data[
            "rolling_std_4"
        ] = shifted_sales.rolling(
            window=4
        ).std()
    
        feature_data[
            "rolling_mean_13"
        ] = shifted_sales.rolling(
            window=13
        ).mean()
    
        feature_data[
            "rolling_max_13"
        ] = shifted_sales.rolling(
            window=13
        ).max()
    
        feature_data[
            "week_of_year"
        ] = feature_data.index.isocalendar().week.astype(
            int
        )
    
        feature_data[
            "month"
        ] = feature_data.index.month
    
        feature_data[
            "quarter"
        ] = feature_data.index.quarter
    
        feature_data[
            "year"
        ] = feature_data.index.year
    
        return feature_data.dropna()
    
    
    feature_data = create_time_series_features(
        series
    )
    
    
    print(
        feature_data.head()
    )
    Critical correction: Rolling features use shift(1) before calculating the window. This prevents the current target from leaking into its own predictors.

    Encode Cyclical Calendar Features

    Week numbers are cyclical. Week 52 is close to week 1, even though their ordinary numerical difference is large.

    Cyclical Encoding
    \[ x_{\sin} = \sin \left( \frac{2\pi x}{m} \right) \] \[ x_{\cos} = \cos \left( \frac{2\pi x}{m} \right) \]
    feature_data[
        "week_sin"
    ] = np.sin(
        2
        * np.pi
        * feature_data["week_of_year"]
        / 52
    )
    
    
    feature_data[
        "week_cos"
    ] = np.cos(
        2
        * np.pi
        * feature_data["week_of_year"]
        / 52
    )
    
    
    feature_data = feature_data.drop(
        columns=[
            "week_of_year"
        ]
    )

    Split the Feature Dataset Chronologically

    target_name = "sales"
    
    
    model_features = [
        column_name
        for column_name
        in feature_data.columns
        if column_name != target_name
    ]
    
    
    training_features = feature_data.loc[
        feature_data.index
        <
        validation_series.index.min(),
        model_features
    ]
    
    
    training_target = feature_data.loc[
        training_features.index,
        target_name
    ]
    
    
    validation_features = feature_data.loc[
        (
            feature_data.index
            >=
            validation_series.index.min()
        )
        &
        (
            feature_data.index
            <
            test_series.index.min()
        ),
        model_features
    ]
    
    
    validation_target = feature_data.loc[
        validation_features.index,
        target_name
    ]
    
    
    test_features = feature_data.loc[
        feature_data.index
        >=
        test_series.index.min(),
        model_features
    ]
    
    
    test_target = feature_data.loc[
        test_features.index,
        target_name
    ]
    
    
    print(
        training_features.shape,
        validation_features.shape,
        test_features.shape
    )

    Train a Gradient Boosting Model

    from sklearn.ensemble import HistGradientBoostingRegressor
    
    
    forecast_model = HistGradientBoostingRegressor(
        learning_rate=0.05,
        max_iter=300,
        max_leaf_nodes=15,
        min_samples_leaf=10,
        l2_regularization=0.1,
        random_state=42
    )
    
    
    forecast_model.fit(
        training_features,
        training_target
    )
    
    
    validation_predictions = forecast_model.predict(
        validation_features
    )
    
    
    test_predictions = forecast_model.predict(
        test_features
    )
    
    
    test_predictions = np.maximum(
        test_predictions,
        0.0
    )
    Nonnegative-target rule: Clipping predictions to zero may be suitable for nonnegative unit sales, but it should be documented. It may be inappropriate if the target legitimately includes negative net adjustments.

    Forecast Evaluation Metrics

    Forecast Error
    \[ e_t = y_t-\hat{y}_t \]
    Mean Absolute Error
    \[ MAE = \frac{1}{n} \sum_{t=1}^{n} |e_t| \]
    Root Mean Squared Error
    \[ RMSE = \sqrt{ \frac{1}{n} \sum_{t=1}^{n} e_t^2 } \]
    Weighted Absolute Percentage Error
    \[ WAPE = \frac{ \sum_{t=1}^{n} |y_t-\hat{y}_t| }{ \sum_{t=1}^{n}|y_t| } \times 100 \]
    Forecast Bias
    \[ Bias = \frac{1}{n} \sum_{t=1}^{n} \left( \hat{y}_t-y_t \right) \]

    Reusable Forecast Evaluation

    from sklearn.metrics import mean_absolute_error
    from sklearn.metrics import mean_squared_error
    
    
    def evaluate_forecast(
        model_name,
        actual_values,
        predicted_values
    ):
        actual_array = np.asarray(
            actual_values,
            dtype=float
        )
    
        predicted_array = np.asarray(
            predicted_values,
            dtype=float
        )
    
        mae = mean_absolute_error(
            actual_array,
            predicted_array
        )
    
        rmse = mean_squared_error(
            actual_array,
            predicted_array
        ) ** 0.5
    
        denominator = np.abs(
            actual_array
        ).sum()
    
        wape = (
            np.abs(
                actual_array
                - predicted_array
            ).sum()
            / denominator
            * 100
            if denominator > 0
            else np.nan
        )
    
        bias = np.mean(
            predicted_array
            - actual_array
        )
    
        return {
            "model": model_name,
            "mae": float(mae),
            "rmse": float(rmse),
            "wape_percent": float(wape),
            "bias": float(bias)
        }
    
    
    model_result = evaluate_forecast(
        model_name="Gradient Boosting",
        actual_values=test_target,
        predicted_values=test_predictions
    )
    
    
    print(
        model_result
    )

    MASE for Baseline-Relative Evaluation

    Mean Absolute Scaled Error compares forecast error with the in-sample error of a naïve reference forecast.

    Mean Absolute Scaled Error
    \[ MASE = \frac{ \frac{1}{n} \sum_{t=1}^{n} |y_t-\hat{y}_t| }{ \frac{1}{T-m} \sum_{t=m+1}^{T} |y_t-y_{t-m}| } \]
    def calculate_mase(
        training_values,
        actual_values,
        predicted_values,
        seasonal_period=1
    ):
        training_array = np.asarray(
            training_values,
            dtype=float
        )
    
        actual_array = np.asarray(
            actual_values,
            dtype=float
        )
    
        predicted_array = np.asarray(
            predicted_values,
            dtype=float
        )
    
        scale_errors = np.abs(
            training_array[
                seasonal_period:
            ]
            -
            training_array[
                :-seasonal_period
            ]
        )
    
        scale = scale_errors.mean()
    
        if scale == 0:
    
            return np.nan
    
        forecast_mae = np.abs(
            actual_array
            - predicted_array
        ).mean()
    
        return float(
            forecast_mae
            / scale
        )
    
    
    mase = calculate_mase(
        training_values=training_target,
        actual_values=test_target,
        predicted_values=test_predictions,
        seasonal_period=1
    )
    
    
    print(
        "MASE:",
        mase
    )
    Interpretation: A MASE below one indicates a lower average absolute error than the selected in-sample naïve scaling reference.

    Compare Models and Baselines

    comparison_rows = []
    
    
    comparison_rows.append(
        evaluate_forecast(
            model_name="Naive",
            actual_values=(
                naive_evaluation["actual"]
            ),
            predicted_values=(
                naive_evaluation["prediction"]
            )
        )
    )
    
    
    if not seasonal_naive_evaluation.empty:
    
        comparison_rows.append(
            evaluate_forecast(
                model_name="Seasonal Naive",
                actual_values=(
                    seasonal_naive_evaluation[
                        "actual"
                    ]
                ),
                predicted_values=(
                    seasonal_naive_evaluation[
                        "prediction"
                    ]
                )
            )
        )
    
    
    comparison_rows.append(
        evaluate_forecast(
            model_name="Moving Average",
            actual_values=(
                moving_average_evaluation[
                    "actual"
                ]
            ),
            predicted_values=(
                moving_average_evaluation[
                    "prediction"
                ]
            )
        )
    )
    
    
    comparison_rows.append(
        evaluate_forecast(
            model_name="Gradient Boosting",
            actual_values=test_target,
            predicted_values=test_predictions
        )
    )
    
    
    comparison_table = pd.DataFrame(
        comparison_rows
    ).sort_values(
        by="mae"
    )
    
    
    print(
        comparison_table
    )
    No fabricated result: The article does not claim that the Machine Learning model wins. Select a model only after calculating results using identical forecast dates, horizons, and information cutoffs.

    Plot Actual vs Forecast Sales

    forecast_results = pd.DataFrame(
        {
            "actual": test_target,
            "forecast": test_predictions
        },
        index=test_target.index
    )
    
    
    plt.figure(
        figsize=(12, 6)
    )
    
    
    plt.plot(
        forecast_results.index,
        forecast_results["actual"],
        marker="o",
        label="Actual Sales"
    )
    
    
    plt.plot(
        forecast_results.index,
        forecast_results["forecast"],
        marker="o",
        label="Forecast Sales"
    )
    
    
    plt.title(
        "Actual vs Forecast Sales"
    )
    
    plt.xlabel(
        "Week"
    )
    
    plt.ylabel(
        "Sales"
    )
    
    plt.legend()
    
    plt.tight_layout()
    
    plt.show()

    Residual Analysis

    forecast_results[
        "error"
    ] = (
        forecast_results["actual"]
        - forecast_results["forecast"]
    )
    
    
    forecast_results[
        "absolute_error"
    ] = (
        forecast_results["error"]
        .abs()
    )
    
    
    plt.figure(
        figsize=(12, 5)
    )
    
    
    plt.plot(
        forecast_results.index,
        forecast_results["error"],
        marker="o"
    )
    
    
    plt.axhline(
        y=0,
        linestyle="--"
    )
    
    
    plt.title(
        "Forecast Errors Over Time"
    )
    
    plt.xlabel(
        "Week"
    )
    
    plt.ylabel(
        "Actual Minus Forecast"
    )
    
    plt.tight_layout()
    
    plt.show()
    Error Pattern Possible Concern
    Mostly positive errors The model tends to underforecast
    Mostly negative errors The model tends to overforecast
    Errors increase over time The model may not represent a changing trend
    Large errors in repeated calendar periods Important seasonality may be missing
    Large errors during promotions Promotion information may be missing or misaligned
    Large errors during stockouts Observed sales may not represent unconstrained demand

    Inspect the Largest Errors

    largest_errors = (
        forecast_results
        .sort_values(
            by="absolute_error",
            ascending=False
        )
        .head(20)
    )
    
    
    print(
        largest_errors
    )

    Review High-Error Periods For

    • Promotions and price changes
    • Stockouts and supply interruptions
    • Store closures and holidays
    • One-time bulk orders
    • Returns and accounting corrections
    • Product introductions or discontinuations
    • Competitor or market events
    • Data-load failures
    • Changes in product or store hierarchy
    • Unusual weather or local events where relevant
    Incorrect Practice Delete high-error periods only because they reduce reported model accuracy.
    Correct Practice Determine whether each period reflects bad data, a rare valid event, or a recurring business condition the model should learn.

    Expanding-Window Validation

    Time-series validation trains on early observations and evaluates on later observations. The training window expands as the forecast origin moves forward.

    Expanding Window
    Train 1 → Validate 1
    Train 1 + 2 → Validate 2
    Train 1 + 2 + 3 → Validate 3
    from sklearn.model_selection import TimeSeriesSplit
    from sklearn.metrics import mean_absolute_error
    
    
    time_series_split = TimeSeriesSplit(
        n_splits=5
    )
    
    
    fold_results = []
    
    
    feature_matrix = feature_data[
        model_features
    ]
    
    
    target_vector = feature_data[
        target_name
    ]
    
    
    for fold_number, (
        train_indices,
        validation_indices
    ) in enumerate(
        time_series_split.split(
            feature_matrix
        ),
        start=1
    ):
        fold_model = HistGradientBoostingRegressor(
            learning_rate=0.05,
            max_iter=300,
            max_leaf_nodes=15,
            min_samples_leaf=10,
            l2_regularization=0.1,
            random_state=42
        )
    
        fold_model.fit(
            feature_matrix.iloc[
                train_indices
            ],
            target_vector.iloc[
                train_indices
            ]
        )
    
        fold_predictions = fold_model.predict(
            feature_matrix.iloc[
                validation_indices
            ]
        )
    
        fold_predictions = np.maximum(
            fold_predictions,
            0.0
        )
    
        fold_mae = mean_absolute_error(
            target_vector.iloc[
                validation_indices
            ],
            fold_predictions
        )
    
        fold_results.append(
            {
                "fold": fold_number,
                "training_end": (
                    feature_matrix.index[
                        train_indices[-1]
                    ]
                ),
                "validation_start": (
                    feature_matrix.index[
                        validation_indices[0]
                    ]
                ),
                "validation_end": (
                    feature_matrix.index[
                        validation_indices[-1]
                    ]
                ),
                "mae": float(
                    fold_mae
                )
            }
        )
    
    
    fold_table = pd.DataFrame(
        fold_results
    )
    
    
    print(
        fold_table
    )
    Validation requirement: Cross-validation must preserve chronological order. Ordinary shuffled K-fold validation is inappropriate for forecasting.

    One-Step vs Multi-Step Forecasting

    Strategy Description Important Consideration
    One-Step Predicts only the next period Uses newly observed actuals before the next forecast
    Recursive Feeds earlier predictions into later lag features Errors may accumulate across the horizon
    Direct Trains a separate model for each forecast horizon Requires more models and enough data per horizon
    Multi-Output Predicts several future values together Requires an algorithm supporting multi-output targets

    Recursive Forecasting Helper

    def create_feature_row(
        history,
        forecast_date
    ):
        history_series = pd.Series(
            history
        ).sort_index()
    
        required_history = 52
    
        if len(history_series) < required_history:
    
            raise ValueError(
                "At least 52 historical observations "
                "are required for the configured features."
            )
    
        week_number = int(
            forecast_date.isocalendar().week
        )
    
        return pd.DataFrame(
            [
                {
                    "lag_1": history_series.iloc[-1],
                    "lag_2": history_series.iloc[-2],
                    "lag_4": history_series.iloc[-4],
                    "lag_13": history_series.iloc[-13],
                    "lag_52": history_series.iloc[-52],
                    "rolling_mean_4": (
                        history_series.iloc[-4:].mean()
                    ),
                    "rolling_std_4": (
                        history_series.iloc[-4:].std()
                    ),
                    "rolling_mean_13": (
                        history_series.iloc[-13:].mean()
                    ),
                    "rolling_max_13": (
                        history_series.iloc[-13:].max()
                    ),
                    "month": forecast_date.month,
                    "quarter": forecast_date.quarter,
                    "year": forecast_date.year,
                    "week_sin": np.sin(
                        2
                        * np.pi
                        * week_number
                        / 52
                    ),
                    "week_cos": np.cos(
                        2
                        * np.pi
                        * week_number
                        / 52
                    )
                }
            ],
            index=[
                forecast_date
            ]
        )
    
    
    def recursive_forecast(
        model,
        historical_series,
        number_of_periods,
        frequency="W-SUN"
    ):
        history = historical_series.copy()
    
        future_dates = pd.date_range(
            start=(
                history.index.max()
                + pd.tseries.frequencies.to_offset(
                    frequency
                )
            ),
            periods=number_of_periods,
            freq=frequency
        )
    
        forecast_values = []
    
        for forecast_date in future_dates:
    
            feature_row = create_feature_row(
                history=history,
                forecast_date=forecast_date
            )
    
            prediction = float(
                model.predict(
                    feature_row[
                        model_features
                    ]
                )[0]
            )
    
            prediction = max(
                prediction,
                0.0
            )
    
            history.loc[
                forecast_date
            ] = prediction
    
            forecast_values.append(
                prediction
            )
    
        return pd.Series(
            forecast_values,
            index=future_dates,
            name="forecast"
        )

    Generate a Future Forecast

    future_forecast = recursive_forecast(
        model=forecast_model,
        historical_series=series,
        number_of_periods=8,
        frequency="W-SUN"
    )
    
    
    print(
        future_forecast
    )
    
    
    plt.figure(
        figsize=(12, 6)
    )
    
    
    plt.plot(
        series.tail(26).index,
        series.tail(26).values,
        marker="o",
        label="Historical Sales"
    )
    
    
    plt.plot(
        future_forecast.index,
        future_forecast.values,
        marker="o",
        label="Future Forecast"
    )
    
    
    plt.title(
        "Historical and Forecast Sales"
    )
    
    plt.xlabel(
        "Week"
    )
    
    plt.ylabel(
        "Sales"
    )
    
    plt.legend()
    
    plt.tight_layout()
    
    plt.show()
    Recursive-risk warning: Later forecasts depend on earlier predictions. Evaluate accuracy separately for every horizon because errors can grow over time.

    Evaluate by Forecast Horizon

    Horizon Meaning Typical Use
    \(h=1\) Next period Immediate operational response
    \(h=2\) Two periods ahead Short lead-time planning
    \(h=4\) Four periods ahead Monthly operational planning for weekly data
    Long Horizon Several future periods Capacity, procurement, or budgeting decisions
    Reporting requirement: An overall average can hide weak long-horizon performance. Report error separately by horizon whenever multiple future periods are forecast.

    Forecast Intervals

    A point forecast provides one value. A forecast interval communicates a range of plausible outcomes under the assumptions of the chosen interval method.

    Generic Forecast Interval
    \[ \left[ \hat{y}_{t+h}^{L}, \hat{y}_{t+h}^{U} \right] \]

    The following educational example creates an empirical interval using validation residual quantiles.

    validation_errors = (
        validation_target.to_numpy()
        - validation_predictions
    )
    
    
    lower_error_quantile = np.quantile(
        validation_errors,
        0.05
    )
    
    
    upper_error_quantile = np.quantile(
        validation_errors,
        0.95
    )
    
    
    forecast_interval = pd.DataFrame(
        {
            "forecast": future_forecast,
            "lower_bound": np.maximum(
                future_forecast
                + lower_error_quantile,
                0.0
            ),
            "upper_bound": np.maximum(
                future_forecast
                + upper_error_quantile,
                0.0
            )
        }
    )
    
    
    print(
        forecast_interval
    )
    Interval limitation: This simple residual-quantile interval assumes validation errors are representative. It may not capture horizon-specific uncertainty, changing variance, or structural breaks.

    Plot the Forecast Interval

    plt.figure(
        figsize=(12, 6)
    )
    
    
    plt.plot(
        series.tail(26).index,
        series.tail(26).values,
        label="Historical Sales"
    )
    
    
    plt.plot(
        forecast_interval.index,
        forecast_interval["forecast"],
        marker="o",
        label="Forecast"
    )
    
    
    plt.fill_between(
        forecast_interval.index,
        forecast_interval["lower_bound"],
        forecast_interval["upper_bound"],
        alpha=0.25,
        label="Empirical Forecast Interval"
    )
    
    
    plt.title(
        "Sales Forecast with Interval"
    )
    
    plt.xlabel(
        "Week"
    )
    
    plt.ylabel(
        "Sales"
    )
    
    plt.legend()
    
    plt.tight_layout()
    
    plt.show()

    Hierarchical Forecasting

    Real sales data often contains several related aggregation levels, such as products within categories, stores within regions, and regions within a total business.

    Aggregation Consistency
    \[ \text{Total Forecast} = \sum_{i=1}^{k} \text{Component Forecast}_i \]
    Approach General Idea
    Bottom-Up Forecast detailed series and aggregate them
    Top-Down Forecast the total and allocate it to lower levels
    Middle-Out Forecast an intermediate level and aggregate or distribute
    Reconciliation Adjust forecasts so related hierarchy levels remain coherent
    Planning requirement: Product, store, category, region, and total forecasts should not contradict one another without an explicit reconciliation policy.

    Intermittent and New-Product Sales

    Situation Challenge Possible Response
    Intermittent Sales Many zero periods with occasional positive demand Evaluate specialized intermittent-demand methods
    New Product No target history exists Use similar products, attributes, plans, or human judgment
    Discontinued Product Historical pattern no longer continues Represent lifecycle status explicitly
    Replacement Product Demand transfers from an older item Model substitution and transition where possible
    Promotion Launch Historical baseline may not represent promoted demand Use promotion plans available at forecast time

    External and Planned Features

    Forecasting models may use prices, promotions, holidays, weather, store openings, marketing plans, and other variables only when their values are available or forecastable at the historical and future prediction times.

    Feature Availability Question
    Promotion Was the promotion plan finalized at the forecast origin?
    Price Was the future price approved or only observed later?
    Holiday Was the calendar known in advance?
    Weather Is the future value an actual observation or a forecast?
    Inventory Is the value current stock, planned supply, or future actual stock?
    Marketing Spend Was the plan committed before the forecast was created?
    Backtest Leakage Use future actual weather, prices, inventory, or promotion execution while claiming the model used information available beforehand.
    Correct Backtest Use the value or forecast that would genuinely have been available at each historical forecast origin.

    Save the Forecasting Model

    import json
    import joblib
    import sklearn
    
    
    model_path = (
        "sales_forecasting_model.joblib"
    )
    
    
    metadata_path = (
        "sales_forecasting_metadata.json"
    )
    
    
    joblib.dump(
        forecast_model,
        model_path
    )
    
    
    metadata = {
        "project": "Sales Forecasting",
        "target": "sales",
        "frequency": "W-SUN",
        "forecast_strategy": "recursive",
        "model_features": model_features,
        "minimum_history": 52,
        "forecast_model": (
            "HistGradientBoostingRegressor"
        ),
        "scikit_learn_version": (
            sklearn.__version__
        ),
        "test_metrics": model_result,
        "training_end": str(
            training_target.index.max()
        ),
        "test_end": str(
            test_target.index.max()
        )
    }
    
    
    with open(
        metadata_path,
        "w",
        encoding="utf-8"
    ) as metadata_file:
    
        json.dump(
            metadata,
            metadata_file,
            indent=4
        )
    
    
    print(
        "Saved model:",
        model_path
    )
    
    print(
        "Saved metadata:",
        metadata_path
    )
    Artifact requirement: Preserve the model, feature definitions, time frequency, horizon, aggregation rules, target definition, training cutoff, library versions, metrics, and known limitations.

    Reload and Forecast

    loaded_model = joblib.load(
        "sales_forecasting_model.joblib"
    )
    
    
    with open(
        "sales_forecasting_metadata.json",
        "r",
        encoding="utf-8"
    ) as metadata_file:
    
        loaded_metadata = json.load(
            metadata_file
        )
    
    
    loaded_forecast = recursive_forecast(
        model=loaded_model,
        historical_series=series,
        number_of_periods=8,
        frequency=loaded_metadata[
            "frequency"
        ]
    )
    
    
    print(
        loaded_forecast
    )
    Artifact security: Load serialized model files only from trusted sources. Preserve the artifact checksum, owner, approval status, data version, and dependency versions.

    Validate Forecast Input

    def validate_sales_history(
        sales_series,
        expected_frequency,
        minimum_history
    ):
        if not isinstance(
            sales_series,
            pd.Series
        ):
            raise TypeError(
                "sales_series must be a pandas Series."
            )
    
        if not isinstance(
            sales_series.index,
            pd.DatetimeIndex
        ):
            raise TypeError(
                "sales_series must use a DatetimeIndex."
            )
    
        if sales_series.index.has_duplicates:
    
            raise ValueError(
                "The sales index contains duplicate dates."
            )
    
        if not sales_series.index.is_monotonic_increasing:
    
            raise ValueError(
                "The sales index must be sorted chronologically."
            )
    
        if sales_series.isna().any():
    
            raise ValueError(
                "The sales history contains missing values."
            )
    
        if len(sales_series) < minimum_history:
    
            raise ValueError(
                "Insufficient history for the configured features."
            )
    
        expected_index = pd.date_range(
            start=sales_series.index.min(),
            end=sales_series.index.max(),
            freq=expected_frequency
        )
    
        missing_dates = expected_index.difference(
            sales_series.index
        )
    
        if len(missing_dates) > 0:
    
            raise ValueError(
                "The time series contains missing periods."
            )
    
        return sales_series.astype(
            float
        )
    
    
    validated_series = validate_sales_history(
        sales_series=series,
        expected_frequency="W-SUN",
        minimum_history=52
    )

    Forecasts Are Conditional Estimates

    A sales forecast is conditional on the historical data, feature assumptions, model, business environment, and future conditions. It is not a guaranteed outcome.

    Unsafe Interpretation Treat the point forecast as the exact amount that will be sold.
    Responsible Interpretation Present the point forecast with its horizon, interval, assumptions, historical cutoff, model version, and known limitations.

    Human Adjustments

    Authorized planners may know about future events that are absent from the model. Adjustments should be recorded separately rather than silently overwriting the statistical forecast.

    Forecast Value Meaning
    Statistical Forecast Output generated by the forecasting model
    Human Adjustment Documented change based on additional approved information
    Final Consensus Forecast Approved value used for planning
    Actual Outcome Observed value after the period is complete
    Audit principle: Store the model forecast, adjustment amount, reason, approver, timestamp, final forecast, and actual outcome separately.

    Forecast Value Added

    Forecast Value Added compares the error before and after a process, adjustment, or model stage.

    Forecast Value Added
    \[ FVA = Error_{\text{Before}} - Error_{\text{After}} \]
    FVA Result Interpretation
    Positive The later process reduced the selected error measure
    Zero The later process did not change the selected error measure
    Negative The later process increased the selected error measure

    Common Project Mistakes

    Mistake Randomly shuffle time-series rows before splitting.
    Correction Preserve chronological order in training, validation, and testing.
    Mistake Use current-period sales inside an unshifted rolling feature.
    Correction Shift the series before calculating rolling statistics.
    Mistake Treat all missing periods as zero sales.
    Correction Distinguish confirmed zero demand, closure, stockout, and missing data.
    Mistake Compare models using different dates or forecast horizons.
    Correction Use identical forecast origins, horizons, targets, and data cutoffs.
    Mistake Ignore naïve and seasonal naïve benchmarks.
    Correction Require sophisticated models to beat strong simple baselines.
    Mistake Report only one average error.
    Correction Report scale-dependent error, percentage or scaled error, bias, horizon performance, and interval coverage.
    Mistake Assume historical sales represent unconstrained demand.
    Correction Investigate stockouts, capacity restrictions, fulfillment failures, and product availability.
    Mistake Overwrite model forecasts with undocumented manual values.
    Correction Store model forecasts and human adjustments separately.

    Responsible Forecasting

    Responsible-Use Checklist

    • Document the intended planning decision
    • Protect confidential commercial and customer data
    • Use only features available at the forecast origin
    • Evaluate errors across products, stores, and regions
    • Review whether low-volume groups receive unstable forecasts
    • Communicate uncertainty and assumptions
    • Preserve model and human-adjustment histories
    • Provide review for consequential planning decisions
    • Document known cold-start and stockout limitations
    • Monitor changing market conditions
    • Define fallback forecasts
    • Maintain rollback and incident-response procedures

    Production Monitoring

    Monitoring Area Example Measure
    Data Freshness Delay between source availability and forecast generation
    Missing Periods Expected dates absent from the input series
    Feature Availability Missing promotions, prices, holidays, or stock information
    Forecast Accuracy MAE, RMSE, WAPE, and MASE after actuals mature
    Forecast Bias Persistent overforecasting or underforecasting
    Horizon Error Error for each forecast distance
    Interval Coverage Proportion of actual values inside forecast intervals
    Group Performance Error by product, store, region, category, and volume band
    Forecast Overrides Frequency, direction, size, and result of manual adjustments
    Structural Change New products, closures, pricing changes, or market shifts
    Pipeline Reliability Failures, latency, schema changes, and incomplete outputs

    Recommended Project Workflow

    1

    Define the Forecast

    Specify the target, units, granularity, frequency, origin, horizon, users, and operational decision.

    2

    Audit and Version the Data

    Review source systems, dates, duplicates, missing periods, returns, stockouts, units, and hierarchy changes.

    3

    Create a Regular Time Series

    Aggregate transactions consistently and distinguish missing observations from confirmed zero values.

    4

    Explore Time-Series Structure

    Inspect trend, seasonality, volatility, anomalies, closures, and structural breaks.

    5

    Create Strong Baselines

    Evaluate naïve, seasonal naïve, and moving-average forecasts.

    6

    Engineer Leakage-Safe Features

    Create shifted lags, rolling statistics, calendar variables, and genuinely available external features.

    7

    Validate Chronologically

    Use expanding-window or rolling-origin validation and preserve an independent final test period.

    8

    Evaluate and Diagnose

    Compare errors, bias, horizon performance, intervals, groups, and important failure periods.

    9

    Save Complete Artifacts

    Preserve the model, features, target, aggregation, frequency, cutoff, metrics, and dependency versions.

    10

    Monitor and Improve

    Track data freshness, drift, error, bias, overrides, interval coverage, planning impact, and model reliability.

    Hands-On Practice

    1

    Define a Forecasting Contract

    Document the target, frequency, horizon, forecast origin, data cutoff, and business decision.

    2

    Audit Missing Periods

    Distinguish confirmed zeros, closures, stockouts, and missing data.

    3

    Compare Three Baselines

    Evaluate naïve, seasonal naïve, and moving-average forecasts.

    4

    Create Lag Features

    Build shifted lag and rolling-window features without target leakage.

    5

    Run Walk-Forward Validation

    Compare model performance across several historical forecast origins.

    6

    Analyze High-Error Periods

    Review promotions, holidays, stockouts, closures, bulk orders, and data corrections.

    7

    Generate Multi-Step Forecasts

    Produce recursive forecasts and compare errors by horizon.

    8

    Track Forecast Value Added

    Compare baseline, model, adjusted, and final forecast errors.

    Knowledge Check

    1

    Why must a sales forecast define its horizon?

    Explain how forecast difficulty and business use change with distance into the future.

    2

    Why should time-series rows not be shuffled?

    Explain chronological order and future-information leakage.

    3

    What is a naïve forecast?

    Explain how it uses the most recent observation as the next prediction.

    4

    What is a seasonal naïve forecast?

    Explain how it uses the corresponding value from an earlier seasonal period.

    5

    Why must rolling features be shifted?

    Explain how an unshifted window can contain the target being predicted.

    6

    What does forecast bias measure?

    Explain persistent overforecasting or underforecasting.

    7

    Why may observed sales underestimate demand?

    Explain stockouts, capacity limits, fulfillment failures, and availability constraints.

    8

    Why should error be evaluated by horizon?

    Explain how longer-range forecasts can be less accurate than near-term forecasts.

    9

    What does a forecast interval communicate?

    Explain the range of plausible outcomes under the interval method's assumptions.

    10

    Why should human adjustments be recorded separately?

    Explain auditability, forecasting improvement, and Forecast Value Added analysis.

    Project Summary

    Project Stage Main Outcome
    Problem Definition Defines the target, frequency, horizon, origin, and business decision
    Data Governance Documents sources, units, timing, corrections, and permitted use
    Aggregation Converts transactions to a consistent time frequency
    Data Audit Checks missing periods, duplicates, returns, and unusual values
    Time-Series Exploration Identifies trend, seasonality, volatility, and structural change
    Baselines Establishes naïve, seasonal, and moving-average references
    Feature Engineering Creates shifted lags, rolling statistics, and calendar features
    Chronological Validation Measures performance without training on future observations
    Modeling Fits a nonlinear regression model to time-safe features
    Evaluation Measures MAE, RMSE, WAPE, MASE, bias, and horizon error
    Forecast Intervals Communicates uncertainty around point predictions
    Persistence Saves the model and complete forecasting metadata
    Monitoring Tracks freshness, drift, error, bias, coverage, and overrides

    Key Takeaway

    A reliable sales forecasting project begins with a precise target, time frequency, forecast origin, horizon, and planning decision. Historical data must be aggregated consistently, missing periods must be distinguished from true zero sales, and every feature must be available at the forecast creation time. Strong naïve baselines, chronological validation, shifted lag and rolling features, multiple error metrics, horizon-level analysis, forecast intervals, and high-error investigation are essential. Production use also requires saved forecasting metadata, documented human adjustments, monitoring for bias and drift, and clear communication that forecasts are conditional estimates rather than guaranteed outcomes.

    Get Machine Learning updates by email

    Short, practical lessons and course updates from Rumman Ansari. No spam.