Module 8 discussion

Part I: ARIMA Model Selection and Forecasting

Data: Annual Flow of the Nile River (1871-1970)

data("Nile")
autoplot(Nile) +
  ggtitle("Annual Flow of the Nile River (1871-1970)") +
  xlab("Year") + ylab("Flow (10^8 cubic metres)") +
  theme_minimal()

Stationarity Check

adf.test(Nile)

    Augmented Dickey-Fuller Test

data:  Nile
Dickey-Fuller = -3.3657, Lag order = 4, p-value = 0.0642
alternative hypothesis: stationary

The ADF test result guides whether differencing is needed. If p > 0.05 the series is non-stationary and we apply first differencing (d = 1). If p <= 0.05 the series is already stationary (d = 0).

ts_diff <- diff(Nile)
adf.test(ts_diff)

    Augmented Dickey-Fuller Test

data:  ts_diff
Dickey-Fuller = -6.5924, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

After first differencing the series becomes stationary, confirming d = 1.

ACF and PACF Analysis

par(mfrow = c(1, 2))
acf(ts_diff, main = "ACF - First Differenced Nile", lag.max = 20)
pacf(ts_diff, main = "PACF - First Differenced Nile", lag.max = 20)

par(mfrow = c(1, 1))

Manual parameter justification:

• ACF cuts off after lag 1 → suggests MA(1), so q = 1
• PACF cuts off after lag 1 → suggests AR(1), so p = 1
• d = 1 confirmed from differencing above

Manually selected: ARIMA(1,1,1)

Fit Manual vs Auto ARIMA

manual_fit <- arima(Nile, order = c(1, 1, 1))
summary(manual_fit)

Call:
arima(x = Nile, order = c(1, 1, 1))

Coefficients:
         ar1      ma1
      0.2544  -0.8741
s.e.  0.1194   0.0605

sigma^2 estimated as 19769:  log likelihood = -630.63,  aic = 1267.25

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE     MASE
Training set -16.06603 139.8986 109.9998 -4.005967 12.78745 0.825499
                    ACF1
Training set -0.03228482

auto_fit <- auto.arima(Nile, seasonal = FALSE)
summary(auto_fit)

Series: Nile 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.2544  -0.8741
s.e.  0.1194   0.0605

sigma^2 = 20177:  log likelihood = -630.63
AIC=1267.25   AICc=1267.51   BIC=1275.04

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE     MASE
Training set -16.06603 139.8986 109.9998 -4.005967 12.78745 0.825499
                    ACF1
Training set -0.03228482

Reconciliation: The manually selected ARIMA(1,1,1) is compared against auto.arima() which minimizes AIC/BIC. Both approaches agree on d=1. Any difference in p or q reflects the algorithm optimizing information criteria versus visual ACF/PACF reading — both are valid and typically yield close results for this dataset.

Forecast — Next 10 Years

fc <- forecast(manual_fit, h = 10, level = 95)

autoplot(fc) +
  ggtitle("ARIMA(1,1,1) Forecast — Nile River Flow, Next 10 Years") +
  xlab("Year") + ylab("Flow (10^8 cubic metres)") +
  theme_minimal()

---

Part II: Seasonal ARIMA

Data: Monthly Gas Production in Australia

data("gas")
autoplot(gas) +
  ggtitle("Monthly Gas Production in Australia") +
  xlab("Year") + ylab("Gas Production") +
  theme_minimal()

Stationarity and Differencing

ts_log <- log(gas)
adf.test(ts_log)

    Augmented Dickey-Fuller Test

data:  ts_log
Dickey-Fuller = -0.515, Lag order = 7, p-value = 0.9811
alternative hypothesis: stationary

ts_sdiff <- diff(diff(ts_log, lag = 12))
adf.test(ts_sdiff)

    Augmented Dickey-Fuller Test

data:  ts_sdiff
Dickey-Fuller = -7.7321, Lag order = 7, p-value = 0.01
alternative hypothesis: stationary

ACF and PACF after Differencing

par(mfrow = c(1, 2))
acf(ts_sdiff, lag.max = 36, main = "ACF after Seasonal + Regular Differencing")
pacf(ts_sdiff, lag.max = 36, main = "PACF after Seasonal + Regular Differencing")

par(mfrow = c(1, 1))

Fit Seasonal ARIMA

sarima_fit <- auto.arima(ts_log, seasonal = TRUE, stepwise = FALSE)
summary(sarima_fit)

Series: ts_log 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.3304  -0.8256
s.e.   0.0462   0.0473

sigma^2 = 0.002553:  log likelihood = 719.45
AIC=-1432.91   AICc=-1432.86   BIC=-1420.49

Training set error measures:
                       ME       RMSE        MAE         MPE      MAPE     MASE
Training set 0.0001938054 0.04972152 0.03559847 0.006812281 0.3798947 0.357045
                     ACF1
Training set -0.003186571

Parameter Interpretation

The fitted model is typically ARIMA(0,1,1)(0,1,1)[12]:

• d = 1: one regular difference removes the upward trend in gas production
• D = 1: one seasonal difference at lag 12 removes the repeating annual pattern
• MA(1) (q = 1): non-seasonal moving average — the coefficient θ₁ weights the previous forecast error to correct short-term autocorrelation. A value between -0.5 and -0.8 is typical and indicates effective error correction
• SMA(1) (Q = 1): seasonal moving average at lag 12 — the coefficient Θ₁ corrects for year-over-year forecast errors, capturing how production in one month relates to the same month the previous year
• [12]: seasonal period is 12 months

The negative signs on both MA coefficients are expected after differencing and indicate the model is appropriately correcting for over-differencing artifacts.

fc_sarima <- forecast(sarima_fit, h = 24, level = 95)
autoplot(fc_sarima) +
  ggtitle("Seasonal ARIMA Forecast — Australian Gas Production (Log Scale)") +
  xlab("Year") + ylab("Log Gas Production") +
  theme_minimal()

---

Part III: White Noise, Random Walk, and Random Walk with Drift

Mathematical Definitions and Equations

White Noise (WN): \[\epsilon_t \sim WN(0, \sigma^2)\] where \(\epsilon_t\) is independently and identically distributed with mean 0 and constant variance \(\sigma^2\). Each observation is purely random with no relationship to past values.

Random Walk (RW): \[Y_t = Y_{t-1} + \epsilon_t, \quad \epsilon_t \sim WN(0, \sigma^2)\] where \(Y_0 = 0\) by convention. Each value equals the previous value plus a random shock. Shocks accumulate permanently with no tendency to revert to a mean.

Random Walk with Drift (RW + Drift): \[Y_t = \mu + Y_{t-1} + \epsilon_t, \quad \epsilon_t \sim WN(0, \sigma^2)\] where \(\mu\) is the drift parameter. Adds a systematic directional trend to the random walk. When \(\mu > 0\) the series trends upward; when \(\mu < 0\) it trends downward.

Relationship via Differencing

Taking the first difference \(\Delta Y_t = Y_t - Y_{t-1}\):

Differencing RW gives WN: \[\Delta Y_t = Y_t - Y_{t-1} = \epsilon_t \sim WN(0, \sigma^2) \checkmark\]

Differencing RW with drift gives WN + constant: \[\Delta Y_t = Y_t - Y_{t-1} = \mu + \epsilon_t \checkmark\]

This confirms all three processes are related through differencing. RW and RWD are integrated processes of order 1 — differencing once reduces them to stationarity.

Simulation and Plots

set.seed(42)
n <- 200
mu <- 0.3
eps <- rnorm(n, 0, 1)

wn  <- eps
rw  <- cumsum(eps)
rwd <- cumsum(mu + eps)

df <- data.frame(
  time    = rep(1:n, 3),
  value   = c(wn, rw, rwd),
  process = factor(rep(c("White Noise", "Random Walk", "RW with Drift"), each = n),
                   levels = c("White Noise", "Random Walk", "RW with Drift"))
)

ggplot(df, aes(x = time, y = value, color = process)) +
  geom_line() +
  facet_wrap(~process, scales = "free_y", ncol = 1) +
  theme_minimal() +
  ggtitle("White Noise, Random Walk, and Random Walk with Drift") +
  xlab("Time") + ylab("Value") +
  theme(legend.position = "none")

Comparison of the Three Processes

Property	White Noise	Random Walk	RW with Drift
Trend	None	None (wanders randomly)	Upward (μ>0) or Downward (μ<0)
Shock	Temporary, dissipates immediately	Permanent, accumulates forever	Permanent, accumulates with trend
Memory	None — each observation independent	Infinite — all past shocks persist	Infinite — all past shocks persist
Stationarity	Yes — constant mean and variance	No — variance grows over time	No — mean and variance both grow
Forecast	Constant at mean = 0	Flat at last observed value	Trending linearly by μ per period

Visual description:

• White Noise oscillates randomly around zero with no pattern — looks like static
• Random Walk wanders up and down without direction — drifts away from zero but has no trend
• RW with Drift wanders like a random walk but with a clear upward slope driven by μ = 0.3