ARIMA Seasonal Simulation

NAMA: Tito Fauzan Putra
NRP: 5003221164
Kelas: Pengantar Analisis Deret Waktu A

ARIMA(0,0,0)(3,1,0) Model

Model Dengan Backshift Operator:

\[ (\Phi_{3})B^{12}(1-B^{12})y_t=a_t \\ (1-\Phi_1B^{12}-\Phi_2B^{24}-\Phi_3B^{36})(1-B^{12})y_t=a_t \\ \]

Jika backshift di konversi menjadi yt maka:

\[ y_{t-n} = B^ny_t \]

Maka persamaan eksplisit menjadi:

\[ (1 - \Phi_1 B^{12} - \Phi_2 B^{24} -\Phi_3 B^{36}) (y_t - B^{12} y_t) = a_t \\ y_t - B^{12} y_t - \Phi_1B^{12} y_t + \Phi_1 B^{24} y_t - \Phi_2 B^{24} y_t + \Phi_2 B^{36} y_t - \Phi_3 B^{36} y_t + \Phi_3 B^{48} y_t = a_t \\ y_t = B^{12} y_t + \Phi_1 B^{12} y_t - \Phi_1B^{24} y_t + \Phi_2B^{24}y_t-\Phi_2B^{36}y_t+\Phi_3B^{36}y_t-\Phi_3B^{48}y_t + a_t \\ y_t = y_{t-12} + \Phi_1 (y_{t-12} - y_{t-24}) + \Phi_2 (y_{t-24} - y_{t-36}) + \Phi_3 (y_{t-36} - y_{t-48}) + a_t \]

Membuat Data SARIMA Dengan Rumus Manual

Menentukan parameter model SARIMA

Model: SARIMA(0,0,0)(3,1,0)^12 SARIMA(0,0,0)(3,1,0)^12

• p=0, d=0, q=0 (orde non-musiman)
• P=3, D=1, Q=0 (orde musiman)
• S=12 (periode musiman)

Menentukan koefisien AR musiman (Phi):

• SAR(1) = 0,7
• SAR(2) = 0,6
• SAR(3) = -0,8

Menentukan Parameter Model

# Parameter Model
sphi1 <- 0.7
sphi2 <- 0.6
sphi3 <- -0.8
sd_a <- 1 # Standar deviasi dari a_t (noise)
phi_seasonal <- c(sphi1, sphi2, sphi3)

# Panjang data
n_desired <- 1000
n_burn_in_manual <- 100 # Periode burn-in
n_total <- n_desired + n_burn_in_manual

Melakukan Simulasi Model Manual

# Inisialisasi vektor y untuk menyimpan semua nilai (termasuk burn-in)
y_manual <- numeric(n_total)

# Menghasilkan data menggunakan loop
set.seed(164) # Untuk reproduktifitas
for (t in 1:n_total) {
  # Hasilkan white noise a_t
  a_t <- rnorm(1, mean = 0, sd = sd_a)
  
  # Ambil nilai lag, jika t-k <= 0, anggap y_{t-k} = 0 (umum untuk inisialisasi)
  y_t_minus_12 <- ifelse(t - 12 > 0, y_manual[t - 12], 0)
  y_t_minus_24 <- ifelse(t - 24 > 0, y_manual[t - 24], 0)
  y_t_minus_36 <- ifelse(t - 36 > 0, y_manual[t - 36], 0)
  y_t_minus_48 <- ifelse(t - 48 > 0, y_manual[t - 48], 0)
  
  # Terapkan formula
  if (t <= 12) { # Tidak ada lag musiman yang cukup
     y_manual[t] <- a_t 
  } else if (t <= 24) { # Hanya bisa menghitung sampai lag 12
     y_manual[t] <- y_t_minus_12 + 
                   a_t 
  } else if (t <= 36) { # Hanya bisa menghitung sampai lag 24
     y_manual[t] <- y_t_minus_12 + 
                   sphi1 * (y_t_minus_12 - y_t_minus_24) + 
                   a_t 
  } else if (t <= 48) { # Hanya bisa menghitung sampai lag 36
    y_manual[t] <- y_t_minus_12 + 
                   sphi1 * (y_t_minus_12 - y_t_minus_24) + 
                   sphi2 * (y_t_minus_24 - y_t_minus_36) + 
                   a_t 
  } else { # t > 48, semua lag tersedia
    y_manual[t] <- y_t_minus_12 + 
                   sphi1 * (y_t_minus_12 - y_t_minus_24) + 
                   sphi2 * (y_t_minus_24 - y_t_minus_36) + 
                   sphi3 * (y_t_minus_36 - y_t_minus_48) + 
                   a_t
  }
}

# Buang periode burn-in untuk mendapatkan Y_t final
data_sarima_manual_Yt <- y_manual[(n_burn_in_manual + 1):n_total]

# Melihat beberapa data pertama dari Y_t
print(head(data_sarima_manual_Yt))

[1]   3.3307716   2.8772205   4.6720470  -3.9119302   0.4191986 -12.6695242

# Memvisualisasikan data Y_t (opsional)
plot(ts(data_sarima_manual_Yt), main = "Data Simulasi Yt SARIMA(0,0,0)(3,1,0)^12 Manual Loop", 
     ylab = "Nilai Yt", col = "red")

Plot Yt ACF & PACF Data Simulai ARIMA(0, 0, 0)(3, 1, 0)^12 Tanpa Library

# Menghitung dan memplot ACF Yt
acf(data_sarima_manual, 
    lag.max = 48, # Tampilkan hingga lag 48 untuk melihat beberapa siklus musiman
    main = "Yt ACF untuk SARIMA(0,0,0)(3,1,0)^12 Tanpa Library",
    ylab = "ACF")

# Menghitung dan memplot PACF Yt
pacf(data_sarima_manual, 
     lag.max = 48, # Tampilkan hingga lag 48
     main = "Yt PACF untuk SARIMA(0,0,0)(3,1,0)^12 Tanpa Library",
     ylab = "PACF")

Menghitung Wt ARIMA(0,0,0)(3,1,0)^12 Tanpa Library

# Hitung Wt = Yt - Y_{t-12}
# data_sarima_manual_Yt adalah seri Yt setelah burn-in
# Wt akan memiliki panjang n_desired - 12
data_sarima_manual_Wt <- diff(data_sarima_manual_Yt, lag = 12)

# Melihat beberapa data pertama dari Wt
cat("\nBeberapa data pertama dari Wt:\n")

Beberapa data pertama dari Wt:

print(head(data_sarima_manual_Wt))

[1] -2.331840  1.751699  2.436126  1.353494 -5.649257 -3.176931

# Memvisualisasikan data Wt
plot(ts(data_sarima_manual_Wt), 
     main = "Data Diferensiasi Musiman Wt = Yt - Y(t-12) Tanpa Library", 
     ylab = "Nilai Wt", col = "purple")

Plot Wt ACF & PACF Data Simulai ARIMA(0, 0, 0)(3, 1, 0)^12 Tanpa Library

# Menghitung dan memplot ACF Wt
acf(data_sarima_manual_Wt, 
    lag.max = 48, # Tampilkan hingga lag 48 untuk melihat beberapa siklus musiman
    main = "Wt ACF untuk SARIMA(0,0,0)(3,1,0)^12 Tanpa Library",
    ylab = "ACF")

# Menghitung dan memplot PACF Wt
pacf(data_sarima_manual_Wt, 
     lag.max = 48, # Tampilkan hingga lag 48
     main = "Wt PACF untuk SARIMA(0,0,0)(3,1,0)^12 Tanpa Library",
     ylab = "PACF")

Kesimpulan:

• Plot ACF mengalami dies down pada lag kelipatan 12
• Plot PACF mengalami cutoff pada lag 36

Secara teori sudah benar

Membuat Data SARIMA(0,0,0)(3,1,0) Dengan Library

# Panjang data yang ingin dihasilkan
n_desired <- 1000

# Minimal 48 untuk model ini, tetapi lebih banyak lebih baik
n_burn_in_arima_sim <- n_burn_in_manual

# Membuat daftar model untuk arima.sim()
model_spec <- list(
  order = c(0, 0, 0), # (p,d,q) non-musiman
  seasonal = list(
    order = c(3, 1, 0), # (P,D,Q) musiman
    period = 12,
    sar = phi_seasonal # Koefisien AR musiman
  )
)

set.seed(1) # Untuk reproduktifitas
data_sarima_sim <- arima.sim(model = model_spec, 
                             n = n_desired, 
                             n.start = n_burn_in_arima_sim,
                             sd = 164) # Standar deviasi dari a_t

# Melihat beberapa data pertama
print(head(data_sarima_sim))

Time Series:
Start = 1 
End = 6 
Frequency = 1 
[1] -101.740135    6.907003 -149.391150   25.916719 -107.351882  289.835112

# Memvisualisasikan data (opsional)
plot(data_sarima_sim, main = "Data Simulasi Yt SARIMA(0,0,0)(3,1,0)^12 via arima.sim", 
     ylab = "Nilai", col = "blue")

Plot Yt ACF & PACF Data Simulai ARIMA(0, 0, 0)(3, 1, 0)^12 Via arima.sim

# Menghitung dan memplot ACF Yt
acf(data_sarima_sim, 
    lag.max = 48, # Tampilkan hingga lag 48 untuk melihat beberapa siklus musiman
    main = "Yt ACF untuk SARIMA(0,0,0)(3,1,0)^12 arima.sim",
    ylab = "ACF")

# Menghitung dan memplot PACF Yt
pacf(data_sarima_sim, 
     lag.max = 48, # Tampilkan hingga lag 48
     main = "Yt PACF untuk SARIMA(0,0,0)(3,1,0)^12 arima.sim",
     ylab = "PACF")

Menghitung Wt ARIMA(0,0,0)(3,1,0)^12 Dengan Library

# Hitung Wt = Yt - Y_{t-12}
# data_sarima_manual_Yt adalah seri Yt setelah burn-in
# Wt akan memiliki panjang n_desired - 12
data_sarima_library_Wt <- diff(data_sarima_sim, lag = 12)

# Melihat beberapa data pertama dari Wt
cat("\nBeberapa data pertama dari Wt:\n")

Beberapa data pertama dari Wt:

print(head(data_sarima_library_Wt))

Time Series:
Start = 13 
End = 18 
Frequency = 1 
[1]  336.63442 -113.62121  115.38071  -90.33722   54.87305 -335.60969

# Memvisualisasikan data Wt
plot(ts(data_sarima_library_Wt), 
     main = "Data Diferensiasi Musiman Wt = Yt - Y(t-12) Via arima.sim", 
     ylab = "Nilai Wt", col = "purple")

Plot Wt ACF & PACF Data Simulai ARIMA(0, 0, 0)(3, 1, 0)^12 Dengan Library

# Menghitung dan memplot ACF Wt
acf(data_sarima_library_Wt, 
    lag.max = 48, # Tampilkan hingga lag 48 untuk melihat beberapa siklus musiman
    main = "Wt ACF untuk SARIMA(0,0,0)(3,1,0)^12 Dengan Library",
    ylab = "ACF")

# Menghitung dan memplot PACF Wt
pacf(data_sarima_library_Wt, 
     lag.max = 60, # Tampilkan hingga lag 48
     main = "Wt PACF untuk SARIMA(0,0,0)(3,1,0)^12 Dengan Library",
     ylab = "PACF")

Kesimpulan:

• Plot ACF mengalami cut off pada lag 12
• Plot PACF mengalami dies down pada kelipatan 12

Secara teori belum benar

Model Tanpa Differencing

\[ \Phi_3 (B^{12}) y_t = a_t \\ y_t - \Phi_1 B^{12} y_t - \Phi_2 B^{24} y_t - \Phi_3 B^{36} y_t = a_t \\ y_t = \Phi_1 B^{12} y_t + \Phi_2 B^{24} y_t + \Phi_3 B^{36} + a_t \\ y_t = \Phi_1 y_{t-12} + \Phi_2 y_{t-24} + \Phi_3 y_{t-36} + a_t \]

Menentukan Parameter

# Parameter Model
phi1_m <- sphi1
phi2_m <- sphi2
phi3_m <- sphi3 # Koefisien yang sama dengan arima.sim di atas
sd_a_m <- sd_a   # Standar deviasi dari a_t (noise)

# Untuk model stasioner (D=0, d=0), kita bisa menambahkan konstanta c
constant_c <- 0

# Panjang data
n_desired_manual_P3D0 <- 1000
n_burn_in_manual_P3D0 <- 100 # Periode burn-in
n_total_manual_P3D0 <- n_desired_manual_P3D0 + n_burn_in_manual_P3D0

Simulasi Model ARIMA(0, 0, 0)(3, 0, 0)^12 Tanpa Library

# Inisialisasi vektor y untuk menyimpan semua nilai (termasuk burn-in)
y_manual_P3D0 <- numeric(n_total_manual_P3D0)

# Menghasilkan data menggunakan loop
set.seed(164) # Untuk reproduktifitas
for (t in 1:n_total_manual_P3D0) {
  # Hasilkan white noise a_t
  a_t <- rnorm(1, mean = 0, sd = sd_a_m)
  
  # Ambil nilai lag, jika t-k <= 0, anggap y_{t-k} = 0 (untuk mean nol proses awal)
  # Ini adalah pendekatan sederhana untuk inisialisasi; periode burn-in membantu.
  y_t_minus_12 <- ifelse(t - 12 > 0, y_manual_P3D0[t - 12], 0)
  y_t_minus_24 <- ifelse(t - 24 > 0, y_manual_P3D0[t - 24], 0)
  y_t_minus_36 <- ifelse(t - 36 > 0, y_manual_P3D0[t - 36], 0)
  
  # Terapkan formula
  # Kita perlu setidaknya 36 observasi sebelumnya untuk menggunakan formula lengkap.
  # Logika if/else ini menangani iterasi awal selama periode burn-in.
  if (t > 36) {
    y_manual_P3D0[t] <- constant_c + 
                        phi1_m * y_t_minus_12 + 
                        phi2_m * y_t_minus_24 + 
                        phi3_m * y_t_minus_36 + 
                        a_t
  } else if (t > 24) { # Hanya bisa menggunakan lag 12 dan 24
    y_manual_P3D0[t] <- constant_c + 
                        phi1_m * y_t_minus_12 + 
                        phi2_m * y_t_minus_24 + 
                        a_t # phi3_m * y_t_minus_36 tidak ada
  } else if (t > 12) { # Hanya bisa menggunakan lag 12
    y_manual_P3D0[t] <- constant_c + 
                        phi1_m * y_t_minus_12 + 
                        a_t # phi2_m dan phi3_m tidak ada
  } else { # Tidak ada lag musiman yang cukup
    y_manual_P3D0[t] <- constant_c + a_t
  }
}

# Buang periode burn-in
data_sarima_manual_P3D0 <- y_manual_P3D0[(n_burn_in_manual_P3D0 + 1):n_total_manual_P3D0]

# Melihat beberapa data pertama
print(head(data_sarima_manual_P3D0))

[1] -4.1002584  2.5593445  2.8826714  0.2394817 -6.4316180 -3.9651008

# Memvisualisasikan data (opsional)
plot(ts(data_sarima_manual_P3D0), 
     main = "Data Simulasi SARIMA(0,0,0)(3,0,0)^12 Manual Loop",
     ylab = "Nilai", col = "purple4")

Plot Yt ACF & PACF Data Simulai ARIMA(0, 0, 0)(3, 0, 0)^12 Tanpa Library

# Menghitung dan memplot ACF
acf(data_sarima_manual_P3D0, 
    lag.max = 48, # Tampilkan hingga lag 48 untuk melihat beberapa siklus musiman
    main = "Yt ACF untuk SARIMA(0,0,0)(3,0,0)^12 Manual",
    ylab = "ACF")

# Menghitung dan memplot PACF
pacf(data_sarima_manual_P3D0, 
     lag.max = 48, # Tampilkan hingga lag 48
     main = "Yt PACF untuk SARIMA(0,0,0)(3,0,0)^12 Manual",
     ylab = "PACF")

Kesimpulan:

• Plot ACF mengalami Dies down pada lag kelipatan 12
• Plot PACF mengalami cut off pada lag 36

Secara teori sudah benar

Menentukan Parameter

# Koefisien AR Musiman (Phi1, Phi2, Phi3)
phi_seasonal_P3D0 <- phi_seasonal

# Panjang data yang ingin dihasilkan (setelah burn-in)
n_desired_P3D0 <- 1000

# Periode burn-in (untuk membiarkan model stabil)
# Minimal 36 untuk model ini, tetapi lebih banyak lebih baik
n_burn_in_arima_sim_P3D0 <- 100

Simulasi Model ARIMA(0, 0, 0)(3, 0, 0)^12 Dengan Library

# Membuat daftar model untuk arima.sim()
model_spec_P3D0 <- list(
  order = c(0, 0, 0), # (p,d,q) non-musiman
  seasonal = list(
    order = c(3, 0, 0), # (P,D,Q) musiman -> D=0 sekarang
    period = 12,
    sar = phi_seasonal_P3D0 # Koefisien AR musiman
  )
)

# Menghasilkan data (secara default mean nol)
set.seed(456) # Untuk reproduktifitas
data_sarima_sim_P3D0 <- arima.sim(model = model_spec_P3D0,
                                  n = n_desired_P3D0,
                                  n.start = n_burn_in_arima_sim_P3D0,
                                  sd = 1) # Standar deviasi dari a_t

# Melihat beberapa data pertama
print(head(data_sarima_sim_P3D0))

Time Series:
Start = 1 
End = 6 
Frequency = 1 
[1]  0.11815133  0.86990262 -0.09193621  0.06889879 -1.68242675  1.11695555

# Memvisualisasikan data (opsional)
plot(data_sarima_sim_P3D0, 
     main = "Data Simulasi SARIMA(0,0,0)(3,0,0)^12 via arima.sim",
     ylab = "Nilai", col = "green4")

Plot ACF & PACF Data Simulai ARIMA(0, 0, 0)(3, 0, 0)^12 Via arima.sim

# Menghitung dan memplot ACF
acf(data_sarima_sim_P3D0, 
    lag.max = 48, # Tampilkan hingga lag 48 untuk melihat beberapa siklus musiman
    main = "ACF untuk SARIMA(0,0,0)(3,0,0)^12",
    ylab = "ACF")

# Menghitung dan memplot PACF
pacf(data_sarima_sim_P3D0, 
     lag.max = 48, # Tampilkan hingga lag 48
     main = "PACF untuk SARIMA(0,0,0)(3,0,0)^12",
     ylab = "PACF")

Kesimpulan:

• Plot ACF tidak ada lag yang signifikan, mungkin hanya pada lag 27
• Plot PACF yang memiliki lag yang signifikan ada pada lag 12, dan 27

Secara teori belum benar