arima (4,1,0)(0,1,4)[12]

set.seed(12)
# Parameter AR dan MA seasonal
phi1 = -0.5; phi2 = 0.7; phi3 = 0.1; phi4 = -0.4
Theta1 = -0.1; Theta2 = 0.3; Theta3 = 0.7; Theta4 = -0.8

n <- 1000
at <- rnorm(n + 50, mean = 0, sd = 2)  # tambah buffer untuk seasonal lag
y <- numeric(n + 50)

# Iterasi dimulai dari t = 49 agar tidak ada index negatif
for (t in 49:(n + 48)) {
  dif = -y[t-1] - y[t-12] + y[t-13]
  
  ar <- phi1*(y[t-1] - y[t-2] - y[t-13] + y[t-14]) +
    phi2*(y[t-2] - y[t-3] - y[t-14] + y[t-15]) +
    phi3*(y[t-3] - y[t-4] - y[t-15] + y[t-16]) +
    phi4*(y[t-4] - y[t-5] - y[t-16] + y[t-17])
  
  ma <- at[t] - Theta1*at[t-12] - Theta2*at[t-24] - Theta3*at[t-36] - Theta4*at[t-48]
  
  y[t] <- ar + ma - dif
}

# Ambil hanya bagian simulasi (buang 48 awal)
simulasi <- y[49:(n+48)]

# Plot tanpa differencing
ts.plot(simulasi, main = "Simulasi ARIMA(4,1,0)(0,1,4)[12]", ylab = "Y", xlab = "Waktu")

# Differencing biasa (non-musiman & musiman)
diff1 <- diff(simulasi, differences = 1)
diff2 <- diff(diff1, lag = 12, differences = 1)

# Plot dengan differencing
ts.plot(diff2, main = "Data setelah differencing biasa dan musiman", ylab = "Y", xlab = "Waktu")

# ACF Sebelum Differencing
acf(simulasi,lag.max=60)

pacf(simulasi,lag.max=60)

acf(diff2,lag.max=60)

pacf(diff2,lag.max=60)