Building year-round store × category demand forecasts and production recommendations for 300+ bakeries

Every morning a bakery decides how much of each product to bake before it knows how many customers will come. Bake too little and customers leave empty-handed; bake too much and the rest is thrown away. Tous les Jours stores differ by location, customer base, day-of-week pattern, seasonality, promotions, weather and events, so no chain-wide rule fits a particular store. Baking what sold on the same day last week, plus a margin, repeats last week's noise and ignores what the calendar knows in advance. And the initial request — a forecast for one Christmas-season item — would not help a store decide what to bake on an ordinary Tuesday.

In this post we walk through the system built instead. We re-scoped the problem to store × mid-level category × day, all year, for more than 300 stores; built data pipelines with OCR and RPA; described stores with engineered features and clustering; compared Temporal Fusion Transformer, PatchTST and DeepAR as global forecasting models; turned forecasts into production recommendations; and delivered them through a Streamlit dashboard in an automated loop. Testing showed a 5.7% revenue lift. Where the live model on this page differs from the project — a regularized log-linear model instead of the deep models, and an explicit newsvendor rule — we say so.

Solution overview

The system is a scheduled batch loop, not an online service. Collection jobs refresh store data; a feature step builds the store × category × day panel and store descriptors; global models forecast the next two weeks for every series with an interval; a recommendation step turns forecasts into quantities; a dashboard publishes them, and actual sales flow back into evaluation.

The numbered steps in the diagram:

1. Data pipelines pull each store's sales mix, trends, customer characteristics, payment methods and discounts; OCR code and RPA gather commercial-district data for all stores.
2. Sales are aggregated to store × mid-level category × day; stores are profiled and grouped with clustering; calendar, promotion and event features are engineered.
3. Temporal Fusion Transformer, PatchTST and DeepAR are trained as global models across all store-category series and compared.
4. Forecasts are scored on two-week windows: accuracy of the aggregate sales count and prediction-interval accuracy over the 14 days.
5. Forecasts and intervals become production recommendations per store and category; the live model makes the rule explicit as a newsvendor quantile.
6. A Streamlit dashboard publishes the results; collection, forecasting, dashboard generation, delivery and evaluation run as one automated loop.

Technology stack

Step 1: Re-scope the forecast to the decision stores make

The first discussion was about forecasting a Christmas-season cake. That is a visible problem — one item, one week, a large spike — but it is not the decision a store faces on most days. Every morning each store decides how much of each category to bake, and ordinary days add up to most of the year's revenue. We redefined the prediction unit around that decision.

Mid-level categories sit between two poor choices: single items are sparse and change with the range, and store totals say nothing about what to bake. A category such as bread or cakes is stable enough to forecast and close enough to the oven to act on.

Why re-scope at all? A seasonal forecast is used a few days a year. A production decision is made every day in every store, so improving it is an always-on operational gain — unlike large event-driven sales campaigns, which typically run around twice a month.

Step 2: Collect store and district data with pipelines, OCR and RPA

Sales totals alone do not explain why two stores of similar size sell different things. We identified the source data that does and built pipelines for it: each store's sales mix and trends, customer characteristics, payment methods and discounts. Commercial-district analysis data — what surrounds each store — was gathered for all stores with OCR code I developed and CJ Foodville's RPA.

Sales are aggregated to store × mid-level category × day, and each store gets a profile from its sales mix, payment and discount patterns and district data. Clustering the profiles groups stores that behave alike; the clusters become model inputs and give a store with a short history similar stores to borrow from.

import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

def build_panel(sales: pd.DataFrame) -> pd.DataFrame:
    """Transactions -> one row per store x category x day, zero days kept."""
    return (sales.groupby(["store", "category", "date"])["qty"].sum()
                 .unstack("date", fill_value=0)
                 .stack().rename("qty").reset_index())

def store_profiles(sales: pd.DataFrame, district: pd.DataFrame) -> pd.DataFrame:
    """Sales mix, payment mix and discount share per store, plus district data."""
    mix = pd.crosstab(sales["store"], sales["category"], values=sales["amount"],
                      aggfunc="sum", normalize="index")
    pay = pd.crosstab(sales["store"], sales["payment_type"], normalize="index")
    disc = sales.groupby("store")["discounted"].mean().rename("discount_share")
    return mix.join(pay, rsuffix="_pay").join(disc).join(district).fillna(0)

def cluster_stores(profiles: pd.DataFrame, k: int) -> pd.Series:
    X = StandardScaler().fit_transform(profiles)
    labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(X)
    return pd.Series(labels, index=profiles.index, name="store_cluster")

Simplified. Column names are generic; district holds the commercial-district features per store, and k is chosen by inspection.

Step 3: Engineer features for one model across all stores

Three hundred stores with a handful of categories each is too many series to model one at a time, and too few days per series to learn rare events such as Christmas. A global model learns the shared structure once — what an office district does on a Saturday, what a holiday does to a residential street — and lets each store-category contribute its own level.

The drivers — day of week, seasonality, promotions, weather, events, location — become features, and their interactions matter as much as their main effects. The live model encodes them as:

• Weekly shape by district type (office, residential, station, campus): an office district and a residential street have opposite weekends.
• Holidays and rain by district type; promotions and temperature by category: a holiday empties an office district and fills a residential one.
• Christmas phase by category: cakes multiply in the days before Christmas while sandwiches dip.
• Campus term transitions: for campus stores, whether the forecast day is in a vacation minus the share of vacation days in the level window, so a vacation lull is not carried into the first week of term.
• Recent level: log mean sales over the 28 and 7 days before the forecast cutoff.

The level is computed at the forecast cutoff, not on the forecast day. Training rows are built for horizons of 1 to 14 days, each using only data available at its cutoff; a model trained on yesterday's level would look excellent in training and fail on day 14.

import pandas as pd

def cross(df, by, value=None, prefix=""):
    """One-hot of the combined `by` columns, optionally scaled by `value`."""
    key = df[by].astype(str).agg(":".join, axis=1)
    d = pd.get_dummies(key, prefix=prefix, dtype=float)
    return d if value is None else d.mul(df[value], axis=0)

def design_matrix(df: pd.DataFrame) -> pd.DataFrame:
    """One row per store x category x target day."""
    xmas = cross(df, ["category", "xmas_phase"], prefix="xmas")
    parts = [
        df[["level28", "level7"]],                     # recent level at the cutoff
        cross(df, ["district", "dow"], prefix="dow"),   # weekly shape by district
        cross(df, ["district"], "holiday", "hol"),      # holidays by district
        cross(df, ["district"], "rain", "rain"),        # rain by district
        cross(df, ["category"], "promo", "promo"),      # promotion lift by category
        cross(df, ["category"], "temp", "temp"),        # temperature by category
        xmas.loc[:, ~xmas.columns.str.endswith(":0")],  # Christmas phase by category
        df[["term_shift"]],                             # campus term transitions
    ]
    return pd.concat(parts, axis=1)

Simplified. The level and term_shift columns are computed at the forecast cutoff before this step.

Step 4: Compare TFT, PatchTST and DeepAR on two-week windows

In the project, three deep time-series architectures were compared for the store × category forecasts. Each can be trained as one global model over all series:

• Temporal Fusion Transformer (TFT) separates static covariates (store descriptors), inputs known in advance (calendar, planned promotions) and observed inputs (past sales), weighs them with variable-selection networks, encodes the past with an LSTM encoder–decoder and adds interpretable attention. Trained with a quantile loss, it outputs several quantiles at once.
• PatchTST cuts each series into patches of consecutive days and feeds them as tokens to a Transformer encoder, one channel at a time; a patch carries local shape such as a weekly cycle.
• DeepAR is an autoregressive recurrent network that outputs a distribution for each next day and samples paths, from which any quantile — or the distribution of a two-week total — can be read.

The project reports two measures on a two-week window: the accuracy of the aggregate sales count, and prediction-interval accuracy over the 14 days. The first matters for planning; the second for whether an interval can feed a production rule.

import pandas as pd

def score_window(actual: pd.DataFrame, fc: pd.DataFrame) -> dict:
    """Both frames: one row per series x day of a 14-day window.
    fc holds the quantiles q10, q50 and q90."""
    df = actual.merge(fc, on=["series", "date"])
    tot = df.groupby("series")[["qty", "q50"]].sum()
    total_acc = 1 - (tot["qty"] - tot["q50"]).abs().sum() / tot["qty"].sum()
    inside = df["qty"].between(df["q10"], df["q90"])
    return {"two_week_total_acc": total_acc,
            "days_in_interval": inside.groupby(df["series"]).sum().mean()}

def compare(candidates: dict, panel: pd.DataFrame, cutoffs) -> pd.DataFrame:
    rows = []
    for name, make_model in candidates.items():
        for cutoff in cutoffs:
            end = cutoff + pd.Timedelta(days=14)
            train = panel[panel["date"] < cutoff]
            test = panel[(panel["date"] >= cutoff) & (panel["date"] < end)]
            model = make_model().fit(train)
            fc = model.predict(cutoff, horizon=14)      # q10 / q50 / q90
            rows.append({"model": name, "cutoff": cutoff,
                         **score_window(test, fc)})
    return pd.DataFrame(rows).groupby("model").mean(numeric_only=True)

Simplified. make_model wraps one of the candidate architectures behind a common fit/predict interface (not shown); the interval here runs from the 10th to the 90th percentile.

Why probabilistic models? A production decision needs to know how wrong the forecast may be, not only its centre. Quantile and distributional forecasts provide the interval that a production rule consumes.

Step 5: Reproduce the global model in the browser

The deep models do not run in a web page, so the live model uses a stand-in with the same shape: shared calendar, event and district effects plus a store-category level. It is a 56-term ridge regression on log(1 + sales), solved in closed form on one store in five and applied to all 300.

log(1 + y_(s,c,d)) = β₀ + β₁·level28_(s,c) + β₂·level7_(s,c)
                   + dow[district_s, weekday_d] + hol[district_s]·1(holiday_d) + rain[district_s]·1(rain_d)
                   + promo[c]·1(promo_d) + temp[c]·temp_d + xmas[c, phase_d] + τ·term_shift_(s,d) + ε

ε ~ N(0, σ_c²)            ridge penalty λ on all terms except β₀

point forecast    ŷ = k_c · (exp(μ̂) − 1)                       k_c: smearing factor per category
80% interval      [ exp(μ̂ − 1.28·σ_c) − 1 ,  exp(μ̂ + 1.28·σ_c) − 1 ]

The log scale turns multiplicative effects into additive ones. Back-transforming a log-scale prediction underestimates the mean, so a per-category smearing factor — actual over back-transformed totals on the training rows — rescales the point forecast. The residual standard deviation per category gives the interval and, in the next step, the production quantity.

import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge

def fit_global(X: pd.DataFrame, y: pd.Series, category: pd.Series):
    """One ridge model on log1p(sales) for every store and category."""
    z = np.log1p(y)
    model = Ridge(alpha=1.0).fit(X, z)              # intercept is not penalized
    z_hat = pd.Series(model.predict(X), index=y.index)
    sigma = (z - z_hat).groupby(category).std()     # residual s.d. per category
    smear = y.groupby(category).sum() / np.expm1(z_hat).groupby(category).sum()
    return model, sigma, smear

def forecast(model, sigma, smear, X: pd.DataFrame, category: pd.Series):
    mu = pd.Series(model.predict(X), index=X.index)
    s, k = category.map(sigma), category.map(smear)
    return pd.DataFrame({
        "mu": mu, "sigma": s,
        "point": k * np.expm1(mu),
        "lo80": np.expm1(mu - 1.2816 * s).clip(lower=0),
        "hi80": np.expm1(mu + 1.2816 * s),
    })

Simplified. A Python equivalent of the in-browser fit, which solves the same regularized normal equations in JavaScript.

Step 6: Turn the forecast into a production quantity

A forecast is not a production plan. A store needs a number, and the right number depends on which mistake is worse: an unsold item costs what it took to make, and a customer turned away costs the margin. This is the newsvendor problem, and its answer is a quantile of the demand distribution, not its mean.

choose q to maximize      E[ p · min(D, q) − c · q ]

underage cost   c_u = p − c          margin lost on a sale that could not be made
overage cost    c_o = c              cost of an unsold unit (no salvage)

q* = F_D⁻¹( c_u / (c_u + c_o) ) = F_D⁻¹( (p − c) / p )

log-normal forecast:   q* = exp( μ̂ + σ_c · Φ⁻¹((p − c) / p) ) − 1

When the margin exceeds the unit cost, the critical ratio is above one half and the rule deliberately bakes above the median forecast — by more for categories with wider intervals. In the project, forecasts were delivered with intervals and production recommendations through the dashboard; the live model makes one standard rule explicit, with invented prices and costs.

import numpy as np
from scipy.stats import norm

def production_quantity(mu, sigma, price, cost):
    """Newsvendor quantity for a log-normal forecast on the log1p scale."""
    critical_ratio = (price - cost) / price        # c_u / (c_u + c_o)
    z = norm.ppf(critical_ratio)
    return np.maximum(0, np.round(np.expm1(mu + z * sigma)))

def rule_of_thumb(same_day_last_week):
    return np.round(1.10 * same_day_last_week)

def outcome(demand, produced, price, cost):
    sold = np.minimum(demand, produced)
    return {
        "revenue": (sold * price).sum(),
        "lost_sales": ((demand - sold) * price).sum(),
        "waste_cost": ((produced - sold) * cost).sum(),
        "profit": (sold * price - produced * cost).sum(),
        "sold_out_share": (demand > produced).mean(),
    }

Simplified. mu and sigma come from the model in Step 5; the rule of thumb is the live model's baseline.

The live model scores both rules on every store, category and day of a two-week test window. It forecasts all 14 days from one cutoff, while the rule of thumb always looks back seven days. In the ordinary fortnight the rule of thumb is slightly better on two-week totals; at Christmas, last week's sales say little about Christmas Eve.

Step 7: Deliver through a dashboard and automate the loop

Users saw the results in a Streamlit dashboard: store characteristics, sales rankings, sales mix, payment and discount information, and store- and category-level forecasts.

import pandas as pd
import streamlit as st

@st.cache_data
def load(run: str):
    return (pd.read_parquet(f"runs/{run}/forecasts.parquet"),
            pd.read_parquet(f"runs/{run}/stores.parquet"))

run = st.sidebar.selectbox("Run", list_runs())
fc, stores = load(run)
store = st.sidebar.selectbox("Store", stores.index)

left, right = st.columns(2)
left.subheader("Store profile")
left.dataframe(stores.loc[[store], profile_columns].T)
right.subheader("Sales mix")
right.bar_chart(stores.loc[store, mix_columns])

category = st.selectbox("Category", sorted(fc["category"].unique()))
view = fc[(fc["store"] == store) & (fc["category"] == category)].set_index("date")
st.line_chart(view[["lo80", "point", "hi80"]])
st.dataframe(view[["point", "lo80", "hi80", "recommended"]])

Simplified. list_runs(), the column lists and the file layout are placeholders for the outputs of the scheduled pipeline.

The larger change was operational. Data collection, forecasting, dashboard generation, result delivery and evaluation were automated as one loop, so forecasts are refreshed and scored without manual work. That is what makes the revenue effect continuous: a better production decision applies every day, in every store.

Results

Testing showed a 5.7% revenue lift. In an early two-week aggregate sales-count forecast, accuracy improved from an 80% baseline to about 90%, and prediction-interval accuracy improved from 10/14 to 11.5/14.

The key contribution was shifting the problem from a seasonal item forecast to a regular demand-forecasting and production-recommendation system used in store operations, with the whole loop — data collection, forecasting, dashboard, delivery and evaluation — automated. The outcome was connected to a technology-licensing agreement.

Lessons learned

• Scope to the decision, not the request. The unit that mattered was store × category × day, all year — not the seasonal item that prompted the project.
• Pool across stores. A global model learns calendar, event and district effects once and leaves each series only its level; stores with short histories borrow from similar ones.
• Forecast a distribution, then decide. A production quantity needs an interval, and the right quantile depends on margin and cost — the same forecast gives different quantities for different categories.
• Automation is part of the model. Without automated collection and delivery a daily forecast is a report; with it, it is part of how stores operate.

Conclusion

Bakery production is a daily decision under uncertainty in every store. Re-scoping a seasonal item forecast into store × category × day forecasts for more than 300 stores, feeding it with automated pipelines, comparing probabilistic deep time-series models and delivering recommendations through a dashboard made forecasting part of store operations, with a 5.7% revenue lift in testing.

The pattern — one global probabilistic model, a decision rule that weighs the two kinds of error, and an automated loop that delivers and scores the result — applies to other perishable goods whose demand varies by location and calendar.

Limitations

• Sales are censored by what was baked: a sold-out day hides how much more could have sold, and forecasts trained on sales inherit that bias unless it is modelled.
• New stores and new products have no history; they borrow from similar stores, which works for districts and less well for genuinely new items.
• The live model's stores, sales and prices are generated, its forecasts come from a ridge log-linear model rather than the deep models compared in the project, and its production rule is a textbook newsvendor quantile.

About the demo and confidentiality

Stores, districts, categories, prices, weather and sales in the embedded model are generated. No store, sales, customer or commercial-district data from Tous les Jours or CJ Foodville, and no model configuration or production-rule parameters, appear in this post; code is simplified and written for illustration.