Exercise solutions: Section 9.11

Author

Rob J Hyndman and George Athanasopoulos

fpp3 9.11, Ex7

Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.

  1. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
aus_airpassengers |> autoplot(Passengers)

fit <- aus_airpassengers |>
  model(arima = ARIMA(Passengers))
report(fit)
Series: Passengers 
Model: ARIMA(0,2,1) 

Coefficients:
          ma1
      -0.8963
s.e.   0.0594

sigma^2 estimated as 4.308:  log likelihood=-97.02
AIC=198.04   AICc=198.32   BIC=201.65
fit |> gg_tsresiduals()

fit |> forecast(h = 10) |> autoplot(aus_airpassengers)

  1. Write the model in terms of the backshift operator.
    year 
4.307764

\[(1-B)^2y_t = (1+\theta B)\varepsilon_t\] where \(\varepsilon\sim\text{N}(0,\sigma^2)\), \(\theta = -0.90\) and \(\sigma^2 = 4.31\).

  1. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 1 + pdq(0,1,0))) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers)

  • Both containing increasing trends, but the ARIMA(0,2,1) model has an implicit trend due to the double-differencing, while the ARIMA(0,1,0) with drift models the trend directly via the trend coefficient.
  • The intervals are narrower when there are fewer differences.
  1. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to part b. Remove the constant and see what happens.
aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 1 + pdq(2,1,2))) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers)

aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 0 + pdq(2,1,2)))
# A mable: 1 x 1
         arima
       <model>
1 <NULL model>
  • There is little difference between ARIMA(2,1,2) with drift and ARIMA(0,1,0) with drift.
  • Removing the constant causes an error because the model cannot be estimated.
  1. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 1 + pdq(0,2,1))) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers)
Warning: Model specification induces a quadratic or higher order polynomial trend. 
This is generally discouraged, consider removing the constant or reducing the number of differences.

The forecast trend is now quadratic, and there is a warning that this is generally a bad idea.

fpp3 9.11, Ex8

For the United States GDP series (from global_economy):

  1. If necessary, find a suitable Box-Cox transformation for the data;
us_economy <- global_economy |>
  filter(Code == "USA")
us_economy |>
  autoplot(GDP)

lambda <- us_economy |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)
lambda
[1] 0.2819443
us_economy |>
  autoplot(box_cox(GDP, lambda))

It seems that a Box-Cox transformation may be useful here.

  1. fit a suitable ARIMA model to the transformed data using ARIMA();
fit <- us_economy |>
  model(ARIMA(box_cox(GDP, lambda)))
report(fit)
Series: GDP 
Model: ARIMA(1,1,0) w/ drift 
Transformation: box_cox(GDP, lambda) 

Coefficients:
         ar1  constant
      0.4586  118.1822
s.e.  0.1198    9.5047

sigma^2 estimated as 5479:  log likelihood=-325.32
AIC=656.65   AICc=657.1   BIC=662.78
  1. try some other plausible models by experimenting with the orders chosen;
fit <- us_economy |>
  model(
    arima010 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 0)),
    arima011 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 1)),
    arima012 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 2)),
    arima013 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 3)),
    arima110 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 0)),
    arima111 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 1)),
    arima112 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 2)),
    arima113 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 3)),
    arima210 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 0)),
    arima211 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 1)),
    arima212 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 2)),
    arima213 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 3)),
    arima310 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 0)),
    arima311 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 1)),
    arima312 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 2)),
    arima313 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 3))
  )
  1. choose what you think is the best model and check the residual diagnostics;
fit |>
  glance() |>
  arrange(AICc) |>
  select(.model, AICc)
# A tibble: 16 × 2
   .model    AICc
   <chr>    <dbl>
 1 arima110  657.
 2 arima011  659.
 3 arima111  659.
 4 arima210  659.
 5 arima012  660.
 6 arima112  661.
 7 arima211  661.
 8 arima310  662.
 9 arima013  662.
10 arima312  663.
11 arima311  664.
12 arima113  664.
13 arima212  664.
14 arima313  665.
15 arima213  666.
16 arima010  668.

The best according to the AICc values is the ARIMA(1,1,0) w/ drift model.

best_fit <- us_economy |>
  model(ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 0)))
best_fit |> report()
Series: GDP 
Model: ARIMA(1,1,0) w/ drift 
Transformation: box_cox(GDP, lambda) 

Coefficients:
         ar1  constant
      0.4586  118.1822
s.e.  0.1198    9.5047

sigma^2 estimated as 5479:  log likelihood=-325.32
AIC=656.65   AICc=657.1   BIC=662.78
best_fit |> gg_tsresiduals()

augment(best_fit) |> features(.innov, ljung_box, dof = 2, lag = 10)
# A tibble: 1 × 4
  Country       .model                                         lb_stat lb_pvalue
  <fct>         <chr>                                            <dbl>     <dbl>
1 United States ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 0))    3.81     0.874

The residuals pass the Ljung-Box test, but the histogram looks like negatively skewed.

  1. produce forecasts of your fitted model. Do the forecasts look reasonable?
best_fit |>
  forecast(h = 10) |>
  autoplot(us_economy)

These look reasonable. Let’s compare a model with no transformation.

fit1 <- us_economy |> model(ARIMA(GDP))
fit1 |>
  forecast(h = 10) |>
  autoplot(us_economy)

Notice the effect of the transformation on the forecasts. Increase the forecast horizon to see what happens. Notice also the width of the prediction intervals.

us_economy |>
  model(
    ARIMA(GDP),
    ARIMA(box_cox(GDP, lambda))
  ) |>
  forecast(h = 20) |>
  autoplot(us_economy)

  1. compare the results with what you would obtain using ETS() (with no transformation).
us_economy |>
  model(ETS(GDP)) |>
  forecast(h = 10) |>
  autoplot(us_economy)

The point forecasts are similar, however the ETS forecast intervals are much wider.

fpp3 9.11, Ex15

Consider the number of Snowshoe Hare furs traded by the Hudson Bay Company between 1845 and 1935 (data set pelt).

  1. Produce a time plot of the time series.
pelt |>
  autoplot(Hare)

  1. Assume you decide to fit the following model: \[ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \phi_4 y_{t-4} + \varepsilon_t, \] where \(\varepsilon_t\) is a white noise series. What sort of ARIMA model is this (i.e., what are \(p\), \(d\), and \(q\))?
  • This is an ARIMA(4,0,0), hence \(p=4\), \(d=0\) and \(q=0\).
  1. By examining the ACF and PACF of the data, explain why this model is appropriate.
pelt |> gg_tsdisplay(Hare, plot="partial")

fit <- pelt |> model(AR4 = ARIMA(Hare ~ pdq(4,0,0)))
fit |> gg_tsresiduals()

  • The significant spike at lag 4 of the PACF indicates an AR(4).
  • The residuals from this model are clearly whhite noise.

fpp3 9.11, Ex16

The population of Switzerland from 1960 to 2017 is in data set global_economy.

  1. Produce a time plot of the data.
swiss_pop <- global_economy |>
  filter(Country == "Switzerland") |>
  select(Year, Population) |>
  mutate(Population = Population / 1e6)

autoplot(swiss_pop, Population)

  1. You decide to fit the following model to the series: \[y_t = c + y_{t-1} + \phi_1 (y_{t-1} - y_{t-2}) + \phi_2 (y_{t-2} - y_{t-3}) + \phi_3( y_{t-3} - y_{t-4}) + \varepsilon_t\] where \(y_t\) is the Population in year \(t\) and \(\varepsilon_t\) is a white noise series. What sort of ARIMA model is this (i.e., what are \(p\), \(d\), and \(q\))?

This is an ARIMA(3,1,0), hence \(p=3\), \(d=1\) and \(q=0\).

  1. Explain why this model was chosen using the ACF and PACF of the differenced series.
swiss_pop |> gg_tsdisplay(Population, plot="partial")

swiss_pop |> gg_tsdisplay(difference(Population), plot="partial")

The significant spike at lag 3 in the PACF, coupled with the exponential decay in the ACF, for the differenced series, signals an AR(3) for the differenced series.