Demand Forecasting with Statistical Models: How Far Can Demand History Alone Take Us?僅憑歷史需求，我們能預測到多準？

This is the statistical deep dive behind Forecasting iPhone 17 Pro Demand: From Statistics to Deep Learning. The broader article compares three generations of forecasting models. Here, we slow down and ask a narrower question:

How much of tomorrow's demand can we forecast from demand history alone?

The answer establishes a baseline that every more complex model must beat.

This article is structured as follows:

1. The decision behind the forecast
2. One series, several behaviours
3. Forward-walking backtesting
4. Traditional statistical time-series forecasting
   • Holt–Winters
   • ARIMA
   • SARIMA
5. How do we choose parameters?
   • STL
   • ADF
   • ACF and PACF
   • AIC and BIC
6. Evaluation metrics and model selection
7. From point forecast to a 95% prediction interval
8. Recommendation for inventory decision

The decision behind the forecast

Every evening, the inventory team in a synthetic London Central Apple Store scenario must decide how many units of one configuration—iPhone 17 Pro 256GB Light—should be available the following day.

Forecast too low and the store risks stockouts, lost sales, and disappointed customers. Forecast too high and expensive devices occupy working capital and storage capacity. At one store, the difference may look modest. Across many stores and configurations, small daily errors accumulate into a material inventory decision.

The modelling target is therefore:

and the operational question is:

How many units should be positioned before tomorrow begins?

A demand forecast does not automatically tell us how many phones to stock. That decision also depends on the cost of having too few phones compared with the cost of having too many. If stockouts are more expensive, the planner should choose a higher point in the forecast distribution. If excess inventory is more expensive, the planner should choose a lower point. The newsvendor model expresses this choice as:

Here, $C_u$ is the cost of understocking, $C_o$ is the cost of overstocking, $F^{-1}$ selects the corresponding quantile from the forecast demand distribution, and $q^*$ is the recommended inventory quantity. Because these costs are not available in this experiment, we cannot calculate an optimal inventory quantile. We therefore report the point forecast and 95% prediction interval as forecasting outputs, leaving the service-level decision to the business.

One series, several behaviours

The synthetic dataset contains 730 simulated daily observations from 19 September 2025 to 18 September 2027. Mean demand is 39.5 units per day. The model receives only the date and historical unit demand for this store-SKU combination.

Daily observations are shown in grey; the seven-day rolling mean exposes the slower movement in demand.

The chart contains several distinct forms of temporal structure:

• Lifecycle: demand begins high after launch and declines toward a more mature level.
• Weekly seasonality: the day of the week changes the expected demand.
• Calendar effects: Christmas creates sharp local peaks, while closure days can produce genuine zero demand.
• Level changes: the local baseline rises and falls rather than remaining constant across two years.
• Irregular shocks: promotions and other commercial events produce spikes that are not explained by trend or weekday alone.

This is a univariate forecasting problem: the history of $y_t$ is the entire information set. Price, promotions, marketing, holidays, weather, and competitor activity are deliberately withheld. Their effects appear only after they reach the series as unexplained demand shocks.

That constraint is useful. Before giving a machine-learning model more features, we first want to know what the demand series can explain by itself.

Forward-walking backtesting

A random train-test split would leak future market conditions into the past. Instead, the observations are divided chronologically:

The evaluation then reproduces the nightly workflow:

1. Fit the model using every observation available at the end of day $t$.
2. Forecast demand for day $t+1$.
3. Reveal the actual demand for $t+1$.
4. Add that observation to the available history.
5. Move forward one day and repeat.

This is variously called walk-forward validation, rolling-origin evaluation, or time-series backtesting. The forecast horizon remains one day, while the forecast origin advances through the validation and test windows.

At each step, the models use only observations available at that point in time. Candidate specifications are selected on the validation window and then locked before evaluation on the test window. This prevents the test data from influencing model design and keeps the reported performance realistic.

Traditional statistical time-series forecasting

Both Holt–Winters and ARIMA assume that the series' own history contains useful information about its future. They differ in what they believe that information looks like.

Holt–Winters: forecast the evolving components

Holt–Winters belongs to the exponential-smoothing family. It represents a time series through interpretable states:

For an additive model with seasonal period $m$, the state updates are

and the $h$-step-ahead forecast is

Each state has a business interpretation:

• Level $\ell_t$: the current baseline after removing seasonal effects.
• Trend $b_t$: the estimated daily change in that baseline.
• Seasonal state $s_t$: the recurring effect for a position in the cycle; with $m=7$, the model maintains a weekday pattern.

Suppose it is Thursday evening and the model is forecasting demand for Saturday, two days ahead. It currently estimates:

• baseline demand: 35 units;
• trend: −0.2 units per day;
• Saturday seasonal effect: +8 units.

The forecast is:

The model therefore forecasts approximately 43 units.

When Saturday's actual demand becomes available, Holt–Winters updates all three estimates. The level parameter $\alpha$ determines how strongly the baseline reacts, $\beta$ controls how quickly the trend changes, and $\gamma$ determines how much of the difference is treated as a change in the Saturday effect.

The smoothing parameters $\alpha$, $\beta$, and $\gamma$ lie between zero and one. A high value reacts strongly to new observations; a low value preserves more of the previous state. In other words, recent observations matter most, but older observations fade gradually rather than disappearing at a fixed cutoff.

Holt–Winters works best when observations are equally spaced, the seasonal period is known and stable, level and trend evolve gradually, and unusual events do not dominate the series. Additive seasonality assumes a weekday effect of roughly constant magnitude. Multiplicative seasonality is excluded because the series contains genuine zero-demand closure days, while the implementation used here requires strictly positive observations.

We validate three weekly additive candidates:

• level plus seasonality;
• level, additive trend, and seasonality;
• level, damped trend, and seasonality.

The component specification is selected on the validation window; within each specification, the smoothing coefficients are estimated numerically from the training history rather than chosen by hand.

The candidates are close, but additive trend has the lowest validation WAPE and is locked as the family winner.

A damped trend replaces the indefinitely extrapolated $h b_t$ term with a trend whose incremental effect decays over the horizon. This often prevents implausible long-range extrapolation, although it is less consequential for a one-day forecast.

ARIMA: model changes, persistence, and surprises

Holt–Winters describes demand through level, trend, and seasonality. ARIMA approaches the same problem differently. It asks:

After making the series sufficiently stable, how does current demand depend on previous values and previous forecasting errors?

The process can be understood as:

Original demand
      ↓ differencing
Approximately stationary series
      ↓ AR and MA modelling
Forecast temporal dependence
      ↓ reverse differencing
Forecast demand in its original units

ARIMA stands for AutoRegressive Integrated Moving Average. Its notation is $\mathrm{ARIMA}(p,d,q)$, where each parameter represents a different source of information.

Integrated (I): model the change instead of the level

ARIMA generally works best when the series it models is approximately stationary. Its mean, variance, and dependence structure should not change substantially over time.

If demand has a changing level or trend, we can take the first difference:

Instead of modelling demand directly, ARIMA then models the daily change. Suppose demand increases from 40 units on Monday to 44 units on Tuesday:

The transformed observation is an increase of four units. If one ordinary difference is applied, then $d=1$. If the original series is already sufficiently stationary, then $d=0$.

Differencing does not mean that the final forecast is expressed only as a change. After forecasting $w_{t+1}$, ARIMA adds the predicted change back to the latest observed demand. For example:

Latest demand:             44 units
Forecasted daily change:   +3 units
Forecasted demand:         47 units

This reversal returns the prediction to the business quantity we care about: units demanded.

Autoregression (AR): use previous values

The autoregressive component assumes that recent values of the stationary or differenced series may contain information about the next value.

An AR(1) relationship is:

Here, $w_{t-1}$ is the previous transformed value, $\phi_1$ measures how strongly it carries forward, $c$ is a constant, and $\varepsilon_t$ is a new, unpredictable shock.

Suppose yesterday's demand increased by five units and the model estimates $\phi_1=0.6$. The autoregressive contribution to today's expected change is:

The model is saying that part of yesterday's movement is likely to continue today. The order $p$ controls how many previous values are included: $p=1$ uses one previous value, while $p=2$ uses the previous two.

Moving average (MA): use previous forecast errors

The moving-average component uses recent forecast errors:

The previous error is:

Suppose yesterday's forecast was 40 units but actual demand was 45:

The model underestimated demand by five units. If recent surprises tend to carry information into the following day, the MA component uses part of this error to adjust the next forecast.

The order $q$ determines how many previous errors are used. Despite its name, this is not an ordinary moving average of previous demand; it is a weighted combination of previous forecasting errors.

Combining the components

An ARIMA$(1,1,1)$ model first calculates the daily difference $w_t=y_t-y_{t-1}$. It then models that difference using one previous value and one previous forecast error:

Finally, it reverses the difference to produce a forecast in units demanded. In plain language, the model asks:

Given the most recent change in demand and how wrong the previous forecast was, what change should we expect next?

The general ARIMA equation summarises the same idea for any orders $p$ and $q$:

Here, $w_t=\Delta^d y_t$ is the series after applying $d$ differences. The first sum combines the previous $p$ transformed values, while the second sum combines the previous $q$ forecast errors. This article follows the plus-sign convention used by statsmodels; some textbooks define the MA coefficients with the opposite sign.

SARIMA: extend the same idea to weekly demand

ARIMA models short-term relationships between recent values and forecast errors. SARIMA extends it to recurring seasonal patterns:

The lowercase parameters describe ordinary ARIMA behaviour. The uppercase parameters describe seasonal behaviour:

• $P$: seasonal autoregressive order;
• $D$: seasonal differencing order;
• $Q$: seasonal moving-average order;
• $m$: number of observations in one seasonal cycle.

For daily demand with a weekly pattern, $m=7$.

Seasonal differencing

A weekly seasonal difference compares demand with the same weekday one week earlier:

Suppose demand is 48 units this Monday and was 42 units the previous Monday:

Seasonal naive assumes this difference will be zero:

SARIMA instead predicts the weekly change and adds it to last week's demand:

In simple terms:

Seasonal naive repeats last week. SARIMA starts from last week and adjusts for recent changes and forecast errors.

The selected specification

For a Monday forecast, the model begins with demand from the previous Monday. It estimates how demand has changed relative to last week, adjusts for recent forecasting errors, and adds the predicted change back to last Monday's demand.

The validation process described next selects:

• $p=1$: use one previous value of the transformed series;
• $d=0$: no ordinary daily differencing;
• $q=1$: use one recent forecast error.
• $P=0$: no seasonal autoregressive term;
• $D=1$: apply weekly seasonal differencing;
• $Q=1$: use a forecast error from seven days earlier;
• $m=7$: use a seven-day seasonal cycle.

How do we choose parameters?

STL, ADF, ACF, PACF, AIC, and BIC help us identify reasonable statistical model specifications. They answer different questions and narrow the candidate set, but they do not select the final model automatically. Chronological validation performance and residual diagnostics make that decision.

STL — Seasonal-Trend decomposition using LOESS

STL explores whether the time series contains trend, seasonality, changing seasonal strength, or unusual events. LOESS means Locally Estimated Scatterplot Smoothing. Seasonal-Trend decomposition using LOESS writes the training series as

where $T_t$ is a smooth trend, $S_t$ is the repeating seasonal component, and $R_t$ is the remainder.

STL is a diagnostic, not a forecasting model. In this plot it provides four useful clues:

1. The trend changes substantially after launch and around holiday periods, so a fixed historical mean is inadequate.
2. A repeating seven-day component supports $m=7$ for both Holt–Winters and SARIMA.
3. Seasonal amplitude varies, making multiplicative seasonality worth considering in theory, but zero-demand days make it invalid here.
4. Large remainders coincide with events that trend and weekday effects cannot explain. A history-only model will continue to miss some of them.

Standard STL, Holt–Winters, and SARIMA each handle one seasonal period cleanly. If both weekly and annual seasonalities needed explicit representation, we would need an extension such as MSTL, TBATS, or dynamic harmonic regression.

ADF — Augmented Dickey–Fuller test

ADF helps assess whether ordinary differencing may be necessary. It tests whether a series contains a unit root, a common form of non-stationarity.

A weakly stationary series has an approximately constant mean, constant variance, and an autocovariance that depends on lag rather than calendar time. ARIMA needs this behaviour in the series it actually models.

The Augmented Dickey–Fuller test uses:

• $H_0$: the series has a unit root and is non-stationary;
• $H_1$: the series does not have a unit root under the test specification.

A common rule is:

The results are:

Because the raw-demand p-value is below 0.05, we reject the unit-root null. That supports including candidates with $d=0$: ordinary first differencing is not automatically required.

It does not prove that the raw series has no trend, seasonality, changing variance, or structural breaks. ADF addresses a specific unit-root question. Weekly differencing is still reasonable because the weekly pattern and correlation plots show a strong lag-seven relationship. We therefore test seasonal orders with both $D=0$ and $D=1$.

ACF (Autocorrelation Function) and PACF (Partial Autocorrelation Function)

ACF identifies persistent temporal relationships, seasonal cycles, and possible moving-average structure.

The autocorrelation function measures the correlation between the series and its previous values:

at each lag $k$. The partial autocorrelation at lag $k$ measures the incremental relationship after controlling for shorter lags.

For example:

• lag 1 compares today with yesterday;
• lag 7 compares today with the same weekday last week;
• lag 14 compares today with the same weekday two weeks ago.

Typical interpretations include:

• slow decay across many lags: differencing may be needed;
• spikes at 7, 14, and 21: weekly seasonality;
• a sharp cutoff after lag $q$: possible MA order $q$.

The PACF measures the direct relationship between $y_t$ and $y_{t-k}$ after controlling for the shorter lags. For example, PACF at lag 7 asks:

Does demand from seven days ago add information after accounting for lags 1 through 6?

Typical interpretations include:

• a strong lag-1 PACF: consider $p=1$;
• strong lags 1 and 2: consider $p=2$;
• a seasonal spike at lag 7: consider a seasonal AR term $P=1$.

The raw ACF is strongly positive at short lags and shows recurring strength at weekly multiples. The PACF after lag-seven differencing retains a clear lag-one relationship and a smaller set of notable seasonal-multiple spikes. Together, these patterns motivate:

• a weekly period $m=7$;
• non-seasonal AR orders $p=1$ and $p=2$ as compact alternatives;
• a non-seasonal MA term $q=1$;
• testing either a seasonal AR term $P=1$ or a seasonally differenced MA structure with $D=1$, $Q=1$.

ACF and PACF heuristics are cleanest for pure AR or pure MA processes. Real series with mixed AR, MA, seasonal terms, shocks, and finite samples rarely produce textbook cutoffs. They suggest a candidate set; they do not read the final orders directly from the chart.

AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion)

AIC and BIC compare fitted candidate specifications while penalising unnecessary complexity.

AIC compares statistical models by balancing likelihood fit with model complexity:

$AIC=-2\log L+2k$

where $L$ is the likelihood and $k$ is the number of estimated parameters.

BIC uses a stronger complexity penalty:

$BIC=-2\log L+k\log n$

where $n$ is the number of observations.

AIC places relatively more emphasis on predictive fit, whereas BIC's stronger penalty generally favours simpler models. Both are in-sample specification diagnostics; neither directly measures operational forecast error.

The tools play complementary roles:

The resulting candidate set is deliberately small:

Lower is preferred among models fitted to the same target, observations, transformation, and comparable likelihood. They are valuable specification diagnostics, but they do not directly measure one-day-ahead operational error. Because differencing can also change the effective likelihood comparison, the chronological validation WAPE remains the primary selection criterion.

Evaluation metrics and model selection

No single metric describes every kind of forecast failure.

MAE: the typical miss in business units

MAE answers the most tangible question: by how many phones is the forecast wrong on an average day? Every unit of error receives equal weight.

RMSE: emphasise large misses

Squaring errors makes RMSE more sensitive to occasional large misses. This is useful when any large forecast miss is disproportionately costly. RMSE penalises overforecasting and underforecasting symmetrically; bias is still needed to describe direction. Comparing RMSE with MAE also reveals whether errors are fairly consistent or dominated by spikes.

WAPE: error relative to total demand

WAPE is the primary selection metric because it is scale-normalised, aggregates naturally across the evaluation window, and remains defined when individual days have zero demand. MAPE would divide by each actual value and therefore break on closure days.

A WAPE of 25% means total absolute error equals roughly 25% of total observed demand. It does not mean the model is “75% accurate.” Whether 25% is useful depends on the cost and tolerance of the inventory process and on performance relative to credible baselines.

Bias: the direction of total error

Negative bias means systematic underforecasting; positive bias means overforecasting. Bias cannot stand alone because positive and negative errors can cancel. A nearly unbiased model may still be highly inaccurate.

Validation selects the family winners

On the validation window, additive-trend Holt–Winters wins its family with 28.51% WAPE. The selected SARIMA specification wins its family with 27.59%. Those choices are then locked before the test window is opened.

Three simple baselines provide context. Historical mean predicts the expanding average of all demand observed so far, naive repeats yesterday's demand, and seasonal naive repeats demand from seven days earlier.

The final test results are:

Holt–Winters lowers WAPE by about 32% relative to seasonal naive and misses by about 8.9 units on an average day. Its aggregate bias is almost zero. SARIMA is only 0.28 WAPE points behind, so the result does not justify a general claim that Holt–Winters is the superior family. It wins this particular test window.

The forecast paths also show why the statistical models win: they adapt more smoothly than seasonal naive while retaining the weekly movement. None of them anticipates every sharp commercial shock.

Residual diagnostics add an important qualification. The Ljung–Box test asks whether a group of residual autocorrelations is jointly different from zero. At lags 7, 14, and 21, Holt–Winters has p-values near 0.005, evidence that temporal structure remains in its residuals. SARIMA's p-values range from 0.31 to 0.89, providing no strong evidence of remaining autocorrelation at those lags.

SARIMA is statistically cleaner; Holt–Winters is marginally more accurate on the test window. These findings are compatible. Ljung–Box asks whether residual dependence remains. WAPE asks how large the point errors were. Operational accuracy and statistical adequacy are related, but they are not the same criterion.

From point forecast to a 95% prediction interval

After evaluation, both selected models are refitted using all 730 observations. Their one-day-ahead forecasts for the next simulated date, 19 September 2027, are:

The point forecasts are nearly identical. The test result makes Holt–Winters' 40 units the statistical recommendation, but the interval is just as important as the centre.

SARIMA's state-space likelihood produces a forecast distribution and a model-based interval. Conceptually,

with $z_{0.975}\approx1.96$, subject to the model's residual assumptions. The Holt–Winters implementation used in the notebook does not expose an analytic interval, so its interval is an explicit Gaussian residual approximation:

Lower bounds are clipped at zero because negative unit demand is not operationally meaningful.

Across the 56 test days, the approximate Holt–Winters interval covers 100% of actuals with a mean width of 59.2 units. SARIMA covers 98.2% with a mean width of 57.7 units. Both coverage rates exceed the nominal 95%, suggesting that the intervals may be conservative. With only 56 observations—and an approximate Holt–Winters interval—this remains a diagnostic rather than proof of calibration.

The bands are also wide relative to mean daily demand of 39.5 units. That is an honest result. A history-only model can estimate the recurring centre of demand, but it cannot know tomorrow's promotion, competitor launch, weather, or other external shock.

These intervals quantify uncertainty; they do not prescribe an inventory target. The final step is to combine the forecast with the current inventory position, replenishment lead time, and the business cost of shortages and surplus stock.

Recommendation for inventory decision

For the next simulated day, use Holt–Winters' 40-unit forecast as the expected-demand baseline—not as the order quantity. Holt–Winters becomes the initial production model, while SARIMA continues running in parallel as a challenger.

The replenishment decision should account for the inventory already available, deliveries confirmed to arrive in time, and any safety buffer required by the chosen service level:

If expected demand is 40 units, the store has 27 units available, five more will arrive before the next selling day, and no safety buffer has yet been added, then:

The immediate replenishment requirement is eight units. This calculation assumes that replenishment can arrive before tomorrow's demand occurs. If the lead time is longer, the inventory target must cover forecast demand across the full lead-time period.

The upper end of the 95% prediction interval—approximately 70 units—should not automatically become the stock target. It describes forecast uncertainty, not the service level the business should choose. The safety buffer should instead reflect the cost of a stockout relative to the cost of carrying surplus stock.

Known promotions, launches, closures, and price changes should be recorded as explicit business adjustments because neither statistical model sees these events in advance. In production, compare Holt–Winters and SARIMA using rolling WAPE and bias alongside stockouts and excess inventory. Promote SARIMA only if it delivers a sustained improvement across these operational measures.