fpp3<\/code> ecosystem, which turns time series analysis from an arcane statistical ritual into something surprisingly readable. And we’ll be honest about where these methods fall flat, because overselling a forecasting technique is how you end up with a warehouse full of ice cream in January.<\/p>\nThe Data: Five Years of Frozen Profits<\/h2>\n
Let’s start with our running example: monthly ice cream shipments (in thousands of cases) from a mid-size manufacturer serving the U.S. Northeast region, January 2020 through December 2024.<\/p>\n
\n\n\nYear<\/th>\n Jan<\/th>\n Feb<\/th>\n Mar<\/th>\n Apr<\/th>\n May<\/th>\n Jun<\/th>\n Jul<\/th>\n Aug<\/th>\n Sep<\/th>\n Oct<\/th>\n Nov<\/th>\n Dec<\/th>\n<\/tr>\n<\/thead>\n \n\n2020<\/td>\n 42<\/td>\n 45<\/td>\n 52<\/td>\n 68<\/td>\n 89<\/td>\n 118<\/td>\n 138<\/td>\n 132<\/td>\n 105<\/td>\n 72<\/td>\n 51<\/td>\n 48<\/td>\n<\/tr>\n \n2021<\/td>\n 44<\/td>\n 47<\/td>\n 55<\/td>\n 72<\/td>\n 93<\/td>\n 124<\/td>\n 145<\/td>\n 137<\/td>\n 108<\/td>\n 75<\/td>\n 53<\/td>\n 50<\/td>\n<\/tr>\n \n2022<\/td>\n 43<\/td>\n 46<\/td>\n 54<\/td>\n 70<\/td>\n 91<\/td>\n 121<\/td>\n 142<\/td>\n 135<\/td>\n 107<\/td>\n 74<\/td>\n 52<\/td>\n 49<\/td>\n<\/tr>\n \n2023<\/td>\n 45<\/td>\n 48<\/td>\n 57<\/td>\n 74<\/td>\n 96<\/td>\n 128<\/td>\n 149<\/td>\n 141<\/td>\n 112<\/td>\n 77<\/td>\n 55<\/td>\n 52<\/td>\n<\/tr>\n \n2024<\/td>\n 46<\/td>\n 49<\/td>\n 58<\/td>\n 76<\/td>\n 99<\/td>\n 131<\/td>\n 153<\/td>\n 145<\/td>\n 115<\/td>\n 79<\/td>\n 56<\/td>\n 53<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\nEven eyeballing this table reveals the two dominant signals: July is always the peak (3.2x the January trough), and demand creeps upward about 2-3% per year. The question is whether we can quantify those signals precisely enough to forecast what comes next \u2014 and, more importantly, to understand how confident we should be in that forecast.<\/p>\n
![\"Monthly](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_demand_lineplot.png\")
<\/p>\n
Decomposition: Taking the Signal Apart<\/h2>\nSTL \u2014 The Swiss Army Knife of Decomposition<\/h3>\n
STL (Seasonal and Trend decomposition using Loess) is the most versatile decomposition method available. Developed by Cleveland et al. in 1990, it uses locally weighted regression (LOESS) to iteratively separate a time series into three additive components:<\/p>\n
y_t = Trend + Seasonal + Remainder<\/strong><\/p>\n\n- Trend<\/strong>: The underlying long-term trajectory. For our ice cream data, this is the slow upward drift from roughly 80,000 cases\/month average in 2020 to about 88,000 in 2024.<\/li>\n
- Seasonal<\/strong>: The repeating calendar-driven pattern. This is the July peak and January trough that we can set our production schedule by.<\/li>\n
- Remainder<\/strong>: Everything else \u2014 the noise, the one-off events, the unexplained variation. This is where surprises live.<\/li>\n<\/ul>\n
What makes STL superior to classical decomposition is flexibility. The seasonal component can evolve over time (useful if ice cream season is gradually starting earlier due to climate change). The trend can bend without breaking. And when you turn on the robust<\/code> option, outliers won’t derail the entire decomposition \u2014 which matters when your 2020 data includes a few pandemic-era anomalies.<\/p>\n![\"STL](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_stl_decomposition.png\")
<\/p>\n
The decomposition tells us something quantitatively useful: the strength of seasonality<\/strong> for this data is approximately 0.95 on a 0-to-1 scale (where 1 means seasonality completely dominates the signal). That’s extremely strong \u2014 but expected for ice cream, where July volumes run more than three times January’s. The strength of trend<\/strong> is around 0.4, confirming that growth exists but is modest compared to the seasonal swing.<\/p>\nThis matters for planning. A product with F_S = 0.95 demands a seasonally differentiated inventory strategy. You need production ramp-ups starting in March, peak warehouse capacity by May, and a rapid drawdown plan by September. Treating each month the same would be like staffing a ski resort identically in July and January.<\/p>\n
Reading the Seasonal Fingerprint<\/h2>\n
Decomposition gives you the broad picture, but two specialized visualizations from the feasts<\/code> package let you inspect the seasonal pattern at a finer resolution.<\/p>\nThe Seasonal Plot: Year-Over-Year Overlay<\/h3>\n
A seasonal plot (gg_season()<\/code>) overlays each year’s data on a single January-through-December axis. For our ice cream data, this produces five nearly parallel curves \u2014 one per year \u2014 that rise from January through July and fall back to December.<\/p>\n![\"Seasonal](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_seasonal_plot.png\")
<\/p>\n
What to look for:<\/p>\n
\n- Consistency<\/strong>: If the curves stack neatly, the seasonal pattern is stable. Our ice cream data is textbook-consistent.<\/li>\n
- Shifting peaks<\/strong>: If July’s peak starts migrating toward June across years, that’s a structural change worth investigating (earlier summers, shifting promotion calendars, changing consumer behavior).<\/li>\n
- Outlier years<\/strong>: A year that breaks the pattern \u2014 say 2020 dipping in April due to lockdown disruptions \u2014 stands out immediately.<\/li>\n<\/ul>\n
The Subseries Plot: Is Each Month Stable?<\/h3>\n
A subseries plot (gg_subseries()<\/code>) shows a separate mini time-series for each month, with a horizontal blue line marking that month’s mean. January gets its own panel showing all five January values over 2020-2024. February gets its own panel. And so on.<\/p>\n![\"Seasonal](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_subseries_plot.png\")
<\/p>\n
This plot answers a question the seasonal plot cannot: Is the seasonal pattern shifting over time?<\/strong> If January’s five data points show a clear upward trend within that panel, it means winter demand is growing faster than you might assume from the overall trend. For production planning, that’s the difference between maintaining a flat winter production schedule and gradually increasing it.<\/p>\nAutocorrelation: The Demand Memory Test<\/h2>\n
The Autocorrelation Function (ACF) measures how strongly today’s demand correlates with demand at various lags. For monthly data, the ACF at lag 1 asks: "Does this month’s demand tell me anything about next month’s?" At lag 12: "Does this month tell me about the same month next year?"<\/p>\n
![\"ACF](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_acf_plot.png\")
<\/p>\n
For ice cream demand, the ACF tells a clear story:<\/p>\n
\n- Strong positive spikes at lags 12, 24, and 36<\/strong>: Annual seasonality is dominant. July 2023 is highly correlated with July 2022 and July 2021. This is your most exploitable pattern.<\/li>\n
- Negative correlation around lag 6<\/strong>: Demand six months ago is anti-correlated<\/em> with current demand. Makes sense \u2014 when July was high, January was low. These are literally opposite seasons.<\/li>\n
- Slow decay at seasonal lags<\/strong>: The correlations at lag 12, 24, 36 remain strong rather than dying off quickly. This tells us the seasonal pattern is stable over the full 5-year window.<\/li>\n<\/ul>\n
The Partial Autocorrelation Function (PACF) strips out indirect effects: it shows the direct<\/em> correlation between y_t and y_{t-k} after removing the influence of intermediate lags. A sharp cutoff in the PACF at lag 12, with the spike at lag 1 also significant, suggests an ARIMA model with both non-seasonal and seasonal autoregressive terms \u2014 something like ARIMA(1,0,0)(1,1,0)_12. But we’ll let the automated model selection sort out the specifics.<\/p>\nForecasting: ETS and ARIMA<\/h2>\n
With the diagnostic work done, we can fit proper forecasting models. Two families dominate supply chain time series forecasting: ETS<\/strong> (Exponential Smoothing State Space) and ARIMA<\/strong> (AutoRegressive Integrated Moving Average). They approach the problem from different angles, and comparing them is standard practice.<\/p>\nETS: Exponential Smoothing, Properly<\/h3>\n
ETS models are named by three components: Error<\/strong> (additive or multiplicative), Trend<\/strong> (none, additive, or additive damped), and Seasonal<\/strong> (none, additive, or multiplicative). For ice cream demand with its constant seasonal amplitude and mild linear trend, automatic model selection typically lands on ETS(M,A,M)<\/strong> \u2014 multiplicative errors and seasonality with an additive trend.<\/p>\nThe model estimates three smoothing parameters:<\/p>\n
\n- Alpha<\/strong> (level): How fast the model adapts to changes in the baseline demand level. Higher alpha = more reactive, but noisier forecasts.<\/li>\n
- Beta<\/strong> (trend): How fast the model adjusts the growth rate. For a slow, steady growth market like ice cream, this tends to be small.<\/li>\n
- Gamma<\/strong> (seasonal): How fast the seasonal pattern can evolve. For a product where seasons are driven by physics (temperature) rather than fashion, a low gamma is appropriate.<\/li>\n<\/ul>\n
ARIMA: The Pattern Matching Approach<\/h3>\n
Where ETS thinks in terms of smoothed levels, ARIMA thinks in terms of differencing and correlations. For strongly seasonal monthly data, automatic selection typically produces a SARIMA<\/strong> model \u2014 something like ARIMA(1,0,1)(0,1,1)_12 \u2014 which means:<\/p>\n\n- Seasonal differencing (D=1): Subtract last year’s same-month value to remove the seasonal pattern<\/li>\n
- A seasonal moving average term (Q=1): Account for shocks that persist across seasonal cycles<\/li>\n
- Non-seasonal AR(1) and MA(1) terms: Capture month-to-month dynamics<\/li>\n<\/ul>\n
The beauty of fable::ARIMA()<\/code> is that it automates this selection process using AICc (corrected Akaike Information Criterion), testing hundreds of candidate specifications and choosing the most parsimonious model that adequately captures the data’s structure.<\/p>\nThe Benchmark: Seasonal Naive<\/h3>\n
Before celebrating any model’s performance, we need a benchmark. The Seasonal Naive<\/strong> method is the simplest possible seasonal forecast: predict that next July’s demand will equal this July’s demand. Period. No smoothing, no parameters, no optimization.<\/p>\nIf your ETS or ARIMA model can’t beat Seasonal Naive, it isn’t worth the complexity. This sounds obvious, but it’s a test that embarrassingly many production forecasting systems fail, particularly for highly seasonal products where last year’s pattern is already an excellent predictor.<\/p>\n
Putting Models Head to Head<\/h3>\n
We split the data into training (2020-2023) and test (2024) sets, fit all three models on training data, and forecast 12 months ahead.<\/p>\n
![\"Forecast](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_forecast_comparison.png\")
<\/p>\n
The accuracy comparison uses four standard metrics:<\/p>\n
\n\n\nModel<\/th>\n RMSE<\/th>\n MAE<\/th>\n MAPE<\/th>\n MASE<\/th>\n<\/tr>\n<\/thead>\n \n\nSeasonal Naive<\/td>\n ~2.4<\/td>\n ~2.2<\/td>\n ~2.3%<\/td>\n 1.00<\/td>\n<\/tr>\n \nETS<\/td>\n ~1.5<\/td>\n ~1.2<\/td>\n ~1.4%<\/td>\n ~0.55<\/td>\n<\/tr>\n \nARIMA<\/td>\n ~1.3<\/td>\n ~1.1<\/td>\n ~1.2%<\/td>\n ~0.50<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\nNote: This synthetic data is very clean \u2014 real-world demand data typically produces larger errors across all models. Run the R code below for exact values on your own data.<\/em><\/p>\n\n- RMSE<\/strong> (Root Mean Squared Error): Penalizes large misses. Lower is better.<\/li>\n
- MAE<\/strong> (Mean Absolute Error): Average miss in the same units as demand. More robust to outliers.<\/li>\n
- MAPE<\/strong> (Mean Absolute Percentage Error): Scale-independent \u2014 useful for comparing across products.<\/li>\n
- MASE<\/strong> (Mean Absolute Scaled Error): Ratio of your model’s MAE to the Naive forecast’s MAE. Below 1.0 means you’re beating the benchmark.<\/li>\n<\/ul>\n
Both ETS and ARIMA achieve MASE values well below 1.0, confirming they add value beyond the naive benchmark. The ARIMA model holds a slight edge here, but with only 12 test observations, the difference is not statistically decisive. In practice, many supply chain teams run both and average the forecasts \u2014 ensemble approaches tend to be more robust than either model alone.<\/p>\n
Prediction Intervals: What You Don’t Know Matters Most<\/h2>\n
Point forecasts are dangerous. A forecast of "148,000 cases in July" sounds precise, but supply chain decisions require understanding the range<\/em> of plausible outcomes. This is where prediction intervals earn their keep.<\/p>\nBoth ETS and ARIMA produce proper probability distributions, not just point estimates. The 80% prediction interval tells you: "There’s an 80% chance actual demand falls in this range." The 95% interval is wider, covering more extreme scenarios.<\/p>\n
For our ice cream data, a 12-month-ahead July forecast might look like:<\/p>\n
\n- Point forecast<\/strong>: 158,000 cases<\/li>\n
- 80% interval<\/strong>: 146,000 to 170,000 cases<\/li>\n
- 95% interval<\/strong>: 139,000 to 177,000 cases<\/li>\n<\/ul>\n
This is vastly more useful for planning than a single number. The lower bound of the 80% interval tells you the minimum you should produce to avoid stockouts in most scenarios. The upper bound tells you the maximum warehouse capacity you need. The gap between them \u2014 about 24,000 cases \u2014 is the quantified cost of uncertainty.<\/p>\n
Notice something important: the intervals widen as the forecast horizon increases. A 1-month-ahead forecast is much tighter than a 12-month-ahead forecast. This is not a weakness of the model \u2014 it’s an accurate reflection of reality. Uncertainty genuinely grows with time. Any forecasting system that doesn’t show this widening is hiding risk from you.<\/p>\n
Cross-Validation: Don’t Trust a Single Test<\/h2>\n
A single train\/test split can be misleading. Maybe 2024 was an unusually easy year to predict. Time series cross-validation uses a rolling origin<\/strong> approach: start with 36 months of training data, forecast 6 months ahead, add one more month of training data, forecast again, and repeat. This produces dozens of forecast-vs-actual comparisons across different time periods.<\/p>\n![\"Cross-validation](\"https:\/\/inphronesys.com\/wp-content\/uploads\/2026\/02\/ts_cv_rmse_horizon.png\")
<\/p>\n
The cross-validation results confirm the patterns we saw in the single test: both ETS and ARIMA beat Seasonal Naive at all horizons, ARIMA holds a marginal advantage, and the accuracy gap narrows as the forecast horizon lengthens. By 6 months out, even the best statistical model is only modestly better than "same as last year." This is useful information \u2014 it tells you exactly how far ahead your sophisticated models provide value.<\/p>\n
Where Time Series Analysis Breaks Down<\/h2>\n
Here’s where most blog posts would wrap up with a triumphant conclusion about the power of statistical forecasting. Instead, let’s talk about the failure modes, because knowing when your tool won’t<\/em> work is at least as important as knowing when it will.<\/p>\nThe Past-Equals-Future Assumption<\/h3>\n
Every model we’ve discussed extrapolates historical patterns. This works beautifully when the future cooperates \u2014 which it usually does for ice cream demand, because the physics of temperature and human preferences for cold treats haven’t changed much. But it fails catastrophically during structural breaks.<\/p>\n
In March 2020, no ARIMA model on earth would have predicted the simultaneous collapse of food-service ice cream demand and spike in retail take-home demand. The model doesn’t know about pandemics. It knows about lags and seasonal patterns and moving averages. When the world changes in ways it has never changed before, the model has nothing useful to say.<\/p>\n
No External Drivers<\/h3>\n
Pure ETS and ARIMA models are univariate<\/strong> \u2014 they only look at the demand series itself. They can’t incorporate temperature forecasts (which would obviously help for ice cream), promotional calendars, pricing changes, or competitor activity. The models will happily forecast "normal" July demand right through a planned July price increase that could boost volumes 20%.<\/p>\nYou can<\/em> add external regressors via dynamic regression (ARIMA(demand ~ temperature + promo)<\/code>), but this requires having those variables available at forecast time \u2014 which for temperature means weather forecasts, not historical weather. It helps, but it’s not a free lunch.<\/p>\nIntermittent Demand: The Spare Parts Problem<\/h3>\n
Time series analysis works best with continuous, relatively smooth demand patterns. Products with intermittent demand<\/strong> \u2014 many zeros punctuated by occasional, unpredictable orders \u2014 violate the basic assumptions. This is common in spare parts, industrial components, and niche products.<\/p>\nFor these items, standard ETS and ARIMA produce forecasts that are always slightly positive and never match either the zeros or the spikes. Specialized methods like Croston’s method or the Syntetos-Boylan Approximation model demand occurrence and demand size separately, which is a fundamentally different approach.<\/p>\n
The Cold-Start Problem<\/h3>\n
New products have no history. Time series models need at least two full seasonal cycles \u2014 that’s 24 months for annual seasonality \u2014 before they can reliably estimate the seasonal pattern. A new ice cream flavor launched in March has exactly zero July data points to learn from. For new product launches, you’re better served by analogous product data, market research, or judgment-based methods until sufficient history accumulates.<\/p>\n
Overfitting: When the Model Memorizes the Noise<\/h3>\n
A model with too many parameters can fit the training data exquisitely \u2014 capturing every wiggle and bump \u2014 while being terrible at forecasting new data. This is called overfitting, and it’s the statistical equivalent of studying the answer key instead of the material. An ARIMA(3,1,3)(2,1,2)_12 with 10 parameters will almost certainly fit 60 data points better than a parsimonious ARIMA(1,0,1)(0,1,1)_12 with 3 parameters. But the simpler model will usually forecast better.<\/p>\n
AICc penalizes model complexity and helps guard against overfitting. Cross-validation is even more reliable \u2014 it directly measures out-of-sample performance, which is what actually matters for supply chain planning.<\/p>\n
Correlation, Not Causation<\/h3>\n
Time series models identify pattern<\/em> \u2014 they don’t explain why<\/em>. The ACF tells you that July demand is correlated with last July’s demand, but it doesn’t tell you whether the driver is temperature, school holidays, Fourth of July marketing campaigns, or all three. If you need causal understanding to design interventions (e.g., "Would moving our summer promotion from July to June shift the seasonal peak?"), you need causal inference methods and domain expertise, not time series analysis.<\/p>\nYour Next Steps<\/h2>\n
Time series analysis won’t solve all your demand planning problems, but it will solve the most common one: extracting actionable structure from historical demand data. Here’s how to start this week:<\/p>\n
\n- \n
Decompose your top 10 SKUs.<\/strong> Run STL decomposition on your highest-volume products and compute the strength of seasonality (F_S) and strength of trend (F_T) for each. Products with F_S > 0.6 need seasonally differentiated safety stock and production plans. Products with F_T > 0.5 need trend-adjusted procurement \u2014 your buyer shouldn’t be ordering the same volume as last year for a product growing 8% annually.<\/p>\n<\/li>\n- \n
Benchmark your current forecasts against Seasonal Naive.<\/strong> Compute the MASE for whatever forecasting method you’re using today. If MASE > 1.0, your current approach is losing<\/em> to "same as last year" and you need to fix that before layering on more complexity. The R code below gives you a direct template.<\/p>\n<\/li>\n- \n
Run ETS and ARIMA on 3-5 product families.<\/strong> Use cross-validation (not a single train\/test split) and compare RMSE by forecast horizon. This tells you not just which<\/em> model is better, but how far ahead<\/em> each model provides meaningful accuracy over the naive baseline.<\/p>\n<\/li>\n- \n
Use prediction intervals for safety stock calculations.<\/strong> Stop treating forecasts as point estimates. The 95% upper bound of your prediction interval gives you a principled, statistically grounded service-level target. This replaces the guesswork of "add 20% buffer" with a calculation tied to actual forecast uncertainty.<\/p>\n<\/li>\n- \n
Document your limitations.<\/strong> For every product family, note whether any of the failure modes apply: new product (cold start), intermittent demand, heavy promotion influence, or known upcoming structural changes. Flag these for methods beyond pure time series analysis \u2014 regression, judgment, or Croston’s method as appropriate.<\/p>\n<\/li>\n<\/ol>\n\nShow R Code<\/strong><\/summary>\n# =============================================================================\n# Time Series Analysis of Ice Cream Demand\n# Complete FPP3 Pipeline\n# =============================================================================\n# Required: install.packages("fpp3")\n# This loads: tsibble, feasts, fable, fabletools, ggplot2, and more\n# =============================================================================\n\nlibrary(fpp3)\n\n# =============================================================================\n# 1. DATA PREPARATION\n# =============================================================================\n\n# Five years of monthly ice cream shipments (thousands of cases)\n# Northeast US manufacturer, 2020-2024\nice_cream <- tibble(\n month = yearmonth(seq(as.Date("2020-01-01"),\n as.Date("2024-12-01"),\n by = "month")),\n demand = c(\n 42, 45, 52, 68, 89, 118, 138, 132, 105, 72, 51, 48, # 2020\n 44, 47, 55, 72, 93, 124, 145, 137, 108, 75, 53, 50, # 2021\n 43, 46, 54, 70, 91, 121, 142, 135, 107, 74, 52, 49, # 2022\n 45, 48, 57, 74, 96, 128, 149, 141, 112, 77, 55, 52, # 2023\n 46, 49, 58, 76, 99, 131, 153, 145, 115, 79, 56, 53 # 2024\n )\n) |>\n as_tsibble(index = month)\n\n# =============================================================================\n# 2. EXPLORATORY VISUALIZATION\n# =============================================================================\n\n# Time series plot\nice_cream |>\n autoplot(demand) +\n labs(title = "Monthly Ice Cream Demand (2020-2024)",\n subtitle = "Thousands of cases \u2014 Northeast US manufacturer",\n y = "Demand (thousand cases)",\n x = NULL) +\n theme_minimal(base_size = 13)\n\n# Seasonal plot \u2014 overlay each year\nice_cream |>\n gg_season(demand) +\n labs(title = "Seasonal Pattern: Ice Cream Demand",\n subtitle = "Each line represents one year \u2014 note the consistent July peak",\n y = "Demand (thousand cases)") +\n theme_minimal(base_size = 13)\n\n# Seasonal subseries \u2014 one panel per month\nice_cream |>\n gg_subseries(demand) +\n labs(title = "Monthly Demand Subseries (2020-2024)",\n subtitle = "Blue lines show monthly means \u2014 all months show slight upward drift",\n y = "Demand (thousand cases)") +\n theme_minimal(base_size = 13)\n\n# ACF \u2014 autocorrelation diagnostics\nice_cream |>\n ACF(demand, lag_max = 36) |>\n autoplot() +\n labs(title = "Autocorrelation: Clear Annual Seasonality",\n subtitle = "Strong spikes at lags 12, 24, 36 confirm the annual cycle") +\n theme_minimal(base_size = 13)\n\n# PACF \u2014 partial autocorrelation\nice_cream |>\n PACF(demand, lag_max = 36) |>\n autoplot() +\n labs(title = "Partial Autocorrelation",\n subtitle = "Sharp cutoff at lag 12 suggests seasonal AR term") +\n theme_minimal(base_size = 13)\n\n# Combined display\nice_cream |>\n gg_tsdisplay(demand, plot_type = "season")\n\n# =============================================================================\n# 3. STL DECOMPOSITION\n# =============================================================================\n\n# STL with fixed seasonality (appropriate for stable seasonal products)\nstl_decomp <- ice_cream |>\n model(STL(demand ~ season(window = "periodic"))) |>\n components()\n\n# Plot all components\nstl_decomp |>\n autoplot() +\n labs(title = "STL Decomposition of Ice Cream Demand",\n subtitle = "Trend + Seasonal + Remainder \u2014 seasonality dominates") +\n theme_minimal(base_size = 13)\n\n# Quantify strength of trend and seasonality\nice_cream |>\n features(demand, feat_stl(s.window = "periodic"))\n# Expected output:\n# trend_strength ~ 0.4 (moderate trend)\n# seasonal_strength_year ~ 0.95 (very strong seasonality)\n# seasonal_peak_year = 7 (July)\n# seasonal_trough_year = 1 (January)\n\n# Seasonal amplitude\nstl_decomp |>\n as_tibble() |>\n summarise(\n seasonal_amplitude = max(season_year) - min(season_year),\n peak_month = month(month[which.max(season_year)]),\n trough_month = month(month[which.min(season_year)])\n )\n\n# STL with flexible seasonality (allows seasonal pattern to evolve)\nstl_flex <- ice_cream |>\n model(STL(demand ~ trend(window = 21) +\n season(window = 13),\n robust = TRUE)) |>\n components()\n\nstl_flex |>\n autoplot() +\n labs(title = "STL Decomposition (Flexible Seasonality)",\n subtitle = "season(window = 13) allows the seasonal shape to evolve over time") +\n theme_minimal(base_size = 13)\n\n# =============================================================================\n# 4. TRAIN \/ TEST SPLIT\n# =============================================================================\n\ntrain <- ice_cream |> filter(year(month) <= 2023)\ntest <- ice_cream |> filter(year(month) == 2024)\n\ncat("Training set:", nrow(train), "observations (2020-2023)\\n")\ncat("Test set: ", nrow(test), "observations (2024)\\n")\n\n# =============================================================================\n# 5. MODEL FITTING\n# =============================================================================\n\n# Fit three models on training data\nfit <- train |>\n model(\n snaive = SNAIVE(demand),\n ets = ETS(demand),\n arima = ARIMA(demand)\n )\n\n# Inspect ETS selection\nfit |> select(ets) |> report()\n# Typically selects ETS(M,A,M) or ETS(A,A,A) for this data\n\n# Inspect ARIMA selection\nfit |> select(arima) |> report()\n# Typically selects ARIMA(1,0,1)(0,1,1)[12] or similar\n\n# Residual diagnostics \u2014 check for remaining autocorrelation\nfit |> select(ets) |> gg_tsresiduals()\nfit |> select(arima) |> gg_tsresiduals()\n\n# Ljung-Box test on ARIMA residuals\naugment(fit) |>\n filter(.model == "arima") |>\n features(.innov, ljung_box, lag = 24, dof = 3)\n# p-value > 0.05 = residuals are white noise (good)\n\n# =============================================================================\n# 6. FORECASTING WITH PREDICTION INTERVALS\n# =============================================================================\n\n# Generate 12-month forecasts\nfc <- fit |> forecast(h = 12)\n\n# Plot forecasts vs actuals\nfc |>\n autoplot(ice_cream, level = c(80, 95)) +\n labs(title = "12-Month Demand Forecast Comparison",\n subtitle = "Shaded regions show 80% and 95% prediction intervals",\n y = "Demand (thousand cases)",\n x = NULL) +\n facet_wrap(~ .model, ncol = 1, scales = "free_y") +\n theme_minimal(base_size = 13)\n\n# All models on one plot\nfc |>\n autoplot(ice_cream, level = NULL) +\n labs(title = "Forecast Comparison: ETS vs ARIMA vs Seasonal Naive",\n y = "Demand (thousand cases)",\n x = NULL) +\n theme_minimal(base_size = 13)\n\n# Extract prediction intervals for a specific month\nfc |>\n hilo(level = c(80, 95)) |>\n filter(month == yearmonth("2024 Jul"))\n\n# =============================================================================\n# 7. FORECAST ACCURACY EVALUATION\n# =============================================================================\n\n# Compare against actual 2024 data\naccuracy_results <- accuracy(fc, test)\nprint(accuracy_results)\n\n# Formatted comparison table\naccuracy_results |>\n select(.model, RMSE, MAE, MAPE, MASE) |>\n arrange(MASE) |>\n mutate(across(where(is.numeric), ~ round(., 2)))\n\n# =============================================================================\n# 8. TIME SERIES CROSS-VALIDATION\n# =============================================================================\n\n# Rolling-origin cross-validation\n# Start with 36 months, add 1 month at a time, forecast 6 months ahead\ncv_data <- ice_cream |>\n stretch_tsibble(.init = 36, .step = 1)\n\ncat("Number of CV folds:", max(cv_data$.id), "\\n")\n\n# Fit models on each expanding training set\ncv_fit <- cv_data |>\n model(\n ets = ETS(demand),\n arima = ARIMA(demand),\n snaive = SNAIVE(demand)\n )\n\n# Forecast 6 months ahead from each origin\ncv_fc <- cv_fit |>\n forecast(h = 6) |>\n group_by(.id) |>\n mutate(h = row_number()) |>\n ungroup() |>\n as_fable(response = "demand", distribution = demand)\n\n# Accuracy by forecast horizon\ncv_accuracy <- cv_fc |>\n accuracy(ice_cream, by = c("h", ".model"))\n\n# Plot RMSE by horizon\ncv_accuracy |>\n ggplot(aes(x = h, y = RMSE, colour = .model)) +\n geom_line(linewidth = 1) +\n geom_point(size = 2.5) +\n scale_color_manual(values = c("arima" = "#e74c3c",\n "ets" = "#2980b9",\n "snaive" = "#95a5a6")) +\n labs(title = "Forecast Accuracy by Horizon (Cross-Validation)",\n subtitle = "Lower RMSE is better \u2014 all models degrade with longer horizons",\n x = "Forecast Horizon (months ahead)",\n y = "RMSE (thousand cases)",\n color = "Model") +\n theme_minimal(base_size = 13)\n\n# MASE by horizon\ncv_accuracy |>\n ggplot(aes(x = h, y = MASE, colour = .model)) +\n geom_line(linewidth = 1) +\n geom_point(size = 2.5) +\n geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +\n annotate("text", x = 5.5, y = 1.05, label = "Naive benchmark",\n color = "grey50", size = 3.5) +\n scale_color_manual(values = c("arima" = "#e74c3c",\n "ets" = "#2980b9",\n "snaive" = "#95a5a6")) +\n labs(title = "MASE by Horizon: Are We Beating Naive?",\n subtitle = "Values below 1.0 = better than Seasonal Naive",\n x = "Forecast Horizon (months ahead)",\n y = "MASE",\n color = "Model") +\n theme_minimal(base_size = 13)\n\n# =============================================================================\n# 9. SPECIFIC ETS MODEL VARIANTS\n# =============================================================================\n\n# Compare ETS specifications explicitly\nfit_ets_variants <- train |>\n model(\n auto = ETS(demand),\n additive = ETS(demand ~ error("A") + trend("A") + season("A")),\n multiplicative = ETS(demand ~ error("M") + trend("A") + season("M")),\n damped = ETS(demand ~ error("M") + trend("Ad") + season("M"))\n )\n\n# Compare AICc values (lower is better)\nglance(fit_ets_variants) |>\n select(.model, AICc, BIC) |>\n arrange(AICc)\n\n# Forecast comparison\nfit_ets_variants |>\n forecast(h = 12) |>\n accuracy(test) |>\n select(.model, RMSE, MAE, MASE) |>\n arrange(MASE)\n\n# =============================================================================\n# 10. APPLY TO YOUR OWN DATA\n# =============================================================================\n#\n# Replace the ice cream data with your own demand data:\n#\n# my_data <- read_csv("my_demand_data.csv") |>\n# mutate(month = yearmonth(date_column)) |>\n# as_tsibble(index = month, key = product_id)\n#\n# # Quick diagnostic\n# my_data |>\n# features(demand, feat_stl(s.window = "periodic")) |>\n# select(product_id, trend_strength, seasonal_strength_year,\n# seasonal_peak_year, seasonal_trough_year)\n#\n# # Fit and forecast\n# my_fit <- my_data |>\n# model(\n# snaive = SNAIVE(demand),\n# ets = ETS(demand),\n# arima = ARIMA(demand)\n# )\n#\n# my_fc <- my_fit |> forecast(h = 12)\n#\n# # Evaluate\n# my_fc |> accuracy(my_test_data) |>\n# select(.model, RMSE, MAE, MASE) |>\n# arrange(MASE)\n#\n# # If MASE > 1.0 for ETS\/ARIMA, your data may have characteristics\n# # that these models can't handle well (intermittent demand, structural\n# # breaks, strong external drivers). Consider:\n# # - Croston's method for intermittent demand\n# # - Dynamic regression for external drivers\n# # - Judgment-based adjustments for known structural changes\n<\/code><\/pre>\n<\/details>\nInteractive Dashboard<\/h2>\n
Explore the data yourself \u2014 adjust the ETS smoothing parameters, switch between model types, and see how decomposition, seasonal patterns, and forecast accuracy change in real time.<\/p>\n
\nInteractive Dashboard<\/strong><\/p>\nExplore the data yourself \u2014 adjust parameters and see the results update in real time.<\/p>\n
Open Interactive Dashboard →<\/a><\/div>\nReferences<\/h2>\n\n- Hyndman, R.J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice<\/em>, 3rd edition. OTexts. https:\/\/otexts.com\/fpp3\/<\/a><\/li>\n
- Cleveland, R.B., Cleveland, W.S., McRae, J.E., & Terpenning, I. (1990). "STL: A Seasonal-Trend Decomposition Procedure Based on Loess." Journal of Official Statistics<\/em>, 6(1), 3-73.<\/li>\n
- Wang, X., Smith, K.A., & Hyndman, R.J. (2006). "Characteristic-based clustering for time series data." Data Mining and Knowledge Discovery<\/em>, 13(3), 335-364.<\/li>\n
- FRED Blog (2024). "Ice cream is a seasonal product, right?" Federal Reserve Bank of St. Louis. https:\/\/fredblog.stlouisfed.org\/2024\/05\/ice-cream-is-a-seasonal-product-right\/<\/a><\/li>\n<\/ol>\n","protected":false},"excerpt":{"rendered":"
Your demand data is trying to tell you something. We use STL decomposition, seasonal diagnostics, and ETS\/ARIMA models to extract trend, seasonality, and noise from ice cream sales data \u2014 then honestly discuss where these methods break down.<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[13,115],"tags":[126,122,125,128,127,119,124,123,26,121],"class_list":["post-1096","post","type-post","status-publish","format-standard","hentry","category-data-science","category-supply-chain-management","tag-arima","tag-demand-forecasting","tag-ets","tag-forecast-accuracy","tag-fpp3","tag-r-programming","tag-seasonality","tag-stl-decomposition","tag-supply-chain-analytics","tag-time-series-analysis"],"_links":{"self":[{"href":"https:\/\/inphronesys.com\/index.php?rest_route=\/wp\/v2\/posts\/1096","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/inphronesys.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/inphronesys.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/inphronesys.com\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/inphronesys.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=1096"}],"version-history":[{"count":2,"href":"https:\/\/inphronesys.com\/index.php?rest_route=\/wp\/v2\/posts\/1096\/revisions"}],"predecessor-version":[{"id":1098,"href":"https:\/\/inphronesys.com\/index.php?rest_route=\/wp\/v2\/posts\/1096\/revisions\/1098"}],"wp:attachment":[{"href":"https:\/\/inphronesys.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1096"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/inphronesys.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1096"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/inphronesys.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1096"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}