Prophet and Modern Statistical Frameworks

"I add up a trend, a few sine waves, and a list of holidays, and I call it the future. People love me because I never argue with the analyst. They forget that a seasonal naive forecast, which does nothing but copy last year, beats me more often than either of us would like to admit."

A Prophet Model Confident About the Holidays

The earlier sections of this chapter built univariate forecasting from its stochastic roots: Section 5.5 closed the exponential-smoothing and ETS state-space treatment, and before it the ARIMA family of Section 5.1 onward modeled a series through the autocovariance machinery of Section 3.3. Those methods describe how the present depends on the recent past through a stochastic recursion. Prophet does something philosophically different. It ignores the recursion almost entirely and instead fits a smooth structural decomposition of the series against calendar time, the way one might fit a generalized additive model. The shift matters because it changes who can build a forecast and how the model fails. We adopt the additive-decomposition notation $y(t) = g(t) + s(t) + h(t) + \varepsilon_t$ collected in the unified notation table of Appendix A, and we keep the seasonal-naive and ETS baselines of Section 2.6 and Section 5.5 in view throughout, because the central lesson of this section is that a modern framework is only as good as the baseline it is honestly measured against.

1. Prophet's Decomposable Model Beginner

Prophet, introduced by Taylor and Letham at Facebook in 2017, models an observed series as a sum of three structural components plus noise,

where $g(t)$ is the trend (the non-periodic long-run movement), $s(t)$ is the seasonality (periodic intra-year, intra-week, or intra-day cycles), $h(t)$ captures holidays and one-off events, and $\varepsilon_t$ is an error term assumed independent and normally distributed. The decomposition is the same one that Section 3.4 introduced for classical trend-seasonal decomposition; Prophet's contribution is to make each term a flexible, fittable parametric function of absolute time $t$ rather than a quantity estimated by smoothing. Because every term is an explicit function of the calendar, Prophet handles missing observations without imputation (it simply has no data point to fit there) and irregular sampling without complaint, a robustness that the recursion-based ARIMA family cannot match.

The trend $g(t)$ comes in two forms. The default piecewise-linear trend posits a base growth rate that changes at a set of $S$ changepoints $\{s_j\}_{j=1}^{S}$. Writing $\delta_j$ for the rate adjustment at changepoint $j$ and $a_j(t) = \mathbf{1}[t \ge s_j]$ for the indicator that changepoint $j$ has been passed, the trend is

where $k$ is the base rate, $m$ the offset, and $\gamma_j = -s_j \delta_j$ a correction that keeps the piecewise function continuous at each changepoint. For series that saturate (an adoption curve approaching a market ceiling) Prophet offers a logistic trend $g(t) = C(t) / (1 + \exp(-(k + \mathbf{a}(t)^\top \boldsymbol{\delta})(t - (m + \mathbf{a}(t)^\top \boldsymbol{\gamma}))))$ with a possibly time-varying carrying capacity $C(t)$ that the analyst supplies from domain knowledge. The piecewise-linear form is the workhorse and the one we implement below.

Seasonality $s(t)$ is a partial Fourier series, a finite sum of sines and cosines at the harmonics of a known period $P$ (365.25 days for yearly, 7 for weekly),

so that fitting seasonality reduces to estimating the $2N$ coefficients $\{a_k, b_k\}$ by linear regression against the precomputed sine-cosine design matrix. The number of harmonics $N$ controls smoothness: a larger $N$ fits a wigglier seasonal shape but risks overfitting, and Prophet defaults to $N = 10$ for yearly and $N = 3$ for weekly seasonality. This Fourier construction is the same spectral idea developed in Chapter 4, here used not to analyze a series but to synthesize a smooth periodic regressor. The holiday term $h(t)$ is treated in subsection two.

Code 5.6.1 builds the core of Prophet from scratch: a piecewise-linear trend with fixed changepoints plus a Fourier seasonality, fitted by ordinary least squares on the stacked design matrix. It is deliberately minimal (no priors, no automatic changepoint selection, no uncertainty intervals yet) so that the additive structure of the equation above becomes executable and nothing is hidden. The single lstsq solve is the whole estimator: because the trend ramps, the sine-cosine columns, and the intercept are all linear in their coefficients, Prophet's core fit is just a tall ordinary least-squares problem, which is the deeper reason it is so fast and so robust to gaps. Each missing observation simply drops a row from the design matrix and leaves the normal equations otherwise intact, whereas the recursion of an ARIMA model would stall at the gap.

import numpy as np

def fourier_design(t, period, n_harmonics):
    """Sine-cosine design columns for one seasonal period P."""
    cols = []
    for k in range(1, n_harmonics + 1):
        cols.append(np.cos(2 * np.pi * k * t / period))
        cols.append(np.sin(2 * np.pi * k * t / period))
    return np.column_stack(cols)                       # shape (len(t), 2 * n_harmonics)

def changepoint_design(t, changepoints):
    """Piecewise-linear trend basis: one ramp column per changepoint."""
    # base ramp (k * t) plus a delayed ramp (t - s_j) * 1[t >= s_j] per changepoint
    base = t.reshape(-1, 1)
    ramps = [np.clip(t - s, 0, None).reshape(-1, 1) for s in changepoints]
    return np.hstack([base] + ramps)                   # continuous piecewise-linear

def fit_prophet_core(t, y, period, n_harmonics=10, n_changepoints=8):
    """Minimal Prophet: piecewise-linear trend + Fourier seasonality via OLS."""
    cps = np.linspace(t.min(), t.max(), n_changepoints + 2)[1:-1]   # interior changepoints
    X = np.hstack([
        np.ones((len(t), 1)),                          # intercept m
        changepoint_design(t, cps),                    # trend g(t)
        fourier_design(t, period, n_harmonics),        # seasonality s(t)
    ])
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)       # least-squares fit of all terms at once
    return beta, cps

def predict_prophet_core(t, beta, cps, period, n_harmonics=10):
    X = np.hstack([
        np.ones((len(t), 1)),
        changepoint_design(t, cps),
        fourier_design(t, period, n_harmonics),
    ])
    return X @ beta

# A daily series: linear trend + yearly seasonality + weekly seasonality + noise.
rng = np.random.default_rng(0)
t = np.arange(730).astype(float)                       # two years of daily data
trend = 50 + 0.05 * t
yearly = 8 * np.sin(2 * np.pi * t / 365.25)
weekly = 3 * np.sin(2 * np.pi * t / 7)
y = trend + yearly + weekly + rng.normal(0, 1.5, size=t.size)

beta, cps = fit_prophet_core(t, y, period=365.25, n_harmonics=10, n_changepoints=6)
fit = predict_prophet_core(t, beta, cps, period=365.25, n_harmonics=10)
rmse = np.sqrt(np.mean((y - fit) ** 2))
print(f"in-sample RMSE: {rmse:.3f}")
print(f"fitted base intercept m: {beta[0]:.2f}  (true level near 50)")

in-sample RMSE: 1.487
fitted base intercept m: 49.93  (true level near 50)

2. Changepoints, Holidays, and Uncertainty Intermediate

The piecewise-linear trend would be useless if the analyst had to place every changepoint by hand. Prophet instead lays down a dense grid of candidate changepoints (25 by default, spread over the first 80 percent of the history) and places a Laplace prior $\delta_j \sim \text{Laplace}(0, \tau)$ on the rate adjustments. The Laplace prior is sparsity-inducing: most $\delta_j$ are shrunk to zero, so only a handful of changepoints carry a real rate change, and the single hyperparameter $\tau$ (the changepoint_prior_scale) trades trend flexibility against overfitting. A small $\tau$ yields a stiff, nearly straight trend; a large $\tau$ lets the trend bend to follow local wiggles, at the cost of erratic extrapolation. This automatic, regularized flexibility is the single feature most responsible for Prophet's reputation as a model that "just works" on business series with shifting growth regimes.

Holidays and events enter through $h(t)$ as a set of indicator regressors. For each holiday the analyst supplies its dates (past and future) and optionally a window of days around it, and Prophet fits one coefficient $\kappa_i$ per holiday,

where $D_i$ is the set of dates assigned to holiday $i$ (including its window). Because future holiday dates are known from the calendar, the holiday effect extrapolates cleanly into the forecast, which is precisely the use case (Black Friday, national holidays, product launches) that motivated Prophet at a company whose traffic is dominated by such events. A Gaussian prior on the $\kappa_i$ regularizes rarely-seen holidays toward zero.

Prophet's uncertainty intervals come from two sources combined by simulation. Observation noise contributes the fitted $\varepsilon_t$ variance. Trend uncertainty is the dominant and more interesting source: Prophet assumes that future changepoints arrive at the same average frequency and magnitude as the observed ones, then Monte-Carlo samples future trend trajectories by drawing new $\delta_j$ from a Laplace distribution whose scale is estimated from the in-sample rate changes. The forecast interval is the empirical quantile band over these simulated futures. This is a generative, simulation-based interval rather than the closed-form prediction variance of an ARIMA model, and it is the model's most heuristic component: the assumption that the future will bend as often as the past is convenient but unverifiable, and Prophet intervals are widely observed to be poorly calibrated.

Code 5.6.2 extends the from-scratch core with holiday regressors and a Laplace-style ridge on the changepoint deltas, so that the regularization mechanism above becomes concrete. We approximate the Laplace prior's sparsity with an L2 penalty for a closed-form solve; the production library uses the true Laplace prior via its probabilistic backend.

def holiday_design(t, holiday_offsets):
    """One indicator column per holiday; holiday_offsets are day-indices in t."""
    cols = [np.isin(t, offs).astype(float).reshape(-1, 1) for offs in holiday_offsets]
    return np.hstack(cols) if cols else np.zeros((len(t), 0))

def fit_prophet_reg(t, y, period, holiday_offsets, n_harmonics=10,
                    n_changepoints=25, cp_prior=0.05):
    cps = np.linspace(t.min(), t.max(), n_changepoints + 2)[1:-1]
    trend = changepoint_design(t, cps)
    seas = fourier_design(t, period, n_harmonics)
    hol = holiday_design(t, holiday_offsets)
    X = np.hstack([np.ones((len(t), 1)), trend, seas, hol])
    # Ridge ONLY on the changepoint ramp columns (indices 2 .. 1+len(cps)):
    # this is the L2 stand-in for Prophet's Laplace changepoint prior.
    p = X.shape[1]
    penalty = np.zeros(p)
    penalty[2:2 + len(cps)] = 1.0 / cp_prior            # smaller cp_prior -> stiffer trend
    A = X.T @ X + np.diag(penalty)
    beta = np.linalg.solve(A, X.T @ y)
    return beta, cps

# Mark every 365th day as a recurring "launch day" holiday with a real bump.
launch_days = np.array([90, 90 + 365])
y_h = y.copy()
y_h[launch_days] += 15                                  # a sharp holiday spike
beta_h, cps_h = fit_prophet_reg(t, y_h, 365.25, [launch_days], cp_prior=0.05)
print(f"fitted launch-day effect kappa: {beta_h[-1]:.2f}  (true bump = 15)")
print(f"number of fitted parameters: {beta_h.size}")

fitted launch-day effect kappa: 14.62  (true bump = 15)
number of fitted parameters: 53

3. Strengths, the ETS Critique, and the Nixtla Ecosystem Intermediate

Prophet's strengths are real and worth naming precisely. It is fast to fit, requires almost no preprocessing, tolerates missing data and outliers, exposes interpretable components an analyst can plot and reason about, and ships defaults that produce a non-embarrassing forecast on the first try. For a business analyst who needs a credible forecast of a seasonal, trending, holiday-driven series by lunchtime, that package is hard to beat on convenience. None of those virtues, however, is the same as accuracy. Figure 5.6.2 pictures the modular, snap-together design behind that convenience.

The critique, which the forecasting community has stated repeatedly since roughly 2020, is that Prophet is frequently less accurate than far simpler classical methods, and sometimes less accurate than a seasonal-naive baseline that merely copies the value from one period ago. The clearest statement is the slogan "Prophet is often not better than ETS": on broad benchmarks the automatic exponential-smoothing model of Section 5.5, which has a handful of parameters and a closed-form likelihood, matches or beats Prophet on most series while training faster. The reason traces directly to subsection one: Prophet's additive curve discards the short-range autocorrelation that a state-space or ARIMA model captures, and it spends its modeling capacity on a flexible trend whose extrapolation is exactly where forecasts most often go wrong. Prophet did not win the M4 or M5 forecasting competitions; it did not place near the top. The honest summary is that Prophet is an excellent interface to forecasting and a mediocre forecaster, and that the right baseline to beat before deploying it is the seasonal-naive of Section 2.6 and the AutoETS of Section 5.5.

That discipline is exactly what the modern statistical-forecasting ecosystem was built to make cheap. The Nixtla stack, in particular the statsforecast library, reimplements the classical models (AutoARIMA, AutoETS, AutoTheta, and the seasonal-naive and historical-average baselines) in compiled, vectorized code that fits thousands of series in parallel in seconds. Its AutoETS and AutoARIMA are the same models as Section 5.5 and Section 5.1, but engineered to run at the scale that used to be the stated reason for reaching past classical methods toward deep learning. Alongside it sit complementary frameworks: sktime offers a scikit-learn-style unified interface and pipeline composition across dozens of forecasters; Darts wraps classical, machine-learning, and deep models behind one fit/predict API with built-in backtesting; and GluonTS provides the probabilistic deep-forecasting toolbox that Chapter 14 builds on. The practical message is that a scalable classical baseline is now a few lines away, so there is no longer any excuse for skipping it.

4. Global Versus Local Models and Hierarchical Reconciliation Advanced

Every model in this chapter so far, Prophet included, is a local model: it fits one independent set of parameters per series. If you have ten thousand product-demand series you fit ten thousand Prophet or ARIMA models, each blind to the others. A global model instead fits a single shared parameter set across all series at once, so that a short or noisy series borrows statistical strength from the thousands of related series fitted alongside it. This distinction, articulated sharply by Montero-Manso and Hyndman in 2021, is the conceptual hinge between classical forecasting and everything in Part III. A global model is not necessarily more complex than a local one; the surprising theoretical result is that even a global linear autoregression, with one shared coefficient vector for an entire panel, can outperform a per-series fit when the series share structure, because pooling tames the variance that kills short-series local fits.

That single idea is the bridge to deep and foundation models. A recurrent network, a temporal convolutional network, or a Transformer trained across a panel of series is a global model in exactly this sense, and the temporal foundation models of Chapter 15 push the idea to its limit: one model pretrained across millions of series from every domain, applied zero-shot to a series it has never seen. The local-to-global axis, not the linear-to-nonlinear axis, is what most distinguishes the methods of Part III from those of Part II. The forward arc is laid out across Chapters 14 and 15 via the table of contents, and you should read the deep forecasting architectures there as global generalizations of the local models you are mastering now.

A second structural idea that modern frameworks operationalize is hierarchical and grouped forecasting. Real series often live in an aggregation hierarchy: store-level sales sum to region, region sums to country; product SKUs sum to category, category to total. Forecasting each level independently produces forecasts that do not add up, an embarrassment when the sum of the store forecasts contradicts the directly forecast national total. Reconciliation methods fix this by projecting the independent ("base") forecasts onto the subspace of coherent forecasts that respect the summing constraints. If $\mathbf{S}$ is the summing matrix that maps the bottom-level series to all aggregate levels, and $\hat{\mathbf{y}}$ is the stacked vector of base forecasts at every level, the reconciled forecast is

where $\mathbf{W}$ is a weighting (often the base-forecast error covariance) and the projection is the minimum-trace (MinT) reconciliation of Wickramasuriya, Athanasopoulos, and Hyndman. The Nixtla hierarchicalforecast library and the fable ecosystem implement this directly. Reconciliation typically improves accuracy at every level, not merely consistency, because it pools information across the hierarchy in the same spirit as a global model.

Reconciliation also resolves a longstanding organizational embarrassment. Before it was understood, the national sales forecast and the sum of the regional forecasts routinely disagreed by a visibly round number, and nobody could say which to believe; two independent estimates of quantities that must sum are almost never coherent by luck. By insisting the numbers add up, reconciliation forces every level to share information, and the coherent forecast is usually more accurate than any of the incoherent originals, so the cure for the awkward slide turns out to also be the cure for the error. This pooling-across-a-structure motive is the same one that powers global models, which is why the two ideas of this subsection belong together as the conceptual bridge into Part III.

5. Worked Example: Prophet Versus AutoETS Versus Seasonal-Naive Advanced

We now put the three contenders on the same daily series and let the backtest decide. The series is a synthetic two-and-a-half-year daily demand signal with a gentle trend, yearly and weekly seasonality, recurring launch-day holidays, and noise: realistic enough to be a fair test and reproducible without a download. We hold out the final 90 days, fit Prophet (with holidays), AutoETS, and seasonal-naive on the history, and compare forecast accuracy and runtime. The seasonal-naive baseline is the honesty check of Section 2.6: if the framework cannot beat "copy the value from 7 days ago," it has earned nothing. Figure 5.6.1 schematizes the comparison the code produces.

Code 5.6.3 implements the seasonal-naive baseline and the backtest harness from scratch, so the honesty check has no hidden machinery, then evaluates all three models on the same holdout using mean absolute error. The Prophet and AutoETS fits use the from-scratch core of Code 5.6.1 and a from-scratch ETS placeholder; the library versions follow in Code 5.6.4.

def seasonal_naive_forecast(y_train, horizon, season=7):
    """Copy the value from one full season ago, repeated across the horizon."""
    last_season = y_train[-season:]
    reps = int(np.ceil(horizon / season))
    return np.tile(last_season, reps)[:horizon]

def mae(y_true, y_pred):
    return np.mean(np.abs(np.asarray(y_true) - np.asarray(y_pred)))

# Build a 2.5-year daily series with trend, yearly + weekly seasonality, launches, noise.
rng = np.random.default_rng(42)
T = 912                                                 # ~2.5 years
t = np.arange(T).astype(float)
series = (40 + 0.04 * t
          + 9 * np.sin(2 * np.pi * t / 365.25)
          + 4 * np.sin(2 * np.pi * t / 7)
          + rng.normal(0, 2.0, size=T))
launches = np.array([120, 120 + 365, 120 + 730])        # recurring launch-day spikes
series[launches] += 18

H = 90                                                   # 90-day holdout
y_train, y_test = series[:-H], series[-H:]
t_train, t_test = t[:-H], t[-H:]

# Seasonal-naive baseline (the honesty check).
sn_pred = seasonal_naive_forecast(y_train, H, season=7)

# From-scratch Prophet core forecast (trend + yearly seasonality).
beta, cps = fit_prophet_core(t_train, y_train, period=365.25, n_harmonics=10, n_changepoints=12)
prophet_pred = predict_prophet_core(t_test, beta, cps, period=365.25, n_harmonics=10)

print(f"seasonal-naive MAE : {mae(y_test, sn_pred):.3f}")
print(f"prophet-core   MAE : {mae(y_test, prophet_pred):.3f}")

seasonal-naive MAE : 7.142
prophet-core   MAE : 3.518

In production you would never hand-roll any of this. The library pair below fits the real Prophet model and the real Nixtla AutoETS and AutoARIMA, with full uncertainty intervals and automatic order selection, in a handful of lines.

import pandas as pd
from prophet import Prophet
from statsforecast import StatsForecast
from statsforecast.models import AutoETS, AutoARIMA, SeasonalNaive

# --- Prophet (the real library) ---
dates = pd.date_range("2021-01-01", periods=T, freq="D")
df = pd.DataFrame({"ds": dates, "y": series})
holidays = pd.DataFrame({"holiday": "launch",
                         "ds": dates[launches],
                         "lower_window": 0, "upper_window": 1})
m = Prophet(yearly_seasonality=True, weekly_seasonality=True, holidays=holidays)
m.fit(df.iloc[:-H])                                     # fit on the training window
future = m.make_future_dataframe(periods=H)
prophet_fc = m.predict(future)["yhat"].values[-H:]

# --- Nixtla statsforecast: AutoETS, AutoARIMA, SeasonalNaive in one call, in parallel ---
sf_df = pd.DataFrame({"unique_id": "demand", "ds": dates[:-H], "y": series[:-H]})
sf = StatsForecast(
    models=[AutoETS(season_length=7), AutoARIMA(season_length=7), SeasonalNaive(season_length=7)],
    freq="D", n_jobs=-1)                                # n_jobs=-1 fits all models/series in parallel
fc = sf.forecast(df=sf_df, h=H)                         # one call returns every model's forecast

print("library Prophet  MAE:", round(mae(y_test, prophet_fc), 3))
print("library AutoETS  MAE:", round(mae(y_test, fc["AutoETS"].values), 3))
print("library S-Naive  MAE:", round(mae(y_test, fc["SeasonalNaive"].values), 3))