# Loading required packages
library(fpp3)DATA 624 Homework 3
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:
global_economy |>
filter(Country == "Australia") |>
autoplot(Population)
The population is yearly and shows a steady upward trend with no seasonal pattern in this case a random walk with drift is the best benchmark
aus_pop <- global_economy |> filter(Country == "Australia")
aus_pop |>
model(RW(Population ~ drift())) |>
forecast(h = 10) |>
autoplot(aus_pop)
As shown in the graph the forecast extends the average historical increase of ~250 thousand people per year
aus_production |>
autoplot(Bricks)Warning: Removed 20 rows containing missing values or values outside the scale range
(`geom_line()`).
The Brick production shows quarterly with a clear repeating seasonal pattern and no stable long run trend
bricks <- aus_production |> filter(!is.na(Bricks))
bricks |>
model(SNAIVE(Bricks)) |>
forecast(h = "3 years") |>
autoplot(bricks)
A seasonal naive forecast was produced for brick production repeating the most recent year’s quarterly pattern for the next three years. The prediction intervals look wide and widen which reflects large variability in the series
lambs <- aus_livestock |>
filter(Animal == "Lambs", State == "New South Wales")
lambs |> autoplot(Count)
lambs |> gg_season(Count)
NSW is a monthly graph with no stable trend and weak season trends. A seasonal naive method would be best so we can retain the mild pattern.
lambs |>
model(SNAIVE(Count)) |>
forecast(h = "3 years") |>
autoplot(lambs)
Seasonal naive forecast for NSW labs repeat the most recent years pattern. The wide prediction intervals show that the series is noisy and it shows a single year of history is a weak basis for forecasting.
hh_budget |> autoplot(Wealth)
hh_budget |>
model(RW(Wealth ~ drift())) |>
forecast(h = 5) |>
autoplot(hh_budget)
Random walk with drift forecast for household wealth extend each country’s average annual change giving modest upward projections for all of the countries. The USA widen more with the USA given that its history is the most volatile.
takeaway <- aus_retail |>
filter(Industry == "Takeaway food services") |>
summarise(Turnover = sum(Turnover))
takeaway |> autoplot(Turnover)
takeaway |>
model(SNAIVE(Turnover)) |>
forecast(h = "3 years") |>
autoplot(takeaway)
For this forecast it shows food turnover repeat the last observed years pattern with tight intervals.
5.2
Use the Facebook stock price (data set gafa_stock) to do the following:
Produce a time plot of the series.
fb <- gafa_stock |> filter(Symbol == "FB")
fb |> autoplot(Close)
Produce forecasts using the drift method and plot them.
fb_ts <- fb |>
mutate(day = row_number()) |>
update_tsibble(index = day, regular = TRUE)
fb_ts |>
model(RW(Close ~ drift())) |>
forecast(h = 30) |>
autoplot(fb_ts)
The Drift forecast for FB’s closing price rise slightly from the last observation. The prediction intervals widen with the horizon and are far larger than the forecast slope.
Show that the forecasts are identical to extending the line drawn between the first and last observations.
first <- fb_ts |> slice(1)
last <- fb_ts |> slice(n())
slope <- (last$Close - first$Close) / (last$day - first$day)
fb_ts |>
model(RW(Close ~ drift())) |>
forecast(h = 30) |>
autoplot(fb_ts) +
geom_abline(intercept = first$Close - slope * first$day,
slope = slope, colour = "red", linetype = "dashed")
The drift forecast coincides with the straight line
Try using some of the other benchmark functions to forecast the same data set. Which do you think is best? Why?
fb_ts |>
model(
Mean = MEAN(Close),
Naive = NAIVE(Close),
Drift = RW(Close ~ drift())
) |>
forecast(h = 30) |>
autoplot(fb_ts, level = NULL)
The mean is the least suitable because it averages the full history and it is below the latest price. Naive and drift give nearly identical short horizon forecast, and naive is preferred because drift slope reflects the earlier upward run rather than the recent decline.
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. The following code will help.
recent_production <- aus_production |>
filter(year(Quarter) >= 1992)
fit <- recent_production |> model(SNAIVE(Beer))
fit |> gg_tsresiduals()Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_bin()`).Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_rug()`).
The seasonal naive residuals are centered near zero, it also shows a normal histogram though the ACF shows a strong negative spike at lag 4.
fit |> forecast() |> autoplot(recent_production)
Given that the ACF has a large negative spike at lag 4 and a marginal one at lag1 it indicates that the method transmits the noise of a single past quarter into each forecast.
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.
aus_exports <- global_economy |> filter(Country == "Australia")
fit_exp <- aus_exports |> model(NAIVE(Exports))
fit_exp |> gg_tsresiduals()Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 1 row containing non-finite outside the scale range
(`stat_bin()`).Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_rug()`).
fit_exp |> forecast(h = 5) |> autoplot(aus_exports)
The forecast for Australian exports are flat at the last observed value, with prediction interval that widen in population.
bricks <- aus_production |> filter(!is.na(Bricks))
fit_bricks <- bricks |> model(SNAIVE(Bricks))
fit_bricks |> gg_tsresiduals()Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_bin()`).Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_rug()`).
fit_bricks |> forecast(h = "3 years") |> autoplot(bricks)
Because the residuals are strongly correlated the prediction intervals probably understate the true uncertainty.
5.7
For your retail time series (from Exercise 7 in Section 2.10):
set.seed(12345678)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
myseries |> autoplot(Turnover)
Create a training dataset consisting of observations before 2011 using
Check that your data have been split appropriately by producing the following plot.
myseries_train <- myseries |>
filter(year(Month) < 2011)
autoplot(myseries, Turnover) +
autolayer(myseries_train, Turnover, colour = "red")
Fit a seasonal naïve model using SNAIVE() applied to your training data
Check the residuals.
fit <- myseries_train |>
model(SNAIVE(Turnover))
fit |> gg_tsresiduals()Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 12 rows containing non-finite outside the scale range
(`stat_bin()`).Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_rug()`).
The residuals look strongly autocorrelated wit an ACF of 0.75 in the lag1 and it goes down slowly.
Produce forecasts for the test data
fc <- fit |>
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 |> 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 Norther… Clothin… SNAIV… Trai… 0.439 1.21 0.915 5.23 12.4 1 1 0.768fc |> 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… Nort… Clothin… Test 0.836 1.55 1.24 5.94 9.06 1.36 1.28 0.601On the training set the seasonal naive model has an RMSE 1.21 and an MAE 0.9, and on the test set it rises to 1.55 and 1.24 so the test error is larger.
How sensitive are the accuracy measures to the amount of training data used?
starts <- c(1989, 1996, 2001, 2006)
results <- lapply(starts, function(s) {
train <- myseries |> filter(year(Month) >= s, year(Month) < 2011)
train |>
model(SNAIVE(Turnover)) |>
forecast(new_data = anti_join(myseries, myseries_train)) |>
accuracy(myseries) |>
mutate(start = s)
})Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`bind_rows(results) |> select(start, RMSE, MAE, MAPE)# A tibble: 4 × 4
start RMSE MAE MAPE
<dbl> <dbl> <dbl> <dbl>
1 1989 1.55 1.24 9.06
2 1996 1.55 1.24 9.06
3 2001 1.55 1.24 9.06
4 2006 1.55 1.24 9.06