\[ Φ_4(B^{12})ϕ_4(B)y_t​=Θ_4(B^{12})θ_4(B)a_t \]

\[ y_t=0.2y_{t-1}-0.2y_{t-2}+0.1y_{t-3}+0.2y_{t-4}+0.2y_{t-12}+0.1y_{t-24}-0.2y_{t-36}+0.3y_{t-48}\] \[+0.3a_{t-1}+0.2a_{t-2}-0.2a_{t-3}+0.1a_{t-4}+0.3a_{t-12}+0.2a_{t-24}-0.1a_{t-36}+0.2a_{t-48} \]

set.seed(123) 
at <- rnorm(1000, mean = 0, sd = 2)
#parameter non-seasonal (4,0,4) 
phi1 <- 0.2 
phi2 <- -0.2 
phi3 <- 0.1 
phi4 <- 0.2 
theta1 <- 0.3 
theta2 <- 0.2 
theta3 <- -0.2 
theta4 <- 0.1

#parameter seasonal (4,0,4)^12 
Phi1 <- 0.2 
Phi2 <- 0.1 
Phi3 <- -0.2 
Phi4 <- 0.3 
Theta1 <- 0.3 
Theta2 <- 0.2 
Theta3 <- -0.1 
Theta4 <- 0.2
n <- length(at); n 
## [1] 1000
y <- numeric(n)
y[1:48] <- 0 

for (t in 49:n) { 

  ar <- phi1*(y[t-1]) + phi2*(y[t-2]) + phi3*(y[t-3]) + phi4*(y[t-4]) 
  ma <- theta1*at[t-1] + theta2*at[t-2] + theta3*at[t-3] + theta4*at[t-4] 
  sar <- Phi1*(y[t-12]) + Phi2*(y[t-24]) + Phi3*(y[t-36]) + Phi4*(y[t-48]) 
  sma <- Theta1*at[t-12] + Theta2*at[t-24] + Theta3*at[t-36] + Theta4*at[t-48] 

  y[t] <- ar + ma + sar + sma + at[t] #model 

  }

y <- y[-(1:100)]
plot.ts(y, main = "Simulasi ARIMA(4,0,4)(4,0,4)[12]")

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