Seasonal ARIMA(0,0,0)(4,1,4)^12 Model

\[ (1-\Phi_1B^{12}-\Phi_2B^{24}-\Phi_3B^{36}-\Phi_4B^{48})(Y_t-Y_{t-12}) = (1-\Theta_1B^{12}-\Theta_2B^{24}-\Theta_3B^{36}-\Theta_4B^{48})a_t \] \[ y_t = y_{t-12} + \Phi_1 (y_{t-12} - y_{t-24}) + \Phi_2 (y_{t-24} - y_{t-36}) + \Phi_3 (y_{t-36} - y_{t-48}) + \Phi_4 (y_{t-48} - y_{t-60}) + a_t - \Theta_1 a_{t-12} - \Theta_2 a_{t-24} - \Theta_3 a_{t-36} - \Theta_4 a_{t-48} \]

set.seed(123)
at <- rnorm(560, mean = 0, sd = 2)
# Seasonal AR coefficients (Phi_1 to Phi_4)
phi_seasonal <- c(0.8, -0.6, 0.4, -0.2) 

# Seasonal MA coefficients (Theta_1 to Theta_4)
theta_seasonal <- c(0.3, -0.1, 0.5, -0.02)
n <- length(at)
y <- numeric(n)
y[1:60] <- 0  # Initialize first 60 lags

for (t in 61:n) {
  
  # Seasonal AR terms (Phi_1 to Phi_4)
  ar_terms <- 
    phi_seasonal[1] * (y[t-12] - y[t-24]) +
    phi_seasonal[2] * (y[t-24] - y[t-36]) +
    phi_seasonal[3] * (y[t-36] - y[t-48]) +
    phi_seasonal[4] * (y[t-48] - y[t-60])
  
  # Seasonal MA terms (Theta_1 to Theta_4)
  ma_terms <- 
    -theta_seasonal[1] * at[t-12] -
    theta_seasonal[2] * at[t-24] -
    theta_seasonal[3] * at[t-36] -
    theta_seasonal[4] * at[t-48]
  
  # Combine all terms
  y[t] <- y[t-12] + ar_terms + at[t] + ma_terms
}

# Discard/nurn in the first 100 observations to avoid initialization artifacts:
y <- y[-(1:60)] # Final simulated series
y_ts <- ts(y, frequency = 12) 
ts.plot(y_ts)

par(mfrow=c(1,2))
acf(y, main = "ACF Original")
pacf(y, main = "PACF Original")

par(mfrow=c(1,2))
# differencing lag = 12
acf(diff(y,12), main = "ACF differencing")
pacf(diff(y,12), main = "PACF differencing")