Retail Demand Forecasting with Prophet: From Model Fundamentals to Engineering Practice

1. Introduction: From Retail Demand Forecasting to Prophet #

Retail demand forecasting does more than estimate future daily demand for each store-item combination. It must also identify temporal patterns and business effects in historical sales while accounting for observational limitations. Consider a Melbourne supermarket planning replenishment for the next four weeks: weekends create recurring patterns, promotions and public holidays create event effects, and long-term changes in demand form a trend. Store closures, stockouts, or periods when an item is not on sale may reduce observed sales or make them zero, without implying an equivalent decline in underlying demand.

Meta’s open-source Prophet is a time-series forecasting method with a Python implementation. Its decomposable model combines trend, seasonality, and holiday effects and can incorporate additional business variables [1]. Prophet is better suited to series with clear seasonal patterns and enough history to cover multiple seasonal cycles, so this article does not assume that it is appropriate for every retail series.

This article constructs a synthetic retail dataset, defines each demand series by a store-SKU (stock-keeping unit) combination, uses daily observations and forecasts, predicts demand for the next 28 days, and discusses how to extend the workflow to large-scale, multi-series forecasting. It combines model formulas, parameter estimation, forecast reconstruction, time-series cross-validation, and candidate-model comparisons to examine when Prophet is appropriate and where its limits lie. All data and business parameters are teaching assumptions and do not represent the operations of any actual retailer.

The article first defines the forecast target and develops the teaching case, then explains Prophet’s model structure, data processing, model configuration, MAP parameter estimation, and daily forecast reconstruction. It then evaluates the model through time-series backtesting, error diagnosis, and parameter experiments, before comparing Prophet with intermittent-demand methods, ARIMA, tree-based models, and cross-series approaches. The final sections cover batch training, forecast delivery, monitoring, and rollback. Complete Python code, the accompanying notebook, and setup instructions are available in the GitHub repository: retail-demand-forecast.

2. Demand Forecasting Project: From Forecast Target to Case Design #

Sales vs. Demand: Define the Forecast Target First #

This article defines the following: store_id uniquely identifies a store, product_id identifies a product concept, and sku_id identifies a stock-keeping unit (SKU) with specific specifications, capacity, packaging, or sales unit. The same product_id can correspond to multiple sku_ids. For example, 500g and 1kg packages of a certain brand of pasta can share one product_id but use different sku_ids.

A demand sequence is jointly defined by store_id and sku_id; each daily observation is identified by the business date business_date, so the daily record key is store_id + sku_id + business_date. The product master data structure may differ across enterprises; when connecting to a real system, the actual definitions should be used.

At this granularity, the team must also clarify the forecast target: should the model predict the sales recorded by the system, or the latent demand that customers would have generated under unconstrained supply? Recorded sales are directly observable. Latent demand is closer to what replenishment decisions require, but it usually needs to be estimated during stockouts. The two targets require different data treatment, model interpretation, and business use.

Observed sales are constrained by both latent demand and available supply. For example, if customers want to buy 100 units but only 60 units are available that day, at most 60 units can be sold. Under the simplifying assumptions that the store is operating normally and that customers do not substitute products, switch channels, or delay purchases:

Notation explanation:

These three quantities are non-negative and use the same product, store, and date scope. This formula simplifies the intra-day supply and demand arrival order, so whether a stockout occurred cannot be inferred solely from the total inventory on that day.

Real-world retail also involves product substitution, channel switching, and delayed purchases. Even so, the relationship above illustrates the central issue: when a stockout occurs, observed sales are a censored observation of latent demand. If supply-constrained sales are treated as complete demand during training, the model may underestimate true demand.

Production data needs to distinguish at least the following statuses:

The teaching benchmark masks stockout days in the model target column while retaining the original sales values and stockout flags. This prevents the model from interpreting supply-constrained sales as low demand. The model then fits demand patterns using sales from the remaining valid operating days, but this does not recover latent demand on stockout days; those masked dates still require a dedicated demand-imputation method.

Teaching Case: Generating Sales Data and Defining the Time Split #

This article constructs a simulated daily sales series for a high-frequency food item in an Australian large supermarket as a teaching case. The data covers January 1, 2023, to December 31, 2025, a total of 1096 days; the model uses data up to December 3, 2025, for training and forecasts the last 28 days. All data and business rules are for teaching purposes and do not represent the actual sales situation of any enterprise.

This case aims to demonstrate how trends, seasonality, holidays, discounts, stockouts, and store closures affect sales. It represents high-frequency demand with many non-zero observations and does not cover intermittent demand with large numbers of zero sales; related models will be discussed in Chapter 7.

Assuming normal store operation, the noise-free mean demand for the day is obtained by summing contributions from trend, weekend, annual seasonality, holidays, and discounts, and adding random noise:

Here, \(\mu_t^{\mathrm{open}}\) represents the noise-free mean demand assuming normal store operation; \(t\) is the number of days since January 1, 2023, and \(W_t\), \(H_t\), and \(q_t\) represent weekend, public-holiday, and discount variables, respectively. \(D_t^{\mathrm{sim}}\) is the simulated latent demand after random noise is added.

\(\epsilon_t\) is independent and identically distributed random noise, normally distributed across dates with a mean of 0 and variance of \(6^2\), thus a standard deviation of 6 units. Superscripts open and sim denote the scenario of assumed normal operation and simulated quantity, respectively.

After latent demand is generated, operating and supply constraints are applied to produce the sales records. On operating days, sales equal the smaller of demand and available supply, rounded to the nearest integer; on closed days, sales are recorded as zero:

The effective supply quantity uses the following simplified rules:

Here, \(O_t\) is the 0/1 operating status, \(V_t\) is the effective supply quantity, and \(y_t^{\mathrm{sales}}\) is the observed sales volume after supply constraint and rounding to units. The code sets closed-day sales to zero via \(O_t\). The fixed 20-unit cap is only used for generating teaching data; this example does not simulate continuously changing inventory trajectories.

The discount \(q_t\) is 0.20 for the first 7 days of every 14-day period, and 0 otherwise; operating days are randomly flagged for stockout with a 3% probability. The teaching store is closed on Good Friday and Christmas, and other dates are treated as full-day operations. In actual use, transaction restrictions need to be verified based on store category and applicable licenses, and a store operating calendar [2] should be provided.

The generated CSV contains date ds, raw observed sales y, operating flag is_open, closed holiday flag is_closed_holiday, discount discount, and stockout flag stockout. These fields can usually be recorded or predetermined in real business systems. The file also retains conditional_mean and latent_demand, used for verifying the synthetic data generation mechanism and diagnosing model errors; these two columns are not used as model input, and are typically not directly observable in real business.

Raw observed sales form the basis for constructing the model target column. Before training, target values corresponding to stockout days and closed days will be set to missing, but the original CSV still retains sales records. Potential demand on closed days represents a “hypothetical operation” counterfactual; this example does not simulate reduced operating hours for ANZAC Day.

The figure below first shows the complete three-year sales, then magnifies a 21-day window including discounts, stockouts, and store closures.

The table below lists daily samples within the same window, showing only the six real business fields mentioned above. On March 30, a stockout occurred, and observed sales were constrained by the 20-unit supply cap; on April 7, the store was closed for Good Friday, and observed sales were recorded as zero.

This set of smooth cycles and stable additive effects makes it easy to compare the generation formula and fitted components item by item. After changing seeds, dates, or formulas, the model needs to be re-run, and component values and evaluation results in the main text updated. For the complete generation process and field descriptions, please refer to the data generation Notebook.

3. Prophet’s Modeling Framework: Representing Temporal and Business Patterns #

This chapter explains how Prophet models time patterns and business influences in demand series, starting with its model structure. Prophet decomposes the demand series into trend, seasonality, holiday effects, and additional business variables, establishing a mathematical form for each component. Readers can refer to the data generation rules in Chapter 2 to understand how these patterns are incorporated into the model. Chapters 4 and 5 will detail input transformation and parameter estimation, respectively.

To avoid confusion between engineering and statistical terminology, this paper first establishes unified conventions. In data tables, y is referred to as the target variable in machine learning and the response variable in statistical modeling. A row of data used for fitting is called a training sample, and in likelihood derivation, it’s an observation. Model input columns are generally called features; those participating in regression calculations are also known as regressors or covariates. Engineering workflows can train models in batches, while a single Prophet model estimates parameters through fitting. In this article, ’label’ refers only to observed target values or their availability status in the data pipeline, and should not be confused with unobservable latent demand.

Understanding Forecasts with the Additive Model #

This paper uses \(t\) to denote the number of days from a fixed start date, \(y(t)\) for the observed target variable on day \(t\), and \(\hat y(t)\) for the point forecast given by the model. For clarity, the formulas in this chapter use business scales such as days and sales units. Prophet internally scales time and target values during fitting, so scales must be unified before comparing formula parameters with internal model parameters. Each symbol is defined upon its first appearance and used consistently thereafter; internally scaled quantities are denoted with a tilde above the symbol.

Prophet’s additive model expresses the daily observed value as the sum of several components: a long-term trend, recurring seasonality, holiday and business event effects, and random error not yet explained by the model. Each component adopts the same additive scale as the target variable:

Variable definitions:

If additional business features like discount rates or display status are provided, the model estimates corresponding coefficients for each feature and adds the contribution of “feature value multiplied by coefficient” to the forecast. These inputs are called extra regressors, and their dot product with the coefficients forms the regression component:

New regression term notations are as follows, with others consistent with the above formula:

In additive mode, the target variable, trend, seasonality, event effects, and regression components all use the same units, e.g., items/day. During prediction, the model sums the fitted trend, seasonality, event effects, and regression components to obtain the point forecast \(\hat y(t)\); the random error term \(\varepsilon_t\) describes unexplained fluctuations of observed values around these systematic components. Prediction intervals also incorporate uncertainties such as observational noise and future trend changes; their meaning will be further discussed in Chapter 6.

Thus, Prophet more closely resembles a regression model with a time structure. Given dates and business inputs for the next 28 days, the model can directly compute daily forecasts; the default model does not automatically leverage temporal dependencies between yesterday’s sales or residuals; such structures require separate modeling. Chapter 5 will use the Boxing Day forecast of approximately 129.62 items as an example to verify how each component sums up to the final result.

Trend: Representing Growth, Changepoints, and Saturation #

Linear Trend: Fixed Daily Change #

A linear trend model assumes that demand changes by a fixed quantity each day, increasing or decreasing from a base level. Multiplying the daily change by the number of elapsed days and adding the starting level yields the linear trend:

Variable descriptions:

According to the data generation rule \(85+0.015t\) from Chapter 2, if short-term factors like weekends, seasonality, holidays, and discounts are not considered, the trend baseline for this product on January 1, 2023, is 85 units per day, subsequently increasing by 0.015 units daily, which equates to an annual increase of approximately 5.48 units. This slow upward trend can be interpreted as a composite reflection of factors such as community population growth, store maturation, or long-term foot traffic increase. It should be noted that these are merely business interpretations for a didactic scenario and do not imply that the model can identify specific causes.

The trend value is not the final observed sales. Actual sales are also influenced by factors such as seasonality, holidays, discounts, random fluctuations, and supply constraints. The fitted trend also absorbs the average levels of other modules, so its parameters do not necessarily have to be identical item-by-item to the generated values. In the didactic model, the trend value on December 4 is approximately 105.50 units, accumulating an increase of about 0.46 units by the end of the month. The method of allocating weekend increments will be explained in the next section.

Changepoints: Changing Slope While Maintaining Continuity #

A linear trend assumes a constant growth rate, but in actual business, long-term changes rarely conform perfectly to this assumption. The Prophet model therefore introduces a set of candidate changepoints, allowing the trend to adjust its slope at these specific times.

For each candidate changepoint, the model constructs an indicator variable: it takes 0 when the date is before the changepoint, and 1 on or after the changepoint. Let the positions of the \(J\) candidate changepoints be \(c_1,\ldots,c_J\). The corresponding indicator variable is:

Notation descriptions:

When time reaches the \(j\)-th changepoint, \(a_j(t)\) changes from 0 to 1, and the trend slope increases by \(\delta_j\). The model simultaneously introduces an intercept adjustment \(\gamma_j\) to prevent abrupt changes in the trend value when the slope changes. This yields the expression for a continuous piecewise linear trend:

Notation descriptions:

\(\gamma_j\) is not an independently estimated free parameter but is determined by the continuity constraint \(\gamma_j=-c_j\delta_j\). At \(t=c_j\), the new slope contribution \(c_j\delta_j\) exactly cancels out the intercept adjustment, thereby maintaining the continuity of the trend value, with only the subsequent slope changing. This relationship applies to piecewise linear trends, while piecewise logistic trends use different continuity rules.

For example, assume that when time reaches \(t=100\), the trend slope increases from 0.02 units/day to 0.05 units/day. The slope change \(\delta=0.03\) units/day, and the intercept adjustment \(\gamma=-100\times0.03=-3\) units. At the changepoint, the new slope contribution of \(0.03\times100=3\) units exactly cancels out the -3 units intercept adjustment, so the trend value does not jump. Thereafter, the trend value increases by an additional 0.03 units per day compared to the extended value under the original slope. This example only serves to illustrate the principle of continuity; the didactic case does not feature such trend shifts.

The positions of candidate changepoints and their slope changes are determined in two different stages. In Prophet version 1.4.0, as verified in this document, if the user does not manually provide changepoint dates, the default configuration generates 25 candidate changepoints. Since the training historical data in this document is sufficiently long, the actual number of candidate changepoints is \(J=25\). The model first excludes records where the target value y is missing, then approximately equally spaced candidate positions are set from the first 80% of valid observations in the training history. For complete daily series, these positions are also approximately equally spaced on the calendar. The candidate dates are then converted into internal time coordinates \(c_1,\ldots,c_J\); these positions are determined before parameter optimization and are not searched for by the optimizer among all dates [3, 4].

During the fitting process, Prophet jointly estimates the initial slope \(k\), intercept \(m\), the slope changes \(\delta_1,\ldots,\delta_J\) for each candidate position, as well as seasonality, holiday effects, and additional regression coefficients, all within the same objective function. In this document, \(J=25\), so model.params["delta"] contains 25 corresponding slope changes. A sparsity-inducing prior encourages most \(\delta_j\) values to shrink close to zero, retaining larger changes only at positions where the data provides sufficient support; a larger changepoint_prior_scale value generally results in greater trend flexibility. The intercept adjustment \(\gamma_j\) is then determined based on the continuity constraint of \(c_j\) and \(\delta_j\).

Although the didactic example does not manually set trend shifts, the model still uses the aforementioned candidate positions and shrinkage mechanism. After model fitting, the slope changes corresponding to each candidate position can be inspected via model.params["delta"]. Candidate positions merely indicate where the model allows the slope to change; only when the estimated \(\delta_j\) deviates significantly from zero does it imply that the fitted result retains a significant slope change at that position.

Trend changes may be related to factors such as foot traffic, product display, or store adjustments. The model can only estimate the location and magnitude of the change; the specific causes still require support from business records.

Logistic Trend: Gradual Slowdown as Saturation is Approached #

Take a newly opened community supermarket as an example: after a high-frequency food item is stocked, as nearby customers gradually become aware of the store and form repeat purchases, daily demand may show rapid growth in the early stages. Once the potential customer base and purchase frequency within the catchment area stabilize, long-term demand may gradually approach 200 units per day, rather than continuing to grow at a fixed rate.

If both historical data and business judgment support such a process of growth gradually stabilizing, then a constant slope or piecewise linear trend may no longer be appropriate. Therefore, the Prophet model also supports logistic growth. When upper and lower bounds are fixed and the growth rate is positive, the trend will rise along an S-shaped curve, gradually slowing down as it approaches the upper limit. Without considering changepoints, its expression is:

Notation descriptions:

When \(C(t)=C\) and \(F(t)=F\) are fixed values and \(\rho>0\), the logistic trend has the following basic characteristics:

• When \(t\) is much less than \(t_0\), \(g(t)\) approaches the lower bound \(F\).
• When \(t=t_0\), \(g(t_0)=(F+C)/2\), and the trend value is at the midpoint between the upper and lower bounds.
• When \(t\) is much greater than \(t_0\), \(g(t)\) approaches the upper bound \(C\).
• \(\rho\) controls the speed at which the curve rises around \(t_0\); a larger value of \(\rho\) indicates a more concentrated change and a steeper curve.

The figure below takes a fixed lower bound \(F=0\) and a fixed upper bound \(C=200\) as an example. The trend reaches the midpoint of 100 units at \(t=t_0\) and grows fastest around this position; thereafter, the growth rate gradually slows down, and the trend continuously approaches but does not cross the upper bound.

The didactic example in this document uses a linear trend; the logistic settings in the figure are only for extended experimental illustration. The upper limit of 200 is solely for explaining the function’s properties and does not represent the actual demand ceiling for the product. Upper and lower bounds in production models need to have business justification and be validated through historical data and time-series cross-validation.

The use of \(\rho\) and \(t_0\) here aims to avoid confusion with the slope \(k\) and intercept \(m\) in linear trends. When upper and lower bounds are fixed, \(t_0\) is also the inflection point of the curve; however, when cap or floor vary over time, the aforementioned symmetry, asymptotic properties, and inflection point position are no longer fully applicable and require interpretation in conjunction with the boundary functions. cap and floor should reflect saturation levels with business justification [5].

Shelf inventory determines how much demand can be met on a given day, while cap describes the saturation upper limit that the trend component can approach. The two have different meanings and must be set separately. Directly using today’s inventory as cap would introduce existing supply constraints into demand forecasting. Furthermore, it should be noted that logistic bounds only constrain the trend component; after superimposing seasonality, holiday effects, and additional regression components, the final forecast may still be higher than cap or lower than floor.

In actual projects, newly opened stores, products entering a new business district for the first time, or new products in their introductory phase may all experience a process from growth to stabilization. AI/ML engineers should combine sequence length, business phase, and market capacity to determine whether both linear and logistic trends should be included as candidate models, and then compare them through time-series cross-validation, rather than choosing solely based on the observed curve shape.

When using a logistic trend, historical data must provide cap for each date, in addition to ds and y; future forecast data must also provide cap using the same methodology. If a non-zero lower bound is required, floor should also be provided for both historical and forecast intervals, ensuring that cap > floor for every date. Bounds can vary over time, but their future values must be known at the time of forecasting or be determinable in advance based on business rules. During model initialization, Prophet(growth="logistic") must also be explicitly set; otherwise, Prophet defaults to a linear trend [5].

In actual projects, cap should not be set solely based on historical maximum sales, but rather by comprehensively considering the performance of similar stores or products in their mature phase, the potential customer base in the business district, product penetration, and business plans, and then testing the sensitivity of the results to boundary settings through cross-validation. For new stores or products lacking their own sales history, it is necessary to leverage similar series, hierarchical information, or cold-start rules; enabling a logistic trend itself does not solve the problem of insufficient data.

Seasonality: Describing Repetitive Patterns with Fourier Series #

Sales exhibit seasonality over a week or a year, manifesting as repeating highs and lows tied to specific dates. Prophet uses Fourier series to describe this periodic structure.

From a time-domain perspective, seasonality is a curve that repeats over time; from a frequency perspective, this curve can be decomposed into multiple sets of sine and cosine waves of different frequencies. Low-frequency waveforms capture broad, gradual changes within a period, while high-frequency waveforms supplement local fluctuations. By superimposing these periodic basis functions with different weights, one can express or approximate smooth weekly and yearly patterns.

Prophet does not extract periods from the frequency domain using Fast Fourier Transform (FFT). Instead, the model constructs Fourier basis functions based on dates and predefined period lengths. During fitting, the coefficients of these basis functions are jointly estimated along with parameters for trend, holidays, and extra regressors. For a seasonality component with period length \(P\), it can be written as:

Notation Explanation:

In this document, time is measured in days, so \(P\) is also calculated in days. When multiple seasonalities are present simultaneously, \(s(t)\) in the base model is the sum of the individual \(s_P(t)\) components.

In Prophet, the period length \(P\) in the formula corresponds to period, and the Fourier order \(N_P\) corresponds to fourier_order. This document sets weekly seasonality to \(P=7\) and \(N_P=3\), and yearly seasonality to \(P=365.25\) and \(N_P=5\). These settings define the periodic functions the model can use but do not directly specify the exact magnitude of sales increases on weekends; individual Fourier coefficients still need to be estimated from training data.

To understand these settings, it’s necessary to distinguish between data sampling frequency, period length, and curve complexity:

How Do Dates Form Features, and When Do Sales Participate? #

For intuition, let Monday be \(t=0\), Tuesday \(t=1\), and so on through Sunday at \(t=6\). Substituting these values into the sine and cosine functions yields the corresponding weekly Fourier features. For the first harmonic:

Monday and the following Monday obtain identical feature values, indicating that this set of functions repeats every seven days. The calculation method for the second and third orders is the same; simply replace \(n\) in the formula with 2 and 3, respectively.

Dates determine the basis functions, and sales are used to estimate the coefficients of these basis functions. Given \(P\) and \(N_P\), the model directly constructs sine and cosine features based on dates, a process that does not involve sales. For a third-order weekly seasonality, each date will generate six determined feature values. If the training data includes 1,000 valid dates, a feature matrix of 1,000 rows and 6 columns will be formed: each row corresponds to a date, and each column corresponds to a basis function. \(u_{P,n}\) and \(v_{P,n}\) are the parameters to be estimated, which, along with the trend, holiday, extra regressor, and noise parameters, are jointly estimated based on sales during the training period.

Overall, this joint estimation problem can be conceptually stated as:

Here, \(\mathbf X\) contains date and business features determined before fitting, \(\mathbf y\) represents observed sales during the training period, and \(\theta\) summarizes all parameters to be estimated. The seasonality coefficients \(u_{P,n}\) and \(v_{P,n}\) are just a part of \(\theta\) and are not solved independently from the trend and other model components. Furthermore, joint estimation does not mean the model automatically identifies business reasons: if promotions are concentrated on weekends, promotional features and weekly seasonality features might jointly explain the same portion of sales variation. Chapter 5 will explain MAP parameter estimation within this overall framework.

What Do Period Length and Fourier Order Control, Respectively? #

Period length \(P\) determines the repetition cycle of the pattern, while Fourier order \(N_P\) determines the amount of detail the model can express within one period. When the order is \(N_P\), the model will generate \(2N_P\) columns of sine and cosine features and estimate the same number of coefficients.

In the training data, \(W_t\) is 1 on Saturdays and Sundays and 0 on other days, so \(16W_t\) indicates that weekend demand is 16 units higher than on weekdays. This pattern resembles a step wave that repeats every seven days. For daily data, there are only seven observation points in a week; after removing the weekly average, six independently variable directions remain. A third-order Fourier series precisely generates six feature columns, thus sufficient to express this weekly pattern across seven discrete positions. The model’s baseline level absorbs the weekly average contribution of \(16W_t\), while the weekly seasonality component expresses the offset of each day of the week relative to the weekly average.

Here, “expressive capability” specifically refers to discrete time points with one observation per day. If time were treated as a continuous variable, the Fourier functions would still form smooth curves between adjacent dates, not continuous square waves with vertical edges.

The yearly periodic term itself is a sine wave, containing only a first-order harmonic with a period of 365.25 days, so \(N_P=1\) already provides sufficient expressive capability. This document’s initial setting of \(N_P=5\) means the model could also use second to fifth-order harmonics to express more subtle intra-year variations; these additional waveforms are not necessary to restore the training data generation rule and might increase the risk of fitting noise. Internal validation ultimately chose a yearly order of 1, consistent with the known data generation mechanism.

Period length and Fourier order define the model’s expressive space, but they do not guarantee that the fitted results will perfectly reproduce the generating formula. In actual training, seasonality coefficients are jointly estimated with trend, holiday, and discount coefficients and are subject to random noise, missing target values, and prior shrinkage. Real-world projects should use time-series cross-validation to select \(N_P\) and seasonality_prior_scale: the former controls the curve’s detailed expressive capability, and the latter controls the shrinkage strength of seasonality coefficients. Reliably estimating yearly seasonality also requires historical data covering a sufficient number of complete years and business changes.

How Should Seasonality Coefficients Be Interpreted? #

The sine and cosine coefficients of the same order need to be interpreted together: they jointly determine the amplitude and phase of that order’s oscillation. Changing the starting point of time might alter individual coefficients but won’t necessarily change the final curve. Therefore, business interpretations should focus on the seasonal contribution generated on specific dates after all orders are superimposed, rather than isolating and interpreting a single Fourier coefficient. Each order describes a mathematical frequency and doesn’t necessarily correspond to an independent business mechanism; factors such as weekend foot traffic, staffing schedules, promotions, and operating hours might all be reflected in the same seasonal curve.

In additive mode, the seasonal contribution shares the same unit as the target variable y. Daily-level data can only identify periodic changes between dates, not hourly patterns within a day, so this document does not model intra-day seasonality. Chapter 5 will explain how to extract coefficients and reconstruct daily contributions.

Holidays and Promotions: Turning Event Calendars and Promotion Plans into Model Inputs #

Holidays and promotions are business inputs beyond trend and seasonality, but their data formats differ: holidays are usually represented as discrete event features, while discounts can be represented as continuous regression variables. Similar to seasonality modeling, the model first constructs features based on the given business information and then estimates the corresponding coefficients using sales from the training period.

Holidays: Constructing Indicator Features Using Event Calendars #

Holidays and business events usually occur on predetermined dates. Prophet generates indicator features for each date based on the event calendar and then estimates the coefficients corresponding to these features. When an event is not active on date \(t\), its feature value and contribution are both zero; when multiple events are active simultaneously, their individual contributions are summed:

Notation Explanation:

Event names, dates, and windows are configurations given before fitting; \(z_\ell(t)\) is determined by these, without needing sales data for calculation. The coefficients \(\kappa_\ell\) for each event feature are then jointly estimated by model.fit() along with trend, seasonality, and other parameters.

Event impacts do not necessarily occur only on the day of the event. Pre-holiday purchases, post-holiday declines, or multi-day activities can be expressed using event windows [6]. Each relative date within the window can correspond to an independent feature and coefficient, for example, “day before holiday” and “holiday itself”; the model does not pre-assume that these dates have the same impact.

For example, if a store is known to be closed on Christmas Day, customers might shift their purchases of holiday goods to the day before Christmas. When modeling, an independent event feature can be generated for “the day before Christmas,” and its corresponding coefficient can estimate the sales change due to early purchases. The closure status on Christmas Day still needs to be handled separately and cannot be replaced by a holiday coefficient. Whether the coefficient for the previous day is positive and by how much should be estimated from historical observations; the event window merely allows the model to express this early purchase effect, it does not pre-assume its existence.

In the teaching data, the holiday term is written as \(28H_t\), where \(H_t=1\) if date \(t\) is a designated public holiday, and 0 otherwise; thus, assuming normal store operation, public holidays add 28 units to the mean demand. The teaching example sets the holiday window to 0, estimating only the effect on the holiday itself; Prophet, however, estimates coefficients separately by holiday name; for instance, the fitted contribution for Boxing Day is approximately +33.06 units.

Target values corresponding to Good Friday and Christmas are excluded due to store closures. Therefore, these event features have no valid observations available for coefficient estimation. The zero coefficients obtained in such cases reflect insufficient training information and prior shrinkage, and cannot be interpreted as a true demand effect of zero. This article retains these columns for inspection, but actual models can remove event features without valid observations.

A public holiday calendar is not equivalent to a store’s complete business calendar. National holidays can be a starting point, but it’s also necessary to check holidays in the store’s state, mall operating schedules, and the company’s own promotional activities. The model can only learn events that are already defined in the data; missing or erroneous calendars will not be automatically corrected by the model.

Promotions: Expressing Discount Impact Using Continuous Business Variables #

The teaching data uses \(130q_t\) to represent the linear contribution of discounts to demand, where \(q_t\) is the discount rate on date \(t\). A 20% discount sale corresponds to \(q_t=0.20\), so this term increases the noise-free mean demand by 26 units on that day. Prophet can treat discounts as continuous regression variables: the variable value is multiplied by the model’s estimated coefficient to get its contribution to the day’s forecast:

Notation Explanation:

The discount rate \(q_t\), additive or multiplicative mode, standardization, and prior scale need to be determined before fitting; the discount coefficient \(\beta_{\text{discount}}\) is then jointly estimated by model.fit() along with other parameters. This article sets standardize=False, using the original discount rate directly. Therefore, when the discount increases by 0.01, i.e., by one percentage point, the predicted contribution increases by \(0.01\beta_{\text{discount}}\). If standardization is enabled, the coefficient corresponds to changes in the standardized variable and cannot be directly interpreted in terms of original percentage points.

The teaching example fits \(\beta_{\text{discount}}\approx132.62\). When \(q_t=0.20\), the discount contribution is approximately +26.52 units, close to the 26 units in the generative formula. This value is only the individual contribution of the discount variable and does not include trend, seasonality, or holiday effects. The discount in the example appears every 14 days, thus always falling on the same day of the week. Historical data must still include both promotional and non-promotional days for that week for the model to have information to distinguish between discount effects and weekly seasonality.

The teaching data only includes two discount levels: \(q_t=0\) and \(q_t=0.20\). Actual stores, for markdown pricing or to accelerate inventory turnover, might apply 40% or 50% discounts, corresponding to \(q_t=0.40\) or \(q_t=0.50\). Directly plugging these values into the current linear model would yield discount contributions of approximately +53.05 units and +66.31 units, respectively; however, this merely extrapolates the linear relationship estimated from 0% to 20% discounts beyond the training range, and does not imply that historical data supports such a magnitude of increase.

Especially when \(q_t=0.50\), deep discounts can quickly lead to stockouts. In such cases, observed sales can at most reach the available inventory; zero sales after a stockout indicate no items to sell, not a lack of customer demand. Therefore, it is crucial to distinguish between observed sales and potential demand: sales records during stockouts are censored observations constrained by supply and cannot be directly used to estimate the full demand increase brought by deep discounts.

Business inputs used for forecasting must cover the entire forecast horizon. Holidays can be generated in advance based on a calendar, while future discounts require an established promotion plan. If the discount plan is not yet finalized, different promotion scenarios need to be constructed, or the discount variable itself needs to be forecasted, rather than directly using historical values.

The current regression form assumes that each additional percentage point of discount yields the same marginal contribution. Actual promotions may involve activation thresholds, diminishing marginal returns, interactions between display and discount, substitution among similar products, and subsequent declines due to advanced purchases; markdown discounts often coincide with product life cycles and inventory pressure. When these mechanisms need to be expressed, discount bins, piecewise linear terms, nonlinear terms, interaction terms, or lagged features can be constructed, and then their improvement to forecasts can be evaluated through backtesting.

Finally, \(\beta_{\text{discount}}\) describes a conditional association within the model. Retailers might schedule promotions when weak demand is expected, and price, promotional activity, and underlying demand might change simultaneously. Therefore, this coefficient cannot be directly interpreted as the incremental sales caused by the promotion. To estimate the causal effect of promotions or price elasticity, appropriate research designs and identification assumptions are also required.

Additive or Multiplicative: Does the Effect Scale in Units or Proportionally? #

Additive effects describe changes in fixed quantities (e.g., “16 units increase on weekends”), while multiplicative effects describe proportional changes relative to the trend level (e.g., “20% increase over the trend on weekends”). These two effects can coexist: a model can first adjust the trend proportionally, then overlay contributions expressed in terms of units. Separating the forecast mean from the observation error, it can be expressed as:

Variable definitions:

The difference between the two modes amplifies with the trend level. For example, if a periodic effect is 20%: when the trend is 100 units, it contributes 20 units; when the trend rises to 150 units, it contributes 30 units. In contrast, if an additive effect is fixed at 20 units, its contribution remains 20 units whether the trend is 100 or 150 units.

Prophet [7] supports combining additive and multiplicative components within the same model. The specific mode needs to be configured before fitting, and the coefficients for each component are still jointly estimated by model.fit(). The trend \(g(t)\) serves as a baseline for other components and is not classified as either additive or multiplicative.

In the tutorial examples of this article, since \(16W_t\), \(28H_t\), and \(130q_t\) represent increments of fixed quantities, weekly seasonality, holidays, and discount regression terms all adopt the additive mode. The previously calculated unit contributions can be directly summed. The multiplicative mode has not been separately validated in the tutorial examples, so its effect on other retail series cannot be inferred from the current experiments.

Although the interface supports mixed modes, this does not mean that this combination should be adopted by default in real-world projects. When choosing, one should judge whether the amplitude of fluctuations expands synchronously with the trend level, and compare the pros and cons of candidate configurations through backtesting. After changing the mode of any component, the model must be re-fitted, and the contribution of the corresponding component re-interpreted based on units or proportions.

Chapter 3 Review #

The tutorial examples in this article all use an additive mode. For comparison, the overall model presented at the beginning of this chapter is listed again here:

Here, \(g(t)\), \(s(t)\), \(h(t)\), and \(\mathbf{x}_t^\top\boldsymbol\beta\) correspond to trend, seasonality, holidays, and extra business variables, respectively, and \(\varepsilon_t\) represents the observation error not yet explained by the model. If some components adopt a multiplicative mode, the aforementioned combined form \(g(t)[1+E_{\mathrm{mult}}(t)]+E_{\mathrm{add}}(t)\) should be used; the equation above specifically corresponds to the purely additive tutorial examples in this article.

The table below summarizes the business problem each component solves, the structure to be determined before fitting, and the parameters to be estimated by model.fit().

It should be noted that the symbols \(\rho\) and \(t_0\) were used previously to explain the shape of a simplified Logistic curve. These are mathematical notations for easier understanding and should not be directly regarded as fields with the same names in model.params. Prophet’s complete Logistic trend also includes growth rate adjustments corresponding to changepoints and continuity constraints.

The model structures and input definitions listed in the table need to be determined before fitting, while the parameters in the last column are jointly estimated by model.fit(). Chapter 4 will elaborate on how to convert historical records and known future information into these model inputs, and Chapter 5 will delve into the MAP estimation process for the parameters.

4. Data Processing and Model Configuration: Turning Business Records into Model Inputs #

Chapter 3 introduced how Prophet utilizes trends, seasonality, holidays, and extra regressors to model demand patterns. This chapter will delve into two engineering problems: first, how to convert Store-SKU business records into consistently defined historical and future inputs; and second, how to determine the model structure and candidate hyperparameters by combining sequence features, store operations, and product attributes. Once the inputs and configurations are ready, Chapter 5 will discuss how the model estimates unknown parameters through MAP optimization.

Constructing Consistently Defined Training and Future Inputs #

Organizing Business Records into Store-SKU Daily Data #

We will define demand sequences by Store-SKU combinations, using daily granularity for observation and forecasting. Before model input, raw transaction records need to be processed and unified into daily data. This involves converting transaction times to the store’s local time zone, then aggregating by store, SKU, and operating day, ensuring each combination has only one record per date. The engineering process also needs to check for duplicate records, date gaps, store operating status, and stockout flags, and construct the training column y according to the target definition specified in Chapter 2.

Prophet’s direct inputs include dates, target values, and registered regressors; holiday calendars are provided through model configuration. Operating status and stockout flags are primarily used to construct targets, filter valid observations, and process forecast results. The table below shows an example of actual records for Boxing Day 2023 in the teaching dataset, where the holiday field comes from the corresponding VIC holiday table.

Business data processing and Prophet’s internal processing have different responsibilities:

Maintaining Consistency Between Training and Future Inputs #

The definition of fields, time conventions, and transformation rules must be consistent between the training and prediction phases. ds represents the store’s local business date in both phases; y is only provided during the training phase; extra regressors like discount must have both historical records and cover the entire forecast horizon; holiday and operating calendars also need to extend into the future.

Taking the discount variable as an example, if training data uses 0.20 to represent a 20% discount, then future promotions must also be entered as 0.20 according to the same definition, not 20; otherwise, the model will receive a value scaled by a hundredfold. Even if there is no promotion on a future date, 0.00 must be explicitly provided; one cannot omit a registered discount just because the target value y has not yet occurred. Similarly, if sales are aggregated by Melbourne local date during the training phase, the prediction phase must use the same time zone and operating day boundaries.

Future inputs must be based on information known when the forecast is issued. Calendars can usually be generated in advance; promotions and prices should use the plan version available at that time, and weather inputs should use the forecasts published at that time. If backtesting uses actual promotion, price, or weather data obtained after the forecast cutoff, it introduces data leakage and overstates expected production performance. Production systems should therefore retain historical plan snapshots so that each backtest origin can reconstruct the information genuinely available at the time.

For example, suppose the replenishment team issues a 28-day forecast on December 1. The promotion plan available then specifies a 20% Boxing Day discount, but the store later changes it to 30%. A backtest originating on December 1 should use the 0.20 known at that time, not the 0.30 ultimately executed. Although 30% is the realized value, it was unavailable when the forecast was issued; using it would make offline evaluation look better than actual production conditions.

Stockout and store closure information play different roles in training and prediction: historical flags are used to determine if sales accurately represent demand, while future operating schedules are used for post-processing of predictions; whether a future stockout will occur is generally not a known fact and needs to be assessed separately by inventory and replenishment processes, thus actual stockout labels from the future cannot be used directly.

How Prophet Converts Business Fields into Model Inputs #

Prophet first parses ds into dates, checks the input, and sorts it by time. Then, it constructs the necessary inputs for different model components. To unify notation, this article uses \(t\) to represent a time point in the model (i.e., a row of sorted data) and \(d_t\) to represent the calendar date corresponding to the ds of that row. Here, \(t\) is an index in the model formula, not equivalent to any internal numerical coordinate.

The same \(d_t\) is converted into different inputs depending on the model component: the trend function \(g(t)\) uses normalized trend coordinates \(\tau_t^{(g)}\); the seasonality function \(s(t)\) uses fixed calendar reference points to generate seasonal coordinates \(\tau_t^{(s)}\) and Fourier features; and the holiday function \(h(t)\) directly uses \(d_t\) to match the event calendar. Training and prediction must follow the same set of conversion rules: the trend continues to use the start point and span saved during fitting, and seasonality continues to use fixed calendar reference points and period configurations. [4]

How Dates Are Converted into Trend Inputs #

The trend function \(g(t)\) uses normalized trend coordinates \(\tau_t^{(g)}\) at time point \(t\). Prophet first calculates how far date \(d_t\) is from the training start point, then divides by the span of the entire training history, instead of directly substituting ds into the trend formula. The training start point corresponds to 0, the end point to 1, and future time points can have coordinates greater than 1:

Therefore, the mapping relationship from ds to the trend component can be summarized as:

ds is first parsed into the date \(d_t\) corresponding to time point \(t\), then converted into the normalized trend coordinate \(\tau_t^{(g)}\); the trend function combines this coordinate with the fitted trend parameters to calculate \(g(t)\).

Notation explanation:

The numerator and denominator use the same time unit. The fitting process saves model.start and model.t_scale, which will be reused during prediction, without rescaling for future windows.

Assume the effective training records are only January 1, 2, 3, 2023. The training span is two days:

Trend coordinates exceeding 1 is normal extrapolation. If future dates were rescaled back to 0 and 1, they would incorrectly be placed back within the training start and end points. The same future date should have a fixed coordinate, regardless of the length of the current prediction window.

The start point and span are taken from the effective training history; if target values for boundary dates are excluded, the attributes saved by the model should be used.

Missing dates do not compress the timeline: even if January 2 is missing, January 1 and 3 are still two days apart.

How Dates Are Converted into Seasonality Inputs #

The seasonality function \(s(t)\) uses calendar coordinates \(\tau_t^{(s)}\) at time point \(t\). Prophet first calculates the number of days date \(d_t\) is from a fixed calendar reference point, then converts \(\tau_t^{(s)}\) into periodic angles \(\alpha_{P,n,t}\) based on the period length \(P\) and harmonic number \(n\), which are used to generate corresponding sine and cosine features. This coordinate is different from \(\tau_t^{(g)}\) used by the trend function. Prophet 1.4.0, as used in this article, sets 1970-01-01 as the fixed reference date:

Therefore, the seasonality formula from Chapter 3 can be written here as:

Therefore, the mapping relationship from ds to the seasonality component can be summarized as:

ds is first parsed as date \(d_t\), then converted into seasonal coordinate \(\tau_t^{(s)}\). For each period \(P\) and harmonic \(n\), the model generates a pair of sine and cosine features. It combines them using the fitted \(u_{P,n}\) and \(v_{P,n}\) coefficients to form the contribution \(s_P(t)\) for one seasonal period, then sums components such as weekly and yearly seasonality to obtain \(s(t)\).

Notation explanation:

Each harmonic number generates two columns: sine and cosine. Seasonal coordinates and periods are both measured in days; the normalized trend coordinate \(\tau_t^{(g)}\) cannot be directly substituted here.

Below, we use two Saturdays one week apart and an adjacent Sunday to illustrate how ds is converted into features used by the seasonality function \(s(t)\). Take the first order weekly period \(P=7,n=1\), and define the position of time point \(t\) within the seven-day period as:

For a weekly period \(P=7\), \(\rho_{7,t}=\tau_t^{(s)}\bmod 7\) represents the position of time point \(t\) within the seven-day cycle. Fourier functions only depend on the angle’s position within a full cycle, so equivalent sine and cosine features can be calculated using \((\alpha_{7,1,t}\bmod 2\pi)\). The fixed reference date 1970-01-01 was a Thursday, corresponding to \(\rho_{7,t}=0\); thus, \(\rho_{7,t}=2\) corresponds to Saturday, and \(\rho_{7,t}=3\) corresponds to Sunday. Two Saturdays seven days apart have the same weekly position and first-order Fourier features.

Only \(n=1\) is shown in the table. When the Fourier order for weekly seasonality is 3, each ds will also calculate \(n=2\) and \(n=3\) using the same method, ultimately generating six columns of sine and cosine features. These feature values are determined solely by the date, period, and order; the sales y is used during fitting to estimate their coefficients. When predicting future dates, the model continues to use the same fixed calendar coordinates, so future Saturdays will obtain the same weekly periodic features as historical Saturdays; trend coordinates, however, continue to extend forward, allowing the same weekly position to correspond to different baseline demand levels.

How Holiday Calendars Are Converted into Event Features #

A time point \(t\) in the holiday function \(h(t)\) corresponds to a row in the input data, and its ds field, once parsed, represents the calendar date \(d_t\). Unlike trend and seasonality, the holiday component does not require converting \(d_t\) into a normalized time or periodic angle. Prophet directly compares \(d_t\) with event dates and their windows, converting the matching results into event features.

Let \(\mathcal H\) be the set of events. For an event \(H\in\mathcal H\), let its calendar date be \(d_H\), and the set of relative dates included in its window be \(\mathcal R_H\). The indicator feature for event \(H\) at relative date \(r\) is:

The contributions of all event features are summed to form the holiday component for time point \(t\):

Therefore, the conversion from ds to the holiday component can be summarized as:

Chapter 3 uses the index \(\ell\) to represent event features. Here, each “event name—relative date” combination \((H,r)\) corresponds to an \(\ell\), so \(z_{H,r}(t)\) and \(\kappa_{H,r}\) are expanded forms of \(z_\ell(t)\) and \(\kappa_\ell\), respectively.

Teaching Case: Window Is 0, Only Marking the Holiday Itself #

The data generation formula uses \(28H_t\) to represent the holiday contribution: when the calendar date \(d_t\) appears in the VIC holiday table, \(H_t=1\), otherwise it is 0. This means that, assuming the store operates normally, public holidays uniformly add 28 units. The teaching configuration sets lower_window and upper_window to 0 for all events. Taking Christmas Day on December 25, 2023, as an example:

Christmas Day is chosen here solely to demonstrate how dates are converted into event features; whether a feature is generated and whether the target value for that day is valid are two independent issues.

Therefore, the contribution of Christmas Day to \(h(t)\) will only take effect on December 25; if that day also matches other events, \(h(t)\) will continue to accumulate corresponding contributions. Prophet does not directly use the known coefficient 28 or a single \(H_t\) from the generation formula. Instead, it generates features separately for events like Christmas Day, Boxing Day, etc., and then estimates their respective coefficients from historical sales. The teaching case uses an additive mode and explicitly sets holidays_prior_scale to 1.0.

The business calendar and holiday features are independent. The teaching store is closed on Good Friday and Christmas Day, so the training targets for these dates are set to missing; Boxing Day operates normally and can provide valid holiday sales observations. Store closures do not remove the holiday features for those days, but dates without valid target values cannot help the model identify the corresponding holiday coefficients.

When Extending Windows, Each Relative Date Forms an Independent Feature #

If Christmas Day is changed to lower_window=-1 and upper_window=0, Prophet will generate two separate columns: “day before” and “day of”:

Each ds still occupies only one row in the feature matrix. \(r=-1\) and \(r=0\) form two columns, “day before” and “day of,” respectively, each with independent coefficients to be estimated. Therefore, the model can learn the different impacts of pre-holiday purchasing and the holiday itself.

The 0s and 1s in the table only indicate whether an event feature is active; they do not pre-define how much sales will increase or decrease. When the same event name and relative dates are used across multiple years, the model reuses the same feature columns and estimates the corresponding coefficients from data spanning multiple years [4, 6].

In real retail, Christmas purchasing might start a week or even earlier. Therefore, lower_window could be expanded to −7. However, lower_window=-7, upper_window=0 would generate eight independent features for the seven days before the holiday and the holiday itself, estimating eight separate coefficients, rather than a single unified “pre-Christmas week” effect.

When only a few complete years of historical data are available, there are very few event observations for each column, and they may overlap with weekdays, promotions, and other activities. In practical projects, time series cross-validation should be used to determine if expanded windows improve forecasting. If the business assumption is that the entire week before the holiday shares the same type of impact, it is more appropriate to construct a pre_christmas_week indicator variable manually and input it as an extra regressor into the model.

What Holiday Aspects Are Handled by Prophet by Default, and What Requires Engineer Configuration? #

Prophet 1.4.0’s division of labor is as follows:

Prophet defaults to handling date matching and feature expansion, but it does not determine which events are applicable to a store, how many days a window should cover, or whether the store is open. The holiday table should also maintain stable event names and store versions visible at each backtesting start point.

Before deployment, cross-reference and fix the calendar version against Business Victoria’s 2025 holiday table and 2026 holiday table. As some non-metropolitan areas may substitute local holidays for Melbourne Cup Day, subdiv="VIC" can only serve as a state-level reference; verification by store is still required.

Therefore, engineers are responsible for defining the event calendar and windows, while Prophet is responsible for converting each ds into event features; the coefficients for these features are then estimated during the fitting process described in Chapter 5.

Regressor Standardization and Target Scaling #

Before passing data to the optimizer, Prophet preprocesses additional regressors and the training target: regressors can be standardized, and the target values are scaled according to the scale saved by the model. Both transformations compute statistics based on training data, but they apply to different objects and are interpreted differently. When making predictions, these statistics saved during the training phase must be reused, rather than being recomputed based on future data.

Additional Regressors: Deciding Whether to Standardize #

Let the original value of the \(j\)-th additional regressor at time \(t\) be \(x_{t,j}\). When standardize=True, Prophet uses the training period mean \(\mu_j\) and standard deviation \(\sigma_j\) for standardization:

The process of computing standardization statistics from training data and then generating regressor contributions under the additive configuration described in this article can be summarized as:

Here, \(\mu_j\) and \(\sigma_j\) are preprocessing statistics calculated from the training data, while \(\beta_j\) is the coefficient estimated jointly with trend, seasonality, and holiday parameters during model.fit(). If standardize=False, Prophet will retain the original input values, which is equivalent to using \(\mu_j=0\) and \(\sigma_j=1\) for standardization. The default setting standardize="auto" keeps binary variables as 0/1 and typically standardizes non-binary variables.

In the tutorial data, the discount contribution is generated by \(130q_t\), where \(q_t\) represents the discount percentage at time \(t\); when input into Prophet, this corresponds to \(x_{t,\mathrm{discount}}=q_t\). Since standardize=False is set in this article, a 20% off sale (\(q_t=0.20\)) will enter the model with its original value of 0.20, and its regressor contribution will be \(\beta_{\mathrm{discount}}q_t\). Prophet will not directly use the known coefficient 130 from the generative formula, but will estimate \(\beta_{\mathrm{discount}}\) via model.fit(). If the estimated value is close to 130, then the contribution corresponding to a 20% discount will be approximately \(130\times0.20=26\) units.

Regardless of whether standardization is applied, the same regressor must use the same definition, units, and transformation rules across both the training and prediction intervals. Holiday indicator features and Fourier features are constructed by Prophet according to established rules, and users do not need to standardize them along with additional regressors.

Target Value: Scaling to the Model’s Internal Scale #

This tutorial case uses historical data up to 2025-12-03 to fit the model, and target values for closed days and out-of-stock days are marked as missing before fitting. After filtering, the training set contains a total of 1,039 valid target values. The configuration in this article uses a linear trend and scaling="absmax", with no lower bound offset. Therefore, Prophet divides the valid training sales \(y_t\) at time \(t\) by the maximum absolute value of the training targets:

The mapping from business sales to the model’s internal target is:

The maximum valid training sales occurred on 2025-04-20, at 193 units, so \(s_y=193\). This value is saved by the model as model.y_scale. For example, the valid training sales of 153 units on 2023-01-01 are internally converted to \(\tilde y_t=153\div193\approx0.793\); while 193 units on 2025-04-20 are converted to 1. Sales on closed days, out-of-stock days, and during the test period do not participate in determining this scale. \(s_y\) is a preprocessing statistic calculated from valid training targets, not a model parameter solved by the optimization algorithm.

During prediction, the inverse scale transformation is performed. For the full additive configuration with no lower bound offset used in this article, if the model obtains a point prediction \(\widehat{\tilde y}_t\) on the internal scale, the demand forecast in the original target units is:

For example, the internal point prediction for 2025-12-06 in the tutorial is approximately \(\widehat{\tilde y}_t=0.69465\); multiplying this by \(s_y=193\) restores it to \(\hat y_t\approx134.07\) units. The yhat value returned by predict() has already completed this restoration step, with units of daily items, and therefore should not be multiplied by model.y_scale again.

Notation #

If all valid training target values are 0, Prophet sets \(s_y\) to 1 to avoid division by zero errors. The formulas above apply only to the configuration in this article; if scaling="minmax" is used or a logistic lower bound is set, target values will also be offset and scaled according to the corresponding lower bound and scale.

The predict() function uses the \(\mu_j\) and \(\sigma_j\) saved during the training phase to transform future regressors, and uses \(s_y\) to restore the trend, additive components, and final prediction result; multiplicative components are still expressed as a proportion relative to the trend. Since this article uses a full additive mode, yhat, the trend, and each additive contribution can be interpreted in terms of daily units. Business processes typically perform non-negative truncation, rounding, or case pack conversion based on this; in the code, forecast_units stores these business processing results, while the original yhat values are also retained.

Determine Model Configuration Based on Store and Product Business Characteristics #

Once dates, holiday calendars, discount variables, and target values have been processed, the data is in the required input format for modeling. Next, we will configure Prophet’s trend form, seasonality, events and regressors, and model complexity, taking into account the historical length, observational granularity, demand patterns, operational rules, and future information availability of the store-SKU sequence. Business judgment can rule out clearly unreasonable options, while the remaining candidate configurations need to be compared via time series cross-validation within the training history.

Determine Model Structure from Teaching Case Evidence #

The generation mechanism of the teaching case is known, so the formulas and fields in Chapter 2 can be directly mapped to modeling decisions. In a real project, since the data generation formula is unavailable, we can only propose candidate model structures based on historical data, business records, and planning information available at the time of prediction, and then validate them through time series cross-validation.

The table above determines the model structure and data processing direction, but not the final values for all parameters. The next section will provide a reproducible teaching baseline, including the number of candidate changepoints, Fourier orders, prior scales, and interval settings. These values still need to be validated in Chapter 6 through time series cross-validation. Specific coefficients for trend, seasonality, holidays, and discounts are not specified during initialization but are jointly estimated by model.fit(), whose solution process will be detailed in Chapter 5.

Complete Configuration for the Teaching Case #

The table below classifies the configuration of the teaching case into three categories: model structure, complexity and regularization, and inference and output. These specific values constitute a reproducible teaching baseline, but it does not mean that all store-SKUs should use the same settings.

The period in add_seasonality() is in days, and fourier_order controls the number of basis functions; if a prior is not specified separately, the global seasonality prior is used. The prior_scale of add_regressor() constrains the coefficient of that variable, standardize determines whether the input is transformed, and mode determines whether its contribution enters the forecast as quantities or proportionally. Both seasonality and regressor variables must be registered before calling fit().

Initialize Model and Generate Forecast #

The previous two sections have determined the model structure and the teaching baseline. Now, we can apply these settings to the Prophet model and prepare the fitting input and prediction input separately. Both stages use the same ds and discount definitions, but only the fitting stage can provide the observed target value y [8, 9]:

The model first registers weekly seasonality, yearly seasonality, and the discount regressor, then calls fit(). Prophet 1.4.0 automatically excludes rows where y is missing, so stockout days and closed days that were masked do not participate in parameter estimation. The remaining valid training rows construct the model input according to the date transformation and scaling rules discussed earlier. The complete model constructor has been saved in the Notebook, with key configurations as follows:

model = Prophet(
    growth="linear",
    holidays=vic_holidays,
    n_changepoints=25,
    changepoint_range=0.80,
    weekly_seasonality=False,
    yearly_seasonality=False,
    daily_seasonality=False,
    seasonality_mode="additive",
    changepoint_prior_scale=0.05,
    seasonality_prior_scale=1.0,
    holidays_mode="additive",
    holidays_prior_scale=1.0,
    mcmc_samples=0,
    interval_width=0.80,
    uncertainty_samples=1000,
)
model.add_seasonality("weekly", period=7, fourier_order=3)
model.add_seasonality("yearly", period=365.25, fourier_order=5)
model.add_regressor(
    "discount",
    prior_scale=1.0,
    standardize=False,
    mode="additive",
)
model.fit(train[["ds", "y", "discount"]], seed=42)
forecast = model.predict(test[["ds", "discount"]])

The predict() function generates trend, seasonality, and holiday features based on future ds, and reads the discount plan for corresponding dates, then outputs point forecast yhat, lower interval yhat_lower, upper interval yhat_upper, and each model component. The actual y retained in the test dataset is only for ex-post evaluation and is not passed to predict(). Here, test refers to the last 28-day holdout window, which does not participate in configuration selection within the training history. Chapter 6 will further clarify the boundaries of the validation set, final evaluation window, and independent test period.

The original model output yhat represents the demand forecast assuming normal store operation. The teaching code first performs non-negative truncation, then sets the sales forecast for closed days to 0 based on is_open, and saves the final result as forecast_units. Rounding or carton conversion is not performed here; these steps should be handled separately according to the actual replenishment process. This field design preserves the distinction between demand forecast and sales forecast considering operational constraints.

The figure below shows yhat from the same fitted model across the training history and the final 28-day holdout window. Black points are valid training targets used in fitting. The blue curve to the left of the red dashed line is the in-sample fit; to the right, it is a genuine out-of-sample forecast for the holdout period. The light-blue band is the model’s 80% prediction interval. The red dashed line marks the forecast origin, 2025-12-03, and the following 28 days were not used in this fit. Historical fit is not a substitute for out-of-sample evaluation; model selection and error comparisons must still rely on the rolling-origin backtests and holdout results in Chapter 6.

Interpreting Prophet Model Components #

While the previous total forecast plot answers “how much to forecast daily,” model.plot_components(forecast) reveals how the model allocates the forecast across its components [9]. These curves are useful for examining the model’s structure and contribution allocation. However, they represent only a statistical decomposition of the fitted results and cannot independently prove a causal effect of a specific business factor on sales.

When interpreting these plots, it’s essential to consider both the y-axis units and the x-axis meanings. This article uses a fully additive model, so weekly, yearly, holidays, and extra_regressors_additive are all expressed in units. Negative values indicate a reduction in demand relative to the baseline, not a negative final forecast. If a multiplicative model were used, these components would represent a relative proportion of the trend, requiring calculation according to the hybrid formula in Chapter 3.

The trend, holidays, and regressor panels are plotted along actual dates, while the weekly and yearly panels display their fitted patterns on standard weekly and annual grids, respectively. The red forecast start line applies only to panels with a date-based x-axis. In the discount panel, a plateau indicates the same discount applied over consecutive days. The diagonal lines at the plateau edges are merely artifacts of the plotting program connecting adjacent daily points and do not signify a gradual intra-day change in discount.

For this fully additive configuration, the total forecast satisfies:

Here, additive_terms already aggregates contributions from weekly, yearly, holidays, and discounts, so these components should not be added again when calculating the total forecast. The Notebook also verifies:

The full historical components plot helps examine long-term trends and recurring structures. In this illustrative case, the yearly cycle peaks around mid-April, and the difference between weekend and weekday contributions is close to the 16 units specified in the generating formula. While these correspondences indicate that the model has captured some known structures, they do not substitute for out-of-sample error validation.

Discount contributions appear dense on a full historical scale, so the zoomed-in plot focuses on the two weeks before the forecast start and the final 28 days to provide a clearer view. Seven-day continuous discounts appear as plateaus, while single-day events like Boxing Day show as spikes, with these two shapes arising from different input mechanisms.

A business review should not merely assess curve smoothness but also verify: whether weekly patterns align with store traffic, whether annual patterns are misled by a few promotions, whether holiday effects are conflated with store closures, and whether the trend has absorbed sustained discounts. The model allocates these contributions during joint fitting; when features are correlated or historical data is limited, the allocation results can be unstable. Forecast accuracy and interval coverage still require validation through time-series backtesting in Chapter 6.

The component plots illustrate demand assuming normal store operations. To demonstrate how operational constraints modify the final sales output, we now compare the original yhat with forecast_units after applying the is_open flag:

The zero sales for Christmas Day in the plot are determined by is_open=False; they do not mean that the underlying demand forecast is zero. Once the operating plan is finalized, sales forecasts and their intervals for closed days can be set to [0,0]. The latent-demand forecast and its uncertainty, assuming the store were open, should still be retained separately.

At this point, data preparation, model initialization, forecast generation, and component inspection form a complete workflow. The model.fit() function in the code has already estimated all parameters; Chapter 5 will delve into the internal process of this step, explaining the statistical model, MAP objective, and optimization solution.

5. Optimization and Interpretability: How Parameters Become Daily Forecasts #

In Chapter 4, we completed data processing and finalized the model configuration. This chapter integrates trend, seasonality, holidays, and extra regressors into a unified statistical model, explaining how the likelihood function and priors together form the Maximum A Posteriori (MAP) optimization objective, and how the optimizer jointly estimates the unknown parameters. Finally, we will trace the path from “fitted parameters → feature contributions → model components → point forecasts” to convert internal calculation results back into daily demand forecasts, and reconstruct specific predictions for Saturdays and Boxing Day.

Building a Unified Statistical Model and Separating Fixed Inputs from Estimated Parameters #

Constructing a Unified Statistical Model from Individual Modules #

After merging all modules, the model first calculates the conditional mean \(\mu_t\) for the \(t\)-th data row, which equals the sum of the trend term and all feature contributions. The actual observed sales \(y_t\) may be higher or lower than this mean. This example uses a normal distribution to describe this fluctuation, with its typical magnitude measured by the standard deviation of the observation error. In additive mode:

Notation Explanation:

By arranging \(\mathbf f_t^\top\) for all training time points row by row, we obtain the design matrix \(\mathbf F\); the corresponding coefficients form the vector \(\mathbf b\). The matrix product \(\mathbf F\mathbf b\) can simultaneously calculate the contributions of seasonality, holidays, and extra regressors for all training time points. This article distinguishes these coefficients according to their business meaning, while Prophet’s Stan implementation stores them uniformly in the beta vector [10].

For any given set of trend parameters, feature coefficients, and observation noise scale, the model can calculate the conditional mean and its observation density for each training time point. The fitting process aggregates these daily densities into a likelihood function for the entire training history, which then, along with the parameter priors, forms the Maximum A Posteriori (MAP) optimization objective. Before delving into the specific objective function, the next section will differentiate between quantities estimated by the optimizer and those determined before fitting.

What Quantities Does the Optimizer Solve For, and Which Are Fixed Before Fitting #

Before formulating Prophet as an optimization problem, we need to clearly distinguish between optimization variables and fixed conditions. Optimization variables are the unknown parameters jointly solved by model.fit(); while model structure, prior hyperparameters, and data transformation quantities are determined before the current solution process begins. Taking the educational case as an example:

The positions of candidate changepoints are generated by configuration rules; the optimizer estimates the slope change \(\delta_j\) at these positions. Fourier basis functions are determined by date and period configurations, and their corresponding coefficients are included in the vector \(\mathbf b\). Similarly, prior_scale determines the strength of the prior penalty; the fit() operation solves for parameters constrained by this prior, not prior_scale itself.

Therefore, with the model structure, training inputs, prior scales, and data transformation quantities fixed, the objective of this chapter’s solution can be stated as: finding the parameter set \(\theta=(k,m,\boldsymbol\delta,\mathbf b,\sigma)\) that maximizes the joint likelihood function of the training data and the prior support for the parameters. The next section will formally express this objective as the MAP formula.

From Joint Likelihood to MAP Numerical Solution #

Defining the MAP Objective with Likelihood and Priors #

We have defined the single-day observation distribution. Given a set of parameters \(\theta\), the trend and feature coefficients determine the conditional mean \(\mu_t\), and the noise scale \(\sigma\) determines the dispersion of observations around that mean. For a set \(\mathcal D\) containing valid training observations along with their dates and inputs, and a set \(\mathcal T_{\mathrm{train}}\) of training time points with valid target values, Prophet obtains the likelihood for the entire training history by multiplying the conditional densities for each day:

Here, the single-day density \(p(y_t\mid t,\mathbf f_t,\theta)\) has a mean of \(\mu_t\) and a standard deviation of \(\sigma\).

This product is based on the assumption of conditional independence: given the date features and parameters, the observation errors for each day are mutually independent and follow a normal distribution. The likelihood measures the compatibility of the current parameters with the entire training data.

If we solely pursue higher training likelihood, the model might use excessively large changepoint magnitudes or feature coefficients to over-explain random fluctuations. To avoid this problem, Prophet introduces parameter priors \(p(\theta)\), to reduce the support for complex parameter values. Maximum a Posteriori Estimation (MAP) strikes a balance between fitting the training data and adhering to prior constraints, selecting the set of parameters that maximizes the posterior density:

By taking the logarithm, the product in the joint likelihood becomes a sum. This allows the optimizer to simultaneously weigh the evidence provided by each day’s data and the constraints imposed by the parameter priors. In the next section, we will switch to Prophet’s internal scale and elaborate on the priors and their corresponding optimization penalties.

How Prior Scales Enter the Optimization Objective #

The prior scale determines the concentration of parameter priors around zero. When the MAP maximization problem is rephrased as a negative log-likelihood minimization problem, the same scale manifests as the penalty strength in the objective function. The following sections sequentially explain Prophet’s internal parameters, the impact of priors on demand forecasting, and the objective that the optimizer ultimately minimizes.

Internal Parameters: What the Optimizer Actually Solves For #

Previously, \(\theta\) denoted the complete set of parameters at the business scale. Prophet scales the target values and transforms time coordinates before optimization, so the optimizer actually solves for the set of internal parameters \(\tilde\theta\). Correspondingly, \(\sigma\) and \(\tilde\sigma\) represent the observation noise standard deviations at the business scale and internal scale, respectively.

In the linear, fully additive tutorial example, the optimizer simultaneously adjusts the initial trend, changepoint shifts, all feature coefficients, and the noise scale:

Here, \(\tilde\theta\) is the internal parameterization of \(\theta\); \(\tilde k\) and \(\tilde m\) are the internal trend slope and intercept; \(\tilde{\boldsymbol\delta}\) collects the \(J\) changepoint slope shifts; \(\tilde{\mathbf b}\) collects the \(K\) Fourier, holiday, and extra regressor coefficients; and \(\tilde\sigma>0\) is the internal observation noise standard deviation.

For the \(i\)-th valid training observation, the internal conditional mean equals the sum of the internal trend and contributions from each feature:

Here, \(\tilde g_i\) is the internal trend, \(f_{i,r}^{\mathrm{model}}\) is the value of the \(r\)-th model feature for the \(i\)-th observation, and \(\tilde b_r\) is the corresponding internal coefficient. The latter half of Chapter 5 will convert these internal results back to sales units.

Parameter Priors: How They Impact Future Demand Forecasts #

Priors do not directly constrain future demand, but rather the trend changes and feature coefficients that generate predictions. The training likelihood requires the model to explain historical sales, while parameter priors suppress excessively large changes that lack data support. Together, they determine the fitting amplitude of each component and ultimately influence future 28-day demand forecasts.

Prophet uses a zero-centered Laplace prior for changepoint shifts and zero-centered normal priors for Fourier, holiday, and extra regressor coefficients:

Here, zero is the center of the prior distribution before fitting, not a value that the coefficients must take after fitting. Taking feature coefficients as an example, \(\mathbb E_{\mathrm{prior}}[\tilde b_r]=0\) indicates that before observing training sales, the model favors smaller positive or negative effects. After incorporating training data, the likelihood and prior jointly form the posterior distribution:

Thus, the posterior distribution can deviate from zero. The MAP estimate used in this text takes the location of highest posterior density, not the expected value of the prior distribution. A zero-centered normal prior will shrink coefficients with insufficient evidence towards zero, but it will not fix all coefficients at zero.

\(\tau\) corresponds to changepoint_prior_scale, and \(s_r\) is determined by the prior_scale of the respective feature. A smaller \(\tau\) will cause more changepoint shifts to shrink towards zero, making future trends closer to a stable extension; a smaller \(s_r\) will strengthen the shrinkage of feature coefficients towards zero, weakening the impact of seasonality, holidays, or promotions on future forecasts. Larger scales allow the model to express stronger changes, but also increase the risk of the model tracking a few anomalous observations.

The tutorial data generates discount contributions via \(130q_t\); when \(q_t=0.20\), the hidden true value in the synthetic data is a demand increase of 26 units. Real-world projects do not know this value before fitting; it is visible here only because the data generation mechanism for the tutorial data was predefined by the author. Therefore, the 26 units are only used for post-fitting recovery checks and do not participate in prior scale configuration.

During fitting, Prophet simultaneously estimates the trend, seasonality, holiday, discount coefficients, and noise scale. If, after controlling for other components, sales on promotional days remain consistently higher than comparable non-promotional days, increasing the discount coefficient can reduce the forecast residuals for these dates; however, a larger coefficient also increases the L2 cost \(\tilde b_{\mathrm{discount}}^2/(2s_{\mathrm{discount}}^2)\) corresponding to the zero-centered normal prior. The MAP optimizer balances likelihood improvement against prior cost.

\(s_{\mathrm{discount}}\) is the prior standard deviation configured before fitting, not a parameter estimated by model.fit() from sales data. The tutorial example sets prior_scale=1.0, so \(s_{\mathrm{discount}}=1.0\). Since y_scale=193, and the discount variable uses its original scale in additive mode, this configuration corresponds to a prior standard deviation for the discount coefficient at the original sales scale of:

For a 20% discount, the prior standard deviation for the discount component is:

This indicates that before observing training sales, the model allows a 20% discount to generate a wide range of positive or negative demand changes. 38.6 units describes the prior scale, not a known true promotional uplift, nor a value the model must ultimately predict.

Promotions in the tutorial data occur every 14 days, and the historical data includes both promotional and non-promotional days for the same day of the week. These repeated comparisons provide data evidence for the discount effect, causing the likelihood gain from a positive discount coefficient to outweigh the increased prior cost from the coefficient deviating from zero. MAP ultimately yields an internal discount coefficient of:

Converted back to the original sales scale, the discount coefficient is approximately:

Thus, the fitted contribution for a future 20% discount day is approximately:

26.52 units is close to the hidden true value of 26 units in the synthetic data, indicating that the model recovered the discount effect under the current tutorial data and configuration; it does not prove that the promotional effect in a real-world project would also equal 26 units. If the training history only contains a few promotional days, or if promotions always coincide with a specific day of the week, holiday, or other business event, the data will struggle to separately identify the discount effect, and the prior will pull coefficients with insufficient evidence more strongly towards zero.

Real-world projects cannot know the true value of a promotional effect beforehand. Engineers can construct multiple candidate prior_scale values based on existing promotional analyses, similar store-items, price elasticity studies, or reasonable ranges provided by business experts available before the forecast origin, and then select among them using rolling time-series cross-validation within the training history. The final evaluation window and ex-post observed promotional effects cannot be used to retrospectively adjust priors. The aforementioned scale conversion depends on target scaling, feature encoding, and standardization methods, and cannot be directly replicated for other store-SKU series. changepoint_prior_scale=0.05 similarly acts on internal slope shifts and cannot be interpreted as an increase or decrease of 0.05 units per day.

\(\tau\) and \(s_r\) are hyperparameters given before fitting, while \(\tilde\delta_j\) and \(\tilde b_r\) are the parameters solved by the optimizer. When training data provides sufficient evidence, these fitted parameters can significantly deviate from the prior center.

Joint Prior: How Constraints on Different Parameters Are Combined #

Prophet sets priors for various parameter types separately; under the assumption of independent priors, multiplying these densities yields the joint prior for the complete parameter set:

Here, \(\mathbf s=(s_1,\ldots,s_K)^\top\) collects the prior standard deviations for each feature coefficient in beta. Generally, different \(s_r\) values can be used for different seasonality, holiday, and extra regressor variables. The tutorial example sets seasonality_prior_scale, holidays_prior_scale, and the discount variable’s prior_scale all to 1.0, hence:

Although these internal prior scales are the same, the encoding and value ranges of Fourier, 0/1 holiday, and continuous discount features differ, so they do not represent the same range of business effects when converted to daily demand contributions. The 1.0 here is merely a tutorial baseline, not a verified optimal value. Real-world projects should combine available business information prior to the forecast origin to set a few candidate scales for seasonality, holiday, and discount coefficients, then compare overall error, errors for respective date groups, and component stability using rolling time-series cross-validation within the training history; if individual experiments show that these modules share the explanation of the same segment of sales changes, a few candidate combinations also need to be jointly validated. A prior that is too tight might suppress real changes, while one that is too loose might interpret random fluctuations as trend, seasonality, or event effects; the final evaluation window cannot participate in hyperparameter selection.

\(\mathbf s\) does not include the changepoint prior scale \(\tau=0.05\), nor the prior scales for \(\tilde k\), \(\tilde m\), and \(\tilde\sigma\) individually. That priors can be written as a product does not mean that the parameters are independent after fitting; different parameters may still exhibit posterior correlation due to jointly explaining sales.

Negative Log-Posterior: How It Transforms into the Optimization Objective #

Numerical optimizers typically seek minima, so Prophet rephrases the MAP maximization problem as a negative log-posterior minimization problem; both formulations yield the same optimal parameters:

Here, \(\mathcal L\) is the negative log-posterior objective, \(\tilde{\mathcal D}\) is the training data represented at the model’s internal scale; \(\tilde\theta\) comprises the parameters to be estimated, and \(\tau\) and \(\mathbf s\) are prior scales fixed before fitting. The posterior normalization constant, which does not vary with parameters, is omitted from the equation.

Expanding this objective, the Laplace prior yields \(|\tilde\delta_j|/\tau\), forming an L1 penalty for changepoint shifts; the normal prior yields \(\tilde b_r^2/(2s_r^2)\), forming an L2 penalty for feature coefficients. Combining these with the negative log-likelihood of a normal observation model, and omitting constants unrelated to the parameters to be estimated, yields:

Prophet 1.4.0 uses a normal prior with standard deviation 5 for internal \(\tilde k\) and \(\tilde m\), and a normal prior with standard deviation 0.5 for the positively constrained \(\tilde\sigma\), which are incorporated into \(\mathcal R\). This formula applies to the linear, fully additive tutorial example; constants unrelated to the parameters to be estimated can be omitted after fixing prior scales [10].

In the tutorial configuration, changepoint_prior_scale=0.05 means \(\tau=0.05\), so the weight for the changepoint L1 penalty is \(1/\tau=20\). If this scale is halved, the penalty for the same changepoint shift will double. The penalty for feature coefficients is controlled by \(s_r\); halving \(s_r\) will quadruple the L2 penalty for the same coefficient. Different feature columns can use different \(s_r\) values, so itemizing the penalties is closer to Prophet’s implementation than using a uniform weight.

Thus, the MAP objective combines four types of information into a single optimization problem: training residuals, observation noise, L1 penalties for trend changes, and L2 penalties for feature coefficients. All parameters need to be estimated jointly, so related features may share the explanation of the same segment of sales changes; prior independence does not imply posterior independence after fitting, nor does it signify that the model has identified business causal effects.

How the Optimizer Finds MAP Parameters #

The previous section defined parameter estimation as a minimization problem: finding a set of \(\tilde\theta\) in the internal parameter space that minimizes \(\mathcal L(\tilde\theta)\). Prophet first calculates initial values from the training data, then the Stan optimizer simultaneously updates \(\tilde k\), \(\tilde m\), \(\tilde{\boldsymbol\delta}\), \(\tilde{\mathbf b}\), and \(\tilde\sigma\). These parameters collectively determine the trend, the contribution of each feature, and observational noise, and thus cannot be solved separately in unrelated steps.

In each iteration, Stan computes the gradient of the objective function with respect to each parameter using automatic differentiation. The gradient reflects how training residuals, noise terms, changepoint L1 penalty, and feature coefficient L2 penalty change when parameters are slightly perturbed; the optimizer uses this to select the next set of parameters, causing \(\mathcal L(\tilde\theta)\) to progressively decrease. Iteration stops when the objective value, parameter updates, or gradients satisfy convergence criteria, and the current parameters are taken as the MAP numerical solution [11].

In the Prophet 1.4.0 Python backend verified in this paper, the MAP path typically uses the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) quasi-Newton algorithm. L-BFGS approximates the curvature of the objective function using limited historical update information, making it suitable for models with many parameters. When valid observations are fewer than 100, the backend defaults to the Newton method; when fallback is enabled, L-BFGS abnormal termination may also cause a switch to Newton. These choices are version-dependent low-level solver implementations and are typically not business hyperparameters that retail modelers need to adjust [12].

After optimization converges, Prophet writes the obtained parameters to model.params. The 6 weekly seasonality coefficients, 10 yearly seasonality coefficients, various holiday coefficients, and discount coefficients from the tutorial example are all stored in beta; \(k\), \(m\), and delta describe the trend, and sigma_obs describes observational noise. Normal termination of numerical optimization only means that the algorithm found a solution satisfying the convergence criteria. In engineering practice, it is still necessary to check for abnormal termination, stability under different initial conditions or configurations, and out-of-sample prediction performance.

The fitted parameters do not necessarily have to be exactly 16, 14, or 130 from the generative formula item by item. This is because weekly and yearly patterns need to be reconstructed collectively through multiple Fourier basis functions, and prior shrinkage, random noise, and effective target filtering will also alter parameter allocation. The next part of this chapter will recombine the internal parameters into predictions for specific dates, and Chapter 6 will then examine whether these predictions remain accurate in future windows.

How to Choose Between MAP and MCMC #

mcmc_samples is a non-negative integer parameter set when initializing a Prophet model. Although its name includes “samples,” its primary role is to determine whether the model performs full Bayesian posterior sampling: a default value of 0 means MCMC is not run, while a positive integer enables MCMC and further influences the scale of sampling iterations.

When mcmc_samples=0, Stan finds a set of maximum a posteriori (MAP) parameters through numerical optimization; when mcmc_samples>0, Stan instead uses Markov Chain Monte Carlo (MCMC) to sample from the parameter posterior distribution. Both paths use the same model structure, training data, and parameter priors; the main difference is whether the final result is a single set of MAP parameters or multiple sets of posterior samples.

The choice of path should be dictated by the intended use of the demand forecast output, rather than considering MCMC as a more advanced and inherently more accurate model. While MCMC can preserve forecast differences arising from various parameter combinations, it cannot fix unaddressed stockouts, incorrect promotion plans, overlooked structural changes, or insufficient training information; nor does it guarantee a reduction in MAE, WAPE, or bias. If the business primarily uses point forecasts for future daily demand, MAP is generally recommended; if replenishment decisions rely on service levels, safety stock, or demand quantiles, and parameter uncertainty might significantly affect decisions, then MCMC can be further evaluated. Point forecasts from both paths should be compared using the same rolling forecast origin and forecast window, while prediction intervals require separate checking of actual coverage [4, 13].

Taking the tutorial example, the MAP baseline uses 1,039 valid training observations up to December 3, 2025, to forecast daily demand for the subsequent 28 days. To evaluate MCMC, the training data, trend form, Fourier order, holiday calendar, discount inputs, and all prior scales should remain unchanged, with mcmc_samples changed from 0 to, for example, 300, and the random seed kept fixed. This way, the differences between MAP and MCMC primarily stem from the parameter inference method, rather than changes in data or model structure.

For this diagnostic experiment, one should first check the convergence, effective sample size, and divergent transitions of the 4 chains, then compare MAE, WAPE, bias, and prediction interval coverage on the same set of rolling forecast origins; it’s also useful to specifically check if the demand distributions for Boxing Day and the 20% discount day become notably wider due to the inclusion of parameter uncertainty. The 300 here is merely an initial sampling configuration, not a verified optimal value, nor does it imply that exactly 300 posterior samples will be retained. Prophet 1.4.0 uses half of this parameter for warmup in each chain and the other half for formal sampling.

The Stan backend uses the No-U-Turn Sampler (NUTS) to generate posterior samples. NUTS is built upon Hamiltonian Monte Carlo (HMC): HMC uses the gradient of the posterior density to construct sampling trajectories in the parameter space, reducing inefficient exploration caused by random walks; NUTS, in turn, automatically stops extending the trajectory when it begins to turn back, thus adaptively determining the trajectory length for each iteration. For this example, each retained sample is a set of possible trend, seasonality, holiday, discount, and noise parameters; substituting these parameters separately into the inputs for the next 28 days yields a set of demand forecast distributions [12, 14]. This paper did not run this MCMC control experiment and therefore does not report improvements in prediction accuracy or interval coverage based on it.

uncertainty_samples is the number of simulations used by Prophet to generate prediction intervals during the predict() phase; it is not a tutorial data field nor a model parameter estimated by fit(). The default value in Prophet 1.4.0 is 1000; this paper explicitly set uncertainty_samples=1000 during model initialization to clearly demonstrate the prediction interval generation configuration. Prophet repeatedly simulates future paths based on the fitted results and then extracts the 80% interval specified by interval_width=0.80 from the simulation results. Increasing the number of simulations typically only reduces Monte Carlo variability in interval quantiles and increases forecast computation, without directly improving the yhat point forecast; setting it to 0 prevents the generation of uncertainty intervals [4, 13].

The uncertainty simulated by uncertainty_samples also depends on the parameter inference method. When mcmc_samples=0, Prophet bases its prediction intervals on a single set of MAP parameters, primarily propagating observational noise and model assumptions about future trend changes, without including the full posterior uncertainty of seasonality, holiday, and discount coefficients; when mcmc_samples>0, the forecasts further aggregate future paths under multiple sets of posterior parameters. In short, mcmc_samples determines how unknown parameters are estimated, while uncertainty_samples determines how fitted results are used to approximate future forecast distributions. This paper requires reconstructing a specific set of fitted parameters item by item, therefore mcmc_samples=0 is used; subsequent parameter reconstructions are all based on a single MAP estimate.

Reconstructing Future Predictions from Fitted Parameters #

The previous section explained how Prophet solves for unknown parameters using MAP or MCMC. This section will demonstrate how to reconstruct daily demand forecasts from the fitted results, following the sequence of “future inputs → model features → component contributions → total forecast.” We will use Saturday and Boxing Day from the tutorial case study to examine the calculation process.

How Future Inputs Combine with Fitted Parameters #

For each future date, Prophet first reuses the transformation rules defined in Chapter 4: mapping dates to trend and seasonality coordinates, matching holiday features, and processing regression variables such as discounts in the same way as during the training phase. The model then multiplies these features by the fitted parameters to obtain contributions from trend, seasonality, holidays, and extra regressors, which are then combined into a point forecast yhat. Observation noise is not included in yhat but is used for simulating prediction intervals.

In the default MAP path, the trend of the point forecast continues to extend at the last estimated growth rate after the last candidate changepoint in the training period. Whether future trend changes will occur again remains unknown; Prophet only simulates possible future trend changes according to model assumptions when generating prediction intervals.

Therefore, Prophet can generate forecasts for the next 28 days at once, but the model still requires known or pre-planned business inputs for each target date. For example, after estimating the parametric relationship between discounts and demand during the training phase, the forecasting phase must provide a discount plan for the next 28 days; the model will not infer unprovided promotional arrangements on its own.

How to Reconstruct Saturday and Boxing Day Predictions #

When reconstructing forecasts, it is necessary to distinguish between model configuration, fitted parameters, and daily component contributions. Periods and Fourier orders are configured by engineers, feature coefficients are fitted by the model, and daily component contributions are obtained by multiplying the feature values of the target date by their corresponding coefficients. Prophet stores the internal coefficients for seasonality, holidays, and extra regressors uniformly in beta. For numerical reconstruction, it is necessary to follow the feature column order stored by the model and scale the internal results back to the sales scale according to the target scaling rules.

The tutorial case study uses MAP, a linear trend, a fully additive mode, absmax scaling, and no lower bound offset. For any target date, the trend and each component’s contribution can first be summed on the internal scale, then multiplied by the target scaling factor saved during the training phase to convert the forecast back to units:

where \(\hat y_t\) is the point forecast expressed in original sales units for target date \(t\); \(s_y\) is the training-target scaling factor; \(\tilde g_t\) is the internal trend for that day; \(\mathbf f_t\) is the complete feature row constructed from date \(d_t\) and business inputs according to the model’s stored rules; and \(\tilde{\mathbf b}\) represents the internal beta coefficients. \(\top\) denotes transpose, and the inner product multiplies feature values by corresponding coefficients term by term and then sums them. This formula does not directly apply to multiplicative mode or configurations with a lower-bound offset.

The trend can also be calculated segment-by-segment from an initial straight line. When a changepoint has not yet been reached, its additional contribution is zero; after passing that changepoint, the new velocity change is multiplied by the elapsed time and added to the original trend:

Here, \(\tau_t^{(g)}\) is the normalized trend coordinate defined in Chapter 4; \(\tilde k\), \(\tilde m\), \(\tilde\delta_j\) are the internal trend parameters stored by the model; \(\tilde c_j\) is the normalized \(j\)-th candidate changepoint, and \(J\) is the number of candidates; \(\max\) takes the larger of the two values. This term is zero before the changepoint and grows with time after the changepoint, which is equivalent to the continuous piecewise linear formula in Chapter 3.

Complete reconstruction requires sequentially reusing the model’s stored input transformations, calculating products for each term according to the training feature column order, and then performing a numerical consistency check against the predict() output. Below, we verify the predictions for an ordinary Saturday (December 6, 2025) and Boxing Day (December 26, 2025).

Prophet 1.4.0’s Fourier features are ordered as sin1, cos1, sin2, cos2, etc., corresponding to \(v\) and \(u\) in this text’s notation, respectively. In this example, \(s_y=193\). After scaling the internal weekly Fourier coefficients back to the sales scale, we get:

The subscript 7 indicates a 7-day period, and \(n\) indicates the Fourier order. During calculation, the six weekly period features for the target date should be multiplied by their corresponding coefficients, following the model’s actual stored sine and cosine column order.

The six weekly period contributions for December 6, 2025, are as follows. The source code arranges sine and cosine in pairs; the column order may differ from the formula in this text but must be consistent with the stored coefficients.

The sum of these six terms is 11.1193 units, corresponding to the Saturday component contribution. The actual discount on December 6 was 0.20, and the complete reconstruction value is 134.0674 units; there was no discount on Boxing Day, and the reconstructed value is 129.6224 units. Both reconstructions passed the numerical consistency assertion against predict(). During calculation, the actual values of the feature columns must be used, and the day’s promotional plan cannot be arbitrarily altered for manual calculation convenience.

The model’s forecast for business demand on December 26, 2025 (Boxing Day, a Friday), can be broken down as follows:

Placing the reconstructed components back on the operating calendar again shows the distinction between the raw demand forecast and the final sales forecast after store-opening constraints are applied:

The store is open on Boxing Day, so the operating calendar does not alter the reconstructed forecast of 129.62 units. Christmas Day is different: Prophet’s forecast of approximately 95.52 units represents demand under the counterfactual assumption that the store is open, whereas the final sales forecast becomes zero after the closure rule is applied. Business systems should retain both the raw demand forecast and its components, as well as the sales output after applying operating rules. Actual sales will also reflect random variation, rounding, and on-site supply constraints. At this point, the fitted parameters have been reconstructed into daily forecasts; Chapter 6 next examines their errors, bias, and interval coverage in future windows.

6. Evaluation, Diagnostics, and Parameter Selection: How to Assess Model Reliability #

Chapter 5 converted the fitted parameters into daily demand forecasts for the next 28 days, but producing forecasts does not prove that the model is reliable. Model evaluation requires chronologically partitioned historical windows and a consistent definition of evaluation targets, benchmarks, and metrics. Within that framework, we examine point-forecast error and bias, stability across forecast origins, and empirical prediction-interval coverage. This chapter first establishes the evaluation protocol, then compares benchmarks and candidate configurations through rolling-origin backtesting, and finally uses prediction intervals, residuals, and the known data-generating process to analyze error sources and model boundaries.

This chapter will use three sets of results. The first set comes from the “teaching baseline model” configured in Chapter 4; the second set comes from the “validation-selected model” chosen through internal validation within the training history; the third set comes from a separately constructed “out-of-stock mechanism experiment,” used to analyze target value bias when high demand is more likely to trigger stockouts. These three sets of results cannot be directly mixed and compared due to their differing data sources or evaluation objectives.

Start by Defining the Evaluation Protocol and Data Boundaries #

Before calculating metrics, it’s crucial to define the purpose of each historical window: which are for configuration selection, which for model performance evaluation, and which experiments are solely for explaining model mechanisms. If the same data is used for both parameter selection and reporting final results, evaluation metrics will typically be overly optimistic. This chapter adopts the following protocol:

Although the final 28 days were not used for configuration selection in the main text, their results were examined while the tutorial was being developed. They therefore serve only as a teaching evaluation and cannot be treated as an independent test set. Real projects should reserve a test period that has not influenced any development decision; comparisons with conditional means and latent demand are possible here only because the data is simulated.

How to Select Evaluation Metrics and Benchmark Methods #

Metrics are only comparable when the evaluation objects and data scope are consistent. Let \(\mathcal T_{\mathrm{eval}}\) be the set of evaluation dates. This example only includes operating, in-stock, and valid target value dates. \(y(t)\) denotes the actual sales on day \(t\), and \(\hat y(t)\) denotes the forecast_units corresponding to these sales after non-negative processing and operating rule constraints.

Given the instability of Mean Absolute Percentage Error (MAPE) near zero sales, this paper primarily uses Mean Absolute Error (MAE), Weighted Absolute Percentage Error (WAPE), and Aggregate Percentage Bias (Bias) for evaluation. Here, Bias is an aggregate metric for forecast error, not the statistical bias of an estimator.

MAE retains the sales unit, answering how many units the model deviates on average per day:

WAPE divides the total absolute error during the evaluation period by the total actual sales, facilitating comparisons between windows or series with different sales magnitudes. For example, a cumulative absolute error of 50 units with actual sales of 1,000 units results in a WAPE of 5%:

Bias retains the sign of “forecast minus actual,” used to determine whether the model generally tends to overestimate or underestimate:

MAE and WAPE measure error magnitude without distinguishing between overestimation or underestimation; a positive Bias indicates overall overestimation, while a negative Bias indicates overall underestimation, thus it should be reported alongside MAE or WAPE. When the evaluation set is empty, none of the three metrics are defined; when the total actual sales are zero, WAPE and Bias are not defined.

In addition to the above metrics, evaluation also requires a sufficiently simple benchmark for reference. This paper uses a seasonal naïve baseline: for each day of the week, it takes the most recent valid sales before the forecast cutoff date and repeats it across the entire forecast window; if the corresponding day of the week in the most recent week lacks a valid target value, it looks back to an earlier week. The benchmark can only use information available before the cutoff date. Complex models are only valuable if they consistently outperform the benchmark in the same backtesting window and with the same evaluation scope, justifying the additional training and maintenance costs. Chapter 7 will compare Prophet, LightGBM, and this baseline on the same evaluation dates.

All candidate models and benchmark methods must adopt the same operating day, out-of-stock filtering, and sales post-processing rules; otherwise, metric differences may stem from the evaluation scope rather than the model itself. For example, in this case, Prophet’s MAE on operating days is 4.81 units; if closed days where both predicted and actual sales are zero are included, the MAE drops to 4.64 units, while WAPE remains unchanged. The decrease in MAE results from adding zero-error days to the evaluation set and does not represent an improvement in the model’s inherent performance.

A single aggregated metric is insufficient to support replenishment decisions. Retail systems should separately observe SKU-level errors, category-level biases, performance of high-volume and long-tail products, and cumulative demand error within the replenishment cycle. Replenishment teams typically focus more on total demand for the next two weeks rather than the forecast error for a single day. Counterfactual demand on out-of-stock days or closed days requires additional validation information and cannot be directly evaluated using recorded sales.

How to Check Stability with Rolling-Origin Backtesting #

Randomly splitting training and test sets disrupts time order and may indirectly allow the model to access future information. Retail forecasting is better suited for rolling-origin backtesting: successively select multiple historical cutoff dates, refit the model using only the data visible before each cutoff date, and then forecast the future 28 days that have not yet occurred. This simulates model performance when deployed at different times and checks whether errors change significantly with the forecast origin or forecast horizon [15].

Prophet’s cross_validation() copies the current model’s structure and configuration at each historical cutoff date and refits it using data prior to that cutoff, rather than allowing the current set of fitted parameters to “look back” at history directly. However, the built-in backtesting reads additional regressors from the passed historical table. For future business inputs like discounts, the values in each backtesting window must represent promotional plans already determined at that time; if the historical table only stores actual discounts that occurred afterward, direct backtesting will leak future information to the model. Only when planned values can indeed be obtained in advance, or when plan snapshots for each cutoff date can be recovered, is it suitable to use the built-in diagnostics below directly:

from prophet.diagnostics import cross_validation, performance_metrics

cv = cross_validation(
    model,
    initial="730 days",
    period="28 days",
    horizon="28 days",
    parallel=None,
)
assert cv["cutoff"].nunique() == 12

scores = performance_metrics(cv, rolling_window=0.1)
print(scores[["horizon", "mae", "rmse", "coverage"]].tail())

This set of settings corresponds to the time structure of the teaching case: initial="730 days" makes the first fit use approximately two years of history, covering two annual cycles; period="28 days" means the forecast origin moves forward every 28 days; horizon="28 days" corresponds to the task of forecasting demand for the next four weeks. This results in 12 forecast origins. parallel=None simply indicates sequential execution of each fit, without altering the backtesting protocol.

cross_validation() generates out-of-sample forecast records corresponding to each cutoff and target date ds; performance_metrics() then calculates built-in metrics such as MAE, Root Mean Squared Error (RMSE), and prediction interval coverage from these records. The WAPE, Bias, operating day filtering, and out-of-stock filtering used in this paper are not automatically performed by this function and need to be calculated separately on the cv details according to the unified scope from the previous section.

performance_metrics(cv, rolling_window=0.1) first arranges records by forecast horizon and then performs a rolling aggregation over approximately 10% of the forecast-actual pairs, so a given horizon in the results does not solely represent that exact step. To compare errors for days 1, 7, and 28 separately, one should first calculate horizon_days=(ds-cutoff).dt.days and then group by the exact number of days [15].

If it’s not possible to recover the prices or promotional plans visible at each cutoff date from the historical table, one should write a custom backtesting loop: load corresponding data and plan snapshots by cutoff date, retrain the model, and then forecast for that window. The information boundaries used for backtesting must be consistent with what would be available at the actual deployment time.

Next, independent controlled experiments will be used to differentiate the effects of candidate changepoint ranges, prior scales, Fourier orders, and seasonal patterns, using only the validation window within the training history to select configurations. The last 28 days evaluation window does not participate in this selection process.

How Configuration and Model Structure Change Prediction Results #

This section explains how different configurations and model structures affect prediction results through four sets of comparative experiments. The seasonality and ablation experiments use two internal prediction start points from the teaching case: October 8, 2025, and November 5, 2025. Each prediction spans 28 days, yielding a total of 56 valid evaluation samples. The trend experiment and the additive/multiplicative experiment, respectively, use independently generated data and their final 28-day validation periods. The evaluation targets for different experiments may be observed sales or known conditional means during the generation process, so numerical values should not be directly compared across experiments. These experiments aim to explain model mechanisms and narrow down candidate options; they do not replace the configuration selection in the next section, nor do they use the final 28-day teaching evaluation window to select parameters.

Trend Changepoint Location and Magnitude of Change #

The trend experiment independently generates 700 days of data, increasing the growth slope from 0.025 units/day to 0.275 units/day on day 600. The first 672 days are used for training, and the subsequent 28 days for validation. changepoint_range (c-range) determines where candidate changepoints can be distributed within the training history, and changepoint_prior_scale (cps) determines how much the model allows the slope to change. When changepoint_range=0.8, the candidate range ends before the actual turning point; when extended to 0.95, candidate changepoints can cover the vicinity of the turning point.

The MAE in the table is calculated relative to observed values after adding random noise; the black line in the figure represents the known trend without noise. Keeping changepoint_range=0.8 and increasing the prior scale from 0.05 to 0.5, the validation MAE is still 5.88 units. Expanding the candidate range to 0.95 and using a prior scale of 0.05, the validation MAE drops to 1.33 units. This result indicates that loosening the range of slope change cannot compensate for the issue of candidate locations not covering the true turning point. The candidate range and prior scale need to be selected in conjunction and judged through historical validation. The gray line in the figure marks the true turning point, the red line marks the training cutoff date, and other curves represent the trends fitted by each configuration.

Seasonal Complexity and Prior Shrinkage #

Fourier order determines how complex a curve seasonality can express, and seasonality_prior_scale determines how strongly the corresponding coefficients are regularized. The teaching case uses the same two internal prediction start points, with the following results:

These four configurations do not cover all combinations of order and prior scale, so when comparing, one parameter needs to be kept constant. When the prior scale is 0.01, annual Fourier orders 1 and 5 can be compared. When the annual Fourier order is 5, prior scales 0.01 and 1.0 can be compared. When the prior scale is 1.0, annual Fourier orders 5 and 15 can be compared. In the last comparison group, the training MAE of the 15-order model is slightly lower, but the validation WAPE is higher than that of the 5-order model, indicating that while a more complex curve can improve historical fit, it does not lead to better out-of-sample predictions. Different products still need to select their order and prior scale based on their respective historical validation results. Noise level, target value scale, and sample size also influence the selection results.

Ablation Experiment Removing Model Components #

Ablation studies start from the teaching benchmark model, removing only one source of information at a time: removing discount means no longer adding the discount variable via add_regressor("discount"); removing holidays means not passing the holiday calendar when initializing the model; removing annual seasonality means no longer adding the custom annual seasonality defined in this article. Other configurations, training data, prediction start points, and evaluation metrics remain consistent, and models are refitted separately.

After removing discounts, holidays, or annual seasonality, validation errors all increase, but to varying degrees. Training MAE and validation WAPE have different units and purposes; models should be compared column by column, and 4.55 units cannot be directly compared with 4.01%.

Overall metrics may also mask large errors on sparse events. The two internal validation windows comprise 56 valid samples in total, but only 1 of them is a holiday sample: the complete teaching benchmark model’s MAE on that day was 2.48 units, increasing to 31.28 units after removing holidays. Meanwhile, the overall WAPE only increased from 4.01% to 4.65%. A single sample is insufficient to stably estimate model performance on all holidays, but this comparison shows that when evaluating sparse events, grouping rules, sample size, and group errors must be reported simultaneously.

Since each variant is refitted, the parameters of the remaining components are also jointly adjusted. Therefore, ablation results measure the change in the entire model’s prediction after removing one information source, rather than simply subtracting the original component’s contribution from the complete model’s yhat. Further examination of errors on promotion days, holidays, and ordinary dates, as well as changes in other components around Boxing Day, can help determine whether missing structures are absorbed by trends or other seasonal components. An increase in validation error indicates that the component provided out-of-sample predictive information in the current data and model, but it cannot be concluded that the corresponding business factor has a causal effect of the same magnitude.

Additive vs. Multiplicative Seasonality #

To compare the “additive” and “multiplicative” modes of seasonality_mode separately, this experiment constructed two daily store-item demand series. Both series use the same date ds, growth trend, and weekly seasonality pattern. The only difference is whether weekly fluctuations enter the conditional mean as a fixed number of units or a fixed proportion. Let the base demand be \(\ell_t=50+0.15t\), and the weekly seasonality pattern be \(W_t=\mathbb{I}(\text{weekend})-2/7\); thus, weekdays take \(-2/7\), weekends take \(5/7\), and the average for a full week is 0.

Both series generate observed targets according to \(y_t=\mu_t+\varepsilon_t\), where \(\varepsilon_t\sim\mathcal N(0,2^2)\). Corresponding to the data expression in this teaching material, ds is the daily date, and y is the target value actually used by the model; \(\mu_t\) serves only as a diagnostic ground truth in the simulation experiment and is not used in model training. To isolate the seasonal form, the experiment does not include discount, holiday, store closure, or stock-out mechanisms, nor does it add annual seasonality or extra regressors.

Each series contains 730 days; the first 702 days are used for training, and the final 28 days for validation. The two fits only change seasonality_mode, keeping other model configurations consistent; therefore, the results mainly reflect whether the seasonal form matches the data generation mechanism.

The table below uses the noiseless conditional mean \(\mu_t\) of the final 28 days to calculate MAE, aiming to check whether the model can recover systematic trend and weekly seasonality structures, avoiding a single instance of random noise masking the differences between the two modes. The evaluation target for this metric differs from the MAE based on observed sales earlier, so the numerical values cannot be directly compared.

The black line in the figure represents the known conditional mean, and other curves represent predictions under the two seasonality_mode settings. For fixed unit fluctuations, the additive model is closer to the conditional mean. For proportional fluctuations, the multiplicative model can proportionally expand the weekly seasonality amplitude with the base demand, which is consistent with the MAE results in the table.

Real business scenarios do not have directly observable generation formulas, so it is only possible to first check whether seasonal fluctuation amplitude changes with demand levels, and then compare the two configurations via rolling origin backtesting. It should also be noted that seasonality_mode is a global configuration, and holidays and extra regressors whose mode is not separately specified will also adopt it. This experiment only includes trend and weekly seasonality, so the differences can be attributed to the seasonal form. In a complete business model, it should also be checked whether holidays and extra regressors use additive or multiplicative modes.

The above experiments explain how parameters and components affect prediction results but do not provide an “optimal configuration” that can be directly deployed. The next section will define a candidate set and formally select the configuration for the teaching case using time windows within the training history.

How to Select Configurations within Training History #

Initial comparative experiments aim to understand parameter effects; for formal configuration selection, we need to limit the candidate set and evaluation rules. This article presupposes three Prophet candidate models and, based on the tutorial case’s training history, selects only October 8, 2025, and November 5, 2025, as forecast origins. Each origin forecasts 28 days, totaling 56 valid samples, and their WAPE is ultimately aggregated. The other structures and evaluation criteria for the three candidate models remain consistent. Note: cps: changepoint_prior_scale, sps: seasonality_prior_scale.

This section adopts the following selection rule: prioritize the candidate model with the lowest internal validation WAPE; if WAPE is identical, process according to a fixed candidate order. This calculation process does not rely on the target values of the last 28 days. Based on this, Candidate C is selected as the “validation-selected model.” After model selection, it will be re-fitted using the complete training history up to December 3, 2025, and used to forecast data for December 4-31.

The internal validation WAPE for all three candidate models is close to 4%, with Candidate C’s advantage being only about 0.03 percentage points. While this result supports selecting C for now, it is insufficient to prove that a lower annual periodicity order and prior scale consistently perform better across other time windows or product series. Candidate C’s first-order annual periodicity happens to be close to the data generation mechanism in the tutorial, but in real business, this ground truth is unknown, so more historical origins and representative store-product series are needed for verification.

After fitting with Candidate C, the WAPE for operating days in the last 28-day tutorial evaluation was 4.21%, almost identical to the tutorial baseline model’s 4.22%. This indicates that the internal selection process standardizes configuration decisions but cannot guarantee significant improvement in a single future forecast window. Given that the results for this window were viewed during the tutorial development process, it serves only as a tutorial evaluation and must not be used again to expand the candidate set or adjust selection rules.

In a real project, candidate configurations, rolling origins, and selection metrics should be fixed before development, and each validation result and selection basis should be saved. After internal selection is complete, an independent test period (holdout period) not involved in any development decisions must be used to evaluate the final model; if the model continues to be modified based on this test period, these data will lose their value for independent testing.

How to Interpret and Test Prediction Intervals #

Point forecasts aim to answer “how many units are expected to be sold,” while prediction intervals describe the range within which future observations might fall under the model’s assumptions. This section will revert to the tutorial baseline model (i.e., annual periodicity order 5, seasonality prior scale 1.0) to maintain consistency with previous component and parameter explanations. The validation-selected model from the previous section is not used here, nor are configurations re-selected.

Under the MAP path with mcmc_samples=0, Prophet’s prediction intervals primarily propagate observation noise and future trend changes under model assumptions, and do not include the full posterior uncertainty of seasonality, holidays, and discount coefficients. Future trend simulation also assumes that historical slope changes can represent the future; intervals do not automatically cover unknown promotions, unusual shocks, stockout demand recovery errors, or omitted model structures [13].

In this example, we re-fit the tutorial baseline model at each historical origin, keeping uncertainty_samples=1000 constant, and only setting interval_width to 0.80 and 0.95 respectively. For the same forecast origin, MAP parameters and yhat are identical; the only change is the upper and lower quantiles extracted from the simulated forecast distribution. While increasing uncertainty_samples can reduce simulation variability in quantile estimation, it does not directly improve the accuracy of point forecasts.

Prediction intervals need to be checked for both actual coverage and average width. Let the lower and upper bounds be \(\hat y_t^{\mathrm{lower}}\) and \(\hat y_t^{\mathrm{upper}}\) respectively. The coverage rate on the evaluation set is:

The average interval width is:

Higher coverage is not necessarily more useful, as infinitely widening the interval can also cover more observations. Effective prediction intervals need to be sufficiently narrow while being close to the nominal coverage to provide discriminatory power for replenishment decisions.

Interval backtesting uses 12 historical origins, evaluating only 323 samples from operating days without stockouts:

The figure shows the results for the last rolling origin, while the table summarizes data from all 12 origins. The actual coverage for the 80% interval is 75.85%, differing by 4.15 percentage points from the nominal level; the 95% interval actually covers 92.88%, differing by 2.12 percentage points. When the nominal coverage increases from 80% to 95%, the average interval width also increases from 14.78 units to 22.73 units, illustrating the trade-off between coverage and interval width. These 323 samples provide empirical calibration results for the current tutorial series, but in actual projects, such discrepancies still need to be continuously monitored across more time windows and store-product series.

To pinpoint calibration discrepancies, coverage and average width can be further calculated by forecast horizon, promotion status, holidays, and demand scale. For example, analyzing the grouped results for the 28th day and sparse holidays can reveal the composition of overall coverage in specific business scenarios.

In inventory decisions, prediction intervals are a key input for calculating service levels. When formulating a 95% service level strategy, it is necessary to estimate the cumulative demand distribution within the replenishment lead time, and comprehensively consider daily error correlation, replenishment batch size, stockout cost, and backlog cost. Business systems can input the single-day yhat_upper along with these constraints into the inventory optimization process to calculate safety stock and order quantities.

How to Diagnose Sources of Forecast Error #

Forecast errors typically arise from three mechanisms: future random fluctuations, systematic structures not yet captured by the model, and restricted target values due to store closures or stockouts. To diagnose error sources, evidence corresponding to each mechanism is needed. This section will analyze the first two types of errors using the tutorial baseline model and simulated ground truth, and analyze restricted sales through a separately constructed stockout mechanism experiment, explaining the data and purpose of both experiments.

Random Noise vs. Systematic Error #

Simulated data allows us to compare the same set of forecasts against noise-free conditional means, latent demand with noise, and final observed sales. In the last 28-day evaluation, we focused only on operating days and obtained the following results:

The model’s error relative to the conditional mean is minimal, but significantly increases relative to latent demand and observed sales. This set of differences indicates that the tutorial baseline model has recovered most of the known systematic structures, and the remaining error primarily stems from random perturbations and sales rounding. This forecast window includes one store closure day, and no stockouts occurred during business hours; the metrics in the table use only operating day data corresponding to sales. Simulated data provides conditional mean and latent demand as diagnostic ground truths; in real projects, sufficient inventory periods, controlled experiments, or demand recovery results can be used to establish approximate references.

Restricted Sales and Selection Bias from Demand-Related Stockouts #

The original stockout flags in the tutorial case were independent of demand, suitable for demonstrating how to mask restricted sales. To further observe the selection bias formed by “higher demand leading to higher likelihood of stockout,” this section additionally generated 730 days of demand data and fixed available supply at 115 units. When latent demand exceeded 115 units, the system still recorded 115 units of sales and flagged it as a stockout; the true demand could only be determined to be at least 115 units, so such records are right-censored observations, not truncated observations excluded from the sample. The first 702 days were used for training, and the last 28 days for evaluation, with the evaluation target being the latent demand known during generation.

When right-censored sales are treated as accurate demand, the model retains all 702 training labels, but compresses latent demand exceeding 115 units to the supply cap, resulting in a 6.74% WAPE and -3.74% Bias. After masking stockout days, the number of available labels drops to 580. At this point, WAPE is 6.89% and Bias is -4.74%. The latter set of samples retains only dates where demand did not exceed supply, thus the data is more concentrated at lower demand levels, leading to a further negative shift in forecast bias.

From a statistical learning perspective, masked target values are Missing Not At Random (MNAR): the higher the latent demand, the more likely the target value is to be removed. The stockout flag actually provides lower-bound information, “latent demand is at least 115 units”; directly deleting observations simultaneously loses this valuable information. If the daily available supply cap is reliable, actual projects can employ censored likelihood, inventory-constrained demand models, or dedicated demand recovery processes, combined with inventory, arrivals, shelf availability, search or order signals, and substitute purchase information, to estimate true demand during stockouts.

This set of results reveals how demand-related stockouts create selection bias. In real projects, repeated validation across multiple supply levels, stockout rates, and historical forecast origins can evaluate the stability of demand recovery methods across different store-product series.

Omitted Structures in Residuals #

A residual is defined as the observed value minus the fitted value. If residuals are consistently positive or negative over time, it usually indicates systematic underestimation or overestimation; if they vary repeatedly at fixed intervals, it might suggest omitted seasonality; if they consistently shift in the same direction after a promotion, it could indicate promotion lag or pull-back effects; if the magnitude of residuals expands with demand level, heteroscedasticity or proportional fluctuations might be present. Furthermore, a few extreme residuals should be reviewed in conjunction with stockouts, unusual orders, data quality, and concurrent business events.

The figure above first shows the in-sample residuals of the tutorial baseline model by date, then calculates the calendar-lag correlation for natural days separated by 1-28 days. Store closures or restricted target values remain missing here, and each lag is paired only with residuals from dates with the same actual calendar interval, thus preserving the original calendar distance. This calculation method is particularly suitable for retail series with date gaps, differing from standard equally-spaced autocorrelation calculated after compressing the time axis.

In this example, the absolute values of the calendar-lag correlation coefficients for 1-28 days are all below approximately 0.08, and the peaks near 7, 14, and 28 days are also weak. This indicates that after fitting the tutorial baseline model, fewer linear periodic signals are retained, and the residual pattern approximates random fluctuations around 0. To further assess the white noise assumption, mean and variance stability, along with confidence intervals, significance tests, and multiple comparison control, still need to be checked.

The value of residual diagnostics lies in guiding verifiable model modifications: for example, a peak near 7 days can prompt engineers to re-examine weekly seasonality or operating patterns; continuous shifts after promotions can be converted into lagged variables or event windows; expanding fluctuations with demand levels can guide comparisons of additive vs. multiplicative modes or observation distributions; and extreme residuals concentrated on stockout days can drive target value filtering and demand recovery. Each modification should be re-backtested using the same rolling origins, and out-of-sample error should be the criterion for whether to retain it.

When reporting grouped errors for promotion days, holidays, and ordinary dates, grouping rules, overlap relationships, and sample sizes should be provided simultaneously. Chapter 7 will combine these groupings and exact forecast steps to compare Prophet, tree models, and seasonal naive baselines.

7. Model Extensions and Comparisons: When to Keep Using Prophet #

The first six chapters introduced Prophet’s fitting, interpretation, and evaluation. This chapter will explore four aspects to determine whether to continue using Prophet: whether a single store-item series provides sufficient signals, whether future business variables can be obtained at the time of forecasting, whether historical sales dependencies require other models to express them, and whether Prophet can outperform simple benchmarks and other candidates under a unified backtesting protocol. After these determinations, Chapter 8 will discuss system delivery.

How to Select Forecasting Methods Based on Data Conditions #

Distinguishing True Zero Demand, Constrained Sales, and Missing Records #

Before assessing data sufficiency, it’s crucial to confirm the business meaning of zero values. When a product is normally for sale, available for purchase, and adequately supplied, no transactions can be recorded as zero sales; under the simplified assumptions of this paper, such records are also used as observations of zero demand. Constrained sales due to unlisted products, discontinued products, store closures, out-of-stock situations, and missing observations due to interface failures should each retain status flags. After excluding these statuses, if non-zero demand still appears at irregular intervals, it constitutes a candidate for intermittent demand.

After cleaning business statuses, if a series still contains a large number of true zero demands, Prophet’s modeling approach for conditional means should be examined to see if it suits the current data. For series with sufficient non-zero events, Prophet can be used as a baseline for the average demand rate and compared with intermittent demand methods such as Croston, SBA, TSB, or two-part models [16, 17, 18, 19]. Appendix A defines metrics describing demand intervals and positive demand volume fluctuations, Appendix B uses these metrics to compare four teaching scenarios, and Appendix C introduces the basic mechanisms and formulas of each candidate method.

Are Effective History and Non-Zero Events Sufficient? #

Prophet can fit data with gaps in dates and can treat zero sales as numerical observations. However, to stably estimate the contributions of trends, seasonality, holidays, and business variables, a sufficiently long effective history, a sufficient number of effective observations, and an out-of-sample evaluation window covering the relevant patterns are required. Successful numerical optimization only indicates that the model has solved for its parameters; the number of effective events and the historical coverage determine whether these parameters have learnable, verifiable business meaning.

A high proportion of zero values and insufficient available data are two related but distinct problems. Even if a series has mostly zero values, as long as the observation period is long enough and non-zero demand reappears repeatedly, it may still provide learnable and verifiable signals; new products or very low-frequency items, even with a similar proportion of zero values, might only have a few non-zero sales. The retail industry includes high-frequency necessities, seasonal products, new products, and long-tail items; model selection should be based on the effective history and number of events provided by each specific store-item series.

How to Choose Forecasting Granularity and Information Sharing Methods #

When non-zero events for a single series are insufficient, daily demand can be aggregated into weekly demand, single-store demand can be aggregated into store-group demand, or a global model capable of sharing information across stores and products can be adopted. New and very low-frequency products can first use category baselines, similar product mappings, or business fallback rules, and then evaluate independent local models after sufficient events have accumulated.

Aggregation also changes the forecast target. If replenishment decisions are made at the store-group-weekly granularity, Prophet, statistical models, and tree models can be directly compared at that granularity; if decisions still require store-daily forecasts, then temporal or store-dimensional disaggregation must be designed, and consistency constraints for hierarchical or grouped forecasts should be considered. The evaluation should cover the entire process of “aggregate forecast—disaggregation or reconciliation—replenishment decision”, not just the forecast error at the aggregated level [20, 21].

Therefore, an initial routing can be performed based on data sufficiency and actual decision granularity:

Modeling begins with examining data status, effective history, and non-zero events, then determining forecasting granularity, information sharing methods, and candidate models. All approaches should be compared at the actual decision granularity, using the same information set and a rolling forecast protocol. This paper has not yet run experiments on intermittent demand methods, global models, or hierarchical reconciliation; therefore, this section only describes the selection paths and applicable boundaries, with performance ranking to be provided by future experiments.

How to Determine if Business Variables Can Be Used for Future Forecasts #

After determining forecasting granularity and candidate models, it’s also necessary to define the information available at the time the forecast is initiated. Prophet’s trend, seasonality, and holiday components can be directly calculated based on future dates; whereas additional variables like price, discount, display, and weather must have their values for each future target day provided before calling predict(). Engineers need to register these variables via add_regressor() before fitting and ensure that training and prediction tables use the same field definitions, numerical standards, and missing value handling rules [6].

The vector \(\mathbf{x}_t\) mentioned earlier is the collection of these additional variables. Whether a variable is suitable for demand forecasting critically depends on whether its values for each future target day can be obtained at the time the forecast is initiated; ex-post records in historical data need to be converted into plan snapshots visible at the time, external forecasts, or scenario inputs.

The discount plan in the teaching case study was already determined when the forecast was initiated, so it can be used as an extra regressor. Both the training and future periods use a 0-1 ratio to represent the discount, with standardize=False. In real projects, prices, displays, and advertising expenditures should also save the planned versions visible at the forecast start point; ex-post execution records can only be used for forecasts initiated after their corresponding time.

Whether new continuous variables are standardized, whether an additive or multiplicative mode is adopted, and how original scale coefficients are interpreted, follow the rules from Chapters 4 and 5. Preprocessing parameters can only be estimated from the corresponding training window; highly correlated variables require checking for identifiability, and non-linear or interactive relationships can be verified by explicitly constructing features or in later comparisons with tree models. All new variables should be tested for incremental value through time-series backtesting and ablation experiments.

If variables like weather or foot traffic come from external forecasts, their predicted values can be used as Prophet inputs, or multiple business scenarios can be set up to calculate demand separately [22]. Prediction errors of external variables will propagate into demand forecasts, so backtesting should cover the entire pipeline of “external variable forecast—demand forecast” and retain the input version actually visible at each forecast start point.

Business plans and external forecasts address the question of where to obtain future inputs; sales lags and rolling statistics also involve another temporal boundary: as the forecast horizon increases, some sales required to construct features for the target day have not yet occurred. The next section will specifically discuss the construction methods and information leakage risks of such features in 28-day multi-step forecasting.

How Sales Lag Features Maintain Temporal Boundaries in Multi-Step Forecasting #

Sales lags and rolling statistics can express recent demand levels, but their availability in multi-step forecasting changes with the forecast horizon. Assume the forecast start point is the end of business on day \(d\), and the target day is \(d+h\). The sales feature for the day before the target day is \(y(d+h-1)\): when \(h=1\), the model uses the observed \(y(d)\); when \(h\geq2\), the required sales occur after the forecast start point and have not yet happened when the forecast is initiated. The average of the seven days before the target day will also progressively include unknown dates as \(h\) increases.

Therefore, lag features must be reconstructed around each forecast start point. If shift(1) or rolling calculations are first applied to a complete data table including the evaluation window, and then the next 28 days are sliced out, farther target days will read true sales from after the forecast start point, leading to temporal information leakage. Training, validation, and online processes should all record the forecast start point and item-by-item check whether each target day’s input was visible at that time.

When using recent sales information, the following designs can be chosen based on business update frequency and forecasting tasks:

Another design is to first use Prophet to generate a baseline for trends, calendar, and business plans, and then have a second model learn the remaining predictable residual structure beyond the baseline. Training input for the second stage should come from out-of-sample residuals generated by historical rolling forecasts, not in-sample fitted residuals from the same training period; a complete backtest should also run both stages simultaneously to reflect actual errors in the future window.

How ARIMA Expresses Temporal Dependencies in Sales and Errors #

The previous section discussed how to explicitly construct sales lag features. ARIMA (AutoRegressive Integrated Moving Average) provides another way to express temporal dependencies: the differencing component handles changes in the series level, the autoregressive component uses historical values of the differenced series, and the moving average component describes how historical innovation terms—errors that are new in the current period and were unpredictable before—affect the current value [23]. The ‘moving average’ here refers to the error structure and is not the same concept as ‘average sales over the past seven days’ in feature engineering.

In multi-step forecasting, ARIMA progressively generates future values according to its estimated temporal dependencies, thus not requiring future true sales as lagged inputs. Autoregressive terms use previously generated forecasts for farther horizons, unobserved future innovation terms are treated according to their conditional expectation, and forecast uncertainty propagates accordingly. ARIMA typically models based on equally spaced time series; store closure days, missing dates, and target value masking in retail data need to be handled first according to a unified calendar and target scope.

If the model also requires business inputs like price, discount, or holidays, a regression model with ARIMA errors can be used: the regression part expresses the contribution of external variables, while the ARIMA structure describes the temporal correlation that still exists after regression [24]. The values of these external variables in the future forecast window still need to be prepared in advance, and historical backtesting should also use the planned or external forecast versions visible at each forecast start point.

When residual checks in Chapter 6 reveal stable temporal correlation, ARIMA-like methods can be added as candidates, and their incremental value tested using the same target variable, forecast start point, future business inputs, and evaluation window. If residuals primarily indicate missed holidays, out-of-stock constraints, or business plan discrepancies, then the corresponding data or structure should be corrected first, before deciding whether to incorporate ARIMA errors.

How to Establish a Unified Protocol for Comparing Candidate Models #

After discussing intermittent demand methods, additional business variables, sales lag, and ARIMA, it is necessary to standardize the experimental protocol when incorporating these methods into the candidate set. This helps differentiate whether error differences stem from the model structure or input information.

While standardizing the information boundary, it is not required that all models use exactly the same feature representation. Prophet can express periodicity through Fourier terms, tree models can use date features like day of the week and month, and neural networks can read historical windows. However, all these inputs must be constructed from the same set of information available at the time of forecasting. When a candidate model incorporates additional business variables, historical windows, or cross-series information, changes to the information set should be reported separately to avoid attributing all gains from new information solely to the algorithm.

If future inputs are reliable, the number of variables is limited, and relationships are suitable for linear expression, Prophet with regression terms can be considered. For a large number of categorical attributes, threshold effects, and non-linear interactions, tree models are more worthy of comparison. When multiple related series can share patterns and historical windows have stable value, corresponding global models or neural network architectures can be further validated.

Model Family Overview: Tree Models and Neural Networks #

Regression trees split samples based on rules like “is discount above a certain value” or “is it a weekend,” and provide a forecast value at each leaf node. Gradient-Boosted Decision Trees (GBDT) combine multiple trees, improving the training objective round by round. Both LightGBM and XGBoost can express piecewise non-linear relationships and interactions: the former uses mechanisms like histograms to improve efficiency, while the latter provides a regularized boosting tree framework [25, 26]. However, actual error, speed, and memory footprint still need to be compared within the target data and operating environment.

Neural network forecasting models, on the other hand, encompass a variety of different architectures. Temporal Fusion Transformer (TFT) is geared towards multi-step forecasting, distinguishing between static variables, known future inputs, and inputs visible only in history [27]; N-BEATS’ original task is univariate point forecasting, primarily relying on historical target windows [28]. Therefore, when evaluating neural networks, the specific architecture, input window, and forecasting method should be clearly defined, rather than treating “neural network” as a single candidate.

Model Family Comparison and Engineering Considerations #

The table below compares model structures and engineering requirements, but this does not constitute a fixed performance ranking:

Model families and training scope are two choices: Both tree models and neural networks can employ local or global training. Cross-series sharing depends on sample and structural design; new products and categories still require separate validation for cold start scenarios.

The teaching experiment in the next section standardized the target variable, time splitting, available inputs, and forecast post-processing, focusing on reporting forecast errors and feature behavior. Therefore, the current conclusions are limited to forecasting performance and feature usage. Overall costs in a production environment (including training time, inference costs, and batch failure rates) will need to be quantified separately.

What the Comparison of Prophet and LightGBM on the Same Data Shows #

Experiment Setup and Overall Results #

This experiment compares Prophet, LightGBM, and the seasonal naïve method [29] on a teaching sequence. This method is used as a baseline: for each day of the week, it reads the most recent valid sales value prior to the forecast origin and uses it as the forecast for the same day of the week in the future. This aims to test whether complex models truly outperform simple methods that merely perpetuate weekly patterns.

Let \(T\) be the training cutoff date, \(h\) be the forecast horizon, and \(\operatorname{dow}(t)\) denote the day of the week for date \(t\). The seasonal naïve forecast in this paper is written as:

Here, \(\tau(T+h)\) is the most recent valid date before the forecast origin that shares the same day of the week as the target date. If the training data has no missing or masked target values, this definition is equivalent to the standard seasonal naïve method with a 7-day period; we look forward for valid observations to skip unavailable target values caused by out-of-stock situations or store closures. This baseline requires at least one valid observation for each day of the week in the training history, which is met by the teaching data; if real data lacks valid history for a certain day of the week, a mean-based baseline or a business fallback value would need to be defined separately.

LightGBM uses date features, VIC public holidays, and discount plans known at the time of forecasting, but no sales lags or rolling statistics; forecasts for the next 28 days are generated directly based on known features for each target date. The model is first validated chronologically within the training history, with the number of trees selected via early stopping, and then refitted using all valid training history. All three methods use the same target variable, evaluation dates, out-of-stock filtering, business calendar, and non-negative post-processing.

Under the current teaching data and feature set, Prophet achieved lower WAPE in both the final evaluation window and rolling backtesting. Both Prophet and LightGBM’s WAPE were significantly lower than the seasonal naïve baseline. The trends, seasonality, holidays, and discount mechanisms in the teaching data are relatively close to Prophet’s model structure, so this set of results suggests that Prophet matches the current data generating mechanism well; its relative performance on real retail data still needs to be re-validated on corresponding business samples.

The final evaluation window includes 28 calendar days, of which 27 were operating days with no out-of-stock. The 12 rolling origins collectively formed 323 “forecast origin—target date” records that were operating days and not out-of-stock; since the same target date might be forecasted repeatedly from different origins, 323 should be understood as the number of forecast records, not the number of mutually independent calendar dates. Summarized WAPE is defined as the sum of absolute errors for all valid records divided by the sum of actual sales. All three methods use the same business calendar and non-negative rules, with sales for closed stores fixed at zero.

Prophet uses the teaching configuration determined previously; LightGBM fixes the learning rate, number of leaves, and regularization parameters, only selecting the number of trees via internal time validation and early stopping. Subsequent holiday diagnostics will additionally compare three pre-specified candidate settings. Therefore, the ranking in the table applies to the current limited configuration range; a more complete model ranking would require designing comparable hyperparameter searches for each candidate. XGBoost and neural networks are reserved for future experiments; the table only reports methods actually run in this instance.

Why LightGBM Did Not Use Holiday Features #

Overall results show LightGBM’s error is slightly higher than Prophet’s, and a grouped check further localized part of the gap to holiday records. In the training data corresponding to the final evaluation window, the original LightGBM adopted a design that can be described in English as named-holiday indicator features, i.e., establishing indicator variables separately for each holiday name; this paper abbreviates it as “separate holiday encoding.” To illustrate this design, the table below uses readable holiday names for conceptual column names:

The first two columns represent “separate holiday encoding”: Boxing Day only activates its corresponding column, and Australia Day only activates its corresponding column. The last column represents the “unified public holiday flag” compared later, which only distinguishes whether a day is a public holiday, no longer distinguishing specific names. To avoid spaces and special characters in names, the actual code uses column names like holiday_00, holiday_01, etc., and separately saves the correspondence between numbers and holiday names.

The original design generated 15 holiday features in total. In the final training data, each feature had only 0-3 valid positive instances, and the final fitted model showed zero splits and gain for these features. min_child_samples=20 requires leaf nodes to obtain more training samples, therefore, these sparse indicator variables are difficult to form effective splits independently; default feature pre-filtering might further limit them from entering the split search [30].

The 12 rolling origins collectively formed 8 public holiday “forecast origin—target date” records where stores were open. In this small sample, LightGBM and Prophet’s WAPE were approximately 17.83% and 3.65% respectively; for regular dates, WAPE was approximately 5.31% and 4.47% respectively. This set of results is used to localize the concentration of errors, and the sample size is not yet sufficient to support stable holiday performance ranking. Indicator variables corresponding to fixed closed holidays like Good Friday and Christmas Day had no positive instances in the valid training samples because target values for these dates were masked; the model cannot estimate corresponding sales effects from these training samples, while zero sales after closure are handled by business calendar rules.

How Results Change After Merging Holiday Features #

The above diagnosis indicates that after being separated by name, each holiday feature received very few effective samples. To examine the impact of feature sparsity, this experiment compared three LightGBM holiday feature settings using two historical forecast origins before the final 28-day evaluation window, with each origin evaluating the next 28 days:

By lowering the minimum leaf sample threshold, separately defined holiday indicator variables gained more opportunities to participate in splits. The unified public holiday flag, on the other hand, consolidated multiple sparse events into a single feature: the teaching data set the same additive effect of \(28H_t\) for all public holidays, thus, the unified flag not only increased the number of positive instances but also more closely approximated the data generating mechanism, achieving lower WAPE in internal validation. Since target values for closed days were masked, this unified flag primarily obtained parameter information from open public holidays in the valid training samples. If different holidays in real business have different directions or magnitudes, the unified flag might obscure these differences, and the scope of merging would still need to be determined based on business mechanisms and historical sample sizes.

Two historical origins are sufficient to demonstrate the process of “first selecting configuration within training history, then evaluating the final window,” but candidate ranking may still be affected by specific origins; in real projects, the number of historical origins should be expanded, and stability across different stores and products should be checked.

Candidate selection was based entirely on these two historical origins, with the final 28-day evaluation window reserved for a single independent evaluation. The selected unified public holiday flag configuration had a WAPE of 4.36% in that window, lower than the original separate holiday encoding configuration of 4.68%, and narrowed the gap with Prophet’s approximately 4.22%.

Feature importance in the figure is summarized by split gain on the training objective, used to describe how the model uses features; causal contributions require analysis through corresponding causal designs. Sparse holiday features missing from the figure can be cross-checked in the full gain table to confirm whether they participated in splits.

This experiment shows that model comparison depends not only on the algorithm family but also on how features represent business mechanisms and how many valid samples each feature receives. Prophet benefited from its explicit seasonality and holiday structure on the current teaching data; the unified public holiday flag reduced LightGBM’s WAPE in the final evaluation window from 4.68% to 4.36%, indicating that part of the gap in the original results stemmed from holiday feature representation and sample support. Real projects still need to be re-compared on representative stores, products, and independent evaluation windows.

How to Position Prophet in a Candidate Model System #

In a candidate model system, Prophet is first and foremost a time series model with a clear structure that is easy to review, and can also serve as an interpretable baseline for other methods. Engineers can explicitly configure trends, seasonality, events, and planned variables, generate forecasts directly based on future dates, and examine the contribution of each component to the results. For sequences with sufficient valid history, relatively stable patterns, and where business teams need to review the forecast structure, Prophet is worth including in formal backtesting.

Different data mechanisms correspond to different candidate directions. The table below is used to determine the priority of experiments and does not represent a fixed performance ranking across all datasets:

This chapter has run Prophet, LightGBM, and the seasonal naïve baseline; for XGBoost, ARIMA, neural networks, and intermittent demand methods, the main text explains their modeling mechanisms and candidate conditions, with performance ranking reserved for future experiments. Regardless of which model path a sequence enters, simple baselines matching the task should be retained, such as seasonal naïve, category baselines, or business rules.

All candidates follow the unified comparison protocol from the previous section. The final selection will integrate rolling backtesting, critical date and product group errors from Chapter 6, as well as training cost, delivery success rate, and rollback conditions discussed in the next chapter; when candidate models use additional business variables, historical windows, or cross-series information, changes in the information set should also be reported separately.

When target-granularity data is sufficient, future inputs are available at launch, and unified backtesting shows Prophet’s error is competitive and its component explanations have business value, Prophet can be considered a primary candidate. If other methods achieve lower and more stable errors, Prophet can still serve as an interpretable baseline or fallback method; when other paths are more suitable for the current sequence in terms of forecasting effectiveness, running cost, and maintainability, that sequence can be assigned to the corresponding model, and Prophet no longer needs to be maintained separately.

After completing model allocation for store-SKU sequences, the target values, plan snapshots, model versions, forecast results, backtesting, and fallback rules need to be integrated into a traceable delivery process. Chapter 8 will discuss this engineering boundary, and Chapter 9 will then provide implementation recommendations based on it.

8. MLOps: From Candidate Models to a Production-Ready Demand Forecasting System #

Chapter 7 identified candidate models based on demand sequence data. This chapter integrates these models into a production-ready batch forecasting pipeline, first by defining the gap between teaching experiments and production systems, and explaining how demand forecasting fits into replenishment decisions. Next, we will discuss daily data construction, task scalability, parallel execution, model version management, and failure fallback strategies. Finally, this chapter will cover continuous evaluation, release, and handover of responsibilities after models go live.

Retail systems typically generate 28-day forecasts in daily batches, but forecasting frequency, model retraining frequency, and model release times can be configured independently. The system needs to orchestrate these three cadences based on data changes, computational costs, and backtesting results.

Defining System Boundaries and Business Outputs #

What Still Needs Validation Before Moving from a Teaching Experiment to Production #

The accompanying Notebook provides evidence of model behavior and experimental processes on teaching data. To integrate candidate models into a replenishment system, it is also necessary to verify real-world interfaces, production scale, and operational responsibilities, and to define the following evidence boundaries:

Model Training and Release Workflow. Candidate models complete backtesting and quality checks using historical data; only candidates that pass the quality gate will replace the current production version:

flowchart LR
    A["Versioned Training Data"] --> B["Candidate Model Training"]
    B --> C["Temporal Backtesting & Group Evaluation"]
    C --> D{"Quality Gate Passed?"}
    D -->|Yes| E["Release New Model Version"]
    D -->|No| F["Retain Current Production Version"]

Daily Batch Forecasting Workflow. Daily tasks read already published models. If input is unqualified, model loading fails, or inference fails, the system executes corresponding cold start or fallback rules, and explicitly logs the status:

flowchart LR
    A["Transaction, Inventory & Plan Snapshots"] --> B["Feature & Quality Checks"]
    B --> C{"Input Qualified?"}
    C -->|Yes| D["Published Model Forecasting"]
    C -->|No| E["Cold Start or Data Fallback"]
    D --> F{"Inference Successful?"}
    F -->|Yes| G["Write to Versioned Forecast Table"]
    F -->|No| H["Model Operation Fallback"]
    E --> G
    H --> G
    G --> I["Replenishment Calculation & Anomaly Review"]

How to Convert Demand Forecasts into Replenishment Recommendations #

Once demand forecasts enter replenishment decision-making, the system first needs to determine how long the current inventory needs to cover. This timeframe is called the protection period, typically determined by the replenishment lead time and the review period between inventory checks. Following previous notation, let \(T\) be the forecast origin, \(h\) be the forecast horizon, and \(\hat y_{T+h\mid T}\) denote the point forecast provided by the model at time \(T\) for the \(h\)-th future day. Summing the daily forecasts within the protection period yields the cumulative demand forecast for that period:

Where \(H_{\mathrm{prot}}\) is the number of forecast steps included in the protection period, and \(\hat y_{\mathrm{prot}\mid T}\) is the corresponding cumulative demand forecast. In the teaching example, Prophet’s yhat output has been converted back to SKU business units, representing daily demand point forecasts. If replenishment only needs to cover actual operating days, then forecast_units should first be derived based on the known business calendar, and then the protection period total should be calculated.

The replenishment system adds safety stock to the cumulative demand forecast to form a target inventory level, then subtracts the current inventory position; if existing inventory is sufficient to cover the target, the recommended order quantity is zero. This simplified rule is written as:

Notation:

Retail Case Studies. For comparison, let’s assume inventory position is calculated as ‘current available inventory + ordered in-transit inventory - committed quantity’, and that forecasts, inventory, and order quantities all use the same SKU unit:

In the egg case, the target inventory level is \(420+60=480\) cartons; the current inventory position is 270 cartons, so the preliminary recommended order is 210 cartons. In the wine case, the target inventory level is 30 bottles, which is lower than the current inventory position of 38 bottles, thus the recommended order quantity is zero. The promotional beverage case yields a preliminary recommendation of 390 cans; if the supplier requires ordering by case of 24 cans, this would need to be rounded up to 17 cases, or 408 cans, followed by checks for minimum order quantity, storage capacity, and shelf life.

These three examples also illustrate that demand forecasting is just one input to replenishment calculations. Quantities in the formula must use the same SKU units, and the scope of inventory and in-transit quantities must match the protection period. \(Q_T\) is a preliminary recommendation before applying business constraints, and actual orders must also consider case packs, minimum order quantities, shelf life, supplier restrictions, and storage capacity.

Safety stock should be determined based on the error distribution of cumulative demand over the protection period and the stockout risk acceptable to the business. Forecast errors for different days may be correlated, so this correlation needs to be preserved during calculation. If the system already uses demand quantiles to set target inventory, it is also necessary to check whether the same risk buffer is already included to avoid redundantly adding safety stock.

Building a Scalable Batch Forecasting Pipeline #

A batch forecasting pipeline first converts transaction records into store-SKU daily time series suitable for modeling. It then estimates the training and prediction scale based on the number of valid series. The system executes model training and prediction using complete series as task units, saving the model, input conditions, prediction results, and failure status. This ensures that every run is traceable and reproducible.

Aggregating Transaction Records into Store-SKU Daily Series #

The data pipeline first aggregates transaction records by (store_id, sku_id, business_date) in a data warehouse or distributed engine to obtain the daily observed sales volume SUM(sold_units). We assume business_date is generated according to the store’s local operating day; sold_units represents the quantity sold as recorded by the system, which may not necessarily equal unconstrained potential demand. Returns need to be reconciled separately: if the goal is customer purchase demand, directly subtracting returns as negative sales might alter the business meaning of the target variable.

Transaction aggregation results typically only include dates when transactions occurred. The system also needs to join product availability periods, store business calendars, inventory, price, and promotion snapshots at the same granularity, and build a complete date index for each valid store-SKU series. Missing transaction records can only be imputed as zero sales when it’s confirmed that the product was on sale, the store was open, and transaction data is complete for that day; situations like store closure, product not listed, or interface data missing should retain their respective business statuses. Before joining dimension tables, it is also necessary to check if the keys are unique to prevent sales from being amplified by many-to-many joins.

A single store-SKU series’ three-year history only contains approximately 1,095 daily observations. The distributed engine is responsible for cleaning, aggregating, joining, and date imputation, while fit() for each series is executed as an independent task.

How to Estimate the Task Scale for Local Models #

Standard Prophet is a local forecasting model, meaning each series is fitted independently; placing multiple store-SKU series in the same table does not automatically share parameters or representations. In contrast, global forecasting models train multiple series jointly, allowing them to share some parameters or representations. Here, “local/global” describes the training scope across series, not local versus global optima.

With 1,000 stores and 10,000 candidate products each, the theoretical number of combinations can reach ten million. The production system will first exclude combinations that are not on sale, have been delisted, or do not meet forecasting conditions. Then, it calculates the actual number of series based on valid on-sale combinations for the target period:

Where \(N_{\mathrm{series}}\) is the number of valid store-SKU combinations requiring forecasting, and \(H\) is the forecast horizon (number of steps) for each series. When all series generate complete \(H\) target date forecasts, \(N_{\mathrm{outputs}}\) is the number of result rows produced by a single batch forecast. Different channels and fulfillment nodes may adopt different criteria for valid combinations. Therefore, the production task scale should be calculated based on the actual forecast targets for the target period; valid combinations also dynamically adjust with product listing, delisting, and changes in store status.

The system should first establish valid on-sale combinations, then compute series profiles such as historical length, number of non-zero events, Average Demand Interval (ADI), and squared coefficient of variation of non-zero demand. Series profiles need to specify both time and organizational granularity: the same product might be sparse at the single-store daily granularity but stable at the regional weekly granularity. The system follows the routing rules established in Chapter 7, assigning Prophet, intermittent demand baselines, global models, or cold start strategies to each series; this chapter only discusses how these tasks are executed in batches. Specific definitions and time scopes for ADI and demand variability can be found in Appendix A.

Parallel Training Using Complete Series as Units #

For series suitable for Prophet, the (store_id, sku_id) pair typically serves as the task unit, with the full training window handed to one execution process or node (worker). Each task retains the complete historical data for a single series, allowing the model to estimate trends, seasonality, and event effects along a continuous timeline; the scheduling system then distributes parallel tasks among different series.

Batch training can utilize process pools, task queues, or distributed frameworks. When Spark performs grouped training, each worker needs to have the same dependencies and be able to accommodate the complete series; varying series lengths can also lead to uneven processing times.

The total fitting workload for batch backtesting depends on the number of series participating in backtesting, the number of configurations to compare per series, the number of backtesting start points used per configuration set, and the measured average time per fit. Multiplying these four yields the cumulative time consumption for all fitting tasks, not the actual waiting time after parallel execution:

Notation:

Here, \(N_{\mathrm{fit}}\) is a subset of \(N_{\mathrm{series}}\). Taking Prophet batch backtesting as an example, it only counts series routed to the Prophet candidate path after model diversion, not those assigned to intermittent demand methods, global models, or cold start rules.

Assuming \(N_{\mathrm{fit}}=100{,}000\) series participating in backtesting for this run, 4 configurations, 5 backtesting start points, and an average fitting time of 2 seconds, the total workload would be approximately 4 million seconds, or about 1,111 hours. Under ideal conditions with 200 compute slots at 100% utilization, the computation time would be about 5.6 hours; actual runtime would also need to consider the long-tail effect of tasks, scheduling overhead, prediction process, final model retraining, data reading, and failure retries. The 2 seconds here is merely a capacity estimate; actual planning should use measured results from the target environment.

The system can periodically filter grouped configuration sets on representative series and reuse validated configurations during planned retraining; daily predictions would then read the currently published model. For high-value series, those with significant drift, or deteriorating forecast performance, a more comprehensive parameter search would be triggered.

Prophet’s diagnostic interface supports process, thread, and Dask parallel backtesting [15]. This addresses the execution method for backtesting tasks. Production environments still need to orchestrate training between series separately and control nested parallelism to prevent each worker from launching another set of processes and exhausting resources.

Saving Models and Their Prediction Conditions #

Prophet models contain Stan backend objects, and it is officially recommended to use JSON format for model serialization and restoration, rather than directly using Python pickle [22].

The accompanying implementation uses model_to_json() / model_from_json() to save and load models, and checks if the loaded point predictions are consistent with those before saving; randomly generated interval endpoints are not used for deterministic consistency checks.

Model files can only restore the fitted parameters. To reproduce a prediction result, the system also needs to record the training data snapshot, training cutoff date, prediction input snapshot, feature definitions, plan version, parameter configuration, code version, dependency version, and backtesting results.

The forecast table is recommended to save the following fields:

The system can use (store_id, sku_id, forecast_origin, ds, run_id) to identify an immutable forecast record, ensuring idempotent writes when retrying the same run. The system retains the immutable results of each run, then publishes forecasts for downstream consumption via a “currently active version” pointer; when some tasks fail, they fall back according to rules, and the source of each forecast is recorded.

How to Handle Failures and Generate Traceable Results #

Batch jobs should isolate exceptions for each series, ensuring that a single model failure only affects its corresponding task. Engineers can define fallback strategies in advance, such as using the last valid and feature-compatible model with the current plan for prediction, seasonal naive forecasting, or category rules. Each strategy should be backtested separately, and anomaly status and fallback forecasts should be saved together.

The system should also distinguish between normal forecasts, unqualified input data, model unavailability, inference failure, successful fallback, and product discontinuation, among other statuses. True zero demand should be retained as a numerical zero; records for which predictions cannot be generated should be retained as null values, with the failure reason and fallback source written. This prevents downstream replenishment systems from interpreting operational anomalies as zero demand.

Continuous Evaluation, Deployment, and Handover of Responsibilities #

After a model is deployed, the system must continuously record actual predictions and evaluate the performance of different forecast horizons once target values mature. Monitoring results can trigger data investigation, retraining of candidate models, and controlled releases, but each stage requires independent quality gates and responsible parties. This section will sequentially explain online evaluation, CI/CT/CD processes, and system handover requirements.

How to Continuously Record and Evaluate Forecasts Post-Deployment #

Post-deployment monitoring needs to observe four dimensions: operational status, data quality, forecast performance, and business results:

Stockout rate, inventory turnover, and shrinkage are also influenced by replenishment strategies, supply execution, and store operations, so these metrics are used to evaluate the complete decision system and are not solely attributable to the forecast model. Forecast error must be evaluated using actual values consistent with the modeling objective. For example, evaluating only in-stock dates might overlook stockout risks; conversely, using constrained sales directly to assess potential demand forecasts might misjudge a reasonably high forecast as an error. Monitoring reports should clearly define the scope of the target variable, valid observation conditions, and target value coverage rate.

Continuous evaluation should verify the immutable forecasts actually issued at the time, rather than recalculating historical forecasts using a new model. The system can link forecasts with actuals by store_id, sku_id, forecast_origin, ds, and run_id, and retain model_version. Here, forecast_origin distinguishes forecasts issued from different origins for the same target date, run_id identifies the specific delivery batch, and model_version records the production model used at that time. Offline backtesting simulates historical points in time, while online evaluation checks actual delivered results; they serve different purposes.

After the target date passes, sales records may still change due to returns, interface re-uploads, or inventory status revisions. The system should define target value maturity rules: only when records reach a specified waiting time, complete necessary revisions, and pass quality checks, can they enter formal performance evaluation. Evaluation data can store actual_available_at, actual_version, stockout_status, and evaluation_status to explain when target values are available, which revision version is used, and whether they can participate in evaluation.

For the same target date, forecasts issued 1 day, 7 days, and 28 days in advance might coexist. For example, the demand for a certain store-SKU on December 20, 2025, could come from three different batch forecasts:

The target date for all three forecasts is the same, but each can only use the model version published at that time and the available business plans. If plan retraining occurs during this period, the training history cutoff date will also change with the model version, potentially leading to three different forecast values. After the target value for December 20 matures, the system will compare the same actual value with each of the three forecasts, thereby evaluating the model’s performance at 1-day, 7-day, and 28-day forecast horizons. The system should retain each forecast and summarize them separately by exact forecast horizon. When rolling windows overlap, the same actual value might enter multiple forecast-actual pairs, so reports need to provide the number of forecast records, target dates, and forecast origins simultaneously.

The statistical window must also be clearly defined as by forecast origin or by target date. Target values that have not yet arrived or matured should remain missing. Metrics from different periods are comparable only when the comparison windows reach the same level of maturity. Metrics can further be grouped by forecast horizon, category, store, sales volume, and promotion status to check whether overall results mask systemic underestimation in important scenarios.

Drift monitoring needs to distinguish three types of changes: distribution changes in input features like price belong to feature drift; distribution changes in sales targets belong to target drift; and changes in the relationship between inputs and targets belong to concept drift. Changes in input and target distributions can provide early warnings; after target values mature, the team can then combine forecast error, bias, and residual structure to judge whether the model relationship has changed. When metrics suddenly improve, it is also important to check whether stockout days have been excluded, whether the product mix has changed, and whether the target value coverage rate has decreased.

Retraining rules should simultaneously specify error thresholds, minimum effective sample sizes, consecutive trigger windows, and data quality conditions. For example, when a certain category shows significant underestimation in multiple consecutive windows with sufficient valid samples, the team should first investigate plan coverage, stockout status, and business changes, then decide whether to initiate candidate model retraining.

How CI, CT, and CD Connect Testing, Training, and Deployment #

Machine Learning Operations (MLOps) needs to distinguish between code and configuration validation, candidate model training, model deployment, and daily batch forecasting. Continuous Integration (CI), Continuous Training (CT), and Continuous Delivery (CD) are responsible for the first three stages, respectively [31]:

Daily batch forecasting reads the currently active model version and is not part of CT, nor does it imply daily retraining. CT is responsible for generating and evaluating candidate versions, and CD decides which version goes into production; therefore, forecasting, evaluation, retraining, and deployment can adopt different frequencies.

Continuous evaluation provides trigger signals for CT. Monitoring alerts should first lead to data and plan input investigation, followed by comparing candidate models against the current version and baseline methods under the same information conditions and evaluation windows. Only after a candidate version passes quality gates can CD deploy it; otherwise, the current qualified version continues to be used or established rollback rules are executed.

For example, the system can generate forecasts daily, initiate candidate model training once a week, and only deploy a new version after the candidate version passes checks for overall metrics, key product groupings, and underestimation risk. When performance continuously deteriorates, additional training can be initiated, but the same evaluation and deployment thresholds should still apply.

What Responsibilities Need to Be Defined During System Handover #

When delivering a system, the team needs to clearly define responsibilities for data quality, model approval, batch operations, and replenishment decisions:

These items describe areas of responsibility and do not require a company to establish four separate positions; the same team can undertake multiple responsibilities. Handover materials should cover data contracts and target variable definitions, model configurations and evaluation reports, operation schedules and deployment standards, alerting and rollback rules, as well as disaster recovery and historical forecast reproduction steps. The receiving team should be able to reproduce forecasts based on saved operational versions and practice recovery processes for data loss, batch failure, or candidate deployment failure.

When continuous monitoring discovers data anomalies, the problem will be returned to the data pipeline for processing; when model error continuously deteriorates on valid samples, the system will initiate retraining or re-compare candidate models; if business objectives, product scope, or replenishment rules change, the team will need to reconfirm forecast targets and acceptance criteria. The purpose of responsibility handover is to ensure that each type of problem has clear evidence, processing procedures, and responsible parties.

9. Summary and Implementation Recommendations #

Core Conclusions #

This series of articles began with a 28-day daily demand forecasting task for a store-SKU, and successively discussed forecast targets, Prophet model structure, data transformation, parameter solving, temporal backtesting, candidate model comparison, and production delivery. The core conclusion is: whether Prophet is suitable for a retail demand series does not depend on the model name, but on whether the target variable is clear, historical information is sufficient, temporal patterns can be identified, future business variables can be obtained, and whether the model can consistently outperform simple benchmarks in temporal backtesting.

Recommended Implementation Sequence #

Real projects can gradually expand their scope according to the following sequence:

Complete Python code, accompanying notebooks, and running instructions are available in the GitHub repository: retail-demand-forecast. The tutorial code is used to reproduce data generation, model fitting, component reconstruction, and backtesting results; when migrating to real business, revalidation is still required based on enterprise data definitions, available information, and business constraints.

Production implementation should follow the chain of evidence established in this article: save truly available data and plan snapshots at the time of forecasting, compare candidate models during temporal backtesting, validate data, forecasting, and fallback processes through small-scale pilot runs, and then integrate versions that meet quality gates into the replenishment system. After deployment, the team should continuously evaluate forecasts using mature and consistently defined target values, and retain model versions, operational status, and business accountability. Prophet can serve as a primary model, an interpretable baseline, or a fallback method; its final position in the system should be jointly determined by data conditions, out-of-sample performance, business value, running costs, and maintenance complexity.

Appendix: Supplementary Demand Mechanisms and Candidate Methods #

This appendix will unfold in the following order: first, defining key metrics; next, comparing teaching scenarios; and finally, introducing candidate methods. Appendix A defines the Average Demand Interval (ADI) and the Coefficient of Variation of non-zero demand quantities (CV+^2). Appendix B compares four scenarios—groceries, red wine, garden supplies, and luxury handbags—using a unified standard. Appendix C introduces candidate methods applicable to intermittent demand. All four sets of teaching data use the same date range, facilitating the comparison of demand characteristics under a consistent framework.

All four datasets follow the same teaching business calendar, setting observed sales for Good Friday and Christmas Day to zero, but retaining the potential demand assuming the stores were operating normally. The y in the figures represents observed sales after applying business status and corresponding supply rules, which is not equivalent to the potential demand in the formulas. Real projects require replacing these teaching rules with actual store operating records.

These four scenarios adopt the data design from the teaching examples in the main text. The calendar date for each record is stored in ds, and \(t\) represents the number of days calculated from the sample start. The generation formula first produces potential demand based on dates and business variables, then applies business and supply rules to obtain observed sales \(y_t\), which are written into the y field. Here, \(t\) is used to express the generation pattern of the teaching data and is not a raw field that needs to be separately provided to Prophet.

Appendix A: Describing Demand Series with ADI and Demand Quantity Variation #

A system can describe a series using the Average Demand Interval (ADI) and the squared coefficient of variation of non-zero demand quantities [32]:

First, identify the set of valid observations from daily records sorted by ds. This paper excludes closed days and retains dates where businesses are operating normally and target values are available. \(T\) is the number of these valid dates, and \(N_+\) is the count of dates among them where \(y_t>0\); \(\bar y_+\) and \(s_+^2\) are the sample mean and sample variance, respectively, of all positive sales \(y_t\).

ADI describes how frequently demand events occur. When \(\widehat{\mathrm{ADI}}=1\), there is, on average, one positive demand event per valid observation period; when \(\widehat{\mathrm{ADI}}=2\), positive demand occurs, on average, approximately once every two observation periods. A larger ADI typically indicates longer intervals between non-zero events, making the series sparser.

\(\widehat{\mathrm{CV}}_+^2\) describes how unstable positive demand quantities are once demand occurs. The coefficient of variation is defined as the standard deviation divided by the mean, and squaring it yields \(\widehat{\mathrm{CV}}_+^2\); since it is unitless, it allows for comparison between series with different average sales volumes. A value close to zero indicates that positive demand quantities are relatively concentrated, while an increasing value indicates stronger fluctuations in positive demand quantities relative to their own mean.

These two metrics describe different dimensions: ADI focuses on the frequency and intervals of positive demand occurrences, while \(\widehat{\mathrm{CV}}_+^2\) focuses solely on the magnitude variation of positive demand quantities. They describe the target series y after processing business statuses, not simulated potential demand; if stockout days are still retained in y, the metrics will also be affected by constrained sales. When selecting a model, both metrics need to be considered simultaneously, in conjunction with factors such as historical data length, temporal patterns, business status, and backtesting.

When calculating these metrics, it is also necessary to clarify whether to use continuous calendar time or operational time that only includes normal business hours and valid sales opportunities, and to maintain consistency in definition across training, backtesting, and production monitoring. Closed days and out-of-stock items can be excluded from operational time based on business definitions; however, truncated demand due to stockouts and cases of missing interface data should retain their original dates and unknown status. If there are no non-zero demand events, these two statistics need to be specially flagged and cannot be calculated directly using the standard formulas.

Appendix B: Data Profiles of Four Teaching Scenarios #

Unified Comparison of Four Scenarios #

The table below compares the results generated from four sets of teaching data using a fixed random seed. For all scenarios, dates with is_open=True are considered effective operating days, and sequence profiles are calculated using observed sales y. Here, \(N_+\) denotes the number of operating days where \(y>0\); the definitions of \(\widehat{\mathrm{ADI}}\) and \(\widehat{\mathrm{CV}}_+^2\) are provided in Appendix A.

The first three scenarios primarily exhibit frequent non-zero sales, not forming distinct intermittent demand. The wine scenario has a small number of operating days with zero sales, so its ADI is slightly above 1; the grocery and yard supplies scenarios have positive sales on every effective operating day. In contrast, the luxury handbag scenario has zero sales on 93.58% of operating days, with an average interval of about 15.57 operating days between positive sales, and significantly higher relative fluctuation in positive sales. The demand type depends on the specific store-SKU-time granularity and cannot be judged solely by product category or the name of the generative distribution.

By combining \(\widehat{\mathrm{ADI}}\) and \(\widehat{\mathrm{CV}}_+^2\), the four scenarios reveal different demand profiles: the former distinguishes the frequency of positive demand occurrences, while the latter distinguishes the stability of purchase quantities once demand occurs. The same method can be extended to hundreds or thousands of SKUs in a store, initially grouping them based on their historical sequences, and then designing common candidate models, feature templates, and parameter search ranges for SKUs with similar characteristics. Such grouping can reduce the workload of repeatedly designing and tuning parameters for individual sequences but cannot directly replace time-series backtesting; the final model and configuration still need to be determined based on the out-of-sample performance of each group and important SKUs.

Grocery Teaching Case: High-Frequency Demand and Out-of-Stock Constraints #

The main text uses the grocery scenario to illustrate Prophet’s model definition, configuration, solving, and evaluation. This sequence has positive sales on every effective operating day and is primarily used to demonstrate how trends, weekly patterns, yearly patterns, holidays, promotions, and out-of-stock situations jointly affect observed sales, thus representing data conditions where Prophet can easily establish an effective baseline.

The current data includes 24 out-of-stock days. The average observed sales on operating days and \(\widehat{\mathrm{CV}}_+^2\) in the table are derived from observed sales y, thus retaining constrained records caused by inventory limits, and do not equate to complete potential demand. The data generation formula and time series plots are presented in Chapter 2 and will not be repeated here.

When fitting Prophet, ds and the y processed according to the business rules in the main text constitute the basic input; discounts can be used as additional regression variables when future plans are known, while holidays and business status are used to construct event features and sales opportunity constraints, respectively. The noise-free demand mean and potential demand are used only to verify the data generation mechanism for teaching purposes.

Wine Scenario: Count Series with Frequent Non-Zero Demand #

For a hypothetical product of a single 750 mL bottle of wine, the daily Poisson intensity is first synthesized, and then the number of potential demand bottles is sampled:

The current scenario does not simulate out-of-stock situations. Potential demand is passed through the business calendar to obtain observed sales:

Notation Explanation:

The \(10^{-6}\) in the formula is a tiny lower bound set to keep the Poisson intensity positive; \(\max\) indicates taking the larger value, and \(\sim\) means “follows the… distribution”.

The formula first calculates \(\lambda_t\) based on the weekday, date, and activity schedule corresponding to ds, then samples potential demand \(D_t^{\mathrm{wine}}\), and finally obtains observed sales \(y_t\) via \(O_t\). Therefore, y on operating days equals the potential demand generated in this scenario, and y on closed days is zero.

When fitting Prophet, ds and the processed y constitute the basic input. Weekly and annual patterns are constructed by Prophet based on ds; party events and carton promotions can be configured as holiday features or additional regression variables only if they can be predetermined in the forecast horizon. Poisson intensity \(\lambda_t\) and potential demand \(D_t^{\mathrm{wine}}\) are simulation ground truths for diagnostic purposes and are not used as model features.

The figure above shows the observed sales along with closed and out-of-stock markers throughout the complete history; the figure below magnifies March 18 to April 7, 2023, and highlights the carton promotion windows. The magnified view is used to verify that promotion markers, sales changes, and supply status are correctly aligned by date. It displays the results of synthetic data generation and does not independently prove that promotions have causal effects of the same magnitude on sales.

Data Characteristics and Comparison Suggestions

The Poisson distribution is used for non-negative integer counts but does not automatically imply intermittent demand. Given the intensity \(\lambda_t\) for a particular day, the probability of a Poisson variable taking a value of zero is:

The lower the intensity, the higher the probability of zero demand; for example, at \(\lambda_t=3\), it’s about 5.0%, and at \(\lambda_t=0.2\), it’s about 81.9%. Whether it constitutes intermittent demand also needs to be judged in conjunction with the actual zero value proportion, number of non-zero events, demand interval, and fluctuation of non-zero demand volume at the chosen granularity.

The unified profile table shows that the current wine scenario has 1,090 effective operating days, with positive sales on 1,076 of them, a zero value proportion of about 1.28% on operating days, \(\widehat{\mathrm{ADI}}\approx1.013\), and \(\widehat{\mathrm{CV}}_+^2\approx0.387\). Although it uses a Poisson observation distribution, it is still a count series with frequent non-zero demand and does not form distinct intermittent demand.

Readers can compare Prophet with event inputs against Poisson regression and check for negative prediction rates, event grouping errors, and interval coverage. To study intermittent demand, one can reduce demand intensity or change the demand generation mechanism, then re-examine the zero-value rate, ADI, and non-zero demand volume fluctuations. The demand type should be determined by the data characteristics at the target granularity, rather than by brand or distribution name.

Garden Consumables Scenario: Continuous Demand with Strong Seasonality #

This scenario models a garden consumable, such as leaf collection bags or fall overseeding grass seeds. Business demand is synthesized and then subjected to independent random noise:

This scenario does not simulate stockouts, so observed sales are:

Notation Explanation:

The formula first generates trend, annual cycle, weekend, and event effects based on the ds date, then adds random noise to derive the latent demand \(D_t^{\mathrm{garden}}\). For operating days, y equals this latent demand; for closed days, \(O_t\) sets y to zero.

When fitting Prophet, ds and the processed y form the core input. Weekends and annual cycles can be constructed from ds; Labour Day, Easter, and clearance events must be converted into holiday features or business variables available during the forecast period. The noise-free demand level \(\mu_t^{\mathrm{garden}}\), random noise \(\eta_t\), and latent demand \(D_t^{\mathrm{garden}}\) are used solely for simulation and diagnosis.

The upper graph displays annual variations, weekly fluctuations, and closed dates throughout the full history. The lower graph magnifies May 18 to June 7, 2023, highlighting the May 25–31 autumn clearance window. This magnified view is used to verify the date alignment of clearance markers with the sales series; sales are still simultaneously influenced by trend, seasonality, weekend effects, and random noise.

Experimental Suggestions

The unified profile table indicates that positive sales occurred on all 1,090 valid operating days in the garden consumables scenario, with an average observed daily sales volume of approximately 64.92 units. \(\widehat{\mathrm{ADI}}=1.000\) and \(\widehat{\mathrm{CV}}_+^2\approx0.235\). This scenario is designed for studying strong seasonality, weekend effects, and event-driven changes in continuous demand, not for intermittent demand cases.

The binary clearance marker corresponds to an additional 12 units, while weather effects would require constructing additional inputs. Readers can compare additive versus multiplicative models, examining off-season errors and non-negative post-processing; floor=0 constrains the trend, but total forecasts and lower interval bounds are still influenced by other components [4, 5].

The four teaching scenarios can be compared using a common forecast cutoff date, multiple-origin backtesting, and several random seeds. The configuration choice uses only the training history, reserving the final evaluation window for a single use after development is complete; real-world projects should also set aside an independent test period that did not participate in development. When interpreting results, WAPE, the proportion of zero-sales operating days, and the number of observed events should be reported concurrently, with model ranking determined by experimental outcomes.

Luxury Handbag Scenario: Intermittent High-Value Demand #

This scenario models a high-end handbag at a Melbourne CBD store, separating the simulation of “whether positive demand occurs on a given day” from “how many units are purchased if demand occurs.” Let \(R_t\) denote whether positive demand occurs on day \(t\), and \(p_t\) its corresponding probability:

When \(R_t=1\), the positive demand quantity \(Z_t\) is drawn from the following discrete distribution:

Notation Explanation:

The theoretical mean and squared coefficient of variation for the positive demand quantity distribution are:

Thus, the conditional mean of latent demand on day \(t\) is \(\operatorname E(D_t^{\mathrm{luxury}}\mid p_t)=1.65p_t\), corresponding to conditional_mean in the data. The teaching store is closed on Good Friday and Christmas Day, and this scenario does not simulate stockouts, so observed sales are:

The final modeling table writes the date of the \(t\)-th record to ds and the observed sales \(y_t\) to y. Weekends and the December peak season can be constructed from ds; VIP events and known discounts can only be configured as additional Prophet inputs if they are already scheduled for the forecast period. occurrence_probability, positive_demand, conditional_mean, and latent_demand record diagnostic ground truth from the synthetic data generation process, which is typically not directly observable in real business and should not be used as model features.

The upper graph illustrates sparse positive sales and VIP event dates throughout the full history. The lower graph magnifies December 1–21, 2023, simultaneously showing sales and the demand occurrence probability from the generation process. Weekends, the December peak season, and VIP events increase \(p_t\), but \(R_t\) is still randomly generated by a Bernoulli distribution, so a higher occurrence probability does not guarantee sales on a given day. This figure is for verifying the generation mechanism and date alignment and does not imply that these factors have the same magnitude of causal effect in a real store.

Data Characteristics and Comparison Suggestions

With the current fixed random seed, 1,090 valid operating days are generated, with positive sales occurring on 70 days. The proportion of zero-sales operating days is approximately 93.58%, with \(\widehat{\mathrm{ADI}}\approx15.571\) and \(\widehat{\mathrm{CV}}_+^2\approx0.949\). The sample mean for positive demand is 1.30 units, showing a sampling difference from the theoretical mean of 1.65 units; theoretical quantities describe the generating distribution, while sample quantities describe this specific generated finite sequence, and the two should not be conflated.

This scenario, characterized by long intervals between positive demands and highly volatile positive demand quantities, is suitable for comparing candidate methods such as seasonal naive benchmarks, Prophet, Croston, SBA, TSB, and hurdle models. Prophet can leverage dates, VIP events, and known discounts to explain the average demand rate, but its continuous curve structure may not be ideal for directly modeling numerous zero values and discrete purchase counts. Experiments should use a consistent forecast cutoff date, future available inputs, and rolling-origin backtesting, separately checking metrics related to zero-sales days, positive demand days, event windows, and inventory decisions; the final selection should be determined by out-of-sample results.

Appendix C: Candidate Methods for Intermittent Demand #

The luxury handbag scenario illustrates that when daily demand is largely zero, models need to simultaneously consider the frequency of positive demand occurrence and the quantity purchased when it does occur. This appendix introduces the intermittent demand methods included in the candidate set in Chapter 7. Such demand is also common in scenarios like industrial spare parts [33], but the specific method should still be chosen based on the data conditions of the target retail series. The main text uses these methods to define the applicability boundaries of Prophet, without pre-ranking their forecasting performance.

This section uses \(Y_t\) to denote the random variable for sales on date \(t\), and \(y_t\) to denote the specific observed sales in the data table, corresponding to Prophet’s y. When comparing different methods, all candidate models use the same ds date range, target value processing rules, forecast cutoff date, and forecast window; each method then uses these time-ordered observations according to its own structure.

Croston and SBA: Smoothing Demand Quantity and Demand Interval Separately #

The Croston method only performs exponential smoothing on non-zero demand quantity \(z_i\) and the interval \(\ell_i\) between two non-zero demands when non-zero demand occurs [16]:

Here, \(z_i\) is the \(i\)-th positive sales in the target column y, \(\ell_i\) is the non-zero demand interval calculated based on valid dates ordered by ds, and \(\alpha\in(0,1]\) is the smoothing parameter. \(\hat y\) represents the average demand rate per unit time, not the time of the next demand occurrence. SBA (Syntetos–Boylan Approximation) adds a bias correction to the same decomposition [17]:

This correction comes from an analysis of the bias in Croston’s estimation under specific assumptions. Actual selection should also be validated in conjunction with product data, evaluation metrics, and inventory targets.

TSB: Updating Demand Occurrence Probability in Every Period #

The TSB (Teunter–Syntetos–Babai) method updates the demand occurrence probability in each period instead. For the target column y ordered by ds, let \(o_t=\mathbf 1(y_t>0)\), smooth the occurrence probability with \(\alpha_p\); non-zero demand quantity uses \(\alpha_z\), and is only updated when \(y_t>0\) [18]:

Where \(\alpha_p,\alpha_z\in(0,1]\). Consecutive zero values cause TSB’s demand occurrence probability to decrease incrementally, thus quickly reflecting fading demand. Product listing, de-listing, and operational status can be used as additional business information in conjunction with this statistical update mechanism.

Zero-Inflated Models and Hurdle Models: Explicitly Modeling Abundant Zero Values #

Zero-inflated models handle zero values using a mixture distribution. For example, a Zero-Inflated Poisson (ZIP) model uses \(\pi_t\) to represent the probability of an additional zero component, and \(\lambda_t\) to represent the Poisson count mean [19]:

Here, \(a!\) is the factorial of the positive integer demand quantity \(a\). Count distributions themselves can also produce zero values; ‘structural zeros’ represent a latent mixture component in the model; business status flags can help determine its relationship with unlisted items, store closures, or missed sales opportunities. Hurdle models (also known as two-part models) adopt another decomposition: the first part predicts whether demand is greater than zero, and the second part fits a zero-truncated distribution only for positive demand quantities. Both types of models link demand occurrence with demand quantity but employ different zero-value generation assumptions.

Method Differences and Experimental Boundaries #

This table compares model mechanisms and does not represent a predetermined performance ranking. Actual experiments still require unifying target variables, time granularity, forecast windows, and available information, and using a rolling-origin backtest consistent with the main text to compare Prophet, intermittent demand methods, and other candidate models.

References #

[1] Taylor, S. J., & Letham, B. Forecasting at Scale. The American Statistician, 72(1), 37–45, 2018. https://doi.org/10.1080/00031305.2017.1380080. Author preprint: https://facebook.github.io/prophet/static/prophet_paper_20170113.pdf

[2] Business Victoria. Operating on a restricted trading day. https://business.vic.gov.au/business-information/public-holidays/operating-on-a-restricted-trading-day

[3] Prophet. Trend Changepoints. https://facebook.github.io/prophet/docs/trend_changepoints.html

[4] Prophet 1.4.0. Date preprocessing and forecast components. forecaster.py. https://github.com/facebook/prophet/blob/v1.4.0/python/prophet/forecaster.py

[5] Prophet. Saturating Forecasts. https://facebook.github.io/prophet/docs/saturating_forecasts.html

[6] Prophet. Seasonality, Holiday Effects, and Regressors. https://facebook.github.io/prophet/docs/seasonality,_holiday_effects,_and_regressors.html

[7] Prophet. Multiplicative Seasonality. https://facebook.github.io/prophet/docs/multiplicative_seasonality.html

[8] Prophet. Installation. https://facebook.github.io/prophet/docs/installation.html

[9] Prophet. Quick Start. https://facebook.github.io/prophet/docs/quick_start.html

[10] Prophet 1.4.0. Stan model source. prophet.stan. https://github.com/facebook/prophet/blob/v1.4.0/python/stan/prophet.stan

[11] Stan. CmdStan User’s Guide: Optimization. https://mc-stan.org/docs/cmdstan-guide/optimize_config.html

[12] Prophet 1.4.0. Python Stan backend source. models.py. https://github.com/facebook/prophet/blob/v1.4.0/python/prophet/models.py

[13] Prophet. Uncertainty Intervals. https://facebook.github.io/prophet/docs/uncertainty_intervals.html

[14] Stan. Reference Manual: MCMC Sampling. https://mc-stan.org/docs/reference-manual/mcmc.html

[15] Prophet. Diagnostics. https://facebook.github.io/prophet/docs/diagnostics.html

[16] Croston, J. D. Forecasting and Stock Control for Intermittent Demands. Operational Research Quarterly, 23(3), 289–303, 1972. https://doi.org/10.1057/jors.1972.50

[17] Syntetos, A. A., & Boylan, J. E. The accuracy of intermittent demand estimates. International Journal of Forecasting, 21(2), 303–314, 2005. https://doi.org/10.1016/j.ijforecast.2004.10.001

[18] Teunter, R. H., Syntetos, A. A., & Babai, M. Z. Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214(3), 606–615, 2011. https://doi.org/10.1016/j.ejor.2011.05.018

[19] Lambert, D. Zero-Inflated Poisson Regression, with an Application to Defects in Manufacturing. Technometrics, 34(1), 1–14, 1992. https://doi.org/10.1080/00401706.1992.10485228

[20] Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice, 3rd edition. Chapter 11: Forecasting hierarchical and grouped time series. OTexts. https://otexts.com/fpp3/hierarchical.html

[21] Prophet. Non-Daily Data. https://facebook.github.io/prophet/docs/non-daily_data.html

[22] Prophet. Additional Topics. https://facebook.github.io/prophet/docs/additional_topics.html

[23] Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice, 3rd edition. Section 9.5: Non-seasonal ARIMA models. OTexts. https://otexts.com/fpp3/non-seasonal-arima.html

[24] Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice, 3rd edition. Section 10.2: Regression with ARIMA errors using fable. OTexts. https://otexts.com/fpp3/regarima.html

[25] LightGBM. Features. https://lightgbm.readthedocs.io/en/v4.6.0/Features.html

[26] XGBoost. Introduction to Boosted Trees. https://xgboost.readthedocs.io/en/stable/tutorials/model.html

[27] Lim, B., Arik, S. O., Loeff, N., & Pfister, T. Temporal Fusion Transformers for Interpretable Multi-horizon Time Series Forecasting. arXiv:1912.09363, 2019; revised 2020. https://arxiv.org/abs/1912.09363

[28] Oreshkin, B. N., Carpov, D., Chapados, N., & Bengio, Y. N-BEATS: Neural basis expansion analysis for interpretable time series forecasting. arXiv:1905.10437, 2019; revised 2020. https://arxiv.org/abs/1905.10437

[29] Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice, 2nd edition. Section 3.1: Seasonal naïve method. OTexts. https://otexts.com/fpp2/simple-methods.html

[30] LightGBM 4.6.0. Parameters: min_data_in_leaf and feature_pre_filter. https://lightgbm.readthedocs.io/en/v4.6.0/Parameters.html

[31] Google Cloud. MLOps: Continuous delivery and automation pipelines in machine learning. https://cloud.google.com/architecture/mlops-continuous-delivery-and-automation-pipelines-in-machine-learning

[32] Syntetos, A. A., Boylan, J. E., & Croston, J. D. On the categorization of demand patterns. Journal of the Operational Research Society, 56(5), 495–503, 2005. https://doi.org/10.1057/palgrave.jors.2601841

[33] Pennings, C. L. P., van Dalen, J., & van der Laan, E. A. Exploiting elapsed time for managing intermittent demand for spare parts. European Journal of Operational Research, 258(3), 958–969, 2017. https://doi.org/10.1016/j.ejor.2016.09.017