1 Introduction

For this week’s discussion I decided to use the aus_retail data set from the fpp3 libary. I decided to use that data set because I had not worked with it yet and I thought it would be interesting to explore a new data set and see the different patterns I can find. What I focus on from the data set is monthly turnover for Cafes, Restaurants and Takeaway Food Services in New South Wales.

2 Visualizing the data

retail <- aus_retail %>%
  filter(State == "New South Wales",
         Industry == "Cafes, restaurants and takeaway food services")

retail %>%
  autoplot(Turnover) +
  labs(title = "NSW Cafe/Restaurant/Takeaway Turnover",
       y = "Turnover ($ million)")

The plot shows a steady upward trend over time along with a repeating seasonal cycle which means that the turnover consistently peaks around the same months each year. Before fitting any models, we need to split the data into a training set and a held-out test set so that forecast accuracy can be evaluated on data the models have not seen.

retail_train <- retail %>% filter(Month <= yearmonth("2016 Dec"))
retail_test  <- retail %>% filter(Month > yearmonth("2016 Dec"))

3 Fitting the Models

retail_fit <- retail_train %>%
  model(
    SES  = ETS(Turnover ~ error("A") + trend("N") + season("N")),
    Holt = ETS(Turnover ~ error("A") + trend("A") + season("N")),
    HW   = ETS(Turnover ~ error("M") + trend("A") + season("M"))
  )

3.1 Smoothing Parameters

tidy(retail_fit)
## # A tibble: 23 × 5
##    State           Industry                                .model term  estimate
##    <chr>           <chr>                                   <chr>  <chr>    <dbl>
##  1 New South Wales Cafes, restaurants and takeaway food s… SES    alpha  6.04e-1
##  2 New South Wales Cafes, restaurants and takeaway food s… SES    l[0]   1.46e+2
##  3 New South Wales Cafes, restaurants and takeaway food s… Holt   alpha  5.77e-1
##  4 New South Wales Cafes, restaurants and takeaway food s… Holt   beta   1.00e-4
##  5 New South Wales Cafes, restaurants and takeaway food s… Holt   l[0]   1.40e+2
##  6 New South Wales Cafes, restaurants and takeaway food s… Holt   b[0]   2.90e+0
##  7 New South Wales Cafes, restaurants and takeaway food s… HW     alpha  7.81e-1
##  8 New South Wales Cafes, restaurants and takeaway food s… HW     beta   5.76e-3
##  9 New South Wales Cafes, restaurants and takeaway food s… HW     gamma  1.00e-4
## 10 New South Wales Cafes, restaurants and takeaway food s… HW     l[0]   1.42e+2
## # ℹ 13 more rows
  • Alpha (\(\alpha\)) controls how much weight is placed on the most recent observation when updating the level. A high alpha makes the model react quickly to recent changes; a low alpha smooths more heavily over history.
  • Beta (\(\beta\)), present in Holt’s and HW, controls how quickly the estimated trend/slope adapts to new information.
  • Gamma (\(\gamma\)), present only in HW, controls how quickly the seasonal component is allowed to change from year to year.

3.2 Model Summaries

report(retail_fit %>% select(SES))
## Series: Turnover 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.6036689 
## 
##   Initial states:
##      l[0]
##  146.0329
## 
##   sigma^2:  1655.062
## 
##      AIC     AICc      BIC 
## 5610.427 5610.485 5622.526
report(retail_fit %>% select(Holt))
## Series: Turnover 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.5770106 
##     beta  = 0.0001000184 
## 
##   Initial states:
##      l[0]     b[0]
##  140.1919 2.904238
## 
##   sigma^2:  1639.102
## 
##      AIC     AICc      BIC 
## 5608.371 5608.517 5628.537
report(retail_fit %>% select(HW))
## Series: Turnover 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.7812533 
##     beta  = 0.005759187 
##     gamma = 0.0001004155 
## 
##   Initial states:
##     l[0]      b[0]      s[0]     s[-1]    s[-2]    s[-3]    s[-4]    s[-5]
##  142.432 0.4738337 0.9947482 0.9275107 1.022848 1.142098 1.027138 1.023615
##      s[-6]     s[-7]     s[-8]     s[-9]    s[-10]    s[-11]
##  0.9884304 0.9852397 0.9828876 0.9390628 0.9860981 0.9803232
## 
##   sigma^2:  0.0015
## 
##      AIC     AICc      BIC 
## 4917.340 4918.874 4985.903

HW’s AIC is much lower than SES and Holt’s, which confirms that the seasonal component is capturing real,systematic structure in the data rather than fitting noise. SES and Holt’s are essentially tied on AIC, meaning the addition of a trend term alone barely improved in-sample fit once seasonality is left unmodeled.

Holt’s estimated beta is essentially 0, meaning the trend component is not updating from its initial estimate at all. The model has locked in a fixed slope of about 2.9 units/month and extrapolates it linearly for the entire forecast horizon. This is worth mentioning, since a frozen trend estimate is a risk when forecasting far into the future. Which again demonstrate that this type of modeling is better for short and medium time horizons.

4 Forecast Visualization

retail_fc <- retail_fit %>% forecast(h = "2 years")

retail_fc %>%
  autoplot(retail, level = NULL) +
  labs(title = "SES vs Holt's vs Holt-Winters: NSW Cafe Turnover",
       y = "Turnover ($ million)")

retail_fc %>%
  autoplot(retail_test, level = NULL) +
  labs(title = "Forecast Comparison Over the Test Period")

SES produces a flat forecast since it has no mechanism for trend or seasonality. Holt’s forecast slopes upward, capturing the trend but missing the seasonal peaks and troughs. HW is the only model that reproduces both the upward trajectory and the repeating seasonal shape.

5 Accuracy Comparison

accuracy(retail_fc, retail_test) %>%
  select(.model, RMSE, MAE, MAPE) %>%
  arrange(RMSE)
## # A tibble: 3 × 4
##   .model  RMSE   MAE  MAPE
##   <chr>  <dbl> <dbl> <dbl>
## 1 HW      29.9  21.7  1.66
## 2 SES     98.4  80.1  6.33
## 3 Holt   112.  102.   8.05

HW is the clear winner on every metric,confirming that modeling both trend and seasonality is essential for this series.

Holt’s actually performs worse than SES on the held-out test set, despite explicitly modeling a trend that SES ignores. Holt’s beta was estimated at essentially zero, its trend component never adapts and instead extrapolates the initial slope in a straight line for the full 2-year forecast horizon.

6 Residual Diagnostics

retail_fit %>% select(HW) %>% gg_tsresiduals()

augment(retail_fit) %>%
  filter(.model == "HW") %>%
  features(.innov, ljung_box, lag = 24)
## # A tibble: 1 × 5
##   State           Industry                              .model lb_stat lb_pvalue
##   <chr>           <chr>                                 <chr>    <dbl>     <dbl>
## 1 New South Wales Cafes, restaurants and takeaway food… HW        34.1    0.0829

For the HW model, the Ljung-Box test on 24 lags gives a statistic of 34.1 with a p-value of 0.083. Since this is above the conventional 0.05 threshold, we fail to reject the null hypothesis that the residuals are uncorrelated which supports the conclusion that HW has adequately captured the systematic (trend and seasonal) structure in the series, leaving behind residuals that are reasonably close to white noise. The p-value isn’t especially large,so there may be a small amount of remaining structure the model isn’t fully capturing, which is consistent with real-world retail data that often has irregular effects (holiday timing shifts, promotions, etc.) beyond a clean seasonal cycle.

7 Scenario Analysis: Sudden Shifts and Disruptions

To examine how each model responds to an sudden, one-off shock, we put an artificial 60% drop in turnover for a single month into the training data and refit all three models.

shocked <- retail_train %>%
  mutate(Turnover = if_else(Month == yearmonth("2015 Jun"),
                             Turnover * 0.4,
                             Turnover))
shock_fit <- shocked %>%
  model(
    SES  = ETS(Turnover ~ error("A") + trend("N") + season("N")),
    Holt = ETS(Turnover ~ error("A") + trend("A") + season("N")),
    HW   = ETS(Turnover ~ error("M") + trend("A") + season("M"))
  )

shock_fit %>%
  forecast(h = "1 year") %>%
  autoplot(shocked, level = NULL) +
  labs(title = "Effect of a Sudden Shock on Each Model's Forecast")

7.1 Discussion

Exponential smoothing models update their level, trend, and seasonal states gradually, based on the smoothing parameters estimated above. This means a single random shock is partially absorbed into the smoothed state rather than recognized as an outlier. The model has no way to distinguish a genuine one-off shock from the start of a real shift in the series.

  • The shock can distort the estimated level (and, depending on timing, the seasonal factor for that month) for several periods afterward.
  • None of these models adapt instantly, the degree of lingering distortion depends on how high alpha (and gamma, for HW) is: higher values recover faster but are noisier in normal periods.
  • For datasets prone to irregular shocks such as demand disruptions, policy changes, and extreme weather practical improvements include:
    • Outlier detection/correction before model fitting, so a known anomaly doesn’t contaminate the estimated states.
    • Damped trend to prevent forecasts from over-extrapolating a trend that may have been distorted by a shock.
    • Dynamic regression / ARIMA with exogenous regressors to explicitly encode known disruptions for example a dummy variable for a lockdown period rather than relying on the smoothing states to infer them after the fact.
    • Shorter or rolling estimation windows so the model re-weights recent behavior more quickly following a known structural change.

8 Conclusion

Holt-Winters was clearly the best-fitting and best-forecasting model for this dataset, both by AIC and by out-of-sample accuracy. This is expected given that NSW cafe turnover has both a real upward trend and a genuine, repeating seasonal cycle which are features that only HW is equipped to model.

Holt’s Linear Trend model performed worse than plain SES on the test set, despite modeling more structure. This is a useful reminder that a more complex model is only better if its extra parameters are well-estimated and genuinely reflect the data-generating process. Otherwise the added structure can actively hurt forecast accuracy.

The shock scenario confirmed a broader limitation of all exponential smoothing methods.This is because they update level/trend/seasonal states gradually by fixed smoothing parameters, none of them can distinguish a genuine one-off shock from the start of a real shift in the series. The fake shock partially contaminates the estimated states for a period afterward, regardless of which model is used. For data prone to sudden disruptions, this suggests exponential smoothing should be paired with outlier detection/correction, a damped trend to avoid over-extrapolation, or exogenous regressors that explicitly flag known disruptions, rather than relying on the smoothing parameters alone to react appropriately.