Controlled comparisons: changing the forecast or the policy with run_comparison¶

An inventory policy has two separate inputs: the target it is given (here, a forecast quantile) and the rule that turns the target into orders. To learn what each one does, change one at a time and keep everything else fixed.

SimulationEngine.run_comparison is built for this. It takes several fitted policies and runs each one on the same demand, from its own copy of the same opening state, with the same windows. Any difference between the branches therefore comes from the policies alone.

This notebook shows:

• how targets from an external forecaster enter Stockcast (here Smooth's ETS models; any forecasting library works);
• one run_comparison call with five branches that vary the forecast probability, the forecast model, and the policy rule;
• how to read controlled contrasts from the results, pooled and per item;
• what the two rules do differently, in the event table.

The data are synthetic: five items with Poisson demand.

In [1]:

import numpy as np
import pandas as pd

# Show every number with two decimals.
pd.options.display.float_format = "{:.2f}".format

import numpy as np import pandas as pd # Show every number with two decimals. pd.options.display.float_format = "{:.2f}".format

1. Demand: history for the forecaster, the future for the simulation¶

The forecast origin is 28 February 2026. The 56 days up to and including it are history: the forecaster sees only these. The 42 days after it are the future demand that the simulation replays. Each item has its own Poisson rate.

In [2]:

origin = pd.Timestamp("2026-02-28")
rates = {"tea": 9.0, "coffee": 7.0, "pasta": 6.0, "cereal": 5.0, "soup": 4.0}
skus = list(rates)
n_history, n_days = 56, 42

rng = np.random.default_rng(6)
history_dates = pd.date_range(end=origin, periods=n_history, freq="D")
future_dates = pd.date_range(origin + pd.Timedelta(days=1), periods=n_days, freq="D")
history = pd.DataFrame(
    [
        {"unique_id": sku, "date": date, "y": float(rng.poisson(rate))}
        for sku, rate in rates.items()
        for date in history_dates
    ]
)
demand = pd.DataFrame(
    [
        {
            "unique_id": sku,
            "period": period,
            "date": date,
            "y": float(rng.poisson(rate)),
        }
        for sku, rate in rates.items()
        for period, date in enumerate(future_dates)
    ]
)

origin = pd.Timestamp("2026-02-28") rates = {"tea": 9.0, "coffee": 7.0, "pasta": 6.0, "cereal": 5.0, "soup": 4.0} skus = list(rates) n_history, n_days = 56, 42 rng = np.random.default_rng(6) history_dates = pd.date_range(end=origin, periods=n_history, freq="D") future_dates = pd.date_range(origin + pd.Timedelta(days=1), periods=n_days, freq="D") history = pd.DataFrame( [ {"unique_id": sku, "date": date, "y": float(rng.poisson(rate))} for sku, rate in rates.items() for date in history_dates ] ) demand = pd.DataFrame( [ { "unique_id": sku, "period": period, "date": date, "y": float(rng.poisson(rate)), } for sku, rate in rates.items() for period, date in enumerate(future_dates) ] )

In [3]:

import matplotlib.dates as mdates
import matplotlib.pyplot as plt

HISTORY_COLOR = "#9CA3AF"
FUTURE_COLOR = "#4C78A8"

fig, axes = plt.subplots(
    len(skus),
    1,
    figsize=(10, 7.5),
    sharex=True,
    sharey=True,
    layout="constrained",
)
for ax, sku in zip(axes, skus):
    past = history.loc[history["unique_id"].eq(sku)]
    future = demand.loc[demand["unique_id"].eq(sku)]
    ax.plot(
        past["date"],
        past["y"],
        marker="o",
        markersize=2,
        linewidth=1,
        color=HISTORY_COLOR,
        label="History",
    )
    ax.plot(
        future["date"],
        future["y"],
        marker="o",
        markersize=2,
        linewidth=1,
        color=FUTURE_COLOR,
        label="Future demand (simulated)",
    )
    ax.axvline(
        origin, color="black", linestyle="--", linewidth=1, label="Forecast origin"
    )
    ax.set_ylabel(sku, rotation=0, ha="right", va="center")
    ax.set_ylim(bottom=0)
    ax.grid(axis="y", alpha=0.2)
    ax.spines[["top", "right"]].set_visible(False)
axes[-1].xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
axes[0].legend(frameon=False, ncol=3, loc="lower left", bbox_to_anchor=(0, 1.0))
fig.suptitle("Daily demand per item")
plt.show()

import matplotlib.dates as mdates import matplotlib.pyplot as plt HISTORY_COLOR = "#9CA3AF" FUTURE_COLOR = "#4C78A8" fig, axes = plt.subplots( len(skus), 1, figsize=(10, 7.5), sharex=True, sharey=True, layout="constrained", ) for ax, sku in zip(axes, skus): past = history.loc[history["unique_id"].eq(sku)] future = demand.loc[demand["unique_id"].eq(sku)] ax.plot( past["date"], past["y"], marker="o", markersize=2, linewidth=1, color=HISTORY_COLOR, label="History", ) ax.plot( future["date"], future["y"], marker="o", markersize=2, linewidth=1, color=FUTURE_COLOR, label="Future demand (simulated)", ) ax.axvline( origin, color="black", linestyle="--", linewidth=1, label="Forecast origin" ) ax.set_ylabel(sku, rotation=0, ha="right", va="center") ax.set_ylim(bottom=0) ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) axes[-1].xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) axes[0].legend(frameon=False, ncol=3, loc="lower left", bbox_to_anchor=(0, 1.0)) fig.suptitle("Daily demand per item") plt.show()

2. Targets from an external forecaster¶

Stockcast does not fit forecasts; it takes the target a forecaster produces. With a lead time of L = 2 days and a review every R = 2 days, an order must protect L + R = 4 days of demand, so the target is a quantile of the 4-day total demand. Smooth returns it directly with cumulative=True. It is a quantile of the total, never a sum of four daily quantiles.

We fit two ETS models per item on history only: ANN (level only) and AAN (level and trend). From ANN we take the 90% and 95% quantiles, from AAN the 95% quantile. Orders are in whole units, so each target is rounded up. Every target is dated: it covers the four days after the origin, up to the end date.

In [4]:

import math

from smooth import ES

lead_time, review_period = 2, 2
horizon = lead_time + review_period
forecast_inputs = [("ANN", 0.90), ("ANN", 0.95), ("AAN", 0.95)]

target_rows = []
for sku in skus:
    y = history.loc[history["unique_id"].eq(sku), "y"].reset_index(drop=True)
    for model_name in ["ANN", "AAN"]:
        model = ES(model_name)
        model.fit(y)
        for probability in [p for m, p in forecast_inputs if m == model_name]:
            forecast = model.predict(
                h=horizon,
                interval="prediction",
                level=probability,
                side="upper",
                cumulative=True,
            )
            target_rows.append(
                {
                    "unique_id": sku,
                    "model": model_name,
                    "probability": probability,
                    "target": float(math.ceil(forecast.upper.iloc[0, 0])),
                    "end": origin + pd.Timedelta(days=horizon),
                }
            )
targets = pd.DataFrame(target_rows)

import math from smooth import ES lead_time, review_period = 2, 2 horizon = lead_time + review_period forecast_inputs = [("ANN", 0.90), ("ANN", 0.95), ("AAN", 0.95)] target_rows = [] for sku in skus: y = history.loc[history["unique_id"].eq(sku), "y"].reset_index(drop=True) for model_name in ["ANN", "AAN"]: model = ES(model_name) model.fit(y) for probability in [p for m, p in forecast_inputs if m == model_name]: forecast = model.predict( h=horizon, interval="prediction", level=probability, side="upper", cumulative=True, ) target_rows.append( { "unique_id": sku, "model": model_name, "probability": probability, "target": float(math.ceil(forecast.upper.iloc[0, 0])), "end": origin + pd.Timedelta(days=horizon), } ) targets = pd.DataFrame(target_rows)

In [5]:

# The 4-day targets, one column per forecast input.
targets.assign(
    forecast=targets["model"] + " " + targets["probability"].map("{:.0%}".format)
).pivot(index="unique_id", columns="forecast", values="target").loc[
    skus, ["ANN 90%", "ANN 95%", "AAN 95%"]
]

# The 4-day targets, one column per forecast input. targets.assign( forecast=targets["model"] + " " + targets["probability"].map("{:.0%}".format) ).pivot(index="unique_id", columns="forecast", values="target").loc[ skus, ["ANN 90%", "ANN 95%", "AAN 95%"] ]

Out[5]:

forecast	ANN 90%	ANN 95%	AAN 95%
unique_id
tea	46.00	48.00	47.00
coffee	36.00	38.00	40.00
pasta	33.00	35.00	36.00
cereal	27.00	29.00	29.00
soup	20.00	22.00	25.00

What to notice: raising the probability from 90% to 95% raises every ANN target, by about two units (tea from 46 to 48). Changing the model moves targets in both directions: AAN gives higher targets for coffee (40 against 38), pasta and soup (25 against 22), a lower one for tea (47 against 48), and the same for cereal.

3. Two rules, fed by the same targets¶

The same target plays a different role in each rule:

• Order-up-to (OrderUpToPolicy): at each review, order the gap between the target S and the inventory position (stock on hand plus stock on order). The order size follows demand.
• Reorder point with a fixed quantity, (s,Q) (ReorderPointPolicy with order_quantity): at each review, if the inventory position is at or below the target s, order exactly Q = 15 units. The order size never changes.

Each policy states the probability of its target with service_level, and checks that the target's end date matches its own L + R window. We fit five policies: order-up-to with each of the three forecast inputs, and (s,Q) with the two 95% inputs.

In [6]:

from stockcast.policies import OrderUpToPolicy, ReorderPointPolicy


def target_table(model_name, probability):
    rows = targets["model"].eq(model_name) & targets["probability"].eq(probability)
    return targets.loc[rows]


fixed_quantity = 15.0
policies = {}
for model_name, probability in forecast_inputs:
    policy = OrderUpToPolicy(
        lead_time=lead_time,
        review_period=review_period,
        freq="D",
        service_level=probability,
        allow_backorders=False,
        date_column="end",
    )
    policy.fit(
        target_table(model_name, probability),
        target_column="target",
        forecast_origin=origin,
    )
    policies[f"Order-up-to, {model_name} {probability:.0%}"] = policy

for model_name in ["ANN", "AAN"]:
    policy = ReorderPointPolicy(
        lead_time=lead_time,
        review_period=review_period,
        freq="D",
        service_level=0.95,
        order_quantity=fixed_quantity,
        allow_backorders=False,
        date_column="end",
    )
    policy.fit(
        target_table(model_name, 0.95),
        reorder_point_column="target",
        forecast_origin=origin,
    )
    policies[f"(s,Q), {model_name} 95%"] = policy

from stockcast.policies import OrderUpToPolicy, ReorderPointPolicy def target_table(model_name, probability): rows = targets["model"].eq(model_name) & targets["probability"].eq(probability) return targets.loc[rows] fixed_quantity = 15.0 policies = {} for model_name, probability in forecast_inputs: policy = OrderUpToPolicy( lead_time=lead_time, review_period=review_period, freq="D", service_level=probability, allow_backorders=False, date_column="end", ) policy.fit( target_table(model_name, probability), target_column="target", forecast_origin=origin, ) policies[f"Order-up-to, {model_name} {probability:.0%}"] = policy for model_name in ["ANN", "AAN"]: policy = ReorderPointPolicy( lead_time=lead_time, review_period=review_period, freq="D", service_level=0.95, order_quantity=fixed_quantity, allow_backorders=False, date_column="end", ) policy.fit( target_table(model_name, 0.95), reorder_point_column="target", forecast_origin=origin, ) policies[f"(s,Q), {model_name} 95%"] = policy

4. One comparison, five branches¶

run_comparison takes the policies as {label: policy}, the demand and the opening state. Every branch starts from the same stock, three days of average demand per item and nothing on order, and replays the same 42 days. Unmet demand is lost.

In [7]:

from stockcast.core import InventoryStateDataFrame, SimulationEngine

opening = InventoryStateDataFrame.from_observed(
    pd.DataFrame({"unique_id": skus, "on_hand": [3 * rates[sku] for sku in skus]}),
    opening_date=origin,
)
comparison = SimulationEngine().run_comparison(policies, demand, opening)

from stockcast.core import InventoryStateDataFrame, SimulationEngine opening = InventoryStateDataFrame.from_observed( pd.DataFrame({"unique_id": skus, "on_hand": [3 * rates[sku] for sku in skus]}), opening_date=origin, ) comparison = SimulationEngine().run_comparison(policies, demand, opening)

Each branch is scored with stockcast.evaluation on the same costs: holding 0.20 per unit and day, a lost sale 5.00 per unit, and 1.50 per order line (one item on one order). Because only the policies differ, pairs of branches that differ in one input are controlled contrasts:

• probability: order-up-to ANN 90% against ANN 95%;
• model: ANN 95% against AAN 95%, within either rule;
• rule: order-up-to against (s,Q), with the same forecast input.

In [8]:

from stockcast.evaluation import (
    InventoryEvaluator,
    avg_on_hand,
    fill_rate,
    lost_sales_units,
    total_cost,
)

costs = {
    "cost_components": ["holding", "shortage", "ordering"],
    "holding_cost_per_unit_period": 0.20,
    "shortage_cost_per_unit": 5.00,
    "order_cost_per_sku_line": 1.50,
    "order_cost_per_unit": 0.0,
}
score_rows = []
for label in policies:
    evaluator = InventoryEvaluator()
    evaluator.fit(comparison[label])
    scores = evaluator.evaluate(
        [lost_sales_units, avg_on_hand, total_cost, fill_rate],
        groupby=["unique_id"],
        context=costs,
    )
    score_rows.append(scores.assign(branch=label))
per_item = pd.concat(score_rows, ignore_index=True)

pooled = per_item.groupby("branch", sort=False).agg(
    lost_sales_units=("lost_sales_units", "sum"),
    avg_on_hand=("avg_on_hand", "mean"),
    total_cost=("total_cost", "sum"),
)

BRANCH_COLORS = ["#72B7B2", "#4C78A8", "#F58518", "#B6D7F0", "#FBD3A6"]
fig, axes = plt.subplots(1, 3, figsize=(13, 3.4), sharey=True, layout="constrained")
for ax, column, title, label_format in [
    (axes[0], "lost_sales_units", "Lost sales (units)", "{:.0f}"),
    (axes[1], "avg_on_hand", "Average stock per item and day", "{:.1f}"),
    (axes[2], "total_cost", "Total cost", "{:.0f}"),
]:
    ax.barh(pooled.index, pooled[column], color=BRANCH_COLORS)
    for y, value in enumerate(pooled[column]):
        ax.text(value, y, " " + label_format.format(value), va="center", fontsize=8)
    ax.set(title=title)
    ax.set_xlim(0, pooled[column].max() * 1.18)
    ax.grid(axis="x", alpha=0.2)
    ax.spines[["top", "right"]].set_visible(False)
axes[0].invert_yaxis()
plt.show()

from stockcast.evaluation import ( InventoryEvaluator, avg_on_hand, fill_rate, lost_sales_units, total_cost, ) costs = { "cost_components": ["holding", "shortage", "ordering"], "holding_cost_per_unit_period": 0.20, "shortage_cost_per_unit": 5.00, "order_cost_per_sku_line": 1.50, "order_cost_per_unit": 0.0, } score_rows = [] for label in policies: evaluator = InventoryEvaluator() evaluator.fit(comparison[label]) scores = evaluator.evaluate( [lost_sales_units, avg_on_hand, total_cost, fill_rate], groupby=["unique_id"], context=costs, ) score_rows.append(scores.assign(branch=label)) per_item = pd.concat(score_rows, ignore_index=True) pooled = per_item.groupby("branch", sort=False).agg( lost_sales_units=("lost_sales_units", "sum"), avg_on_hand=("avg_on_hand", "mean"), total_cost=("total_cost", "sum"), ) BRANCH_COLORS = ["#72B7B2", "#4C78A8", "#F58518", "#B6D7F0", "#FBD3A6"] fig, axes = plt.subplots(1, 3, figsize=(13, 3.4), sharey=True, layout="constrained") for ax, column, title, label_format in [ (axes[0], "lost_sales_units", "Lost sales (units)", "{:.0f}"), (axes[1], "avg_on_hand", "Average stock per item and day", "{:.1f}"), (axes[2], "total_cost", "Total cost", "{:.0f}"), ]: ax.barh(pooled.index, pooled[column], color=BRANCH_COLORS) for y, value in enumerate(pooled[column]): ax.text(value, y, " " + label_format.format(value), va="center", fontsize=8) ax.set(title=title) ax.set_xlim(0, pooled[column].max() * 1.18) ax.grid(axis="x", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) axes[0].invert_yaxis() plt.show()

What to notice:

• Probability (order-up-to, ANN 90% → 95%): lost sales fall from 5 to 2 units, while average stock rises from 11.3 to 13.2 and total cost from 657 to 720. A higher quantile buys service with stock.
• Model (ANN 95% → AAN 95%): under order-up-to, lost sales fall from 2 to 0 for 0.9 more units of average stock. Under (s,Q), lost sales stay at 41: the higher reorder points add stock and cost (1,004 → 1,034) without serving more demand.
• Rule (order-up-to → (s,Q), both ANN 95%): lost sales rise from 2 to 41 units and average stock rises from 13.2 to 16.0. The fixed quantity is short for some items and too large for others; the next two sections show where.

5. Pooled and per-item service¶

The totals above pool all five items. A pooled number can hide one item that is served badly, so the fill rate (fulfilled demand as a share of demand) is shown per item. groupby=["unique_id"] in evaluate gives one row per item.

In [9]:

from IPython.display import display

# Fill rates sit close to 1, so this table shows three decimals.
fill_rates = per_item.pivot(index="branch", columns="unique_id", values="fill_rate")
with pd.option_context("display.float_format", "{:.3f}".format):
    display(fill_rates.loc[list(policies), skus])

from IPython.display import display # Fill rates sit close to 1, so this table shows three decimals. fill_rates = per_item.pivot(index="branch", columns="unique_id", values="fill_rate") with pd.option_context("display.float_format", "{:.3f}".format): display(fill_rates.loc[list(policies), skus])

unique_id	tea	coffee	pasta	cereal	soup
branch
Order-up-to, ANN 90%	1.000	0.986	1.000	1.000	0.994
Order-up-to, ANN 95%	1.000	0.993	1.000	1.000	1.000
Order-up-to, AAN 95%	1.000	1.000	1.000	1.000	1.000
(s,Q), ANN 95%	0.898	0.986	1.000	1.000	1.000
(s,Q), AAN 95%	0.898	0.986	1.000	1.000	1.000

All three order-up-to branches serve at least 98.6% of every item's demand. Both (s,Q) branches serve 89.8% of tea demand while pasta, cereal and soup are fully served: the pooled fill rate hides that one item carries almost all the lost sales.

6. The two rules in the event table¶

The event table shows why (s,Q) falls behind on tea. Both branches below use the same ANN 95% target, 48 units, as S or as s. The dots are the inventory position at each review, before ordering; the bars are the orders placed at that review.

In [10]:

OUT_COLOR = "#4C78A8"
SQ_COLOR = "#F58518"

fig, axes = plt.subplots(
    2,
    1,
    figsize=(10, 5.5),
    sharex=True,
    layout="constrained",
)
for label, color, shift in [
    ("Order-up-to, ANN 95%", OUT_COLOR, -0.2),
    ("(s,Q), ANN 95%", SQ_COLOR, 0.2),
]:
    events = comparison[label].to_event_frame()
    tea = events.loc[events["unique_id"].eq("tea") & events["is_review_period"]]
    axes[0].plot(
        tea["date"],
        tea["decision_inventory_position"],
        marker="o",
        color=color,
        label=label,
    )
    axes[1].bar(
        tea["date"] + pd.Timedelta(days=shift),
        tea["order_quantity"],
        width=0.4,
        color=color,
        label=label,
    )
tea_target = target_table("ANN", 0.95).set_index("unique_id").loc["tea", "target"]
axes[0].axhline(
    tea_target,
    color="black",
    linestyle="--",
    linewidth=1,
    label=f"Target: S or s = {tea_target:.0f}",
)
axes[0].set(title="Tea: inventory position at each review", ylabel="Units")
axes[1].set(title="Tea: order placed at each review", ylabel="Units")
for ax in axes:
    ax.set_ylim(bottom=0)
    ax.grid(axis="y", alpha=0.2)
    ax.spines[["top", "right"]].set_visible(False)
    ax.legend(frameon=False, fontsize=8, loc="lower left")
axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
plt.show()

OUT_COLOR = "#4C78A8" SQ_COLOR = "#F58518" fig, axes = plt.subplots( 2, 1, figsize=(10, 5.5), sharex=True, layout="constrained", ) for label, color, shift in [ ("Order-up-to, ANN 95%", OUT_COLOR, -0.2), ("(s,Q), ANN 95%", SQ_COLOR, 0.2), ]: events = comparison[label].to_event_frame() tea = events.loc[events["unique_id"].eq("tea") & events["is_review_period"]] axes[0].plot( tea["date"], tea["decision_inventory_position"], marker="o", color=color, label=label, ) axes[1].bar( tea["date"] + pd.Timedelta(days=shift), tea["order_quantity"], width=0.4, color=color, label=label, ) tea_target = target_table("ANN", 0.95).set_index("unique_id").loc["tea", "target"] axes[0].axhline( tea_target, color="black", linestyle="--", linewidth=1, label=f"Target: S or s = {tea_target:.0f}", ) axes[0].set(title="Tea: inventory position at each review", ylabel="Units") axes[1].set(title="Tea: order placed at each review", ylabel="Units") for ax in axes: ax.set_ylim(bottom=0) ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) ax.legend(frameon=False, fontsize=8, loc="lower left") axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) plt.show()

What to notice:

• Order-up-to reviews at a position of about 30 (the target minus two days of demand) and orders the gap back to 48, so the order size follows demand, here between 12 and 25 units.
• (s,Q) finds the position below s = 48 at every review, so it orders 15 units every time. Fifteen units per two days is 7.5 a day, less than tea's average demand of 9. The position sinks until lost sales stop it, near 15 units from 21 March on.
• For the slow items the same 15 units are more than they sell between reviews (soup sells 4 a day), which is where the extra average stock of (s,Q) comes from. One fixed Q cannot fit items with different demand rates, while order-up-to adapts the order to each item.

All numbers come from one demand path; with run_comparison the same contrasts can be repeated on other paths or other targets.