Forecasting U.S. Retail Alcohol Sales Using Exponential Smoothing
Luke Volm
2025-11-20
1 Introduction
Forecasting plays a central role in understanding how economic indicators evolve over time, especially for industries influenced by both long-term trends and seasonal consumer behavior. In this analysis, we examine more than thirty years of monthly U.S. retail sales for beer, wine, and liquor stores using a range of exponential smoothing models. Because the dataset displays clear upward growth as well as pronounced seasonal patterns, it provides an ideal setting for comparing simple, trend-based, and fully seasonal smoothing approaches. Our goal is to evaluate how well each model captures the underlying structure of the series and to determine which forecasting method delivers the most accurate predictions for the 12-month holdout period. By systematically fitting and assessing all relevant smoothing models, this study identifies the most effective approach for producing reliable short-term forecasts of retail alcohol sales.
2 Description of Data Set
The dataset used in this analysis comes from the Federal Reserve Bank of St. Louis FRED database and contains monthly U.S. retail sales for beer, wine, and liquor stores (NAICS 4453). Each observation reports the total dollar value of nationwide sales for these retail establishments and is measured in millions of dollars. The series begins in January 1992 and extends through 2025, providing over thirty years of monthly data with clear long-term trends and strong seasonal patterns related to consumer spending. Because the dataset is comprehensive, up-to-date, and reflects national retail activity across the alcoholic beverage sector, it is well-suited for comparing exponential smoothing models and evaluating forecasting performance.
# Create the full time series object
sales.ts <- ts(sales$MRTSSM4453USN,
start = c(1992, 1),
frequency = 12)
# ---- HOLD OUT LAST 12 MONTHS ----
n <- length(sales.ts)
h <- 12 # monthly data → hold out 12
# Training data (all but last 12 points)
train.sales <- sales.ts[1:(n - h)]
# Testing data (the last 12 months)
test.sales <- sales.ts[(n - h + 1):n]
# Create TS object for the training portion for modeling
sales.train.ts <- ts(train.sales,
start = c(1992, 1),
frequency = 12)First, the retail sales data were converted into a monthly time series object starting in January 1992 with a frequency of 12 to reflect the monthly observation pattern. To evaluate forecasting performance, the most recent 12 months of the series were set aside as a test set, while all earlier observations were used as training data. The training portion was then stored as its own time series object and served as the basis for fitting the exponential smoothing models, while the held-out 12 months were kept for comparing forecast accuracy across models.
3 Smoothing Models
Now, we will create and plot three different types of smoothing models: simple exponential smoothing, Holt models, and Holt–Winters models. We will report accuracy measures based on the training data and forecast errors. For all modeling, we will be using the 12 most recent observations that we created above.
3.1 Simple Exponential Model
The first model we’ll create is the simple exponential model (SES).
Forecast method: Simple exponential smoothing
Model Information:
Simple exponential smoothing
Call:
ses(y = sales.train.ts, h = 12)
Smoothing parameters:
alpha = 0.1294
Initial states:
l = 1614.0721
sigma: 469.8404
AIC AICc BIC
7148.931 7148.993 7160.837
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 85.74636 468.6373 301.8894 1.033189 8.872947 2.050725 -0.02974129
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 5951.93 5349.806 6554.055 5031.060 6872.801
Sep 2024 5951.93 5344.787 6559.074 5023.384 6880.476
Oct 2024 5951.93 5339.809 6564.052 5015.771 6888.089
Nov 2024 5951.93 5334.871 6568.990 5008.220 6895.641
Dec 2024 5951.93 5329.973 6573.888 5000.728 6903.133
Jan 2025 5951.93 5325.112 6578.748 4993.295 6910.566
Feb 2025 5951.93 5320.290 6583.571 4985.919 6917.941
Mar 2025 5951.93 5315.503 6588.357 4978.599 6925.261
Apr 2025 5951.93 5310.753 6593.108 4971.334 6932.527
May 2025 5951.93 5306.037 6597.823 4964.122 6939.739
Jun 2025 5951.93 5301.356 6602.505 4956.962 6946.898
Jul 2025 5951.93 5296.708 6607.153 4949.854 6954.007The simple exponential model estimates a baseline level of ℓ = 1614.07, which represents the constant smoothed value the model uses since SES does not capture trend or seasonality. The error measures show that SES does not fit the data very well: the model’s RMSE (468.64) and MAPE (8.87%) indicate relatively large forecast errors for monthly retail sales, which is expected because the series has clear trend and seasonal patterns that SES cannot model. The information criteria (AIC = 7148.93) are higher than what we typically see with more flexible models, reinforcing that SES is too simple for this dataset. Overall, the SES model provides a weak baseline and serves mainly as a comparison point for more appropriate trend- and seasonality-based smoothing models.
3.2 Holt Model
Now, we will create and plot our Holt Model.
# Holt linear trend (additive trend, optimal alpha and beta)
fit2 <- holt(sales.train.ts,
initial = "optimal",
h = 12)
summary(fit2)
Forecast method: Holt's method
Model Information:
Holt's method
Call:
holt(y = sales.train.ts, h = 12, initial = "optimal")
Smoothing parameters:
alpha = 0.0608
beta = 0.0013
Initial states:
l = 1558.9422
b = 5.0016
sigma: 456.8164
AIC AICc BIC
7128.932 7129.088 7148.776
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 28.99928 454.4737 289.174 -0.6699864 8.571987 1.964349 0.007190073
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 6182.523 5597.089 6767.957 5287.180 7077.867
Sep 2024 6202.399 5615.838 6788.960 5305.332 7099.467
Oct 2024 6222.275 5634.542 6810.009 5323.414 7121.137
Nov 2024 6242.151 5653.198 6831.105 5341.425 7142.878
Dec 2024 6262.027 5671.807 6852.247 5359.364 7164.691
Jan 2025 6281.903 5690.369 6873.438 5377.229 7186.578
Feb 2025 6301.780 5708.881 6894.678 5395.020 7208.539
Mar 2025 6321.656 5727.345 6915.966 5412.736 7230.575
Apr 2025 6341.532 5745.758 6937.305 5430.375 7252.688
May 2025 6361.408 5764.121 6958.694 5447.937 7274.879
Jun 2025 6381.284 5782.433 6980.135 5465.420 7297.147
Jul 2025 6401.160 5800.693 7001.627 5482.825 7319.495 ME RMSE MAE MPE MAPE MASE ACF1
Training set 28.99928 454.4737 289.174 -0.6699864 8.571987 1.964349 0.007190073
Forecast method: Damped Holt's method
Model Information:
Damped Holt's method
Call:
holt(y = sales.train.ts, h = 12, damped = TRUE)
Smoothing parameters:
alpha = 0.0201
beta = 0.0201
phi = 0.9139
Initial states:
l = 1500.1718
b = 39.0147
sigma: 455.9848
AIC AICc BIC
7128.496 7128.714 7152.308
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 45.93341 453.0599 288.3275 -0.1409079 8.483295 1.958599 0.01724913
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 5925.458 5341.090 6509.826 5031.744 6819.171
Sep 2024 5941.323 5356.522 6526.123 5046.948 6835.698
Oct 2024 5955.822 5370.131 6541.514 5060.085 6851.560
Nov 2024 5969.074 5381.931 6556.218 5071.116 6867.033
Dec 2024 5981.186 5391.963 6570.409 5080.047 6882.325
Jan 2025 5992.255 5400.286 6584.224 5086.916 6897.593
Feb 2025 6002.371 5406.974 6597.768 5091.790 6912.953
Mar 2025 6011.617 5412.113 6611.122 5094.754 6928.480
Apr 2025 6020.067 5415.793 6624.341 5095.910 6944.225
May 2025 6027.790 5418.111 6637.469 5095.367 6960.214
Jun 2025 6034.849 5419.165 6650.532 5093.242 6976.455
Jul 2025 6041.299 5419.052 6663.547 5089.654 6992.945 ME RMSE MAE MPE MAPE MASE ACF1
Training set 45.93341 453.0599 288.3275 -0.1409079 8.483295 1.958599 0.01724913# Holt exponential (multiplicative) damped trend
fit4 <- holt(sales.train.ts,
exponential = TRUE,
damped = TRUE,
h = 12)
summary(fit4)
Forecast method: Damped Holt's method with exponential trend
Model Information:
Damped Holt's method with exponential trend
Call:
holt(y = sales.train.ts, h = 12, damped = TRUE, exponential = TRUE)
Smoothing parameters:
alpha = 0.0204
beta = 0.0204
phi = 0.9041
Initial states:
l = 1527.3352
b = 1.0198
sigma: 0.1362
AIC AICc BIC
7038.291 7038.510 7062.103
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 46.92706 453.3145 288.3949 -0.0825909 8.472369 1.959057 0.0171936
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 5928.279 4893.998 6946.561 4361.387 7508.470
Sep 2024 5945.000 4899.169 7000.031 4365.605 7631.198
Oct 2024 5960.157 4892.634 7011.720 4331.739 7519.816
Nov 2024 5973.893 4896.469 7020.690 4367.942 7568.188
Dec 2024 5986.339 4920.628 7015.249 4360.534 7611.775
Jan 2025 5997.613 4939.476 7044.425 4384.034 7618.656
Feb 2025 6007.823 4964.482 7101.966 4392.766 7685.402
Mar 2025 6017.069 4923.645 7106.157 4372.110 7658.595
Apr 2025 6025.440 4961.296 7067.327 4406.143 7594.932
May 2025 6033.018 4926.393 7134.589 4382.143 7717.067
Jun 2025 6039.877 4951.451 7142.745 4390.314 7771.460
Jul 2025 6046.084 4929.498 7198.133 4354.391 7787.227 ME RMSE MAE MPE MAPE MASE ACF1
Training set 46.92706 453.3145 288.3949 -0.0825909 8.472369 1.959057 0.0171936All three Holt models improve on the simple exponential smoothing approach by capturing the upward trend in monthly retail alcohol sales. The standard Holt linear trend model provides a moderate improvement, with an RMSE of 454.47 and a MAPE of 8.57%, but it still leaves fairly large errors because it cannot represent the strong seasonal structure in the data.
The additive damped Holt model performs slightly better, lowering the RMSE to 453.06 and reducing the MAPE to 8.48%. The damping parameter (φ ≈ 0.91) slightly slows the trend over time, which helps the model follow the data more closely, but the improvement is modest.
The multiplicative damped Holt model yields the lowest MAPE among the three (8.47%) and similar RMSE (453.31), indicating that allowing the trend to grow proportionally rather than additively offers a small advantage. However, all Holt variations produce very similar error levels, and none of them match the accuracy expected from models that explicitly include seasonality.
Overall, the Holt models capture trend effectively but still fall short because they cannot model the repeated seasonal peaks and dips in the series. As a result, they are expected to be outperformed by the Holt–Winters seasonal models in later comparisons.
3.3 Holt-Winters Model
While the Holt model is able to capture the underlying trend in the retail sales data, it still lacks a seasonal component, which limits its ability to fully represent the strong repeating patterns observed in this monthly series. To address this, the Holt-Winters family of models incorporates both trend and seasonality, allowing the forecasts to adjust to regular fluctuations across the year. In the following subsection, we fit four Holt-Winters variations—additive, multiplicative, and their damped counterparts—to determine which seasonal structure best reflects the behavior of U.S. retail alcohol sales and provides the most accurate forecasts.
# Additive Holt-Winters model
fit5 <- hw(sales.train.ts,
h = 12,
seasonal = "additive")
summary(fit5)
Forecast method: Holt-Winters' additive method
Model Information:
Holt-Winters' additive method
Call:
hw(y = sales.train.ts, h = 12, seasonal = "additive")
Smoothing parameters:
alpha = 0.3574
beta = 0.0046
gamma = 0.5626
Initial states:
l = 1798.0845
b = 0.5834
s = 1129.675 42.5114 -60.7972 -64.8807 80.9583 156.8818
31.0221 80.4876 -228.7128 -180.6954 -489.6278 -496.8218
sigma: 112.0801
AIC AICc BIC
6041.853 6043.494 6109.321
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 8.312837 109.763 82.16355 0.2179468 2.746772 0.5581343 -0.06640942
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 6035.333 5891.697 6178.970 5815.660 6255.006
Sep 2024 5944.454 5791.692 6097.215 5710.825 6178.082
Oct 2024 5954.666 5793.079 6116.252 5707.541 6201.791
Nov 2024 6323.481 6153.320 6493.642 6063.242 6583.720
Dec 2024 7933.441 7754.917 8111.965 7660.412 8206.469
Jan 2025 5072.617 4885.912 5259.322 4787.076 5358.158
Feb 2025 5320.824 5126.093 5515.554 5023.009 5618.639
Mar 2025 5847.534 5644.913 6050.154 5537.652 6157.415
Apr 2025 5741.624 5531.231 5952.017 5419.855 6063.392
May 2025 6411.535 6193.473 6629.597 6078.038 6745.032
Jun 2025 6331.304 6105.664 6556.944 5986.218 6676.391
Jul 2025 6446.920 6213.782 6680.057 6090.366 6803.473 ME RMSE MAE MPE MAPE MASE ACF1
Training set 8.312837 109.763 82.16355 0.2179468 2.746772 0.5581343 -0.06640942# Multiplicative Holt-Winters model
fit6 <- hw(sales.train.ts,
h = 12,
seasonal = "multiplicative")
summary(fit6)
Forecast method: Holt-Winters' multiplicative method
Model Information:
Holt-Winters' multiplicative method
Call:
hw(y = sales.train.ts, h = 12, seasonal = "multiplicative")
Smoothing parameters:
alpha = 0.37
beta = 0.0036
gamma = 0.2244
Initial states:
l = 1691.0689
b = 2.8371
s = 1.4039 1.0185 0.9963 0.9622 1.0181 1.0619
0.9859 0.9914 0.9282 0.9078 0.8626 0.863
sigma: 0.0277
AIC AICc BIC
5805.827 5807.467 5873.295
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 7.290776 93.92774 69.49604 0.208493 2.174438 0.4720843 -0.1501059
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 6065.625 5850.595 6280.655 5736.764 6394.485
Sep 2024 5859.936 5638.235 6081.637 5520.874 6198.998
Oct 2024 5957.132 5718.218 6196.046 5591.744 6322.519
Nov 2024 6268.488 6003.424 6533.551 5863.108 6673.868
Dec 2024 7990.113 7635.465 8344.761 7447.726 8532.500
Jan 2025 5155.568 4916.255 5394.882 4789.570 5521.566
Feb 2025 5260.904 5006.321 5515.487 4871.552 5650.255
Mar 2025 5886.849 5590.669 6183.030 5433.880 6339.819
Apr 2025 5781.146 5479.448 6082.844 5319.739 6242.553
May 2025 6377.419 6032.916 6721.922 5850.548 6904.290
Jun 2025 6314.623 5962.177 6667.069 5775.603 6853.643
Jul 2025 6478.886 6105.861 6851.912 5908.393 7049.380 ME RMSE MAE MPE MAPE MASE ACF1
Training set 7.290776 93.92774 69.49604 0.208493 2.174438 0.4720843 -0.1501059# Additive Holt-Winters model with damping
fit7 <- hw(sales.train.ts,
h = 12,
seasonal = "additive",
damped = TRUE)
summary(fit7)
Forecast method: Damped Holt-Winters' additive method
Model Information:
Damped Holt-Winters' additive method
Call:
hw(y = sales.train.ts, h = 12, seasonal = "additive", damped = TRUE)
Smoothing parameters:
alpha = 0.35
beta = 0.0115
gamma = 0.542
phi = 0.98
Initial states:
l = 1819.6951
b = -3.4086
s = 1138.32 78.6614 -40.9844 -122.0709 45.8019 161.7
41.0247 83.0873 -201.6927 -182.9848 -493.5815 -507.2805
sigma: 113.4338
AIC AICc BIC
6052.197 6054.036 6123.634
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set 12.60166 110.9404 83.33934 0.3361054 2.789608 0.5661213
ACF1
Training set -0.05389815
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 6016.537 5871.166 6161.909 5794.211 6238.864
Sep 2024 5915.211 5760.645 6069.778 5678.822 6151.600
Oct 2024 5922.545 5758.778 6086.313 5672.085 6173.006
Nov 2024 6277.807 6104.829 6450.785 6013.260 6542.354
Dec 2024 7875.396 7693.196 8057.596 7596.745 8154.047
Jan 2025 5017.234 4825.799 5208.669 4724.459 5310.008
Feb 2025 5249.625 5048.942 5450.309 4942.707 5556.544
Mar 2025 5772.831 5562.885 5982.776 5451.747 6093.914
Apr 2025 5662.927 5443.707 5882.147 5327.659 5998.194
May 2025 6315.292 6086.785 6543.799 5965.821 6664.763
Jun 2025 6230.131 5992.326 6467.935 5866.440 6593.821
Jul 2025 6337.977 6090.864 6585.089 5960.051 6715.903 ME RMSE MAE MPE MAPE MASE
Training set 12.60166 110.9404 83.33934 0.3361054 2.789608 0.5661213
ACF1
Training set -0.05389815# Multiplicative Holt-Winters model with damping
fit8 <- hw(sales.train.ts,
h = 12,
seasonal = "multiplicative",
damped = TRUE)
summary(fit8)
Forecast method: Damped Holt-Winters' multiplicative method
Model Information:
Damped Holt-Winters' multiplicative method
Call:
hw(y = sales.train.ts, h = 12, seasonal = "multiplicative", damped = TRUE)
Smoothing parameters:
alpha = 0.3819
beta = 0.014
gamma = 0.2121
phi = 0.98
Initial states:
l = 1691.7117
b = 0.6692
s = 1.3902 1.0177 1.0029 0.9666 1.0074 1.0464
0.9875 1.0021 0.9371 0.9255 0.8496 0.8669
sigma: 0.0277
AIC AICc BIC
5806.177 5808.016 5877.614
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 9.640629 94.19309 70.05605 0.285802 2.187223 0.4758884 -0.1683167
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2024 6049.702 5835.237 6264.167 5721.706 6377.698
Sep 2024 5832.631 5610.272 6054.990 5492.563 6172.700
Oct 2024 5924.748 5683.109 6166.387 5555.193 6294.304
Nov 2024 6224.343 5953.965 6494.720 5810.836 6637.849
Dec 2024 7930.380 7564.888 8295.873 7371.408 8489.353
Jan 2025 5111.643 4862.535 5360.752 4730.665 5492.622
Feb 2025 5206.584 4939.083 5474.085 4797.476 5615.692
Mar 2025 5819.929 5505.533 6134.325 5339.102 6300.757
Apr 2025 5707.464 5384.060 6030.868 5212.860 6202.068
May 2025 6287.010 5914.151 6659.868 5716.772 6857.247
Jun 2025 6215.693 5830.636 6600.750 5626.800 6804.586
Jul 2025 6370.842 5959.334 6782.349 5741.495 7000.188 ME RMSE MAE MPE MAPE MASE ACF1
Training set 9.640629 94.19309 70.05605 0.285802 2.187223 0.4758884 -0.1683167The Holt–Winters models dramatically improve the fit compared to the SES and Holt models by explicitly capturing both trend and seasonality in monthly retail alcohol sales. The additive Holt–Winters model already reduces the training RMSE to about 109.76 and MAPE to 2.75%, which is a huge improvement over the non-seasonal models. However, the multiplicative Holt–Winters model performs even better, with the lowest RMSE (93.93) and lowest MAPE (2.17%), indicating that modeling the seasonal pattern as proportional to the level of the series is more appropriate for this data.
The damped Holt–Winters versions do not provide meaningful gains. The additive damped model has slightly higher error (MAPE ≈ 2.79%), and the damped multiplicative model is very close to, but slightly worse than, the non-damped multiplicative model (MAPE ≈ 2.19%). Overall, the Holt–Winters multiplicative model without damping offers the best in-sample fit among all smoothing models considered, and it is the strongest candidate for forecasting the monthly retail sales series.
3.4 Accuracy Comparison of All Smoothing Models
After fitting all eight exponential smoothing models, we next compare their in-sample performance using standard accuracy measures. This summary table allows us to evaluate how well each model captures the underlying structure of the training data and to determine which smoothing approach provides the best overall fit. By examining metrics such as RMSE, MAE, and MAPE across all models, we can identify the most effective candidate before moving on to forecasting and evaluating performance on the held-out test data.
accuracy.table <- round(
rbind(
accuracy(fit1), # SES
accuracy(fit2), # Holt Linear
accuracy(fit3), # Holt Add. Damped
accuracy(fit4), # Holt Exp. Damped
accuracy(fit5), # HW Additive
accuracy(fit6), # HW Multiplicative
accuracy(fit7), # HW Additive Damped
accuracy(fit8) # HW Multiplicative Damped
),
4
)
row.names(accuracy.table) <- c(
"SES",
"Holt Linear",
"Holt Add. Damped",
"Holt Exp. Damped",
"HW Add.",
"HW Mult.",
"HW Add. Damp",
"HW Mult. Damp"
)
kable(accuracy.table,
caption = "Accuracy Measures for Exponential Smoothing Models (Training Data)")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| SES | 85.7464 | 468.6373 | 301.8894 | 1.0332 | 8.8729 | 2.0507 | -0.0297 |
| Holt Linear | 28.9993 | 454.4737 | 289.1740 | -0.6700 | 8.5720 | 1.9643 | 0.0072 |
| Holt Add. Damped | 45.9334 | 453.0599 | 288.3275 | -0.1409 | 8.4833 | 1.9586 | 0.0172 |
| Holt Exp. Damped | 46.9271 | 453.3145 | 288.3949 | -0.0826 | 8.4724 | 1.9591 | 0.0172 |
| HW Add. | 8.3128 | 109.7630 | 82.1636 | 0.2179 | 2.7468 | 0.5581 | -0.0664 |
| HW Mult. | 7.2908 | 93.9277 | 69.4960 | 0.2085 | 2.1744 | 0.4721 | -0.1501 |
| HW Add. Damp | 12.6017 | 110.9404 | 83.3393 | 0.3361 | 2.7896 | 0.5661 | -0.0539 |
| HW Mult. Damp | 9.6406 | 94.1931 | 70.0561 | 0.2858 | 2.1872 | 0.4759 | -0.1683 |
The accuracy results for all eight smoothing models show a dramatic improvement when moving from non-seasonal to seasonal methods. The simple exponential smoothing (SES) and Holt models have relatively large errors, with MAPE values between 8.47% and 8.87%, confirming that models without a seasonal component cannot adequately capture the strong seasonal structure in the retail sales data. Among the non-seasonal models, the Holt exponential damped model performs slightly best, but the improvement is minimal.
Once seasonality is introduced, the Holt-Winters models outperform all other approaches by a wide margin. The additive Holt-Winters model reduces the RMSE to 109.76 and MAPE to 2.75%, while the multiplicative Holt-Winters model performs best overall with the lowest RMSE (93.93) and lowest MAPE (2.17%) of all eight models. The damped Holt-Winters versions do not improve accuracy, and in both the additive and multiplicative cases produce slightly larger errors than their non-damped counterparts.
Overall, the Holt-Winters multiplicative model provides the best in-sample performance and is the strongest candidate for forecasting out-of-sample values.
4 Forecast Plots
To visually compare the performance of the smoothing models, we next examine forecast plots for both the non-seasonal and seasonal approaches. These plots allow us to assess how well each model follows the observed training data and how effectively they project into the 12-month forecast horizon. By overlaying model predictions on the historical series, we can clearly see the limitations of simpler methods and the improvements gained by including trend and seasonal components, providing an intuitive complement to the numerical accuracy measures reported earlier.
4.1 Non-Seasonal Smoothing Models
The first set of forecast plots focuses on the non-seasonal smoothing models, which include the simple exponential smoothing method and the various Holt trend-based models. These methods account for the overall level and, in the case of the Holt models, the underlying trend in the series, but they do not incorporate any seasonal structure. Plotting these forecasts against the historical data allows us to see how well these simpler models track the upward movement in retail sales and highlights the limitations that arise when seasonality is ignored.
# 4.1 Non-seasonal smoothing models plot
n.train <- length(train.sales)
pred.id <- (n.train + 1):(n.train + h)
# Set y-limits to cover training data, forecasts, and test data
y.min <- min(c(train.sales, test.sales,
fit1$mean, fit2$mean, fit3$mean, fit4$mean))
y.max <- max(c(train.sales, test.sales,
fit1$mean, fit2$mean, fit3$mean, fit4$mean))
plot(1:n.train, train.sales,
type = "o", lwd = 1.5, cex = 0.6,
xlim = c(1, n.train + h),
ylim = c(y.min, y.max),
xlab = "Time (Months)",
ylab = "Retail Sales (Millions of Dollars)",
main = "Non-Seasonal Smoothing Models")
# Forecast lines for SES and Holt models
lines(pred.id, fit1$mean, col = "red", lwd = 1.5)
lines(pred.id, fit2$mean, col = "blue", lwd = 1.5)
lines(pred.id, fit3$mean, col = "purple", lwd = 1.5)
lines(pred.id, fit4$mean, col = "darkgreen", lwd = 1.5)
# Forecast points
points(pred.id, fit1$mean, col = "red", pch = 16, cex = 0.7)
points(pred.id, fit2$mean, col = "blue", pch = 17, cex = 0.7)
points(pred.id, fit3$mean, col = "purple", pch = 19, cex = 0.7)
points(pred.id, fit4$mean, col = "darkgreen", pch = 15, cex = 0.7)
# Actual test data
points(pred.id, test.sales,
col = "black", pch = 4, cex = 0.8)
legend("topleft",
legend = c("Training Data",
"SES",
"Holt Linear",
"Holt Add. Damped",
"Holt Exp. Damped",
"Actual Test Data"),
col = c("black", "red", "blue", "purple", "darkgreen", "black"),
pch = c(1, 16, 17, 19, 15, 4),
lty = c(1, 1, 1, 1, 1, NA),
bty = "n",
cex = 0.8)
The non-seasonal smoothing models all capture the overall upward trend in retail alcohol sales, but none of them are able to track the repeating seasonal fluctuations present in the data. The simple exponential smoothing model produces a very smoothed forecast that underestimates the upward momentum of the most recent observations. The Holt models perform noticeably better because they incorporate a trend component, and both the additive and damped versions follow the training data more closely than SES. However, all four non-seasonal models still produce forecasts that smooth over the seasonal peaks and dips, causing their projected values to fall short of the true variability seen in the series. This visual pattern reinforces the accuracy results: models without seasonality cannot adequately represent monthly retail sales and are expected to perform worse than the Holt–Winters models in the next section.
4.2 Holt-Winters Seasonal Smoothing Models Plot
The Holt-Winters models extend the non-seasonal approaches by incorporating both trend and seasonality, allowing them to more accurately reflect the recurring patterns present in the monthly retail sales data. Plotting these seasonal models alongside the training series provides a clear visual comparison of how well each specification captures the annual fluctuations and overall growth of the series. By examining these forecasts, we can see the benefits of modeling seasonality directly and identify which Holt-Winters variation provides the best alignment with the observed data.
# 4.2 Holt-Winters seasonal smoothing models plot (refined)
n.train <- length(train.sales)
pred.id <- (n.train + 1):(n.train + h)
# Set y-limits to cover training data, forecasts, and test data
y.min <- min(c(train.sales, test.sales,
fit5$mean, fit6$mean, fit7$mean, fit8$mean))
y.max <- max(c(train.sales, test.sales,
fit5$mean, fit6$mean, fit7$mean, fit8$mean))
plot(1:n.train, train.sales,
type = "o",
lwd = 1.2,
cex = 0.4, # smaller points for training data
xlim = c(1, n.train + h),
ylim = c(y.min, y.max),
xlab = "Time (Months)",
ylab = "Retail Sales (Millions of Dollars)",
main = "Holt-Winters Trend and Seasonal Models")
# Forecast lines for Holt-Winters models
lines(pred.id, fit5$mean, col = "red", lwd = 1.6) # HW Add.
lines(pred.id, fit6$mean, col = "blue", lwd = 1.6) # HW Mult.
lines(pred.id, fit7$mean, col = "purple", lwd = 1.6) # HW Add. Damp
lines(pred.id, fit8$mean, col = "darkgreen", lwd = 1.6) # HW Mult. Damp
# Forecast points
points(pred.id, fit5$mean, col = "red", pch = 16, cex = 0.7)
points(pred.id, fit6$mean, col = "blue", pch = 17, cex = 0.7)
points(pred.id, fit7$mean, col = "purple", pch = 19, cex = 0.7)
points(pred.id, fit8$mean, col = "darkgreen", pch = 15, cex = 0.7)
# Actual test data
points(pred.id, test.sales,
col = "black", pch = 4, cex = 0.8)
legend("topleft",
legend = c("Training Data",
"HW Add.",
"HW Mult.",
"HW Add. Damp",
"HW Mult. Damp",
"Actual Test Data"),
col = c("black", "red", "blue", "purple", "darkgreen", "black"),
pch = c(1, 16, 17, 19, 15, 4),
lty = c(1, 1, 1, 1, 1, NA),
bty = "n",
cex = 0.8)
The Holt-Winters models provide a much closer match to the structure of the retail sales data than the non-seasonal methods. All four versions successfully capture both the long-term upward trend and the recurring seasonal peaks present throughout the series. The additive and multiplicative models track the training data especially well, aligning closely with the monthly fluctuations that the SES and Holt models were unable to represent. The forecasts extend this seasonal pattern into the 12-month hold-out period, showing smooth and realistic seasonal waves rather than the flattened predictions produced by trend-only models. Among the four models, the multiplicative versions appear to fit slightly more tightly during months with larger seasonal swings, which is consistent with the multiplicative structure allowing the seasonal effect to scale with the level of the series. Overall, the Holt-Winters forecasts visually demonstrate a superior fit, reinforcing the accuracy results in the previous section and indicating that seasonal models are far more appropriate for this type of monthly retail data.
5 Testing Accuracy
To evaluate how well each exponential smoothing model performs on unseen data, the final twelve months of the series were held out and used as a test set. Unlike the training accuracy measures, which reflect how closely each model fits data it has already seen, the testing accuracy provides a more realistic assessment of forecasting performance. By comparing the models’ errors on the withheld observations, we can determine which smoothing method generalizes best and produces the most reliable forecasts for future retail sales. The following subsection computes MSE and MAPE for all eight models using the holdout period and summarizes their comparative performance.
5.1 Testing Accuracy Table
Using the held-out 12 months of data, we calculated out-of-sample forecasting errors for each of the eight smoothing models. Because these observations were not used during model fitting, the resulting metrics provide an unbiased assessment of predictive performance. Following the structure of the case study, we report two key measures for comparison—mean squared error (MSE) and mean absolute percentage error (MAPE). These values allow us to identify which models most accurately capture the true patterns in the retail sales data when forecasting beyond the training period.
# Function to compute MSE and MAPE on test data
acc.fun <- function(test.data, mod.obj) {
PE <- 100 * (test.data - mod.obj$mean) / mod.obj$mean
MAPE <- mean(abs(PE))
E <- test.data - mod.obj$mean
MSE <- mean(E^2)
c(MSE = MSE, MAPE = MAPE)
}
# Apply to all eight models
pred.accuracy <- rbind(
SES = acc.fun(test.sales, fit1),
Holt.Linear = acc.fun(test.sales, fit2),
Holt.Add.Dmp = acc.fun(test.sales, fit3),
Holt.Exp.Dmp = acc.fun(test.sales, fit4),
HW.Add = acc.fun(test.sales, fit5),
HW.Mult = acc.fun(test.sales, fit6),
HW.Add.Dmp2 = acc.fun(test.sales, fit7),
HW.Mult.Dmp2 = acc.fun(test.sales, fit8)
)
pred.accuracy <- round(pred.accuracy, 4)
kable(pred.accuracy,
caption = "Forecast Accuracy (MSE and MAPE) for the 12-Month Holdout Set")| MSE | MAPE | |
|---|---|---|
| SES | 507089.28 | 8.5308 |
| Holt.Linear | 605816.36 | 9.5086 |
| Holt.Add.Dmp | 512087.63 | 8.6187 |
| Holt.Exp.Dmp | 512242.54 | 8.6158 |
| HW.Add | 40707.59 | 2.8287 |
| HW.Mult | 39531.50 | 2.9293 |
| HW.Add.Dmp2 | 25086.73 | 2.1821 |
| HW.Mult.Dmp2 | 24150.35 | 2.2911 |
The testing accuracy table summarizes how well each smoothing model predicts the 12 months of data that were intentionally withheld from the training process. The SES and Holt models produce relatively large errors, with MSE values above 500,000 and MAPE values between 8% and 10%. These high error levels indicate that models without seasonality struggle to forecast monthly retail sales, which contain strong seasonal patterns.
In contrast, all Holt-Winters models achieve much lower forecast errors. The additive and multiplicative Holt-Winters methods reduce MSE to around 40,000 and bring MAPE down to roughly 2–3%, demonstrating a substantial improvement in predictive accuracy. The damped Holt-Winters variants perform even better, with the additive damped model achieving the lowest MAPE (2.18%) and one of the lowest MSE values. Overall, the testing accuracy results show a clear advantage for the Holt-Winters models, confirming that including both trend and seasonality is essential for accurately forecasting this dataset.
5.2 Identifying the Best Forecasting Model
Based on the testing accuracy results, the damped Holt-Winters additive model emerges as the best-performing approach for forecasting monthly retail alcohol sales. This model produced the lowest MAPE (2.18%) among all eight smoothing models, indicating that its percentage errors were smallest when predicting the holdout period. It also achieved one of the lowest MSE values, suggesting strong accuracy in absolute forecasting terms as well. The superior performance of this model is consistent with the structure of the dataset, which contains both a long-term upward trend and pronounced seasonal patterns. The damped trend component helps prevent the forecast from growing too aggressively into the future, while the additive seasonal component accurately reproduces the recurring monthly fluctuations in the data.
Overall, the testing results confirm that incorporating both trend and seasonality is essential for reliable forecasting, and the damped additive Holt-Winters method provides the most accurate and stable predictions for the 12-month horizon considered in this analysis.
6 Final Model Refit Using the Full Dataset
Once the best forecasting model has been identified, the last step is to refit that model using the entire time series. This updates the smoothing parameters using all available information and produces the final version of the model that would be used for real forecasting.
The damped additive Holt-Winters model was selected as the best performer, so we refit it below using the full series from 1992–2025.
sales.full.ts <- ts(sales$MRTSSM4453USN,
start = c(1992, 1),
frequency = 12)
final.model <- hw(sales.full.ts,
h = 12,
seasonal = "additive",
damped = TRUE)
# Extract smoothing parameters
final.params <- final.model$model$par[1:3]
kable(final.params,
caption = "Estimated Smoothing Parameters for the Final Holt-Winters Damped Additive Model")| x | |
|---|---|
| alpha | 0.3256501 |
| beta | 0.0111128 |
| gamma | 0.5634339 |
The estimated smoothing parameters show that the final damped Holt-Winters additive model updates the overall level of the series at a moderate rate (α = 0.3257), maintains a very stable and slowly changing trend component (β = 0.0111), and adjusts the seasonal pattern relatively quickly (γ = 0.5634). Together, these values indicate that the model responds appropriately to both long-term growth and strong recurring seasonal patterns in the retail sales data.
7 Conclusion
This analysis compared a comprehensive set of exponential smoothing models to forecast monthly U.S. retail alcohol sales. The results consistently showed that models lacking a seasonal component—such as simple exponential smoothing and the various Holt trend models—were unable to capture the strong recurring seasonal patterns present in the data, leading to relatively large forecast errors. Once seasonality was incorporated through the Holt-Winters framework, forecasting performance improved dramatically.
Among all eight models evaluated, the damped additive Holt-Winters method provided the most accurate predictions for the 12-month holdout period, producing the lowest MAPE and one of the lowest MSE values. Re-estimating this model on the full dataset yielded smoothing parameters consistent with the characteristics of the series: a moderately updated level, a slowly evolving trend, and a quickly adapting seasonal component. These features make the model well-suited for forecasting a series influenced by long-term economic growth and strong, recurring seasonal shifts.
Overall, the study demonstrates that incorporating both trend and seasonality is essential for accurately modeling monthly retail sales. The damped additive Holt-Winters model stands out as the most reliable and stable forecasting approach for this dataset and provides a strong foundation for future prediction and decision-making.