SARIMA Model Simulation

ARIMA(3,1,3)(3,1,3)^12 Model

\[ (1-\phi_1(B)-\phi_2(B^{2})-\phi_3(B^{3})) (1-\Phi_1(B^{12})-\Phi_2(B^{24})-\Phi_3(B^{36}))(1 - B)(1 - B^{12})y_t = (1+\theta_1(B)+\theta_2(B^{2})+\theta_3(B^{3})) (1+\Theta_1(B^{12})+\Theta_2(B^{24})+\Theta_3(B^{36}))a_t\\ \]

\[ y_t - y_{t-1} - y_{t-12} + y_{t-13} = a_t + \theta_1 a_{t-1} + \theta_2 a_{t-2} + \theta_3 a_{t-3} + \Theta_1 a_{t-12} + \Theta_2 a_{t-24} + \Theta_3 a_{t-36} \]

\[ + \phi_1 (y_{t-1}- y_{t-2} - y_{t-13} + y_{t-14})+ \phi_2 (y_{t-2}- y_{t-3} - y_{t-14} + y_{t-15}) + \phi_3 (y_{t-3}- y_{t-4} - y_{t-15} + y_{t-16}) \\ + \Phi_1 (y_{t-12}- y_{t-13} - y_{t-24} + y_{t-25})+ \Phi_2 (y_{t-24}- y_{t-25} - y_{t-36} + y_{t-37}) + \Phi_3 (y_{t-36}- y_{t-37} - y_{t-48} + y_{t-49}) \]

\[ + \phi_1 \Phi_1 (y_{t-13} - y_{t-14} -y_{t-25} -y_{t-26}) + \phi_1 \Phi_2 (y_{t-25} - y_{t-26} -y_{t-37} -y_{t-38})+ \phi_1 \Phi_3 (y_{t-37} - y_{t-38} -y_{t-49} -y_{t-50})\\ + \phi_2 \Phi_1 (y_{t-14} - y_{t-15} -y_{t-26} -y_{t-27})+ \phi_2 \Phi_2 (y_{t-26} - y_{t-27} -y_{t-38} -y_{t-39}) + \phi_2 \Phi_3 (y_{t-38} - y_{t-39} -y_{t-50} -y_{t-51})\\ + \phi_3 \Phi_1 (y_{t-15} - y_{t-16} -y_{t-27} -y_{t-28}) + \phi_3 \Phi_2 (y_{t-27} - y_{t-28} -y_{t-39} -y_{t-40}) + \phi_3 \Phi_3 (y_{t-39} - y_{t-40} -y_{t-51} -y_{t-52})\\ + \theta_1 \Theta_1 a_{t-13} + \theta_1 \Theta_2 a_{t-25} + \theta_1 \Theta_3 a_{t-37} + \theta_2 \Theta_1 a_{t-14} + \theta_2 \Theta_2 a_{t-26} + \theta_2 \Theta_3 a_{t-38} + \theta_3 \Theta_1 a_{t-15} + \theta_3 \Theta_2 a_{t-27} + \theta_3 \Theta_3 a_{t-39}\\ \\ \]

set.seed(123)
at <- rnorm(1052, mean = 0, sd = 2)

phi  <- c(0.4, -0.2, 0.1)
Phi  <- c(0.3, -0.15, 0.1)
theta <- c(0.5, -0.3, 0.2)
Theta <- c(0.4, -0.2, 0.1)

n <- length(at)
y <- numeric(n)

###### Inisialisasi nilai awal 
y[1:52] <- 0

for (t in 53:n) {
  # Differencing ganda
  diff_y <- y[t-1] + y[t-12] - y[t-13]
  
  # Komponen AR
  ar_non_seasonal <- 
    phi[1] * (y[t-1] - y[t-2] - y[t-13] + y[t-14]) +
    phi[2] * (y[t-2] - y[t-3] - y[t-14] + y[t-15]) +
    phi[3] * (y[t-3] - y[t-4] - y[t-15] + y[t-16])
  
  ar_seasonal <- 
    Phi[1] * (y[t-12] - y[t-13] - y[t-24] + y[t-25]) +
    Phi[2] * (y[t-24] - y[t-25] - y[t-36] + y[t-37]) +
    Phi[3] * (y[t-36] - y[t-37] - y[t-48] + y[t-49])
  
  # Interaksi Cross AR (φΦ)
  cross_ar <- 
    phi[1]*Phi[1] * (y[t-13] - y[t-14] - y[t-25] + y[t-26]) +
    phi[1]*Phi[2] * (y[t-25] - y[t-26] - y[t-37] + y[t-38]) +
    phi[1]*Phi[3] * (y[t-37] - y[t-38] - y[t-49] + y[t-50]) +
    phi[2]*Phi[1] * (y[t-14] - y[t-15] - y[t-26] + y[t-27]) +
    phi[2]*Phi[2] * (y[t-26] - y[t-27] - y[t-38] + y[t-39]) +
    phi[2]*Phi[3] * (y[t-38] - y[t-39] - y[t-50] + y[t-51]) +
    phi[3]*Phi[1] * (y[t-15] - y[t-16] - y[t-27] + y[t-28]) +
    phi[3]*Phi[2] * (y[t-27] - y[t-28] - y[t-39] + y[t-40]) +
    phi[3]*Phi[3] * (y[t-39] - y[t-40] - y[t-51] + y[t-52])
  
  # Komponen MA
  ma_non_seasonal <- 
    theta[1] * at[t-1] + 
    theta[2] * at[t-2] + 
    theta[3] * at[t-3]
  
  ma_seasonal <- 
    Theta[1] * at[t-12] + 
    Theta[2] * at[t-24] + 
    Theta[3] * at[t-36]
  
  # Interaksi cross terms MA (θΘ)
  cross_ma <- 
    theta[1] * Theta[1] * at[t-13] +
    theta[1] * Theta[2] * at[t-25] +
    theta[1] * Theta[3] * at[t-37] +
    theta[2] * Theta[1] * at[t-14] +
    theta[2] * Theta[2] * at[t-26] +
    theta[2] * Theta[3] * at[t-38] +
    theta[3] * Theta[1] * at[t-15] +
    theta[3] * Theta[2] * at[t-27] +
    theta[3] * Theta[3] * at[t-39]
  
  y[t] <- diff_y + ar_non_seasonal + ar_seasonal + cross_ar + at[t] + ma_non_seasonal + ma_seasonal + cross_ma }

y_ts <- y[101:length(y)]
### Plot
ts.plot(y_ts, main = "Original Time  Series", ylab = "y")

### ACF dan PACF plot Original Time  Series
par(mfrow = c(1,2))
acf(y_ts, main = "ACF of y_ts")
pacf(y_ts, main = "PACF of y_ts")

diff1 <- diff(y_ts, 1)
ts.plot(diff1, main = "Time Series - Non-Seasonal Differencing",ylab= "y")

### ACF dan PACF setelah differencing non-seasonal
par(mfrow = c(1,2))
acf(diff1, main = "ACF of diff(y_ts, 1)")
pacf(diff1, main = "PACF of diff(y_ts, 1)")

diff2 <- diff(y_ts, lag =12)
ts.plot(diff2, main = "Time Series - Seasonal Differencing",ylab= "y")

### ACF dan PACF setelah differencing seasonal
par(mfrow = c(1,2))
acf(diff2, main = "ACF of after diff(y_ts, D=1, s=12)")
pacf(diff2, main = "PACF of after diff(y_ts, D=1, s=12)")

# Plot time series setelah differencing non seasonal dan seasonal
diff3 <- diff(diff(y_ts,1), lag =12,1)
ts.plot(diff3, main = "Time Series - After Combined Differencing",ylab= "y")

### ACF dan PACF setelah differencing non seasonal dan seasonal
par(mfrow = c(1,2))
acf(diff3, main = "ACF of after diff d=1 D=1, s=12")
pacf(diff3, main = "PACF of after diff d=1 D=1, s=12")