Time series regression, autocorrelation, moving averages.

Section 5.11 Problem #8

8. Consider the number of pigs slaughtered in New South Wales (data set aus_livestock).

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

nsw_pigs <- aus_livestock |>
  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.

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

nsw_pigs |> gg_subseries(Count)

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

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

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

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

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

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

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

Section 7.10 Problem #5

5. The us_gasoline series consists of weekly data for supplies of US finished motor gasoline product, from 2 February 1991 to 20 January 2017. The units are in “million barrels per day”. Consider only the data to the end of 2004.

a.Fit a harmonic regression with trend to the data. Experiment with changing the number Fourier terms. Plot the observed gasoline and fitted values and comment on what you see.

gas <- us_gasoline |> filter(year(Week) <= 2004)
gas |> autoplot(Barrels)

fit <- gas |>
  model(
    fourier1 = TSLM(Barrels ~ trend() + fourier(K = 1)),
    fourier2 = TSLM(Barrels ~ trend() + fourier(K = 2)),
    fourier3 = TSLM(Barrels ~ trend() + fourier(K = 3)),
    fourier4 = TSLM(Barrels ~ trend() + fourier(K = 4)),
    fourier5 = TSLM(Barrels ~ trend() + fourier(K = 5)),
    fourier6 = TSLM(Barrels ~ trend() + fourier(K = 6)),
    fourier7 = TSLM(Barrels ~ trend() + fourier(K = 7)),
    fourier8 = TSLM(Barrels ~ trend() + fourier(K = 8)),
    fourier9 = TSLM(Barrels ~ trend() + fourier(K = 9)),
    fourier10 = TSLM(Barrels ~ trend() + fourier(K = 10)),
    fourier11 = TSLM(Barrels ~ trend() + fourier(K = 11)),
    fourier12 = TSLM(Barrels ~ trend() + fourier(K = 12)),
    fourier13 = TSLM(Barrels ~ trend() + fourier(K = 13)),
    fourier14 = TSLM(Barrels ~ trend() + fourier(K = 14)),
    fourier15 = TSLM(Barrels ~ trend() + fourier(K = 15))
  )

glance(fit) |>
  arrange(AICc) |>
  select(.model, AICc, CV)

Best model has order 7

gas |>
  autoplot(Barrels, col = "gray") +
  geom_line(
    data = augment(fit) |> filter(.model == "fourier7"),
    aes(y = .fitted), col = "blue"
  )

This seems to have captured the annual seasonality quite well.

b. Select the appropriate number of Fourier terms to include by minimising the AICc or CV value.

See above

c. Plot the residuals of the final model using the gg_tsresiduals() function and comment on these. Use a Ljung-Box test to check for residual autocorrelation.

fit |>
  select(fourier7) |>
  gg_tsresiduals()

There is some small remaining autocorrelation at lag 1.

augment(fit) |>
  filter(.model == "fourier7") |>
  features(.innov, ljung_box, lag = 52)

The Ljung-Box test is not significant. Even if the Ljung-Box test had given a significant result here, the correlations are very small, so it would not make much difference to the forecasts and prediction intervals.

d. Generate forecasts for the next year of data and plot these along with the actual data for 2005. Comment on the forecasts.

fit |>
  select("fourier7") |>
  forecast(h = "1 year") |>
  autoplot(filter_index(us_gasoline, "2000" ~ "2006"))

The forecasts fit pretty well for the first six months, but not so well after that.