Three ways to build a cumulative target¶

Notebook 04c took the target S straight from Smooth's cumulative upper bound. S is the stock level that should cover the total demand over the protection horizon, here lead time plus review period, with 95% probability. There is more than one way to get that number from the same forecast, and the ways differ in how they treat the relation between one day's demand and the next. This notebook builds three targets from one forecast and replays each of them on the same held-out demand.

• A, independent-normal: Stockcast adds up the daily means and the daily variances. This assumes that forecast errors on different days are independent.
• B, direct approximate: Smooth calculates the upper bound of the four-day total from its own model.
• C, simulated: Smooth simulates many future demand paths, totals each path, and takes the 95% quantile of the totals.

In [1]:

import pandas as pd

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

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

In [2]:

# Reuse the 04c tea scenario so this comparison starts from a known handoff.
sku = "tea_250g"
forecast_origin = pd.Timestamp("2026-02-14")
period_frequency = "D"
lead_time = 2
review_period = 2
service_level = 0.95
# An (R, S) target protects lead-time demand plus the next review interval.
protection_horizon = lead_time + review_period

# Reuse the 04c tea scenario so this comparison starts from a known handoff. sku = "tea_250g" forecast_origin = pd.Timestamp("2026-02-14") period_frequency = "D" lead_time = 2 review_period = 2 service_level = 0.95 # An (R, S) target protects lead-time demand plus the next review interval. protection_horizon = lead_time + review_period

1. The data¶

The same ten observed days and twelve held-out days as Notebook 04c. The history fits the one shared forecast, and the held-out days become the realized demand for the replay.

In [3]:

# Only these observations are available when the forecast is created.
history_values = [5.0, 6.0, 4.0, 7.0, 6.0, 5.0, 8.0, 6.0, 5.0, 7.0]
history = pd.DataFrame(
    {
        "unique_id": sku,
        "date": pd.date_range(
            end=forecast_origin,
            periods=len(history_values),
            freq=period_frequency,
        ),
        "y": history_values,
    }
)

# The holdout is kept separate and later becomes realized simulation demand.
holdout_values = [7.0, 5.0, 8.0, 6.0, 4.0, 7.0, 6.0, 9.0, 5.0, 6.0, 7.0, 5.0]
holdout = pd.DataFrame(
    {
        "unique_id": sku,
        "period": range(len(holdout_values)),
        "date": pd.date_range(
            forecast_origin + pd.Timedelta(days=1),
            periods=len(holdout_values),
            freq=period_frequency,
        ),
        "y": holdout_values,
    }
)

# The holdout begins the day after the origin and never overlaps history.
assert history["date"].max() == forecast_origin
assert holdout["date"].min() == forecast_origin + pd.Timedelta(days=1)

# Only these observations are available when the forecast is created. history_values = [5.0, 6.0, 4.0, 7.0, 6.0, 5.0, 8.0, 6.0, 5.0, 7.0] history = pd.DataFrame( { "unique_id": sku, "date": pd.date_range( end=forecast_origin, periods=len(history_values), freq=period_frequency, ), "y": history_values, } ) # The holdout is kept separate and later becomes realized simulation demand. holdout_values = [7.0, 5.0, 8.0, 6.0, 4.0, 7.0, 6.0, 9.0, 5.0, 6.0, 7.0, 5.0] holdout = pd.DataFrame( { "unique_id": sku, "period": range(len(holdout_values)), "date": pd.date_range( forecast_origin + pd.Timedelta(days=1), periods=len(holdout_values), freq=period_frequency, ), "y": holdout_values, } ) # The holdout begins the day after the origin and never overlaps history. assert history["date"].max() == forecast_origin assert holdout["date"].min() == forecast_origin + pd.Timedelta(days=1)

In [4]:

import matplotlib.pyplot as plt

HISTORY_COLOR = "#4C78A8"
HOLDOUT_COLOR = "#54A24B"


def style_axis(ax):
    ax.set_ylim(bottom=0)
    ax.grid(axis="y", alpha=0.2)
    ax.spines[["top", "right"]].set_visible(False)


# Plot the data boundary before fitting the forecasting model.
fig, ax = plt.subplots(figsize=(10, 3.5))
ax.plot(
    history["date"],
    history["y"],
    marker="o",
    linewidth=1.8,
    color=HISTORY_COLOR,
    label="History",
)
ax.plot(
    holdout["date"],
    holdout["y"],
    marker="o",
    linewidth=1.8,
    color=HOLDOUT_COLOR,
    label="Held-out demand",
)
ax.axvline(forecast_origin, color="0.35", linestyle="--", label="Forecast origin")
ax.set(title="History and held-out demand", xlabel="Date", ylabel="Demand (units)")
ax.legend(frameon=False, loc="lower center", ncol=3)
style_axis(ax)
plt.show()

import matplotlib.pyplot as plt HISTORY_COLOR = "#4C78A8" HOLDOUT_COLOR = "#54A24B" def style_axis(ax): ax.set_ylim(bottom=0) ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) # Plot the data boundary before fitting the forecasting model. fig, ax = plt.subplots(figsize=(10, 3.5)) ax.plot( history["date"], history["y"], marker="o", linewidth=1.8, color=HISTORY_COLOR, label="History", ) ax.plot( holdout["date"], holdout["y"], marker="o", linewidth=1.8, color=HOLDOUT_COLOR, label="Held-out demand", ) ax.axvline(forecast_origin, color="0.35", linestyle="--", label="Forecast origin") ax.set(title="History and held-out demand", xlabel="Date", ylabel="Demand (units)") ax.legend(frameon=False, loc="lower center", ncol=3) style_axis(ax) plt.show()

2. One forecast, three targets¶

Smooth fits one ES(model="ANN") model to the history. This optional example needs smooth>=1.0.7 on Python 3.11+. Smooth's ES class fixes the error distribution to Gaussian (dnorm), and that assumption needs checking for real demand. All three methods start from this one fit.

First the ordinary daily forecast: for each of the next four days a mean and a 95% upper bound. These daily bounds are the starting point of method A only. They are not the target. Adding four daily upper bounds would give a number that is too high, because it is unlikely that every day reaches its own 95% bound at once.

In [5]:

from smooth import ES

# One univariate fit on the observed history; every method reuses this object.
model = ES(model="ANN")
model.fit(history["y"])
print(f"Fitted smoothing parameter alpha: {float(model.coef[0]):.2f}")

daily_forecast = model.predict(
    h=protection_horizon,
    interval="prediction",
    level=service_level,
    side="upper",
    cumulative=False,
)

forecast_dates = pd.date_range(
    forecast_origin + pd.Timedelta(days=1),
    periods=protection_horizon,
    freq=period_frequency,
)
daily_forecast_output = pd.DataFrame(
    {
        "date": forecast_dates,
        "forecast_mean": daily_forecast.mean.to_numpy(),
        "daily_upper_95": daily_forecast.upper.iloc[:, 0].to_numpy(),
        "heldout_demand": holdout.loc[: protection_horizon - 1, "y"].to_numpy(),
    }
)

# Four daily means with separate upper bounds, beside demand later observed.
daily_forecast_output

from smooth import ES # One univariate fit on the observed history; every method reuses this object. model = ES(model="ANN") model.fit(history["y"]) print(f"Fitted smoothing parameter alpha: {float(model.coef[0]):.2f}") daily_forecast = model.predict( h=protection_horizon, interval="prediction", level=service_level, side="upper", cumulative=False, ) forecast_dates = pd.date_range( forecast_origin + pd.Timedelta(days=1), periods=protection_horizon, freq=period_frequency, ) daily_forecast_output = pd.DataFrame( { "date": forecast_dates, "forecast_mean": daily_forecast.mean.to_numpy(), "daily_upper_95": daily_forecast.upper.iloc[:, 0].to_numpy(), "heldout_demand": holdout.loc[: protection_horizon - 1, "y"].to_numpy(), } ) # Four daily means with separate upper bounds, beside demand later observed. daily_forecast_output

Fitted smoothing parameter alpha: 0.00

Out[5]:

	date	forecast_mean	daily_upper_95	heldout_demand
0	2026-02-15	5.90	7.99	7.00
1	2026-02-16	5.90	7.99	5.00
2	2026-02-17	5.90	7.99	8.00
3	2026-02-18	5.90	7.99	6.00

Method A: independent-normal, inside Stockcast¶

For a normal forecast, each daily upper bound reveals the daily standard deviation: sigma = (q95 - mean) / z, where z is the 95% point of the standard normal curve. Stockcast then calculates

S = sum(mean) + z * sqrt(sum(sigma^2))

This is correct only when the daily forecast errors are independent. Normal marginals alone do not establish that, so Stockcast records the result as independent_normal_marginals and the reader decides whether the assumption is defensible. The inputs go to OrderUpToPolicy from stockcast.policies as a table with one row per forecast day.

In [6]:

from statistics import NormalDist

# z converts each upper bound back into the sigma it implies under normality.
z = NormalDist().inv_cdf(service_level)
marginal_inputs = pd.DataFrame(
    {
        "unique_id": sku,
        "fh": range(1, protection_horizon + 1),
        "date": forecast_dates,
        "mean": daily_forecast_output["forecast_mean"].to_numpy(),
        "std": (
            (
                daily_forecast_output["daily_upper_95"]
                - daily_forecast_output["forecast_mean"]
            )
            / z
        ).to_numpy(),
    }
)
marginal_inputs

from statistics import NormalDist # z converts each upper bound back into the sigma it implies under normality. z = NormalDist().inv_cdf(service_level) marginal_inputs = pd.DataFrame( { "unique_id": sku, "fh": range(1, protection_horizon + 1), "date": forecast_dates, "mean": daily_forecast_output["forecast_mean"].to_numpy(), "std": ( ( daily_forecast_output["daily_upper_95"] - daily_forecast_output["forecast_mean"] ) / z ).to_numpy(), } ) marginal_inputs

Out[6]:

	unique_id	fh	date	mean	std
0	tea_250g	1	2026-02-15	5.90	1.27
1	tea_250g	2	2026-02-16	5.90	1.27
2	tea_250g	3	2026-02-17	5.90	1.27
3	tea_250g	4	2026-02-18	5.90	1.27

In [7]:

from stockcast.policies import OrderUpToPolicy

policy_independent = OrderUpToPolicy(
    lead_time=lead_time,
    review_period=review_period,
    freq=period_frequency,
    service_level=service_level,
    allow_backorders=False,
    date_column="date",
)
policy_independent.fit(
    marginal_inputs,
    mean_column="mean",
    std_column="std",
    forecast_origin=forecast_origin,
)

# The fitted target, with how it was calculated.
independent_target = policy_independent.get_target_levels().copy()
independent_target["calculation_method"] = policy_independent.get_target_metadata()[
    "calculation_method"
]
independent_target

from stockcast.policies import OrderUpToPolicy policy_independent = OrderUpToPolicy( lead_time=lead_time, review_period=review_period, freq=period_frequency, service_level=service_level, allow_backorders=False, date_column="date", ) policy_independent.fit( marginal_inputs, mean_column="mean", std_column="std", forecast_origin=forecast_origin, ) # The fitted target, with how it was calculated. independent_target = policy_independent.get_target_levels().copy() independent_target["calculation_method"] = policy_independent.get_target_metadata()[ "calculation_method" ] independent_target

Out[7]:

	unique_id	target_level	calculation_method
0	tea_250g	27.78	independent_normal

Methods B and C: Smooth calculates the total¶

With cumulative=True, Smooth returns the distribution of the four-day total directly, so Stockcast needs no independence assumption. Its approximate route estimates an upper bound for the total from the model's own variance. Its simulated route generates many joint future paths, adds up each path, and takes the quantile of those totals.

The simulated result has Monte Carlo variability. Smooth (checked with 1.0.7 and 1.0.8) ignores its seed argument for interval="simulated", so repeated runs can differ. Keep the generated targets and the forecasting settings if you need to replay a run; the engine's own seed does not seed the external forecast.

In [8]:

# Two calls, same question: the 95% upper bound of the four-day total.
cumulative_approx = model.predict(
    h=protection_horizon,
    interval="approximate",
    level=service_level,
    side="upper",
    cumulative=True,
)
cumulative_sim = model.predict(
    h=protection_horizon,
    interval="simulated",
    level=service_level,
    side="upper",
    cumulative=True,
    nsim=5000,
)
approx_target = float(cumulative_approx.upper.iloc[-1, 0])
sim_target = float(cumulative_sim.upper.iloc[-1, 0])

# Two calls, same question: the 95% upper bound of the four-day total. cumulative_approx = model.predict( h=protection_horizon, interval="approximate", level=service_level, side="upper", cumulative=True, ) cumulative_sim = model.predict( h=protection_horizon, interval="simulated", level=service_level, side="upper", cumulative=True, nsim=5000, ) approx_target = float(cumulative_approx.upper.iloc[-1, 0]) sim_target = float(cumulative_sim.upper.iloc[-1, 0])

3. Compare the three targets¶

All three share the horizon and the requested probability. Any difference comes from the independence assumption, the way the total is calculated, or Monte Carlo error. covers_realized_total checks each target against the demand that was later observed over the same four days. One realization cannot establish 95% coverage, so read it as an illustration.

In [9]:

import numpy as np

LABELS = ["A: independent-normal", "B: direct approximate", "C: simulated"]
independent_level = float(independent_target.loc[0, "target_level"])
realized_total = float(holdout.loc[: protection_horizon - 1, "y"].sum())

comparison_targets = pd.DataFrame(
    {
        "method": LABELS,
        "protection_target": [independent_level, approx_target, sim_target],
        "construction": [
            "sum(mean) + z * sqrt(sum(sigma^2))",
            "model variance of the total",
            "quantile of summed joint paths (5000)",
        ],
    }
)

# External targets must be finite and non-negative before they reach the policy.
assert np.isfinite(comparison_targets["protection_target"]).all()
assert comparison_targets["protection_target"].ge(0).all()
comparison_targets["covers_realized_total"] = comparison_targets[
    "protection_target"
].ge(realized_total)

comparison_targets

import numpy as np LABELS = ["A: independent-normal", "B: direct approximate", "C: simulated"] independent_level = float(independent_target.loc[0, "target_level"]) realized_total = float(holdout.loc[: protection_horizon - 1, "y"].sum()) comparison_targets = pd.DataFrame( { "method": LABELS, "protection_target": [independent_level, approx_target, sim_target], "construction": [ "sum(mean) + z * sqrt(sum(sigma^2))", "model variance of the total", "quantile of summed joint paths (5000)", ], } ) # External targets must be finite and non-negative before they reach the policy. assert np.isfinite(comparison_targets["protection_target"]).all() assert comparison_targets["protection_target"].ge(0).all() comparison_targets["covers_realized_total"] = comparison_targets[ "protection_target" ].ge(realized_total) comparison_targets

Out[9]:

	method	protection_target	construction	covers_realized_total
0	A: independent-normal	27.78	sum(mean) + z * sqrt(sum(sigma^2))	True
1	B: direct approximate	27.78	model variance of the total	True
2	C: simulated	28.25	quantile of summed joint paths (5000)	True

A and B agree here. The fitted smoothing parameter printed in section 2 is zero, so the model forecasts a constant level, the history mean of 5.9, with the same spread on every day. The daily errors then really are independent, which is exactly what method A assumes. With a model that smooths with a larger parameter, their errors on neighbouring days would be related, and A and B would differ. Method C comes from random paths, so it varies a little from run to run and can sit slightly above or below B. Notice also that none of the three targets is anywhere near the sum of the four daily upper bounds, 31.96, which is why that sum is never used.

4. Replay the targets on the same demand¶

Each target becomes a fitted (R, S) policy with the same lead time, review period, and probability. SimulationEngine.run_comparison from stockcast.core then replays all three policies on identical held-out demand and the same opening stock, so any difference in the results comes from the target alone. The targets stay fixed here; Notebook 04d shows them updating.

In [10]:

def direct_policy(target_value):
    """Fit an (R, S) policy from one externally calculated cumulative target."""
    policy = OrderUpToPolicy(
        lead_time=lead_time,
        review_period=review_period,
        freq=period_frequency,
        service_level=service_level,
        allow_backorders=False,
        date_column="target_end_date",
    )
    policy.fit(
        pd.DataFrame(
            {
                "unique_id": [sku],
                "protection_target": [target_value],
                "target_end_date": [forecast_dates[-1]],
            }
        ),
        target_column="protection_target",
        forecast_origin=forecast_origin,
    )
    return policy


policy_approx = direct_policy(approx_target)
policy_simulated = direct_policy(sim_target)

# Each fitted policy holds exactly the target listed above.
policies = [policy_independent, policy_approx, policy_simulated]
fitted_levels = [float(p.get_target_levels().loc[0, "target_level"]) for p in policies]
assert np.allclose(fitted_levels, comparison_targets["protection_target"])

def direct_policy(target_value): """Fit an (R, S) policy from one externally calculated cumulative target.""" policy = OrderUpToPolicy( lead_time=lead_time, review_period=review_period, freq=period_frequency, service_level=service_level, allow_backorders=False, date_column="target_end_date", ) policy.fit( pd.DataFrame( { "unique_id": [sku], "protection_target": [target_value], "target_end_date": [forecast_dates[-1]], } ), target_column="protection_target", forecast_origin=forecast_origin, ) return policy policy_approx = direct_policy(approx_target) policy_simulated = direct_policy(sim_target) # Each fitted policy holds exactly the target listed above. policies = [policy_independent, policy_approx, policy_simulated] fitted_levels = [float(p.get_target_levels().loc[0, "target_level"]) for p in policies] assert np.allclose(fitted_levels, comparison_targets["protection_target"])

In [11]:

# Opening stock is a declared inventory input, independent of the forecast.
opening_stock = pd.DataFrame({"unique_id": [sku], "observed_on_hand": [12.0]})
opening_stock

# Opening stock is a declared inventory input, independent of the forecast. opening_stock = pd.DataFrame({"unique_id": [sku], "observed_on_hand": [12.0]}) opening_stock

Out[11]:

	unique_id	observed_on_hand
0	tea_250g	12.00

In [12]:

from stockcast.core import InventoryStateDataFrame

# The engine deep-copies this state per branch; the caller-owned object is preserved.
empty_state = InventoryStateDataFrame(
    [sku],
    max_lead_time=lead_time,
    allow_backorders=False,
)
opening_inventory = empty_state.initialize_from_observed(
    opening_stock,
    on_hand_column="observed_on_hand",
    opening_date=forecast_origin,
)
opening_inventory

from stockcast.core import InventoryStateDataFrame # The engine deep-copies this state per branch; the caller-owned object is preserved. empty_state = InventoryStateDataFrame( [sku], max_lead_time=lead_time, allow_backorders=False, ) opening_inventory = empty_state.initialize_from_observed( opening_stock, on_hand_column="observed_on_hand", opening_date=forecast_origin, ) opening_inventory

Out[12]:

InventoryStateDataFrame(n_skus=1, total_on_hand=12, total_backorders=0, has_stockout=False, has_backorder=False)

In [13]:

from stockcast.core import SimulationEngine

# Three branches, one shared demand calendar and timing contract.
comparison = SimulationEngine().run_comparison(
    policies=policies,
    labels=LABELS,
    demand_source=holdout,
    inventory=opening_inventory,
)

from stockcast.core import SimulationEngine # Three branches, one shared demand calendar and timing contract. comparison = SimulationEngine().run_comparison( policies=policies, labels=LABELS, demand_source=holdout, inventory=opening_inventory, )

5. Read the inventory outcomes¶

stockcast.evaluation applies one contract to every branch: the same scoring window and the same pooled metrics. These are unit and service measures on declared synthetic demand, so no cost or operational claim follows.

In [14]:

from stockcast.evaluation import (
    InventoryEvaluator,
    avg_on_hand,
    demand_period_service_level,
    fill_rate,
    order_event_count,
    order_units,
)

metric_rows = []
for label, branch in comparison.results.items():
    evaluator = InventoryEvaluator()
    evaluator.fit(branch, window="scoring")
    row = evaluator.evaluate(
        metrics=[
            fill_rate,
            demand_period_service_level,
            avg_on_hand,
            order_units,
            order_event_count,
        ],
        groupby=[],
    )
    metric_rows.append(row.assign(method=label))
pd.concat(metric_rows).set_index("method")

from stockcast.evaluation import ( InventoryEvaluator, avg_on_hand, demand_period_service_level, fill_rate, order_event_count, order_units, ) metric_rows = [] for label, branch in comparison.results.items(): evaluator = InventoryEvaluator() evaluator.fit(branch, window="scoring") row = evaluator.evaluate( metrics=[ fill_rate, demand_period_service_level, avg_on_hand, order_units, order_event_count, ], groupby=[], ) metric_rows.append(row.assign(method=label)) pd.concat(metric_rows).set_index("method")

Out[14]:

	fill_rate	demand_period_service_level	avg_on_hand	order_units	order_event_count
method
A: independent-normal	1.00	1.00	5.31	78.78	6
B: direct approximate	1.00	1.00	5.31	78.78	6
C: simulated	1.00	1.00	5.71	79.25	6

See the targets and the stock they produce¶

The left panel compares the three targets with the demand that was later observed over the same four days, and with the sum of the daily upper bounds, drawn hatched because it is not a valid target. The right panel shows end-of-day stock for each branch. A and B overlap exactly, so B is dashed on top of A. In this run C's target is a little higher, so it holds a little more stock. Every branch served all demand, so the extra stock buys no better service here.

One realization cannot show calibration: a 95% demand target does not guarantee a 95% fill rate, and the suitability of the uncertainty model needs separate held-out evidence. For a bootstrap workflow, generate coherent future-demand paths with the dependence you want to keep, total each path over the protection horizon, take the 95% quantile, and hand the result over as in method B.

In [15]:

import matplotlib.dates as mdates

A_COLOR = "#F58518"
B_COLOR = "#4C78A8"
C_COLOR = "#54A24B"
branch_colors = dict(zip(LABELS, [A_COLOR, B_COLOR, C_COLOR]))
daily_upper_sum = float(daily_forecast_output["daily_upper_95"].sum())

fig, axes = plt.subplots(1, 2, figsize=(13, 4.6), layout="constrained")

# Left: the targets, the demand later observed, and the sum of daily bounds.
bars = axes[0].bar(
    comparison_targets["method"],
    comparison_targets["protection_target"],
    color=list(branch_colors.values()),
)
axes[0].bar(
    ["Sum of daily\nupper bounds"],
    [daily_upper_sum],
    color="none",
    edgecolor="0.5",
    hatch="//",
)
axes[0].axhline(
    realized_total,
    color="0.35",
    linestyle="--",
    linewidth=1.4,
    label=f"Realized four-day total ({realized_total:.1f})",
)
for bar in [*bars, *axes[0].patches[3:]]:
    axes[0].text(
        bar.get_x() + bar.get_width() / 2,
        bar.get_height() + 0.3,
        f"{bar.get_height():.1f}",
        ha="center",
        fontsize=9,
    )
axes[0].set(
    title="Same forecast, different cumulative targets", ylabel="Target S (units)"
)
axes[0].tick_params(axis="x", rotation=12)
axes[0].legend(frameon=False, loc="lower right")

# Right: end-of-day stock; receipts arrive first, then demand is served.
for label, branch in comparison.results.items():
    branch_events = branch.to_event_frame(window="all")
    branch_periods = branch_events.loc[branch_events["event_type"].eq("period")]
    axes[1].step(
        branch_periods["date"],
        branch_periods["ending_on_hand"],
        where="mid",
        color=branch_colors[label],
        linewidth=2.4 if label.startswith("A") else 1.8,
        linestyle="--" if label.startswith("B") else "-",
        label=label,
    )
axes[1].set(
    title="End-of-day stock under the three targets",
    xlabel="Date",
    ylabel="On hand (units)",
)
axes[1].set_xticks(list(branch_periods["date"])[::2])
axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))
axes[1].tick_params(axis="x", rotation=30)
axes[1].legend(frameon=False, loc="upper left")
for ax in axes:
    style_axis(ax)
plt.show()

import matplotlib.dates as mdates A_COLOR = "#F58518" B_COLOR = "#4C78A8" C_COLOR = "#54A24B" branch_colors = dict(zip(LABELS, [A_COLOR, B_COLOR, C_COLOR])) daily_upper_sum = float(daily_forecast_output["daily_upper_95"].sum()) fig, axes = plt.subplots(1, 2, figsize=(13, 4.6), layout="constrained") # Left: the targets, the demand later observed, and the sum of daily bounds. bars = axes[0].bar( comparison_targets["method"], comparison_targets["protection_target"], color=list(branch_colors.values()), ) axes[0].bar( ["Sum of daily\nupper bounds"], [daily_upper_sum], color="none", edgecolor="0.5", hatch="//", ) axes[0].axhline( realized_total, color="0.35", linestyle="--", linewidth=1.4, label=f"Realized four-day total ({realized_total:.1f})", ) for bar in [*bars, *axes[0].patches[3:]]: axes[0].text( bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.3, f"{bar.get_height():.1f}", ha="center", fontsize=9, ) axes[0].set( title="Same forecast, different cumulative targets", ylabel="Target S (units)" ) axes[0].tick_params(axis="x", rotation=12) axes[0].legend(frameon=False, loc="lower right") # Right: end-of-day stock; receipts arrive first, then demand is served. for label, branch in comparison.results.items(): branch_events = branch.to_event_frame(window="all") branch_periods = branch_events.loc[branch_events["event_type"].eq("period")] axes[1].step( branch_periods["date"], branch_periods["ending_on_hand"], where="mid", color=branch_colors[label], linewidth=2.4 if label.startswith("A") else 1.8, linestyle="--" if label.startswith("B") else "-", label=label, ) axes[1].set( title="End-of-day stock under the three targets", xlabel="Date", ylabel="On hand (units)", ) axes[1].set_xticks(list(branch_periods["date"])[::2]) axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %d")) axes[1].tick_params(axis="x", rotation=30) axes[1].legend(frameon=False, loc="upper left") for ax in axes: style_axis(ax) plt.show()