Rpubs URL: https://rpubs.com/umaisabdullah/DATA624_HW3

Objective

Complete and submit exercises 5.1, 5.2, 5.3, 5.4 and 5.7 for Homework 3 from the Hyndman online Forecasting book. We have to submit both our Rpubs link as well as attach the .pdf file with our code.

Overview

Homework 1 was about looking at a time series and naming the patterns in it. Homework 2 was about preparing a series for modelling by transforming it and splitting it into parts. This homework is the first one where I actually forecast something.

The methods here are deliberately simple. The naive method says the next value equals the last value. The seasonal naive method says the next value equals the value from the same period last year. The drift method draws a straight line from the first observation to the last one and continues it. The mean method says every future value equals the average of the past.

These are called benchmark methods, and the point of them is not that they are good. The point is that they give me something to beat. If a complicated model cannot outperform a method that just repeats last year’s number, then the complicated model is not earning its keep. So the real skill in this chapter is not fitting the models, which takes one line each. It is checking whether a model is any good, which means looking at the residuals and measuring accuracy on data the model never saw.

Two ideas run through every exercise below.

The first is residual diagnostics. A residual is the gap between what actually happened and what the model said would happen. If a model has captured everything useful in the data, the residuals should look like random noise with no pattern left in them. Any pattern still sitting in the residuals is information the model failed to use.

The second is the difference between training accuracy and test accuracy. A model always looks better on the data it was fitted to than on new data. So a number that comes from the training set tells me how well the model memorised the past, not how well it will forecast the future. Only a test set held back from fitting can tell me that.

Pre-Requisite

All of the data for this homework comes from packages, so nothing needs to be downloaded into the working directory. Only one package is required:

  1. fpp3 loads the core forecasting functions and all of the datasets used below.

Exercise 5.1

Question. Produce forecasts for the following series using whichever of NAIVE(y), SNAIVE(y) or RW(y ~ drift()) is more appropriate in each case:

My thought process. The exercise gives me three methods and five series, and the whole task is matching them up. I do not want to guess, so I use two questions to decide each one.

The first question is whether the series has a seasonal pattern. If it does, seasonal naive is the only one of the three that can reproduce it, because the other two carry forward a single number and a flat line cannot have a season in it.

The second question is whether the series has a steady trend. If it does and there is no season, drift is right, because drift continues the slope while naive flattens out. If there is no trend and no season, naive is the honest choice.

So I plot each series before choosing, and I say why the method fits rather than just naming it.

Australian Population

aus_pop <- global_economy |>
  filter(Country == "Australia") |>
  select(Year, Population)

aus_pop |>
  autoplot(Population) +
  labs(title = "Australian population",
       x = "Year", y = "Population")

This is close to a straight upward line with no season, because population data is annual so there is no within year pattern to have. Drift is the right choice. It continues the climb, while naive would predict the population stops growing the moment the data ends, which is clearly wrong.

aus_pop |>
  model(Drift = RW(Population ~ drift())) |>
  forecast(h = 10) |>
  autoplot(aus_pop) +
  labs(title = "Australian population, 10 year drift forecast",
       x = "Year", y = "Population")

Method chosen: drift. The forecast continues the upward slope at roughly 251,000 people a year, which is the average yearly growth across the whole history. The forecast for 2018 is about 24.85 million, rising to about 25.35 million by 2020.

Bricks

bricks <- aus_production |>
  filter(!is.na(Bricks)) |>
  select(Bricks)

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

Two things matter here. The data is quarterly and there is a visible pattern that repeats each year, so the series is seasonal. There is also no single consistent trend. Production rises until the early 1980s and then drifts down and sideways for the rest of the series.

That combination rules out drift. Drift would draw a line from roughly 190 in 1956 to 435 in 2005 and forecast continued growth, even though production has been flat to declining for 25 years. It would also ignore the seasonal swing entirely. Seasonal naive handles both problems by carrying forward the last four quarters.

bricks |>
  model(SNAIVE(Bricks)) |>
  forecast(h = "3 years") |>
  autoplot(bricks) +
  labs(title = "Brick production, 3 year seasonal naive forecast",
       x = "Quarter", y = "Bricks (millions)")

Method chosen: seasonal naive. It repeats the last full year, which was 428, 397, 355 and 435 for Q3 2004 through Q2 2005, and keeps repeating it. The forecast is flat from year to year but keeps the quarterly shape.

Notice the series ends in 2005 Q2 rather than Q4. That matters because seasonal naive always looks back exactly four quarters, so the forecast for 2005 Q3 comes from 2004 Q3, not from the most recent observation.

NSW Lambs

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

lambs |>
  autoplot(Count) +
  labs(title = "Lambs slaughtered in New South Wales",
       x = "Month", y = "Count")

This is monthly data with a repeating within year pattern, so it is seasonal. The level moves around over the decades but there is no straight trend to extend, so seasonal naive fits better than drift.

lambs |>
  model(SNAIVE(Count)) |>
  forecast(h = "2 years") |>
  autoplot(lambs) +
  labs(title = "NSW lambs, 2 year seasonal naive forecast",
       x = "Month", y = "Count")

Method chosen: seasonal naive. The forecast repeats the last twelve months of data. The prediction intervals here are very wide, which is worth noticing rather than passing over. This series is noisy from month to month, and seasonal naive builds its intervals from how badly the same month differed from a year earlier in the past. Noisy history means wide intervals, and that is the model being honest about how uncertain it is.

Household wealth

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

This is annual data for four countries, so there is no seasonality to capture. All four countries trend upward over the period despite a dip around the 2008 financial crisis. Drift handles that because it links the first and last observations for each country separately.

hh_budget |>
  model(Drift = RW(Wealth ~ drift())) |>
  forecast(h = 5) |>
  autoplot(hh_budget) +
  labs(title = "Household wealth, 5 year drift forecast",
       x = "Year", y = "Wealth (% of income)")

Method chosen: drift. Each of the four countries gets its own slope, which is the point of using drift on a keyed dataset rather than fitting one line to everything. All four forecasts continue upward from their 2016 values of roughly 422 for Australia, 565 for Canada, 602 for Japan and 609 for the USA.

Australian takeaway food turnover

takeaway <- aus_retail |>
  filter(Industry == "Takeaway food services") |>
  summarise(Turnover = sum(Turnover))

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

This one has both a strong upward trend and a clear monthly seasonal pattern, which is the awkward case. None of the three methods handles both at once.

I choose seasonal naive because the seasonal pattern is the part that a flat forecast would get badly wrong every single month, while the trend is gradual enough that ignoring it costs less over a short horizon. Drift would capture the growth but would produce a smooth line with no December peak at all, which is a worse mistake for a food retail series.

takeaway |>
  model(SNAIVE(Turnover)) |>
  forecast(h = "2 years") |>
  autoplot(takeaway) +
  labs(title = "Takeaway food turnover, 2 year seasonal naive forecast",
       x = "Month", y = "Turnover ($ million)")

Method chosen: seasonal naive. The forecast repeats 2018 for both 2019 and 2020.

The weakness is visible in the plot. The forecast is flat from one year to the next while the history has been climbing steadily, so seasonal naive will systematically forecast too low. That is a real limitation and it is exactly the kind of gap a proper model is meant to close later. For now the value of this forecast is as a benchmark to beat, not as a forecast I would actually use.

Exercise 5.2

Question. Use the Facebook stock price (data set gafa_stock) to do the following:

My thought process. There is a practical problem to solve before any forecasting happens. Stock data only exists on trading days, so the series skips weekends and holidays. A tsibble indexed by calendar date therefore has gaps in it, and the forecasting functions will complain. The standard fix is to re-index by trading day number so the series becomes regular.

Part c is the interesting one. It asks me to prove a claim rather than just run a function, so I check the arithmetic by hand and compare it against what the model returns.

(a) Time plot

fb_stock <- gafa_stock |>
  filter(Symbol == "FB") |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)

fb_stock |>
  autoplot(Close) +
  labs(title = "Facebook closing stock price, 2014 to 2018",
       x = "Trading day", y = "Closing price (USD)")

The price climbs from around 55 dollars at the start of 2014 to a peak above 215 in mid 2018, then falls sharply in the second half of 2018 and ends the year around 131. That late drop matters for the rest of this exercise.

(b) Drift forecasts

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

fb_fc <- fb_fit |>
  forecast(h = 63)

fb_fc |>
  autoplot(fb_stock) +
  labs(title = "Facebook closing price, drift forecast for 63 trading days",
       x = "Trading day", y = "Closing price (USD)")

(c) Showing the forecast is the line from first to last observation

The claim. The drift method estimates its slope as the total change divided by the number of steps, so the slope is the last observation minus the first, divided by one less than the number of observations. That is exactly the slope of a straight line joining the first point to the last point. If that is true, I should be able to reproduce the model’s forecast with arithmetic.

first_close <- fb_stock$Close[1]
last_close  <- fb_stock$Close[nrow(fb_stock)]
n_obs       <- nrow(fb_stock)

slope <- (last_close - first_close) / (n_obs - 1)

c(first = first_close, last = last_close, n = n_obs, slope = slope)
##        first         last            n        slope 
## 5.471000e+01 1.310900e+02 1.258000e+03 6.076372e-02
# Forecast by hand for each horizon, then compare with what the model produced.
manual <- tibble(
  h = 1:63,
  by_hand = last_close + (1:63) * slope
)

comparison <- fb_fc |>
  as_tibble() |>
  mutate(h = row_number()) |>
  select(h, from_model = .mean) |>
  left_join(manual, by = "h") |>
  mutate(difference = from_model - by_hand)

head(comparison, 3)
## # A tibble: 3 × 4
##       h from_model by_hand difference
##   <int>      <dbl>   <dbl>      <dbl>
## 1     1       131.    131.          0
## 2     2       131.    131.          0
## 3     3       131.    131.          0
tail(comparison, 3)
## # A tibble: 3 × 4
##       h from_model by_hand difference
##   <int>      <dbl>   <dbl>      <dbl>
## 1    61       135.    135.   0       
## 2    62       135.    135.   0       
## 3    63       135.    135.  -2.84e-14
# The largest gap anywhere in the 63 forecasts.
max(abs(comparison$difference))
## [1] 2.842171e-14

The proof. The first close is 54.71, the last is 131.09, and there are 1258 observations. The slope works out to 0.0608 dollars per trading day. Forecasting by hand gives 131.15 at horizon 1 and 134.92 at horizon 63, and the model returns the same values. The largest difference anywhere across all 63 forecasts is zero apart from floating point rounding.

Drawing the line makes the same point visually.

fb_fc |>
  autoplot(fb_stock) +
  geom_segment(aes(x = 1, y = first_close,
                   xend = n_obs + 63,
                   yend = last_close + 63 * slope),
               colour = "red", linetype = "dashed") +
  labs(title = "Drift forecast lies exactly on the line from first to last observation",
       subtitle = "Red dashed line joins day 1 to the final day and continues",
       x = "Trading day", y = "Closing price (USD)")

The forecast sits exactly on the dashed line. This is worth understanding rather than just verifying, because it tells me what drift actually pays attention to. Drift uses two numbers, the first observation and the last one. Everything in between, including the entire 2018 peak and crash, has no effect on the slope at all.

(d) Comparing the other benchmark methods

A plot alone will not settle which method is best, so I split the data and measure. The training set is everything before 2018 and the test set is all of 2018, which the models never see.

fb_train <- fb_stock |> filter(Date < as.Date("2018-01-01"))
fb_test  <- fb_stock |> filter(Date >= as.Date("2018-01-01"))

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

fb_bench_fc <- fb_models |>
  forecast(h = nrow(fb_test))

fb_bench_fc |>
  autoplot(fb_train, level = NULL) +
  autolayer(fb_test, Close, colour = "black") +
  labs(title = "Three benchmark methods against the actual 2018 prices",
       subtitle = "Black line shows what actually happened",
       x = "Trading day", y = "Closing price (USD)")

accuracy(fb_bench_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

Which is best, and why. Naive wins on every measure, with an RMSE of 20.5 against 33.1 for drift and 66.8 for mean. The ranking is naive first, drift second, mean a distant third.

The mean method losing is no surprise. It forecasts the average of the last four years, which is far below where the price actually was at the end of 2017, so it is wrong from the first day.

Drift losing to naive is the result worth explaining. Drift fits a line from the 2014 price of 55 to the end of 2017 price of about 176, which is a steep upward slope, so it forecasts the price continuing to rise through 2018. The price actually fell. Naive predicts no change at all, and in a year where the price ended lower than it started, predicting no change turned out to be closer than predicting a rise.

There is a general lesson underneath this specific result. Stock prices behave much like a random walk, where the best available estimate of tomorrow’s price genuinely is today’s price. Drift only beats naive when a trend keeps going, and assuming a past trend continues is precisely the assumption that fails for stock prices. The methods with the fewest assumptions tend to do best on this kind of data.

Exercise 5.3

Question. Apply a seasonal naive method to the quarterly Australian beer production data from 1992. Check if the residuals look like white noise, and plot the forecasts. What do you conclude?

My thought process. The exercise hands me the code, so the work is in reading the output rather than writing it. Checking for white noise means asking three separate questions, and I want to answer each one deliberately rather than glancing at the plots and calling it fine.

The questions are whether the residuals are uncorrelated, whether they have a mean near zero, and whether their spread stays steady. I also run a Ljung-Box test, because the autocorrelation plot shows me spikes but does not tell me whether they are large enough to take seriously.

recent_production <- aus_production |>
  filter(year(Quarter) >= 1992)

fit <- recent_production |>
  model(SNAIVE(Beer))
fit |> gg_tsresiduals()

augment(fit) |>
  features(.innov, ljung_box, lag = 8)
## # A tibble: 1 × 3
##   .model       lb_stat lb_pvalue
##   <chr>          <dbl>     <dbl>
## 1 SNAIVE(Beer)    32.3 0.0000834
resid_beer <- augment(fit)$.innov

c(mean = mean(resid_beer, na.rm = TRUE),
  sd   = sd(resid_beer, na.rm = TRUE))
##      mean        sd 
## -1.571429 16.174386
fit |>
  forecast() |>
  autoplot(recent_production) +
  labs(title = "Australian beer production, seasonal naive forecast",
       x = "Quarter", y = "Beer (megalitres)")

Are the residuals uncorrelated? No. The autocorrelation plot shows a clear spike at lag 4 that reaches outside the significance bounds. The Ljung-Box test gives a statistic of 32.3 with a p-value of 0.00008, so I reject the idea that the residuals are independent.

A spike at lag 4 specifically is informative rather than just a failure. Lag 4 on quarterly data means one year, so what is left in the residuals is a year to year relationship. Seasonal naive uses only last year’s value for the same quarter and nothing else, so it cannot pick up the fact that the last several years together carry information about the next one. That leftover structure shows up exactly where the plot puts it.

Is the mean near zero? Close but not quite. The residual mean is -1.57 megalitres against a residual standard deviation of 16.2, so the bias is small relative to the noise. The negative sign means the forecasts run slightly high on average, which makes sense because beer production drifts gently downward over this period and seasonal naive always forecasts last year’s level.

Is the variance steady? Yes. The time plot of residuals shows a band of roughly constant width across the whole period with no fanning out at either end.

Is the distribution normal? Roughly. The histogram is centred near zero and has a single peak with no severe skew. It is not a textbook bell curve, but it is close enough that prediction intervals built on a normal assumption will not be badly wrong.

What I conclude. The residuals are not white noise, because the correlation at lag 4 is real and statistically significant. So the model has not extracted everything useful from this series.

That said, the failure is a mild one and the model is still usable. The forecasts track the seasonal shape well, the MAPE on the training data is 3.16 percent, and the bias is small. What the diagnostics tell me is that a better model exists, not that this one is broken. Since the residuals are not white noise, the prediction intervals are probably a little narrower than they should be, because the interval calculation assumes the errors are independent and they are not.

For a benchmark this is a good result. Seasonal naive sets a bar at about 3 percent error, and any more sophisticated model I fit later has to clear that bar to justify itself.

Exercise 5.4

Question. 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.

My thought process. The method choice is the first decision and it follows the same rule as Exercise 5.1. Exports is annual data, so there is no season for seasonal naive to use and naive is the only sensible choice of the two. Bricks is quarterly with a clear yearly pattern, so seasonal naive it is.

What makes this exercise worth doing is that the two series give opposite answers on the residual check, so I get to see what passing and failing each look like.

Australian Exports

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

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

Annual data, so there is no season for seasonal naive to use, and naive is the only sensible choice of the two. The series does trend upward, which naive will not capture.

fit_exports <- aus_exports |>
  model(NAIVE(Exports))

fit_exports |> gg_tsresiduals()

augment(fit_exports) |>
  features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 NAIVE(Exports)    16.4    0.0896
accuracy(fit_exports) |>
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 1 × 4
##   .model          RMSE   MAE  MAPE
##   <chr>          <dbl> <dbl> <dbl>
## 1 NAIVE(Exports)  1.24 0.985  5.83
fit_exports |>
  forecast(h = 10) |>
  autoplot(aus_exports) +
  labs(title = "Australian exports, 10 year naive forecast",
       x = "Year", y = "Exports (% of GDP)")

Residual check. This one passes, though only just. The Ljung-Box test gives a statistic of 16.4 with a p-value of 0.090, which is above 0.05, so there is no significant evidence against independence overall. The autocorrelation plot is mostly consistent with that: lag 1 sits just outside the bounds at about -0.3, but every other spike is inside them. The residual mean is 0.15, which is small against a series that sits between about 12 and 23, so the bias is minor. It is not zero by accident, though. It equals the average yearly increase in the series, (21.3 - 13.0) / 57 about 0.15, so it is the upward trend that naive ignores. The histogram is roughly symmetric around zero.

Conclusion for exports. The residuals are close to white noise, which means naive has captured what there is to capture in this series. That is not because naive is clever. It is because the series behaves like a random walk, where each year’s value is last year’s value plus an unpredictable shock. When that is genuinely how a series works, naive is not just a benchmark, it is close to the correct model.

The forecast is a flat line at the last observed value of 21.3 percent, with intervals that widen as the horizon grows. The widening is right, since uncertainty about a random walk really does accumulate the further out you look.

Bricks

bricks_series <- aus_production |>
  filter(!is.na(Bricks)) |>
  select(Bricks)

fit_bricks <- bricks_series |>
  model(SNAIVE(Bricks))
fit_bricks |> gg_tsresiduals()

augment(fit_bricks) |>
  features(.innov, ljung_box, lag = 8)
## # A tibble: 1 × 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 SNAIVE(Bricks)    274.         0
accuracy(fit_bricks) |>
  select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 1 × 5
##   .model          RMSE   MAE  MAPE  MASE
##   <chr>          <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(Bricks)  48.3  35.5  8.84     1
fit_bricks |>
  forecast(h = "3 years") |>
  autoplot(bricks_series) +
  labs(title = "Brick production, 3 year seasonal naive forecast",
       x = "Quarter", y = "Bricks (millions)")

Residual check. This one fails, and badly. The Ljung-Box statistic is 274 with a p-value indistinguishable from zero. The autocorrelation plot shows large spikes at the early lags that decay slowly rather than dropping away, which is the signature of a series where consecutive residuals are strongly related.

The residual time plot adds a second problem that the beer series did not have. The spread is not constant. Residuals are small in the 1950s and 1960s when production was low, much larger through the 1970s and 1980s when production was high, and the swings around the 1982 downturn are the largest anywhere in the plot. The residual mean is 4.2, which is a mild upward bias.

Conclusion for bricks. The residuals are clearly not white noise and this is a much worse result than either of the other two series. There are two separate reasons.

The first is that brick production goes through long cycles tied to the construction industry, which rise and fall over several years. Seasonal naive only ever looks back exactly four quarters, so it has no way to know the series is in the middle of a multi year downturn. When production falls for three years running, the model forecasts last year’s higher numbers every single time and is wrong in the same direction throughout. That sustained one directional error is what produces the slowly decaying autocorrelation.

The second is the changing variance. Seasonal naive produces errors roughly proportional to the level of the series, so as production grew the errors grew with it. A Box-Cox transformation of the kind used in Homework 2 would help with that part.

Comparing the three series. Putting the exercises side by side makes the point clearly. Exports passed the residual test, beer failed mildly, and bricks failed badly. The difference is not the method. It is how well each series matches what the method assumes. Naive suits exports because exports behave like a random walk. Seasonal naive suits beer reasonably because beer has a stable season and little trend. Seasonal naive suits bricks poorly because bricks has long cycles and growing variance, and neither of those is something the method can represent.

The residual diagnostics found all of this without me needing to know anything about the construction industry in advance, which is the reason to run them.

Exercise 5.7

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

My thought process. This exercise continues the retail series from Homework 1 and Homework 2, so I use the same seed to draw the same series. Keeping it consistent matters, because the answers would not line up with my earlier homework otherwise.

Part g is the one I want to spend real effort on. The obvious reading is that more training data gives better forecasts, and I could write that sentence without running anything. Instead I am going to actually vary the training length and look at what comes out, because with seasonal naive specifically I suspect the obvious answer is wrong.

(a) Create the training set

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# Confirm this is the same series used in Homework 1 and Homework 2.
myseries |>
  as_tibble() |>
  distinct(State, Industry, `Series ID`)
## # A tibble: 1 × 3
##   State           Industry               `Series ID`
##   <chr>           <chr>                  <chr>      
## 1 New South Wales Takeaway food services A3349792X
myseries_train <- myseries |>
  filter(year(Month) < 2011)

c(full = nrow(myseries),
  train = nrow(myseries_train),
  test = nrow(myseries) - nrow(myseries_train))
##  full train  test 
##   441   345    96

The series is New South Wales takeaway food services, the same one as in the previous two homeworks. The full series has 441 monthly observations running from April 1982 to December 2018. The split gives 345 months of training data and 96 months, eight years, of test data.

(b) Check the split

autoplot(myseries, Turnover) +
  autolayer(myseries_train, Turnover, colour = "red") +
  labs(title = "Retail turnover, full series in black and training portion in red",
       x = "Month", y = "Turnover ($ million)")

The red section covers everything up to the end of 2010 and the black section continues from 2011 to the end. The split is where it should be.

The plot also shows something that matters for what comes later. The training period ends at a turnover level of around 400 to 460 million, while the test period climbs well beyond that, reaching over 700 by 2018. The model is being asked to forecast a period where the series goes considerably higher than anything it was trained on.

(c) Fit the seasonal naive model

fit_retail <- myseries_train |>
  model(SNAIVE(Turnover))

fit_retail
## # A mable: 1 x 3
## # Key:     State, Industry [1]
##   State           Industry               `SNAIVE(Turnover)`
##   <chr>           <chr>                             <model>
## 1 New South Wales Takeaway food services           <SNAIVE>

(d) Check the residuals

fit_retail |> gg_tsresiduals()

augment(fit_retail) |>
  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 Takeaway food services SNAIVE(Turnover)   1101.         0
resid_retail <- augment(fit_retail)$.innov

c(mean = mean(resid_retail, na.rm = TRUE),
  sd = sd(resid_retail, na.rm = TRUE))
##     mean       sd 
## 11.49940 23.48467

Are the residuals uncorrelated? No, and this is the clearest failure of the three residual checks in this homework. The Ljung-Box statistic at lag 24 is 1101 with a p-value of zero. The autocorrelation plot shows nearly every lag sitting well outside the significance bounds, starting at about 0.89 at lag 1 and decaying slowly rather than dropping off.

Autocorrelation that high at lag 1 means each residual is strongly predictable from the one before it. The cause is the trend. This series grows almost every year, so seasonal naive, which forecasts last year’s value, is too low nearly every month. A long run of errors that are all positive is by definition strongly autocorrelated.

Are they normally distributed? Roughly in shape but not in centre. The histogram has a single peak and is not far from symmetric, with a mild right skew of about 0.36. The real problem is location, not shape. The residual mean is 11.5 rather than zero, against a standard deviation of 23.5. That is a bias of about half a standard deviation, which is substantial.

That number quantifies exactly what I said above. On average the model forecasts 11.5 million dollars too low every month, because the series grows by roughly that much year over year and the model assumes no growth at all.

Answering the question directly. No on both counts. The residuals are strongly correlated and they are not centred on zero. The practical consequence is that the prediction intervals from this model will be too narrow and positioned too low, so they will cover the true value less often than the stated 80 or 95 percent.

(e) Forecast the test data

fc_retail <- fit_retail |>
  forecast(new_data = anti_join(myseries, myseries_train,
                                by = c("State", "Industry", "Series ID",
                                       "Month", "Turnover")))

fc_retail |>
  autoplot(myseries) +
  labs(title = "Seasonal naive forecast against the actual test period",
       x = "Month", y = "Turnover ($ million)")

The shape of the failure is easy to see. The forecast is a flat repetition of 2010 stretched across eight years, while the actual series keeps climbing away from it. The gap widens the whole way, and by 2018 the actual values are far above even the upper 95 percent interval.

(f) Compare accuracy

# Accuracy on the data the model was fitted to.
accuracy(fit_retail) |>
  select(.model, RMSE, MAE, MAPE, MASE, ACF1)
## # A tibble: 1 × 6
##   .model            RMSE   MAE  MAPE  MASE  ACF1
##   <chr>            <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(Turnover)  26.1  19.2  9.59     1 0.890
# Accuracy on the test data the model never saw.
accuracy(fc_retail, myseries) |>
  select(.model, RMSE, MAE, MAPE, MASE, ACF1)
## # A tibble: 1 × 6
##   .model            RMSE   MAE  MAPE  MASE  ACF1
##   <chr>            <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(Turnover)  96.8  79.5  16.3  4.14 0.964

What the comparison shows. Every measure gets substantially worse on the test set.

The RMSE rises from 26.1 on the training data to 96.8 on the test data, which is nearly four times larger. The MAPE goes from 9.6 percent to 16.3 percent. The MASE goes from 1 to 4.14, and that one is worth reading carefully. MASE is scaled so that a value of 1 means the model matches a naive forecast on the training data, which is by definition true for the model the scaling is based on. A test MASE of 4.14 means the test errors are more than four times the size of the typical training error.

Why the gap is this large. The usual reason training accuracy beats test accuracy is that a model fits noise in the data it was trained on. That is not what is happening here, because seasonal naive has no parameters to overfit with.

The real reason is the eight year horizon combined with the trend. On the training data the model is only ever forecasting one year ahead, so it is wrong by about one year of growth. On the test data it is forecasting up to eight years ahead while still repeating 2010, so by the end it is wrong by eight years of accumulated growth. The error grows with the horizon because the thing the model is ignoring, the trend, keeps accumulating.

The test ACF1 of 0.96 confirms the errors are not random. They are almost perfectly correlated, which is what happens when the forecast is too low every single month by a steadily growing amount.

(g) Sensitivity to the amount of training data

Setting up the question properly. To test this fairly I have to change one thing at a time. I keep the test set fixed at 2011 onward and vary only where the training data starts, so any change in accuracy comes from the training length and nothing else.

test_fixed <- myseries |> filter(year(Month) >= 2011)

sensitivity <- purrr::map_dfr(c(1995, 2000, 2005, 2008, 2010), function(start_yr) {
  tr <- myseries |>
    filter(year(Month) >= start_yr, year(Month) < 2011)

  fc <- tr |>
    model(SNAIVE(Turnover)) |>
    forecast(new_data = test_fixed)

  acc <- accuracy(fc, myseries)

  tibble(start_year = start_yr,
         train_months = nrow(tr),
         RMSE = round(acc$RMSE, 2),
         MAE = round(acc$MAE, 2),
         MAPE = round(acc$MAPE, 2))
})

sensitivity
## # A tibble: 5 × 5
##   start_year train_months  RMSE   MAE  MAPE
##        <dbl>        <int> <dbl> <dbl> <dbl>
## 1       1995          192  96.8  79.6  16.3
## 2       2000          132  96.8  79.6  16.3
## 3       2005           72  96.8  79.6  16.3
## 4       2008           36  96.8  79.6  16.3
## 5       2010           12  96.8  79.6  16.3

The result is not what I expected. The accuracy does not change at all. Training on 192 months gives exactly the same RMSE of 96.8 as training on 12 months. Cutting the training data to one sixteenth of its size costs nothing.

Why this happens. Once I saw it, the reason is obvious in hindsight. Seasonal naive forecasts each month by copying the same month from the previous year. For a test set starting in 2011, every forecast traces back to a value in 2010. Nothing before 2010 is ever used.

So the extra years of training data are not being ignored because they are unhelpful. They are ignored because the method never looks at them. Training on 16 years and training on 1 year produce an identical model, because both end with the same twelve values.

# These twelve numbers are the entire model, regardless of training length.
myseries |>
  filter(year(Month) == 2010) |>
  as_tibble() |>
  select(Month, Turnover)
## # A tibble: 12 × 2
##       Month Turnover
##       <mth>    <dbl>
##  1 2010 Jan     413.
##  2 2010 Feb     346 
##  3 2010 Mar     368.
##  4 2010 Apr     395.
##  5 2010 May     389.
##  6 2010 Jun     393.
##  7 2010 Jul     431.
##  8 2010 Aug     411.
##  9 2010 Sep     421.
## 10 2010 Oct     425.
## 11 2010 Nov     403 
## 12 2010 Dec     464.

Checking whether this is a general result or a quirk of the method. A finding like this is only useful if I know how far it extends, so I ran the same test on drift and mean, which do use the whole training set.

sens_others <- purrr::map_dfr(c(1995, 2000, 2005, 2008, 2010), function(start_yr) {
  tr <- myseries |>
    filter(year(Month) >= start_yr, year(Month) < 2011)

  fc <- tr |>
    model(Drift = RW(Turnover ~ drift()),
          Mean  = MEAN(Turnover)) |>
    forecast(new_data = test_fixed)

  acc <- accuracy(fc, myseries)

  tibble(start_year = start_yr,
         train_months = nrow(tr),
         model = acc$.model,
         RMSE = round(acc$RMSE, 1))
})

sens_others |>
  tidyr::pivot_wider(names_from = model, values_from = RMSE)
## # A tibble: 5 × 4
##   start_year train_months Drift  Mean
##        <dbl>        <int> <dbl> <dbl>
## 1       1995          192  98.5  227.
## 2       2000          132 129.   204.
## 3       2005           72 152.   164.
## 4       2008           36 253.   135.
## 5       2010           12 242.   100.

These two are very sensitive to the training length, and interestingly they move in opposite directions.

Drift gets worse with less data, from an RMSE of 98.5 using 192 months to 242 using 12 months. That makes sense, because drift estimates its slope from the first and last observations only. With a short window those two points are close together and the slope estimate becomes unstable.

Mean gets better with less data, from 227 down to 100. That is not because short windows suit the mean method. It is because the series has a strong trend, so a mean taken over 16 years sits far below the 2010 level, while a mean taken over 2010 alone is close to it. The short window happens to land nearer the truth by accident.

Answering the question. The honest answer is that it depends entirely on the method, and for seasonal naive on this series the accuracy is completely insensitive to the amount of training data.

That is the useful conclusion here. The instinct that more training data always helps is wrong, because it assumes the model actually uses the extra data. Seasonal naive uses only the final season no matter how much history it is given.

What genuinely does move the accuracy for this series is the forecast horizon rather than the training size. Forecasting two years ahead from 2013 gives an RMSE of about 39, while forecasting eight years ahead from 2011 gives 96.8. When a model ignores the trend, the error grows with how far ahead it is asked to look, and the length of its history has nothing to do with it.