All exercises are from Chapter 9 of Hyndman and Athanasopoulos, Forecasting: Principles and Practice (3rd ed.), on ARIMA models.

Exercise 9.1

Figure 9.32 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

(a) Explain the differences among these figures. Do they all indicate that the data are white noise?

In the book’s Figure 9.32, all three ACFs are consistent with white noise: the bars are small, show no pattern or slow decay, and nearly all of them sit inside the dashed lines. The figures differ in two ways. In the left panel (36 numbers) the bars are larger and the dashed lines are far from zero. In the middle (360 numbers) and right (1,000 numbers) panels the bars are tiny and the dashed lines sit close to zero. The panels differ in sample size, and since they are random samples the exact bars differ too. All three ACFs are consistent with white noise.

(b) Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

The dashed lines in an ACF plot are placed at plus or minus 1.96 divided by the square root of the series length T. The table computes them for the three series lengths in the figure.

# the dashed lines are at plus/minus 1.96 / sqrt(T) for a series of length T
tibble(T = c(36, 360, 1000)) |>
  mutate(critical_value = 1.96 / sqrt(T)) |>
  knitr::kable(digits = 3)
T critical_value
36 0.327
360 0.103
1000 0.062

The lines are about 0.327 from zero for 36 numbers, 0.103 for 360 and 0.062 for 1,000, so they move toward zero as T grows. With more data an autocorrelation is estimated more precisely, so a smaller bar is already distinguishable from zero.

The autocorrelations differ between the panels because each panel is a different random sample. White noise has zero population autocorrelation at every nonzero lag, but sample autocorrelations fluctuate around zero. In a sample, the bars scatter around zero by an amount that shrinks as the sample grows, which is why the bars in the 1,000-number panel are the smallest. The dashed lines are drawn where about 95% of the bars of a white noise series should fall, so about 1 bar in 20 can cross them by chance without any pattern in the data.

Exercise 9.2

A classic example of a non-stationary series is stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

# gafa_stock only has trading days, so the calendar has gaps on weekends and holidays;
# indexing by trading day gives a regular series, as in the book's stock price examples
amazon <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  mutate(trading_day = row_number()) |>
  update_tsibble(index = trading_day, regular = TRUE)

amazon |> gg_tsdisplay(Close, plot_type = "partial")

The closing price rises from the start of the series to a peak near the end and then falls, so it has no constant mean and no stable level, which already points to a non-stationary series. The ACF supports this: the bars start close to 1 and fall very slowly, staying far above the dashed lines at every lag shown, which is the pattern a series with a trend produces. The PACF has one very large bar at lag 1, close to 1, and only small bars after that, so the largest partial autocorrelation is at lag 1. Together these plots are consistent with a highly persistent, non-stationary series, and the standard remedy is to difference it, which models the day-to-day change instead of the price level. The x-axis is the trading day number, not the calendar date, because of the re-indexing in the code.

Exercise 9.3

For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

(a) Turkish GDP from global_economy.

turkey <- global_economy |> filter(Country == "Turkey")
turkey |> autoplot(GDP)

# Guerrero's method picks lambda; a value near 0 would mean a log
lambda_turkey <- turkey |> features(GDP, guerrero) |> pull(lambda_guerrero)
lambda_turkey
## [1] 0.1572187
turkey |> autoplot(box_cox(GDP, lambda_turkey))

# annual data, so no seasonal difference: ask how many first differences KPSS needs
turkey |> features(box_cox(GDP, lambda_turkey), unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1
turkey |>
  gg_tsdisplay(difference(box_cox(GDP, lambda_turkey)), plot_type = "partial")
## 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()`).

turkey |>
  mutate(diff_gdp = difference(box_cox(GDP, lambda_turkey))) |>
  features(diff_gdp, unitroot_kpss)
## # A tibble: 1 × 3
##   Country kpss_stat kpss_pvalue
##   <fct>       <dbl>       <dbl>
## 1 Turkey     0.0889         0.1

Turkish GDP has an overall upward trend with some years of decline, and the year-to-year swings get larger in absolute size later in the series. Guerrero’s method gives a lambda of about 0.16, between a log and a square root, and the transformed series grows more steadily than the raw one. It still has a clear trend, so it is not stationary. unitroot_ndiffs() asks for one first difference. After taking it, the series moves around a fairly constant level, the ACF and PACF show no strong bars (only lag 4 comes close to the dashed line), and KPSS reports a p-value of 0.1 (the book reports anything above 0.1 as 0.1), so the test does not reject stationarity. The annual data have no seasonal pattern, so the answer is a Box-Cox transformation with lambda 0.16 and one first difference.

(b) Accommodation takings in the state of Tasmania from aus_accommodation.

tasmania <- aus_accommodation |> filter(State == "Tasmania")
tasmania |> autoplot(Takings)

lambda_tas <- tasmania |> features(Takings, guerrero) |> pull(lambda_guerrero)
lambda_tas
## [1] 0.001819643
tasmania |> autoplot(box_cox(Takings, lambda_tas))

# seasonal strength from the STL features: above 0.64 the book suggests one seasonal difference
tasmania |> features(box_cox(Takings, lambda_tas), feat_stl) |> select(State, seasonal_strength_year)
## # A tibble: 1 × 2
##   State    seasonal_strength_year
##   <chr>                     <dbl>
## 1 Tasmania                  0.992
# quarterly data: ask for seasonal differences first, then first differences on top
tasmania |> features(box_cox(Takings, lambda_tas), unitroot_nsdiffs)
## # A tibble: 1 × 2
##   State    nsdiffs
##   <chr>      <int>
## 1 Tasmania       1
tasmania |>
  mutate(sdiff_takings = difference(box_cox(Takings, lambda_tas), 4)) |>
  features(sdiff_takings, unitroot_ndiffs)
## # A tibble: 1 × 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      0
tasmania |>
  gg_tsdisplay(difference(box_cox(Takings, lambda_tas), 4), plot_type = "partial")
## 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()`).

tasmania |>
  mutate(sdiff_takings = difference(box_cox(Takings, lambda_tas), 4)) |>
  features(sdiff_takings, unitroot_kpss)
## # A tibble: 1 × 3
##   State    kpss_stat kpss_pvalue
##   <chr>        <dbl>       <dbl>
## 1 Tasmania     0.198         0.1

Tasmanian takings are quarterly with a strong yearly pattern whose swings grow with the level, and Guerrero’s lambda is almost 0, so this is effectively a log transformation. After it, the swings are about the same size throughout. The seasonal strength is 0.992, well above the 0.64 cutoff, and unitroot_nsdiffs() asks for one seasonal difference (lag 4). After taking it, unitroot_ndiffs() asks for no further first difference. KPSS reports a p-value of 0.1 for the seasonally differenced series (anything above 0.1 is reported as 0.1), so the test does not reject stationarity. The ACF is large over the first three lags and then drops inside the dashed lines, so I stop with only the seasonal difference. The answer is a log-type transformation (lambda close to 0) and one seasonal difference, with no first difference.

(c) Monthly sales from souvenirs.

souvenirs |> autoplot(Sales)

lambda_sou <- souvenirs |> features(Sales, guerrero) |> pull(lambda_guerrero)
lambda_sou
## [1] 0.002118221
souvenirs |> autoplot(box_cox(Sales, lambda_sou))

# seasonal strength from the STL features: above 0.64 the book suggests one seasonal difference
souvenirs |> features(box_cox(Sales, lambda_sou), feat_stl) |> select(seasonal_strength_year)
## # A tibble: 1 × 1
##   seasonal_strength_year
##                    <dbl>
## 1                  0.940
souvenirs |> features(box_cox(Sales, lambda_sou), unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
souvenirs |>
  mutate(sdiff_sales = difference(box_cox(Sales, lambda_sou), 12)) |>
  features(sdiff_sales, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0
souvenirs |>
  gg_tsdisplay(difference(box_cox(Sales, lambda_sou), 12), plot_type = "partial", lag_max = 36)
## 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()`).

souvenirs |>
  mutate(sdiff_sales = difference(box_cox(Sales, lambda_sou), 12)) |>
  features(sdiff_sales, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.240         0.1

Souvenir sales have a large spike every December, and the spikes generally get larger over the years, so the variance rises with the level. Guerrero’s lambda is again close to 0, so a log transformation is appropriate, and it makes the December spikes more similar in size. The seasonal strength is 0.940, above the 0.64 cutoff, and unitroot_nsdiffs() asks for one seasonal difference (lag 12), and unitroot_ndiffs() then asks for no first difference. KPSS reports a p-value of 0.1 after the seasonal difference, which does not reject stationarity. The ACF of the differenced series is not perfectly clean: it stays positive for the first ten lags and goes negative later, so some autocorrelation is left. A stationary series can still be autocorrelated, so that alone is not a reason to difference again. I still choose one seasonal difference only, because KPSS does not reject stationarity after the seasonal difference, unitroot_ndiffs() asks for no first difference, and, as the book warns in Section 9.1, unnecessary differences can add autocorrelation that is not in the data. That leaves a log-type transformation and one seasonal difference.

Exercise 9.5

For your retail data (from Exercise 7 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

set.seed(219562)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
myseries |> distinct(State, Industry)
## # A tibble: 1 × 2
##   State      Industry                                     
##   <chr>      <chr>                                        
## 1 Queensland Cafes, restaurants and takeaway food services
myseries |> autoplot(Turnover)

lambda_retail <- myseries |> features(Turnover, guerrero) |> pull(lambda_guerrero)
lambda_retail
## [1] 0.1422668
myseries |> autoplot(box_cox(Turnover, lambda_retail))

# seasonal strength from the STL features: above 0.64 the book suggests one seasonal difference
myseries |> features(box_cox(Turnover, lambda_retail), feat_stl) |> select(State, Industry, seasonal_strength_year)
## # A tibble: 1 × 3
##   State      Industry                                     seasonal_strength_year
##   <chr>      <chr>                                                         <dbl>
## 1 Queensland Cafes, restaurants and takeaway food servic…                  0.753
# seasonal difference first, then check whether a first difference is still needed
myseries |> features(box_cox(Turnover, lambda_retail), unitroot_nsdiffs)
## # A tibble: 1 × 3
##   State      Industry                                      nsdiffs
##   <chr>      <chr>                                           <int>
## 1 Queensland Cafes, restaurants and takeaway food services       1
myseries |>
  mutate(sdiff = difference(box_cox(Turnover, lambda_retail), 12)) |>
  features(sdiff, unitroot_ndiffs)
## # A tibble: 1 × 3
##   State      Industry                                      ndiffs
##   <chr>      <chr>                                          <int>
## 1 Queensland Cafes, restaurants and takeaway food services      1
myseries |>
  gg_tsdisplay(difference(box_cox(Turnover, lambda_retail), 12), plot_type = "partial", 
               lag_max = 36)
## 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()`).

myseries |>
  gg_tsdisplay(difference(difference(box_cox(Turnover, lambda_retail), 12)),
               plot_type = "partial", lag_max = 36)
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

myseries |>
  mutate(ddiff = difference(difference(box_cox(Turnover, lambda_retail), 12))) |>
  features(ddiff, unitroot_kpss)
## # A tibble: 1 × 4
##   State      Industry                                      kpss_stat kpss_pvalue
##   <chr>      <chr>                                             <dbl>       <dbl>
## 1 Queensland Cafes, restaurants and takeaway food services    0.0132         0.1

My series is Queensland cafes, restaurants and takeaway food turnover. It has an upward trend, a yearly seasonal pattern and variance that grows with the level, so I transformed it first with Guerrero’s lambda (about 0.14). The seasonal strength is 0.753, above the 0.64 cutoff, and unitroot_nsdiffs() asks for one seasonal difference. After that, the plot still drifts and the ACF decays slowly, with the ACF starting near 0.8 at lag 1 and still above the dashed lines at lag 9, and unitroot_ndiffs() asks for one more first difference. After both differences, the series is centred on zero and the spread is more stable overall, with a few large deviations (the biggest early in the series), the ACF no longer decays slowly, and KPSS reports a p-value of 0.1, which does not reject stationarity. The answer for this series is a Box-Cox transformation with lambda 0.14, one seasonal difference and one first difference. The large negative bar at lag 12 in the final ACF suggests a seasonal moving average term for modelling later.

Exercise 9.6

Simulate and plot some data from simple ARIMA models.

(a) Use the following R code to generate data from an AR(1) model with phi_1 = 0.6 and sigma^2 = 1. The process starts with y_1 = 0.

set.seed(624)
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

(b) Produce a time plot for the series. How does the plot change as you change phi_1?

sim |> autoplot(y)

# same errors e each time, so only phi_1 differs between the panels
ar1 <- function(phi) {
  y <- numeric(100)
  for(i in 2:100)
    y[i] <- phi*y[i-1] + e[i]
  tibble(idx = seq_len(100), y = y, phi = paste("phi_1 =", phi))
}
bind_rows(ar1(0), ar1(0.6), ar1(0.95), ar1(-0.6)) |>
  as_tsibble(index = idx, key = phi) |>
  autoplot(y) +
  facet_wrap(vars(phi), ncol = 1, scales = "free_y") +
  theme(legend.position = "none")

With phi_1 = 0.6 the series moves up and down more smoothly than pure noise, because each value carries over part of the one before it. Changing phi_1 (same errors in each panel) shows how much carry-over there is. At 0 the series is just the white noise. At 0.6 it drifts in longer swings. At 0.95 it wanders far from zero in slow, long movements and looks like a random walk, which makes sense since a phi_1 of 1 would be exactly a random walk. At -0.6 the series zigzags around zero, because each value tends to partly reverse the one before it, although not every step changes sign. All four values (0, 0.6, 0.95 and -0.6) are inside the range -1 to 1 that keeps an AR(1) stationary, so even the 0.95 series is stationary, only very persistent.

(c) Write your own code to generate data from an MA(1) model with theta_1 = 0.6 and sigma^2 = 1.

# MA(1): each value is the current error plus 0.6 times the previous error
y_ma <- numeric(100)
for(i in 2:100)
  y_ma[i] <- e[i] + 0.6*e[i-1]
sim_ma <- tsibble(idx = seq_len(100), y = y_ma, index = idx)

(d) Produce a time plot for the series. How does the plot change as you change theta_1?

sim_ma |> autoplot(y)

ma1 <- function(theta) {
  y <- numeric(100)
  for(i in 2:100)
    y[i] <- e[i] + theta*e[i-1]
  tibble(idx = seq_len(100), y = y, theta = paste("theta_1 =", theta))
}
bind_rows(ma1(0), ma1(0.6), ma1(0.95), ma1(-0.6)) |>
  as_tsibble(index = idx, key = theta) |>
  autoplot(y) +
  facet_wrap(vars(theta), ncol = 1, scales = "free_y") +
  theme(legend.position = "none")

With theta_1 = 0.6 the series looks smoother than white noise, because each value includes part of the previous error, so neighbouring values share some of the same shock. Compared with the AR(1) plots, the effect lasts only one period. Changing theta_1 (same errors in each panel) works the same way: at 0 the series is white noise, at 0.6 and 0.95 neighbouring values are similar and the series is smoother, and at -0.6 consecutive values tend to move in opposite directions, which makes the line more jagged. Unlike the AR(1) with phi_1 = 0.95, the MA(1) with theta_1 = 0.95 stays near zero and does not wander, because a shock affects only two values.

(e) Generate data from an ARMA(1,1) model with phi_1 = 0.6, theta_1 = 0.6 and sigma^2 = 1.

y_arma <- numeric(100)
for(i in 2:100)
  y_arma[i] <- 0.6*y_arma[i-1] + e[i] + 0.6*e[i-1]
sim_arma <- tsibble(idx = seq_len(100), y = y_arma, index = idx)

(f) Generate data from an AR(2) model with phi_1 = -0.8, phi_2 = 0.3 and sigma^2 = 1. (Note that these parameters will give a non-stationary series.)

y_ar2 <- numeric(100)
for(i in 3:100)
  y_ar2[i] <- -0.8*y_ar2[i-1] + 0.3*y_ar2[i-2] + e[i]
sim_ar2 <- tsibble(idx = seq_len(100), y = y_ar2, index = idx)

(g) Graph the latter two series and compare them.

bind_rows(
  tibble(idx = seq_len(100), y = y_arma, model = "ARMA(1,1)"),
  tibble(idx = seq_len(100), y = y_ar2,  model = "AR(2)")
) |>
  as_tsibble(index = idx, key = model) |>
  autoplot(y) +
  facet_wrap(vars(model), ncol = 1, scales = "free_y") +
  theme(legend.position = "none")

The ARMA(1,1) series is stationary: it moves up and down around zero in smooth swings with a fairly stable spread, much like the AR(1) with phi_1 = 0.6. The AR(2) series with phi_1 = -0.8 and phi_2 = 0.3 is not stationary, and the plot shows it: the values alternate in sign from one period to the next, and the size of the swings keeps growing until they reach the thousands by the end. The parameters break the AR(2) stationarity condition phi_2 - phi_1 < 1, since here phi_2 - phi_1 = 0.3 - (-0.8) = 1.1, so the series is explosive rather than settling around a mean. The two panels use very different y-axis scales for that reason.

Exercise 9.7

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

(a) 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.

The exercise text says the data cover 1970 to 2011, but the aus_airpassengers data I have run through 2016, so I use the full supplied series.

# first and last year in the data
range(aus_airpassengers$Year)
## [1] 1970 2016
aus_airpassengers |> autoplot(Passengers)

fit_air <- aus_airpassengers |> model(auto = ARIMA(Passengers))
report(fit_air)
## 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_air |> gg_tsresiduals()

# K = p + q = 1 parameters, so dof = 1
augment(fit_air) |> features(.innov, ljung_box, lag = 10, dof = 1)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 auto      6.70     0.669
fit_air |> forecast(h = 10) |> autoplot(aus_airpassengers)

The passenger series has an overall upward trend and no seasonality, as expected for annual data. ARIMA() chose an ARIMA(0,2,1) with no constant. This means it differenced twice and then fitted a single moving average term, with ma1 = -0.8963. The residuals look like white noise: the ACF bars are all inside the dashed lines and the Ljung-Box p-value is 0.669. I used 10 lags and set dof = 1 to account for the one estimated MA coefficient, which leaves 9 degrees of freedom for the test. The residual histogram is roughly centred on zero, with a few large residuals, the largest in 2010. The forecasts continue the upward path as a nearly straight line, and the prediction intervals widen as the horizon grows, as the book says to expect when d is greater than 0.

(b) Write the model in terms of the backshift operator.

For an ARIMA(0,2,1) with no constant, the model is (1 - B)^2 y_t = (1 + theta_1 B) epsilon_t. With the ma1 estimate of -0.8963 printed by report() in part (a), that is (1 - B)^2 y_t = (1 - 0.8963 B) epsilon_t. Written out without the backshift operator, it is y_t = 2 y_{t-1} - y_{t-2} + epsilon_t - 0.8963 epsilon_{t-1}.

(c) Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.

fit_air_c <- aus_airpassengers |>
  model(auto = ARIMA(Passengers),
        drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)))
report(fit_air_c |> select(drift))
## Series: Passengers 
## Model: ARIMA(0,1,0) w/ drift 
## 
## Coefficients:
##       constant
##         1.4191
## s.e.    0.3014
## 
## sigma^2 estimated as 4.271:  log likelihood=-98.16
## AIC=200.31   AICc=200.59   BIC=203.97
fit_air_c |> forecast(h = 10) |> autoplot(aus_airpassengers, level = NULL)

The ARIMA(0,1,0) with drift is a random walk with drift: the forecast goes up by the estimated drift, 1.42 million passengers, each year, so it is a straight line from the last observation with that slope (the actual yearly changes also contain random errors). Compared with part (a), its forecast line is less steep, so by the end of the horizon it is lower than the ARIMA(0,2,1) forecast.

(d) Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.

fit_air_d <- aus_airpassengers |>
  model(auto = ARIMA(Passengers),
        drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)),
        arima212 = ARIMA(Passengers ~ 1 + pdq(2, 1, 2)))
report(fit_air_d |> select(arima212))
## Series: Passengers 
## Model: ARIMA(2,1,2) w/ drift 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2  constant
##       -0.5518  -0.7327  0.5895  1.0000    3.2216
## s.e.   0.1684   0.1224  0.0916  0.1008    0.7224
## 
## sigma^2 estimated as 4.031:  log likelihood=-96.23
## AIC=204.46   AICc=206.61   BIC=215.43
fit_air_d |> forecast(h = 10) |> autoplot(aus_airpassengers, level = NULL)

# the two d = 1 models can be compared by AICc; the (0,2,1) from part (a) has d = 2 and cannot
fit_air_d |>
  select(-auto) |>
  glance() |>
  select(.model, AIC, AICc, BIC) |>
  knitr::kable(digits = 2)
.model AIC AICc BIC
drift 200.31 200.59 203.97
arima212 204.46 206.61 215.43

With the constant, the ARIMA(2,1,2) forecast is almost identical to the drift model from part (c); the two lines overlap in the plot, and both are less steep than the ARIMA(0,2,1) from part (a). Both models have d = 1, so their AICc values can be compared, and the AICc table favours the drift model (200.59, against 206.61 for the ARIMA(2,1,2) with a constant).

Now I remove the constant from the ARIMA(2,1,2), as the exercise asks.

# same order with the constant removed (0 in the formula, as in the book's ~ 0 syntax)
fit_air_d2 <- aus_airpassengers |>
  model(arima212_noc = ARIMA(Passengers ~ 0 + pdq(2, 1, 2)))
## Warning: 1 error encountered for arima212_noc
## [1] non-stationary AR part from CSS
fit_air_d2
## # A mable: 1 x 1
##   arima212_noc
##        <model>
## 1 <NULL model>

With the constant removed, the default estimation fails. The output gives the error “non-stationary AR part from CSS” and the model comes back as a NULL model, so there is no fitted model and no forecast to plot or compare.

(e) Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?

fit_air_e <- aus_airpassengers |>
  model(auto = ARIMA(Passengers),
        arima021c = ARIMA(Passengers ~ 1 + pdq(0, 2, 1)))
## 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.
report(fit_air_e |> select(arima021c))
## Series: Passengers 
## Model: ARIMA(0,2,1) w/ poly 
## 
## Coefficients:
##           ma1  constant
##       -1.0000    0.0571
## s.e.   0.0585    0.0213
## 
## sigma^2 estimated as 3.855:  log likelihood=-95.1
## AIC=196.21   AICc=196.79   BIC=201.63
fit_air_e |> forecast(h = 10) |> autoplot(aus_airpassengers, level = NULL)

The ARIMA(0,2,1) with a constant gives the warning that the model “induces a quadratic or higher order polynomial trend”, which the book says is not recommended. The forecasts show why: the line is steeper than the one from part (a) and curves upward, so the gap between the two grows with the horizon. The report prints the ma1 estimate as -1.0000. Both models have d = 2, so the AICc values can be compared, and the constant model scores slightly better (AICc 196.79 in the report above, against 198.32 for the model in part (a)). I would still be cautious about using it, because the quadratic trend extrapolates an accelerating rise in passengers, which the book warns against, and a slightly lower AICc is not enough to justify that.

Exercise 9.8

For the United States GDP series (from global_economy):

us_gdp <- global_economy |> filter(Code == "USA")
us_gdp |> autoplot(GDP)

(a) if necessary, find a suitable Box-Cox transformation for the data;

lambda_us <- us_gdp |> features(GDP, guerrero) |> pull(lambda_guerrero)
lambda_us
## [1] 0.2819443
us_gdp |> autoplot(box_cox(GDP, lambda_us))

US GDP curves upward at an increasing rate, so a transformation is worthwhile. Guerrero’s method gives lambda of about 0.28. After the Box-Cox transformation the series is close to a straight upward line, with a small dip around 2009.

(b) fit a suitable ARIMA model to the transformed data using ARIMA();

fit_us <- us_gdp |> model(auto = ARIMA(box_cox(GDP, lambda_us)))
report(fit_us)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda_us) 
## 
## 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

ARIMA() selected an ARIMA(1,1,0) with drift on the transformed data. It takes one difference, adds one autoregressive term (ar1 = 0.4586), and includes a constant, which gives the forecasts an upward drift on the transformed scale.

(c) try some other plausible models by experimenting with the orders chosen;

us_gdp |>
  gg_tsdisplay(difference(box_cox(GDP, lambda_us)), plot_type = "partial")
## 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()`).

fit_us_all <- us_gdp |>
  model(auto     = ARIMA(box_cox(GDP, lambda_us)),
        search   = ARIMA(box_cox(GDP, lambda_us), stepwise = FALSE, approx = FALSE),
        arima110 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(1, 1, 0)),
        arima011 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(0, 1, 1)),
        arima210 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(2, 1, 0)))
fit_us_all |> pivot_longer(!Country, names_to = "Model name", values_to = "Orders")
## # A mable: 5 x 3
## # Key:     Country, Model name [5]
##   Country       `Model name`                  Orders
##   <fct>         <chr>                        <model>
## 1 United States auto         <ARIMA(1,1,0) w/ drift>
## 2 United States search       <ARIMA(1,1,0) w/ drift>
## 3 United States arima110     <ARIMA(1,1,0) w/ drift>
## 4 United States arima011     <ARIMA(0,1,1) w/ drift>
## 5 United States arima210     <ARIMA(2,1,0) w/ drift>
glance(fit_us_all) |> arrange(AICc) |> select(.model, AIC, AICc, BIC) |> knitr::kable(digits = 2)
.model AIC AICc BIC
auto 656.65 657.10 662.78
search 656.65 657.10 662.78
arima110 656.65 657.10 662.78
arima011 658.73 659.18 664.86
arima210 658.64 659.41 666.82

The ACF and the PACF of the differenced series each have a significant bar at lag 1 and nothing outside the dashed lines after that, so an AR(1) and an MA(1) are both reasonable candidates, and the automatic choice is an AR(1). The plots can only identify a pure AR or a pure MA, not a mixed model, so beyond that I tried neighbouring orders, as the book’s stepwise search does by adding or removing one AR or MA term. I compared that with a harder search (stepwise = FALSE, approx = FALSE), which also returned ARIMA(1,1,0) with drift, and with an ARIMA(0,1,1) and an ARIMA(2,1,0), both with drift. The AICc values in the table are 657.10 for the (1,1,0), 659.18 for the (0,1,1) and 659.41 for the (2,1,0), so the (1,1,0) has the lowest AICc of the candidates I tried (tied with the automatic searches, which chose the same model), and extra terms do not improve it. All models have d = 1, so the AICc comparison is valid.

(d) choose what you think is the best model and check the residual diagnostics;

fit_us_all |> select(auto) |> gg_tsresiduals()

augment(fit_us_all) |>
  filter(.model == "auto") |>
  features(.innov, ljung_box, lag = 10, dof = 1)
## # A tibble: 1 × 4
##   Country       .model lb_stat lb_pvalue
##   <fct>         <chr>    <dbl>     <dbl>
## 1 United States auto      3.81     0.923

I chose the ARIMA(1,1,0) with drift, which has the lowest AICc among the candidates. Its residuals look like white noise: all ACF bars are inside the dashed lines and the Ljung-Box p-value is 0.923. I used 10 lags and set dof = 1 to account for the one estimated AR coefficient, which leaves 9 degrees of freedom for the test. The residual plot has one large negative value near 2009, and the histogram shows it as a lone outlier on the left; the data alone do not tell me why. Apart from that, the residuals are centred near zero and fairly evenly spread.

(e) produce forecasts of your fitted model. Do the forecasts look reasonable?

fit_us_all |>
  forecast(h = 10) |>
  filter(.model == "auto") |>
  autoplot(us_gdp)

The forecasts continue the upward curve of GDP for the next ten years, with prediction intervals that widen with the horizon but stay fairly narrow around the central line. They look reasonable to me: the central line follows the steady growth of the last few decades. The intervals depend on the model and on the historical pattern continuing, and they leave out some sources of uncertainty, such as the estimated parameters and the choice of model, so they are probably narrower than the true uncertainty, which the book notes is a general weakness of ARIMA intervals.

(f) compare the results with what you would obtain using ETS() (with no transformation).

fit_us_ets <- us_gdp |> model(ets = ETS(GDP))
report(fit_us_ets)
## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9990876 
##     beta  = 0.5011949 
## 
##   Initial states:
##          l[0]        b[0]
##  448093333334 64917355687
## 
##   sigma^2:  7e-04
## 
##      AIC     AICc      BIC 
## 3190.787 3191.941 3201.089
fit_us_ets |> forecast(h = 10) |> autoplot(us_gdp)

# AICc cannot compare ETS with ARIMA, so I use a training and test set as in Section 5.8:
# train through 2007, forecast the 10 years 2008-2017, and compare against what happened.
# lambda is re-estimated on the training years only, so no later data leaks into the model
train_us <- us_gdp |> filter(Year <= 2007)
lambda_train <- train_us |> features(GDP, guerrero) |> pull(lambda_guerrero)
lambda_train
## [1] 0.1897405
fit_us_test <- train_us |>
  model(arima = ARIMA(box_cox(GDP, lambda_train)),
        ets = ETS(GDP))
fit_us_test
## # A mable: 1 x 3
## # Key:     Country [1]
##   Country                         arima          ets
##   <fct>                         <model>      <model>
## 1 United States <ARIMA(1,1,0) w/ drift> <ETS(M,A,N)>
fc_us_test <- fit_us_test |> forecast(h = 10)

# training accuracy and test accuracy side by side, as in the book's ARIMA vs ETS example
bind_rows(accuracy(fit_us_test), accuracy(fc_us_test, us_gdp)) |>
  select(.model, .type, RMSE, MAE, MAPE)
## # A tibble: 4 × 5
##   .model .type       RMSE     MAE  MAPE
##   <chr>  <chr>      <dbl>   <dbl> <dbl>
## 1 arima  Training 7.92e10 5.61e10  1.47
## 2 ets    Training 9.64e10 7.11e10  1.97
## 3 arima  Test     2.75e12 2.54e12 14.8 
## 4 ets    Test     1.72e12 1.66e12  9.93
fc_us_test |> autoplot(us_gdp |> filter(Year >= 1990), level = NULL)

On the full series the two models give similar forecasts in direction, but the ETS intervals are much wider and widen faster. For accuracy, I trained both on 1960 to 2007 and forecast 2008 to 2017. The table shows training and test accuracy. On the training years the ARIMA has the lower RMSE (7.92e10 against 9.64e10 for ETS), but on the test set the ETS model is more accurate on all three measures. This is one ten-year test period, and it includes the 2009 drop in the series, so it shows how the two methods did over that particular decade and does not settle which is better in general. It also answers a different question from comparing the two models by their AICc, which is not valid across model classes.