Rpubs URL: https://rpubs.com/umaisabdullah/DATA624_HW6

Objective

Complete and submit exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7 and 9.8 for Homework 6 from the Hyndman online Forecasting book (Chapter 9, ARIMA models). We have to submit both our Rpubs link as well as the .pdf file with our code.

Overview

Homework 5 used exponential smoothing, which describes a series through its level, trend and season. ARIMA models take a different route: they describe a series through its autocorrelations, the way each observation is related to the ones before it. Everything in this chapter hangs off three ideas.

The first is stationarity. ARIMA modelling assumes the series has a stable mean and variance and no predictable long-term pattern. Real series usually do not, so we transform them (Box-Cox to stabilize variance) and difference them (to remove trend and seasonality) until they do. The exercises in 9.2, 9.3 and 9.5 are all about deciding how much of that to do, using the ACF plot and the KPSS unit root test as evidence.

The second is the building blocks: an AR(p) term regresses the series on its own past values, an MA(q) term regresses it on past forecast errors, and the I(d) term is the number of differences taken. Exercise 9.6 simulates each of these so their behavior can be seen directly, and 9.7 and 9.8 fit them to real data.

The third is that ARIMA() automates the order selection, but the result is a starting point rather than an answer. Including a constant changes the long-run forecast shape dramatically (a constant with d = 1 gives a drift, with d = 2 gives a quadratic trend), and models with different amounts of differencing cannot be compared by AIC. Exercises 9.7 and 9.8 spend most of their effort on those two points.

Pre-Requisite

All of the data comes from the fpp3 package, so nothing needs to be downloaded. The only package required is fpp3, which loads the forecasting functions and every dataset used below.

Exercise 9.1

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

  1. Explain the differences among these figures. Do they all indicate that the data are white noise?
  2. 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?

My thought process. The figure is just three ACF plots of independent random numbers at three sample sizes. Rather than describe a picture I cannot see here, I regenerate the same experiment with a fixed seed and compute the critical values directly, since the answer to both parts comes down to how sample size affects them.

set.seed(624)
wn <- list(
  n36   = rnorm(36),
  n360  = rnorm(360),
  n1000 = rnorm(1000)
)

acf_one <- function(x, n) {
  tsibble(t = seq_along(x), y = x, index = t) |>
    ACF(y, lag_max = 20) |>
    autoplot() +
    labs(title = paste0("ACF of ", n, " random numbers"))
}

acf_one(wn$n36, 36)

acf_one(wn$n360, 360)

acf_one(wn$n1000, 1000)

tibble(T = c(36, 360, 1000)) |>
  mutate(critical_value = 1.96 / sqrt(T))
## # A tibble: 3 × 2
##       T critical_value
##   <dbl>          <dbl>
## 1    36         0.327 
## 2   360         0.103 
## 3  1000         0.0620
# how many of the first 20 autocorrelations fall outside the bounds
count_out <- function(x) {
  r <- acf(x, lag.max = 20, plot = FALSE)$acf[-1]
  sum(abs(r) > 1.96 / sqrt(length(x)))
}
c(n36 = count_out(wn$n36), n360 = count_out(wn$n360), n1000 = count_out(wn$n1000))
##   n36  n360 n1000 
##     0     0     2

(a) Differences among the figures

All three are the same underlying process, independent draws from a normal distribution, so all three are white noise, and in all three almost every spike sits inside the dashed blue bounds. What changes with sample size is how tightly the spikes cluster around zero. With 36 observations the spikes are comparatively large, reaching a few tenths, with 360 they are visibly smaller, and with 1,000 they are barely off the axis. The bounds themselves also pull in toward zero as the sample grows.

Yes, all three indicate white noise. A properly behaved white noise series is expected to have about 5 percent of its autocorrelations outside the 95 percent bounds purely by chance. In my simulation none of the first 20 autocorrelations crossed the bounds for the 36 and 360 point series, and two of the 20 did for the 1,000 point series, which is in line with the 1 in 20 that chance alone produces and is not evidence of structure. What would indicate a problem is many spikes outside the bounds, or a clear pattern in them.

(b) Why the critical values and the autocorrelations differ

The critical values are \(\pm 1.96/\sqrt{T}\), which is about 0.327 for T = 36, 0.103 for T = 360 and 0.062 for T = 1,000. The bounds are the 95 percent range of the sampling distribution of an autocorrelation when the series really is white noise, and that distribution has a standard deviation of roughly \(1/\sqrt{T}\). A larger sample estimates each autocorrelation more precisely, so the range of values that could arise by chance shrinks, and the bounds move closer to zero.

The autocorrelations differ between the figures for the same reason. The true autocorrelation of white noise is exactly zero at every lag, but what is plotted is an estimate from a finite sample, and an estimate carries sampling error. The smaller the sample, the larger that error, which is why the 36 point plot has the biggest spikes even though nothing is going on. The bounds and the spikes shrink together, so the conclusion drawn from each plot is the same.

Exercise 9.2

Question. A classic example of a non-stationary series is the daily closing IBM stock price series (data set gafa_stock). Use R to 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. Produce the same plots again after differencing the data.

My thought process. The prompt text mentions IBM and then asks for Amazon, and gafa_stock contains Google, Amazon, Facebook and Apple but no IBM, so I follow the explicit instruction and use Amazon. Stock prices are only recorded on trading days, so the raw dates have gaps every weekend. I re-index by trading day number so the time series functions see a regular series, the same approach the book takes.

amzn <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)

amzn |>
  autoplot(Close) +
  labs(title = "Amazon daily closing price", x = "Trading day", y = "Close (USD)")

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

Why this series is non-stationary. The time plot shows a strong upward trend with a large, sustained rise and no return to any fixed mean, so the mean is clearly not constant. The ACF starts near 1 and decays very slowly, staying significantly positive across the whole lag range, which is the signature of a trending series: each price is closely tied to the ones near it. The PACF has one enormous spike at lag 1, close to 1, and almost nothing after that, which is what a random walk looks like, where the current value is essentially the previous value plus noise. All three plots point the same way, and the standard remedy is to difference once.

amzn |>
  features(Close, list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 4
##   Symbol kpss_stat kpss_pvalue ndiffs
##   <chr>      <dbl>       <dbl>  <int>
## 1 AMZN        14.0        0.01      1

The KPSS test agrees. It tests the null hypothesis that the series is stationary, and the KPSS statistic of about 14 (p-value of 0.01 or smaller) rejects that very firmly, and unitroot_ndiffs recommends one difference.

amzn |>
  mutate(diff_close = difference(Close)) |>
  gg_tsdisplay(diff_close, plot_type = "partial")

amzn |>
  mutate(diff_close = difference(Close)) |>
  features(diff_close, unitroot_kpss)
## # A tibble: 1 × 3
##   Symbol kpss_stat kpss_pvalue
##   <chr>      <dbl>       <dbl>
## 1 AMZN       0.149         0.1

After differencing. The differenced series moves around zero with no trend, and the ACF and PACF collapse to almost nothing, with nearly all spikes inside the bounds. The KPSS statistic falls to about 0.15 with a p-value above 0.1, so we no longer reject stationarity. In practical terms, the day to day change in the Amazon price looks like noise, which is the classic random walk behavior: yesterday’s price is the best guide to today’s, and the change from one day to the next is not predictable from earlier changes.

Exercise 9.3

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

  1. Turkish GDP from global_economy.
  2. Accommodation takings in the state of Tasmania from aus_accommodation.
  3. Monthly sales from souvenirs.

My thought process. The same recipe applies to each series. First choose the Box-Cox lambda with the Guerrero method to stabilize the variance. Then, on the transformed series, check for seasonality first with unitroot_nsdiffs (the seasonal difference is taken before any ordinary difference, since a seasonal difference often removes the trend as well) and then for further differencing with unitroot_ndiffs, confirming the final result with the KPSS test.

(a) Turkish GDP

turkey <- global_economy |>
  filter(Country == "Turkey") |>
  select(Year, GDP)

turkey |>
  autoplot(GDP) +
  labs(title = "Turkish GDP", x = "Year", y = "GDP (current USD)")

lambda_turkey <- turkey |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)
lambda_turkey
## [1] 0.1572187
turkey |>
  features(box_cox(GDP, lambda_turkey), list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 3
##   kpss_stat kpss_pvalue ndiffs
##       <dbl>       <dbl>  <int>
## 1      1.50        0.01      1
turkey |>
  mutate(bc_diff = difference(box_cox(GDP, lambda_turkey))) |>
  features(bc_diff, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0889         0.1
turkey |>
  mutate(bc_diff = difference(box_cox(GDP, lambda_turkey))) |>
  autoplot(bc_diff) +
  labs(title = "Turkish GDP, Box-Cox transformed and differenced once",
       x = "Year", y = "Differenced Box-Cox GDP")

The Guerrero lambda is about 0.157, a value close to a log transformation, which suits a series whose growth is roughly proportional to its level. The transformed series is still strongly non-stationary (KPSS statistic about 1.5, p-value 0.01), and unitroot_ndiffs recommends one difference. After a single first difference the KPSS test no longer rejects, and the plot fluctuates around a stable mean. Since this is annual data there is no seasonal difference to consider. Transformation: Box-Cox with lambda of about 0.16. Differencing: d = 1.

(b) Tasmanian accommodation takings

tas <- aus_accommodation |>
  filter(State == "Tasmania") |>
  select(Date, Takings)

tas |>
  autoplot(Takings) +
  labs(title = "Tasmanian accommodation takings", x = "Quarter", y = "Takings ($ million)")

lambda_tas <- tas |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)
lambda_tas
## [1] 0.001819643
tas |>
  features(box_cox(Takings, lambda_tas),
           list(unitroot_kpss, unitroot_nsdiffs, unitroot_ndiffs))
## # A tibble: 1 × 4
##   kpss_stat kpss_pvalue nsdiffs ndiffs
##       <dbl>       <dbl>   <int>  <int>
## 1      1.84        0.01       1      1
tas |>
  mutate(bc_sdiff = difference(box_cox(Takings, lambda_tas), lag = 4)) |>
  features(bc_sdiff, list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 3
##   kpss_stat kpss_pvalue ndiffs
##       <dbl>       <dbl>  <int>
## 1     0.198         0.1      0
tas |>
  mutate(bc_sdiff = difference(box_cox(Takings, lambda_tas), lag = 4)) |>
  autoplot(bc_sdiff) +
  labs(title = "Tasmanian takings, Box-Cox transformed and seasonally differenced",
       x = "Quarter", y = "Seasonally differenced Box-Cox takings")

The seasonal swing grows with the level of the series, so a transformation is needed, and the Guerrero lambda comes out at essentially 0 (about 0.002), which is a log transformation. The transformed series has a strong quarterly seasonal pattern, and unitroot_nsdiffs asks for one seasonal difference (lag 4). Once that seasonal difference is taken, the KPSS statistic falls to about 0.20 with a p-value above 0.1 and unitroot_ndiffs asks for no further differencing, so the seasonal difference alone has been enough to remove the trend too. Transformation: log (lambda of about 0). Differencing: one seasonal difference at lag 4, and no ordinary difference, so D = 1 and d = 0.

(c) Monthly souvenir sales

souvenirs |>
  autoplot(Sales) +
  labs(title = "Monthly souvenir shop sales", x = "Month", y = "Sales ($)")

lambda_souv <- souvenirs |>
  features(Sales, features = guerrero) |>
  pull(lambda_guerrero)
lambda_souv
## [1] 0.002118221
souvenirs |>
  features(box_cox(Sales, lambda_souv), list(unitroot_kpss, unitroot_nsdiffs))
## # A tibble: 1 × 3
##   kpss_stat kpss_pvalue nsdiffs
##       <dbl>       <dbl>   <int>
## 1      1.79        0.01       1
souvenirs |>
  mutate(bc_sdiff = difference(box_cox(Sales, lambda_souv), lag = 12)) |>
  features(bc_sdiff, list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 3
##   kpss_stat kpss_pvalue ndiffs
##       <dbl>       <dbl>  <int>
## 1     0.240         0.1      0
souvenirs |>
  mutate(bc_sdiff = difference(box_cox(Sales, lambda_souv), lag = 12)) |>
  autoplot(bc_sdiff) +
  labs(title = "Souvenir sales, Box-Cox transformed and seasonally differenced",
       x = "Month", y = "Seasonally differenced Box-Cox sales")

The souvenir series has the clearest case for a transformation of the three. The big December spike and the general level both grow over time, so the seasonal swing is multiplicative, and the Guerrero lambda is again essentially 0 (about 0.002), a log transformation. There is one seasonal difference recommended at lag 12. After it, the KPSS p-value is above 0.1 and no further ordinary difference is needed. Transformation: log (lambda of about 0). Differencing: one seasonal difference at lag 12, D = 1 and d = 0.

Exercise 9.5

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

Find an appropriate order of differencing (after transformation if necessary) to obtain stationary data.

My thought process. This is the same NSW takeaway food series I have used in every homework so far, drawn with the same seed. It has the same shape as the Tasmanian and souvenir series above, strong trend plus seasonality that grows with the level, so I expect the same recipe to apply, but I run the tests rather than assume.

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

myseries |>
  as_tibble() |>
  distinct(State, Industry, `Series ID`)
## # A tibble: 1 × 3
##   State           Industry               `Series ID`
##   <chr>           <chr>                  <chr>      
## 1 New South Wales Takeaway food services A3349792X
myseries |>
  autoplot(Turnover) +
  labs(title = "NSW takeaway food services turnover",
       x = "Month", y = "Turnover ($ million)")

lambda_retail <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)
lambda_retail
## [1] 0.002144737
myseries |>
  features(box_cox(Turnover, lambda_retail),
           list(unitroot_kpss, unitroot_nsdiffs, unitroot_ndiffs))
## # A tibble: 1 × 6
##   State           Industry               kpss_stat kpss_pvalue nsdiffs ndiffs
##   <chr>           <chr>                      <dbl>       <dbl>   <int>  <int>
## 1 New South Wales Takeaway food services      7.07        0.01       1      1
myseries |>
  mutate(bc_sdiff = difference(box_cox(Turnover, lambda_retail), lag = 12)) |>
  features(bc_sdiff, list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 5
##   State           Industry               kpss_stat kpss_pvalue ndiffs
##   <chr>           <chr>                      <dbl>       <dbl>  <int>
## 1 New South Wales Takeaway food services    0.0657         0.1      0
myseries |>
  mutate(bc_sdiff = difference(box_cox(Turnover, lambda_retail), lag = 12)) |>
  gg_tsdisplay(bc_sdiff, plot_type = "partial", lag_max = 36)

Result. The Guerrero lambda is about 0.002, effectively a log transformation, matching the earlier finding that the seasonal amplitude grows with the level. The transformed series fails the KPSS test badly (statistic about 7.07), and unitroot_nsdiffs asks for one seasonal difference. After taking it (lag 12), the KPSS statistic falls to about 0.066 with a p-value above 0.1, and unitroot_ndiffs asks for no further differencing, so by the KPSS test one seasonal difference is enough. The differenced series has no remaining trend or obvious seasonal swing, but it is not a clean noise series either. It makes long excursions away from zero, and the ACF confirms this: the first several autocorrelations start near 0.9 and decay only slowly, and there is a run of negative autocorrelations around lags 24 to 36, with a PACF spike near lag 13. That is a lot of persistence left over.

So this is a borderline case. The formal test accepts stationarity at D = 1, d = 0, and I go with that, but the slowly decaying ACF means an extra first difference (D = 1, d = 1) would also be a defensible choice. Either way, the remaining short-run persistence and the seasonal structure around lag 12 and 24 are what the non-seasonal and seasonal AR and MA terms of an ARIMA model would then have to explain. Transformation: log (lambda of about 0). Differencing: one seasonal difference at lag 12, D = 1 and d = 0, with d = 1 as a reasonable alternative.

One more thing worth being careful about: unitroot_ndiffs run on the transformed series before any seasonal difference (the first table above) gives its own answer, and it does not know a seasonal difference is coming. The seasonal difference is taken first, which is the standard order, and once it has been taken unitroot_ndiffs asks for no further ordinary difference.

Exercise 9.6

Question. Simulate and plot some data from simple ARIMA models.

  1. 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\).
  2. Produce a time plot for the series. How does the plot change as you change \(\phi_1\)?
  3. Write your own code to generate data from an MA(1) model with \(\theta_1 = 0.6\) and \(\sigma^2 = 1\).
  4. Produce a time plot for the series. How does the plot change as you change \(\theta_1\)?
  5. Generate data from an ARMA(1,1) model with \(\phi_1 = 0.6\), \(\theta_1 = 0.6\) and \(\sigma^2 = 1\).
  6. 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.)
  7. Graph the latter two series and compare them.

My thought process. Simulation is the clearest way to see what each parameter does, because the truth is known. I write one small generator per model, use the same seed so different parameter values are compared on the same underlying shocks, and keep 100 observations as in the book.

(a) AR(1) with phi = 0.6

set.seed(624)
e <- rnorm(100)

sim_ar1 <- function(phi, e) {
  y <- numeric(length(e))
  for (i in 2:length(e)) y[i] <- phi * y[i - 1] + e[i]
  tsibble(idx = seq_along(y), y = y, index = idx)
}

ar_06 <- sim_ar1(0.6, e)

(b) Time plot, and varying phi

ar_06 |>
  autoplot(y) +
  labs(title = "AR(1), phi = 0.6", x = "t", y = "y")

bind_rows(
  sim_ar1(0.0,  e) |> as_tibble() |> mutate(phi = "phi = 0 (white noise)"),
  sim_ar1(0.6,  e) |> as_tibble() |> mutate(phi = "phi = 0.6"),
  sim_ar1(0.95, e) |> as_tibble() |> mutate(phi = "phi = 0.95"),
  sim_ar1(-0.8, e) |> as_tibble() |> mutate(phi = "phi = -0.8")
) |>
  ggplot(aes(idx, y)) +
  geom_line() +
  facet_wrap(~ phi, ncol = 1, scales = "free_y") +
  labs(title = "AR(1) simulations with the same shocks and different phi", x = "t", y = "y")

With phi of 0.6 the series is smooth-looking: a high value tends to be followed by another high value, so it makes slow excursions away from zero and drifts back. Changing phi controls how much memory the series has. With phi of 0 there is no memory and the series is plain white noise. As phi approaches 1 the excursions become longer and larger and the series starts to look like a random walk (phi of exactly 1 is one). With a negative phi, such as -0.8, the series flips sign nearly every step and looks jagged, because a positive value is followed by a negative one. All of these use the same underlying shocks, so the only thing changing between panels is how the shocks are propagated.

(c) MA(1) with theta = 0.6

sim_ma1 <- function(theta, e) {
  y <- numeric(length(e))
  for (i in 2:length(e)) y[i] <- e[i] + theta * e[i - 1]
  tsibble(idx = seq_along(y), y = y, index = idx)
}

ma_06 <- sim_ma1(0.6, e)

(d) Time plot, and varying theta

ma_06 |>
  autoplot(y) +
  labs(title = "MA(1), theta = 0.6", x = "t", y = "y")

bind_rows(
  sim_ma1(0.0,  e) |> as_tibble() |> mutate(theta = "theta = 0 (white noise)"),
  sim_ma1(0.6,  e) |> as_tibble() |> mutate(theta = "theta = 0.6"),
  sim_ma1(0.95, e) |> as_tibble() |> mutate(theta = "theta = 0.95"),
  sim_ma1(-0.8, e) |> as_tibble() |> mutate(theta = "theta = -0.8")
) |>
  ggplot(aes(idx, y)) +
  geom_line() +
  facet_wrap(~ theta, ncol = 1, scales = "free_y") +
  labs(title = "MA(1) simulations with the same shocks and different theta", x = "t", y = "y")

The MA(1) series looks much less persistent than the AR(1) one. Each value depends on only the current and previous shock, so the memory is exactly one step long and the series never wanders far from zero for long. A larger positive theta gives a slightly smoother look, since consecutive values share more of the same shock, and a negative theta gives a more jagged look, as consecutive shocks partly cancel. Compared with the AR(1) panels, even theta near 0.95 does not produce the long slow excursions that phi near 0.95 does, and that difference (a shock that fades after one step versus a shock that decays geometrically forever) is the essential difference between the two model types.

(e) ARMA(1,1) with phi = 0.6 and theta = 0.6

sim_arma11 <- function(phi, theta, e) {
  y <- numeric(length(e))
  for (i in 2:length(e)) y[i] <- phi * y[i - 1] + theta * e[i - 1] + e[i]
  tsibble(idx = seq_along(y), y = y, index = idx)
}

arma_11 <- sim_arma11(0.6, 0.6, e)

arma_11 |>
  autoplot(y) +
  labs(title = "ARMA(1,1), phi = 0.6, theta = 0.6", x = "t", y = "y")

(f) AR(2) with phi1 = -0.8 and phi2 = 0.3

sim_ar2 <- function(phi1, phi2, e) {
  y <- numeric(length(e))
  for (i in 3:length(e)) y[i] <- phi1 * y[i - 1] + phi2 * y[i - 2] + e[i]
  tsibble(idx = seq_along(y), y = y, index = idx)
}

ar_2 <- sim_ar2(-0.8, 0.3, e)

ar_2 |>
  autoplot(y) +
  labs(title = "AR(2), phi1 = -0.8, phi2 = 0.3", x = "t", y = "y")

The stationarity conditions for an AR(2) are \(\phi_2 - \phi_1 < 1\), \(\phi_2 + \phi_1 < 1\) and \(|\phi_2| < 1\). Here \(\phi_2 - \phi_1 = 0.3 + 0.8 = 1.1\), which violates the first condition, so this is a non-stationary process, as the exercise warns.

(g) Comparing the two series

bind_rows(
  arma_11 |> as_tibble() |> mutate(model = "ARMA(1,1): phi = 0.6, theta = 0.6"),
  ar_2 |> as_tibble() |> mutate(model = "AR(2): phi1 = -0.8, phi2 = 0.3")
) |>
  ggplot(aes(idx, y)) +
  geom_line() +
  facet_wrap(~ model, ncol = 1, scales = "free_y") +
  labs(title = "ARMA(1,1) (stationary) versus AR(2) (non-stationary)", x = "t", y = "y")

Note that the two panels have different y axes on purpose, otherwise the ARMA(1,1) would be a flat line. The ARMA(1,1) series is stationary: it bounces around zero with a roughly constant spread for the whole 100 observations, and its values stay within a few units of zero. The AR(2) series looks entirely different. It oscillates, flipping sign almost every step, and the size of the oscillations grows steadily over time, staying near zero for the first 40 or so observations and then reaching the thousands by the end, while the ARMA series never leaves the range of about -5 to 3. That is what non-stationarity from an explosive root looks like: the negative phi1 produces the sign flipping, and the parameter combination makes each swing bigger than the last instead of dying out. The ARMA model would be fine to forecast with a constant mean and variance. The AR(2) model would not, since it has no stable mean or variance to forecast.

Exercise 9.7

Question. 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.
  2. Write the model in terms of the backshift operator.
  3. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
  4. 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.
  5. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?

My thought process. This exercise is a tour of how the constant and the amount of differencing shape a forecast. I start with what ARIMA() selects, then vary the constant and the orders one at a time against that baseline so the effect of each change is isolated.

aus_airpassengers |>
  autoplot(Passengers) +
  labs(title = "Australian air passengers", x = "Year", y = "Passengers (millions)")

(a) Automatic selection, residuals and forecasts

fit_auto <- aus_airpassengers |>
  model(auto = ARIMA(Passengers))

report(fit_auto)
## 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_auto |>
  gg_tsresiduals()

augment(fit_auto) |>
  features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 auto      6.70     0.754
fit_auto |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers) +
  labs(title = "Australian air passengers, ARIMA(0,2,1) forecast",
       x = "Year", y = "Passengers (millions)")

ARIMA() selects an ARIMA(0,2,1) with no constant, with an MA coefficient of about -0.896. The residual time plot has no pattern, the ACF is inside the bounds and the Ljung-Box test at lag 10 has a p-value of about 0.75, so there is no evidence against white noise. The forecast continues the upward trend in a straight line, with prediction intervals that widen quickly with the horizon.

(b) Backshift notation

An ARIMA(0,2,1) with no constant is

\[(1 - B)^2 y_t = (1 + \theta_1 B)\varepsilon_t\]

and with the estimated coefficient, \(\theta_1 = -0.896\),

\[(1 - B)^2 y_t = (1 - 0.896B)\varepsilon_t.\]

The left side is the second difference, the change in the change, so this model says that the slope of the series follows a random walk with a little smoothing from the MA term. Since the MA coefficient is close to -1, the model is close to what would be obtained by differencing twice something that only needed to be differenced once, a sign of mild over-differencing that shows up again in part (e).

(c) ARIMA(0,1,0) with drift

fit_c <- aus_airpassengers |>
  model(drift = ARIMA(Passengers ~ pdq(0, 1, 0)))

report(fit_c)
## 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
bind_rows(
  fit_auto |> forecast(h = 10),
  fit_c    |> forecast(h = 10)
) |>
  autoplot(aus_airpassengers, level = 80) +
  labs(title = "ARIMA(0,2,1) versus ARIMA(0,1,0) with drift",
       x = "Year", y = "Passengers (millions)")

The ARIMA(0,1,0) with drift is a random walk with drift: each year the series moves by a constant of about 1.42 million passengers plus noise. Its point forecast is a straight line, like part (a), but the two differ in what they extrapolate. The drift model uses the average growth of the whole history, while the ARIMA(0,2,1) follows the recent slope, so the (0,2,1) line is steeper. The drift model’s prediction intervals are narrower at long horizons, since a model with one difference accumulates uncertainty more slowly than one with two.

The two models cannot be ranked by AIC. The (0,2,1) reports an AIC of about 198.0 and the drift model about 200.3, but the first is computed from the likelihood of the twice-differenced series and the second from the once-differenced series, so the numbers measure fit to two different data sets. The choice between them is really a choice about the trend: whether next decade’s growth should follow the long-run average (drift) or the recent slope (d = 2).

(d) ARIMA(2,1,2) with drift, and without the constant

fit_d <- aus_airpassengers |>
  model(
    with_drift = ARIMA(Passengers ~ 1 + pdq(2, 1, 2)),
    no_constant = ARIMA(Passengers ~ 0 + pdq(2, 1, 2), method = "ML")
  )

report(fit_d |> select(with_drift))
## 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
report(fit_d |> select(no_constant))
## Series: Passengers 
## Model: ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       1.8081  -0.8111  -1.9642  1.0000
## s.e.  0.1121   0.1126   0.1091  0.1102
## 
## sigma^2 estimated as 3.9:  log likelihood=-97.19
## AIC=204.38   AICc=205.88   BIC=213.52
# roots of the AR polynomial for the no-constant model;
# a modulus close to 1 means the model is close to needing another difference
ar_nc <- tidy(fit_d) |>
  filter(.model == "no_constant", grepl("^ar", term)) |>
  pull(estimate)
ar_root_modulus <- Mod(polyroot(c(1, -ar_nc)))
ar_root_modulus
## [1] 1.017281 1.211981
min_root <- min(ar_root_modulus)

bind_rows(
  fit_auto |> forecast(h = 10),
  fit_c    |> forecast(h = 10),
  fit_d    |> select(with_drift) |> forecast(h = 10)
) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(0,2,1), ARIMA(0,1,0) with drift and ARIMA(2,1,2) with drift",
       x = "Year", y = "Passengers (millions)")

fit_d |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(2,1,2), with and without the constant",
       x = "Year", y = "Passengers (millions)")

glance(fit_d) |> select(.model, AIC, AICc, BIC)
## # A tibble: 2 × 4
##   .model        AIC  AICc   BIC
##   <chr>       <dbl> <dbl> <dbl>
## 1 with_drift   204.  207.  215.
## 2 no_constant  204.  206.  214.

With drift. The ARIMA(2,1,2) with drift gives a forecast that is a straight upward line, and it sits almost on top of the ARIMA(0,1,0) with drift forecast from part (c), below the steeper ARIMA(0,2,1) line. The drift term of about 3.2 sets the slope, and the AR and MA terms add little to the point forecast. Its fair comparison is the (0,1,0) with drift, since both have d = 1, and it loses: its AICc is about 206.6 against about 200.6 for the random walk with drift. (It cannot be ranked against the (0,2,1) by AICc, for the reason given in part (c).) The extra AR and MA terms are not earning their place, and the fitted MA(2) coefficient is estimated at exactly 1.0, on the boundary of invertibility, a warning sign that the model is over-parameterized for 42 observations.

Without the constant. I had to use maximum likelihood (method = "ML") to fit this one, since the default estimation could not find a stationary AR part without the constant. I expected removing the constant to flatten the forecast, and it did not. The forecast still rises, and more steeply than the version with drift, reaching roughly 100 million passengers by the end of the horizon against about 87 million with drift. Without a constant the model has no built-in average growth rate, so a stationary ARMA part on the differenced series would eventually pull the forecast to a flat line. The reason it does not is in the AR coefficients reported above. The smallest modulus of the AR polynomial’s roots is about 1.017. That is very close to 1, the boundary of stationarity, which means the AR part is behaving almost like another unit root. In effect the model has found a way to difference the series a second time through its AR terms, so it carries the recent slope forward much like the ARIMA(0,2,1) does, rather than settling down to a flat forecast. The lesson is that removing the constant does not guarantee a flat forecast when d = 1. What the constant does is anchor the slope to the long-run average growth, and removing it hands the slope over to the AR and MA terms and the most recent data, which makes the forecast more sensitive to where the sample happens to end.

(e) ARIMA(0,2,1) with a constant

fit_e <- aus_airpassengers |>
  model(const = 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_e)
## 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
bind_rows(
  fit_auto |> forecast(h = 10),
  fit_e    |> forecast(h = 10)
) |>
  autoplot(aus_airpassengers, level = 80) +
  labs(title = "ARIMA(0,2,1) with and without a constant",
       x = "Year", y = "Passengers (millions)")

Including a constant in a model with d = 2 produces a quadratic trend, and fable warns about this (the warning is shown for this chunk only). The fitted constant is about 0.057 and the MA coefficient moves to -1.0000, the edge of the invertible region, both signs that the model is being pushed to an extreme. The forecast curves upward, above the no-constant version and with a wider band, so passenger numbers are projected to grow faster and faster each year. Interestingly the fit statistics look better (AIC of about 196.2 against 198.0, a fair comparison this time since both models have d = 2), but this is a case where I would not trust the criterion over judgment: an ever-accelerating forecast for passenger numbers is not believable, and the MA coefficient of exactly -1 says the second difference is over-differencing the data. The book’s general advice is to avoid a constant when d = 2 and I agree.

Exercise 9.8

Question. For the United States GDP series (from global_economy):

  1. if necessary, find a suitable Box-Cox transformation for the data;
  2. fit a suitable ARIMA model to the transformed data using ARIMA();
  3. try some other plausible models by experimenting with the orders chosen;
  4. choose what you think is the best model and check the residual diagnostics;
  5. produce forecasts of your fitted model. Do the forecasts look reasonable?
  6. compare the results with what you would obtain using ETS() (with no transformation).

My thought process. US GDP is a smooth, fast-growing series, so I expect a transformation to be needed and a single difference with drift to be enough. I will let ARIMA() make its choice and then test neighbors to see whether the automatic pick holds up. The important caution for part (c) is that AIC can only be compared across models with the same amount of differencing.

us_gdp <- global_economy |>
  filter(Code == "USA") |>
  select(Year, GDP)

us_gdp |>
  autoplot(GDP) +
  labs(title = "United States GDP", x = "Year", y = "GDP (current USD)")

(a) Box-Cox transformation

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

The series curves upward, so the variance grows with the level and a transformation helps. The Guerrero lambda is about 0.282, which is between a log and a square root, so I use that value directly.

(b) ARIMA on the transformed data

fit_us_auto <- us_gdp |>
  model(auto = ARIMA(box_cox(GDP, lambda_us)))

report(fit_us_auto)
## 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
us_gdp |>
  features(box_cox(GDP, lambda_us), list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 3
##   kpss_stat kpss_pvalue ndiffs
##       <dbl>       <dbl>  <int>
## 1      1.55        0.01      1
us_gdp |>
  mutate(bc_diff = difference(box_cox(GDP, lambda_us))) |>
  features(bc_diff, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.208         0.1

ARIMA() selects an ARIMA(1,1,0) with drift: one difference to remove the trend, an AR(1) coefficient of about 0.46 and a drift constant of about 118 on the transformed scale. The unit root tests in the second and third tables agree that one difference is what the series needs: the KPSS test rejects stationarity for the transformed series, unitroot_ndiffs recommends d = 1, and after one difference the KPSS test no longer rejects.

(c) Other plausible models

fit_us <- us_gdp |>
  model(
    auto  = ARIMA(box_cox(GDP, lambda_us)),
    ma1  = ARIMA(box_cox(GDP, lambda_us) ~ pdq(0, 1, 1)),
    arma  = ARIMA(box_cox(GDP, lambda_us) ~ pdq(1, 1, 1)),
    ar2   = ARIMA(box_cox(GDP, lambda_us) ~ pdq(2, 1, 0)),
    d2    = ARIMA(box_cox(GDP, lambda_us) ~ pdq(0, 2, 2))
  )

glance(fit_us) |>
  select(.model, sigma2, log_lik, AIC, AICc, BIC) |>
  arrange(AICc)
## # A tibble: 5 × 6
##   .model sigma2 log_lik   AIC  AICc   BIC
##   <chr>   <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 d2      5743.   -321.  648.  649.  654.
## 2 auto    5479.   -325.  657.  657.  663.
## 3 ma1     5689.   -326.  659.  659.  665.
## 4 arma    5580.   -325.  659.  659.  667.
## 5 ar2     5580.   -325.  659.  659.  667.

The automatic model is already the ARIMA(1,1,0) with drift, so I do not refit that order separately; the alternatives are its neighbors.

A caution about this table. The (0,2,2) model has the lowest AICc, about 649 against 657 for the rest, but I do not take that at face value. It has d = 2 while the others have d = 1, and information criteria are computed from the likelihood of the differenced data, so they are not comparable between models with different amounts of differencing. Comparing them would mean comparing likelihoods of two different series. Among the models that can be compared (all with d = 1), the automatic choice ARIMA(1,1,0) with drift has the lowest AICc, about 657.1. The alternatives, MA(1), ARMA(1,1) and AR(2), all come in about 2 AICc points higher, and adding a parameter to the automatic model (the ARMA(1,1) and the AR(2)) does not improve the fit. The simplest model is also the best of the comparable ones.

(d) Best model and residual diagnostics

fit_us |>
  select(auto) |>
  gg_tsresiduals()

augment(fit_us |> select(auto)) |>
  features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 auto      3.81     0.955

I choose the ARIMA(1,1,0) with drift. The residuals show no pattern in the time plot, the ACF stays inside the bounds and the histogram is roughly bell shaped. The Ljung-Box test at lag 10 gives a p-value of about 0.96, so there is no evidence of leftover autocorrelation and the model has captured the structure in the data.

(e) Forecasts

fit_us |>
  select(auto) |>
  forecast(h = 10) |>
  autoplot(us_gdp) +
  labs(title = "US GDP, ARIMA(1,1,0) with drift and Box-Cox, 10 year forecast",
       x = "Year", y = "GDP (current USD)")

The forecasts look reasonable. They continue the historical upward growth in a smooth curve (a straight line on the transformed scale, which becomes a gently accelerating curve after back-transforming), and the prediction intervals widen with the horizon in the way that they should, reflecting increasing uncertainty rather than pretending the next decade is as knowable as next year. The drift is estimated from the whole 1960 to 2017 history, so the forecast assumes growth continues at the historical average rate, and that is the assumption to keep in mind. It would not anticipate a recession, which is a limit of all of these models and not a flaw in this one.

(f) Comparison with ETS

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
bind_rows(
  fit_us |> select(auto) |> forecast(h = 10),
  fit_us_ets |> forecast(h = 10)
) |>
  autoplot(us_gdp, level = 80) +
  labs(title = "US GDP, ARIMA (Box-Cox) versus ETS (no transformation)",
       x = "Year", y = "GDP (current USD)")

bind_rows(
  fit_us |> select(auto) |> accuracy(),
  fit_us_ets |> accuracy()
) |>
  select(.model, RMSE, MAE, MAPE, MASE, ACF1)
## # A tibble: 2 × 6
##   .model          RMSE           MAE  MAPE  MASE   ACF1
##   <chr>          <dbl>         <dbl> <dbl> <dbl>  <dbl>
## 1 auto   147232166279.  87539655357.  1.52 0.257 0.0232
## 2 ets    166906260796. 104653662551.  1.91 0.307 0.153
bind_rows(
  fit_us |> select(auto) |> forecast(h = 10) |> filter(Year == 2027),
  fit_us_ets |> forecast(h = 10) |> filter(Year == 2027)
) |>
  as_tibble() |>
  select(.model, Year, .mean)
## # A tibble: 2 × 3
##   .model  Year   .mean
##   <chr>  <dbl>   <dbl>
## 1 auto    2027 2.82e13
## 2 ets     2027 2.61e13

ETS() on the untransformed data selects ETS(M,A,N), multiplicative error with an additive trend and no season, with alpha near 1 (about 0.999) and beta of about 0.50. The alpha near 1 means that the level essentially tracks the latest observation, so the model behaves like a random walk with a trend that adapts.

On training accuracy the ARIMA model wins on every measure: its RMSE is about 1.47e11 against about 1.67e11 for the ETS, and its MAPE is 1.5 percent against 1.9 percent. The ARIMA residuals are also essentially uncorrelated (ACF1 of 0.02) while the ETS residuals still carry some (ACF1 of 0.15), meaning ETS leaves a little predictable structure behind. (AICs cannot be compared here, since the ARIMA is fitted to the Box-Cox scale and the ETS to the original scale.) The forecasts for 2027 are close, about 28.2 trillion dollars for the ARIMA model against about 26.1 trillion for the ETS. The ETS forecast is lower, and I read that as a difference in trend shape: ETS extends a straight line on the original scale, while the ARIMA’s straight line on the transformed scale curves upward when converted back. The ETS intervals are also markedly wider, which is consistent with a multiplicative error term letting the uncertainty grow with the level. Given the better training fit and cleaner residuals I would use the ARIMA(1,1,0) with drift, though the honest summary is that on this series both are reasonable and the choice matters less over a short horizon than it would over a long one.

Summary

Across the seven exercises, one theme returned repeatedly: the constant, the transformation and the amount of differencing each have a large effect on the forecast, and they should be chosen on evidence rather than by default.