What it actually does
Forecast with a weighted average of the past where recent observations count for more, and the weights fall away geometrically. One parameter controls how fast the memory fades: near one and the model chases every wiggle, near zero and it barely moves. That is simple exponential smoothing, which handles a level. Add a second component that tracks how fast the level is changing and you have a trend. Add a third that tracks the repeating pattern and you have Holt-Winters. Each component is updated by the same rule — blend what you just observed with what you already believed — so the entire model is three numbers of state and three parameters, and it never looks at more than the last observation.
- The series decomposes into level, trend and a single seasonal cycle — a structural commitment made up front rather than learned, and the source of both the method’s strength and its limits
- Recent observations are more informative than distant ones, with importance decaying geometrically
- The seasonal period is known and fixed, and there is only one of them
- Sampling is regular with no gaps, since the seasonal state is indexed by position rather than timestamp
- Errors are additive with roughly constant variance, if the ETS prediction intervals are to be believed
Loss function
Choose the smoothing parameters that would have produced the smallest one-step forecast errors over the history. Two things are worth noting. Written this way it looks like a heuristic, and for decades it was treated as one — the ETS state-space formulation later showed that each variant corresponds to a proper statistical model with an actual likelihood, which is what gave the method prediction intervals and automatic model selection rather than a rule of thumb. And every parameter is bounded to [0, 1] with a direct interpretation: it is the fraction of the new observation absorbed into that component, so a fitted alpha of 0.9 tells you the series is close to a random walk without any further analysis.
- the one-step-ahead forecast made before seeing y_t — the model is scored on what it did not know
- level smoothing: how much of the new observation replaces the current level
- trend smoothing: how quickly the estimated slope is allowed to change
- seasonal smoothing: how quickly the seasonal pattern is allowed to change
How it is minimized
Bounded numerical optimization of the smoothing parameters and initial state, with model selection by AIC over the ETS family
- the level: where the series is now, with the seasonal effect removed
- the trend: how fast the level is moving, estimated from how much the level just changed
- the seasonal component for this position in the cycle, updated once per period
- the seasonal period; the seasonal state is m numbers, and the rest of the model is two
Every line has the same shape — a convex blend of what was just observed and what was already believed — which is what makes the whole model O(1) in state and one pass in time. There is no matrix, no lag window, and nothing that grows with the length of the history. Fitting means choosing three bounded parameters plus the initial state, which is a small box-constrained optimization solved by L-BFGS-B or Nelder-Mead in milliseconds. Two decisions matter more than the optimizer. Damping the trend — multiplying the slope by a factor below one at each horizon step — is the single most valuable variant, because an undamped linear trend extrapolated far enough is always wrong and damping is what stops it. And the initial state genuinely matters on short series, which is why modern implementations estimate it jointly with the parameters rather than seeding it from the first few observations.
alpha
Level smoothing, and the parameter that matters most. Near 1 the model is a random walk; near 0 it is a global mean (0.05 to 0.4 on stable series, fitted rather than chosen)
beta
Trend smoothing. Often fits near zero, which is the model saying the slope is stable — and is a useful diagnostic in itself
gamma
Seasonal smoothing. Usually small, because a seasonal pattern that changes quickly is generally a sign of a modelling problem elsewhere
phi (damping)
Trend damping factor. The most valuable single addition to the method, and it should be the default rather than an option (0.80 to 0.98; below 0.8 the trend dies almost immediately)
seasonality type
Additive or multiplicative. Multiplicative when the seasonal swing grows with the level, which on sales data it usually does
error type (ETS)
Additive or multiplicative errors. Does not change the point forecast and does change the prediction intervals, which is the whole reason ETS exists
The optimization is small and well behaved — three bounded parameters — but the surface is not always informative. A common outcome is a fitted alpha at the boundary: alpha = 1 means the model has decided the best forecast is the last observation, which is a random walk and is sometimes correct and sometimes a sign that the series has a structural break the model is chasing. Beta fitting to zero is the analogous message about the trend. The named failure mode is the undamped trend: Holt’s linear method extrapolates a straight line forever, so a series with a temporary upswing produces a forecast that keeps climbing, and at long horizons the error is unbounded. Damping is the fix, and the M-competition evidence for making it the default rather than an option is about as strong as forecasting evidence gets.
One pass over the series per parameter evaluation — O(n) time, O(m) memory for the seasonal state and O(1) for everything else. A fit is tens of evaluations over three parameters, so milliseconds. That cost profile is the point: with constant state per series, a million series is a million cheap independent recursions rather than a million optimizations over a design matrix, which is why this and not ARIMA sits underneath large-scale demand planning.
How, where, and why it is used
Time-Series Forecasting
Primary fitChoose the ETS variant — which components exist, and whether errors and seasonality are additive or multiplicative — then fit the smoothing parameters and initial state by minimizing one-step error, and forecast by projecting the components forward. The seasonal component repeats, the level stays put, and the damped trend decays toward zero, which is what makes the long-horizon forecast approach a sensible constant instead of a straight line.
- Large-scale demand and inventory forecasting, where the number of series makes per-series cost the binding constraint
- Any series with clear level, trend and seasonal structure and no useful exogenous drivers
- The benchmark in forecasting competitions, where it repeatedly finishes at or near the top
- Short series, where there is not enough history to identify ARIMA orders but enough to estimate three parameters
It commits to a structure instead of learning one, and on real business data that constraint helps rather than hurts — there is rarely enough signal to identify an arbitrary lag structure, and pretending otherwise fits noise. It is also fast enough to run everywhere and simple enough that a planner can read the fitted alpha and understand what the model believes. The M-competition results are the honest summary: it is competitive with or better than fitted ARIMA on business data, and the field has largely accepted that generality is not accuracy. What it cannot do is use exogenous information, which is the clear case for SARIMAX or a regression instead.
- Use multiplicative seasonality when the seasonal swing grows with the level, which on sales data is the common case
- Damp the trend by default; an undamped linear extrapolation is wrong at long horizons and the damping parameter is nearly free
- Estimate the initial state jointly with the parameters rather than seeding from the first cycle, which matters on short series
- Aggregate or switch to a Croston-type method for intermittent demand, where most periods are zero and the Gaussian error assumption fails outright
Rolling-origin backtesting with MASE against seasonal-naive. Inspect the fitted parameters as a diagnostic: alpha pinned at 1 means the model has become a random walk, and beta at 0 means it found no trend — both are informative before any accuracy number is read.
- An undamped trend extrapolating a temporary upswing indefinitely, which is the classic failure of Holt’s method
- Multiplicative seasonality on a series that touches zero, where the model divides by a near-zero level
- Fitting parameters on the whole series before splitting, which leaks and is easy to do because the fit feels like preprocessing
- Multiple seasonalities — daily and weekly together — which the classical form cannot represent at all
Anomaly Detection
AdaptedThe EWMA control chart is exponential smoothing used for monitoring, and it predates the forecasting application. Smooth the series, compare each new observation against the smoothed value, and signal when the deviation exceeds control limits derived from the smoothed variance. Because the smoothed statistic has far lower variance than the raw series, the chart detects small sustained shifts that a threshold on individual observations misses entirely.
- Statistical process control, where the EWMA chart is standard for detecting small persistent shifts in a process mean
- Infrastructure and business metric monitoring against a trend-and-season-aware baseline
- Change-point screening, where a run of same-signed deviations indicates a level shift rather than an outlier
It is the right detector for a small sustained shift, which is precisely what a per-observation threshold cannot see: averaging is what makes a shift of half a standard deviation visible, and the geometric weighting is what keeps the average responsive. It is also O(1) state, so it runs on a stream at any volume. Against it: the smoothing that gives sensitivity to sustained shifts costs sensitivity to single large spikes, so it is the complement to a threshold rather than a replacement, and the control limits assume a stable variance that many real metrics do not have.
- Fit the smoothing parameters on a confirmed-clean period, since a contaminated fit sets both the baseline and the control limits
- Model seasonality explicitly rather than smoothing across it, or every Monday reads as an anomaly
- Choose alpha for the shift size you care about — smaller alpha detects smaller sustained shifts more slowly, and that trade is the design of the chart
Average run length to detection for a given shift size, which is the process-control metric and is more informative here than precision@k, alongside the false-alarm rate on clean data at the chosen limits.
- A single large spike absorbed into the smoothed level, shifting the baseline and masking what follows
- Control limits computed from a period whose variance differs from the current one, which mis-sets the whole chart
- Treating a run of same-signed deviations as many independent anomalies when it is one level shift
Optimization
AdaptedIt is the forecasting engine inside large-scale inventory and capacity optimization, and the reason is cost rather than accuracy. Replenishment decisions need a demand distribution per SKU per location, which is millions of forecasts refreshed nightly; a model with constant state and a one-pass fit makes that affordable, and a model requiring a per-series optimization over a design matrix does not.
- Replenishment and safety-stock systems across large retail and distribution catalogues
- Capacity and staffing planning where the forecast has to be refreshed frequently across many resources
- Any planning system whose binding constraint is the number of series rather than the difficulty of any one of them
The interesting content here is that the model was chosen for its cost profile and turned out to be competitive on accuracy as well — which is the opposite of how these decisions usually go, and worth internalizing. Its state is three numbers plus the seasonal cycle, its fit is one pass, and it parallelizes perfectly across series because the series are independent. What it does not supply on its own is the distribution the optimizer needs; the ETS state-space formulation is what provides intervals, and using the point forecast with an assumed variance is the shortcut that quietly breaks the service-level guarantee.
- Forecast over the lead-time window rather than one period, since the relevant distribution is demand until the next replenishment arrives
- Use the ETS interval rather than a point forecast plus an assumed variance, or the achieved service level will not match the target
- Handle intermittent demand separately; a Gaussian interval on a mostly-zero series is meaningless and over-stocks systematically
Realized service level and holding cost under the plan, not forecast error. A change that improves MASE and worsens achieved service level has optimized the wrong quantity, which happens whenever the interval degrades while the mean improves.
- Scaling the one-step forecast variance by the horizon, which is wrong for any model with trend or seasonality
- Ignoring lead time so the distribution covers the wrong window
- Treating a forecast refresh as free when the true cost is the downstream replanning it triggers
Control & Operations
Primary fitThe default forecaster in operational planning at scale. One model per SKU, per location, per resource, refitted on a schedule and consumed by whatever orders, schedules or provisions. Its constant state also makes it natural on a stream: the recursion updates in place as each observation arrives, with no window to maintain.
- Demand planning across retail and distribution catalogues in the millions of series
- Workforce and capacity scheduling against forecast arrival rates
- EWMA-based process control on manufacturing lines
- Streaming metric smoothing where a moving average would need a window buffer and this does not
Operational forecasting is dominated by the number of series rather than the difficulty of any single one, and this is the method whose cost profile matches that — constant state, one-pass fit, perfectly parallel across series, and interpretable enough that a planner can override it with a reason. The M-competition evidence that it is also competitive on accuracy is what settles the argument. Its limits show up wherever the structure it commits to is wrong: multiple seasonalities, promotions and holidays, or any case where an exogenous driver is the real explanation.
- Multiplicative seasonality for sales-like series whose swing scales with volume
- Damped trend as the default across a catalogue, since undamped extrapolation on thousands of series will produce a few spectacular failures
- Route intermittent and slow-moving items to a Croston-type method rather than forcing them through the same model
- Handle promotions outside the model — as an override or an adjustment — because the classical form has nowhere to put a regressor
Rolling-origin MASE at the aggregation level decisions are made at, plus interval coverage at the operating service level. Track the share of series where the fitted alpha sits at a boundary, which is a fleet-level diagnostic that individual accuracy numbers will not surface.
- Promotions learned as seasonality because there is no way to declare them as interventions
- Stockouts treated as low demand, teaching the model to forecast the constraint rather than the demand
- One model per series becoming unmonitorable at catalogue scale, so failures are found by the business rather than by the system
What it costs to run
Milliseconds per series: one pass per parameter evaluation, tens of evaluations over three bounded parameters. Millions of series is an embarrassingly parallel batch job, which is the property that put this method underneath large-scale planning rather than any accuracy argument.
Forecasting is projecting three components forward — O(h) arithmetic with no history required, because the state IS the summary of the history. The model is a handful of numbers per series, which makes storing millions of them trivial.
The state updates online with every observation at no cost, so the level, trend and seasonal components stay current continuously. Only the smoothing parameters need periodic refitting, and they move slowly — which is why production systems refresh state nightly and parameters monthly or less.
- Track the distribution of fitted alpha across the fleet — a shift toward 1 means series are becoming less predictable and the model is degrading toward a random walk
- Watch the share of series where a parameter sits at a boundary, which is a fleet-level signal individual accuracy numbers never surface
- Monitor interval coverage rather than point error, since the intervals degrade first when the error assumption stops holding
- Alert on forecasts that diverge over the horizon, which on an undamped trend is the failure that produces the spectacular outliers in a large catalogue
- The state is the model: level, trend and the full seasonal cycle must be persisted and restored exactly, and reinitializing from recent history silently produces a different forecaster
- Multiplicative seasonality breaks on a series that reaches zero, and a catalogue always contains some — the variant must be chosen per series, not globally
- An undamped trend extrapolates without bound, so across thousands of series a few will produce forecasts that are obviously absurd and need clamping
- Missing periods break the seasonal indexing: the recursion assumes regular sampling, so a gap must be imputed rather than skipped or the whole cycle shifts
- Fitted parameters at a boundary are informative rather than broken, and pipelines that reject them as invalid discard the model’s clearest diagnostic
Pros and cons, in context
Constant state and a one-pass fit
Three numbers plus the seasonal cycle, updated in place, with no history required at forecast time. This is why it sits underneath planning systems with millions of series, and it is a cost argument rather than an accuracy one.
Repeatedly competitive with far more general models
The M-competitions found it at or above fitted ARIMA on real business data. Committing to a structure helps when there is not enough signal to identify an arbitrary one, which on business series is most of the time.
Parameters are interpretable on sight
Each is the fraction of a new observation absorbed into a component, so a fitted alpha of 0.9 says the series is near a random walk without any further analysis. Genuinely rare, and useful for triage across a large fleet.
The ETS formulation supplies proper prediction intervals
What turned a heuristic into a statistical model with a likelihood, AIC-based variant selection, and calibrated intervals. Without it the method is a smoothing rule; with it the quantile-shaped decisions downstream are defensible.
No exogenous regressors
Price, promotion, weather and holidays have nowhere to go, so their effects are either absorbed into the seasonal component or missed. This is the clearest case for SARIMAX or a regression instead, and it is a structural limit rather than a gap in the implementation.
Exactly one seasonality
A series with daily and weekly cycles cannot be expressed. TBATS and structural time-series models exist for this, and both give up the simplicity that motivated the method in the first place.
The undamped trend extrapolates without bound
Holt’s linear method turns a temporary upswing into a forecast that climbs forever. Damping fixes it and should be the default; leaving it off across a large catalogue guarantees a handful of spectacular failures.
Fails on intermittent demand
A mostly-zero series breaks the Gaussian error assumption and the intervals with it, which systematically over-stocks. Croston-type methods exist precisely for this, and the routing decision has to be made per series.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Holt-Winters additive - the three update equations, transcribed.
Every line has the same shape: blend what was just observed with what was
already believed. That is why the whole model is three numbers of state plus
the seasonal cycle, and why it never looks further back than one observation.
"""
def forecast_one_step(level, trend, seasonal, period_index):
"""The forecast made BEFORE seeing y_t - which is what is scored."""
return level + trend + seasonal[period_index]
def fit(series, m, alpha, beta, gamma):
"""m is the seasonal period. Returns the state and the sum of squared
one-step errors, which is the objective the parameters are chosen on."""
# Initialization: level from the first cycle, trend from the difference
# between the first two cycles, seasonality from the deviations. Crude,
# and on a short series it visibly matters - which is why real
# implementations estimate the initial state jointly with the parameters.
level = sum(series[:m]) / m
trend = (sum(series[m : 2 * m]) - sum(series[:m])) / (m * m)
seasonal = [series[i] - level for i in range(m)]
total = 0.0
for t in range(m, len(series)):
index = (t - m) % m
prediction = forecast_one_step(level, trend, seasonal, index)
error = series[t] - prediction
total += error * error
previous_level = level
# level: blend the deseasonalized observation with where the level
# was heading.
level = alpha * (series[t] - seasonal[index]) + (1.0 - alpha) * (level + trend)
# trend: blend how much the level just moved with the previous slope.
trend = beta * (level - previous_level) + (1.0 - beta) * trend
# season: blend this period's deviation with the stored one.
seasonal[index] = gamma * (series[t] - level) + (1.0 - gamma) * seasonal[index]
return level, trend, seasonal, total
def select_parameters(series, m, step=0.1):
"""Grid search over the unit cube. Every parameter is a FRACTION, so the
search space is bounded by construction - which is what makes a grid
viable here and not for most models."""
best = (None, None, None, float("inf"))
value = step
while value < 1.0:
beta_value = step
while beta_value < 1.0:
gamma_value = step
while gamma_value < 1.0:
_, _, _, total = fit(series, m, value, beta_value, gamma_value)
if total < best[3]:
best = (value, beta_value, gamma_value, total)
gamma_value += step
beta_value += step
value += step
return best[0], best[1], best[2]
def forecast(level, trend, seasonal, m, horizon):
"""Project forward. The seasonal component repeats and the trend is
extended linearly - which is exactly the behaviour damping exists to
limit, since an undamped slope is extrapolated forever."""
out = []
for h in range(1, horizon + 1):
out.append(level + h * trend + seasonal[(h - 1) % m])
return outO(n) per parameter evaluation, and a three-dimensional grid at step 0.1 means about 700 passes over the series.
How this language actually wants it written.
"""Holt-Winters - typed, damped, additive or multiplicative, bounded fit."""
from dataclasses import dataclass
from typing import Literal
import numpy as np
from numpy.typing import NDArray
from scipy.optimize import minimize
Vector = NDArray[np.float64]
Seasonality = Literal["additive", "multiplicative"]
@dataclass(frozen=True)
class SmoothingState:
"""The state IS the model: level, trend, and one value per seasonal slot.
Constant in the length of the history, which is the property that makes
this method viable across millions of series.
"""
level: float
trend: float
seasonal: Vector
period: int
@dataclass(frozen=True)
class HoltWinters:
state: SmoothingState
alpha: float
beta: float
gamma: float
phi: float # trend damping
seasonality: Seasonality
sse: float
def forecast(self, horizon: int) -> Vector:
"""Damped trend: the slope contributes phi + phi^2 + ... + phi^h, which
converges rather than growing linearly. An undamped Holt forecast
extrapolates a straight line forever, and at long horizons that is
always wrong."""
steps = np.arange(1, horizon + 1)
damping = np.cumsum(self.phi ** steps)
base = self.state.level + damping * self.state.trend
season = self.state.seasonal[(steps - 1) % self.state.period]
return base + season if self.seasonality == "additive" else base * season
def _run(
series: Vector,
period: int,
alpha: float,
beta: float,
gamma: float,
phi: float,
seasonality: Seasonality,
) -> tuple[SmoothingState, float]:
"""One pass. Sequential in t by construction - each update reads the state
the previous step produced."""
level = float(series[:period].mean())
trend = float((series[period : 2 * period].sum() - series[:period].sum()) / period**2)
seasonal = (
series[:period] - level if seasonality == "additive" else series[:period] / level
).astype(np.float64)
sse = 0.0
for t in range(period, series.size):
index = (t - period) % period
if seasonality == "additive":
prediction = level + phi * trend + seasonal[index]
else:
prediction = (level + phi * trend) * seasonal[index]
error = series[t] - prediction
sse += error * error
previous = level
if seasonality == "additive":
level = alpha * (series[t] - seasonal[index]) + (1 - alpha) * (level + phi * trend)
seasonal[index] = gamma * (series[t] - level) + (1 - gamma) * seasonal[index]
else:
level = alpha * (series[t] / seasonal[index]) + (1 - alpha) * (level + phi * trend)
seasonal[index] = gamma * (series[t] / level) + (1 - gamma) * seasonal[index]
trend = beta * (level - previous) + (1 - beta) * phi * trend
return SmoothingState(level, trend, seasonal, period), sse
def fit(
series: Vector,
period: int,
seasonality: Seasonality = "additive",
damped: bool = True,
) -> HoltWinters:
"""Fit by bounded optimization. Raises ValueError on malformed input."""
if series.ndim != 1:
raise ValueError(f"series must be 1-D, got shape {series.shape}")
if series.size < 2 * period:
raise ValueError(f"need at least two full cycles, got {series.size} for period {period}")
if not np.isfinite(series).all():
raise ValueError("series contains NaN or inf; the recursion cannot skip a period")
if seasonality == "multiplicative" and (series <= 0).any():
raise ValueError(
"multiplicative seasonality requires strictly positive values; this series "
"reaches zero, which would divide the level by a near-zero seasonal factor"
)
def objective(packed: Vector) -> float:
alpha, beta, gamma, phi = packed
return _run(series, period, alpha, beta, gamma, phi, seasonality)[1]
bounds = [(1e-4, 0.9999)] * 3 + [(0.80, 0.9999) if damped else (1.0, 1.0)]
result = minimize(objective, np.array([0.3, 0.1, 0.1, 0.98]),
method="L-BFGS-B", bounds=bounds)
alpha, beta, gamma, phi = result.x
state, sse = _run(series, period, alpha, beta, gamma, phi, seasonality)
return HoltWinters(
state=state, alpha=float(alpha), beta=float(beta), gamma=float(gamma),
phi=float(phi), seasonality=seasonality, sse=float(sse),
)What changed: Two modelling additions and one optimizer change. Trend damping enters as a real parameter, because an undamped Holt forecast extrapolates a straight line forever and across a catalogue that guarantees occasional absurd forecasts — the damped sum converges instead. Multiplicative seasonality is supported for series whose swing scales with the level, which on sales data is the common case, and the positivity precondition it requires becomes an explicit error rather than a division by a near-zero factor. The grid search is replaced by bounded L-BFGS-B, which turns roughly seven hundred passes over the series into a few dozen.
O(n) per evaluation with tens of evaluations under L-BFGS-B instead of hundreds under a grid.
Faster, and the justification for each change.
"""Holt-Winters - vectorized ACROSS series, not across time."""
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
def fit_batch(
panel: Matrix,
period: int,
alpha: Vector,
beta: Vector,
gamma: Vector,
phi: Vector,
) -> tuple[Matrix, Vector]:
"""Run the recursion for MANY series at once.
The recursion is irreducibly sequential in time - each update reads the
state the previous step produced, and no rewrite changes that. But series
are completely INDEPENDENT of one another, so the loop over series can
disappear entirely while the loop over time remains.
That inversion is the whole optimization, and it is the shape production
demand-planning systems actually use: one Python-level iteration per
TIMESTEP for the entire catalogue, rather than one per series per step.
panel is (n_series, n_timesteps), and every parameter is a per-series
vector so a whole catalogue can be fitted with different smoothing
constants in one pass.
"""
if panel.ndim != 2:
raise ValueError(f"panel must be (series, time), got shape {panel.shape}")
if panel.shape[1] < 2 * period:
raise ValueError("every series needs at least two full cycles")
# One C-contiguous float64 block; every timestep slices a column out of it.
data = np.ascontiguousarray(panel, dtype=np.float64)
n_series, n_steps = data.shape
# State for the whole catalogue, allocated once.
level = data[:, :period].mean(axis=1)
trend = (
data[:, period : 2 * period].sum(axis=1) - data[:, :period].sum(axis=1)
) / period**2
seasonal = data[:, :period] - level[:, None] # (n_series, period)
sse = np.zeros(n_series, dtype=np.float64)
previous = np.empty(n_series, dtype=np.float64) # reused every step
prediction = np.empty(n_series, dtype=np.float64)
error = np.empty(n_series, dtype=np.float64)
for t in range(period, n_steps):
index = (t - period) % period
column = data[:, t]
season = seasonal[:, index]
# prediction = level + phi * trend + season, built in place so no
# temporary of length n_series is created per step.
np.multiply(phi, trend, out=prediction)
prediction += level
prediction += season
np.subtract(column, prediction, out=error)
sse += error * error
np.copyto(previous, level)
# level <- alpha * (y - s) + (1 - alpha) * (level + phi * trend)
np.subtract(column, season, out=prediction)
prediction *= alpha
level *= 1.0 - alpha
level += prediction
level += (1.0 - alpha) * phi * trend
# season <- gamma * (y - level) + (1 - gamma) * season
np.subtract(column, level, out=prediction)
prediction *= gamma
seasonal[:, index] *= 1.0 - gamma
seasonal[:, index] += prediction
# trend <- beta * (level - previous) + (1 - beta) * phi * trend
np.subtract(level, previous, out=prediction)
prediction *= beta
trend *= (1.0 - beta) * phi
trend += prediction
return seasonal, sseWhat changed: The recursion is sequential in time and no rewrite changes that — each update reads the state the previous step produced. What is available instead is the other axis: series are completely independent, so the loop over series disappears and only the loop over time remains. One interpreter iteration per timestep for an entire catalogue, rather than one per series per timestep, is the inversion that makes this the shape production demand-planning systems actually use. Every per-step temporary is written through preallocated buffers, because at a million series each intermediate is an eight-megabyte allocation repeated once per timestep.
Eliminate Python-level loops over samples
The loop over series is removed entirely by treating the catalogue as a matrix, so the interpreter runs once per timestep rather than once per series per timestep.
Trade-off: Every series must share the same length, seasonal period and model variant — a real constraint, since a catalogue always contains series that need multiplicative seasonality or a different period, and those have to be batched separately.
Pre-allocate output arrays and use in-place operations
The per-step prediction, error and previous-level buffers are allocated once and rewritten, avoiding several n_series-length allocations on every timestep.
Trade-off: The update chain now reuses one buffer for four different quantities in sequence, which is fast and materially harder to read — reordering any line corrupts the state silently rather than raising.
Ensure contiguous memory layout and a single dtype
A C-contiguous float64 panel means each timestep reads a column with a fixed stride and the state vectors stay cache-resident across the whole pass.
Trade-off: Series-major layout makes the per-timestep column access strided rather than sequential; the transposed layout would fix that and break the seasonal slice, so one of the two always pays.
O(n_series * n_steps) with one interpreter iteration per timestep for the whole catalogue. Illustrative, not a measured benchmark.