From Forecasts to Decisions: Inventory Planning with Probabilistic Programming

Mean Forecast Plan

We replace the distribution by the mean path of the scenarios, \(\bar{D} = \frac{1}{N} \sum_{n} D^{(n)}\), and solve the deterministic problem

\[\min_{a \geq 0} \; L(a, \bar{D}).\]

This is the plan that a point forecast gives. Against a single path, the cheapest plan brings the stock available on each day to the demand of that day and leaves nothing on the shelf, so the orders follow the mean path once the initial stock is used up.

Chance-Constrained Plan

The chance-constrained plan does not use the costs. It fixes a service level instead, which is how often we are willing to run out. A stockout on day \(t\) is the event \(D_t > y_t\), and the fraction of scenarios in which it happens is our estimate of the stockout probability \(\text{P}(D_t > y_t)\). We choose a small tail probability \(\alpha\), here \(\alpha = 0.05\), and require

\[\text{P}(D_t > y_t) \leq \alpha \quad \text{for every } t = 1, \dots, H,\]

and we look for the plan with the smallest cumulative orders that satisfies it.

Note that the constraint holds day by day. It bounds the chance of a lost sale on any single day, not the chance of at least one lost sale somewhere in the \(15\) days. That probability is larger, because the plan has \(15\) chances to fail: with independent days it would be \(1 - (1 - \alpha)^{15} \approx 54\%\) for \(\alpha = 0.05\).

The days are not independent, so we report both numbers in the evaluation: the stockout probability on a given day, which the constraint controls, and the probability of at least one stockout in the horizon, which it does not.

We build the plan with a forward pass over the days of the horizon, inside all the scenarios at once. Given the orders of days \(1\) to \(t - 1\), every scenario \(n\) carries a simulated stock \(I_{t-1}^{(n)}\) into day \(t\), and the shortfall that the order has to cover in that scenario is \(D_t^{(n)} - I_{t-1}^{(n)}\). The smallest order that leaves at most a fraction \(\alpha\) of the scenarios short is the \((1 - \alpha)\) quantile of that shortfall across the scenarios,

\[a_t = \left( Q_{1-\alpha}\!\left(D_t - I_{t-1}\right) \right)^{+}.\]

Here is an important observation: the forward pass does not place an order each day after observing the demand. It walks through the horizon inside the scenarios, which we have today, and it never sees the real demand. Its output is the full vector \(a\), and we place all \(15\) orders today, as the problem statement requires.

One forward pass over the days builds the whole plan, and no optimizer is needed. The constraint counts how often we run out and ignores by how much.

Remark: This plan leaves the costs out, but the next two plans depend on them, so they need a good estimate of the cost of a lost sale. In many cases, not serving the demand is more costly than having leftovers. For example, there is a cost of customers having a bad experience, and therefore churning to competitors and never coming back. In real applications, it is important to account for this asymmetry in the cost function with estimates coming from customer surveys and other sources (or causal models).

Expected Cost Plan

We minimize the sample average of the loss,

\[\min_{a \geq 0} \; \frac{1}{N} \sum_{n=1}^{N} L(a, D^{(n)}),\]

which is the sample average approximation of the Bayes decision, see Shapiro, Dentcheva and Ruszczyński (2021) for the general theory. To see what this plan does, let’s look at one day in isolation. The stock available is \(y_t\), the demand is \(D_t\), a unit left over costs \(h\) and a unit short costs \(p\). This single-period problem has a name, the newsvendor problem, after the vendor who buys newspapers in the morning without knowing how many will sell. Its solution goes back to Arrow, Harris and Marschak (1951).

Its solution has a simple form. Suppose we raise the stock by one unit. With probability \(\text{P}(D_t \leq y_t)\) the unit is left over and costs \(h\), and with the remaining probability it covers a unit of demand and saves \(p\). The expected cost stops decreasing when the two balance:

\[h \, \text{P}(D_t \leq y_t) = p \, \left[1 - \text{P}(D_t \leq y_t)\right] \quad \Longleftrightarrow \quad \text{P}(D_t \leq y_t) = \frac{p}{p + h}.\]

The optimal stock is the \(p / (p + h)\) quantile of the demand distribution, and the ratio \(p / (p + h)\) is called the critical fractile.

Note that with \(p = h\) the critical fractile is \(1 / 2\), so the optimal stock is the median of the demand. For a symmetric demand distribution, such as the Normal of the example below, the median is the mean. So with equal costs and a symmetric demand the mean forecast is the optimal single-period stock. With our costs, \(p = 2\) and \(h = 0.1\), the fractile is \(0.952\), so the optimal stock is a high quantile of the demand, far above the mean.

Let’s work a simple example with a Normal demand with mean \(\mu = 60\) and standard deviation \(\sigma = 15\). Let \(\Phi\) and \(\varphi\) be the distribution function and the density of the standard Normal. We write the demand as \(D = \mu + \sigma Z\) with \(Z\) standard Normal.

The stock \(q\) is a number that we choose, not a random variable, and we let \(z = (q - \mu) / \sigma\) be the stock in standard units. The gap between the stock and the demand is then \(q - D = (q - \mu) - \sigma Z = \sigma \, (z - Z)\), so the leftover is \((q - D)^{+} = \sigma \, (z - Z)^{+}\) and the expected leftover is \(\sigma \, \text{E}[(z - Z)^{+}]\), which is the integral of \((z - u) \, \varphi(u)\) over \(u \leq z\):

\[\text{E}[(q - D)^{+}] = \sigma \int_{-\infty}^{z} (z - u) \, \varphi(u) \, du = \sigma \left[ z \, \Phi(z) + \varphi(z) \right],\]

where we used \(\varphi'(u) = -u \, \varphi(u)\) for the second term.

The expected shortage follows from the identity \((D - q)^{+} = (q - D)^{+} + (D - q)\). Its expectation is the expected leftover minus \(q - \mu = \sigma z\):

\[\text{E}[(D - q)^{+}] = \sigma \left[ \varphi(z) - z \, \left(1 - \Phi(z)\right) \right].\]

We compute both with jax.scipy.stats.norm, add the two parts of the expected cost and plot it as a function of the stock, with the critical fractile of our costs.

critical_fractile = config.lost_sales_cost / (
    config.lost_sales_cost + config.holding_cost
)

demand_mean, demand_sd = 60.0, 15.0
stock_grid = jnp.linspace(20, 120, 201)
z = (stock_grid - demand_mean) / demand_sd
expected_leftover = demand_sd * (z * norm.cdf(z) + norm.pdf(z))
expected_shortage = demand_sd * (norm.pdf(z) - z * (1 - norm.cdf(z)))
holding_part = config.holding_cost * expected_leftover
lost_sales_part = config.lost_sales_cost * expected_shortage
optimal_stock = demand_mean + demand_sd * norm.ppf(critical_fractile)

fig, ax = plt.subplots()
ax.plot(
    stock_grid,
    holding_part + lost_sales_part,
    color="C0",
    linewidth=2,
    label="expected cost",
)
ax.plot(stock_grid, holding_part, color="C2", linestyle=":", label="holding part")
ax.plot(stock_grid, lost_sales_part, color="C3", linestyle=":", label="lost sales part")
ax.axvline(demand_mean, color="black", linewidth=1, label="mean demand")
ax.axvline(
    optimal_stock,
    color="C1",
    linewidth=2,
    label=f"{critical_fractile:.3f} quantile (critical fractile)",
)
fig.legend(loc="outside lower center", ncol=3)
ax.set(xlabel="stock $q$", ylabel="expected cost of one period")
fig.suptitle("The Newsvendor Trade-Off", fontsize=18, fontweight="bold");

The expected cost has its minimum at the \(0.952\) quantile, about \(85\) units, and not at the mean of \(60\). The lost sales part makes the curve steep on the left and the holding part makes it flat on the right, so a stock above the mean is cheaper than the mean itself.

Our problem has \(15\) days, and the days are linked by the stock that carries over. The single-period rule charges a leftover unit one day of holding, while in the multi-period problem the unit pays \(h\) on every later day until it is sold. It is also left over in the scenarios where day \(t\) was slow, which carry more stock into day \(t + 1\) and tend to have a low level, so the unit can take several days to sell.

So the marginal unit costs more than the single-period rule assumes. The plan that minimizes the expected loss keeps the probability \(\text{P}(D_t \leq y_t)\) of covering the demand of day \(t\) below \(p / (p + h)\) on every day but the last, so the single-period rule, applied day by day, overstocks. We check this in the Out of Sample Evaluation section, where the expected cost plan runs out on \(10.4\%\) of the days against the single-period \(4.8\%\).

The general lost-sales problem, with delivery lead times and orders that react to the stock, is hard, see Zipkin (2008) and the review of Bijvank and Vis (2011). Our plan is fixed in advance, and this makes the problem a convex program: the loss is a convex function of the orders \(a\), so every local minimum is the global one and a gradient method cannot get trapped in a local minimum.

The loss is convex because it can be written as a positive combination of the daily stocks plus a term that is linear in the orders, and each stock is a convex piecewise linear function of the orders. We show this in the Appendix, and we solve the problem numerically.

We keep the single-period rule as a reference, the myopic newsvendor plan: on each day, order up to the \(p / (p + h)\) quantile of the shortfall \(D_t - I_{t-1}\), with the same forward pass as the chance-constrained plan.

The expected cost plan has no service target of its own, i.e. we do not set a stockout probability. The stockout frequency it accepts is a consequence of the costs, and we read it from the outputs below. On the last day of the horizon the problem is a single-period one, since a leftover there is charged holding only once, so the optimal stockout probability on that day is \(h / (p + h)\).

The myopic rule uses the same critical fractile as the order rule of the winning VN2 entry, where Szabłowski (2026) maps \(c_s / (c_s + c_h)\) to a Normal quantile around a point forecast and re-plans every week on the observed inventory. We read the quantile from the scenarios instead, which the report itself recommends:

In production, \(q^\star\) would be implemented via empirical quantiles or a probabilistic forecast, avoiding any normality assumption.

CVaR Plan

The CVaR plan adds a penalty on the worst outcomes. Let \(\text{VaR}_\alpha(L) = Q_{1-\alpha}(L)\) be the \((1 - \alpha)\) quantile of the loss, the value at risk. For a continuous loss, the conditional value at risk is the expected loss beyond that quantile:

\[\text{CVaR}_\alpha(L) = \text{E}\left[ L \mid L \geq \text{VaR}_\alpha(L) \right].\]

With \(N\) equally weighted scenarios, the sample estimate is the mean of the \(k = \alpha N\) largest losses, and this is what the function cvar computes below.

Again, let’s work a simple example. We draw the two quantities on a skewed loss density, a LogNormal with log-mean \(\mu = 5\) and log-standard deviation \(\sigma = 0.6\), for the tail probability \(\alpha = 0.05\). If \(L\) is LogNormal, then \(W = \log L\) is \(\text{Normal}(\mu, \sigma)\), and every quantity we need has a closed form. Let’s derive them one at a time.

First, let’s derive the density. Let \(F\) be the distribution function of \(L\) and \(f = F'\) its density. For \(x > 0\) the event \(L \leq x\) is the event \(W \leq \log x\), so

\[F(x) = \text{P}(L \leq x) = \text{P}(W \leq \log x) = \Phi\!\left(\frac{\log x - \mu}{\sigma}\right).\]

The density is the derivative of \(F\), and the chain rule gives

\[f(x) = \frac{1}{\sigma x} \, \varphi\!\left(\frac{\log x - \mu}{\sigma}\right).\]

Next, we turn to the value at risk. Let \(z_{1-\alpha} = \Phi^{-1}(1 - \alpha)\) be the \((1 - \alpha)\) quantile of the standard Normal. The value at risk is the \(x\) with \(F(x) = 1 - \alpha\), so \((\log x - \mu) / \sigma = z_{1-\alpha}\) and

\[\text{VaR}_\alpha(L) = \exp\left(\mu + \sigma \, z_{1-\alpha}\right).\]

Now we compute the mean. We write \(W = \mu + \sigma Z\) with \(Z\) standard Normal, so

\[\text{E}[L] = \text{E}\left[e^{W}\right] = \int_{-\infty}^{\infty} e^{\mu + \sigma z} \, \varphi(z) \, dz.\]

Completing the square in the exponent, \(\sigma z - z^2 / 2 = \sigma^2 / 2 - (z - \sigma)^2 / 2\), so the integrand is \(e^{\mu + \sigma^2 / 2} \, \varphi(z - \sigma)\), and the density \(\varphi(z - \sigma)\) integrates to one:

\[\text{E}[L] = \exp\left(\mu + \sigma^2 / 2\right).\]

Finally, we compute the conditional value at risk. It is the same integral restricted to the tail \(Z \geq z_{1-\alpha}\), divided by the probability \(\alpha\) of that tail:

\[\text{CVaR}_\alpha(L) = \frac{1}{\alpha} \int_{z_{1-\alpha}}^{\infty} e^{\mu + \sigma z} \, \varphi(z) \, dz = \frac{e^{\mu + \sigma^2 / 2}}{\alpha} \int_{z_{1-\alpha}}^{\infty} \varphi(z - \sigma) \, dz.\]

The last integral is \(1 - \Phi(z_{1-\alpha} - \sigma) = \Phi(\sigma - z_{1-\alpha})\), so

\[\text{CVaR}_\alpha(L) = \frac{1}{\alpha} \exp\left(\mu + \sigma^2 / 2\right) \Phi\left(\sigma - z_{1-\alpha}\right).\]

In the code, dist.LogNormal from NumPyro gives the density, the quantile and the mean, and we compute the conditional value at risk with the last formula. We fix \(\alpha\) here, since the chance constraint and the CVaR share it.

ALPHA = 0.05  # tail probability of the chance constraint and of the CVaR

log_mean, log_sd = 5.0, 0.6
lognormal = dist.LogNormal(log_mean, log_sd)
loss_grid = jnp.linspace(1, 1_000, 1_000)
loss_density = jnp.exp(lognormal.log_prob(loss_grid))
z_alpha = norm.ppf(1 - ALPHA)
expected_loss_example = lognormal.mean
var_example = lognormal.icdf(1 - ALPHA)
cvar_example = expected_loss_example * norm.cdf(log_sd - z_alpha) / ALPHA

fig, ax = plt.subplots()
ax.plot(loss_grid, loss_density, color="C0", linewidth=2, label="loss density")
ax.fill_between(
    loss_grid,
    loss_density,
    where=loss_grid >= var_example,
    color="C3",
    alpha=0.3,
    label=f"worst {ALPHA:.0%} of outcomes",
)
ax.axvline(expected_loss_example, color="black", linewidth=1, label="expected loss")
ax.axvline(var_example, color="C1", linewidth=2, label="VaR")
ax.axvline(cvar_example, color="C3", linewidth=2, label="CVaR")
fig.legend(loc="outside lower center", ncol=3)
ax.set(xlabel="loss", ylabel="density", yticks=[])
fig.suptitle(
    "Value at Risk and Conditional Value at Risk", fontsize=18, fontweight="bold"
);

The value at risk is the point that the loss exceeds with probability \(5\%\). The conditional value at risk is the average loss in the shaded region beyond it, so it sits further to the right, and well above the expected loss because the density has a long right tail.

Having understood the intuition of the CVaR, we can write the objective that the CVaR plan minimizes in practice:

\[\min_{a \geq 0} \; \text{E}[L(a, D)] + \lambda \, \text{CVaR}_\alpha(L(a, D)),\]

where \(\lambda \geq 0\) is the tail weight, the weight we put on the worst \(\alpha\) fraction of the outcomes. With \(\lambda = 0\) we recover the expected cost plan. The CVaR term counts how bad the worst outcomes are, not how often they happen.

Now that we have these concepts in place, let’s put them into practice on a synthetic data set.

Prior Elicitation

The model has five global parameters. For three of them we state what we believe on the scale of the data and use PreliZ to turn the statement into a distribution: pz.maxent returns the member of a family that puts a given mass between two bounds (the one with maximum entropy, so it adds nothing beyond our constraints). Here is the reasoning for each of the three.

• intercept: \(\exp(\beta_0)\) is the demand per day at the start of the history at the list price, with the weekly pattern averaged out on the log scale (the Fourier terms sum to zero over a week). We believe it lies between \(10\) and \(300\) units, a wide range for a product that sells tens of units per day. So we ask for a Normal on \(\beta_0\) with \(94\%\) of its mass between \(\log 10\) and \(\log 300\).

• drift_scale: the level is a random walk, so after \(T\) days the standard deviation of the log level is \(\sigma_\delta \sqrt{T}\), and \(\exp(\sigma_\delta \sqrt{T})\) is the typical factor by which the level has moved. Over the \(364\) days of history we believe that factor lies between \(1.04\) and \(2\), that is, from a \(4\%\) move to a doubling. So we ask for a LogNormal on \(\sigma_\delta\) with \(94\%\) of its mass between \(\log(1.04) / \sqrt{T}\) and \(\log(2) / \sqrt{T}\). The LogNormal density vanishes at zero, which fits a level that moves at least a little over a year.

• concentration: the negative binomial adds \(\mu_t^2 / \phi\) to the Poisson variance, so the extra day-to-day noise is a fraction \(1 / \sqrt{\phi}\) of the level. We believe that fraction lies between a tenth and a half, which is \(\phi\) between \(4\) and \(100\). So we ask for a Gamma with \(94\%\) of its mass in that range.

The other two priors we state directly.

• Each Fourier coefficient in beta gets a \(\text{Normal}(0, 0.3)\), so two standard deviations multiply the demand by \(e^{\pm 0.6}\), between \(0.55\) and \(1.8\).

• The elasticity gets a \(\text{Normal}(-1, 1)\): a higher price should lower the demand, so we center the prior on a negative value, and the standard deviation of \(1\) leaves a positive value possible if the data ask for it.

We run the three elicitations and plot the densities in two rows. The top row is the parameter scale, with the elicited density and the two bounds we asked for. The bottom row is the data scale: the density of the quantity we reasoned about (demand per day, level multiplier after one year, noise fraction), with the bounds mapped through the same transform.

The bottom row needs no draws. If \(\psi = g(\theta)\) is a monotone transform of the parameter \(\theta\) with density \(f\), the density of \(\psi\) is \(f(g^{-1}(\psi)) \, |dg^{-1} / d\psi|\) by the change of variables formula. We write the inverse transforms in jax.numpy and take the derivative with jax.grad.

intercept_bounds = (float(np.log(10)), float(np.log(300)))
drift_bounds = (
    float(np.log(1.04) / np.sqrt(setup.t_obs)),
    float(np.log(2) / np.sqrt(setup.t_obs)),
)
concentration_bounds = (4.0, 100.0)
# the prior mass between the bounds, fixed here and not tied to the plotted bands
ELICITATION_MASS = 0.94

intercept_prior = pz.maxent(
    pz.Normal(), *intercept_bounds, mass=ELICITATION_MASS, plot=False
)
drift_prior = pz.maxent(
    pz.LogNormal(), *drift_bounds, mass=ELICITATION_MASS, plot=False
)
concentration_prior = pz.maxent(
    pz.Gamma(), *concentration_bounds, mass=ELICITATION_MASS, plot=False
)

priors = {
    "intercept": (intercept_prior, intercept_bounds),
    "drift_scale": (drift_prior, drift_bounds),
    "concentration": (concentration_prior, concentration_bounds),
}
# the quantity on the data scale that each prior is about, with the inverse transform
data_scale = {
    "intercept": ("demand per day at the list price", jnp.exp, jnp.log),
    "drift_scale": (
        "level multiplier after one year",
        lambda x: jnp.exp(x * jnp.sqrt(setup.t_obs)),
        lambda y: jnp.log(y) / jnp.sqrt(setup.t_obs),
    ),
    "concentration": (
        "extra noise as a fraction of the level",
        lambda x: 1 / jnp.sqrt(x),
        lambda y: 1 / y**2,
    ),
}


def density_on_data_scale(
    prior: pz.distributions.distributions.Continuous,
    inverse: Callable[[Float[Array, " grid"]], Float[Array, " grid"]],
    grid: Float[Array, " grid"],
) -> Float[Array, " grid"]:
    """Density of a monotone transform of the prior, by change of variables."""
    jacobian = jnp.abs(jax.vmap(jax.grad(inverse))(grid))
    return prior.pdf(inverse(grid)) * jacobian


fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 8), layout="constrained")

for j, (name, (prior, bounds)) in enumerate(priors.items()):
    label, transform, inverse = data_scale[name]
    # parameter scale: the elicited density and the bounds
    ax = axes[0, j]
    x = np.linspace(*prior.ppf([0.001, 0.999]), 500)
    ax.plot(x, prior.pdf(x), color="C0", label="elicited density")
    ax.axvline(bounds[0], color="black", linestyle="--", label="elicited bounds")
    ax.axvline(bounds[1], color="black", linestyle="--")
    ax.set(title=f"{name}\n{prior}", yticks=[])
    # data scale: the same density after the transform, from near zero to the
    # larger of the 99th percentile and the upper bound
    ax = axes[1, j]
    bounds_data = transform(jnp.asarray(bounds))
    upper = 1.05 * max(
        float(transform(prior.ppf([0.01, 0.99])).max()), float(bounds_data.max())
    )
    grid = jnp.linspace(upper / 500, upper, 500)
    ax.plot(grid, density_on_data_scale(prior, inverse, grid), color="C0")
    for bound in bounds_data:
        ax.axvline(bound, color="black", linestyle="--")
    ax.set(title=label, yticks=[])

axes[0, 0].set(ylabel="parameter scale")
axes[1, 0].set(ylabel="data scale")
fig.legend(*axes[0, 0].get_legend_handles_labels(), loc="outside lower center", ncol=2)
fig.suptitle("Elicited Priors", fontsize=18, fontweight="bold");

The elicited priors are a \(\text{Normal}(4.0, 0.9)\) for the intercept, a \(\text{LogNormal}(-4.3, 0.6)\) for the drift scale and a \(\text{Gamma}\) with shape \(2.8\) and rate \(0.06\) for the concentration.

The top row shows the three elicited densities with their bounds, and the mass between the bounds is \(94\%\) in each panel. The bottom row shows what these priors say on the data scale. The demand per day spreads between a few units and a few hundred, with most of the mass below \(100\) units. The level multiplier after one year sits between no move and a doubling, with a tail beyond. The noise fraction sits mostly between \(0.1\) and \(0.3\), inside its bounds of \(0.1\) and \(0.5\).

Note that pz.maxent fixes the total mass between the bounds and not the two tails, so the Gamma puts almost all of its remaining \(6\%\) above \(100\), which is the tail below \(0.1\) in the bottom right panel. The LogNormal does the same: almost all of its remaining \(6\%\) sits above the upper bound, beyond a doubling of the level. The true values of the generator, \(\log 60 \approx 4.1\), \(0.01\) and \(30\), sit inside all three priors.

Model Specification

Recall that a numpyro_forecast model is a plain function with signature (covariates, data=None). In addition, it has the following components:

• Horizon.from_data reads the observed length from data and the full duration from covariates.

• innovations samples the level increments under a time plate over the history and a time_future plate over the horizon.

• predict registers obs over the training window and forecast over the horizon with the negative binomial likelihood.

We define the model with these priors and render it on the training window to see its structure. The last column of covariates is the log relative price.

def demand_model(
    covariates: Float[Array, "duration features"],
    data: Int[Array, "t_obs 1"] | None = None,
) -> None:
    """Local level plus weekly seasonality plus price, negative binomial likelihood."""
    horizon = Horizon.from_data(covariates, data)
    fourier, log_price = covariates[..., :-1], covariates[..., -1]
    num_features = fourier.shape[-1]
    intercept = numpyro.sample(
        "intercept", dist.Normal(intercept_prior.mu, intercept_prior.sigma)
    )
    beta = numpyro.sample(
        "beta", dist.Normal(0.0, 0.3).expand([num_features]).to_event(1)
    )
    elasticity = numpyro.sample("elasticity", dist.Normal(-1.0, 1.0))
    drift_scale = numpyro.sample(
        "drift_scale", dist.LogNormal(drift_prior.mu, drift_prior.sigma)
    )
    concentration = numpyro.sample(
        "concentration",
        dist.Gamma(concentration_prior.alpha, concentration_prior.beta),
    )
    # non-centered increments: the sampler sees standardized drifts
    drift = innovations(
        horizon, "drift", dist.Normal(0.0, drift_scale), reparam=LocScaleReparam(0)
    )
    level = jnp.cumsum(drift, axis=-2)
    regression = (fourier @ beta + elasticity * log_price)[..., None]
    log_rate = intercept + level + regression

    def likelihood(eta: Float[Array, "..."]) -> Distribution:
        return dist.NegativeBinomial2(jnp.exp(eta), concentration)

    predict(horizon, likelihood, log_rate)


numpyro.render_model(
    demand_model,
    model_args=(covariates[: setup.t_obs], data),
    render_distributions=True,
    render_params=True,
)

Inside the time plate we see drift_decentered, the standardized increments that the LocScaleReparam samples, and the deterministic drift that rescales them by drift_scale. On the training window there is no horizon, so the time_future plate does not appear.

Data Preparation

The function to_datatree stores the forecast draws as obs in the predictions group, shaped (chain, draw, time, obs_dim). We flatten the sample dimensions into a scenario matrix of shape \((N, H)\), where each row is one joint draw of the \(15\) days. These rows are the scenarios of the Decision Problem section.

We then shuffle the rows and split them into two sets. The first \(1{,}000\) scenarios are the optimization set, and every plan is built on them. The remaining \(3{,}000\) scenarios are the evaluation set, which no plan sees.

We need the split because a plan fitted to a finite set of scenarios adapts to that set, so its loss on the same scenarios is usually optimistic. The held-out scenarios tell us how the plan does on new draws from the same forecast. We also keep the mean demand path of the optimization set, which several figures use.

Note that the plots keep the day index of the series, so the horizon runs from day \(364\) to day \(378\). This corresponds to the \(15\)-day horizon.

N_SCENARIOS_OPT = 1_000

scenarios = jnp.asarray(
    idata["predictions"]["obs"].values.reshape(-1, setup.horizon), dtype=jnp.float32
)
rng_key, key_split = random.split(rng_key)
scenarios = scenarios[random.permutation(key_split, scenarios.shape[0])]
demand_opt, demand_eval = scenarios[:N_SCENARIOS_OPT], scenarios[N_SCENARIOS_OPT:]
mean_demand_opt = demand_opt.mean(axis=0)

Scenario Visualization

We plot the forecast as a fan over the horizon, next to the last \(60\) days of history. The shaded column is the promotion week in the horizon.

def shade_promotion(ax: plt.Axes) -> None:
    """Shade the promotion days of the horizon."""
    ax.axvspan(
        promo_days_future[0] - 0.5,
        promo_days_future[-1] + 0.5,
        color="C5",
        alpha=0.25,
        label="Promotion week",
    )


forecast_tree = az.from_dict(
    {
        "posterior_predictive": {
            "obs": idata["predictions"]["obs"].sel(obs_dim=0).values
        },
        "observed_data": {"obs": demand_future_true},
        "constant_data": {"day": days_future},
    },
    coords={"time": days_future},
    dims={"obs": ["time"], "day": ["time"]},
)

pc = az.plot_lm(
    forecast_tree,
    x="day",
    y="obs",
    y_obs="obs",
    smooth=False,
    ci_prob=(0.5, CI_PROB),
    point_estimate="mean",
    visuals={
        "ci_band": {"color": "C2"},
        "pe_line": {"color": "C2", "alpha": 1.0, "linewidth": 2},
        "observed_scatter": False,
    },
    figure_kwargs={"figsize": (12, 6)},
)
lighten_bands(pc, "day")
ax = pc.viz["plot"]["day"].item()
pc.viz["ci_band"]["day"].sel(prob=CI_PROB).item().set_label(f"{CI_PROB:.0%} ETI")
pc.viz["ci_band"]["day"].sel(prob=0.5).item().set_label("50% ETI")
pc.viz["pe_line"]["day"].item().set_label("Forecast mean")
last_two_months = slice(setup.t_obs - 60, setup.t_obs)
ax.plot(
    days[last_two_months],
    data[last_two_months, 0],
    "o-",
    markersize=3,
    color="C0",
    label="Observed history",
)
ax.plot(
    days_future,
    demand_future_true,
    "o-",
    markersize=3,
    color="C1",
    label="True future",
)
ax.axvline(setup.t_obs - 0.5, color="black", linestyle="--", linewidth=1)
shade_promotion(ax)
ax.set(xlabel="day", ylabel="units", title="")
fig = pc.viz["figure"].item()
fig.legend(loc="outside lower center", ncol=3)
fig.suptitle("Forecast Fan for the Planning Horizon", fontsize=18, fontweight="bold");

The fan follows the weekly cycle and jumps in the promotion week because the model knows the price path. It widens slowly with the horizon because the level keeps drifting. The true future stays inside the wide band. Next, we look at \(30\) individual scenarios from the optimization set, the paths that the optimizer sees.

fig, ax = plt.subplots(figsize=(12, 6))
lines = ax.plot(days_future, demand_opt[:30].T, color="C2", alpha=0.3, linewidth=1)
lines[0].set_label("Scenario")
ax.plot(
    days_future,
    mean_demand_opt,
    color="C2",
    linewidth=3,
    label="Forecast mean",
)
shade_promotion(ax)
fig.legend(loc="outside lower center", ncol=3)
ax.set(xlabel="day", ylabel="units")
fig.suptitle("Joint Demand Scenarios", fontsize=18, fontweight="bold");

Each path is one joint draw of the \(15\) days. A plan is scored against each path, and the decision rule aggregates the scores.

The paths have to be joint draws for two reasons. First, the stock that a scenario carries into day \(t\) is what is left of the earlier days of the same scenario. Second, the days of one path share the same level and the same parameters.

Two consecutive days are only weakly linked, but the stock carried into day \(t\) depends on the demand of all the earlier days, and a path with high demand on the earlier days tends to have high demand on day \(t\) too. In the Out of Sample Evaluation section we check what happens to a plan built on scenarios without this link.

Chance-Constrained Orders

We run the forward pass on the joint scenarios with the quantile level \(1 - \alpha\). Then we check the guarantee of the quantile on the \(N = 1{,}000\) optimization scenarios: on every day, the number of scenarios with a lost sale is at most \(\alpha N = 50\).

Demand comes in whole units, while the orders of the expected cost plan below land up to a few ten-thousandths of a unit away from whole numbers. Without a tolerance, a day where the demand matches the intended stock can count as a stockout, so we count a stockout only when the lost sales exceed \(10^{-3}\) units.

STOCKOUT_TOL = 1e-3  # lost sales below this are rounding noise of the solver


def is_stockout(sim: Simulation) -> Bool[Array, "*batch horizon"]:
    """Days with lost sales above the tolerance."""
    return sim.lost_sales > STOCKOUT_TOL


orders_cc = forward_pass(demand_opt, quantile_shortfall(1 - ALPHA), config)

stockout_scenarios_cc = is_stockout(
    simulate_scenarios(orders_cc, demand_opt, config)
).sum(axis=0)
assert int(stockout_scenarios_cc.max()) <= tail_size(N_SCENARIOS_OPT)

Expected Cost Orders

The expected cost plan and the CVaR plan minimize the same objective, \(\text{E}[L] + \lambda \, \text{CVaR}_\alpha(L)\), with \(\lambda = 0\) for the first and \(\lambda = 1\) for the second, so one solver covers both. The objective is convex in \(a\), which we show in the Appendix, and the only constraint is \(a \geq 0\).

We use projected gradient descent with Adam. Each step moves the orders a step down the gradient of the objective and then clips them at zero, which is the Euclidean projection onto the feasible set. The loop is a jax.lax.scan inside jax.jit, so later on a sweep over several values of \(\lambda\) is a single jax.vmap.

The objective is piecewise linear, so its gradient is a subgradient, and a step size that does not shrink keeps crossing the kinks. We use a cosine decay schedule that takes the step size from one unit to zero over the \(1{,}500\) steps of the run, so the early steps move fast and the late steps settle the iterates.

The function returns the orders, the trace of the objective and the step size at every step. We start from the mean forecast plan and solve with \(\lambda = 0\).

def objective(
    orders: Float[Array, " horizon"],
    demand: Float[Array, "scenarios horizon"],
    lam: ScalarLike,
    config: InventoryConfig,
) -> ScalarLike:
    """Expected loss plus lam times the CVaR of the loss over the scenarios."""
    losses = simulate_scenarios(orders, demand, config).loss
    return losses.mean() + lam * cvar(losses)


def solve(
    demand: Float[Array, "scenarios horizon"],
    lam: ScalarLike,
    orders_init: Float[Array, " horizon"],
    config: InventoryConfig,
    n_steps: int = 1_500,
    learning_rate: float = 1.0,
) -> tuple[
    Float[Array, " horizon"], Float[Array, " n_steps"], Float[Array, " n_steps"]
]:
    """Projected Adam on the expected cost plus lam times CVaR, as a lax.scan."""
    # the step size decays from learning_rate to zero along a cosine curve
    schedule = optax.cosine_decay_schedule(
        init_value=learning_rate, decay_steps=n_steps
    )
    optimizer = optax.adam(learning_rate=schedule)
    value_and_grad = jax.value_and_grad(objective)

    def step(
        carry: tuple[Float[Array, " horizon"], optax.OptState], _: None
    ) -> tuple[tuple[Float[Array, " horizon"], optax.OptState], ScalarLike]:
        orders, opt_state = carry
        value, grads = value_and_grad(orders, demand, lam, config)
        updates, opt_state = optimizer.update(grads, opt_state, orders)
        orders = jnp.maximum(orders + updates, 0.0)  # projection onto a >= 0
        return (orders, opt_state), value

    (orders, _), values = jax.lax.scan(
        step, (orders_init, optimizer.init(orders_init)), None, length=n_steps
    )
    step_sizes = jax.vmap(schedule)(jnp.arange(n_steps))
    return orders, values, step_sizes


solve_jit = jax.jit(solve, static_argnames=("config", "n_steps"))

orders_saa, trace_saa, step_sizes_saa = solve_jit(demand_opt, 0.0, orders_mean, config)

We also build the myopic newsvendor plan with the forward pass at the critical fractile and compare the two plans.

orders_newsvendor = forward_pass(
    demand_opt, quantile_shortfall(critical_fractile), config
)
cost_saa = objective(orders_saa, demand_opt, 0.0, config)
cost_newsvendor = objective(orders_newsvendor, demand_opt, 0.0, config)

fig, ax = plt.subplots()
ax.plot(
    days_future,
    orders_saa,
    "o-",
    color="C2",
    label=f"Expected cost plan, in-sample expected loss {cost_saa:.0f}",
)
ax.plot(
    days_future,
    orders_newsvendor,
    "s--",
    color="black",
    label=f"Myopic newsvendor plan, in-sample expected loss {cost_newsvendor:.0f}",
)
shade_promotion(ax)
fig.legend(loc="outside lower center", ncol=2)
ax.set(xlabel="day", ylabel="units ordered")
fig.suptitle(
    "Expected Cost Plan Against the Myopic Newsvendor Plan",
    fontsize=18,
    fontweight="bold",
);

The two plans follow the same shape, and the myopic plan orders more on \(12\) of the \(15\) days and in total, \(1{,}410\) units against \(1{,}355\). The extra stock costs \(11\%\) in expected loss, \(158\) against \(143\) on the optimization scenarios.

Here is the intuition for the larger orders of the myopic plan. The single-period rule charges a leftover unit one day of holding cost. In the \(15\)-day problem a leftover unit stays on the shelf and pays \(h\) on every later day until it sells. It is also left over exactly in the slow scenarios, where the level tends to be low, so it takes several days to sell. The optimizer sees this extra cost and orders less.

The optimizer orders more than the myopic rule on the last two days only, by \(5\) and \(7\) units. There, few later days remain, so both plans want a similar stock available, and the optimizer carries less stock into those days.

CVaR Orders

Now we solve the CVaR plan with \(\lambda = 1\): one unit of average loss in the worst \(5\%\) of scenarios counts as much as one unit of expected loss.

orders_cvar, trace_cvar, _ = solve_jit(demand_opt, 1.0, orders_mean, config)

Here are the optimizer traces for both solves: the objective on top and the step size of the schedule below it.

steps = np.arange(1, trace_saa.shape[0] + 1)

fig, axes = plt.subplots(nrows=2, figsize=(10, 7), sharex=True, layout="constrained")
ax = axes[0]
ax.plot(steps, trace_saa, color="C2", label=r"$\lambda = 0$ (expected cost)")
ax.plot(steps, trace_cvar, color="C3", label=r"$\lambda = 1$ (expected cost + CVaR)")
ax.set(ylabel="objective", xscale="log")
ax = axes[1]
ax.plot(steps, step_sizes_saa, color="black")  # the schedule is the same for both
ax.set(xlabel="step", ylabel="step size (both solves)")
fig.legend(loc="outside lower center", ncol=2)
fig.suptitle("Optimizer Traces", fontsize=18, fontweight="bold");

Both objectives drop within the first \(50\) steps and change by less than \(1\%\) after that, settling at \(143\) for the expected cost and \(469\) for the expected cost plus CVaR. The step size decays to zero over the run, so the iterates stop moving while the objective is flat, and the Appendix checks that these values are the optima.

We now have the four plans, so let’s collect them in a list of Plan tuples and compare them. We fix one color per plan and use it in every figure that follows. The shaded column is again the promotion week.

class Plan(NamedTuple):
    """One order plan with its name."""

    name: str
    orders: Float[Array, " horizon"]


plans = [
    Plan("Mean forecast", orders_mean),
    Plan("Chance constrained", orders_cc),
    Plan("Expected cost", orders_saa),
    Plan("CVaR", orders_cvar),
]
plan_names = [plan.name for plan in plans]
PLAN_COLORS = {name: f"C{i}" for i, name in enumerate(plan_names)}
orders_long = pl.DataFrame(
    {
        "plan": np.repeat(plan_names, setup.horizon),
        "day": np.tile(days_future, len(plans)),
        "orders": np.concatenate([plan.orders for plan in plans]),
    }
)

fig, ax = plt.subplots(figsize=(12, 7))
sns.lineplot(
    data=orders_long,
    x="day",
    y="orders",
    hue="plan",
    hue_order=plan_names,
    palette=PLAN_COLORS,
    marker="o",
    ax=ax,
)
ax.plot(
    days_future,
    mean_demand_opt,
    color="black",
    linestyle="--",
    label="Mean demand",
)
shade_promotion(ax)
ax.get_legend().remove()
fig.legend(loc="outside lower center", ncol=3)
ax.set(xlabel="day", ylabel="units ordered")
fig.suptitle("Order Plans", fontsize=18, fontweight="bold");

Here is what we see.

• Mean forecast. The plan follows the mean demand. Its first order is \(23\) units, because the \(50\) units of initial stock cover most of the mean demand of \(73\) on day one.

• Buffers. The three other plans order a buffer above the mean demand, and every plan orders more in the promotion week. On day one the three plans order \(44\) to \(54\) units, which is below the mean demand of \(73\) and looks like no buffer at all. The buffer is in the stock available, which is the initial stock plus the order: \(94\) to \(104\) units on day one, so \(21\) to \(31\) units above the mean demand. On day one the initial stock provides part of the buffer. From day two on, the chance-constrained plan orders above the mean demand on every day, the expected cost and CVaR plans do so on all days but one, and the stock carried over adds to the buffer.

• Chance constrained. This plan orders the most in the first week and carries the largest cumulative order through day \(11\).

• Expected cost. This plan orders the least of the three buffered plans, \(1{,}355\) units against \(1{,}406\) and \(1{,}431\), and stays closest to the mean demand.

• CVaR. This plan orders the most in the second week, where the high-demand scenarios of the tail separate from the rest, and from day \(12\) on its cumulative orders are the largest of the four. Its path is less smooth than the chance-constrained one, because the \(50\) scenarios of the tail shape it.

Key Findings

• Rule and forecast. The forecast fixes the scenarios, and the decision rule applied to them is a separate choice: expected cost, a service level or a weight on the worst outcomes. Plans built from the same scenarios under different rules differ by tens of percent in cost and in tail loss.

• Distribution, not point. The loss has a kink at the stock, so the cost on the mean path underestimates the expected cost, and when the two errors have different prices the best stock is not the mean. A plan built on a point forecast is the most expensive plan on every metric we report, and only a probabilistic forecast gives us the alternatives.

• Multi-period dynamics. Carried stock couples the days, so single-period formulas applied day by day miss the multi-period optimum (here they overstock, at an \(11\%\) higher expected loss on the held-out scenarios). A loss written as a simulator and optimized directly does not depend on closed forms that hold only for one period.

• Frequency and severity. A chance constraint controls how often we run out, and CVaR controls how bad the worst outcomes are. A service level needs the constraint, a bound on bad outcomes needs the tail weight, and the frontier over \(\lambda\) shows the price of one in terms of the other. Beyond some weight the stockout frequency stops improving, because the tail fills with low-demand scenarios that pay holding cost.

• Joint scenarios. Anything with a term that depends on the demand of several days or products at once (carried stock, or a CVaR or a chance constraint on the total loss) needs joint draws, not per-day quantiles. A posterior predictive gives them directly, and a plan built on shuffled days landed at \(7\%\) per day against a \(5\%\) target.

• Out-of-sample scores. Held-out scenarios score a plan on draws it was not fitted to, which is the fair comparison between plans built with different rules. A cross check against an exact solver, in the Appendix, confirms that the gradient solver reaches the optimum.

Model Limitations

• No recourse. All orders are fixed at the start of the horizon. A rolling horizon that re-forecasts and re-solves each day, and commits only the next order, is the practical version of this problem.

• No lead time. Orders arrive on the day we planned them for. The VN2 setting has a two-week delivery lead time, which moves the decision to the stock projected at the delivery date.

• Unconstrained demand history. The forecasting model sees the demand itself. In a real lost-sales setting the history is sales, which are censored on the days the shelf was empty. The VN2 report handles this with stockout-aware features, Caron (2023) with a censored likelihood, and the censored demand example of numpyro_forecast shows the same likelihood in a forecasting model.

• Same posterior throughout. The held-out scenarios come from the posterior that produced the optimization set, so they check the plans against new draws from the same forecast (the expected cost plan scores \(143\) in sample and \(142\) held out) but not against model misspecification. A backtest on real data would show it.

• Linear costs. A lost sale costs the same margin on every day, and leftovers keep their value at the end of the horizon. Perishable goods, order minimums and fixed ordering costs change the simulator, and integer order quantities change the solver.

• Small tail sample. With \(\alpha = 0.05\) and \(1{,}000\) scenarios, the CVaR objective is the mean of \(50\) losses, and those \(50\) scenarios shape the plan. A real application needs more scenarios, or a check of the plan across resampled scenario sets.

Recommendations

1. Elicit the costs. The holding cost and the margin set the critical fractile and every cost-based plan with it. Agree on them with the business before any optimization, and state what a leftover is worth at the end of the horizon.

2. Pick from the frontier. Do not fix a service level or a tail weight by default. Report the expected loss, the tail loss and the stockout frequency of several plans and let the business pick the point.

3. Backtest the plan. Score the forecast and the plan together on a rolling window of real data, with the same loss, so that model misspecification shows up as cost.

Convex Solver Cross Check

The projected gradient solver is convenient because it stays in JAX, but it is a first-order method on a non-smooth objective. The loss has a kink wherever a positive part switches, that is, wherever the stock available equals the demand of a day in one of the scenarios. The gradient jumps at these kinks, so a first-order method can stop short of the exact optimum. In this appendix we confirm that it does not, with a solver that reaches the exact optimum of a linear program. Two things make the objective a linear program:

• The loss is convex and piecewise linear in the plan, since it is a positive combination of the daily stocks plus a linear term, and each stock is a maximum of affine functions of the plan.

• The CVaR has a representation with a threshold variable: an extra decision variable \(\zeta\), a loss threshold, such that the CVaR is the minimum over \(\zeta\) of a convex piecewise linear function of \(\zeta\). Its positive parts become linear constraints once we add one auxiliary variable per scenario.

We start with the second one. The representation of the conditional value at risk of Rockafellar and Uryasev (2000) is

\[\text{CVaR}_\alpha(L) = \min_{\zeta} \left\{ \zeta + \frac{1}{\alpha} \, \text{E}\left[(L - \zeta)^{+}\right] \right\},\]

where a minimizing \(\zeta\) is the value at risk. Let’s verify this representation on the losses of the expected cost plan over the optimization scenarios, in three steps.

• We write the function inside the braces as ru_objective, with the expectation replaced by the average over the scenarios.

• We evaluate it on a grid of thresholds, with jax.vmap over the grid.

• We compare its minimum with the value that cvar computes from the \(k\) largest losses, and the location of the minimum with the value at risk.

def ru_objective(
    threshold: ScalarLike, losses: Float[Array, " scenarios"]
) -> ScalarLike:
    """The threshold zeta plus the scenario average of (L - zeta)^+ over alpha."""
    return threshold + jnp.maximum(losses - threshold, 0.0).mean() / ALPHA


losses_saa_opt = simulate_scenarios(orders_saa, demand_opt, config).loss
threshold_grid = jnp.linspace(0.0, losses_saa_opt.max(), 400)
ru_values = jax.vmap(ru_objective, in_axes=(0, None))(threshold_grid, losses_saa_opt)
k = tail_size(N_SCENARIOS_OPT)

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(
    threshold_grid,
    ru_values,
    color="C0",
    linewidth=2,
    label=r"$\zeta + \text{E}[(L - \zeta)^+] / \alpha$",
)
ax.axhline(
    float(cvar(losses_saa_opt)),
    color="C3",
    linewidth=2,
    label=f"CVaR (mean of the {k} largest losses)",
)
ax.axvline(
    float(jnp.quantile(losses_saa_opt, 1 - ALPHA, method="inverted_cdf")),
    color="C1",
    linewidth=2,
    label="VaR",
)
fig.legend(loc="outside lower center", ncol=3)
ax.set(
    xlabel=r"threshold $\zeta$",
    ylabel="value",
    ylim=(0, 3 * float(cvar(losses_saa_opt))),
)
fig.suptitle(
    "The Rockafellar and Uryasev Representation", fontsize=18, fontweight="bold"
);

The curve is convex in \(\zeta\), its minimum sits at the value at risk, and the minimum value is the conditional value at risk that cvar computes with top_k. So the linear program minimizes the same objective as solve. It minimizes over \(\zeta\) and over the plan at the same time, with one auxiliary variable \(u_n \geq (L_n - \zeta)^{+}\) per scenario, where \(L_n = L(a, D^{(n)})\) is the loss of scenario \(n\).

Now we turn to the loss. The lost sales of a day are the demand minus the stock available, plus the leftover:

\[U_t = (D_t - y_t)^{+} = D_t - y_t + (y_t - D_t)^{+} = D_t - y_t + I_t.\]

With \(y_t = I_{t-1} + a_t\), we sum \(h \, I_t + p \, U_t\) over the days. The stocks \(I_{t-1}\) and \(I_t\) telescope, and we get

\[L(a, D) = h \sum_{t=1}^{H} I_t + p \, I_H + p \left( C_H - I_0 - A_H \right),\]

where \(A_H = \sum_t a_t\) is the total order and \(C_H = \sum_t D_t\) the total demand. The last term is linear in the plan. Each stock is the positive part of an affine function of the previous stock and the order,

\[I_t = \max\left(I_{t-1} + a_t - D_t, \, 0\right),\]

so by induction on \(t\) it is convex in \(a\), and the loss is convex in \(a\). The mean of the \(k\) largest losses is convex and non-decreasing in the losses, so the objective is convex too.

In the linear program we replace every stock \(I_t^{(n)}\) by a variable with two constraints,

\[I_t^{(n)} \geq I_{t-1}^{(n)} + a_t - D_t^{(n)}, \qquad I_t^{(n)} \geq 0.\]

Every stock variable has a positive coefficient in the objective, so the solver pushes each one down to the larger of the two bounds, which is the positive part we want. We write the program in three pieces: stock_constraints builds the stock bounds, scenario_losses the loss of every scenario, and cvar_epigraph the CVaR term with its auxiliary variables and constraint. Then solve_lp assembles them and hands the problem, with \(N \times H\) stock variables, to HiGHS through cvxpy, and we solve it for the two weights.

def stock_constraints(
    stock: cp.Variable,
    orders: cp.Variable,
    demand: Float[np.ndarray, "scenarios horizon"],
    config: InventoryConfig,
) -> list[cp.Constraint]:
    """I_t >= I_{t-1} + a_t - D_t for every scenario and day, tight at the optimum."""
    horizon = demand.shape[1]
    first_day = stock[:, 0] >= config.initial_inventory + orders[0] - demand[:, 0]
    later_days = [
        stock[:, t] >= stock[:, t - 1] + orders[t] - demand[:, t]
        for t in range(1, horizon)
    ]
    return [first_day, *later_days]


def scenario_losses(
    stock: cp.Variable,
    orders: cp.Variable,
    demand: Float[np.ndarray, "scenarios horizon"],
    config: InventoryConfig,
) -> cp.Expression:
    """The loss h sum_t I_t + p I_H + p (C_H - I_0 - A_H) of every scenario."""
    return (
        config.holding_cost * cp.sum(stock, axis=1)
        + config.lost_sales_cost * stock[:, -1]
        + config.lost_sales_cost
        * (demand.sum(axis=1) - config.initial_inventory - cp.sum(orders))
    )


def cvar_epigraph(
    losses: cp.Expression, num_scenarios: int
) -> tuple[cp.Expression, list[cp.Constraint]]:
    """The Rockafellar and Uryasev term zeta + sum_n u_n / k with its constraint."""
    threshold = cp.Variable()  # zeta
    excess = cp.Variable(num_scenarios, nonneg=True)  # u_n >= (L_n - zeta)^+
    # the same integer tail size as `cvar`, so both objectives agree
    term = threshold + cp.sum(excess) / tail_size(num_scenarios)
    return term, [excess >= losses - threshold]


def solve_lp(
    demand: Float[Array, "scenarios horizon"], lam: float, config: InventoryConfig
) -> Float[Array, " horizon"]:
    """The same objective as `solve`, as a linear program."""
    demand = np.asarray(demand, dtype=np.float64)
    num_scenarios, horizon = demand.shape
    orders = cp.Variable(horizon, nonneg=True)  # the plan a
    stock = cp.Variable((num_scenarios, horizon), nonneg=True)  # I_t of every scenario

    # `cp.sum` checks its shape by summing an uninitialized array, which can warn.
    # Such atoms are built here and again when `solve` canonicalizes the problem.
    with warnings.catch_warnings():
        warnings.filterwarnings(
            "ignore",
            message=".*encountered in reduce",
            category=RuntimeWarning,
        )
        constraints = stock_constraints(stock, orders, demand, config)
        losses = scenario_losses(stock, orders, demand, config)
        cvar_term, cvar_constraints = cvar_epigraph(losses, num_scenarios)
        lp = cp.Problem(
            cp.Minimize(cp.sum(losses) / num_scenarios + lam * cvar_term),
            [*constraints, *cvar_constraints],
        )
        lp.solve(solver=cp.HIGHS)

    assert lp.status == cp.OPTIMAL, lp.status
    return jnp.asarray(orders.value, dtype=jnp.float32)


orders_lp_saa = solve_lp(demand_opt, 0.0, config)
orders_lp_cvar = solve_lp(demand_opt, 1.0, config)

We overlay the order plans of both solvers and evaluate objective at both solutions, so that the two values come from the same function. The panel titles carry the relative gap and the largest difference between the two order vectors.

fig, axes = plt.subplots(ncols=2, figsize=(14, 5), sharey=True, layout="constrained")

for ax, (name, lam, orders_jax, orders_lp, color) in zip(
    axes,
    [
        ("Expected cost", 0.0, orders_saa, orders_lp_saa, "C2"),
        ("CVaR", 1.0, orders_cvar, orders_lp_cvar, "C3"),
    ],
    strict=True,
):
    value_jax = float(objective(orders_jax, demand_opt, lam, config))
    value_lp = float(objective(orders_lp, demand_opt, lam, config))
    relative_gap = value_jax / value_lp - 1
    max_order_diff = float(jnp.abs(orders_jax - orders_lp).max())
    ax.plot(
        days_future,
        orders_jax,
        color=color,
        linewidth=6,
        alpha=0.5,
        label="Projected Adam (JAX)",
    )
    ax.plot(
        days_future,
        orders_lp,
        "o",
        color="black",
        markersize=4,
        label="Linear program (cvxpy)",
    )
    shade_promotion(ax)
    ax.set(
        xlabel="day",
        title=f"{name}\nobjective gap {relative_gap:.1e}, "
        f"max order difference {max_order_diff:.1f} units",
    )
    ax.legend(loc="upper left")

axes[0].set(ylabel="units ordered")
fig.suptitle("JAX Solver Against the Linear Program", fontsize=18, fontweight="bold");

The two solvers agree! The relative objective gap is \(1 \times 10^{-7}\) for the expected cost plan and zero to single precision for the CVaR plan, and the orders differ by at most \(0.11\) units on orders between \(44\) and \(152\) units, so the two solvers land on the same plan.

The LP orders of the expected cost plan are whole numbers on every day, so requiring whole-unit orders would not change that plan. We use the JAX solver to vectorize over weights or over products, and a linear or mixed-integer program when the problem gains coupling constraints or integer decisions.

PyMC and PyTensor Implementation

In this appendix we reproduce the held-out scores of the Out of Sample Evaluation section with a second implementation: the forecast with the functional API of pymc_forecast, and the forward pass, the simulator and the solver in PyTensor, with the optimizers of pytensor-ml. We use the same data, the same priors, the same costs and the same split into \(1{,}000\) optimization and \(3{,}000\) evaluation scenarios. We score the plans with the summarize function of this notebook and check the solver against its linear program.

The concepts are the ones of the main text, so we do not repeat them here. The full walkthrough of this implementation, with the prior and posterior predictive checks, the parameter recovery and every figure of the main text, is in the notebook forecast_to_decision_pymc.ipynb. Here we show the key ingredients only.

We import the PyMC stack here and not at the top of the notebook, so that the main text runs without it. We also silence the sampler logs, since we read the diagnostics from the outputs below. The pymc_forecast functions Horizon, innovations and predict share their names with the numpyro_forecast ones, so we keep the module prefix pmf.

import logging
import math

import pymc as pm
import pymc_forecast as pmf
import pytensor
import pytensor.tensor as pt
from pymc_forecast.data import FUTURE_DIM, TIME_DIM
from pytensor_ml.optim import adam, chain, compile_train, cosine_schedule, to_updates
from pytensor_ml.params import trainable

logging.getLogger("pymc").setLevel(logging.ERROR)
logging.getLogger("pytensor").setLevel(logging.ERROR)

rng = np.random.default_rng(seed=seed)

A pymc_forecast model takes its data and covariates as labeled xarray arrays with a time dim, so we wrap the arrays of the Generate Data section. The covariate matrix keeps the Fourier columns and the log relative price as the last column.

FOURIER_LABELS = [f"fourier_{i}" for i in range(covariates.shape[1] - 1)]

covariates_pymc = xr.DataArray(
    np.asarray(covariates, dtype="float64"),
    dims=(TIME_DIM, "covariate"),
    coords={TIME_DIM: days, "covariate": [*FOURIER_LABELS, "log_price"]},
)
covariates_train_pymc = covariates_pymc.isel({TIME_DIM: slice(None, setup.t_obs)})
data_pymc = xr.DataArray(
    np.asarray(data[:, 0], dtype="float64"),
    dims=(TIME_DIM,),
    coords={TIME_DIM: days[: setup.t_obs]},
)

The model is a plain function (covariates, data=None) with the same structure and the same priors as the NumPyro model of the Forecasting Model section. Three functions of pymc_forecast do the same job as their numpyro_forecast namesakes: pmf.Horizon.from_data reads the observed length from data and the full duration from covariates, pmf.innovations registers the level increments over the history and the horizon, and pmf.predict registers the observed and the forecast variables with the negative binomial likelihood.

def demand_model_pymc(
    covariates: xr.DataArray, data: xr.DataArray | None = None
) -> None:
    """Local level plus weekly seasonality plus price, negative binomial likelihood."""
    h = pmf.Horizon.from_data(covariates, data)
    fourier_design = np.ascontiguousarray(
        covariates.sel(covariate=FOURIER_LABELS).values
    )
    log_price = covariates.sel(covariate="log_price").values
    pm.modelcontext(None).add_coord("fourier", FOURIER_LABELS)

    intercept = pm.Normal("intercept", intercept_prior.mu, intercept_prior.sigma)
    beta = pm.Normal("beta", 0.0, 0.3, dims="fourier")
    elasticity = pm.Normal("elasticity", -1.0, 1.0)
    drift_scale = pm.LogNormal("drift_scale", drift_prior.mu, drift_prior.sigma)
    concentration = pm.Gamma(
        "concentration", alpha=concentration_prior.alpha, beta=concentration_prior.beta
    )
    # non-centered increments: the sampler sees standardized drifts
    drift_raw = pmf.innovations(h, "drift_raw", pm.Normal.dist(0.0, 1.0))
    level = pt.cumsum(drift_raw * drift_scale)
    log_rate = intercept + level + pt.dot(fourier_design, beta) + elasticity * log_price

    def likelihood(name: str, eta: pt.TensorVariable, dims: tuple[str, ...], observed):
        return pm.NegativeBinomial(
            name, mu=pt.exp(eta), alpha=concentration, dims=dims, observed=observed
        )

    pmf.predict(h, likelihood, log_rate, expected_observation=pt.exp(log_rate))


pm.model_to_graphviz(
    pmf.build_model(demand_model_pymc, data_pymc, covariates_train_pymc)
)

We fit the model with NUTS through pmf.fit_mcmc, with the same number of chains and draws as the NumPyro fit. Three keyword arguments go to pm.sample through sample_kwargs. We raise the target acceptance rate to \(0.95\), since the PyMC sampler reports tens of divergent transitions on this model at the \(0.9\) of the NumPyro fit. We start the chains in fresh processes (mp_ctx="forkserver") because the JAX runtime of the main text is already active in this kernel and must not be forked. We also turn off the convergence messages of the sampler (compute_convergence_checks=False), because we read the summary table and the divergence count below.

%%time

fit_pymc = pmf.fit_mcmc(
    demand_model_pymc,
    data_pymc,
    covariates_train_pymc,
    draws=1_000,
    tune=1_000,
    chains=4,
    random_seed=rng,
    sample_kwargs={
        "target_accept": 0.95,
        "mp_ctx": "forkserver",
        "compute_convergence_checks": False,
    },
)
idata_pymc = fit_pymc.idata

CPU times: user 1.25 s, sys: 254 ms, total: 1.5 s
Wall time: 12.1 s

Let’s check the sampler diagnostics of the global parameters and the number of divergent transitions.

az.summary(idata_pymc, var_names=global_params)

int(idata_pymc["sample_stats"]["diverging"].sum())

3

The \(\hat{R}\) values are \(1.00\) and the effective sample sizes are above \(1{,}000\) for every parameter, and the sampler reports \(3\) divergent transitions in \(4{,}000\) draws. The non-centered random walk has \(364\) standardized increments, and at target_accept=0.95 the sampler reports a handful of divergences, between \(1\) and \(4\) across the seeds we tried. The posterior means agree with the NumPyro fit of the Model Fitting section to two significant figures on every parameter, and the standard deviations are within a few percent of each other, for example an elasticity of \(-1.50 \pm 0.09\) in both.

The function pmf.forecast rebuilds the model on the full covariates and draws the horizon from every posterior draw. We flatten the draws into the scenario matrix, shuffle the rows and split them as in the Demand Scenarios section.

forecast_pymc = pmf.forecast(
    demand_model_pymc, idata_pymc, data_pymc, covariates_pymc, random_seed=rng
)
scenarios_pymc = (
    forecast_pymc["predictions"]["forecast"]
    .stack(sample=("chain", "draw"))
    .transpose("sample", FUTURE_DIM)
    .values.astype("float64")
)
scenarios_pymc = scenarios_pymc[rng.permutation(scenarios_pymc.shape[0])]
demand_opt_pymc, demand_eval_pymc = (
    scenarios_pymc[:N_SCENARIOS_OPT],
    scenarios_pymc[N_SCENARIOS_OPT:],
)
mean_demand_opt_pymc = demand_opt_pymc.mean(axis=0)

Next, we write the simulator of the Inventory Simulator section as a PyTensor graph. The function simulation_graph is a pytensor.scan over the days with the stock of every scenario as the carry, so the loop is vectorized over scenarios by construction. We compile it once with pytensor.function, with the demand matrix as an input, and wrap the three outputs in the Simulation tuple of the main text as JAX arrays, so that summarize scores the PyTensor plans with the same code as the JAX plans.

def simulation_graph(
    orders: pt.TensorVariable, demand: pt.TensorVariable, config: InventoryConfig
) -> tuple[pt.TensorVariable, pt.TensorVariable, pt.TensorVariable]:
    """The inventory dynamics of one plan against every scenario, as a graph."""

    def step(
        order_t: pt.TensorVariable,
        demand_t: pt.TensorVariable,
        inventory: pt.TensorVariable,
    ) -> tuple:
        available = inventory + order_t  # y_t
        sales_t = pt.minimum(available, demand_t)
        lost_t = demand_t - sales_t
        inventory = available - sales_t
        cost_t = config.holding_cost * inventory + config.lost_sales_cost * lost_t
        return inventory, lost_t, cost_t

    initial_inventory = pt.full(
        (demand.shape[0],), config.initial_inventory, dtype="float64"
    )
    inventory_path, lost_path, cost_path = pytensor.scan(
        step,
        sequences=[orders, demand.T],
        outputs_info=[initial_inventory, None, None],
        return_updates=False,
    )
    return cost_path.sum(axis=0), inventory_path.T, lost_path.T


orders_sym = pt.dvector("orders")
demand_sym = pt.dmatrix("demand")
simulate_fn = pytensor.function(
    [orders_sym, demand_sym], list(simulation_graph(orders_sym, demand_sym, config))
)


def simulate_scenarios_pt(orders: np.ndarray, demand: np.ndarray) -> Simulation:
    """Run the compiled simulator of one plan against a matrix of demand paths."""
    loss, inventory, lost_sales = simulate_fn(
        np.asarray(orders, dtype="float64"), np.asarray(demand, dtype="float64")
    )
    return Simulation(
        jnp.asarray(loss), jnp.asarray(inventory), jnp.asarray(lost_sales)
    )

The forward pass of the Decision Rules section is a second pytensor.scan, over the days, with an order rule applied to the shortfall of every scenario. The quantile rule reads the ceil(prob * N)-th smallest shortfall off pt.sort, which is the inverted_cdf quantile of the JAX version.

def follow_path_pt(shortfall: pt.TensorVariable) -> pt.TensorVariable:
    """Order the shortfall of a single demand path, here the mean path."""
    return shortfall[0]


def quantile_shortfall_pt(prob: float, num_scenarios: int) -> OrderRule:
    """The rule that orders the `prob` quantile of the shortfall over the scenarios."""
    index = math.ceil(prob * num_scenarios) - 1

    def rule(shortfall: pt.TensorVariable) -> pt.TensorVariable:
        return pt.sort(shortfall)[index]

    return rule


def forward_pass_pt(
    demand: np.ndarray, order_rule: OrderRule, config: InventoryConfig
) -> np.ndarray:
    """Daily orders from an order rule applied to the shortfall of every day."""
    demand_sym = pt.dmatrix("demand")

    def step(demand_t: pt.TensorVariable, inventory: pt.TensorVariable) -> tuple:
        shortfall = demand_t - inventory
        order_t = pt.maximum(order_rule(shortfall), 0.0)
        inventory = pt.maximum(inventory + order_t - demand_t, 0.0)
        return inventory, order_t

    initial_inventory = pt.full(
        (demand_sym.shape[0],), config.initial_inventory, dtype="float64"
    )
    _, orders = pytensor.scan(
        step,
        sequences=[demand_sym.T],
        outputs_info=[initial_inventory, None],
        return_updates=False,
    )
    return pytensor.function([demand_sym], orders)(np.asarray(demand, dtype="float64"))

The objective of the Decision Plans section is the same graph with the mean of the losses plus \(\lambda\) times the mean of the \(k\) largest ones, which pt.sort picks out. The solver is projected Adam from pytensor-ml: the orders are a trainable shared variable, adam with a cosine_schedule returns the updates that move them, the transform project_nonnegative clips the updated orders at zero, and compile_train compiles the loss, its gradient and the updates into one training step. The function build_solver resets the orders and the optimizer state before each run and loops over the \(1{,}500\) steps.

def objective_pt(
    orders: pt.TensorVariable,
    demand: pt.TensorVariable,
    lam: pt.TensorVariable,
    config: InventoryConfig,
    k: int,
) -> pt.TensorVariable:
    """Expected loss plus lam times the CVaR of the loss over the scenarios."""
    losses, _, _ = simulation_graph(orders, demand, config)
    return losses.mean() + lam * pt.sort(losses)[-k:].mean()


lam_sym = pt.dscalar("lam")
objective_fn_pt = pytensor.function(
    [orders_sym, demand_sym, lam_sym],
    objective_pt(orders_sym, demand_sym, lam_sym, config, tail_size(N_SCENARIOS_OPT)),
)


def project_nonnegative(loss_gradients_or_updates, parameters: list) -> dict:
    """Clip the updated parameters at zero, the projection onto a >= 0."""
    updates = to_updates(loss_gradients_or_updates, parameters)
    return updates.replacing({p: pt.maximum(updates[p], 0.0) for p in parameters})


def build_solver(
    config: InventoryConfig, k: int, n_steps: int = 1_500, learning_rate: float = 1.0
) -> Callable[[np.ndarray, np.ndarray, float], tuple[np.ndarray, np.ndarray]]:
    """Compile projected Adam on the expected cost plus lam times CVaR."""
    orders = trainable(np.zeros(setup.horizon), name="orders")
    demand = pt.dmatrix("demand")
    lam = pt.dscalar("lam")
    loss = objective_pt(orders, demand, lam, config, k)
    schedule = cosine_schedule(learning_rate, total_steps=n_steps)
    updates = chain(adam(learning_rate=schedule), project_nonnegative)(loss, [orders])
    # the updates carry the Adam moments and the step counter of the schedule
    initial_state = {
        var: var.get_value().copy() for var in updates if var is not orders
    }
    step = compile_train(loss, updates, parameters=[orders], inputs=[demand, lam])

    def solve(
        orders_init: np.ndarray, demand: np.ndarray, lam: float
    ) -> tuple[np.ndarray, np.ndarray]:
        orders.set_value(np.asarray(orders_init, dtype="float64"))
        for var, value in initial_state.items():
            var.set_value(value)
        demand = np.asarray(demand, dtype="float64")
        values = np.array([step(demand, lam) for _ in range(n_steps)])
        return orders.get_value(), values

    return solve


solve_pt = build_solver(config, tail_size(N_SCENARIOS_OPT))

We build the four plans and the myopic newsvendor reference as in the Decision Plans section, and we check the two solves against the linear program solve_lp of the previous appendix.

orders_mean_pymc = forward_pass_pt(
    mean_demand_opt_pymc[None, :], follow_path_pt, config
)
orders_cc_pymc = forward_pass_pt(
    demand_opt_pymc, quantile_shortfall_pt(1 - ALPHA, N_SCENARIOS_OPT), config
)
orders_saa_pymc, _ = solve_pt(orders_mean_pymc, demand_opt_pymc, 0.0)
orders_cvar_pymc, _ = solve_pt(orders_mean_pymc, demand_opt_pymc, 1.0)
orders_newsvendor_pymc = forward_pass_pt(
    demand_opt_pymc, quantile_shortfall_pt(critical_fractile, N_SCENARIOS_OPT), config
)

plans_pymc = [
    Plan("Mean forecast", orders_mean_pymc),
    Plan("Chance constrained", orders_cc_pymc),
    Plan("Expected cost", orders_saa_pymc),
    Plan("CVaR", orders_cvar_pymc),
]

gaps_pymc = []
for orders_pt, lam in [(orders_saa_pymc, 0.0), (orders_cvar_pymc, 1.0)]:
    orders_lp = np.asarray(solve_lp(jnp.asarray(demand_opt_pymc), lam, config))
    value_pt = float(objective_fn_pt(orders_pt, demand_opt_pymc, lam))
    value_lp = float(objective_fn_pt(orders_lp, demand_opt_pymc, lam))
    gaps_pymc.append(value_pt / value_lp - 1)

assert max(abs(gap) for gap in gaps_pymc) < 1e-5
gaps_pymc

[3.5033775880322082e-09, 2.931576685849535e-07]

Finally, we score every plan on the held-out scenarios with the summarize function of the main text and draw the held-out scores figure of the Out of Sample Evaluation section.

demand_eval_pymc_jax = jnp.asarray(demand_eval_pymc)
scores_pymc = {
    plan.name: simulate_scenarios_pt(plan.orders, demand_eval_pymc)
    for plan in plans_pymc
}
results_pymc = pl.DataFrame(
    [
        {"plan": name, **summarize(sim, demand_eval_pymc_jax)}
        for name, sim in scores_pymc.items()
    ]
)
newsvendor_scores_pymc = summarize(
    simulate_scenarios_pt(orders_newsvendor_pymc, demand_eval_pymc),
    demand_eval_pymc_jax,
)

fig, axes = plt.subplots(
    nrows=2, ncols=2, figsize=(12, 8), sharex=True, layout="constrained"
)

for ax, metric in zip(axes.flat, metrics, strict=True):
    sns.barplot(
        data=results_pymc,
        x="plan",
        y=metric,
        hue="plan",
        hue_order=plan_names,
        palette=PLAN_COLORS,
        legend=False,
        ax=ax,
    )
    reference = ax.axhline(
        newsvendor_scores_pymc[metric],
        color="black",
        linestyle="--",
        label="Myopic newsvendor plan",
    )
    if metric in PROBABILITY_METRICS:
        ax.yaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))
    ax.set(title=metric, xlabel="", ylabel="")

for ax in axes[-1]:
    ax.set_xticks(range(len(plans_pymc)), plan_labels)

fig.legend(handles=[reference], loc="outside lower center")
fig.suptitle(
    "Held-Out Scores by Plan (PyMC and PyTensor)", fontsize=18, fontweight="bold"
);

The PyMC and PyTensor pipeline reproduces the figure of the main text. Here is the comparison, plan by plan, with the main text numbers in parentheses.

• Mean forecast. Expected loss \(170\) (\(172\)), tail loss \(508\) (\(508\)), stockout probability per day \(24.6\%\) (\(24.8\%\)) and in the horizon \(84\%\) (\(82\%\)).

• Chance constrained. Expected loss \(152\) (\(156\)), tail loss \(321\) (\(330\)), stockout probability per day \(5.0\%\) (\(4.7\%\)) and in the horizon \(27\%\) (\(25\%\)).

• Expected cost. Expected loss \(140\) (\(142\)), tail loss \(364\) (\(370\)), stockout probability per day \(10.5\%\) (\(10.4\%\)) and in the horizon \(47\%\) (\(47\%\)).

• CVaR. Expected loss \(148\) (\(153\)), tail loss \(312\) (\(313\)), stockout probability per day \(7.8\%\) (\(7.0\%\)) and in the horizon \(42\%\) (\(39\%\)).

• Myopic newsvendor. Expected loss \(153\) (\(157\)) and stockout probability per day \(4.8\%\) (\(4.4\%\)).

The order of the plans is the same on every metric, the myopic newsvendor line sits next to the chance-constrained bar on every panel, and the PyTensor solver lands on the linear program optimum with relative gaps of \(4 \times 10^{-9}\) and \(3 \times 10^{-7}\). The losses differ by at most \(4\%\) and the stockout probabilities by at most \(3\) percentage points, because the posterior draws and the scenario split differ between the two pipelines, not the model or the decision rules.