SARIMA(4,0,0)(4,0,0)^12

set.seed(123)
a_t <- rnorm(1000, mean = 0, sd = 1)

ϕ_1 <- 0.7
ϕ_2 <- -0.5
ϕ_3 <- 0.3
ϕ_4 <- -0.05

Φ_1 <- 0.4
Φ_2 <- 0.1
Φ_3 <- 0.2
Φ_4 <- -0.05

n <- length(a_t)
z_ <- numeric(n)

# Initialize first fifty two values of z (e.g., zero or any arbitrary values)
z_[1:52] <- 0  

# Simulasi SARIMA(4,0,0)(4,0,0)[12]
for (t in 53:n) {
  z_[t] <- ϕ_1*z_[t-1] + ϕ_2*z_[t-2] + ϕ_3*z_[t-3] + 
    ϕ_4*z_[t-4] + Φ_1*z_[t-12] + Φ_2*z_[t-24] + 
    Φ_3*z_[t-36] + Φ_4*z_[t-48] - ϕ_1*Φ_1*z_[t-13] -
    ϕ_1*Φ_2*z_[t-25] - ϕ_1*Φ_3*z_[t-37] - ϕ_1*Φ_4*z_[t-49] -
    ϕ_2*Φ_1*z_[t-14] - ϕ_2*Φ_2*z_[t-26] - ϕ_2*Φ_3*z_[t-38] -
    ϕ_2*Φ_4*z_[t-50] - ϕ_3*Φ_1*z_[t-15] - ϕ_3*Φ_2*z_[t-27] -
    ϕ_3*Φ_3*z_[t-39] - ϕ_3*Φ_4*z_[t-51] - ϕ_4*Φ_1*z_[t-16] -
    ϕ_4*Φ_2*z_[t-28] - ϕ_4*Φ_3*z_[t-40] - ϕ_4*Φ_4*z_[t-52] +
    a_t[t]
}

# Hilangkan burn-in 
z_ <- z_[-(1:100)]

# Plot time series
ts.plot(z_, main = "Simulasi SARIMA(4,0,0)(4,0,0)[12]", ylab = "Z", xlab = "Time")

# ACF dan PACF
par(mfrow=c(1,2))
acf(z_, main = "ACF", lag.max = 60)
pacf(z_, main = "PACF",lag.max = 60)