library(tidyverse)
library(fpp3)
library(fable)Assignment 6: ARIMA
Do the exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8 in Hyndman. Link to Chapter 9 exercises for reference.
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?
Each of these are indicative of white noise data. For white noise, ACF values should be close to zero at all lags, staying within the blue significance bounds. None of the ACF lags visibly go beyond the threshold that represents the 95% range for in being in scope for being considered white noise. They are different because the number of observations differs; the first plot is for 36 numbers, the middle is for 360, the right is for 1000. With a smaller sample size (left for 36), random fluctuations cause the ACF values to vary more, while larger samples smooth out these fluctuations more.
Figure 9.32: Left: ACF for a white noise series of 36 numbers. Middle: ACF for a white noise series of 360 numbers. Right: ACF for a white noise series of 1,000 numbers.
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 critical values (blue dashed lines) are at different distances from the mean of zero because they depend on the number of observations. As the number of observations increases, the confidence bounds narrow. That’s why the leftmost plot with the fewest observations has a wider critical value range compared to the plot with the most observations which has the narrowest range between its critical values.
The autocorrelations are different in each figure when they each refer to whie noise because of sample variability. While all three represent white noise, the smaller a sample size the more susceptible it is to random variation, causing individual ACF values to deviate more from 0 in contrast to a larger sample size which shows AF values closer to zero, which is the expected behavior of white noise. Fluctuations stabilize as the sample size increases.
2. A classic example of a non-stationary series are 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 |>
filter(Symbol == "AMZN") |>
gg_tsdisplay(Close, plot_type='partial') +
labs(title = "AMZN (Amazon) Daily Closing Stock Price")
The plots indicate that the Amazon (AMZN) daily closing stock prices are non-stationary.
For the time-series plot, we can see an upward trend over time (suggests a lack of stationarity in the mean) and there are periods of more rapid growth and fluctuation, which indicates potential for changes in variation over time. Even without looking at the ACF and PACF plots, we can see that our data is non-stationary based on these observations.
For the ACF plot, we see strong (high indicators of) positive correlations that decay over time as the lag increases (looking at the plot from left to right we see a steady but slow decay). This is another indicator of non-stationarity.
For the PACF plot, we see a significant spike at lag 1, and some significance in further trnds like at lag 5, lag 19, and lag 25. If we were looking at something that exhibited stationarity we would see a cutoff after the first few lags rather than a recurring significance in lags beyond the blue dotted line boundary.
We can expect that the data would need to be differenced in order to become stationary.
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_gdp <- global_economy |> filter(Country == "Turkey")turkey_gdp |> autoplot(GDP) +
labs(title = "Turkish GDP (No Transformations or Differencing Applied)")
Let’s find out the optimal lambda for the Box-Cox transformation using the guerrero feature.
lambda_turkey <- turkey_gdp |> features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda_turkey[1] 0.1572187And now we apply the transformation to the data using the optimal lambda
turkey_gdp <- turkey_gdp |>
mutate(GDP_transformed = box_cox(GDP, lambda_turkey))Now we have our Box-Cox transformed data. Let’s check check for stationarity using a KPSS test and then, as needed, check the order of differencing needed using unitroot_ndiffs.
kpss_turkey <- unitroot_kpss(turkey_gdp$GDP_transformed)
kpss_turkey kpss_stat kpss_pvalue
1.498625 0.010000 Since a high p-value would suggest stationarity (no differencing needed) this indicates a need for differencing. Let’s find out what it should be.
ndiffs_turkey <- unitroot_ndiffs(turkey_gdp$GDP_transformed)
ndiffs_turkeyndiffs
1 Now we know we need to apply first differencing.
if (ndiffs_turkey >= 1) {
turkey_gdp <- turkey_gdp |>
mutate(GDP_diff = difference(GDP_transformed, differences = ndiffs_turkey))
}
turkey_gdp |>
gg_tsdisplay(GDP_diff, plot_type = 'partial') +
labs(title = paste("Box-Cox Transformed and 1st Differenced Turkish GDP"))
The ACF and PACF plots, along with the series showing the differenced GDP, indicate we now have stationary data.
b. Accommodation takings in the state of Tasmania from aus_accommodation.
tasmania_accomodation <- aus_accommodation |>
filter(State == "Tasmania")tasmania_accomodation |> autoplot(Takings) +
labs(title = "Accomodation Takings in Tasmania (No Transformations or Differencing Applied)")
Now we’ll find the optimal lambda for Box-Cox transformation and then apply it to transform the data.
lambda_tasmania <- tasmania_accomodation |>
features(Takings, features = guerrero) |> pull(lambda_guerrero)
lambda_tasmania[1] 0.001819643Now we apply the transformation to the data using the optimal lambda
tasmania_accomodation <- tasmania_accomodation |>
mutate(Takings_transformed = box_cox(Takings, lambda_tasmania))Now that have our Box-Cox transformed data, let’s check check for stationarity using a KPSS test and then, as needed, check the order of differencing needed using unitroot_ndiffs.
kpss_tasmania <- unitroot_kpss(tasmania_accomodation$Takings_transformed)
kpss_tasmania kpss_stat kpss_pvalue
1.837734 0.010000 Since a high p-value would suggest stationarity (no differencing needed) this indicates a need for differencing. Let’s find out what it should be.
ndiffs_tasmania <-
unitroot_ndiffs(tasmania_accomodation$Takings_transformed)
ndiffs_tasmaniandiffs
1 Again we see we need to apply first differencing.
if (ndiffs_tasmania >= 1) {
tasmania_accomodation <- tasmania_accomodation |>
mutate(Takings_diff = difference(Takings_transformed, differences = ndiffs_tasmania))
}
tasmania_accomodation |>
gg_tsdisplay(Takings_diff, plot_type = 'partial') +
labs(title = paste("Box-Cox Transformed and 1st Differenced Tasmanian GDP"))
That didn’t quite work - it still looks like we have non-stationary data. It looks like we need to do seasonal differencing based on the consistent lags shown in the ACF plot above.
tasmania_accomodation <- tasmania_accomodation |>
mutate(
Takings_diff_seasonal = difference(Takings_diff, lag = 4)
)
tasmania_accomodation |>
gg_tsdisplay(Takings_diff_seasonal, plot_type = 'partial') +
labs(title = "Box-Cox Transformed, Differenced, & Seasonally Differenced Tasmania Takings")
unitroot_kpss(tasmania_accomodation$Takings_diff_seasonal) kpss_stat kpss_pvalue
0.04736784 0.10000000 We can confirm that the data is now stationary, seeing as the kpss p-value is greater than 0.05, which indicates stationarity.
Going forward, we should keep in mind that it’s advised to do seasonal differencing first, not second.
c. Monthly sales from souvenirs.
souvenirs |> autoplot(Sales) +
labs(title = "Monthly Sales from Souvenirs (No Transformations or Differencing Applied)")
Now we’ll find the optimal lambda for Box-Cox transformation and then apply it to transform the data.
lambda_souvenirs <- souvenirs |>
features(Sales, features = guerrero) |> pull(lambda_guerrero)
lambda_souvenirs[1] 0.002118221Now we apply the transformation to the data using the optimal lambda
souvenirs_transformed <- souvenirs |>
mutate(Sales_transformed = box_cox(Sales, lambda_souvenirs))Now that have our Box-Cox transformed data, let’s check check for stationarity using a KPSS test and then, as needed, check the order of differencing needed using unitroot_ndiffs.
kpss_souvenirs <- unitroot_kpss(souvenirs_transformed$Sales_transformed)
kpss_souvenirs kpss_stat kpss_pvalue
1.790425 0.010000 Since a high p-value would suggest stationarity (no differencing needed) this p-value 0.01 indicates a need for differencing. Let’s find out what it should be.
ndiffs_souvenirs <-
unitroot_ndiffs(souvenirs_transformed$Sales_transformed)
ndiffs_souvenirsndiffs
1 Again we see we need to apply first differencing. We also know we have seasonality in our data (based on the prior visual inspection) so let’s apply lag=12 to this differencing. We have learned that seasonal differencing should be applied first, as it may resolve the need for further differencing.
if (ndiffs_souvenirs >= 1) {
souvenirs_transformed <- souvenirs_transformed |>
mutate(Sales_diff_seasonal = difference(Sales_transformed, lag = 12))
}
souvenirs_transformed |>
gg_tsdisplay(Sales_diff_seasonal, plot_type = 'partial') +
labs(title = paste("Box-Cox Transformed & Seasonally Differenced Souvenir Sales"))
unitroot_kpss(souvenirs_transformed$Sales_diff_seasonal) kpss_stat kpss_pvalue
0.2397881 0.1000000 This data now appears stationary: a p-value > 0.05 means we fail to reject the null hypothesis, suggesting that the series is now stationary - and we no longer see seasonality in the data or the lags.
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(5889)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))Autoplot of myseries for reference:
autoplot(myseries)Plot variable not specified, automatically selected `.vars = Turnover`
Now we’ll find the optimal lambda for Box-Cox transformation and then apply it to transform the data.
lambda_myseries <- myseries |>
features(Turnover, features = guerrero) |> pull(lambda_guerrero)
lambda_myseries[1] 0.08542275Now we apply the transformation to the data using the optimal lambda
myseries <- myseries |>
mutate(Turnover_transformed = box_cox(Turnover, lambda_myseries))Now that have our Box-Cox transformed data, let’s check check for stationarity using a KPSS test and then, as needed, check the order of differencing needed using unitroot_ndiffs.
kpss_myseries <- unitroot_kpss(myseries$Turnover_transformed)
kpss_myseries kpss_stat kpss_pvalue
6.697358 0.010000 Since a high p-value would suggest stationarity (no differencing needed) this indicates a need for differencing. Let’s find out what it should be.
ndiffs_myseries <-
unitroot_ndiffs(myseries$Turnover_transformed)
ndiffs_myseriesndiffs
1 Once again we need to apply first differencing, and based on visual inspection of the data previously we saw seasonality, so let’s apply lag=12 knowing the data is monthly.
if (ndiffs_myseries >= 1) {
myseries <- myseries |>
mutate(Turnover_diff_seasonal = difference(Turnover_transformed, lag=12))
}
myseries |>
gg_tsdisplay(Turnover_diff_seasonal, plot_type = 'partial') +
labs(title = paste("Box-Cox Transformed and 1st Differenced MySeries"))
unitroot_kpss(myseries$Turnover_diff_seasonal) kpss_stat kpss_pvalue
0.54781259 0.03089807 This data does not yet appear stationary: a p-value > 0.05 means we fail to reject the null hypothesis, suggesting that the series is stationary. But ours is 0.03, so we do not fail to reject the null hypothesis. We need to apply the first order differencing in addition to the seasonal differencing.
myseries <- myseries |>
mutate(Turnover_diff = difference(Turnover_diff_seasonal, differences = 1))
myseries |>
gg_tsdisplay(Turnover_diff, plot_type = 'partial') +
labs(title = paste("Box-Cox Transformed and Differenced MySeries"))
unitroot_kpss(myseries$Turnover_diff) kpss_stat kpss_pvalue
0.01382537 0.10000000 Great! We now have stationarity.
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 ϕ1 = 0.6 and σ^2 = 1. The process starts with y1 = 0.
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 ϕ1?
ggplot(sim, aes(x = idx, y = y)) +
geom_line() +
labs(title = "Simulated AR(1) with ϕ1 = 0.6", y = "y_t")
# phi_1 = 1
for(i in 2:100)
y[i] <- 1*y[i-1] + e[i]
tsibble(idx = seq_len(100), y = y, index = idx) %>%
autoplot(y) +
labs(title = "AR(1) model with ϕ1 = 1 (Random Walk)", y = "y_t")
# phi_1 = 0
for(i in 2:100)
y[i] <- 0*y[i-1] + e[i]
tsibble(idx = seq_len(100), y = y, index = idx) %>%
autoplot(y) +
labs(title = "AR(1) model with ϕ1 = 0 (White Noise)", y = "y_t") 
# phi_1 = -0.6
for(i in 2:100)
y[i] <- -0.6*y[i-1] + e[i]
tsibble(idx = seq_len(100), y = y, index = idx) %>%
autoplot(y) +
labs(title = "AR(1) model with ϕ1 = -0.6 (Oscillates)",
y = "y_t") 
We can see that as ϕ1 changes, the appearance of the plot changes: the simulated AR(1) with ϕ1 = 0.6 shows stationarity; the simulated AR(1) model with ϕ1 = 1 reveals what looks like a Random Walk; the simulated AR(1) model with ϕ1 = 0 looks like White Noise; and the AR(1) model with ϕ1 = -0.6 oscillates.
c. Write your own code to generate data from an MA(1) model with θ1 = 0.6 and σ^2 = 1.
See below.
d. Produce a time plot for the series. How does the plot change as you change θ1?
ma1_white_noise <- arima.sim(n = 100, list(ma = 0))
ma1_positive <- arima.sim(n = 100, list(ma = 0.6, sd = 1))
ma1_negative <- arima.sim(n = 100, list(ma = -0.6))
ma1_noninvertible <- arima.sim(n = 100, list(ma = 1.2))
ma1_white_noise_tbl <- tibble(time = 1:100, value = as.numeric(ma1_white_noise)) |>
as_tsibble(index = time)
ma1_positive_tbl <- tibble(time = 1:100, value = as.numeric(ma1_positive)) |>
as_tsibble(index = time)
ma1_negative_tbl <- tibble(time = 1:100, value = as.numeric(ma1_negative)) |>
as_tsibble(index = time)
ma1_noninvertible_tbl <- tibble(time = 1:100, value = as.numeric(ma1_noninvertible)) |>
as_tsibble(index = time)
autoplot(ma1_white_noise_tbl, value) + ggtitle("MA(1) with θ1 = 0 (White Noise)")
autoplot(ma1_positive_tbl, value) + ggtitle("MA(1) with θ1 = 0.6 (Weighted influence of past errors)")
autoplot(ma1_negative_tbl, value) + ggtitle("MA(1) with θ1 = -0.6 (Oscillating)")
autoplot(ma1_noninvertible_tbl, value) + ggtitle("MA(1) with θ1 = 1.2 (Non-Invertible)")
We can see that as θ1 changes, the appearance of the plot changes: the simulated MA(1) with θ1 = 0.6 shows the weighted influence of past errors in effect; the simulated with θ1 = 0 appears like White Noise (no dependency on past errors); the MA(1) with θ1 = -0.6 appears like an oscillating pattern (negative autocorrelation); and the MA(1) with θ1 = 1.2 is non-invertible (causes unstable estimation).
e. Generate data from an ARMA(1,1) model with ϕ1 = 0.6, θ1 = 0.6 and σ^2 = 1.
y <- numeric(100)
for(i in 2:100)
y[i] <- e[i] + 0.6 * e[i-1] + 0.6 * y[i-1]
arma_11 <- tsibble(idx = seq_len(100), y = y, index = idx)f. Generate data from an AR(2) model with ϕ1 = −0.8, ϕ2 = 0.3 and σ2 = 1. (Note that these parameters will give a non-stationary series.)
y <- numeric(100)
for(i in 3:100)
y[i] <- -0.8 * y[i-1] + 0.3 * y[i-2] + e[i]
ar_2 <- tsibble(idx = seq_len(100), y = y, index = idx)g. Graph the latter two series and compare them.
arma_11 |> gg_tsdisplay(y,plot_type = "partial") +
labs(title = "Simulated ARMA(1,1) Series (ϕ1 = 0.6, θ1 = 0.6)", y = "Value")
unitroot_kpss(arma_11$y) kpss_stat kpss_pvalue
0.1108754 0.1000000 ar_2 |> gg_tsdisplay(y,plot_type = "partial") +
labs(title = "Simulated AR(2) Series (ϕ1 = −0.8, ϕ2 = 0.3 and σ2 = 1)")
We can see that the non-stationary AR(2) exhibits what we learned about when phi is negative; oscillatory behavior. Comparatively, the ARMA(1,1) series does appear stationary.
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.
aus_airpassengers_fit <- aus_airpassengers |>
filter(Year < 2012) |>
model(ARIMA(Passengers))
report(aus_airpassengers_fit)Series: Passengers
Model: ARIMA(0,2,1)
Coefficients:
ma1
-0.8756
s.e. 0.0722
sigma^2 estimated as 4.671: log likelihood=-87.8
AIC=179.61 AICc=179.93 BIC=182.99We can see that an ARIMA(0,2,1) model was selected.
aus_airpassengers_fit |>
gg_tsresiduals() +
labs(title = "Residuals for selected ARIMA(0,2,1) model")
We can confirm that yes, the residuals look like white noise.
aus_airpassengers_fit |>
forecast(h=10) |>
autoplot(aus_airpassengers) +
labs(title = "Passengers on Australian air carriers 1970-2011, ARIMA(0,2,1)",
y = "Passengers (millions)")
Above are the forecasts for the next 10 periods.
b. Write the model in terms of the backshift operator.
y_t = -0.8756_t-1 + E_t
(1-B)^2y_t = (1 - 0.8756B)E_t
c. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
aus_airpassengers_fit2 <- aus_airpassengers |>
filter(Year < 2012) |>
model(ARIMA(Passengers ~ pdq(0,1,0)))
aus_airpassengers_fit2 |>
gg_tsresiduals() +
labs(title = "Residuals for ARIMA(0,1,0)")
aus_airpassengers_fit2 |>
forecast(h=10) |>
autoplot(aus_airpassengers) +
labs(title = "Australian Aircraft Passengers with ARIMA(0,1,0)",
y = "Passengers (millions)")
We can see the results are quite different in terms of the forecast. The ARIMA(0,2,1) forecasts values that are higher than the actual values from 2012 onward, while ARIMA(0,1,0) with drift forecasts values that are lower than the actual values from 2012 onward. The ARIMA(0,1,0) which includes drift seems more reasonable compared to the ARIMA(0,2,1).
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.
aus_airpassengers_fit3 <- aus_airpassengers |>
filter(Year < 2012) |>
model(ARIMA(Passengers ~ pdq(2,1,2)))
aus_airpassengers_fit3 |>
gg_tsresiduals() +
labs(title = "Residuals for ARIMA(2,1,2)")
aus_airpassengers_fit3 |>
forecast(h=10) |>
autoplot(aus_airpassengers) +
labs(title = "Australian Aircraft Passengers with ARIMA(2,1,2)",
y = "Passengers (millions)")
The ARIMA(2,1,2) model appears to fall in the ‘middle’ of the prior two options. Similarly also we see that the residuals appear to be white noise. The slope is more intensified in this model compared to the previous one.
Now if we remove the constant…
aus_airpassengers_fit3 <- aus_airpassengers |>
filter(Year < 2012) |>
model(ARIMA(Passengers ~ 0 + pdq(2,1,2)))Warning: 1 error encountered for ARIMA(Passengers ~ 0 + pdq(2, 1, 2))
[1] non-stationary AR part from CSSreport(aus_airpassengers_fit3)Series: Passengers
Model: NULL model
NULL model… we get a NULL model.
e. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
aus_airpassengers_fit4 <- aus_airpassengers |>
filter(Year < 2012) |>
model(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.We get an interesting warning message here encouraging us to consider removing the constant or reducing the number of differences.
aus_airpassengers_fit4 |>
gg_tsresiduals() +
labs(title = "Residuals for ARIMA(0,2,1) with constant")
aus_airpassengers_fit4 |>
forecast(h=10) |>
autoplot(aus_airpassengers) +
labs(title = "Australian Aircraft Passengers with ARIMA(0,2,1) with constant",
y = "Passengers (millions)")
We see the result of the constant here. The slope is more intensified, which makes sense considering the warning we were given above with implementing this model. That said, the residuals are still appearing to be white noise.
8. For the United States GDP series (from global_economy):
us_gdp <- global_economy |> filter(Country == "United States")us_gdp |> autoplot(GDP) +
labs(title = "United States GDP (No Transformations or Differencing Applied)")
us_gdp |> gg_tsdisplay(GDP, plot_type = "partial")
a. if necessary, find a suitable Box-Cox transformation for the data;
Is it necessary? Well, we can tell by looking at the output above that the data is definitely not stationary, so let’s apply a box-cox transformation first. We’ll find the optimal lambda for Box-Cox transformation and then apply it to transform the data.
lambda_us_gdp <- us_gdp |>
features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda_us_gdp[1] 0.2819443Now we apply the transformation to the data using the optimal lambda
us_gdp <- us_gdp |>
mutate(GDP_Transformed = box_cox(GDP, lambda_us_gdp))Now that have our Box-Cox transformed data, let’s check check for stationarity using a KPSS test and then, as needed, check the order of differencing needed using unitroot_ndiffs.
kpss_us_gdp <- unitroot_kpss(us_gdp$GDP_Transformed)
kpss_us_gdp kpss_stat kpss_pvalue
1.54817 0.01000 Since a high p-value would suggest stationarity (no differencing needed) this p-value 0.01 indicates a need for differencing. Let’s find out what it should be.
ndiffs_us_gdp <-
unitroot_ndiffs(us_gdp$GDP_Transformed)
ndiffs_us_gdpndiffs
1 We won’t yet apply differencing – we’ll leave that to the ARIMA model selection to determine the proper d, but it’s helpful to know this information regardless.
b. fit a suitable ARIMA model to the transformed data using ARIMA();
us_gdp_fit <- us_gdp |>
model(ARIMA(GDP_Transformed))
report(us_gdp_fit)Series: GDP_Transformed
Model: ARIMA(1,1,0) w/ drift
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.78As expected, a d = 1 is what the ARIMA model suggests, along with a p = 1 with drift. The model selected is ARIMA(1,1,0) with drift.
c. try some other plausible models by experimenting with the orders chosen;
We know based on the prior steps taken that differencing of 1 is suggested, so let’s look at some variation on other options.
us_gdp_fit_all <- us_gdp |>
model(
autoarima110 = ARIMA(GDP_Transformed),
arima011 = ARIMA(GDP_Transformed ~ pdq(0,1,1)),
arima111 = ARIMA(GDP_Transformed ~ pdq(1,1,1)),
arima210 = ARIMA(GDP_Transformed ~ pdq(2,1,0)),
arima012 = ARIMA(GDP_Transformed ~ pdq(0,1,2))
)
model_comparison <- glance(us_gdp_fit_all) |>
arrange(AICc) |>
select(.model, sigma2, log_lik, AIC, AICc, BIC)
model_comparison# A tibble: 5 × 6
.model sigma2 log_lik AIC AICc BIC
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 autoarima110 5479. -325. 657. 657. 663.
2 arima011 5689. -326. 659. 659. 665.
3 arima111 5580. -325. 659. 659. 667.
4 arima210 5580. -325. 659. 659. 667.
5 arima012 5634. -326. 659. 660. 667.d. choose what you think is the best model and check the residual diagnostics;
The model selected by the ARIMA() process, called autoarima110 above, is the best model. It has the lowest AICc. Let’s look at its residual diagnostics:
us_gdp_fit_all |>
select(autoarima110) |>
gg_tsresiduals()
The residual time series plot looks like white noise which is great - there is no observable pattern or trend in the data. The ACF plot does not show any spikes above the blue due dashed lines, which also indicates white noise. However the histogram of residuals appears left skewed rather than approximately normally distributed. Let’s look at a Ljung-Box test to confirm that what we’re seeing is in fact white noise residuals.
augment(us_gdp_fit_all) |>
filter(.model == "autoarima110") |>
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 autoarima110 3.81 0.923Here we have a high p-value which strongly validates that the residuals are white noise and this is a good model selection. Despite the slight skewness we see in the histogram of residuals, we can again confirm that the residuals behave like white noise and are not correlated. The ARIMA(1,1,0) model is the best selection.
e. produce forecasts of your fitted model. Do the forecasts look reasonable?
us_gdp_forecast <- us_gdp_fit_all |>
select(autoarima110) |>
forecast(h = "10 years")
us_gdp_forecast |>
autoplot(level = NULL) +
autolayer(us_gdp, GDP_Transformed) +
labs(
title = "Forecasted US GDP (transformed) using ARIMA(1,1,0)",
y = "GDP (Transformed)"
)
These forecasts look reasonable to me.
f. compare the results with what you would obtain using ETS() (with no transformation).
First we’ll fit the ETS model to the original GDP data (no transformation, per the prompt above). This will make it more strange to compare to the transformed ARIMA model, but we’ll move forward regardless.
us_gdp_ets_fit <- us_gdp |>
model(ETS(GDP))Now we’ll produce forecasts also for the next 10 years like we did with ARIMA() except this time using ETS() on the not-transformed (original) data.
us_gdp_ets_forecast <- us_gdp_ets_fit |>
forecast(h = 10)
autoplot(us_gdp_ets_forecast, level = NULL) +
autolayer(us_gdp, GDP, series = "Original Data (Untransformed)") +
labs(
title = "Forecasted US GDP (not transformed) using ETS",
y = "GDP (Not transformed)"
)Warning in geom_line(eval_tidy(expr(aes(!!!aes_spec))), data = object, ..., :
Ignoring unknown parameters: `series`
arima_metrics <- glance(us_gdp_fit_all) |>
select(.model, AICc, BIC)
ets_metrics <- glance(us_gdp_ets_fit) |>
select(.model, AICc, BIC)
bind_rows(arima_metrics, ets_metrics)# A tibble: 6 × 3
.model AICc BIC
<chr> <dbl> <dbl>
1 autoarima110 657. 663.
2 arima011 659. 665.
3 arima111 659. 667.
4 arima210 659. 667.
5 arima012 660. 667.
6 ETS(GDP) 3192. 3201.Visually, the ARIMA(1,1,0) and ETS are “similar” in the sense that it appears the forecast continues the existing trend.
The selected “best” ARIMA model, autoarima110 for ARIMA(1,1,0) has a lowest AICc of 657.1 and a BIC of 662.8. The ETS model without transformation has a AICc of 3191.9 and BIC of 3201.1.
It’s a bit of a challenge to directly compare these two especially because ARIMA is based on a transformed series with differencing while the ETS model is on non-transformed data. These two models are fundamentally different in their approaches.
The AICc and BIC values can’t be compared directly because they’re calculated on different data series (transformed vs original). That said, the values for the ETS model are much higher which overall indicates that we should prefer the ARIMA model in this instance.
While the ETS model could be adjusted to improve performance, it is not as well-suited for this data without preprocessing. Overall, the ARIMA(1,1,0) model on the transformed series appears to be the best fit among the ARIMA models and compared to the ETS model on non-transformed data.