library(tidyverse)
library(fpp3)
library(fable)Assignment 3: Forecasting
Assignment 3: Do exercises 5.1, 5.2, 5.3, 5.4 and 5.7 in the Hyndman book. Link to Chapter 5 exercises for reference.
Loading the Data
5.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
australian_population <- global_economy |>
filter(Country == 'Australia') |>
select(Year, Population)
autoplot(australian_population, Population)
Based on the data and our choice of methods, the random walk with drift will be our most appropriate method for the Australian population data. This is because random walk with drift can incorporate the constant upward trend we see (and we do not see seasonality in the data.)
aus_population_forecast <- australian_population |>
model(RW(Population ~ drift())) |>
forecast(h = 6) |>
autoplot(australian_population)
aus_population_forecast +
labs(
title = "Random walk with Drift Method for Australian Population Forecast")
Bricks (aus_production)
australian_bricks <- aus_production |>
select (Quarter, Bricks) |>
filter(!is.na(Bricks))
autoplot(australian_bricks, Bricks)
There is clearly seasonality in this data. SNAIVE will be most appropriate because we can see seasonal patterns in the quarterly data. The specification of lag is, as we know from the example in section 5.2, optional here, but we will specify the one-year lag nonetheless, with the example in 5.2 as guidance.
australian_bricks_forecast <- australian_bricks |>
model(SNAIVE(Bricks ~ lag("year"))) |>
forecast(h = 12) |>
autoplot(australian_bricks)
australian_bricks_forecast +
labs(title = "Seasonal Naive Method for Australian Bricks Production")
NSW Lambs aus_livestock
NSW_Lambs <- aus_livestock |>
filter(State == "New South Wales", Animal == "Lambs")
autoplot(NSW_Lambs, Count)
This one is a bit trickier to determine: is Seasonal Naive or Random Walk with Drift more appropriate? It seems like there might be some seasonality happening in this monthly tsibble. Let’s look into it more, and if it doesn’t seem like there’s relevant seasonality we will use the Random Walk with Drift method.
gg_season(NSW_Lambs)Plot variable not specified, automatically selected `y = Count`
gg_subseries(NSW_Lambs)Plot variable not specified, automatically selected `y = Count`
We used gg_subseries and gg_season to determine if SNAIVE will be a valid method to model and produce a forecast. From what I can see, the seasonality is definitely ‘noisy’, but it’s there. There is a visible seasonal pattern in terms of certain months showing consistently higher values across different years. While it’s not as visible as the Australian bricks production, for example, we can se that there are months with higher or lower rates of lamb slaughter over time. If we used Random Walk with Drift, the model wouldn’t account for this seasonality. So, due to the clear (while imperfect and noisy) existence of some seasonal patterns over time, we will use SNAIVE as the most appropriate.
NSW_Lambs_forecast <- NSW_Lambs |>
model(SNAIVE(Count)) |>
forecast(h = 12) |>
autoplot(NSW_Lambs)
NSW_Lambs_forecast +
labs(title = "Seasonal Naive Method for Australian Lamb Slaughter")
Household wealth hh_budget
household_wealth <- hh_budget |>
select(Year, Wealth)
autoplot(household_wealth, Wealth)
Since this is a dataset we haven’t looked much into in previous coursework, for more context on what we’re looking at here, I ran ?hhbudget to ind out that hh_budget is an annual tsibble with Wealth as a percentage of net disposable income, from 1995-2016.
Based on the autoplot above and this information, we know that SNAIVE will not be reasonable. Additionally, NAIVE will not be reasonable because there is an upward trend visible for each of the countries in the dataset. So, we will choose Random Walk with Drift as the most appropriate of our choices.
household_wealth_forecast <- household_wealth |>
model(RW(Wealth ~ drift())) |>
forecast(h = 3) |>
autoplot(household_wealth)
household_wealth_forecast +
labs(title = "Random walk with drift method for Household wealth forecast")
Australian takeaway food turnover aus_retail
aus_retail_takeaway_food <- aus_retail |>
filter(Industry == "Cafes, restaurants and takeaway food services")
autoplot(aus_retail_takeaway_food, Turnover)
Like the NSW_Lambs data, we aren’t immediately sure if SNAIVE or Random Walk with Drift is most appropriate in this instance. (NAIVE would not be appropriate; we see there is an increasing trend.) Let’s look more closely.
gg_subseries(aus_retail_takeaway_food)Plot variable not specified, automatically selected `y = Turnover`
Here it appears that the trend is stronger than the seasonality – the turnover is growing over time with a more visible pattern than if we were looking at seasonality alone with SNAIVE. So, we will decide to use the random walk with drift method to model the data and forecast. Seasonality alone does not appear strong enough to choose it over the random walk with drift in this particular case.
aus_retail_takeaway_food_forecast <- aus_retail_takeaway_food |>
model(RW(Turnover ~ drift())) |>
forecast(h = 6) |>
autoplot(aus_retail_takeaway_food) +
facet_wrap(~ State, scales = "free_y")
aus_retail_takeaway_food_forecast +
labs(title = "Random walk with drift method for Australian takeaway forecast")
5.2
Use the Facebook stock price (data set gafa_stock) to do the following:
- Produce a time plot of the series.
facebook_stock_price <- gafa_stock |>
filter(Symbol == "FB") |>
mutate(trading_day = row_number()) |>
update_tsibble(index = trading_day, regular = TRUE)
autoplot(facebook_stock_price, Close) +
labs(title = "Facebook daily closing stock price from `gafa_stock`")
Like the example in section 5.2 of the book, we added a column for the row number because there isn’t data for trading on weekends and other non-trading days like holidays, and udpated the tsibble to have the index = trading_day knowing this step would be important for modeling and forecasting in following steps.
- Produce forecasts using the drift method and plot them.
I did this slightly different than in question 5.1 because of the next steps required. We’ll use this more ‘correct’ way to do the autoplot with autolayer going forward. But, my previous methods still worked for what they were needed for. I am now doing the ‘autoplot’ step separately from creating the forecast.
facebook_stock_price_forecast <- facebook_stock_price |>
model(RW(Close ~ drift())) |>
forecast(h = 36)
autoplot(facebook_stock_price, Close) +
autolayer(facebook_stock_price_forecast, series = "Forecast") +
labs(title = "Random Walk with Drift: Facebook Closing Stock Price Forecast",
x = "Trading Day", y = "Stock Price") +
theme_minimal()Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
Ignoring unknown parameters: `series`Warning in geom_line(mapping = without(mapping, "shape"), data =
unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
`series`
- Show that the forecasts are identical to extending the line drawn between the first and last observations.
First, we’ll compute the first and last observations of the Close closing stock price. Then, we’ll calculate the slope. Then, we’ll add a “drift slope” line by calculating the first value plus the (row number - 1) times the slope as a column in the tsibble. Then, we’ll plot.
first_obs <- facebook_stock_price |> slice(1)
last_obs <- facebook_stock_price |> slice(n())
x_0 <- first_obs$trading_day
y_0 <- first_obs$Close
x_n <- last_obs$trading_day
y_n <- last_obs$Close
drift_slope <- (y_n - y_0) / (x_n - x_0)
extended_x <- seq(x_0, x_n + 36, by = 1)
extended_y <- y_0 + drift_slope * (extended_x - x_0)
autoplot(facebook_stock_price_forecast) +
autolayer(facebook_stock_price, Close, color = "black") +
geom_line(aes(x = extended_x, y = extended_y), color = "violet") +
labs(title = "Drift Method Forecast:
Extending Line Between First and Last Observations")
- Try using some of the other benchmark functions to forecast the same data set. Which do you think is best? Why?
facebook_additional_benchmarks <- facebook_stock_price |>
model(
Mean = MEAN(Close),
Naive = NAIVE(Close),
Drift = RW(Close ~ drift())
) |>
forecast(h = 36) |>
autoplot(facebook_stock_price, level = NULL) +
labs(title = "Facebook closing stock price - multiple benchmark functions") +
guides(color = guide_legend(title = "Forecast"))
facebook_additional_benchmarks
Compared to the Naive and Mean model benchmarks, I think that the Drift method is the best - it provides the best forecast for the Facebook closing stock price. The Mean method (green line) is a mean of all past values, which doesn’t take in to account the upward or downward changes of the stock, making it pretty disconnected from recent changes in closing stock value. The Naive method (blue line) assumes the most recent observed value is best for the future closing price forecast, which would work well for series where there’s little change, but doesn’t make it effective for stock price fluctuations that are highly variable over time. So, the drift is best - it extends the trend established by the first and last observations (as we saw illustrated with the addition of the line in the previous part of this question), and allows for change over time.
5.3
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.
Extracting data of interest and defining and estimating a model:
aus_beer_production <- aus_production |>
filter(year(Quarter) >= 1992)
fit <- aus_beer_production |> model(SNAIVE(Beer))Looking at the residuals:
fit |> gg_tsresiduals()
Looking at some forecasts:
fit |> forecast() |> autoplot(aus_beer_production)
What do you conclude?
Based on the residual time series plot, residuals fluctuate randomly but do show some patterns. The ACF of residuals has some significant autocorrelations (bars crossing the blue threshold), which suggests the model hasn’t captured all the data structure. The histogram of residuals does not appear normal though it’s somewhat symmetric. This makes me conclude that the residuals do not look like white noise, so the SNAIVE model hasnt’ captured all the patterns in the data. Another model might do better.
Based on the forecast analysis, the SNAIVE model is appropriate, though, since beer production follows a strong seasonal pattern and the forecast plot shows the seasonality repeating into the future. However, the confidence intervals grow much taller in their coverage, indicating that the uncertainty increases over time. So, while SNAIVE captures the seasonality we can see, the correlated residuals indicates we may do better with a different model outside the scope of this chapter (as suggested when viewing a similar example in the book).
In general, the seasonal naive method (SNAIVE) provides a reasonable forecast for beer production in Australia, but the residuals analysis shows that we are not looking at white noise, which shows that the model doesn’t fully account for or explain the data.
5.4
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.
Australian Exports series from global_economy
global_economy_aus_exports <- global_economy |>
filter(Country == "Australia") |>
select(Year, Exports)
fit_global_economy_aus_exports <- global_economy_aus_exports |>
model(NAIVE(Exports))Looking at the residuals:
fit_global_economy_aus_exports |> gg_tsresiduals()
Looking at some forecasts:
fit_global_economy_aus_exports |> forecast() |>
autoplot(global_economy_aus_exports)
I’ll be a bit more brief with this overview of what I conclude. We used NAIVE rather than SNAIVE because there is no obious seasonality in this annual data. The residuals look like white noise, so the Naive method appears reasonable - there’s no strong evidence of trends not captured. The residuals appear to fluctuate randomy and most lags in the ACF plot stay within limits, and the residuals histogram appears normal centered around 0. The forecast analysis, however, shows that NAIVE is not the best choice. There is fluctuation in the data, and NAIVE method doesn’t account for this (it shows the forecast as a flat line into the future). A different model than NAIVE or SNAIVE would likely lead to improved forecasts.
Bricks series from aus_production
We already created a tsibble for the Bricks series from aus_production and will use it below. We can remember by looking at question 5.1 that there is seasonality in this data, so we’ll use SNAIVE.
fit_australian_bricks <- australian_bricks |>
model(SNAIVE(Bricks))Looking at the residuals:
fit_australian_bricks |> gg_tsresiduals()
Looking at some forecasts:
fit_australian_bricks |> forecast() |> autoplot(australian_bricks)
Here I would conclude based on the residuals that there appear to be patterns and that it is not white noise, but the SNAIVE method provides what appears to be a decent seasonal forecast. The SNAIVE model doesn’t appear to have fully explained the data and the presence of autocorrelation would indicate that an alternative model outside the scope of our options (NAIVE vs SNAIVE) would be an improvement. The SNAIVE method doesn’t fully account for the variability and trend patterns in the data.
5.7
For your retail time series (from Exercise 7 in Section 2.10):
- Create a training dataset consisting of observations before 2011 using
set.seed(5889)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))Autoplot of myseries for reference:
autoplot(myseries)Plot variable not specified, automatically selected `.vars = Turnover`
myseries_train <- myseries |>
filter(year(Month) < 2011)- Check that your data have been split appropriately by producing the following plot.
autoplot(myseries, Turnover) +
autolayer(myseries_train, Turnover, color = "red")
We can confirm the data have been split appropriately.
- Fit a seasonal naïve model using
SNAIVE()applied to your training data (myseries_train).
fit_myseries <- myseries_train |>
model(SNAIVE(Turnover))- Check the residuals.
fit_myseries |> gg_tsresiduals()
Do the residuals appear to be uncorrelated and normally distributed?
Looking at the above, the residuals do not appear to be uncorrelated – a larger-than-expeted amount of lags in the ACF plot show there is significant correlation. If we were looking at white noise, we wouldn’t see a ACF plot like this, which shows more than expected lags above the boundaries (dashed blue lines). This tells us that the residuals are correlated. The histogram does not appear with expected normality, as it is somewhat bell shaped but not symmetric in its tails and not centered around 0. Also, the innovation residuals plot doesn’t appear to be fully ‘random’ - it appears there are still some remaining patterns.
- Produce forecasts for the test data
fc <- fit_myseries |>
forecast(new_data = anti_join(myseries, myseries_train))Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc |> autoplot(myseries)
- Compare the accuracy of your forecasts against the actual values.
fit_myseries |> accuracy()# 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 Austral… Other r… SNAIV… Trai… 1.15 2.52 1.96 4.96 7.75 1 1 0.719fc |> accuracy(myseries)# 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… Aust… Other r… Test 1.56 7.83 6.33 0.902 14.0 3.24 3.11 0.930The model appears to have a reasonable accuracy - the MPE (mean percentage error) is 0.9%, which tells us that on average the forecast is just slightly over-predicting the actual values (less than 1%). We learned that the low percentage error would suggest that the model is making accurate predictions. The MAPE of 14% indicates that the forecasts could be improved - it shows that the forecast values deviate from the actual values by 14% on average (MAPE stands for Mean Absolute Percentage Error). However, the MASE (Mean Absolute Scaled Error) value = 3.2 - and we’d seek for this to be not greater than 1 if the model were more performant than the naive benchmark.
- How sensitive are the accuracy measures to the amount of training data used?
Using more training data improves forecast accuracy in general, but including outdated historical data can introduce patterns that no longer reflect current trends, potentially reducing forecast reliability. This dataset provides a clear example of this, because the predictions initially overshoot actual values and later fall significantly below them. This suggests that a smaller, more recent subset of training data might yield better forecasts, particularly if the downturn around 2013 is not expected to repeat in the future.
The accuracy of forecasts is highly sensitive to the amount of training data used, especially in datasets with shifting trends. When the overall trend is stable, the choice of training data length has a smaller impact on forecast accuracy. However, in datasets like this one, where trends fluctuate over time, selecting an appropriate training period becomes crucial. If the training and test data exhibit similar patterns, accuracy measures remain stable. But when trends vary significantly across time, the sensitivity of accuracy measures to training data selection becomes more notable.