About this document

This notebook builds and compares forecasting models for a monthly sales series from Statistics South Africa (series code MSV30000, January 1998 to June 2026, in thousands of rand). The goal is a twelve-month-ahead forecast (July 2026 to June 2027) and an honest view of how far it can be trusted.

The workflow is the same in every step: look at the data, hold back the last 12 months as a test set, fit fourteen candidate models to the earlier data only, keep one model per family (plus all the benchmarks), rank the retained models on the test set, check their residuals, and refit the winner to the full series for the final forecast.

Headline result. The seasonal naive benchmark (each month forecast to equal the same month last year) beat every more elaborate model on the 12-month test set, with a MASE of 0.443 and a MAPE of 2.09%. It is therefore the model used for the final forecast.

1. Data

The data comes from Stats SA Excel workbook in which each row is a series and each column named MO + month + year holds one month. The code below keeps the row where H03 == "MSV30000", turns the month columns into a long table, decodes the month and year from the column names (characters 3 to 4 and 5 to 8), and stores the result as a monthly tsibble.

The stopifnot() line is a safety check. It stops the run if there are gaps in the monthly index, any missing values, or any values that are zero or negative. Logs are taken later, so the last condition matters.

raw_data <- read_excel("Excel table from 1998.xlsx")

sales_ts <- raw_data |>
  filter(H03 == "MSV30000") |>
  select(starts_with("MO")) |>
  pivot_longer(cols = everything(), names_to = "month_col", values_to = "sales_value_r1000") |>
  mutate(year = as.integer(substr(month_col, 5, 8)),
         mth  = as.integer(substr(month_col, 3, 4)),
         month_index = yearmonth(paste(year, mth, sep = "-"))) |>
  select(month_index, sales_value_r1000) |>
  as_tsibble(index = month_index)

stopifnot(!has_gaps(sales_ts)$.gaps, !anyNA(sales_ts$sales_value_r1000), all(sales_ts$sales_value_r1000 > 0))
c(start = format(min(sales_ts$month_index)), end = format(max(sales_ts$month_index)), n = nrow(sales_ts))
##      start        end          n 
## "1998 Jan" "2026 Jun"      "342"

The three checks passed. The series starts in January 1998 and ends in June 2026, which is 342 monthly observations with no gaps. The values are in thousands of rand (R’000).

2. Exploration

Before fitting anything we look at the series from several angles. Each plot answers a specific question: how does the level change, are the swings constant in size, is there a seasonal pattern, and how strongly is each month related to earlier months?

p <- sales_ts |> autoplot(sales_value_r1000) + labs(y = "R'000"); p; save_plot(p, "time_plot")

p <- sales_ts |> autoplot(log(sales_value_r1000)) + labs(y = "log(R'000)"); p; save_plot(p, "log_plot")

p <- sales_ts |> gg_season(sales_value_r1000, labels = "both"); p; save_plot(p, "season_plot", h = 4.5)

p <- sales_ts |> gg_subseries(sales_value_r1000); p; save_plot(p, "subseries_plot")

p <- sales_ts |> gg_lag(sales_value_r1000, geom = "point", lags = 1:12); p; save_plot(p, "lag_plot", h = 6)

p <- sales_ts |> ACF(sales_value_r1000, lag_max = 36) |> autoplot(); p; save_plot(p, "acf_plot")

What the six plots show, in the order they appear:

  1. Time plot. Sales rise from roughly R40 billion a month in 1998 to roughly R300 billion by 2025. The seasonal swings grow with the level, which points to a multiplicative pattern. There is a sharp one-off drop in early 2020.
  2. Log plot. After taking logs the swings are roughly constant in size and growth is close to a straight line until about 2023, after which the series moves sideways. This is the reason the models below work on the log scale. Logs also keep forecasts positive and turn regression coefficients into approximate percentage effects.
  3. Seasonal plot. One line per year. The shape repeats: low in January, a peak in March, a dip in April, a build through mid-year, a second peak in October and November, and a fall in December.
  4. Sub-series plot. One panel per calendar month, with the blue line showing the month’s average. Every month trends upward, and the April panel contains the 2020 collapse as an isolated low point.
  5. Lag plots. Every panel is a tight upward line, so sales are strongly related to sales at each of the last twelve lags. Much of this is the trend.
  6. ACF. The autocorrelation starts near 0.97 and only falls to about 0.64 by lag 36, with small bumps at lags 12 and 24. This slow decay is typical of a non-stationary series and is the first sign that differencing will be needed for ARIMA.

The axis labels use scientific notation, so 3e+08 means 3 x 10^8 thousand rand, which is R300 billion.

# Descriptive only (full series). Not used in any evaluation.
stl_full <- sales_ts |> model(STL(log(sales_value_r1000) ~ season(window = "periodic"), robust = TRUE))
p <- components(stl_full) |> autoplot(); p; save_plot(p, "stl_plot", h = 5.5)

comp <- components(stl_full) |> as_tibble()
outliers <- comp |>
  filter(abs(remainder - median(remainder)) > 3 * IQR(remainder)) |>
  transmute(month = format(month_index), remainder)
outliers; save_tbl(outliers, "outliers")
## # A tibble: 13 × 2
##    month     remainder
##    <chr>         <dbl>
##  1 2008 Apr     0.122 
##  2 2008 May     0.0922
##  3 2008 Jun     0.123 
##  4 2008 Jul     0.168 
##  5 2008 Aug     0.138 
##  6 2008 Sept    0.109 
##  7 2008 Oct     0.115 
##  8 2020 Apr    -0.673 
##  9 2020 May    -0.357 
## 10 2020 Jun    -0.176 
## 11 2020 Jul    -0.0899
## 12 2021 Jul    -0.131 
## 13 2021 Oct    -0.0951

3. Train/test split and transformation

The models are evaluated on data that they have not seen, where the last 12 observations (July 2025 to June 2026) are used as the test data. Then the remaining data (January 1998 to June 2025, 330 months) is used to train the forecast models.

test_start <- max(sales_ts$month_index) - 11
train_data <- sales_ts |> filter(month_index <  test_start)
test_data  <- sales_ts |> filter(month_index >= test_start)
stopifnot(nrow(test_data) == 12)

# Guerrero lambda on TRAINING data only, to justify using log (lambda = 0)
train_data |> features(sales_value_r1000, features = guerrero)
## # A tibble: 1 × 1
##   lambda_guerrero
##             <dbl>
## 1          -0.131

The Guerrero method runs on the training data only and gives a Box-Cox parameter of -0.131. A value of 0 is the natural log, and -0.131 is close to 0, which means that the two transformation will behave similar. The log is preferred because the seasonal swings are proportional to the level, coefficients become percentage effects, and back-transformed forecasts cannot be negative.

4. Candidate models (training data)

Fourteen models are fitted to the training data, grouped into the following five families:

The function fit_specs() holds all the model definitions so the same specifications can be refitted to the full series at the end. Within the regression, ETS and ARIMA families the model with the lowest AICc is retained. AICc is only comparable inside a family fitted to the same response, so it is never used across families.

fit_specs <- function(dat) {
  dat |> model(
    # Benchmarks
    Mean   = MEAN(log(sales_value_r1000)),
    Naive  = NAIVE(log(sales_value_r1000)),
    SNaive = SNAIVE(log(sales_value_r1000)),
    Drift  = RW(log(sales_value_r1000) ~ drift()),
    # Regression (EDIT: knot must be justified from history alone)
    TSLM_Trend = TSLM(log(sales_value_r1000) ~ trend() + season()),
    TSLM_Knot  = TSLM(log(sales_value_r1000) ~ trend(knots = yearmonth("2020 Apr")) + season()),
    # ETS (additive on the log scale)
    ETS_AAA  = ETS(log(sales_value_r1000) ~ error("A") + trend("A")  + season("A")),
    ETS_AAdA = ETS(log(sales_value_r1000) ~ error("A") + trend("Ad") + season("A")),
    ETS_Auto = ETS(log(sales_value_r1000)),
    # Decomposition: STL, forecast the SA component, reseasonalise
    STL_ARIMA = decomposition_model(
      STL(log(sales_value_r1000) ~ season(window = "periodic"), robust = TRUE),
      ARIMA(season_adjust ~ PDQ(0, 0, 0))),
    STL_ETS = decomposition_model(
      STL(log(sales_value_r1000) ~ season(window = "periodic"), robust = TRUE),
      ETS(season_adjust ~ season("N"))),
    # ARIMA (EDIT the two manual ones after identification below; "0 +" drops the constant when d + D = 2)
    ARIMA_Auto = ARIMA(log(sales_value_r1000)),
    ARIMA_BJ1  = ARIMA(log(sales_value_r1000) ~ 0 + pdq(0, 1, 1) + PDQ(0, 1, 1)),
    ARIMA_BJ2  = ARIMA(log(sales_value_r1000) ~ 0 + pdq(1, 1, 0) + PDQ(0, 1, 1))
  )
}

fit_cand <- fit_specs(train_data)

fam <- list(
  Benchmark     = c("Mean", "Naive", "SNaive", "Drift"),
  Regression    = c("TSLM_Trend", "TSLM_Knot"),
  ETS           = c("ETS_AAA", "ETS_AAdA", "ETS_Auto"),
  Decomposition = c("STL_ARIMA", "STL_ETS"),
  ARIMA         = c("ARIMA_Auto", "ARIMA_BJ1", "ARIMA_BJ2")
)
stopifnot(all(unlist(fam) %in% names(fit_cand)))

# AICc: within-family comparison only
aicc_tbl <- function(cols) fit_cand |> select(all_of(cols)) |> glance() |> select(model = .model, AIC, AICc, BIC)
aicc_all <- bind_rows(lapply(c("Regression", "ETS", "ARIMA"), function(f) aicc_tbl(fam[[f]]) |> mutate(family = f)))
aicc_all; save_tbl(aicc_all, "aicc")
## # A tibble: 8 × 5
##   model          AIC    AICc     BIC family    
##   <chr>        <dbl>   <dbl>   <dbl> <chr>     
## 1 TSLM_Trend -1459.  -1458.  -1406.  Regression
## 2 TSLM_Knot  -1481.  -1480.  -1424.  Regression
## 3 ETS_AAA      -25.6   -23.6    39.0 ETS       
## 4 ETS_AAdA     -19.3   -17.1    49.1 ETS       
## 5 ETS_Auto     -32.4   -30.5    32.2 ETS       
## 6 ARIMA_Auto  -888.   -888.   -858.  ARIMA     
## 7 ARIMA_BJ1   -938.   -938.   -927.  ARIMA     
## 8 ARIMA_BJ2   -930.   -930.   -918.  ARIMA

Reading the information-criteria table within each family:

The numbers across families are on different scales and should not be compared with each other.

Parameter estimates

This chunk prints the fitted coefficients for the regression, ETS and ARIMA families. Only TSLM_Trend is printed for the regression, not TSLM_Knot. To see the knot model’s coefficients, run tidy(select(fit_cand, TSLM_Knot)) or report() on it.

tslm_coef <- tidy(select(fit_cand, TSLM_Trend)); tslm_coef; save_tbl(tslm_coef, "tslm_coef")
## # A tibble: 13 × 6
##    .model     term           estimate std.error statistic   p.value
##    <chr>      <chr>             <dbl>     <dbl>     <dbl>     <dbl>
##  1 TSLM_Trend (Intercept)    17.7     0.0226       781.   0        
##  2 TSLM_Trend trend()         0.00533 0.0000620     86.0  1.05e-221
##  3 TSLM_Trend season()year2   0.126   0.0286         4.39 1.56e-  5
##  4 TSLM_Trend season()year3   0.204   0.0286         7.10 8.00e- 12
##  5 TSLM_Trend season()year4   0.0857  0.0286         2.99 2.98e-  3
##  6 TSLM_Trend season()year5   0.163   0.0286         5.69 2.89e-  8
##  7 TSLM_Trend season()year6   0.182   0.0287         6.34 8.04e- 10
##  8 TSLM_Trend season()year7   0.172   0.0289         5.94 7.61e-  9
##  9 TSLM_Trend season()year8   0.194   0.0289         6.72 8.51e- 11
## 10 TSLM_Trend season()year9   0.214   0.0289         7.39 1.31e- 12
## 11 TSLM_Trend season()year10  0.262   0.0289         9.06 1.36e- 17
## 12 TSLM_Trend season()year11  0.283   0.0289         9.78 6.34e- 20
## 13 TSLM_Trend season()year12  0.0917  0.0289         3.17 1.66e-  3
ets_par <- bind_rows(lapply(fam$ETS, function(nm)
  tidy(select(fit_cand, all_of(nm))) |> filter(term %in% c("alpha", "beta", "gamma", "phi")) |> mutate(model = nm)))
ets_par; save_tbl(ets_par, "ets_params")
## # A tibble: 10 × 4
##    .model   term  estimate model   
##    <chr>    <chr>    <dbl> <chr>   
##  1 ETS_AAA  alpha 0.646    ETS_AAA 
##  2 ETS_AAA  beta  0.000206 ETS_AAA 
##  3 ETS_AAA  gamma 0.000200 ETS_AAA 
##  4 ETS_AAdA alpha 0.686    ETS_AAdA
##  5 ETS_AAdA beta  0.00169  ETS_AAdA
##  6 ETS_AAdA gamma 0.000100 ETS_AAdA
##  7 ETS_AAdA phi   0.980    ETS_AAdA
##  8 ETS_Auto alpha 0.619    ETS_Auto
##  9 ETS_Auto beta  0.000101 ETS_Auto
## 10 ETS_Auto gamma 0.000103 ETS_Auto
arima_coef <- bind_rows(lapply(fam$ARIMA, function(nm) tidy(select(fit_cand, all_of(nm))) |> mutate(model = nm)))
arima_coef; save_tbl(arima_coef, "arima_coef")
## # A tibble: 11 × 7
##    .model     term     estimate std.error statistic  p.value model     
##    <chr>      <chr>       <dbl>     <dbl>     <dbl>    <dbl> <chr>     
##  1 ARIMA_Auto ar1        0.902    0.0344      26.3  1.24e-81 ARIMA_Auto
##  2 ARIMA_Auto ma1       -0.229    0.0646      -3.54 4.67e- 4 ARIMA_Auto
##  3 ARIMA_Auto ma2       -0.133    0.0626      -2.12 3.44e- 2 ARIMA_Auto
##  4 ARIMA_Auto ma3        0.107    0.0585       1.83 6.85e- 2 ARIMA_Auto
##  5 ARIMA_Auto sar1      -0.694    0.0545     -12.7  2.54e-30 ARIMA_Auto
##  6 ARIMA_Auto sar2      -0.321    0.0525      -6.11 2.91e- 9 ARIMA_Auto
##  7 ARIMA_Auto constant   0.0132   0.00241      5.47 9.24e- 8 ARIMA_Auto
##  8 ARIMA_BJ1  ma1       -0.345    0.0603      -5.72 2.51e- 8 ARIMA_BJ1 
##  9 ARIMA_BJ1  sma1      -0.979    0.105       -9.28 2.64e-18 ARIMA_BJ1 
## 10 ARIMA_BJ2  ar1       -0.242    0.0544      -4.45 1.18e- 5 ARIMA_BJ2 
## 11 ARIMA_BJ2  sma1      -1.000    0.300       -3.33 9.60e- 4 ARIMA_BJ2

Regression (TSLM_Trend). The fitted equation on the log scale is roughly 17.7 + 0.00533 t + seasonal effect, with January as the baseline month. The trend coefficient implies growth of about 0.53% a month, or around 6.6% a year. All seasonal effects are statistically significant at the 1% level. After allowing for the trend, November is about 33% above January and October about 30% above, while April and December are only about 9 to 10% above. Be careful with the p-values: the tests assume independent errors, and the regression residuals are strongly autocorrelated (see the diagnostics section), so the standard errors are too small and the t-statistics too large.

ETS. The level parameter alpha is large in all three models (0.62 to 0.69), so the level reacts quickly to new data. The trend and seasonal parameters (beta and gamma) are at or near the lower bound of about 0.0001, meaning the slope and the seasonal pattern are treated as nearly fixed. In the damped model phi is 0.980, a mild damping.

ARIMA.

  • ARIMA_BJ1 has an MA(1) coefficient of -0.345 and a seasonal MA(1) coefficient of -0.979, both strongly significant. The seasonal term sits almost on the boundary of -1, which is common when the seasonal pattern is very stable but also means the model is close to non-invertible.
  • ARIMA_BJ2 has an AR(1) coefficient of -0.242 and a seasonal MA(1) of -1.000 with a large standard error (0.300).
  • ARIMA_Auto has a non-seasonal AR(1) of 0.902, three MA terms, two seasonal AR terms and a constant of 0.0132. The AR coefficient near 1 shows the series sits close to the line between needing a further difference and not needing one.

Box–Jenkins identification

The manual ARIMA models are chosen the Box-Jenkins way: decide how many differences are needed, then read the autocorrelation (ACF) and partial autocorrelation (PACF) of the differenced series to pick the orders. Everything here uses the training data only.

The first table counts the seasonal and ordinary differences suggested by unit-root rules on the log series, and gives a KPSS test on the seasonally differenced series. The KPSS null hypothesis is that the series is stationary, so a large p-value means no evidence against stationarity.

train_log <- train_data |> mutate(log_sales = log(sales_value_r1000))

ur1 <- train_log |> features(log_sales, list(nsdiffs = unitroot_nsdiffs, ndiffs = unitroot_ndiffs))
ur2 <- train_log |> mutate(d12 = difference(log_sales, 12)) |> filter(!is.na(d12)) |>
  features(d12, list(kpss = unitroot_kpss, ndiffs_after_D1 = unitroot_ndiffs))
unitroot <- bind_cols(ur1, ur2); unitroot; save_tbl(unitroot, "unitroot")
## # A tibble: 1 × 5
##   nsdiffs_nsdiffs ndiffs_ndiffs kpss_kpss_stat kpss_kpss_pvalue
##             <int>         <int>          <dbl>            <dbl>
## 1               1             1          0.276              0.1
## # ℹ 1 more variable: ndiffs_after_D1_ndiffs <int>
p <- train_log |> gg_tsdisplay(difference(log_sales, 12), plot_type = "partial", lag_max = 36)
p; save_plot(p, "tsdisplay_D1", h = 5)

p <- train_log |> gg_tsdisplay(difference(difference(log_sales, 12)), plot_type = "partial", lag_max = 36)
p; save_plot(p, "tsdisplay_D1d1", h = 5)

How many differences? Both rules suggest one difference on the log series. After one seasonal difference the KPSS statistic is 0.276 with a p-value of 0.1 (the upper limit the test reports), so stationarity cannot be rejected. The first plot, though, shows an ACF that decays slowly from about 0.73 at lag 1 and a PACF dominated by a spike at lag 1, which is what an autoregressive coefficient close to one looks like. The tests and the plots therefore disagree a little. The two manual models follow standard practice and take both a first and a seasonal difference (\(d = D = 1\)). The last column of the unit-root table (further differences suggested after the seasonal difference) may be cut off in the printed output, so check it directly if you need it.

Which orders? After both differences the second plot shows:

  • one large negative ACF spike at lag 12 (about -0.5) with nothing similar at lag 24, and PACF spikes at lags 12 and 24 that shrink gradually. That points to a seasonal MA(1) term;
  • a negative ACF spike at lag 1 (about -0.2) and negative PACF values at lags 1 to 3 that decay. That is consistent with either a non-seasonal MA(1) or an AR(1).

These readings give the two candidates in fit_specs(): ARIMA_BJ1 is ARIMA(0,1,1)(0,1,1)[12], the classic airline model, and ARIMA_BJ2 is ARIMA(1,1,0)(0,1,1)[12]. The 0 + in each formula drops the constant, because with \(d + D = 2\) a constant would imply a quadratic trend. The 2020 and 2021 spikes show up in both plots, because differencing turns a one-off drop into a spike followed by a bounce.

5. Retained models and test-set accuracy

Now the retained models are chosen and tested. The code keeps all four benchmarks, both STL models, and the lowest-AICc model from the regression, ETS and ARIMA families. That gives nine retained models.

Every model forecasts 12 months from the same origin (June 2025), and accuracy() compares the forecasts with the held-out data. The measures are RMSE, MAE, MAPE and MASE. MASE divides the test MAE by the in-sample MAE of the seasonal naive method, so a value below 1 means the model does better than the usual in-sample error of “same month last year”. It has no units and does not blow up when values are small, which is why the models are ranked by it. The last column rescales each model’s MASE by that of SNaive.

# One model per family via AICc (valid within a family only; check d and D match for ARIMA).
# Both STL variants are kept (no comparable AICc). EDIT `keep` if you disagree.
best_by_aicc <- function(cols) aicc_tbl(cols) |> slice_min(AICc, n = 1) |> pull(model)
keep <- c(fam$Benchmark, best_by_aicc(fam$Regression), best_by_aicc(fam$ETS),
          fam$Decomposition, best_by_aicc(fam$ARIMA))
non_bench <- setdiff(keep, fam$Benchmark)
fit_all <- fit_cand |> select(all_of(keep))

fc_all <- fit_all |> forecast(h = 12)   # same origin and horizons 1-12 for every model

# accuracy() needs the FULL series so MASE is scaled by the in-sample seasonal naive MAE
acc <- fc_all |> accuracy(sales_ts) |>
  select(model = .model, RMSE, MAE, MASE, MAPE) |>
  mutate(MASE_rel_SNaive = MASE / MASE[model == "SNaive"]) |>
  arrange(MASE)
acc; save_tbl(acc, "accuracy")
## # A tibble: 9 × 6
##   model           RMSE        MAE   MASE  MAPE MASE_rel_SNaive
##   <chr>          <dbl>      <dbl>  <dbl> <dbl>           <dbl>
## 1 SNaive      7208302.   5913106.  0.443  2.09            1   
## 2 TSLM_Knot  13294214.  11066735.  0.828  3.89            1.87
## 3 STL_ARIMA  15193312.  12191521.  0.913  4.32            2.06
## 4 STL_ETS    15403076.  12339814.  0.924  4.37            2.09
## 5 ARIMA_BJ1  15881314.  12950068.  0.969  4.57            2.19
## 6 ETS_Auto   16065632.  13230832.  0.990  4.67            2.24
## 7 Naive      29221266.  23519331.  1.76   8.54            3.98
## 8 Drift      40515218.  32484060.  2.43  11.8             5.49
## 9 Mean      139129267. 137298061. 10.3   46.9            23.2
p <- fc_all |> autoplot(sales_ts |> filter_index("2022 Jan" ~ .), level = NULL)
p; save_plot(p, "test_forecasts", h = 4.5)

The ranking is clear:

Rank Model MASE MAPE (%)
1 SNaive 0.443 2.09
2 TSLM_Knot 0.828 3.89
3 STL_ARIMA 0.913 4.32
4 STL_ETS 0.924 4.37
5 ARIMA_BJ1 0.969 4.57
6 ETS_Auto 0.990 4.67
7 Naive 1.76 8.54
8 Drift 2.43 11.8
9 Mean 10.3 46.9

SNaive has the lowest error on every measure, and the models that estimate a trend are about twice as inaccurate. The most likely reason is the shape of the history. For about 25 years the log series rose steadily, and every trend-carrying model absorbed that growth. Since roughly 2023 the series has moved sideways, so these models start above the data and stay there. The forecast plot supports this reading: the trend-based models get the shape of the year right but sit above the observed line, and Drift keeps climbing. This explanation comes from reading the plot rather than from a month-by-month error table, so hold it with some caution. Mean and Naive do worst because one ignores both level and season and the other ignores the seasonal pattern.

6. Residual diagnostics (retained non-benchmark models)

Good point forecasts are not enough. If a model’s residuals still contain structure, it has missed something. For each retained non-benchmark model the chunk below:

The benchmarks are not put through these checks, so there are no diagnostics for SNaive.

# dof: ARIMA = p+q+P+Q, ETS = number of smoothing parameters, others 0.
# CHECK STL_ARIMA and TSLM by hand (read report()) and override, e.g. lb_dof["STL_ARIMA"] <- 3
dof_for <- function(nm) {
  tt <- tryCatch(tidy(select(fit_all, all_of(nm))), error = function(e) NULL)
  if (is.null(tt)) return(0L)
  if (grepl("^ARIMA", nm)) return(sum(grepl("^(ar|ma|sar|sma)[0-9]+$", tt$term)))
  if (grepl("^ETS", nm))   return(sum(tt$term %in% c("alpha", "beta", "gamma", "phi")))
  0L
}
lb_lag <- 24
lb_dof <- setNames(vapply(non_bench, dof_for, integer(1)), non_bench)

lb_all <- bind_rows(lapply(non_bench, function(nm)
  fit_all |> select(all_of(nm)) |> augment() |> filter(!is.na(.innov)) |>
    features(.innov, ljung_box, lag = lb_lag, dof = lb_dof[[nm]]) |>
    mutate(lag = lb_lag, dof = lb_dof[[nm]])))
lb_all; save_tbl(lb_all, "ljung_box")
## # A tibble: 5 × 5
##   .model    lb_stat lb_pvalue   lag   dof
##   <chr>       <dbl>     <dbl> <dbl> <int>
## 1 TSLM_Knot  1979.      0        24     0
## 2 ETS_Auto     24.8     0.254    24     3
## 3 STL_ARIMA    18.8     0.760    24     0
## 4 STL_ETS      22.0     0.576    24     0
## 5 ARIMA_BJ1    25.0     0.298    24     2
resid_sum <- fit_all |> select(all_of(non_bench)) |> augment() |> as_tibble() |>
  group_by(model = .model) |>
  summarise(n_missing = sum(is.na(.innov)), mean = mean(.innov, na.rm = TRUE), sd = sd(.innov, na.rm = TRUE),
            beyond_3sd = sum(abs(.innov - mean) > 3 * sd, na.rm = TRUE))
resid_sum; save_tbl(resid_sum, "resid_summary")
## # A tibble: 5 × 5
##   model     n_missing      mean      sd beyond_3sd
##   <chr>         <int>     <dbl>   <dbl>      <int>
## 1 ARIMA_BJ1         0 -1.47e- 3 0.0509           3
## 2 ETS_Auto          0 -8.27e- 6 0.00266          3
## 3 STL_ARIMA        12  1.92e- 4 0.0509           1
## 4 STL_ETS          12  2.90e- 5 0.0513           3
## 5 TSLM_Knot         0 -6.20e-18 0.101            4
for (nm in non_bench) {
  p <- gg_tsresiduals(select(fit_all, all_of(nm)))
  print(p); save_plot(p, paste0("resid_", nm), h = 4.5)
}

Ljung-Box results.

Model Q* p-value
TSLM_Knot 1979 < 0.001
ETS_Auto 24.8 0.254
STL_ARIMA 18.8 0.760
STL_ETS 22.0 0.576
ARIMA_BJ1 25.0 0.298

Things to check. The STL models were tested with zero degrees of freedom because the code cannot count the coefficients inside the decomposition model, which makes their p-values slightly generous. The comment in the chunk suggests inspecting them by hand with report(). The ETS_Auto residual standard deviation (about 0.0027) is much smaller than the others (about 0.05) because it appears to have multiplicative errors, so its residuals are relative errors and are not on the same scale. The two STL models have 12 missing innovation residuals, so their plots start about a year after the data do.

7. Final model and forecasts

The final model is the retained model with the lowest test MASE, unless it is overridden after weighing the diagnostics, parsimony and plausibility. The chosen specification is then refitted to all 342 observations, so the forecasts use the most recent data, and used to forecast July 2026 to June 2027.

The forecasts are back-transformed from the log scale, so two point forecasts are produced: the median (the plain back-transform, \(e^{\mu}\)) and the bias-adjusted mean (about \(e^{\mu}(1 + \sigma^2/2)\)), which fable reports by default. The 80% and 95% prediction intervals are built on the log scale and then back-transformed.

# Default = best MASE. EDIT after weighing diagnostics, parsimony and plausibility.
final_name <- acc$model[1]; final_name
## [1] "SNaive"
# Refit the same specification to all observations, forecast 12 months ahead
fit_full <- sales_ts |> fit_specs() |> select(all_of(final_name))
fc_final <- fit_full |> forecast(h = 12)

# Original scale. fable's default is bias_adjust = TRUE, so .mean is the bias-adjusted mean;
# `median` is the plain back-transform.
fc_tab <- fc_final |>
  mutate(median = median(sales_value_r1000)) |>
  hilo(level = c(80, 95)) |>
  unpack_hilo(c(`80%`, `95%`)) |>
  as_tibble() |>
  transmute(month = format(month_index), mean = .mean, median,
            lo80 = `80%_lower`, hi80 = `80%_upper`, lo95 = `95%_lower`, hi95 = `95%_upper`)
fc_tab; save_tbl(fc_tab, "final_forecasts")
## # A tibble: 12 × 7
##    month           mean     median       lo80       hi80       lo95       hi95
##    <chr>          <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
##  1 2026 Jul  298642331. 296413445  253304874. 346858427. 233083685. 376950152.
##  2 2026 Aug  300895567. 298649864. 255216042. 349475450. 234842285. 379794214.
##  3 2026 Sept 311995070. 309666527. 264630509. 362366979. 243505200. 393804148.
##  4 2026 Oct  327060349. 324619368. 277408699. 379864562. 255263315. 412819735.
##  5 2026 Nov  323462164. 321048038. 274356762. 375685447. 252455012. 408278061.
##  6 2026 Dec  262964117. 261001512. 223043038. 305419931. 205237634. 331916656.
##  7 2027 Jan  246728119. 244886690. 209271858. 286562616. 192565800. 311423373.
##  8 2027 Feb  274777634. 272726860. 233063122. 319140752. 214457821. 346827828.
##  9 2027 Mar  303573868. 301308176. 257487745. 352586166. 236932640. 383174800.
## 10 2027 Apr  282034622. 279929686  239218413. 327569388. 220121738. 355987690.
## 11 2027 May  289264485. 287105590. 245350697. 335966520. 225764485. 365113315.
## 12 2027 Jun  291399065. 289224239. 247161223. 338445730. 227430478. 367807610.
p <- fc_final |> autoplot(sales_ts |> filter_index("2022 Jan" ~ .), level = c(80, 95)) +
  labs(x = "Month", y = "R'000")
p; save_plot(p, "final_forecast", h = 4.5)

The selected model is SNaive. It won on accuracy, it is the simplest model on the list, and the runner-up (TSLM_Knot) fails its residual test. The remaining models have better known diagnostics but are about twice as inaccurate on the test year and each carries a trend the recent data do not support.

For a seasonal naive model the median forecast for each month is simply the value from the same month a year earlier, so the median column repeats the last twelve observations. The mean is about 0.75% higher. The forecasts show the same shape as the last year: a peak of about R325 billion in October 2026, about R321 billion in November, a drop to about R261 billion in December, and a low of about R245 billion in January 2027, before recovering to about R300 billion in March.

The intervals are wide. A 95% interval spans roughly plus or minus 27% around the median (for July 2026, about R233 billion to R377 billion). The implied standard deviation on the log scale is about 0.123, against about 0.051 for the ARIMA and STL residuals. Part of the width comes from the 2020 and 2021 shocks, which inflate the standard deviation of the seasonal differences, so the true uncertainty is probably smaller than the intervals suggest. This has not been checked.

Limitations to keep in mind.