# Setup
set.seed(123)

# Simulasi Data SARIMA
n <- 300
at <- rnorm(n + 60, 0, 2)

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

y <- numeric(n + 60) #lag 48 + buffer 12 karena adanya seasonal differencing
y[1:60] <- 0 #Inisialisasi

for (t in 61:(n + 60)) {
  y[t] <- y[t - 12] +
    phi1 * (y[t - 12] - y[t - 24]) +
    phi2 * (y[t - 24] - y[t - 36]) +
    phi3 * (y[t - 36] - y[t - 48]) +
    phi4 * (y[t - 48] - y[t - 60]) +
    at[t]
}

#Burn-in / menghapus buffer awal
y <- y[-(1:60)]

y <- ts(y, frequency = 12)

# Seasonal differencing lag 12
z <- diff(y, lag = 12)

par(mfrow = c(1,1))

# Plot y (before differencing)
plot(y, main = "Data Yt (Sebelum Differencing)", col = "blue", lwd = 2)

# Plot z (after seasonal differencing)
plot(z, main = "Data Setelah Seasonal Differencing ∇12 Yt", col = "red", lwd = 2)

# ACF dan PACF y (before differencing)
acf(y, lag.max = 100, main = "ACF Yt (Before Differencing)")

pacf(y, lag.max = 100, main = "PACF Yt (Before Differencing)")

# ACF dan PACF z (after seasonal differencing)
acf(z, lag.max = 100, main = "ACF Zt (After Seasonal Differencing)")

pacf(z, lag.max = 100, main = "PACF Zt (After Seasonal Differencing)")