Midterm Forecasting Project
Comparing a seasonal naive benchmark, ETS, ARIMA, time-series regression and an equal-weight forecast combination
Rudy Lazaro | ADAN/ADEC 7406 Predictive Analytics and Forecasting | Woods College of Advancing Studies, Boston College
Rendered 23 October 2026
A note on the numbers. This document pulls the series live from FRED each time it is rendered. FRED revises and adds months, so a later run will shift every number below. All figures quoted in the poster come from the render whose output appears here.
1. Series selection
The series measures the retail sales of building material, garden equipment, and supplies dealers. The audience would be the companies that sell these goods. The forecast would help a company decide how much of these goods to stock given the demand for them. We chose a monthly series for more observations, and a multi-year horizon for the same reason.
Getting this forecast wrong could be disastrous for these companies. If we overshoot it, the companies will lose money because they will have excess goods. However, if we undershoot it, they will lose money because they could have sold more. An accurate forecast here will be optimal for the company to make the most profit possible.
| Series | MRTSSM444USN |
| Definition | Retail Sales: Building Material, Garden Equipment and Supplies Dealers |
| Units | Millions of dollars, not seasonally adjusted |
| Frequency | Monthly |
| Source | U.S. Census Bureau, Monthly Retail Trade Survey, via FRED |
| Direct link | https://fred.stlouisfed.org/series/MRTSSM444USN |
| Retrieval | Keyless public CSV endpoint, no API key in this script |
We chose to start in 2005, a few years before the Great Recession, so that we capture both that recession and COVID in our data set. The window therefore spans the 2008-09 housing downturn and the 2020-21 pandemic, which is what makes this a demanding forecasting problem rather than a routine one.
library(fpp3)
library(ggplot2)
set.seed(7406)
bm <- read.csv("https://fred.stlouisfed.org/graph/fredgraph.csv?id=MRTSSM444USN")
names(bm) <- c("date", "value")
bm <- bm |>
dplyr::mutate(date = as.Date(date),
value = suppressWarnings(as.numeric(value))) |>
dplyr::filter(!is.na(value), date >= as.Date("2005-01-01")) |>
dplyr::mutate(month = yearmonth(date)) |>
as_tsibble(index = month) |>
dplyr::select(month, value)
stopifnot(nrow(scan_gaps(bm)) == 0) # no missing months2. Protecting the test period
Given the nature of the data, we decided a chronological split would be best. The forecast origin was April 2022, and no test data touched estimation or model selection.
Every transformation, differencing decision, feature calculation and
model fit below is computed from train only.
test appears for the first time in section 10, where it is
used once, to score forecasts.
n <- nrow(bm)
k <- floor(0.8 * n)
train <- bm |> dplyr::slice(1:k)
test <- bm |> dplyr::slice((k + 1):n)
h <- n - k| Item | Value |
|---|---|
| Series | MRTSSM444USN |
| Units | Millions of dollars, not seasonally adjusted |
| Frequency | Monthly |
| Full sample | 2005 Jan - 2026 Aug |
| Observations | 260 |
| Training set | 2005 Jan - 2022 Apr |
| Training n | 208 |
| Test set | 2022 May - 2026 Aug |
| Test n (h) | 52 |
3. Diagnosing the training series
Figure 1. The full series
The data stays somewhat steady at first. After the 2008 recession we see a dip, and from there it shows an upward trend. We then see a second surge right before our test period. This may be because of a recovery period after COVID.
bm |> ggplot(aes(month, value)) +
geom_line(colour = "grey30") +
geom_vline(xintercept = as.Date(train$month[k]),
linetype = "dashed", colour = "firebrick") +
annotate("text", x = as.Date(train$month[k]), y = max(bm$value),
label = " test period >", hjust = 0, colour = "firebrick", size = 3.2) +
labs(title = "Figure 1. Retail sales, building material and garden supply dealers",
subtitle = paste0("Monthly, NSA. Dashed line = end of training set (",
format(train$month[k]), ")"),
x = NULL, y = "Millions of dollars") +
theme_minimal()
Figure 1. Retail sales, with the end of the training set marked.
Figure 2. Seasonal plot
February is the lowest month and May is the highest within a given year.
train |> gg_season(value, labels = "none") +
labs(title = "Figure 2. Seasonal plot by year",
x = NULL, y = "Millions of dollars") +
theme_minimal()
Figure 2. Seasonal plot by year.
Figure 3. Subseries plot
Even within a given month, we see a similar trend across all of the months. This shows that the seasonality of a given month is persistent regardless of the series’ overall trend.
train |> gg_subseries(value) +
labs(title = "Figure 3. Subseries plot: each month over time",
subtitle = "Blue line = mean for that month",
x = NULL, y = "Millions of dollars") +
theme_minimal()
Figure 3. Subseries plot, each month over time.
Figure 4. STL decomposition
We see in the STL decomposition more confirmation that the series is seasonal. In the remainder we see a surge between 2020 and 2021.
train |> model(STL(value ~ trend(window = 21) + season(window = "periodic"),
robust = TRUE)) |>
components() |> autoplot() +
labs(title = "Figure 4. STL decomposition", x = NULL) +
theme_minimal()
Figure 4. STL decomposition.
Figure 5. Levels, ACF and PACF
The ACF shows a staircase-like pattern, and the staircase gets smaller as the lag increases. This tells us that the farther out a value is from the current one, the less of an effect it has on the current value.
train |> gg_tsdisplay(value, plot_type = "partial", lag_max = 36) +
labs(title = "Figure 5. Training series with ACF and PACF (levels)")
Figure 5. Training series with ACF and PACF, in levels.
Figure 6. After first and seasonal differencing
We see in the differenced series the innovation residual did not change and the ACF and PACF don’t have the same lag structure as the normal series. However, we do see some of the seasonality in them both. Figure 6 shows that even numbered lags explain more about the current value than other lags.
train |> gg_tsdisplay(difference(difference(value, 12)),
plot_type = "partial", lag_max = 36) +
labs(title = "Figure 6. After seasonal and first differencing",
y = "Differenced value")
Figure 6. Series, ACF and PACF after one first and one seasonal difference.
Features, unit-root tests and the variance-level relationship
We decided to use one seasonal difference plus one first difference because Table T2 gives nsdiffs = 1 and ndiffs = 1.
We decided to use log(value) for the regression because we see that the annual level and within-year SD correlation is 0.887 and the Guerrero lambda is 0.58.
t2 <- dplyr::bind_cols(
train |> features(value, feat_stl) |>
dplyr::select(trend_strength, seasonal_strength_year),
train |> features(value, list(unitroot_nsdiffs, unitroot_ndiffs, unitroot_kpss))
)
vl <- train |> as_tibble() |>
dplyr::mutate(yr = lubridate::year(month)) |>
dplyr::group_by(yr) |>
dplyr::summarise(level = mean(value), sd = sd(value), .groups = "drop")
lambda <- train |> features(value, guerrero) |> dplyr::pull(lambda_guerrero)| Statistic | Value |
|---|---|
| Trend strength (STL) | 0.9719 |
| Seasonal strength, yearly (STL) | 0.9416 |
| nsdiffs (seasonal differences) | 1 |
| ndiffs (first differences) | 1 |
| KPSS statistic | 2.2588 |
| KPSS p-value | 0.01 (table floor) |
| corr(annual level, within-year SD) | 0.887 |
| Guerrero Box-Cox lambda | 0.5819 |
4. Benchmark
For the benchmark, we decided a seasonal naive would work best given that the series has strong seasonal strength.
A forecasting model is only worth using if it outperforms a benchmark, which is why every model here is measured against the seasonal naive (Hyndman & Athanasopoulos, 2021).
5. ETS
For the ETS model we chose ETS(M,N,A). This means a multiplicative error, no trend component, and an additive seasonal component. Figure 1 and the trend strength of 0.9719 in Table T2 both show a clear upward trend, so the N does not mean the series has no trend. It means the fitted model carries no separate trend term, and the level component absorbs that growth. Because there is no trend term there is no damping parameter to report.
Having no trend implies that the overall length should not affect the forecast as much. This means that the overall level and the season will have a stronger effect.
6. ARIMA
For the ARIMA we selected the form ARIMA(2,1,2)(0,1,2)[12]. Here p = 2 is the number of non-seasonal autoregressive terms, d = 1 the order of non-seasonal differencing, q = 2 the order of the non-seasonal moving average, P = 0 indicates no seasonal autoregressive term, D = 1 the degree of seasonal differencing, Q = 2 the order of the seasonal moving average, and 12 is the seasonal period, since the data is monthly.
The d = 1 and D = 1 the model chose match the ndiffs = 1 and nsdiffs = 1 reported in Table T2, so the differencing is consistent with the observed temporal structure.
7. Time-series regression
For the regression, we used the log of the value as our dependent variable; for our independent variables we used a trend variable and seasonal dummies. The assumption behind a log-linear model is that the data will keep increasing nominally, so it makes more sense to look at percent changes, which the log does.
Both predictors are justified by the structure of the series rather than by fit alone: the trend term is there because Table T2 reports a trend strength of 0.9719, and the monthly dummies are there because the seasonal strength is 0.9416 and Figures 2 and 3 show that pattern repeating every year.
On future predictor availability. Both predictors, the trend index and the month, are known for any future date. Neither has to be forecast before the regression can be used, so this model does not create a second forecasting problem for itself.
8. Fitting the four models
All four are fitted to train only.
fit <- train |>
model(
bench = SNAIVE(value), # seasonal naive benchmark
ets = ETS(value), # automatic ETS
arima = ARIMA(value), # automatic ARIMA
reg = TSLM(log(value) ~ trend() + season()) # log-linear trend + month dummies
)| Model | Specification |
|---|---|
| Benchmark | SNAIVE(value) |
| ETS | ETS(M,N,A) |
| ARIMA | ARIMA(2,1,2)(0,1,2)[12] |
| Regression | TSLM(log(value) ~ trend() + season()) |
| r_squared | adj_r_squared | sigma2 | statistic | p_value |
|---|---|---|---|---|
| 0.66896 | 0.64859 | 0.01817 | 32.83781 | 0 |
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 9.75857 | 0.03559 | 274.16131 | 0.00000 |
| trend() | 0.00215 | 0.00016 | 13.83044 | 0.00000 |
| season()year2 | -0.02228 | 0.04494 | -0.49592 | 0.62051 |
| season()year3 | 0.23809 | 0.04494 | 5.29849 | 0.00000 |
| season()year4 | 0.35487 | 0.04494 | 7.89693 | 0.00000 |
| season()year5 | 0.43752 | 0.04559 | 9.59631 | 0.00000 |
| season()year6 | 0.36806 | 0.04559 | 8.07312 | 0.00000 |
| season()year7 | 0.29257 | 0.04559 | 6.41735 | 0.00000 |
| season()year8 | 0.25383 | 0.04559 | 5.56741 | 0.00000 |
| season()year9 | 0.20622 | 0.04559 | 4.52305 | 0.00001 |
| season()year10 | 0.24469 | 0.04559 | 5.36673 | 0.00000 |
| season()year11 | 0.18342 | 0.04560 | 4.02270 | 0.00008 |
| season()year12 | 0.12032 | 0.04560 | 2.63876 | 0.00899 |
The adjusted R-squared was 0.649, which means the model explained about 65% of the variance. That figure describes the training fit only; it is not evidence about forecast accuracy, which is settled in section 10.
| .model | AICc |
|---|---|
| bench | NA |
| ets | 4056.578 |
| arima | 3341.109 |
| reg | -816.880 |
The regression AICc is not on the same scale as the ETS or ARIMA numbers, so it can’t be compared to them. However, the ARIMA model did score lower than the ETS model.
9. Residual diagnostics
fit |> dplyr::select(ets) |> gg_tsresiduals(lag_max = 36) +
labs(title = "Figure 7. ETS residual diagnostics")
Figure 7. ETS residual diagnostics.
The innovation residuals are mostly stable, with the exception of a sudden surge around 2022. The ACF produces an unusual pattern, but it appears to show some seasonality. The histogram looks good, since most of the values are very close to zero.
fit |> dplyr::select(arima) |> gg_tsresiduals(lag_max = 36) +
labs(title = "Figure 8. ARIMA residual diagnostics")
Figure 8. ARIMA residual diagnostics.
For the ARIMA, we see the innovation residuals also follow a stable pattern with a surge around 2022. We also see a similarly unusual pattern in the ACF. The histogram does not look good at all; a lot of the values are nowhere near zero.
fit |> dplyr::select(reg) |> gg_tsresiduals(lag_max = 36) +
labs(title = "Figure 9. Regression residual diagnostics")
Figure 9. Regression residual diagnostics.
The regression residuals are not good at all. The innovation residuals show a clear pattern, and the ACF also shows a downward trend. The histogram is the only saving grace, with all of the values being near zero.
t5 <- augment(fit) |>
features(.innov, ljung_box, lag = 24) |>
dplyr::select(.model, lb_stat, lb_pvalue) |>
dplyr::arrange(lb_stat)| .model | lb_stat | lb_pvalue |
|---|---|---|
| arima | 46.19262 | 0.0042 |
| ets | 170.12915 | 0.0000 |
| bench | 927.96134 | 0.0000 |
| reg | 2266.45984 | 0.0000 |
The ARIMA model comes closest to white noise. There isn’t as much seasonality in the ACF compared to the other two models, especially in the lags closer to the current value. All four models reject the null hypothesis of the Ljung-Box test, which means that all of the models carry some sort of temporal structure.
10. Forecast combination and the test period
The combination is the equal-weight mean of the ETS, ARIMA and regression point forecasts, which is the classical construction of Bates & Granger (1969). The benchmark is deliberately excluded from it: the combination is being tested as a competitor to the benchmark, not built from it.
All five approaches forecast the same 52 test months from the same forecast origin.
fc <- fit |> forecast(h = h)
combo <- fc |> as_tibble() |>
dplyr::filter(.model %in% c("ets", "arima", "reg")) |>
dplyr::group_by(month) |>
dplyr::summarise(.mean = mean(.mean), .groups = "drop") |>
dplyr::mutate(.model = "combo")
actual <- test |> as_tibble() |> dplyr::select(month, actual = value)
all_fc <- dplyr::bind_rows(
fc |> as_tibble() |> dplyr::select(.model, month, .mean),
combo) |>
dplyr::left_join(actual, by = "month")ggplot() +
geom_line(data = all_fc, aes(month, .mean, colour = .model), linewidth = 0.6) +
geom_line(data = bm |> dplyr::filter(month >= yearmonth("2017 Jan")),
aes(month, value), colour = "black", linewidth = 0.8) +
geom_vline(xintercept = as.Date(train$month[k]),
linetype = "dashed", colour = "grey50") +
scale_colour_brewer(palette = "Set1", name = "Model") +
labs(title = "Figure 10. Test-period forecasts against actuals",
subtitle = paste0("Black = observed. Forecast origin ", format(train$month[k]),
", horizon h = ", h, " months"),
x = NULL, y = "Millions of dollars") +
theme_minimal()
Figure 10. Test-period forecasts against actuals. This is the primary forecast-versus-actual figure.
All of the models follow the same trends as the actual data, though some track the actuals more closely than others.
Test-set accuracy
# MASE scale: in-sample mean absolute seasonal-naive error on the training set
mase_scale <- mean(abs(diff(train$value, lag = 12)))
acc <- all_fc |>
dplyr::group_by(.model) |>
dplyr::summarise(MAE = mean(abs(actual - .mean)),
RMSE = sqrt(mean((actual - .mean)^2)),
MASE = mean(abs(actual - .mean)) / mase_scale,
MAPE = mean(abs((actual - .mean) / actual)) * 100,
Bias = mean(actual - .mean),
.groups = "drop") |>
dplyr::arrange(MASE)| .model | MAE | RMSE | MASE | MAPE | Bias |
|---|---|---|---|---|---|
| ets | 1541.446 | 1903.114 | 0.82461 | 3.78570 | -742.6050 |
| combo | 1741.894 | 2188.171 | 0.93184 | 3.99790 | 765.3957 |
| bench | 2264.058 | 2769.196 | 1.21117 | 5.06649 | 1306.9423 |
| arima | 3945.799 | 4555.370 | 2.11083 | 9.47397 | -3465.1980 |
| reg | 6503.990 | 6809.534 | 3.47935 | 15.09552 | 6503.9902 |
MAE, RMSE and Bias are in millions of dollars. The MASE denominator is the in-sample mean absolute seasonal-naive error on the training set, 1869.31. Bias is mean(actual minus forecast), so positive means under-forecast and negative means over-forecast.
| .model | Role | MAE | RMSE | MASE |
|---|---|---|---|---|
| ets | Component | 1541.446 | 1903.114 | 0.82461 |
| combo | Combination | 1741.894 | 2188.171 | 0.93184 |
| arima | Component | 3945.799 | 4555.370 | 2.11083 |
| reg | Component | 6503.990 | 6809.534 | 3.47935 |
11. Interpreting the results
Which model performs best. Only the ETS and combination models have lower MAE, RMSE, MASE, and MAPE values than the benchmark, while the ARIMA and regression models do not. This means that those two models did not earn their place.
Does the ranking change across metrics. The rankings do not change across MAE, RMSE, MASE and MAPE. The ranking is the same across all metrics, which means the ETS is the model that clearly outperforms the others.
How much the best model improves on the benchmark. Comparing the benchmark to the ETS we see an MAE of 1541.4 against 2264.1 and a MASE of 0.825 against 1.211, a 32% cut, which means the ETS model outperformed the benchmark by almost a third.
Systematic bias. The ETS and ARIMA models over-forecast, while the regression, seasonal naive, and combination models under-forecast. The log-linear assumption was strained because the data grew too much since the other values before the pandemic carried a lot more of the weight for the regression. Since the other models carry more weight to the more recent values, it makes sense why they did a better job at forecasting.
Did the combination help. The combination model gave every component equal weight, but weight must be earned rather than given, which is why we see the ETS model outperform it. It beat two of its three components and lost to the third.
What characteristics of the series explain the winner. The ETS outperformed the other models because it was able to accurately carry both the multiplicative error, shown by the 0.887 correlation between annual level and within-year SD, and the seasonality, shown by the seasonal strength of 0.9416, without needing to carry a trend alongside them. In the ARIMA and the regression the trend was carried, and both deviated from the actuals by a larger margin. It also beat the seasonal naive since the seasonal naive only carried the seasonality of the data without the multiplicative error.
Did residual diagnostics correspond with out-of-sample performance. No. The ARIMA fit best in sample but finished fourth out of sample, because while it may be best when only looking at one step ahead, it may not be best when looking further ahead. This is the distinction FPP3 draws between residuals, which are based on one-step forecasts, and forecast errors, which can involve multi-step forecasts. AICc and the Ljung-Box statistic are both built on one-step residuals; the test period here is 52 steps.
12. Final recommendation
The model I would use is the ETS model, because it forecasts best compared to the other models. The ARIMA and regression did not outperform the benchmark and the combination still did worse than the ETS model. That said, the models used in this exercise may not be optimal, given the limitations of the project.
I would still ship the ETS because the model is simple, the model is easy to interpret, and it needs no future predictor values at all.
Limitations and what could cause the forecast to fail
The 52-month horizon is not a separate choice: 52 is 20% of our 260 observations, so it follows directly from the chronological 80/20 split. A horizon that long makes the later months much less reliable, because the models receive no new information after the April 2022 forecast origin.
We keep this horizon to keep the 80/20 split. This is because using more test data will give the models bias since the models may not be able to predict other data sets than this one. By not shortening it we make the model more usable for other data sets.
We also see from the Ljung-Box results that all of the models still carry some temporal structure, which would make their forecasts less trustworthy.
The pandemic was an abnormal event which means that the data is working with something abnormal which would affect the accuracy of its forecasts.
I think I would try to use more of the historical data to keep the 80/20 split while still having more data.
References
Bates, J.M., & Granger, C.W.J. (1969). The combination of forecasts. Operational Research Quarterly, 20(4), 451-468.
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, 6(1), 3-73.
Hyndman, R.J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. https://otexts.com/fpp3/
Hyndman, R.J., & Koehler, A.B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4), 679-688.
Ljung, G.M., & Box, G.E.P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297-303.
Data citation. U.S. Census Bureau, Monthly Retail Trade Survey. Retail Sales: Building Material, Garden Equipment and Supplies Dealers [MRTSSM444USN], retrieved from FRED, Federal Reserve Bank of St. Louis, https://fred.stlouisfed.org/series/MRTSSM444USN.
Reproducibility
The analysis runs from a clean session end to end with no manual steps. The data is pulled from the keyless FRED CSV endpoint, so no API key appears anywhere in this file. The random seed is set at the top; the chronological split, forecast origin and horizon are all derived from the data rather than hard-coded.
## R version 4.4.3 (2025-02-28 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] fable_0.5.0 feasts_0.5.0 fabletools_0.8.0 ggtime_1.0.0 tsibbledata_0.4.1
## [6] tsibble_1.2.0 ggplot2_4.0.3 lubridate_1.9.5 tidyr_1.3.1 dplyr_1.1.4
## [11] tibble_3.2.1 fpp3_1.0.3
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.4 sass_0.4.10 generics_0.1.3 anytime_0.3.13
## [5] lattice_0.22-6 digest_0.6.37 magrittr_2.0.3 evaluate_1.0.5
## [9] grid_4.4.3 timechange_0.4.0 RColorBrewer_1.1-3 mixtime_0.3.0
## [13] fastmap_1.2.0 jsonlite_2.0.0 purrr_1.0.4 scales_1.4.0
## [17] numDeriv_2016.8-1.1 jquerylib_0.1.4 cli_3.6.4 rlang_1.1.5
## [21] crayon_1.5.3 vecvec_1.3.0 withr_3.0.2 cachem_1.1.0
## [25] yaml_2.3.12 otel_0.2.0 tools_4.4.3 tzdb_0.5.0
## [29] vctrs_0.6.5 R6_2.6.1 lifecycle_1.0.4 urca_1.3-4
## [33] pkgconfig_2.0.3 progressr_1.0.0 pillar_1.10.1 bslib_0.12.0
## [37] gtable_0.3.6 glue_1.8.0 Rcpp_1.1.1 xfun_0.57
## [41] tidyselect_1.2.1 rstudioapi_0.19.0 knitr_1.52 farver_2.1.2
## [45] nlme_3.1-167 htmltools_0.5.9 labeling_0.4.3 rmarkdown_2.32
## [49] compiler_4.4.3 S7_0.2.2 distributional_0.8.1