Simulasi Arima (3,0,3)(3,0,3)^12

model: ARIMA (3,0,3)(3,0,3)^12

\[ (1-\phi_1B - \phi_2 B^2 -\phi_3 B^3 )(1-\Phi_1 B^{12}- \Phi_2 B^{24} - \Phi_3 B^{36}) 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 \]

atau dalam model yt:

\[ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \Phi_1 y_{t-12} + \Phi_2 y_{t-24} + \Phi_3 y_{t-36} - \phi_1 \Phi_1 y_{t-13} - \phi_1 \Phi_2 y_{t-25} - \phi_1 \Phi_3 y_{t-37} \\ - \phi_2 \Phi_1 y_{t-14} - \phi_2 \Phi_2 y_{t-26} - \phi_2 \Phi_3 y_{t-38} - \phi_3 \Phi_1 y_{t-15} - \phi_3 \Phi_2 y_{t-27} - \phi_3 \Phi_3 y_{t-39} \\ +a_t -\theta_1a_{t-1} - \theta_2 a_{t-2} - \theta_3 a_{t-3} - \Theta_1a_{t-12} - \Theta_2 a_{t-24} - \Theta_3 a_{t-36} + \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(999)
a <- rnorm(2000, mean = 0, sd = 0.5)

#AR 
phi <- c(0.65, 0.35, -0.35)
Phi <- c(0.8, 0.4, -0.3) 

#MR 
theta <- c(0.4, -0.3, 0.35)
Theta <- c(0.5,-0.55,-0.4)

#seasonal 
s <- 12

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

#inisialisasi nilai y awal dengan random normal 
y[1:39] = rnorm(39, mean = 7, sd=4)

#Membuat model yt dengan pola ARIMA (3,0,3)(3,0,3)^12 

for(t in 40:n){
  ar = 0
  ma = a[t]
  for(i in 1:3){
    ar = ar + phi[i]*y[t-i] + Phi[i]*y[t-i*s]
    ma = ma - theta[i]*a[t-i] - Theta[i]*a[t-i*s]
    for(j in 1:3){
      ar = ar - phi[i]*Phi[j]*y[t-(s*j+i)]
      ma = ma + theta[i]*Theta[j]*a[t-(s*j+i)]
    }
  }
  y[t] = ar+ma
}

#menghapus data buffer
y = y[-(1:250)]

head(y, 100)

  [1] -1.83730790  3.28900239 -2.11157285  3.82624956 -0.75731167  1.22921301  0.86745084 -1.36266582
  [9] -1.39815992 -3.06576101 -3.54631596 -1.10842292 -2.28459175  2.42466061 -2.22454162  3.48092041
 [17]  0.17345192  2.73390891  2.30373891 -0.69267886 -0.35414076 -2.77836759 -3.18589739 -3.19630868
 [25] -2.61974028  1.36576536 -2.40657028  3.69788162  0.95227433  3.23741326  2.82736459 -0.40975926
 [33] -1.24910876 -3.20954834 -4.57234414 -2.90742484 -3.11894154 -0.04651996 -3.45206171  0.94773863
 [41] -1.04799685  1.30629775  2.54021101 -1.04199212  0.94239156 -2.55603824 -3.01653250 -3.88977777
 [49] -2.38671799  0.19028946 -1.48706747  1.97618047  1.70960250  2.25405629  3.57102213 -1.12622253
 [57]  0.12832559 -2.89752460 -3.34438345 -3.04763895 -2.63995668 -1.35712164 -3.41729507  0.29308779
 [65] -0.18777257  2.05951277  2.43648786 -1.38865835  0.45412401 -2.46839697 -3.20221810 -2.82913894
 [73] -3.15716154 -0.85056381 -1.92960846  0.58079832  1.42142767  1.98882257  3.29817837 -1.34290008
 [81]  0.67942708 -3.03485011 -3.03442352 -2.01163924 -3.36140553 -0.80435735 -2.41432449  0.52166606
 [89]  0.19903634  2.25770162  2.42203756 -1.09009992  1.12987525 -0.97620869 -1.54636397 -1.39235509
 [97] -3.35495384 -2.54339089 -2.95001343 -1.31237487

#plot data time series secara keseluruhan
ts.plot(y)

#plotting subset y untuk memperjelas pola data
plot(1:100, y[1:100], type = "l", main = "y[1:100]", xlab = "Time", ylab = "y")+
plot(300:400, y[300:400], type = "l", main = "y[300:400]", xlab = "Time", ylab = "y")+
plot(750:800, y[750:800], type = "l", main = "y[750:800]", xlab = "Time", ylab = "y")

integer(0)

acf(y, lag.max = 50)

pacf(y, lag.max = 50)