# Loading required packages
library(fpp3)DATA 624 Homework 6
9.1
Figure 9.32 shows the ACF for 36 random numbers, 360 random numbers and 1,000 random numbers.
The figure in the book is a picture, so I simulate the same three cases with random normal numbers.
set.seed(624)
wn36 <- tsibble(idx = 1:36, y = rnorm(36), index = idx)
wn360 <- tsibble(idx = 1:360, y = rnorm(360), index = idx)
wn1000 <- tsibble(idx = 1:1000, y = rnorm(1000), index = idx)
wn36 |> ACF(y) |> autoplot() + labs(title = "36 random numbers")
wn360 |> ACF(y) |> autoplot() + labs(title = "360 random numbers")
wn1000 |> ACF(y) |> autoplot() + labs(title = "1000 random numbers")
a. Explain the differences among these figures. Do they all indicate that the data are white noise?
Yes, all three look like white noise. Almost every spike stays inside the blue bounds, and the few that cross (2 out of 30 lags for the 1000 numbers, ~ 7%) are close to the 5% we expect by chance. The difference is the scale: the largest spike is 0.288 with 36 numbers, 0.090 with 360 and 0.085 with 1000, and the bounds get tighter as the sample grows.
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 bounds depend only on the sample size. They are about 1.96 divided by the square root of the number of observations, which gives 0.327 for 36 numbers, 0.103 for 360 and 0.062 for 1000. For white noise the sample autocorrelation has a spread of about one over the square root of the sample size, and the 95% band is 1.96 of those. More data means a less noisy estimate, so the band shrinks. The autocorrelations also differ because each figure is a different random sample, and a sample estimate is never exactly 0. Its typical size gets smaller as the sample grows.
9.2
A classic example of a nonstationary series is the daily closing IBM stock price series (data set gafa_stock). Use R to plot the daily closing prices for IBM and compute the ACF and PACF. Explain how each plot shows that the series is nonstationary and should be differenced.
The gafa_stock data has no IBM series, only AAPL, AMZN, FB and GOOG, so I use Google instead. Stock prices skip weekends and holidays, so the dates are not regular. I index by trading day so the ACF can be computed.
goog <- gafa_stock |>
filter(Symbol == "GOOG") |>
mutate(trading_day = row_number()) |>
update_tsibble(index = trading_day, regular = TRUE)
goog |> gg_tsdisplay(Close, plot_type = "partial")
The time plot of Google’s daily closing price rises over time with no fixed mean, and it also jumps up around trading day 400, so the level depends on when you look. The ACF starts near 1 and decays very slowly, still about 0.95 at lag 30, far outside the blue bounds at every lag shown. The PACF has one spike of about 1 at lag 1 and almost nothing after. A slowly decaying ACF with a single PACF spike near 1 is the signature of a unit root, so the series is not stationary and should be differenced.
goog |> gg_tsdisplay(difference(Close), 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()`).
After one difference the series moves around a constant mean and almost all ACF and PACF spikes sit inside the bounds, so one difference is enough.
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.
tr <- global_economy |> filter(Country == "Turkey")
tr |> autoplot(GDP)
GDP grows slowly at first and then climbs faster, so the spread grows with the level and a transformation should help.
lambda_tr <- tr |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda_tr[1] 0.1572187tr <- tr |> mutate(gdp_bc = box_cox(GDP, lambda_tr))
d_tr <- tr |> features(gdp_bc, unitroot_ndiffs) |> pull(ndiffs)
d_tr[1] 1The Guerrero lambda is 0.157 and the unit root test asks for 1 difference(s) of the transformed series.
tr |> gg_tsdisplay(difference(gdp_bc), 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()`).
b. Accommodation takings in the state of Tasmania from aus_accommodation.
tas <- aus_accommodation |> filter(State == "Tasmania")
tas |> autoplot(Takings)
Takings rise over time and follow a clear yearly pattern, so we need a transformation and probably a seasonal difference.
lambda_tas <- tas |>
features(Takings, features = guerrero) |>
pull(lambda_guerrero)
tas <- tas |> mutate(tak_bc = box_cox(Takings, lambda_tas))
lambda_tas[1] 0.001819643D_tas <- tas |> features(tak_bc, unitroot_nsdiffs) |> pull(nsdiffs)
D_tas[1] 1The seasonal difference comes first, because it can also remove the need for a regular difference. Then I test the seasonally differenced series.
tas <- tas |> mutate(tak_sd = difference(tak_bc, 4))
d_tas <- tas |> features(tak_sd, unitroot_ndiffs) |> pull(ndiffs)
d_tas[1] 0tas |> gg_tsdisplay(difference(tak_sd), plot_type = "partial")Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_point()`).
The Guerrero lambda is 0.002, the tests ask for 1 seasonal difference(s) at lag 4 and then 0 regular difference(s). The plot shows the series after both differences.
c. Monthly sales from souvenirs.
souvenirs |> autoplot(Sales)
The December peaks get bigger as sales grow, so the seasonal size depends on the level and a transformation is needed.
lambda_sv <- souvenirs |>
features(Sales, features = guerrero) |>
pull(lambda_guerrero)
sv <- souvenirs |> mutate(sales_bc = box_cox(Sales, lambda_sv))
lambda_sv[1] 0.002118221D_sv <- sv |> features(sales_bc, unitroot_nsdiffs) |> pull(nsdiffs)
D_sv[1] 1sv <- sv |> mutate(sales_sd = difference(sales_bc, 12))
d_sv <- sv |> features(sales_sd, unitroot_ndiffs) |> pull(ndiffs)
d_sv[1] 0sv |> gg_tsdisplay(difference(sales_sd), plot_type = "partial")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()`).
The Guerrero lambda is 0.002, so the transformation is close to a log. The tests ask for 1 seasonal difference(s) at lag 12 and then 0 regular difference(s), and the plot shows the series after both.
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(12345678)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
myseries |> autoplot(Turnover)
The seasonal swings grow with the level, so I start with a Box Cox transformation.
lambda_rt <- myseries |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
myseries <- myseries |> mutate(turn_bc = box_cox(Turnover, lambda_rt))
lambda_rt[1] 0.08303631D_rt <- myseries |> features(turn_bc, unitroot_nsdiffs) |> pull(nsdiffs)
D_rt[1] 1myseries <- myseries |> mutate(turn_sd = difference(turn_bc, 12))
d_rt <- myseries |> features(turn_sd, unitroot_ndiffs) |> pull(ndiffs)
d_rt[1] 0myseries |> gg_tsdisplay(difference(turn_sd), plot_type = "partial")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()`).
The Guerrero lambda is 0.083. The tests ask for 1 seasonal difference(s) at lag 12 and then 0 regular difference(s). The plot shows the transformed series after both differences.
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 phi1 = 0.6 and sigma squared = 1. The process starts with y1 = 0.
set.seed(123)
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 phi1?
sim |> autoplot(y)
To see the effect of phi1 I reuse the same errors e and only change phi1.
ar1 <- function(phi, e) {
y <- numeric(length(e))
for (i in 2:length(e)) {
y[i] <- phi * y[i - 1] + e[i]
}
y
}
ar_all <- tibble(
idx = rep(1:100, 3),
phi = rep(c("phi = 0.2", "phi = 0.6", "phi = 0.9"), each = 100),
y = c(ar1(0.2, e), ar1(0.6, e), ar1(0.9, e))
) |>
as_tsibble(index = idx, key = phi)
ar_all |> autoplot(y) + facet_grid(vars(phi))
Each value keeps a share phi1 of the previous value, so the memory fades quickly when phi1 is small and slowly when it is large. With 0.2 the series jumps around zero like white noise. With 0.6 it is smoother and stays on one side of zero a little longer. With 0.9 it makes long slow swings and looks closer to a random walk, which is what phi1 = 1 would give.
c. Write your own code to generate data from an MA(1) model with theta1 = 0.6 and sigma squared = 1.
ma1 <- function(theta, e) {
y <- numeric(length(e))
for (i in 2:length(e)) {
y[i] <- e[i] + theta * e[i - 1]
}
y
}
ma_sim <- tsibble(idx = 1:100, y = ma1(0.6, e), index = idx)d. Produce a time plot for the series. How does the plot change as you change theta1?
ma_sim |> autoplot(y)
ma_all <- tibble(
idx = rep(1:100, 3),
theta = rep(c("theta = 0.2", "theta = 0.6", "theta = 0.9"), each = 100),
y = c(ma1(0.2, e), ma1(0.6, e), ma1(0.9, e))
) |>
as_tsibble(index = idx, key = theta)
ma_all |> autoplot(y) + facet_grid(vars(theta))
An MA(1) only remembers one step back, so the plots look much more alike than the AR(1) ones. A larger theta1 smooths the series a bit more, but the effect stops at lag 1.
e. Generate data from an ARMA(1,1) model with phi1 = 0.6, theta1 = 0.6 and sigma squared = 1.
arma <- numeric(100)
for (i in 2:100) {
arma[i] <- 0.6 * arma[i - 1] + 0.6 * e[i - 1] + e[i]
}
arma11 <- tsibble(idx = 1:100, y = arma, index = idx)
arma11 |> autoplot(y)
f. Generate data from an AR(2) model with phi1 = -0.8, phi2 = 0.3 and sigma squared = 1. (Note that these parameters will give a nonstationary series.)
ar2 <- numeric(100)
for (i in 3:100) {
ar2[i] <- -0.8 * ar2[i - 1] + 0.3 * ar2[i - 2] + e[i]
}
ar2_sim <- tsibble(idx = 1:100, y = ar2, index = idx)
ar2_sim |> autoplot(y)
g. Graph the latter two series and compare them.
both <- tibble(
idx = rep(1:100, 2),
model = rep(c("ARMA(1,1)", "AR(2)"), each = 100),
y = c(arma, ar2)
) |>
as_tsibble(index = idx, key = model)
both |>
autoplot(y) +
facet_grid(vars(model), scales = "free_y") +
guides(colour = "none")
The ARMA(1,1) is stationary, and it stays between about negative 3.7 and 4.3 around zero. The AR(2) is not stationary because its parameters break the stationarity conditions. The values flip sign at almost every step and keep growing, reaching about negative 2177 and 2018 after 100 steps, which is why the AR(2) plot needs its own scale.
9.7
Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970 to 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.
aus_airpassengers |> autoplot(Passengers)
fit_a <- aus_airpassengers |>
model(auto = ARIMA(Passengers))
report(fit_a)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.65The automatic search selected an ARIMA(0,2,1). The series has a strong trend and no seasonality, so the double difference makes sense.
fit_a |> gg_tsresiduals()
lb <- augment(fit_a) |> features(.innov, ljung_box, lag = 10, dof = 1)
lb# A tibble: 1 × 3
.model lb_stat lb_pvalue
<chr> <dbl> <dbl>
1 auto 6.70 0.669The residuals are centered on 0 with no pattern, and the ACF spikes stay inside the bounds. The Ljung Box p value is 0.669, which is above 0.05, so the residuals look like white noise.
fc_a <- fit_a |> forecast(h = 10)
fc_a |> autoplot(aus_airpassengers)
b. Write the model in terms of the backshift operator.
The ARIMA(0,2,1) has no constant, so with the backshift operator B it is:
\[(1 - B)^2 y_t = (1 + \theta_1 B)\varepsilon_t\]
tidy(fit_a)# A tibble: 1 × 6
.model term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 auto ma1 -0.896 0.0594 -15.1 3.24e-19The estimate of theta1 is the ma1 value above, which is -0.896.
c. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
fit_c <- aus_airpassengers |>
model(drift010 = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)))
fc_c <- fit_c |> forecast(h = 10)
bind_rows(fc_a, fc_c) |>
autoplot(aus_airpassengers, level = NULL)
The ARIMA(0,1,0) with drift is a random walk with drift, so its forecast is a straight line with the average historical change as the slope. The ARIMA(0,2,1) is also close to a line, but its slope follows the most recent data instead of the full history. The intervals of the (0,2,1) also widen faster, because a double difference adds more uncertainty at each step.
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_d <- aus_airpassengers |>
model(
drift212 = ARIMA(Passengers ~ 1 + pdq(2, 1, 2)),
noconst212 = ARIMA(Passengers ~ 0 + pdq(2, 1, 2))
)
fc_d <- fit_d |> forecast(h = 10)
bind_rows(fc_a, fc_c, fc_d |> filter(.model == "drift212")) |>
autoplot(aus_airpassengers, level = NULL)
With drift, the ARIMA(2,1,2) forecast also rises along a straight line, close to the ones from parts a and c.
fc_d |>
filter(.model == "noconst212") |>
autoplot(aus_airpassengers)Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-InfWarning: Removed 10 rows containing missing values or values outside the scale range
(`geom_line()`).
Without the constant, a model with one difference has no drift, so the forecasts level off and stay flat near the last observed value instead of continuing to climb. That is a poor fit for a series that keeps growing.
e. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
fit_e <- aus_airpassengers |>
model(const021 = ARIMA(Passengers ~ 1 + pdq(0, 2, 1)))
fc_e <- fit_e |> forecast(h = 10)
bind_rows(fc_a, fc_e) |>
autoplot(aus_airpassengers, level = NULL)
With two differences, adding a constant puts a quadratic trend in the forecasts, so they curve upward faster than the straight lines from the earlier models. R warns that this is generally discouraged, since a quadratic trend over a long horizon is rarely realistic.
9.8
For the United States GDP series (from global_economy):
us <- global_economy |> filter(Code == "USA")
us |> autoplot(GDP)
a. If necessary, find a suitable Box Cox transformation for the data.
lambda <- us |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda[1] 0.2819443GDP grows faster over time, so a transformation helps. The Guerrero lambda is 0.282, and I use it in all the models below.
us |> autoplot(box_cox(GDP, lambda))
b. Fit a suitable ARIMA model to the transformed data using ARIMA().
fit_auto <- us |>
model(auto = ARIMA(box_cox(GDP, lambda)))
report(fit_auto)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.78The automatic search picked the model shown above, and its AICc is 657.1.
c. Try some other plausible models by experimenting with the orders chosen.
fit <- us |>
model(
auto = ARIMA(box_cox(GDP, lambda)),
search = ARIMA(box_cox(GDP, lambda), stepwise = FALSE, approximation = FALSE),
arima110 = ARIMA(box_cox(GDP, lambda) ~ pdq(1, 1, 0)),
arima011 = ARIMA(box_cox(GDP, lambda) ~ pdq(0, 1, 1)),
arima111 = ARIMA(box_cox(GDP, lambda) ~ pdq(1, 1, 1)),
arima210 = ARIMA(box_cox(GDP, lambda) ~ pdq(2, 1, 0))
)
glance(fit) |> select(.model, AICc) |> arrange(AICc)# A tibble: 6 × 2
.model AICc
<chr> <dbl>
1 auto 657.
2 search 657.
3 arima110 657.
4 arima011 659.
5 arima111 659.
6 arima210 659.All models use the same transformed data, so their AICc values can be compared. The lower the AICc, the better the balance between fit and the number of parameters.
d. Choose what you think is the best model and check the residual diagnostics.
best <- "arima110"fit |> select(all_of(best)) |> gg_tsresiduals()
lb_us <- augment(fit) |>
filter(.model == best) |>
features(.innov, ljung_box, lag = 10, dof = 1)
lb_us# A tibble: 1 × 4
Country .model lb_stat lb_pvalue
<fct> <chr> <dbl> <dbl>
1 United States arima110 3.81 0.923I chose the ARIMA(1,1,0) because it ties for the lowest AICc (657.10) with the auto and search models, which are the same model, and it is simpler than the ARIMA(0,1,1), (1,1,1) and (2,1,0) alternatives, whose AICc values are higher by 2.08, 2.31 and 2.31. The Ljung Box statistic is 3.81 and the p value is 0.923, far above 0.05, so there is no evidence of autocorrelation left in the residuals and they behave like white noise.
e. Produce forecasts of your fitted model. Do they look reasonable?
fc_best <- fit |>
forecast(h = 10) |>
filter(.model == best)fc_best |> autoplot(us)
The forecasts continue the upward path of the series and the interval widens as the horizon grows, which is what we expect. The interval is not symmetric because the forecasts are transformed back to the original scale. This looks reasonable for the next 10 years.
f. Compare the results with what you would obtain using ETS() (with no transformation).
fit_ets <- us |> model(ETS = ETS(GDP))
report(fit_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 fc_ets <- fit_ets |> forecast(h = 10)
bind_rows(fc_best, fc_ets) |>
autoplot(us, level = NULL)
acc <- bind_rows(
accuracy(fit) |> filter(.model == best),
accuracy(fit_ets)
) |>
select(.model, RMSE, MAE, MAPE)
acc# A tibble: 2 × 4
.model RMSE MAE MAPE
<chr> <dbl> <dbl> <dbl>
1 arima110 147232166284. 87539655369. 1.52
2 ETS 166906260796. 104653662551. 1.91AICc cannot be compared between the two, because the ARIMA is fitted to the transformed data and the ETS to the original GDP. The training RMSE, MAE and MAPE are on the original scale, so they can be compared. On the training data arima110 has the lower RMSE.