`library(fpp3)`

# Exercise Week 5: Solutions

# fpp3 5.10, Ex 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`

)- Bricks (
`aus_production`

)- NSW Lambs (
`aus_livestock`

)- Household wealth (
`hh_budget`

)- Australian takeaway food turnover (
`aus_retail`

)

### Australian population

```
|>
global_economy filter(Country == "Australia") |>
autoplot(Population)
```

Data has trend and no seasonality. Random walk with drift model is appropriate.

```
|>
global_economy filter(Country == "Australia") |>
model(RW(Population ~ drift())) |>
forecast(h = "10 years") |>
autoplot(global_economy)
```

### Australian clay brick production

```
|>
aus_production filter(!is.na(Bricks)) |>
autoplot(Bricks) +
labs(title = "Clay brick production")
```

This data appears to have more seasonality than trend, so of the models available, seasonal naive is most appropriate.

```
|>
aus_production filter(!is.na(Bricks)) |>
model(SNAIVE(Bricks)) |>
forecast(h = "5 years") |>
autoplot(aus_production)
```

### NSW Lambs

```
<- aus_livestock |>
nsw_lambs filter(State == "New South Wales", Animal == "Lambs")
|>
nsw_lambs autoplot(Count)
```

This data appears to have more seasonality than trend, so of the models available, seasonal naive is most appropriate.

```
|>
nsw_lambs model(SNAIVE(Count)) |>
forecast(h = "5 years") |>
autoplot(nsw_lambs)
```

### Household wealth

```
|>
hh_budget autoplot(Wealth)
```

Annual data with trend upwards, so we can use a random walk with drift.

```
|>
hh_budget model(RW(Wealth ~ drift())) |>
forecast(h = "5 years") |>
autoplot(hh_budget)
```

### Australian takeaway food turnover

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

This data has strong seasonality and strong trend, so we will use a seasonal naive model with drift.

```
|>
takeaway model(SNAIVE(Turnover ~ drift())) |>
forecast(h = "5 years") |>
autoplot(takeaway)
```

This is actually not one of the four benchmark methods discussed in the book, but is sometimes a useful benchmark when there is strong seasonality and strong trend.

The corresponding equation is \hat{y}_{T+h|T} = y_{T+h-m(k+1)} + \frac{h}{T-m}\sum_{t=m+1}^T(y_t - y_{t-m}), where m=12 and k is the integer part of (h-1)/m (i.e., the number of complete years in the forecast period prior to time T+h).

# fpp3 5.10, Ex 2

Use the Facebook stock price (data set

`gafa_stock`

) to do the following:

- Produce a time plot of the series.

```
<- gafa_stock |>
fb_stock filter(Symbol == "FB")
|>
fb_stock autoplot(Close)
```

An upward trend is evident until mid-2018, after which the closing stock price drops.

- Produce forecasts using the drift method and plot them.

The data must be made regular before it can be modelled. We will use trading days as our regular index.

```
<- fb_stock |>
fb_stock mutate(trading_day = row_number()) |>
update_tsibble(index = trading_day, regular = TRUE)
```

Time to model a random walk with drift.

```
|>
fb_stock model(RW(Close ~ drift())) |>
forecast(h = 100) |>
autoplot(fb_stock)
```

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

Prove drift methods are extrapolations from the first and last observation. First, we will demonstrate it graphically.

```
|>
fb_stock model(RW(Close ~ drift())) |>
forecast(h = 100) |>
autoplot(fb_stock) +
geom_line(
aes(y = Close),
linetype = "dashed", colour = "blue",
data = fb_stock |> filter(trading_day %in% range(trading_day))
)
```

To prove it algebraically, note that \begin{align*} \hat{y}_{T+h|T} = y_T + h\left(\frac{y_T-y_1}{T-1}\right) \end{align*} which is a straight line with slope (y_T-y_1)/(T-1) that goes through the point (T,y_T).

Therefore, it must also go through the point (1,c) where (y_T-c)/(T-1) = (y_T - y_1) / (T-1), so c=y_1.

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

Use other appropriate benchmark methods. The most appropriate benchmark method is the naive model. The mean forecast is terrible for this type of data, and the data is non-seasonal.

```
|>
fb_stock model(NAIVE(Close)) |>
forecast(h = 100) |>
autoplot(fb_stock)
```

The naive method is most appropriate, and will also be best if the efficient market hypothesis holds true.

# fpp3 5.10, Ex 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.

```
# Extract data of interest
<- aus_production |>
recent_production filter(year(Quarter) >= 1992)
# Define and estimate a model
<- recent_production |> model(SNAIVE(Beer))
fit # Look at the residuals
|> gg_tsresiduals() fit
```

- The residuals are not centred around 0 (typically being slightly below it), this is due to the model failing to capture the negative trend in the data.
- Peaks and troughs in residuals spaced roughly 4 observations apart are apparent leading to a negative spike at lag 4 in the ACF. So they do not resemble white noise. Lags 1 and 3 are also significant, however they are very close to the threshold and are of little concern.
- The distribution of the residuals does not appear very normal, however it is probably close enough for the accuracy of our intervals (it being not centred on 0 is more concerning).

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

The forecasts look reasonable, although the intervals may be a bit wide. This is likely due to the slight trend not captured by the model (which subsequently violates the assumptions imposed on the residuals).

# fpp3 5.10, Ex 4

Repeat the exercise for 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

The data does not contain seasonality, so the naive model is more appropriate.

```
# Extract data of interest
<- filter(global_economy, Country == "Australia")
aus_exports # Define and estimate a model
<- aus_exports |> model(NAIVE(Exports))
fit # Check residuals
|> gg_tsresiduals() fit
```

- The ACF plot reveals that the first lag exceeds the significance threshold.
- This data may still be white noise, as it is the only lag that exceeds the blue dashed lines (5% of the lines are expected to exceed it). However as it is the first lag, it is probable that there exists some real auto-correlation in the residuals that can be modelled.
- The distribution appears normal.
- The residual plot appears mostly random, however more observations appear to be above zero. This again, is due to the model not capturing the trend.

```
# Look at some forecasts
|>
fit forecast() |>
autoplot(aus_exports)
```

- The forecasts appear reasonable as the series appears to have flattened in recent years.
- The intervals are also reasonable — despite the assumptions behind them having been violated.

## Australian brick production

The data is seasonal, so the seasonal naive model is more appropriate.

```
# Remove the missing values at the end of the series
<- aus_production |>
tidy_bricks filter(!is.na(Bricks))
# Define and estimate a model
<- tidy_bricks |>
fit model(SNAIVE(Bricks))
# Look at the residuals
|> gg_tsresiduals() fit
```

The residual plot does not appear random. Periods of low production and high production are evident, leading to autocorrelation in the residuals.

The residuals from this model are not white noise. The ACF plot shows a strong sinusoidal pattern of decay, indicating that the residuals are auto-correlated.

The histogram is also not normally distributed, as it has a long left tail.

```
# Look at some forecasts
|>
fit forecast() |>
autoplot(tidy_bricks)
```

- The point forecasts appear reasonable as the series appears to have flattened in recent years.
- The intervals appear much larger than necessary.

# fpp3 5.10, Ex 5

Produce forecasts for the 7 Victorian series in

`aus_livestock`

using`SNAIVE()`

. Plot the resulting forecasts including the historical data. Is this a reasonable benchmark for these series?

```
|>
aus_livestock filter(State == "Victoria") |>
model(SNAIVE(Count)) |>
forecast(h = "5 years") |>
autoplot(aus_livestock)
```

- Most point forecasts look reasonable from the seasonal naive method.
- Some series are more seasonal than others, and for the series with very weak seasonality it may be better to consider using a naive or drift method.
- The prediction intervals in some cases go below zero, so perhaps a log transformation would have been better for these series.

# fpp3 5.10, Ex 6

Are the following statements true or false? Explain your answer.

- Good forecast methods should have normally distributed residuals.

- False.
- Although many good forecasting methods produce normally distributed residuals this is not required to produce good forecasts.
- Other forecasting methods may use other distributions, it is just less common as they can be more difficult to work with.

- A model with small residuals will give good forecasts.

- False.
- It is possible to produce a model with small residuals by making a highly complicated (overfitted) model that fits the data extremely well.
- This highly complicated model will often perform very poorly when forecasting new data.

- The best measure of forecast accuracy is MAPE.

- False.
- There is no single best measure of accuracy - often you would want to see a collection of accuracy measures as they can reveal different things about your residuals. MAPE in particular has some substantial disadvantages - extreme values can result when y_t is close to zero, and it assumes that the unit being measured has a meaningful zero.

- If your model doesn’t forecast well, you should make it more complicated.

- False.
- There are many reasons why a model may not forecast well, and making the model more complicated can make the forecasts worse.
- The model specified should capture structures that are evident in the data. Although adding terms that are unrelated to the structures found in the data will improve the model’s residuals, the forecasting performance of the model will not improve.
- Adding missing features relevant to the data (such as including a seasonal pattern that exists in the data) should improve forecast performance.

- Always choose the model with the best forecast accuracy as measured on the test set.

- False.
- There are many measures of forecast accuracy, and the appropriate model is the one which is best suited to the forecasting task. For instance, you may be interested in choosing a model which forecasts well for predictions exactly one year ahead. In this case, using cross-validated accuracy could be a more useful approach to evaluating accuracy.

# fpp3 5.10, Ex 8

Consider the number of pigs slaughtered in New South Wales (data set

`aus_livestock`

).

- Produce some plots of the data in order to become familiar with it.

```
<- aus_livestock |>
nsw_pigs filter(State == "New South Wales", Animal == "Pigs")
|>
nsw_pigs autoplot(Count)
```

Data generally follows a downward trend, however there are some periods where the amount of pigs slaughtered changes rapidly.

`|> gg_season(Count, labels = "right") nsw_pigs `

`|> gg_subseries(Count) nsw_pigs `

Some seasonality is apparent, with notable increases in December and decreases during January, February and April.

- Create a training set of 486 observations, withholding a test set of 72 observations (6 years).

`<- nsw_pigs |> slice(1:486) nsw_pigs_train `

- Try using various benchmark methods to forecast the training set and compare the results on the test set. Which method did best?

```
<- nsw_pigs_train |>
fit model(
mean = MEAN(Count),
naive = NAIVE(Count),
snaive = SNAIVE(Count),
drift = RW(Count ~ drift())
)|>
fit forecast(h = "6 years") |>
accuracy(nsw_pigs)
```

```
# A tibble: 4 × 12
.model Animal State .type ME RMSE MAE MPE MAPE MASE RMSSE
<chr> <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 drift Pigs New South … Test -4685. 8091. 6967. -7.36 10.1 0.657 0.557
2 mean Pigs New South … Test -39360. 39894. 39360. -55.9 55.9 3.71 2.75
3 naive Pigs New South … Test -6138. 8941. 7840. -9.39 11.4 0.740 0.615
4 snaive Pigs New South … Test -5838. 10111. 8174. -8.81 11.9 0.771 0.696
# ℹ 1 more variable: ACF1 <dbl>
```

The drift method performed best for all measures of accuracy (although it had a larger first order auto-correlation)

- Check the residuals of your preferred method. Do they resemble white noise?

```
|>
fit select(drift) |>
gg_tsresiduals()
```

The residuals do not appear to be white noise as the ACF plot contains many significant lags. It is also clear that the seasonal component is not captured by the drift method, as there exists a strong positive auto-correlation at lag 12 (1 year). The histogram appears to have a slightly long left tail.

# fpp3 5.10, Ex 11

We will use the bricks data from

`aus_production`

(Australian quarterly clay brick production 1956–2005) for this exercise.

- Use an STL decomposition to calculate the trend-cycle and seasonal indices. (Experiment with having fixed or changing seasonality.)

```
<- aus_production |>
tidy_bricks filter(!is.na(Bricks))
|>
tidy_bricks model(STL(Bricks)) |>
components() |>
autoplot()
```

Data is multiplicative, and so a transformation should be used.

```
<- tidy_bricks |>
dcmp model(STL(log(Bricks))) |>
components()
|>
dcmp autoplot()
```

Seasonality varies slightly.

```
<- tidy_bricks |>
dcmp model(stl = STL(log(Bricks) ~ season(window = "periodic"))) |>
components()
|> autoplot() dcmp
```

The seasonality looks fairly stable, so I’ve used a periodic season (window). The decomposition still performs well when the seasonal component is fixed. The remainder term does not appear to contain a substantial amount of seasonality.

- Compute and plot the seasonally adjusted data.

```
|>
dcmp as_tsibble() |>
autoplot(season_adjust)
```

- Use a naïve method to produce forecasts of the seasonally adjusted data.

```
<- dcmp |>
fit select(-.model) |>
model(naive = NAIVE(season_adjust)) |>
forecast(h = "5 years")
|>
dcmp as_tsibble() |>
autoplot(season_adjust) + autolayer(fit)
```

- Use
`decomposition_model()`

to reseasonalise the results, giving forecasts for the original data.

```
<- tidy_bricks |>
fit model(stl_mdl = decomposition_model(STL(log(Bricks)), NAIVE(season_adjust)))
|>
fit forecast(h = "5 years") |>
autoplot(tidy_bricks)
```

- Do the residuals look uncorrelated?

`|> gg_tsresiduals() fit `

The residuals do not appear uncorrelated as there are several lags of the ACF which exceed the significance threshold.

- Repeat with a robust STL decomposition. Does it make much difference?

```
<- tidy_bricks |>
fit_robust model(stl_mdl = decomposition_model(STL(log(Bricks)), NAIVE(season_adjust)))
|> gg_tsresiduals() fit_robust
```

The residuals appear slightly less auto-correlated, however there is still significant auto-correlation at lag 8.

- Compare forecasts from
`decomposition_model()`

with those from`SNAIVE()`

, using a test set comprising the last 2 years of data. Which is better?

```
<- tidy_bricks |>
tidy_bricks_train slice(1:(n() - 8))
<- tidy_bricks_train |>
fit model(
stl_mdl = decomposition_model(STL(log(Bricks)), NAIVE(season_adjust)),
snaive = SNAIVE(Bricks)
)
<- fit |>
fc forecast(h = "2 years")
|>
fc autoplot(tidy_bricks, level = NULL)
```

The decomposition forecasts appear to more closely follow the actual future data.

```
|>
fc accuracy(tidy_bricks)
```

```
# 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 snaive Test 2.75 20 18.2 0.395 4.52 0.504 0.407 -0.0503
2 stl_mdl Test 0.368 18.1 15.1 -0.0679 3.76 0.418 0.368 0.115
```

The STL decomposition forecasts are more accurate than the seasonal naive forecasts across all accuracy measures.

# fpp3 5.10, Ex 12

`tourism`

contains quarterly visitor nights (in thousands) from 1998 to 2017 for 76 regions of Australia.

- Extract data from the Gold Coast region using
`filter()`

and aggregate total overnight trips (sum over`Purpose`

) using`summarise()`

. Call this new dataset`gc_tourism`

.

```
<- tourism |>
gc_tourism filter(Region == "Gold Coast") |>
summarise(Trips = sum(Trips))
gc_tourism
```

```
# A tsibble: 80 x 2 [1Q]
Quarter Trips
<qtr> <dbl>
1 1998 Q1 827.
2 1998 Q2 681.
3 1998 Q3 839.
4 1998 Q4 820.
5 1999 Q1 987.
6 1999 Q2 751.
7 1999 Q3 822.
8 1999 Q4 914.
9 2000 Q1 871.
10 2000 Q2 780.
# ℹ 70 more rows
```

- Using
`slice()`

or`filter()`

, create three training sets for this data excluding the last 1, 2 and 3 years. For example,`gc_train_1 <- gc_tourism |> slice(1:(n()-4))`

.

```
<- gc_tourism |> slice(1:(n() - 4))
gc_train_1 <- gc_tourism |> slice(1:(n() - 8))
gc_train_2 <- gc_tourism |> slice(1:(n() - 12)) gc_train_3
```

- Compute one year of forecasts for each training set using the seasonal naïve (
`SNAIVE()`

) method. Call these`gc_fc_1`

,`gc_fc_2`

and`gc_fc_3`

, respectively.

```
<- bind_cols(
gc_fc |> model(gc_fc_1 = SNAIVE(Trips)),
gc_train_1 |> model(gc_fc_2 = SNAIVE(Trips)),
gc_train_2 |> model(gc_fc_3 = SNAIVE(Trips))
gc_train_3 |> forecast(h = "1 year") )
```

`|> autoplot(gc_tourism) gc_fc `

- Use
`accuracy()`

to compare the test set forecast accuracy using MAPE. Comment on these.

`|> accuracy(gc_tourism) gc_fc `

```
# 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 gc_fc_1 Test 75.1 167. 154. 6.36 15.1 2.66 2.36 -0.410
2 gc_fc_2 Test 12.0 43.1 39.5 1.14 4.32 0.670 0.599 -0.792
3 gc_fc_3 Test 35.8 91.4 83.9 3.56 9.07 1.46 1.30 0.239
```

The second set of forecasts are most accurate (as can be seen in the previous plot), however this is likely due to chance.