Chapter 5, Hyndman and Athanasopoulos

  1. 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)

Looking at the data there is clearly a linear trend, so the drift model would show the continuation of the trend that is inherent in the data.

global_economy %>%
  filter(Country == "Australia") %>%
  autoplot(Population) +
  labs(title = "Australian Population")

The drift model forecast is constructed based on a random walk (RW) along the drift trend inherent in the data. The plot belows shows the forecast out 10 years based on the drift along 80% and 95% cofidence intervals.

global_economy %>%
  filter(Country == "Australia") %>%
  model(RW(Population ~ drift())) %>%
  forecast(h = 10) %>%
  autoplot(global_economy) +
  labs(title = "10-year Forecast of Australian Population")

Bricks (aus_production)

Looking at the data there is clear seasonality along with some trend, so the seasonal naive forecast would be best in this case.

aus_production %>%
  filter(!is.na(Bricks)) %>%
  autoplot(Bricks) +
  labs(title = "Autralian Bricks")

The seasonal naive forecast, show in blue below, demonstrates the seasonality but without any trend. Again the darker blue is within the 80% confidence interval and the lighter blue is showing the 95% confidence interval in the forecast.

aus_production %>%
  filter(!is.na(Bricks)) %>%
  model(SNAIVE(Bricks ~ lag("year"))) %>%
  forecast(h = 20) %>%
  autoplot(aus_production) +
  labs(title = "5-year forecast of Autralian Bricks")

NSW Lambs (aus_livestock)

Looking at the historical data for lamb slaughter in Australia, there is no clear trend or seasonality.

aus_livestock %>%
  filter(State == "New South Wales", Animal == "Lambs") %>%
  autoplot(Count) +
  labs(title = "NSW Lambs")

Checking the residuals of a NAIVE STL decomposition shows no clear trend or seasonality, just white noise but with some disorganized significance in the acf lags.

lamb_dcmp <- aus_livestock %>%
  filter(State == "New South Wales", Animal == "Lambs") %>%
  model(stlf = decomposition_model(
    STL(Count ~ trend(window = 60), robust = TRUE),
    NAIVE(season_adjust)
  ))

lamb_dcmp %>%
  gg_tsresiduals()

Therefore, the Naive model should be sufficient to forecast at this time. But the STL decomposition can also “reseasonalise” the forecast data, which appears to be a little more refined.

aus_livestock %>%
  filter(State == "New South Wales", Animal == "Lambs") %>%
  model(NAIVE(Count)) %>%
  forecast(h = 60) %>%
  autoplot(aus_livestock) +
  labs(y = "Lambs", title = "Monthly NSW Lambs slaughtered Naive method")

lamb_dcmp %>%
  forecast() %>%
  autoplot(aus_livestock) +
  labs(y = "Lambs", title = "Monthly NSW Lambs slaughtered STL decomp 'reseasonalised' method")

Household wealth (hh_budget).

With no seasonality and no consistent trend in household wealth across four countries, the Naive method appears to be the best choice.

hh_budget %>%
  select(Country, Year, Wealth) %>%
  autoplot(Wealth) +
  labs(title = "Household wealth")

The naive forecast uses the last value for each respective country as its forecast, using appropriate confidence intervals for each country based on its wealth scale.

hh_budget %>%
  select(Country, Year, Wealth) %>%
  model(NAIVE(Wealth)) %>%
  forecast(h = 5) %>%
  autoplot(hh_budget) +
  labs(title = "5-year Forecast of Household wealth by Country") +
  facet_wrap(Country ~ ., scales = "free_y")

Australian takeaway food turnover (aus_retail).

Each state looks like it has an upward trend with seasonality included, but the trend is more prominent than the seasonality, which leans towards utilizing a drift method to capture the trend in each State.

aus_retail %>%
  filter(Industry == "Cafes, restaurants and takeaway food services") %>%
  autoplot(Turnover) +
  facet_wrap(State ~ ., scales = "free_y") + 
  theme(legend.position = "none") +
  labs(title = "Australian takeaway food turnover by State")

It appears that the drift forecasts are trending for each territory.

aus_retail %>%
  filter(Industry == "Cafes, restaurants and takeaway food services") %>%
  model(RW(Turnover ~ drift())) %>%
  forecast(h = 60) %>%
  autoplot(aus_retail) +
  facet_wrap(State ~ ., nrow = 4, ncol = 2, scales = "free_y") + 
  theme(legend.position = "none") +
  labs(title = "Australian takeaway food turnover by State")

Looking more closely just at Victoria, it is easiser to see that the drift method does indeed reflect the trend.

aus_retail %>%
  filter(State == "Victoria", Industry == "Cafes, restaurants and takeaway food services") %>%
  model(RW(Turnover ~ drift())) %>%
  forecast(h = 60) %>%
  autoplot(aus_retail) +
  facet_wrap(State ~ ., nrow = 4, ncol = 2, scales = "free_y") + 
  theme(legend.position = "none") +
  labs(title = "Australian takeaway food turnover by State")

  1. Use the Facebook stock price (data set gafa_stock) to do the following:
  1. Produce a time plot of the series.

Facebook (FB) stock price (Close price) peaked on July 25, 2018, after which it slid downhill.

ggplotly(gafa_stock %>%
  filter(Symbol == "FB") %>%
  autoplot(Close) +
  labs(title = "Facebook (FB) stock price"))
  1. Produce forecasts using the drift method and plot them.

I had to reindex the days since the Date index was irregular. Then I was able to plot a 30-day forecast using the drift method.

fb <- gafa_stock %>%
  filter(Symbol == "FB") %>%
  mutate(regular_day = row_number()) %>%
  update_tsibble(index = regular_day, regular = TRUE)

fb %>%
  model(RW(Close ~ drift())) %>%
  forecast(h = 30) %>%
  autoplot(fb) +
  labs(title = "Facebook (FB) stock price with 30-day drift forecast")

  1. Show that the forecasts are identical to extending the line drawn between the first and last observations.

It is possible to add a linear segment between the first and last observations in the original data, which aligns with the forecast trend line.

fb %>%
  model(RW(Close ~ drift())) %>%
  forecast(h = 30) %>%
  autoplot(fb) +
  labs(title = "Facebook (FB) stock price with 30-day drift forecast") +
  geom_segment(aes(x=first(regular_day), y=first(Close), xend=last(regular_day), yend=last(Close)))

  1. Try using some of the other benchmark functions to forecast the same data set. Which do you think is best? Why?

It is “cool” that we can graph multiple forecast types on a single plot. The Mean model (in red) is taking the mean of the entire series, which rose then fell, so the mean forecast is not even close to where the last recorded price falls. The Naive model (in green) takes the last recorded point as the forecast and extends it into the future, which at least continues the plot from where it left off, but does not capture any of the trend. The Seasonal Naive (Random Walk version, in blue) as noted above captures the trend and actually extends it from where it leaves off in the observed values. Therefore, I think the seasonal naive/random walk is the best in this case, since it captures the historical trend and extends it.

fb %>%
  model(Mean = MEAN(Close),
        `Naive` = NAIVE(Close),
        `Seasonal Naive` = RW(Close ~ drift())) %>%
  forecast(h = 90) %>%
  autoplot(fb, level=NULL) +
  labs(title = "Facebook (FB) stock price with Mean, Naive and Seasonal Naive") +
  guides(colour = guide_legend(title = "Forecast"))

  1. 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. The following code will help.

The ACF lag 4 is significant, which makes sense given that the time period is quarterly but the data is recorded monthly.

# 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()

To check if the residuals look like white noise we can use the Ljung-Box test.

fit_resid <- fit %>%
  augment() %>% 
  features(.innov, ljung_box, lag = 8, dof = 1)

fit_resid
## # A tibble: 1 x 3
##   .model       lb_stat lb_pvalue
##   <chr>          <dbl>     <dbl>
## 1 SNAIVE(Beer)    32.3 0.0000362

Here is a plot of the forecast.

# Look a some forecasts
fit %>% forecast() %>% 
  autoplot(recent_production) +
  labs(title = "Forecast of Quarterly Australian beer production")

What do you conclude?

The pvalue from the Ljung-Box test is small, meaning it is significant. This suggests that the residuals do not look like white noise and the forecast model can possibly be improved.

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

Looking at Australian exports first, here is a plot of both methods. Given the trend in the Australian exports data it appears that the Seasonal Naive forecast is most appropriate to capture and reflect the trend.

# Look a some forecasts
global_economy %>%
  filter(Country == "Australia") %>%
  model(`Naive` = NAIVE(Exports),
        `Seasonal Naive` = RW(Exports ~ drift())) %>%
  forecast(h = 10) %>%
  autoplot(global_economy, level = NULL) +
  labs(title = "Australian exports with Naive and Seasonal Naive methods") +
  guides(colour = guide_legend(title = "Forecast"))

Checking the residuals using the Seasonal Naive method only the first lag is significant, but that is probably because the series is annual.

# Extract data of interest
Aus_exports <- global_economy %>%
  filter(Country == "Australia") %>%
  select(-Growth) %>%
  fill_gaps

# Define and estimate a model
fit_aus_exp <- Aus_exports %>% model(RW(Exports ~ drift()))

# Look at the residuals
fit_aus_exp %>% gg_tsresiduals()

Checking for further significance, the p-value is insignificant in that it is > 0.05, meaning the residuals are more like white noise.

aus_export_resid <- Aus_exports %>% 
  model(RW(Exports ~ drift())) %>%
  augment()

aus_export_resid %>%
  features(.innov, ljung_box, lag=10, dof=0)
## # A tibble: 1 x 4
##   Country   .model                lb_stat lb_pvalue
##   <fct>     <chr>                   <dbl>     <dbl>
## 1 Australia RW(Exports ~ drift())    16.4    0.0896

Next looking at Australian Bricks using the NAIVE method, again the Seasonal Naive forecast line aligns better with the overall trend from the first data point to the last in the observed data.

aus_production %>%
  filter(!is.na(Bricks)) %>%
  model(`Naive` = NAIVE(Bricks),
        `Seasonal Naive` = RW(Bricks ~ drift())) %>%
  forecast(h = 10) %>%
  autoplot(aus_production, level = NULL) +
  labs(title = "Australian brick production with Naive and Seasonal Naive methods") +
  guides(colour = guide_legend(title = "Forecast"))

Checking the residuals using the Seasonal Naive method with a random walk for drift, all the even lags are significant, and the residuals appear bimodal and not centered around 0.

# Extract data of interest
aus_bricks <- aus_production %>%
  filter(!is.na(Bricks))

# Define and estimate a model
fit_aus_bricks <- aus_bricks %>% model(RW(Bricks ~ drift()))

# Look at the residuals
fit_aus_bricks %>% gg_tsresiduals()

Checking for further significance, the p-value is significant in that it is 0, meaning the model can probably be further improved with more advanced methods.

aus_brick_resid <- aus_bricks %>% 
  model(RW(Bricks ~ drift())) %>%
  augment()

aus_brick_resid %>%
  features(.innov, ljung_box, lag=10, dof=0)
## # A tibble: 1 x 3
##   .model               lb_stat lb_pvalue
##   <chr>                  <dbl>     <dbl>
## 1 RW(Bricks ~ drift())    360.         0
  1. For your retail time series (from Exercise 8 in Section 2.10):

This again refers to the aus_retail dataset showing retail turnover from 1982 to 2019. The original series is shown once more below.

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

autoplot(myseries,Turnover) +
  labs(title = "Turnover by Month")

  1. Create a training dataset consisting of observations before 2011 using

The training set will contain all data prior to 2011.

myseries_train <- myseries %>%
  filter(year(Month) < 2011)
  1. Check that your data have been split appropriately by producing the following plot.

The training data is shown in red (prior to 2011) and the test data is what remains in black.

autoplot(myseries, Turnover) +
  autolayer(myseries_train, Turnover, colour = "red")

  1. Fit a seasonal naïve model using SNAIVE() applied to your training data (myseries_train).

Only the training data is included in the model, seen as the dataset myseries_train.

fit <- myseries_train %>%
  model(SNAIVE(Turnover))
  1. Check the residuals. Do the residuals appear to be uncorrelated and normally distributed?

The residuals appear to be correlated according to the lag plot, and the residuals appear somewhat normally distributed with an extreme outiler on the upper tail.

fit %>% gg_tsresiduals()

  1. Produce forecasts for the test data

The forecast produces sesonality but without trend or drift. This locks in the forecast at the level of the last data in the training set.

fc <- fit %>%
  forecast(new_data = anti_join(myseries, myseries_train))
## Joining, by = c("State", "Industry", "Series ID", "Month", "Turnover")
fc %>% autoplot(myseries)

  1. Compare the accuracy of your forecasts against the actual values.

Without accounting for trend, the divergence between the forecast and the actual values increases over time both in terms of fit…

fit %>% accuracy()
## # A tibble: 1 x 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 South A~ Clothin~ SNAIV~ Trai~  1.91  5.64  3.96  3.48  11.6     1     1 0.773

and even more so in terms of accuracy.

fc %>% accuracy(myseries)
## # A tibble: 1 x 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~ Sout~ Clothin~ Test  -12.7  15.6  13.2 -24.0  24.5  3.32  2.77 0.762
  1. How sensitive are the accuracy measures to the amount of training data used?

The amount of training data used is not necessarily as important as the properties of the training set compared to the testing set. If the training set has far less variance than the testing sets then it is likely that the model based on the training data will be overfit and fail when applied to the test data. Perhaps this is where cross-validation methods could prove more useful in training and testing the models.