Model: ARIMA(4,0,0)(4,0,0)¹² \[ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \phi_4 y_{t-4} + \Phi_1 y_{t-12} + \Phi_2 y_{t-24} + \Phi_3 y_{t-36} + \Phi_4 y_{t-48} + \epsilon_t \]

set.seed(123)

n <- 1000
at <- rnorm(n, mean = 0, sd = 2)

phi1 <- 0.5
phi2 <- -0.3
phi3 <- 0.2
phi4 <- -0.1

Phi1 <- 0.4
Phi2 <- -0.3
Phi3 <- -0.2
Phi4 <- 0.01

s_lag <- 12
y <- numeric(n)

for (t in (max(4, 4*s_lag) + 1):(n)) {
  ar_nonseasonal <- phi1*y[t-1] + phi2*y[t-2] + phi3*y[t-3] + phi4*y[t-4]
  ar_seasonal <- Phi1*y[t - s_lag] + Phi2*y[t - 2*s_lag] + Phi3*y[t - 3*s_lag] + Phi4*y[t - 4*s_lag]
  
  y[t] <- ar_nonseasonal + ar_seasonal + at[t]
}

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

Time-Series Plot

ts.plot(y)

ACF and PACF

par(mfrow=c(1,2))
acf(y)
pacf(y)