SARIMA Model Simulation

\[\begin{equation*} \begin{split} &(1 - \Phi_1 B^{12} - \Phi_2 B^{24} - \Phi_3 B^{36} - \Phi_4 B^{48}) (1 - B)(1 - B^{12}) y_t = a_t \\ &y_t = y_{t-1} + y_{t-12} - y_{t-13} \\ &\quad + \Phi_1 \left(y_{t-12} - y_{t-13} - y_{t-24} + y_{t-25}\right) \\ &\quad + \Phi_2 \left(y_{t-24} - y_{t-25} - y_{t-36} + y_{t-37}\right) \\ &\quad + \Phi_3 \left(y_{t-36} - y_{t-37} - y_{t-48} + y_{t-49}\right) \\ &\quad + \Phi_4 \left(y_{t-48} - y_{t-49} - y_{t-60} + y_{t-61}\right) + a_t. \end{split} \end{equation*}\]

library(forecast)
library(TSA)

# Seasonal ARIMA(0,1,0)(4,1,0)12
# Generating white noise residual
set.seed(1)
at     <- rnorm(1000, mean = 0, sd = 2)
PHI1   <- 0.7
PHI2   <- -0.8
PHI3   <- 0.65
PHI4   <- -0.2

n <- length(at)
y <- numeric(n)
# Initialize first sixty one values of y (e.g., zero or any arbitrary values)
y[1:61] <- 0

for (t in 62:n) {
  y[t] <- y[t-1] + y[t-12] - y[t-13] +
    PHI1 * (y[t-12] - y[t-13] - y[t-24] + y[t-25]) +
    PHI2 * (y[t-24] - y[t-25] - y[t-36] + y[t-37]) +
    PHI3 * (y[t-36] - y[t-37] - y[t-48] + y[t-49]) +
    PHI4 * (y[t-48] - y[t-49] - y[t-60] + y[t-61]) +
    at[t]
}

y <- y[-(1:100)]

ts.plot(y)

# ACF PACF of The Non-Differenced Simulation Data
par(mfrow=c(1,2))
Acf(y, lag.max=100)
Pacf(y, lag.max=100)

# ACF PACF of The Non-Seasonally Differenced Simulation Data
diff_y <- diff(y)
par(mfrow=c(1,2))
Acf(diff_y, lag.max=100)
Pacf(diff_y, lag.max=100)

# ACF PACF of The Non-Seasonally and Seasonally Differenced Simulation Data
diff_y2 <- diff(diff(y), lag=12)
par(mfrow=c(1,2))
Acf(diff_y2, lag.max=100)
Pacf(diff_y2, lag.max=100)