ARIMA models
Kerwin Charles
library(ggplot2)
library(tidyverse)
library(fpp3)#Question 9.1
##Example a
All three plots indicate white noise. In each case, the autocorrelation spikes at all lags fall within the blue dashed critical bounds, meaning none are statistically significant. However the plots look visually different because of sample size. The left plot (n=36) shows more variation in the autocorrelations with some spikes approaching the bounds, the middle plot (n=360) shows smaller and more uniform spikes, and the right plot (n=1,000) shows spikes that are almost completely flat and near zero. This is expected — with more data, the sample autocorrelations converge more closely to their true value of zero.
##Example b
The critical values are calculated as plus or minus 1.96 divided by the square root of n, where n is the sample size. As n increases, this value gets smaller, so the blue bounds narrow. This is why the bounds are wide for n=36, moderate for n=360, and very narrow for n=1,000.
The autocorrelations themselves also differ because they are estimated from the data. With only 36 observations there is more sampling variability, so the estimated autocorrelations fluctuate more around zero even though the true value is zero. With 1,000 observations the estimates are much more stable and precise, clustering tightly around zero. This is simply the effect of larger sample sizes producing more reliable estimates.
#Question 9.2
amazon_stock <- gafa_stock |>
filter(Symbol == "AMZN")
# Plot closing prices
amazon_stock |>
autoplot(Close) +
labs(y = "Closing Price (USD)", x = "Date",
title = "Amazon Daily Closing Stock Price")
# ACF and PACF
amazon_stock |>
ACF(Close) |>
autoplot() +
labs(title = "ACF: Amazon Closing Price")
amazon_stock |>
PACF(Close) |>
autoplot() +
labs(title = "PACF: Amazon Closing Price")
The time plot shows a clear and sustained upward trend in Amazon’s closing stock price over time, with no tendency to revert to a constant mean. This is a hallmark of non-stationarity — the statistical properties of the series, particularly its mean, change over time. The ACF plot reinforces this conclusion. Rather than dropping to zero quickly as would be expected for a stationary series, the autocorrelations decay very slowly and remain large and significant across many lags. The textbook notes that for non-stationary data, the ACF decreases slowly and the value at lag 1 is often large and positive, both of which are evident here. The PACF shows a single dominant spike at lag 1 close to 1, with all subsequent lags insignificant, which is consistent with a random walk process. As the textbook explains, a random walk is a classic non-stationary model where future movements are unpredictable and equally likely to be up or down. Together, these three plots provide strong evidence that the Amazon stock price series is non-stationary and requires first differencing to achieve stationarity.
#Question 9.3
##Example a
turkish_gdp <- global_economy |>
filter(Country == "Turkey")
# Plot original
turkish_gdp |>
autoplot(GDP) +
labs(y = "GDP (USD)", x = "Year", title = "Turkish GDP")
# Find optimal lambda
lambda_turkey <- turkish_gdp |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda_turkey## [1] 0.1572187# Test for differencing
turkish_gdp |>
features(box_cox(GDP, lambda_turkey), unitroot_ndiffs)# Plot transformed and differenced
turkish_gdp |>
autoplot(difference(box_cox(GDP, lambda_turkey))) +
labs(y = "Differenced Box-Cox GDP", title = "Turkish GDP: Transformed and Differenced")
Turkish GDP displays a clear exponential upward trend with increasing variance over time, making it strongly non-stationary. The Guerrero method suggests an optimal Box-Cox lambda of 0.157, which is close to zero and indicates a near-logarithmic transformation is appropriate to stabilise the variance. Following the Box-Cox transformation, the KPSS-based unitroot_ndiffs function suggests one first difference is required to achieve stationarity. The resulting transformed and differenced series shown in the second plot fluctuates around zero with no obvious trend, indicating that the combination of a Box-Cox transformation with lambda of 0.157 and one first difference is sufficient to produce a stationary series.
##Example b
tasmania_acc <- aus_accommodation |>
filter(State == "Tasmania")
# Plot original
tasmania_acc |>
autoplot(Takings) +
labs(y = "Takings (AUD)", x = "Date", title = "Tasmania: Accommodation Takings")
# Find optimal lambda
lambda_tas <- tasmania_acc |>
features(Takings, features = guerrero) |>
pull(lambda_guerrero)
lambda_tas## [1] 0.001819643# Test for seasonal differencing first
tasmania_acc |>
features(box_cox(Takings, lambda_tas), unitroot_nsdiffs)# Test for first differencing after seasonal
tasmania_acc |>
mutate(diff_takings = difference(box_cox(Takings, lambda_tas), 4)) |>
features(diff_takings, unitroot_ndiffs)# Plot transformed and differenced
tasmania_acc |>
autoplot(difference(difference(box_cox(Takings, lambda_tas), 4))) +
labs(y = "Differenced Box-Cox Takings",
title = "Tasmania Accommodation: Transformed and Differenced")
Tasmania’s accommodation takings display both a clear upward trend and strong quarterly seasonality, with the seasonal swings growing larger as the level of takings increases over time — both features of a non-stationary series. The Guerrero method suggests an optimal Box-Cox lambda of approximately 0.002, which is essentially zero and indicates a log transformation is appropriate to stabilise the growing variance. The unitroot_nsdiffs function confirms that one seasonal difference is required, consistent with the strong quarterly pattern visible in the data. The unitroot_ndiffs result after seasonal differencing returned no further differencing requirement, indicating that the seasonal difference alone is sufficient after the Box-Cox transformation. The resulting transformed and differenced series fluctuates around zero with no obvious trend or seasonal pattern remaining, confirming that a log transformation combined with one seasonal difference of lag 4 is sufficient to produce a stationary series.
##Example c
# Plot original
souvenirs |>
autoplot(Sales) +
labs(y = "Sales", x = "Date", title = "Monthly Souvenir Sales")
# Find optimal lambda
lambda_souv <- souvenirs |>
features(Sales, features = guerrero) |>
pull(lambda_guerrero)
lambda_souv## [1] 0.002118221# Test for seasonal differencing first
souvenirs |>
features(box_cox(Sales, lambda_souv), unitroot_nsdiffs)# Test for first differencing after seasonal
souvenirs |>
mutate(diff_sales = difference(box_cox(Sales, lambda_souv), 12)) |>
features(diff_sales, unitroot_ndiffs)# Plot transformed and differenced
souvenirs |>
autoplot(difference(difference(box_cox(Sales, lambda_souv), 12))) +
labs(y = "Differenced Box-Cox Sales",
title = "Souvenir Sales: Transformed and Differenced")
The monthly souvenir sales series displays a strong upward trend with seasonal spikes that grow dramatically in magnitude over time, particularly the December peak which increases from around 20,000 in the late 1980s to over 100,000 by 1994. This combination of trend and increasing seasonal variation confirms the series is non-stationary. The Guerrero method suggests an optimal Box-Cox lambda of approximately 0.002, essentially zero, indicating a log transformation is appropriate to stabilise the growing variance. The unitroot_nsdiffs function confirms one seasonal difference of lag 12 is required to remove the strong monthly seasonality. The unitroot_ndiffs result after seasonal differencing returned an empty tibble, indicating no further first differencing is needed. The resulting transformed and differenced series fluctuates around zero with no obvious trend or growing seasonal pattern, confirming that a log transformation combined with one seasonal difference of lag 12 is sufficient to achieve stationarity.
#Question 9.5
myseries <- aus_retail |>
filter(`Series ID` == "A3349872X")
# Plot original
myseries |>
autoplot(Turnover) +
labs(y = "Turnover ($ Million)", x = "Month",
title = "NSW: Other Retailing Turnover")
# Find optimal lambda
lambda_retail <- myseries |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
lambda_retail## [1] 0.04872512# Test for seasonal differencing first
myseries |>
features(box_cox(Turnover, lambda_retail), unitroot_nsdiffs)# Test for first differencing after seasonal
myseries |>
mutate(diff_turnover = difference(box_cox(Turnover, lambda_retail), 12)) |>
features(diff_turnover, unitroot_ndiffs)# Plot transformed and differenced
myseries |>
autoplot(difference(difference(box_cox(Turnover, lambda_retail), 12))) +
labs(y = "Differenced Box-Cox Turnover",
title = "NSW Retail: Transformed and Differenced")
The NSW Other Retailing Turnover series displays a strong upward trend with seasonal spikes that grow proportionally over time, confirming it is non-stationary. The Guerrero method suggests an optimal Box-Cox lambda of 0.049, which is close to zero and indicates a near-log transformation is appropriate to stabilise the growing variance. The unitroot_nsdiffs function confirms that one seasonal difference of lag 12 is required to remove the strong monthly seasonal pattern. The unitroot_ndiffs result after seasonal differencing returned no further differencing requirement, indicating that the seasonal difference alone is sufficient after the Box-Cox transformation. The resulting transformed and differenced series fluctuates around zero with no obvious trend or seasonal pattern remaining, confirming that a Box-Cox transformation with lambda of 0.049 combined with one seasonal difference of lag 12 is sufficient to produce a stationary series.
#Question 9.6
##Example a
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
sim_ar1 <- tsibble(idx = seq_len(100), y = y, index = idx)##Example b
sim_ar1 |>
autoplot(y) +
labs(title = "AR(1) Model: phi = 0.6", y = "y", x = "Time")
The AR(1) series with phi of 0.6 fluctuates around zero with moderate persistence and consecutive values tend to move in the same direction before reverting to the mean. As phi increases toward 1 the series becomes more persistent with longer runs above or below zero, resembling a random walk. As phi decreases toward 0 the series becomes more erratic and resembles white noise with less carry-over from one period to the next.
##Example c
y_ma <- numeric(100)
e_ma <- rnorm(100)
for(i in 2:100)
y_ma[i] <- 0.6*e_ma[i-1] + e_ma[i]
sim_ma1 <- tsibble(idx = seq_len(100), y = y_ma, index = idx)##Example d
sim_ma1 |>
autoplot(y) +
labs(title = "MA(1) Model: theta = 0.6", y = "y", x = "Time")
The MA(1) series also fluctuates around zero but with less persistence than the AR(1) model. The current value depends only on the current and one previous error term, so the series reverts to zero more quickly. As theta increases the series shows slightly smoother short-term movements, while as theta decreases toward 0 the series approaches pure white noise.
##Example e
y_arma <- numeric(100)
e_arma <- rnorm(100)
for(i in 2:100)
y_arma[i] <- 0.6*y_arma[i-1] + 0.6*e_arma[i-1] + e_arma[i]
sim_arma <- tsibble(idx = seq_len(100), y = y_arma, index = idx)##Example f
y_ar2 <- numeric(100)
e_ar2 <- rnorm(100)
for(i in 3:100)
y_ar2[i] <- -0.8*y_ar2[i-1] + 0.3*y_ar2[i-2] + e_ar2[i]
sim_ar2 <- tsibble(idx = seq_len(100), y = y_ar2, index = idx)##Example g
sim_arma |>
autoplot(y) +
labs(title = "ARMA(1,1): phi = 0.6, theta = 0.6", y = "y", x = "Time")
sim_ar2 |>
autoplot(y) +
labs(title = "AR(2): phi1 = -0.8, phi2 = 0.3", y = "y", x = "Time")
The ARMA(1,1) series fluctuates around zero with moderate persistence, combining the autocorrelation structure of both an AR and MA component. It behaves similarly to the AR(1) series but with somewhat larger swings due to the added MA component. The AR(2) series with phi1 of -0.8 and phi2 of 0.3 behaves dramatically differently — it exhibits explosive oscillating behaviour with the amplitude growing rapidly over time, confirming the non-stationary nature of these parameter values. The oscillations are caused by the negative phi1 which forces alternating signs, while the overall magnitude grows unboundedly since the parameters fall outside the stationarity region.
#Question 9.7
##Example a
fit_air <- aus_airpassengers |>
model(ARIMA(Passengers ~ pdq(d=1)))
report(fit_air)## 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# Base R residual diagnostics
resid_vec <- augment(fit_air) |> pull(.innov)
resid_vec <- resid_vec[!is.na(resid_vec)]
par(mfrow = c(1,3))
plot(resid_vec, type = "l", main = "Residuals", ylab = "Residuals")
acf(resid_vec, main = "ACF of Residuals")
hist(resid_vec, main = "Histogram", xlab = "Residuals")
par(mfrow = c(1,1))
# Ljung-Box test
Box.test(resid_vec, lag = 10, type = "Ljung-Box")##
## Box-Ljung test
##
## data: resid_vec
## X-squared = 6.7718, df = 10, p-value = 0.7468# Forecast
fit_air |>
forecast(h = 10) |>
autoplot(aus_airpassengers) +
labs(y = "Passengers (millions)", x = "Year",
title = "Australian Air Passengers: ARIMA(0,1,0) w/ drift Forecast")
The ARIMA() function with d=1 constrained selected an ARIMA(0,1,0) model with drift, which is equivalent to a random walk with drift. The drift constant of 1.419 indicates that passenger numbers increase by approximately 1.42 million per year on average. The residual plots show mostly encouraging signs — the ACF has no significant spikes beyond the blue bounds, and the histogram is approximately bell-shaped, though slightly right-skewed due to the large spike visible around index 40 in the residuals plot. The Ljung-Box test returns a p-value of 0.747, which is well above 0.05, confirming that the residuals are consistent with white noise. The forecast plot shows a linear upward trend continuing from the last observed value, with prediction intervals widening over the 10-year horizon as expected for a model with d=1.
##Example b
The ARIMA(0,1,0) model with drift can be written in backshift notation as: (1 - B)y_t = c + epsilon_t where c = 1.419 is the drift constant and epsilon_t is white noise. This simplifies to y_t = c + y_{t-1} + epsilon_t, confirming the random walk with drift structure.
##Example c
fit_drift <- aus_airpassengers |>
model(ARIMA(Passengers ~ 1 + pdq(0,1,0)))
fit_drift |>
forecast(h = 10) |>
autoplot(aus_airpassengers) +
labs(y = "Passengers (millions)", x = "Year",
title = "ARIMA(0,1,0) with Drift")
The ARIMA(0,1,0) model with drift specified manually in part c produces identical forecasts to the model automatically selected in part a. This is because the ARIMA() function with d=1 constrained already selected an ARIMA(0,1,0) with drift as the best model. Both models are random walks with drift, projecting a linear upward trend of approximately 1.42 million passengers per year with widening prediction intervals over the forecast horizon.
##Example d
# ARIMA(2,1,2) with drift
fit_212 <- aus_airpassengers |>
model(ARIMA(Passengers ~ 1 + pdq(2,1,2)))
fit_212 |>
forecast(h = 10) |>
autoplot(aus_airpassengers) +
labs(y = "Passengers (millions)", x = "Year",
title = "ARIMA(2,1,2) with Drift")
# Remove constant
fit_212_nc <- aus_airpassengers |>
model(ARIMA(Passengers ~ 0 + pdq(2,1,2)))
fit_212_nc |>
forecast(h = 10) |>
autoplot(aus_airpassengers) +
labs(y = "Passengers (millions)", x = "Year",
title = "ARIMA(2,1,2) without Constant")
The ARIMA(2,1,2) model with drift produces forecasts very similar to the ARIMA(0,1,0) with drift from parts a and c, projecting a continued linear upward trend with prediction intervals that are noticeably narrower. This is because the additional AR and MA parameters in the ARIMA(2,1,2) model capture more of the short-term dynamics in the data, reducing the residual variance and tightening the prediction intervals. The point forecasts are comparable across all three models, suggesting the drift term is the dominant driver of the forecast trajectory.
Removing the constant from the ARIMA(2,1,2) model has a dramatic effect — the forecast disappears entirely and only the historical data is shown. As the textbook explains, when c=0 and d=1, the long-term forecasts will go to a non-zero constant, meaning the model projects that passenger numbers will flatten out and stop growing. Without the drift constant there is no mechanism to continue the upward trend, so the forecasts essentially become flat at the last observed value. This confirms that the drift constant is essential for this series to produce meaningful long-run forecasts that reflect the sustained growth in Australian air travel.
##Example e
fit_021 <- aus_airpassengers |>
model(ARIMA(Passengers ~ 1 + pdq(0,2,1)))
fit_021 |>
forecast(h = 10) |>
autoplot(aus_airpassengers) +
labs(y = "Passengers (millions)", x = "Year",
title = "ARIMA(0,2,1) with Constant")
The ARIMA(0,2,1) model with a constant produces notably different forecasts compared to the previous models. As the textbook explains, when c is not zero and d=2, the long-term forecasts will follow a quadratic trend. This is clearly visible in the plot — the forecasts curve upward more steeply than the linear projections from the ARIMA(0,1,0) and ARIMA(2,1,2) models, with passenger numbers projected to reach around 100 million by the end of the forecast horizon compared to around 85 million for the d=1 models. The prediction intervals also widen more rapidly, reflecting the greater uncertainty that comes with double differencing. The textbook specifically notes that including a constant with d=2 induces a quadratic trend in forecasts and is generally not recommended, as it tends to produce overly optimistic long-run projections. In this case the quadratic growth assumption is likely unrealistic for air passenger numbers over the long run.
#Question 9.8
##Example a
us_gdp <- global_economy |>
filter(Country == "United States")
us_gdp |>
autoplot(GDP) +
labs(y = "GDP (USD)", x = "Year", title = "United States GDP")
lambda_us <- us_gdp |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda_us## [1] 0.2819443The United States GDP series displays a clear exponential upward trend with increasing variance over time, confirming that a transformation is necessary before fitting an ARIMA model. The Guerrero method suggests an optimal Box-Cox lambda of 0.282, which is between 0 and 0.5, indicating a transformation somewhere between a log and a square root is appropriate to stabilise the growing variance in the series.
##Example b
fit_us <- us_gdp |>
model(ARIMA(box_cox(GDP, lambda_us) ~ pdq(d=1)))
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.78The ARIMA() function with d=1 constrained selected an ARIMA(1,1,0) model with drift applied to the Box-Cox transformed GDP series. The AR(1) coefficient of 0.459 indicates moderate positive autocorrelation in the differenced series, meaning that above-average growth in one year tends to be followed by above-average growth in the next. The drift constant of 118.18 captures the underlying upward trend in the transformed series. The model can be written as:
(1 - 0.459B)(1 - B)y*_t = 118.18 + epsilon_t
where y*_t denotes the Box-Cox transformed GDP series.
##Example c
fit_us_models <- us_gdp |>
model(
auto = ARIMA(box_cox(GDP, lambda_us)),
arima011 = ARIMA(box_cox(GDP, lambda_us) ~ pdq(0,1,1)),
arima110 = ARIMA(box_cox(GDP, lambda_us) ~ pdq(1,1,0)),
arima111 = ARIMA(box_cox(GDP, lambda_us) ~ pdq(1,1,1)),
arima012 = ARIMA(box_cox(GDP, lambda_us) ~ pdq(0,1,2))
)
glance(fit_us_models) |> arrange(AICc) |> select(.model, AICc, BIC)Several alternative ARIMA models were fitted to the Box-Cox transformed US GDP series and compared using AICc. The results show that the ARIMA(1,1,0) model with drift has the lowest AICc of 657.10 and BIC of 662.78, confirming it as the best model among those considered. The ARIMA(0,1,1) model is a close second with an AICc of 659.18, followed by ARIMA(1,1,1) and ARIMA(0,1,2) which both have higher AICc values. The additional parameters in the higher order models do not improve the fit sufficiently to justify their inclusion, so the ARIMA(1,1,0) with drift remains the preferred model on the basis of parsimony and lowest AICc.
##Example d
us_gdp_bc <- global_economy |>
filter(Country == "United States") |>
mutate(GDP_bc = box_cox(GDP, lambda_us))
best_fit <- us_gdp_bc |>
model(ARIMA(GDP_bc ~ pdq(1,1,0)))
resid_vec <- augment(best_fit) |>
filter(!is.na(.innov)) |>
pull(.innov)
par(mfrow = c(1,3))
plot(resid_vec, type = "l", main = "Residuals", ylab = "Residuals")
acf(resid_vec, main = "ACF of Residuals")
hist(resid_vec, main = "Histogram", xlab = "Residuals")
par(mfrow = c(1,1))
Box.test(resid_vec, lag = 10, type = "Ljung-Box")##
## Box-Ljung test
##
## data: resid_vec
## X-squared = 3.8137, df = 10, p-value = 0.9554The residual diagnostics for the ARIMA(1,1,0) model look good. The residuals plot shows fluctuation around zero with no obvious trend or pattern, though there is a notable large negative residual around the 2008 financial crisis period which is an expected structural shock rather than a model deficiency. The ACF plot shows all spikes well within the blue bounds with no significant autocorrelation remaining. The histogram is approximately bell-shaped though slightly left-skewed due to the 2008 outlier. The Ljung-Box test returns a p-value of 0.955, which is well above 0.05, confirming that the residuals are consistent with white noise. The ARIMA(1,1,0) model with drift applied to the Box-Cox transformed series passes all residual checks and is suitable for forecasting.
##Example e
best_fit |>
forecast(h = 10) |>
autoplot(us_gdp_bc) +
labs(y = "Box-Cox GDP", x = "Year",
title = "United States GDP: ARIMA(1,1,0) Forecast")
The forecasts from the ARIMA(1,1,0) model with drift look reasonable. The model projects a continued upward trend in the Box-Cox transformed GDP series, consistent with the long-run growth pattern observed in the historical data. The prediction intervals are relatively narrow, reflecting the good fit of the model and the low residual variance. The forecasts connect smoothly from the last observed value and the trajectory is plausible given the sustained growth trend in US GDP over the entire sample period. The slight dip around 2008 in the historical data does not appear to have distorted the forecasts, which correctly extrapolate the underlying long-run growth trend rather than the short-term fluctuation.
##Example f
fit_ets <- us_gdp |>
model(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.089bind_rows(
best_fit |> accuracy(),
fit_ets |> accuracy()
) |> select(.model, RMSE, MAE)fit_ets |>
forecast(h = 10) |>
autoplot(us_gdp) +
labs(y = "GDP (USD)", x = "Year",
title = "United States GDP: ETS Forecast")
The ETS() function without transformation selected an ETS(M,A,N) model — a Holt linear trend model with multiplicative errors. The forecast plot shows a continued upward trend but with dramatically wider prediction intervals compared to the ARIMA model, reflecting the multiplicative error structure which scales uncertainty proportionally with the level of GDP. The 95% prediction interval fans out very widely, reaching nearly 4e+13 USD at the upper bound by the end of the forecast horizon, which seems overly uncertain.
Comparing the two approaches directly is not straightforward because the ARIMA model was fitted to the Box-Cox transformed series while ETS was fitted to the original untransformed GDP. This is reflected in the RMSE values — the ARIMA RMSE of 72.08 is on the transformed scale while the ETS RMSE of 1.67e+11 is on the original USD scale, so they cannot be directly compared. However, visually the ARIMA model produces tighter and more conservative prediction intervals on the transformed scale, while the ETS model produces very wide intervals on the original scale due to the exponential growth in GDP. The ARIMA model with Box-Cox transformation is arguably more appropriate here as it stabilizes the variance and produces more reliable prediction intervals.