9 ARIMA and SARIMA Models

Learning Goals

  • Explain and illustrate how first and second order differences remove linear and quadratic trends.
  • Explain and illustrate how seasonal differences remove seasonality.
  • Fit ARIMA and SARIMA models to real data and use model selection tools to decide on a final model.
  • Explain and illustrate the general theoretical procedure for doing forecasting.


Backshift Operators

It can be useful to use the backshift operator, \(B\), to simplify the notation of our models, such that

\[B^k Y_t = Y_{t-k}\text{ for }k = 0,1,2,3,4,...\]

Our \(AR(p)\) can be written as

\[(1 - \phi_1B - \cdots-\phi_pB^p)Y_t = W_t\]

Our \(MA(q)\) can be written as

\[Y_t = (1 + \theta_1B + \cdots+\theta_qB^q)W_t\]

Our Seasonal \(AR(P)\) can be written as

\[(1 - \Phi_1B^S - \cdots-\Phi_PB^{PS})Y_t = W_t\]

Our Seasonal \(MA(Q)\) can be written as

\[Y_t = (1 + \Theta_1B^S + \cdots+\Theta_QB^{QS})W_t\]

All together, an ARMA(p,q) + Seasonal ARMA(P,Q)

\[(1 - \Phi_1B^S - \cdots-\Phi_PB^{PS})(1 - \phi_1B - \cdots-\phi_pB^p)Y_t = (1 + \Theta_1B^S + \cdots+\Theta_QB^{QS})(1 + \theta_1B + \cdots+\theta_qB^q)W_t\]

Another fun fact about the backshift operator:

  • A first order difference to remove a linear trend:

\[Y_t - Y_{t-1} = (1-B)Y_t\]

  • A second order difference to remove a quadratic trend:

\[(Y_t - Y_{t-1}) - (Y_{t-1} - Y_{t-2}) = (1-B)^2Y_t\]

  • A first order seasonal difference to remove seasonality:

\[Y_t - Y_{t-S} = (1-B^S)Y_t\]

  • A combination of first order for trend and seasonality:

\[(1-B)(1-B^S)Y_t = (1 - B - B^{S} + B^{S+1}) Y_t = (Y_t - Y_{t-1}) - ( Y_{t-S}-Y_{t-S-1}) \]

Simulate that Model

Download a template RMarkdown file to start from here.

Below, there is three R code chunks to generate three time series and show the theoretical and estimated ACF.

  1. For each, write out the data generation model in random process notation like above (using \(\phi\)’s and \(\theta\)’s).
  2. For each, compare and contrast the estimated ACF & PACF and the theoretical ACF & PACF.
  3. For each, add R code to fit the data generation model to the data using sarima.
  4. For each, add R code to forecast into the future 6 time units using sarima.for.
library(dplyr)
library(astsa)

y1 <- sarima.sim(ar = c(0.9, -0.1))
plot(y1, type='l')

#estimated ACF and PACF
acf2(y1)

#theoretical ACF and PACF
ARMAacf(ar = c(0.9, -0.1),lag.max = 20) %>% plot(type='h',ylim=c(-1,1))
ARMAacf(ar = c(0.9, -0.1), pacf=TRUE,lag.max = 20) %>% plot(type='h',ylim=c(-1,1))

#fit model to data

#forecast
y2 <- sarima.sim(ma = c(0.7), d = 1, D = 1, S = 6)
plot(y2, type='l')

#estimated ACF and PACF
y2 %>% diff(lag = 1) %>% diff(lag = 6) %>% acf2()

#theoretical ACF and PACF
ARMAacf(ma = 0.7, lag.max = 20) %>% plot(type='h',ylim=c(-1,1))
ARMAacf(ma = 0.7, pacf=TRUE,lag.max = 20) %>% plot(type='h',ylim=c(-1,1))

#fit model to data

#forecast
y3 <- sarima.sim(ar = c(0.9), sar = 0.5, S = 12)
plot(y3, type='l')

#estimated ACF and PACF
acf2(y3)

#theoretical ACF and PACF
ARMAacf(ar = c(0.9, rep(0, 10), 0.5, -0.9*0.5)) %>% plot(type='h',ylim=c(-1,1))
ARMAacf(ar = c(0.9, rep(0, 10), 0.5, -0.9*0.5), pacf=TRUE) %>% plot(type='h',ylim=c(-1,1))

#fit model to data

#forecast