ARIMA(0,0,0)(4,0,0)^12 model: \[y_t=\Phi_1y_{t-12}+\Phi_2y_{t-24}+\Phi_3y_{t-36}+\Phi_4y_{t-48}+a_t\]

set.seed(1234)
n <- 1000
sigma <- 2
Phi <- c(0.3, 0.25, 0.15, 0.1)
at <- rnorm(n, mean = 0, sd = sigma)

y <- numeric(n)
y[1:48] <- rnorm(48, mean = 0, sd = sigma)

for (t in 49:n) {
  y[t] <- Phi[1] * y[t - 12] + 
    Phi[2] * y[t - 24] + 
    Phi[3] * y[t - 36] + 
    Phi[4] * y[t - 48] + 
    at[t]
}
ts.plot(y)

par(mfrow=c(1,2))
acf(y, lag.max = 144)
pacf(y, lag.max = 144)