DATA 624 Homework: Chapter 5 Exercises

Exercises 5.1, 5.2, 5.3, 5.4 and 5.7

library(fpp3)   # loads tsibble, tsibbledata, feasts, fable, dplyr, ggplot2, ...

# small helper I'll reuse a few times, it takes fitted model (mable) and
# runs the Ljung-Box test on innovation residuals, same as the book does
# for section 5.4. The book suggests lag = 10 for non-seasonal data and
# lag = 2m for seasonal data (so 8 for quarterly, 24 for monthly).
lb_test <- function(fit, lag) {
  augment(fit) |>          # augment() adds .fitted, .resid and .innov columns
    features(.innov, ljung_box, lag = lag)
}

Exercise 5.1 — Pick the right benchmark for each series

Produce forecasts for the following series using whichever of NAIVE(y), SNAIVE(y) or RW(y ~ drift()) is more appropriate in each case: Australian Population (global_economy), Bricks (aus_production), NSW Lambs (aus_livestock), Household wealth (hh_budget), Australian takeaway food turnover (aus_retail).

What I’m using from section 5.2 is if the data are seasonal, use SNAIVE(), if there’s no seasonality but a clear trend then use the drift method, and if it just wanders around with no trend and no season use NAIVE().

(a) Australian population (global_economy)

aus_pop <- global_economy |>
  filter(Country == "Australia") |>
  mutate(Pop_m = Population / 1e6)          # millions so the axis is readable

aus_pop |>
  autoplot(Pop_m) +
  labs(title = "Australian population", y = "Millions of people", x = "Year")

# annual data -> no seasonality, but a steady upward trend -> drift
aus_pop |>
  model(Drift = RW(Pop_m ~ drift())) |>
  forecast(h = 10) |>                        # 10 years ahead
  autoplot(aus_pop) +
  labs(title = "Australian population: 10 year drift forecast",
       y = "Millions of people", x = "Year")

Answer. Since the data is exclusively yearly, there’s no seasonality and since it keeps going up at an average historical rate lets stick with drift

(b) Bricks (aus_production)

bricks <- aus_production |>
  filter(!is.na(Bricks)) |>       # Bricks stops at 2005 Q2, the rest is NA
  select(Quarter, Bricks)

bricks |>
  autoplot(Bricks) +
  labs(title = "Australian quarterly clay brick production",
       y = "Millions of bricks", x = "Quarter")

# quarterly with a clear seasonal pattern -> seasonal naive
bricks |>
  model(SNAIVE(Bricks)) |>
  forecast(h = "3 years") |>
  autoplot(bricks |> filter_index("1995 Q1" ~ .)) +  # zoom in on recent years
  labs(title = "Bricks: 3 year seasonal naive forecast",
       y = "Millions of bricks", x = "Quarter")

Answer. Bricks is quarterly and there’s definitely a seasonal pattern. Trend is not clear here so in this case I would go with SNAIVE since it just repeats last year’s quarterly pattern.

(c) NSW Lambs (aus_livestock)

nsw_lambs <- aus_livestock |>
  filter(State == "New South Wales", Animal == "Lambs")

nsw_lambs |>
  autoplot(Count / 1e3) +
  labs(title = "Lambs slaughtered in New South Wales",
       y = "Thousands of head", x = "Month")

# monthly with a yearly pattern -> seasonal naive
nsw_lambs |>
  model(SNAIVE(Count)) |>
  forecast(h = "2 years") |>
  autoplot(nsw_lambs |> filter_index("2010 Jan" ~ .)) +
  labs(title = "NSW lambs: 2 year seasonal naive forecast",
       y = "Head", x = "Month")

Answer. There’s no real trend here but we do have monthly data that sort of repeats itself every year. We mainly see some very noisy up and downs over the years, out of the 3 SNAIVE would make the most sense here.

(d) Household wealth (hh_budget)

hh_budget |>
  autoplot(Wealth) +
  labs(title = "Household wealth (% of net disposable income)",
       y = "% of net disposable income", x = "Year")

# annual data, 4 countries (the key), each trending up -> drift
# model() fits one drift model per country automatically (Country is the key)
hh_budget |>
  model(Drift = RW(Wealth ~ drift())) |>
  forecast(h = 5) |>
  autoplot(hh_budget) +
  facet_wrap(vars(Country), scales = "free_y") +
  theme(legend.position = "none") +
  labs(title = "Household wealth: 5 year drift forecast by country",
       y = "% of net disposable income", x = "Year")

Answer. hh_budget is anual so no seasonality but we do see a long term upward trend so I would go with drift so it can keep that up instead of going flat like NAIVE() would.

(e) Australian takeaway food turnover (aus_retail)

# aus_retail is keyed by State, so I add up all the states to get the
# national total. (Northern Territory only starts in 1988 but it's tiny,
# so it doesn't really move the total.)
takeaway <- aus_retail |>
  filter(Industry == "Takeaway food services") |>
  summarise(Turnover = sum(Turnover))

takeaway |>
  autoplot(Turnover) +
  labs(title = "Australian takeaway food turnover",
       y = "$ million AUD", x = "Month")

# strong monthly seasonality -> seasonal naive
takeaway |>
  model(SNAIVE(Turnover)) |>
  forecast(h = "2 years") |>
  autoplot(takeaway |> filter_index("2012 Jan" ~ .)) +
  labs(title = "Takeaway food turnover: 2 year seasonal naive forecast",
       y = "$ million AUD", x = "Month")

Answer. Seasonality is obvious here with peaks at december, although there’s an upward trend I went with SNAIVE because it’s the only one that keeps the seasonal shape, so you might get a lower forecasting due to it just repeating and not following the upward trend.

Series	Frequency	Seasonal?	Trend?	Method
Australian population	Annual	No	Strong, steady	RW(y ~ drift())
Bricks	Quarterly	Yes	Up then down	SNAIVE(y)
NSW lambs	Monthly	Yes	Weak	SNAIVE(y)
Household wealth	Annual	No	Upward	RW(y ~ drift())
Takeaway food	Monthly	Yes	Upward	SNAIVE(y)

Exercise 5.2 — Facebook stock price (gafa_stock)

(a) Time plot

Stock prices were recorded only on trading days, so the date index has gaps for weekends and holidays. Same as the Google example in the book section 5.2 so I re-indexed them by trading day so the tsibble is regular and the benchmark methods work.

fb_stock <- gafa_stock |>
  filter(Symbol == "FB") |>
  mutate(day = row_number()) |>                       # trading day 1, 2, 3...
  update_tsibble(index = day, regular = TRUE)   # day is the new regular index

fb_stock |>
  autoplot(Close) +
  labs(title = "Facebook daily closing stock price (2014-2018)",
       y = "US$", x = "Trading day")

(b) Drift forecasts

fb_drift <- fb_stock |>
  model(Drift = RW(Close ~ drift()))

report(fb_drift)        # shows the estimated slope (drift) per day

## Series: Close 
## Model: RW w/ drift 
## 
## Drift: 0.0608 (se: 0.0681)
## sigma^2: 5.8301

fb_drift_fc <- fb_drift |>
  forecast(h = 60)                              # about 3 months of trading days

fb_drift_fc |>
  autoplot(fb_stock) +
  labs(title = "Facebook closing price: drift forecast (60 trading days)",
       y = "US$", x = "Trading day")

(c) The drift forecast is the line between the first and last observations

From section 5.2, the drift forecast is

\[\hat{y}_{T+h|T} = y_T + h\left(\frac{y_T - y_1}{T - 1}\right)\]

which is just the straight line through the first and last points extended into the future. I’ll compute that line by hand and check it against forecast.

n_obs   <- nrow(fb_stock)
y_first <- first(fb_stock$Close)
y_last  <- last(fb_stock$Close)
slope   <- (y_last - y_first) / (n_obs - 1)   # rise over run from first point to last

c(first = y_first, last = y_last, T = n_obs, slope = slope)

##        first         last            T        slope 
## 5.471000e+01 1.310900e+02 1.258000e+03 6.076372e-02

# line through (1, y_first) and (T, y_last), evaluated at the forecast days
line_check <- fb_drift_fc |>
  as_tibble() |>
  mutate(line = y_first + slope * (day - 1),
         difference = .mean - line) |>
  select(day, drift_forecast = .mean, line, difference)

head(line_check)

## # A tibble: 6 × 4
##     day drift_forecast  line difference
##   <dbl>          <dbl> <dbl>      <dbl>
## 1  1259           131.  131.   0       
## 2  1260           131.  131.   0       
## 3  1261           131.  131.   0       
## 4  1262           131.  131.   0       
## 5  1263           131.  131.   2.84e-14
## 6  1264           131.  131.   0

max(abs(line_check$difference))    # basically 0 (just rounding error)

## [1] 2.842171e-14

# Visual check --> dashed red line from the first to the last observation, extended
fb_drift_fc |>
  autoplot(fb_stock, level = NULL) +     # level = NULL hides the intervals
  annotate("segment",
           x = 1, y = y_first,
           xend = n_obs + 60, yend = y_first + slope * (n_obs + 60 - 1),
           colour = "red", linetype = "dashed") +
  labs(title = "Drift forecast = line from the first to the last observation",
       subtitle = "Dashed red: line from day 1 to the last day, extended 60 days",
       y = "US$", x = "Trading day")

(d) Other benchmark methods, which one is best?

SNAIVE() wont be used here because trading days don’t have a seasonal period, so I compare MEAN(), NAIVE() and drift. First on the full data then with a train and test splits like the Google example in the book section 5.8

fb_stock |>
  model(Mean  = MEAN(Close),
        Naive = NAIVE(Close),
        Drift = RW(Close ~ drift())) |>
  forecast(h = 60) |>
  autoplot(fb_stock, level = NULL) +
  labs(title = "Facebook: MEAN vs NAIVE vs drift forecasts",
       y = "US$", x = "Trading day", colour = "Method")

fb_train <- fb_stock |> filter(year(Date) <= 2017)   # 2014-2017
fb_test  <- fb_stock |> filter(year(Date) == 2018)   # all of 2018

fb_fit <- fb_train |>
  model(Mean  = MEAN(Close),
        Naive = NAIVE(Close),
        Drift = RW(Close ~ drift()))

# new_data = fb_test so the forecasts line up with the actual 2018 trading days
fb_fc <- fb_fit |> forecast(new_data = fb_test)

fb_fc |>
  autoplot(fb_stock |> filter(year(Date) >= 2017), level = NULL) +
  labs(title = "Forecasts for 2018 using 2014-2017 as training",
       y = "US$", x = "Trading day", colour = "Method")

accuracy(fb_fc, fb_stock) |>
  select(.model, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)

## # A tibble: 3 × 5
##   .model  RMSE   MAE  MAPE  MASE
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 Naive   20.5  16.3  10.2  14.0
## 2 Drift   33.1  24.5  16.0  21.0
## 3 Mean    66.8  63.8  36.3  54.7

# Do the naive residuals look like white noise? (book does this for Google)
fb_stock |> model(NAIVE(Close)) |> gg_tsresiduals()

fb_stock |> model(NAIVE(Close)) |> lb_test(lag = 10)

## # A tibble: 1 × 4
##   Symbol .model       lb_stat lb_pvalue
##   <chr>  <chr>          <dbl>     <dbl>
## 1 FB     NAIVE(Close)    12.1     0.276

Answer. I’d go with naive as the best one as it has the lowest error by a large margin. RMSE is 20.5 compared to drift’s 33.1 and mean 66.8. Mean is sort of useless because the stock trended up since its forecast however naive over drift because drift projected the run up from previous years during the July 2018 crash.

Exercise 5.3 — Seasonal naive on Australian beer production

Apply a seasonal naïve method to the quarterly Australian beer production data from 1992. Check if the residuals look like white noise, and plot the forecasts.

# Extract data of interest
recent_production <- aus_production |>
  filter(year(Quarter) >= 1992)

# Define and estimate a model
fit <- recent_production |> model(SNAIVE(Beer))

# Look at the residuals
fit |> gg_tsresiduals()

# Look at some forecasts
fit |> forecast() |> autoplot(recent_production) +
  labs(title = "Australian beer production: seasonal naive forecasts",
       y = "Megalitres", x = "Quarter")

Ljung-Box test (lag = 8 since it’s quarterly, \(2m = 2 \times 4\)) plus the residual mean to back up what I see in the plots.

lb_test(fit, lag = 8)

## # A tibble: 1 × 3
##   .model       lb_stat lb_pvalue
##   <chr>          <dbl>     <dbl>
## 1 SNAIVE(Beer)    32.3 0.0000834

# the actual autocorrelations as  bounds are about +/- 2/sqrt(T)
augment(fit) |> ACF(.innov, lag_max = 8)

## # A tsibble: 8 x 3 [1Q]
## # Key:       .model [1]
##   .model            lag      acf
##   <chr>        <cf_lag>    <dbl>
## 1 SNAIVE(Beer)       1Q -0.237  
## 2 SNAIVE(Beer)       2Q -0.00987
## 3 SNAIVE(Beer)       3Q  0.234  
## 4 SNAIVE(Beer)       4Q -0.531  
## 5 SNAIVE(Beer)       5Q  0.132  
## 6 SNAIVE(Beer)       6Q -0.0747 
## 7 SNAIVE(Beer)       7Q -0.0927 
## 8 SNAIVE(Beer)       8Q  0.0140

2 / sqrt(nrow(recent_production))

## [1] 0.2324953

# mean of the residuals which should be around 0 if the forecasts are truly unbiased
augment(fit) |>
  as_tibble() |>
  summarise(mean_resid = mean(.innov, na.rm = TRUE))

## # A tibble: 1 × 1
##   mean_resid
##        <dbl>
## 1      -1.57

Answer. So residuals are not white noise. When looking at the ACF at lag 4, its a large negative spike but the Ljung-Box p-value tells us that this is not by chance and that the pattern is indeed real. Since the forecast keeps its the seasonal shape of high 4Q and low 2Q, it sems reasonable, its just not account for the decline as SNAIVE() is prone to do because it’s just repeating the previous year.

Exercise 5.4 — Australian Exports and Bricks

Repeat the previous exercise using the Australian Exports series from global_economy and the Bricks series from aus_production. Use whichever of NAIVE() or SNAIVE() is more appropriate in each case.

(a) Australian Exports (global_economy)

aus_exports <- global_economy |>
  filter(Country == "Australia") |>
  select(Year, Exports)

aus_exports |>
  autoplot(Exports) +
  labs(title = "Australian exports", y = "% of GDP", x = "Year")

# annual data -> no seasonality, so NAIVE
fit_exports <- aus_exports |> model(NAIVE(Exports))

fit_exports |> gg_tsresiduals()

lb_test(fit_exports, lag = 10)                   # non-seasonal -> lag 10

## # A tibble: 1 × 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 NAIVE(Exports)    16.4    0.0896

fit_exports |>
  forecast(h = 10) |>
  autoplot(aus_exports) +
  labs(title = "Australian exports: naive forecasts",
       y = "% of GDP", x = "Year")

Answer. Exports is annual so there’s no seasonality which essentiallyy makes NAIVE() the right pick over SNAIVE(). The Ljung-Box p-value is about 0.09, sol we can’t reject white noise, and the residuals do mostly look random. The forecast is a flat line for the 2017 value (~21% of GDP) with intervals that fan out. It seems fair for the next few years, since exports have mostly bounced between those intervals.

(b) Bricks (aus_production)

# same bricks object as 5.1(b), quarterly -> SNAIVE
fit_bricks <- bricks |> model(SNAIVE(Bricks))

fit_bricks |> gg_tsresiduals()

lb_test(fit_bricks, lag = 8)                     # quarterly -> lag 2m = 8

## # A tibble: 1 × 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 SNAIVE(Bricks)    274.         0

fit_bricks |>
  forecast(h = "3 years") |>
  autoplot(bricks) +
  labs(title = "Bricks: seasonal naive forecasts",
       y = "Millions of bricks", x = "Quarter")

Answer. Bricks are quarterly and seasonal so like before, SNAIVE is the most likely choice. The residuals are not just noise as observed by lag1 at 0.8 then slowly decaying and the Ljung-Box statistic being huge at 274 with a p value of 0 which means that is not by chance.

Exercise 5.7 — Retail series: training/test split with SNAIVE()

For your retail time series (from Exercise 7 in Section 2.10):

set.seed(12345678)                          # same seed as HW 1 and HW 2
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# Same series as before
myseries |> as_tibble() |> distinct(State, Industry, `Series ID`)

## # A tibble: 1 × 3
##   State              Industry                                        `Series ID`
##   <chr>              <chr>                                           <chr>      
## 1 Northern Territory Clothing, footwear and personal accessory reta… A3349767W

(a) Training set before 2011

myseries_train <- myseries |>
  filter(year(Month) < 2011)

# how many months in each piece? (NT series starts in April 1988)
c(train = nrow(myseries_train), test = nrow(myseries) - nrow(myseries_train))

## train  test 
##   273    96

(b) Check the split

autoplot(myseries, Turnover) +
  autolayer(myseries_train, Turnover, colour = "red") +
  labs(title = "NT clothing & footwear turnover: training data in red",
       y = "$ million AUD", x = "Month")

The red covers April 1988 through December 2010 and the black that’s left is 2011 through 2018, so the split worked.

(c) Fit a seasonal naive model

fit <- myseries_train |>
  model(SNAIVE(Turnover))                   # I name the variable to be explicit
fit

## # A mable: 1 x 3
## # Key:     State, Industry [1]
##   State              Industry                                 `SNAIVE(Turnover)`
##   <chr>              <chr>                                               <model>
## 1 Northern Territory Clothing, footwear and personal accesso…           <SNAIVE>

(d) Check the residuals

fit |> gg_tsresiduals()

lb_test(fit, lag = 24)                        # monthly -> lag 2m = 24

## # A tibble: 1 × 5
##   State              Industry                           .model lb_stat lb_pvalue
##   <chr>              <chr>                              <chr>    <dbl>     <dbl>
## 1 Northern Territory Clothing, footwear and personal a… SNAIV…    746.         0

# mean of the residuals
augment(fit) |>
  as_tibble() |>
  summarise(mean_resid = mean(.innov, na.rm = TRUE))

## # A tibble: 1 × 1
##   mean_resid
##        <dbl>
## 1      0.439

Answer. No, the residuals are not uncorrelated the main reason the upward trend while turnover keeps growing, but SNAIVE() just repeats last year, so most months come in higher than what was forecasted and the errors stay positive for longer streches which is why the ACF is high at lag 1 and Ljung-Box rejects white noise . They’re also not the most normal, with a longer left tail that comes from a few big negative residuals. ## (e) Forecasts for the test data

fc <- fit |>
  forecast(new_data = anti_join(myseries, myseries_train))
fc |> autoplot(myseries) +
  labs(title = "SNAIVE forecasts for 2011-2018 vs actual turnover",
       y = "$ million AUD", x = "Month")

(f) Accuracy

fit |> accuracy()             # training set

## # A tibble: 1 × 12
##   State    Industry .model .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>    <chr>    <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Norther… Clothin… SNAIV… Trai… 0.439  1.21 0.915  5.23  12.4     1     1 0.768

fc |> accuracy(myseries)      # test set

## # A tibble: 1 × 12
##   .model    State Industry .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>     <chr> <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(T… Nort… Clothin… Test  0.836  1.55  1.24  5.94  9.06  1.36  1.28 0.601

Answer. The forecast did worst on test compared to training which is expected since test data is data that hasn’t been seen before. RSME which is the typical error went from 1.2 to 1.6 while MASE was 1 exactly on the training but that’s because it measure agaist SNAIVE() itself. The test is 1.36 so error where roughly 36% bigger. MAPE which is the percentage error went down 12.4% to 9.1%. Because sales were larger in 2011 - 2018, an error of the same size is smaller in percent. ME was 0.84 meaning the forecast were too low in average but that’s because SNAIVE() doesn’t take the upward trend into consideration, it just repeats meaning its forecasting lower.

(g) How sensitive are the accuracy measures to the amount of training data?

I tried this two ways.

1. Keep the test set fixed here from 2011–2018 and change where the training starts.

start_years <- c(1988, 1995, 2000, 2005, 2009)
test_2011   <- myseries |> filter(year(Month) >= 2011)

# map_dfr() runs the same fit, forecast and accuracy steps once per start year
# and stacks the rows into one table
sens_start <- purrr::map_dfr(start_years, function(s) {
  train_s <- myseries |> filter(year(Month) >= s, year(Month) < 2011)
  fit_s   <- train_s |> model(SNAIVE(Turnover))
  fc_s    <- fit_s |> forecast(new_data = test_2011)
  bind_rows(
    accuracy(fit_s)            |> mutate(set = "Training"),
    accuracy(fc_s, myseries)   |> mutate(set = "Test")
  ) |>
    mutate(train_start = s, n_train = nrow(train_s))
})

sens_start |>
  select(train_start, n_train, set, RMSE, MAE, MAPE, MASE) |>
  arrange(set, train_start)

## # A tibble: 10 × 7
##    train_start n_train set       RMSE   MAE  MAPE  MASE
##          <dbl>   <int> <chr>    <dbl> <dbl> <dbl> <dbl>
##  1        1988     273 Test     1.55  1.24   9.06  1.36
##  2        1995     192 Test     1.55  1.24   9.06  1.36
##  3        2000     132 Test     1.55  1.24   9.06  1.36
##  4        2005      72 Test     1.55  1.24   9.06  1.36
##  5        2009      24 Test     1.55  1.24   9.06  1.36
##  6        1988     273 Training 1.21  0.915 12.4   1   
##  7        1995     192 Training 1.16  0.902 11.0   1   
##  8        2000     132 Training 0.971 0.778  7.60  1   
##  9        2005      72 Training 1.08  0.895  7.78  1   
## 10        2009      24 Training 1.12  0.817  6.46  1

2. Move the end of the training set which could have more or less history, and forecast the next 2 years after each cutoff so the horizon is the same each time.

cutoffs <- c(2005, 2008, 2011, 2014, 2016)

sens_cutoff <- purrr::map_dfr(cutoffs, function(e) {
  train_e <- myseries |> filter(year(Month) < e)
  test_e  <- myseries |> filter(year(Month) >= e, year(Month) < e + 2)  # 24 mo
  fit_e   <- train_e |> model(SNAIVE(Turnover))
  fc_e    <- fit_e |> forecast(new_data = test_e)
  bind_rows(
    accuracy(fit_e)            |> mutate(set = "Training"),
    accuracy(fc_e, myseries)   |> mutate(set = "Test (next 2 yrs)")
  ) |>
    mutate(train_end = e - 1, n_train = nrow(train_e))
})

sens_cutoff |>
  select(train_end, n_train, set, RMSE, MAE, MAPE, MASE) |>
  arrange(set, train_end)

## # A tibble: 10 × 7
##    train_end n_train set                RMSE   MAE  MAPE  MASE
##        <dbl>   <int> <chr>             <dbl> <dbl> <dbl> <dbl>
##  1      2004     201 Test (next 2 yrs) 0.777 0.608  5.94 0.635
##  2      2007     237 Test (next 2 yrs) 0.987 0.883  7.46 0.940
##  3      2010     273 Test (next 2 yrs) 0.765 0.6    4.59 0.656
##  4      2013     309 Test (next 2 yrs) 0.832 0.708  5.51 0.804
##  5      2015     333 Test (next 2 yrs) 1.52  1.38   9.68 1.61 
##  6      2004     201 Training          1.28  0.959 14.4  1    
##  7      2007     237 Training          1.25  0.940 13.4  1    
##  8      2010     273 Training          1.21  0.915 12.4  1    
##  9      2013     309 Training          1.18  0.881 11.5  1    
## 10      2015     333 Training          1.15  0.860 11.0  1

sens_cutoff |>
  ggplot(aes(x = train_end, y = RMSE, colour = set)) +
  geom_line() +
  geom_point() +
  labs(title = "RMSE vs. where the training set ends",
       x = "Last year of training data", y = "RMSE", colour = NULL)

Answer. With SNAIVE(), there is no difference in test accuracy based upon the number of years used as training data. Since SNAIVE() only uses the past 12 months before the cut off date, so beginning the training period in either 1988 or 2009 gives the same predictions and the same test RMSE of 1.55. However there is a large change in the training errors when I cut off the noisy early 1990’s from the in-sample. Moving the cut off point to another year changed the test accuracy by a large amount looking at RMSE move from approximately 0.78 to 1.52. This was not due to the amount of data used as training data but due to the fact that a different set of years were being forecasted each time, and some years are more difficult to forecast than others. For instance sales increased significantly in 2016 and last year’s trend did not anticipate this increase, therefore the results obtained using SNAIVE() depend heavily upon where you divide your data so this is why it is suggested to use time series cross validation in the book section 5.10 rather than relying on one train/test division.

Session info

# Package versions for reproducibility
sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Sequoia 15.7.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## 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     
##  [5] tsibbledata_0.4.1 tsibble_1.2.0     ggplot2_4.0.3     lubridate_1.9.5  
##  [9] tidyr_1.3.2       dplyr_1.2.1       tibble_3.3.1      fpp3_1.0.3       
## 
## loaded via a namespace (and not attached):
##  [1] ggdist_3.3.3         utf8_1.2.6           rappdirs_0.3.4      
##  [4] sass_0.4.10          generics_0.1.4       anytime_0.3.13      
##  [7] digest_0.6.39        magrittr_2.0.5       evaluate_1.0.5      
## [10] grid_4.6.1           timechange_0.4.0     RColorBrewer_1.1-3  
## [13] mixtime_0.3.0        fastmap_1.2.0        jsonlite_2.0.0      
## [16] purrr_1.2.2          scales_1.4.0         jquerylib_0.1.4     
## [19] cli_3.6.6            rlang_1.3.0          crayon_1.5.3        
## [22] vecvec_1.3.0         withr_3.0.3          cachem_1.1.0        
## [25] yaml_2.3.12          otel_0.2.0           tools_4.6.1         
## [28] tzdb_0.5.0           vctrs_0.7.3          R6_2.6.1            
## [31] lifecycle_1.0.5      pkgconfig_2.0.3      progressr_1.0.0     
## [34] pillar_1.11.1        bslib_0.12.0         gtable_0.3.6        
## [37] glue_1.8.1           Rcpp_1.1.2           xfun_0.60           
## [40] tidyselect_1.2.1     rstudioapi_0.19.0    knitr_1.51          
## [43] farver_2.1.2         htmltools_0.5.9      labeling_0.4.3      
## [46] rmarkdown_2.32       compiler_4.6.1       S7_0.2.2            
## [49] distributional_0.8.1