library(fpp3)
## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──
## ✔ tibble      3.2.1     ✔ tsibble     1.1.4
## ✔ dplyr       1.1.3     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.0     ✔ feasts      0.3.1
## ✔ lubridate   1.9.3     ✔ fable       0.3.3
## ✔ ggplot2     3.4.4     ✔ fabletools  0.3.4
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
library(tsibble)

Do exercises 5.1, 5.2, 5.3, 5.4 and 5.7 in the Hyndman book.

  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)

global_economy |> filter(Country=='Australia') |>
    autoplot(Population)

Given the lack of seasonality in the plot it doesn’t make sense to apply a seasonal benchmark method.

econ_aus_sub <- global_economy |> dplyr::filter(Country=='Australia' & Year<=2015)
econ_aus_fit <- econ_aus_sub |>
    model(naive = NAIVE(Population), drift = RW(Population ~ drift()))

econ_fc <- econ_aus_fit |>
    forecast(h = 2)

econ_fc |>
  autoplot(global_economy |> dplyr::filter(Country=='Australia' & Year>'2000'),level=NULL) +
  labs(y = "Population",title = "Forecasts for annual Australian Population")

A visual examination of the data indicates that the drift method appears to be a better method to forecast with this time series.

Let’s see if the accuracy metrics concur with this eye test:

accuracy(econ_fc,global_economy |> dplyr::filter(Country=='Australia'))
## # A tibble: 2 × 11
##   .model Country   .type      ME    RMSE     MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <fct>     <chr>   <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 drift  Australia Test  183879. 196987. 183879. 0.751 0.751 0.745 0.764  -0.5
## 2 naive  Australia Test  554087  587088. 554087  2.26  2.26  2.25  2.28   -0.5

Consistent with that expectation the varied score metrics all rate the drift method to be a better means of forecasting for this particular dataset; however, there could definitely be improvement by trying to find alternatives.

Bricks (aus_production)

aus_production
## # A tsibble: 218 x 7 [1Q]
##    Quarter  Beer Tobacco Bricks Cement Electricity   Gas
##      <qtr> <dbl>   <dbl>  <dbl>  <dbl>       <dbl> <dbl>
##  1 1956 Q1   284    5225    189    465        3923     5
##  2 1956 Q2   213    5178    204    532        4436     6
##  3 1956 Q3   227    5297    208    561        4806     7
##  4 1956 Q4   308    5681    197    570        4418     6
##  5 1957 Q1   262    5577    187    529        4339     5
##  6 1957 Q2   228    5651    214    604        4811     7
##  7 1957 Q3   236    5317    227    603        5259     7
##  8 1957 Q4   320    6152    222    582        4735     6
##  9 1958 Q1   272    5758    199    554        4608     5
## 10 1958 Q2   233    5641    229    620        5196     7
## # ℹ 208 more rows
aus_production |> autoplot(Bricks)
## Warning: Removed 20 rows containing missing values (`geom_line()`).

There is clearer seasonality for brick production than the prior dataset and the seasonal naive forecast is what is expected to do the best from these benchmark methods.

aus_prod_fit <- aus_production |> filter(year(Quarter)<='2000') |> model(naive = NAIVE(Bricks), drift = RW(Bricks ~ drift()),snaive= SNAIVE(Bricks))

aus_prod_fc <- aus_prod_fit |> forecast(h=20)

aus_prod_fc |>
    autoplot(aus_production |> filter(year(Quarter)<='2005'),level=NULL)
## Warning: Removed 2 rows containing missing values (`geom_line()`).

It’s rather hard to tell visually which benchmark method is doing a decent job forecasting this time series.

accuracy(aus_prod_fc,aus_production)
## # A tibble: 3 × 10
##   .model .type     ME  RMSE   MAE    MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 drift  Test    5.14  33.8  28.0  0.450  7.38 0.789 0.698 0.292
## 2 naive  Test   15.2   39.4  34.7  2.95   8.92 0.978 0.814 0.350
## 3 snaive Test  -23.3   48.7  36.7 -6.90  10.1  1.04  1.01  0.332

After comparing the metrics across the 3 different forecasting methods it appears that the drift method has the lowest scores across most of the prevalent error statistics and therefore we would choose to use that out of the three options.

NSW Lambs (aus_livestock)

lambs <- aus_livestock|>  filter(Animal=='Lambs' & State=='New South Wales')
lambs |> autoplot(Count)

lambs |> ACF(Count) |> autoplot()

pre_lambs <- lambs |> filter(year(Month)<=2015)
   
lamb_fit <- pre_lambs |>
    model(naive = NAIVE(Count), drift = RW(Count ~ drift()),snaive= SNAIVE(Count))

lamb_fc <- lamb_fit |>
    forecast(h=36)

lamb_fc |>
    autoplot(lambs |> filter(year(Month)>='2010'),level=NULL)

Based on the forecasted graphic it appears that seasonal naive method would be best of the 3 basic approaches to forecast the data.

Do the accuracy statistics agree with the eye test:

accuracy(lamb_fc,lambs)
## # A tibble: 3 × 12
##   .model Animal State  .type     ME   RMSE    MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <chr>  <fct>  <fct>  <chr>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 drift  Lambs  New S… Test  21655. 45021. 37173.  4.21  8.81 0.878 0.805 0.0335
## 2 naive  Lambs  New S… Test  15242. 42843. 34008.  2.65  8.17 0.804 0.766 0.0580
## 3 snaive Lambs  New S… Test  -2650. 46290. 36472. -1.65  9.02 0.862 0.828 0.0398

Interestingly enough the naive method which is equal to the last observation is best! This might be due to the fact that the seaonality and cycles are moving with somewhat high variance and perhaps the stability of the last forecast is the least of the worst.

Household wealth (hh_budget).

hh_budget
## # A tsibble: 88 x 8 [1Y]
## # Key:       Country [4]
##    Country    Year  Debt     DI Expenditure Savings Wealth Unemployment
##    <chr>     <dbl> <dbl>  <dbl>       <dbl>   <dbl>  <dbl>        <dbl>
##  1 Australia  1995  95.7 3.72          3.40   5.24    315.         8.47
##  2 Australia  1996  99.5 3.98          2.97   6.47    315.         8.51
##  3 Australia  1997 108.  2.52          4.95   3.74    323.         8.36
##  4 Australia  1998 115.  4.02          5.73   1.29    339.         7.68
##  5 Australia  1999 121.  3.84          4.26   0.638   354.         6.87
##  6 Australia  2000 126.  3.77          3.18   1.99    350.         6.29
##  7 Australia  2001 132.  4.36          3.10   3.24    348.         6.74
##  8 Australia  2002 149.  0.0218        4.03  -1.15    349.         6.37
##  9 Australia  2003 159.  6.06          5.04  -0.413   360.         5.93
## 10 Australia  2004 170.  5.53          4.54   0.657   379.         5.40
## # ℹ 78 more rows
hh_budget |> autoplot(Wealth)

Let’s focus on Australia to be consistent with other time series reviewed in this homework.

gg_lag(hh_budget |> filter(Country=='Australia'),Wealth,geom='point',lags=1:20)

hh_budget |> filter(Country=='Australia') |> ACF(Wealth,lag_max=30) |>autoplot()

There are definitely stronger correlations every 5 to 10 years, but there isn’t much of seasonality. Therefore we will exclude seasonal naive from this analysis.

hh_sub <- hh_budget |> filter(Country=='Australia' & Year<='2006')
hh_fit <- hh_sub |>
    model(naive=NAIVE(Wealth),RW(Wealth ~ drift()))

hh_fc <- hh_fit |>
    forecast(h=10)

hh_fc |>
    autoplot(hh_budget |> filter(Country=='Australia'),level=NULL)

The naive method would appear to be the best forecast method out of the basic benchmarks.

accuracy(hh_fc,hh_budget |> filter(Country=='Australia'))
## # A tibble: 2 × 11
##   .model          Country .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <chr>           <chr>   <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 RW(Wealth ~ dr… Austra… Test  -87.7  92.5  87.7 -24.7  24.7  9.12  7.71 0.0874
## 2 naive           Austra… Test  -41.8  52.2  44.9 -12.3  13.0  4.67  4.36 0.276

The various accuracy metrics are consistent with the eye test that show the naive method is much preferred to the drift one.

Australian takeaway food turnover (aus_retail).

aus_retail
## # A tsibble: 64,532 x 5 [1M]
## # Key:       State, Industry [152]
##    State                        Industry           `Series ID`    Month Turnover
##    <chr>                        <chr>              <chr>          <mth>    <dbl>
##  1 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Apr      4.4
##  2 Australian Capital Territory Cafes, restaurant… A3349849A   1982 May      3.4
##  3 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Jun      3.6
##  4 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Jul      4  
##  5 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Aug      3.6
##  6 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Sep      4.2
##  7 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Oct      4.8
##  8 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Nov      5.4
##  9 Australian Capital Territory Cafes, restaurant… A3349849A   1982 Dec      6.9
## 10 Australian Capital Territory Cafes, restaurant… A3349849A   1983 Jan      3.8
## # ℹ 64,522 more rows
aus_retail |> filter(Industry=='Takeaway food services' & State=='New South Wales') |> autoplot(Turnover) + labs(title='New South Wales Takeaway Food')

aus_turn_sub <- aus_retail |> filter(Industry=='Takeaway food services' & State=='New South Wales')
aus_turn_sub |> ACF(Turnover) |> autoplot()

There is a fairly strong correlation across lagged months, but the standard yearly period seems to have a slightly higher correlation.

turn_fit <- aus_turn_sub |> filter(year(Month)<='2015') |>
    model(naive=NAIVE(Turnover),snaive=SNAIVE(Turnover),drift=RW(Turnover~drift()))

turn_fc <- turn_fit |> forecast(h=36)

turn_fc |> autoplot(aus_turn_sub,level=NULL)

Given the upward trend it looks as though the drift method may be best at approximating the patterns in the New South Wales takeaway turnover.

accuracy(turn_fc,aus_turn_sub)
## # A tibble: 3 × 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 drift  New So… Takeawa… Test  -11.1   36.4  27.6 -2.50  5.24  1.25  1.22 0.447
## 2 naive  New So… Takeawa… Test    9.61  40.8  32.8  1.22  6.02  1.49  1.37 0.533
## 3 snaive New So… Takeawa… Test   76.3   78.8  76.3 13.9  13.9   3.45  2.64 0.685

The most popular accuracy metrics to evaluate forecasts appear to be at their lowest for the drift method which is consistent with the graphic representation of the point forecasts.

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

Produce a time plot of the series.

gafa_stock |> filter(Symbol=='FB') |> autoplot(Close)

Produce forecasts using the drift method and plot them.

fb_ts <- gafa_stock |> filter(Symbol=='FB') |>
    mutate(trading_day = row_number()) |>
    update_tsibble(index = trading_day, regular = TRUE)

fb_jan_2018 <- fb_ts |>
  filter(yearmonth(Date) >= yearmonth("2017 Dec"))

fb_pre_2018 <- fb_ts |> filter(Symbol=='FB' & year(Date)<= 2017) 
    
fb_pre_2018 |>
    model(RW(Close ~ drift())) |> 
    forecast(h=30) |> 
    autoplot(fb_pre_2018) +
  labs(y = "Price ($)",title = "Facebook closing stock prices",subtitle='30 day forecast')

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

fb_pre_2018 |>
    model(RW(Close ~ drift())) |> 
    forecast(h=30) |> 
    autoplot(fb_pre_2018) +
    #autolayer()
    geom_segment(aes(x = 0, y = pull(fb_pre_2018 |> filter(trading_day==1),Close), xend = 1007, yend = pull(fb_pre_2018 |> filter(trading_day==1007),Close),linetype=2 ),,colour='red',linetype=2) +
  labs(y = "Price ($)",title = "Facebook closing stock prices",subtitle='30 day forecast')

The difference from the first to last points has the exact same slope as the forecasted point as shown in the textbook.

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

fb_pre_2018 |>
    model(naive=NAIVE(Close),mean=MEAN(Close)) |> 
    forecast(h=30) |>
    autoplot(fb_pre_2018 |> filter(trading_day>500),level=NULL)

Given the longer term trend of the data it would appear that the best method from a visual analysis is the drift method as it captures the change that occurs from the beginning and the end.

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

# Look a some forecasts
fit |> forecast() |> autoplot(recent_production)

fit |> gg_tsresiduals()
## Warning: Removed 4 rows containing missing values (`geom_line()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).
## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

What do you conclude?

There appears to be limited autocorrelation among the residuals given that there is only one value at lag 4 that has some significance. In some ways this is expected to be highest lagged correlation based on the seasonal time periods and Rob Hyndman mentioned 1 in 20 having some correlation is not that unusual given it equates to 5%. The histogram seems somewhat normally distributed although the mean value is not as common as some higher and lower residuals that may imply some bias and it is not in fact centered at zero. The innovation residuals show some variation across time although only have two large variants at -40.

  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.
aus_ex <- global_economy |> filter(Country=='Australia')
aus_ex |> autoplot(Exports)

aus_ex |> ACF(Exports) |> autoplot()

Since there isn’t a clear season based on the year time frame it makes sense to use the naive method as a forecast.

aus_fit_pre_2015 <- aus_ex |> filter(Year<'2010')
aus_fit <- aus_fit_pre_2015 |> model(naive=NAIVE(Exports))
aus_ex_fc <- aus_fit |> forecast(h=7)

aus_ex_fc |>
    autoplot(aus_ex,level=NULL)

aus_fit |> gg_tsresiduals()
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

There doesn’t seem to be much autocorrelation between the residuals across the lagged time frame in the ACF plot. The histogram of residual values looks fairly normallly distributed although there isn’t much spread overall. The residuals look to be mostly homoskedastic across the time period although the later time frames look to be slightly more variable.

aus_production |> autoplot(Bricks)
## Warning: Removed 20 rows containing missing values (`geom_line()`).

There is seasonality in this data set with quarterly periods.

brick_fit <- aus_production |> filter(year(Quarter)<'1995') |>
    model(naive = NAIVE(Bricks),snaive=SNAIVE(Bricks))

brick_fc <- brick_fit |> forecast(h=40)

brick_fc |>
    autoplot(aus_production,level=NULL)
## Warning: Removed 20 rows containing missing values (`geom_line()`).

It’s hard to say from the graphic that the forecasts are either that good, but the seasonal naive appears better.

Let’s compare accuracy metrics for each method:

accuracy(brick_fc,aus_production)
## # A tibble: 2 × 10
##   .model .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 naive  Test  -81.8  89.3  81.8 -21.9  21.9  2.30  1.84 0.513
## 2 snaive Test  -73.1  79.9  73.5 -19.2  19.3  2.07  1.65 0.648

This results in a similar conclusion that the seasonal naive is preferred to the naive method for Australian Brick time series.

brick_fit %>% dplyr::select(naive) |> gg_tsresiduals()
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

A brief review of the residuals would appear to indicate that there is significant autocorrelation between the residual values as most of the lagged values have some significance. The spread of the residuals is a little skewed on the left ned and is also not centered at zero. Lastly, the variance is fairly herteroskedastic with a mixture of values above and below zero but appears to be a little unusual from 1975 - 1985. This likely speaks to the poor results from the forecasted metric overall and why it is likely not practical to use either method for predictive purposes.

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

Create a training dataset consisting of observations before 2011 using

myseries <- aus_retail |> filter(`Series ID` == 'A3349335T')
myseries_train <- aus_retail |>
  filter(year(Month) < 2011 & `Series ID` == 'A3349335T')

Check that your data have been split appropriately by producing the following plot.

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

The training set appears to be accurate based on the execution of the prior code and review of the graphic.

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

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

Check the residuals.

fit |> gg_tsresiduals()
## Warning: Removed 12 rows containing missing values (`geom_line()`).
## Warning: Removed 12 rows containing missing values (`geom_point()`).
## Warning: Removed 12 rows containing non-finite values (`stat_bin()`).

The residuals appear to violate some of the ideal conditions for optimal forecasts. There is heteroskedastic variance across the time frame, non-zero correlation among the residuals, and does not have a normally distributed shape.

Do the residuals appear to be uncorrelated and normally distributed?

The residuals appear highly correlated with most of the lagged correlations in the ACF as they are above the lack of significance line. The shape of the residuals is not normally distributed either as there is considerable right skew and it is centered around 50 rather than zero. The variance is also not consistent across the range of values and is more likely heteroskedastic particularly during the financial crisis in 2008 - 2009.

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 New Sou… Superma… SNAIV… Trai…  61.6  72.2  61.7  6.39  6.40     1     1 0.602
fc |> 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… New … Superma… Test   409.  480.  409.  15.9  15.9  6.64  6.65 0.928

The seasonal naive forecast are substantially higher values across the board and from the prior plot it appears that a prior year trend is not sufficient given the longe term trend growth cycle that occurs in this time series.

How sensitive are the accuracy measures to the amount of training data used?

The accuracy is very sensitive to the training data although this benchmark method using prior season’s data may not be the best type of forecasting strategy given the fairly consistent growth trend in the data.