Simulasi Estimasi Parameter SARIMA pada Data Konsumsi Energi Data Center AI

Import Library dan Menentukan Set Seed

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

set.seed(027)

Pada tahap awal ini dilakukan persiapan dengan memanggil beberapa package yang dibutuhkan. Package forecast digunakan untuk melakukan pemodelan ARIMA/SARIMA, sedangkan ggplot2 digunakan untuk membuat visualisasi data, dan dplyr untuk manipulasi data. Selain itu, fungsi set.seed(27) digunakan agar proses simulasi menghasilkan output yang konsisten setiap kali kode dijalankan, sehingga hasil penelitian dapat direplikasi.

Penentuan Parameter Asli

true_phi   <- 0.6
true_theta <- -0.5
true_sigma <- 1

Bagian ini menetapkan parameter asli (true parameter) dari model yang akan disimulasikan. Parameter \(\phi\) (phi) sebesar 0.6 merepresentasikan komponen autoregressive orde satu (AR(1)), sedangkan \(\theta\) (theta) sebesar -0.5 merupakan komponen moving average musiman orde satu (SMA(1)). Nilai standar deviasi error ditetapkan sebesar 1. Parameter-parameter ini nantinya akan menjadi acuan dalam mengevaluasi hasil estimasi model, khususnya dalam menghitung bias.

Fungsi Pembangkit Data Simulasi

generate_data <- function(n, add_outlier = FALSE) {
  burn_in <- 100
  total <- n + burn_in
  
  e <- rnorm(total, 0, true_sigma)
  w <- numeric(total)
  
  for (t in 8:total) {
    w[t] <- true_phi * w[t-1] + e[t] + true_theta * e[t-7]
  }
  
  y <- numeric(total)
  y[1:7] <- 50 + rnorm(7, 0, 2)
  
  for (t in 8:total) {
    y[t] <- y[t-7] + w[t]
  }
  
  y <- y[(burn_in+1):total]
  
  # Tambah outlier
  if (add_outlier) {
    idx <- sample(1:n, size = floor(0.05*n))
    y[idx] <- y[idx] + 5*true_sigma
  }
  
  return(y)
}

Fungsi ini digunakan untuk membangkitkan data time series yang mengikuti struktur model SARIMA. Proses dimulai dengan membangkitkan error acak (white noise), kemudian digunakan untuk membentuk komponen proses autoregressive dan seasonal moving average pada variabel \(w_t\)​. Selanjutnya, nilai \(w_t\)​ diintegrasikan secara musiman untuk menghasilkan data akhir \(y_t\)​, yang merepresentasikan data konsumsi energi. Bagian burn-in digunakan untuk menghilangkan pengaruh nilai awal agar data lebih stabil. Selain itu, fungsi ini juga menyediakan opsi untuk menambahkan outlier berupa lonjakan ekstrem pada sebagian kecil data, yang bertujuan untuk mensimulasikan kondisi data yang tidak bersih.

Output dari fungsi ini adalah sebuah vektor data time series sepanjang \(n\).

Menampilkan Data Hasil Simulasi

data_contoh <- generate_data(30, add_outlier = TRUE)

tabel_data <- data.frame(
  Waktu = 1:30,
  Nilai = data_contoh
)

head(tabel_data, 10)

##    Waktu    Nilai
## 1      1 51.79698
## 2      2 51.89604
## 3      3 54.15316
## 4      4 49.12209
## 5      5 49.18368
## 6      6 53.76931
## 7      7 56.69035
## 8      8 52.12547
## 9      9 51.78915
## 10    10 53.87331

Pada bagian ini ditampilkan contoh data hasil simulasi dalam bentuk tabel. Data terdiri dari dua kolom, yaitu waktu dan nilai observasi. Tampilan ini bertujuan untuk memberikan gambaran awal mengenai struktur data yang dihasilkan, termasuk kemungkinan adanya lonjakan nilai akibat outlier.

Output yang dihasilkan berupa tabel beberapa observasi pertama dari data simulasi.

Visualisasi Data Time Series

ggplot(tabel_data, aes(x=Waktu, y=Nilai)) +
  geom_line(color="blue") +
  geom_point() +
  labs(title="Contoh Data Simulasi SARIMA",
       x="Waktu",
       y="Nilai") +
  theme_minimal()

Visualisasi ini digunakan untuk melihat pola data terhadap waktu. Melalui grafik ini dapat diamati adanya pola musiman serta fluktuasi nilai dari waktu ke waktu. Jika terdapat outlier, maka akan terlihat sebagai lonjakan tajam pada titik tertentu dalam grafik.

Output yang dihasilkan adalah grafik time series yang menggambarkan perubahan nilai observasi sepanjang waktu. Grafik tersebut menunjukkan pola data time series hasil simulasi SARIMA selama 30 periode. Terlihat adanya fluktuasi nilai yang relatif stabil dengan pola yang berulang, yang mengindikasikan adanya komponen musiman (periode 7). Selain itu, beberapa titik yang melonjak lebih tinggi dapat menunjukkan adanya outlier atau gangguan pada data.

ggplot(tabel_data, aes(x=Nilai)) +
  geom_histogram(aes(y=..density..), bins=10, fill="skyblue", alpha=0.7) +
  geom_density(color="red", linewidth=1) +
  labs(title="Distribusi Data Simulasi") +
  theme_minimal()

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Bagian ini bertujuan untuk melihat distribusi data yang dihasilkan dari simulasi. Histogram digunakan untuk menggambarkan frekuensi data, sedangkan kurva density memberikan gambaran bentuk distribusi secara halus. Keberadaan outlier biasanya akan menyebabkan distribusi menjadi tidak simetris atau memiliki ekor yang lebih panjang.

Output berupa grafik histogram yang dilengkapi dengan kurva density. Grafik tersebut menunjukkan distribusi data hasil simulasi dalam bentuk histogram yang dilengkapi kurva density. Terlihat bahwa data menyebar di sekitar nilai tengah tertentu dengan bentuk distribusi yang cukup mendekati normal, meskipun terdapat sedikit kemencengan yang kemungkinan disebabkan oleh adanya outlier.

Analisis ACF dan PACF

ts_data <- ts(data_contoh, frequency=7)

acf(ts_data, main="ACF")

pacf(ts_data, main="PACF")

Analisis ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function) dilakukan untuk mengidentifikasi pola ketergantungan dalam data time series. ACF menunjukkan hubungan antar observasi pada berbagai lag, sedangkan PACF membantu dalam menentukan orde dari komponen autoregressive. Pola yang muncul pada grafik ini dapat digunakan sebagai dasar dalam pemilihan model SARIMA yang sesuai.

Berdasarkan grafik ACF dan PACF, terlihat adanya pola musiman yang ditunjukkan oleh lonjakan signifikan pada lag ke-7. Pola pada ACF menunjukkan cut off pada lag musiman tersebut, sementara PACF cenderung menurun secara bertahap (tailing off). Hal ini mengindikasikan adanya komponen Seasonal Moving Average orde satu (SMA(1)).

Pada bagian non-musiman, terdapat lonjakan signifikan pada lag awal, yang menunjukkan kemungkinan adanya komponen Autoregressive orde satu (AR(1)). Dengan demikian, model yang sesuai adalah SARIMA dengan kombinasi AR(1) pada bagian non-musiman dan SMA(1) pada bagian musiman.

Simulasi Estimasi Parameter

n_rep <- 100
ukuran <- c(60, 180, 365)
outlier_flag <- c(FALSE, TRUE)

hasil <- data.frame()

for (n in ukuran) {
  for (out in outlier_flag) {
    
    kategori <- case_when(
      n == 60 ~ "Kecil",
      n == 180 ~ "Sedang",
      TRUE ~ "Besar"
    )
    
    for (i in 1:n_rep) {
      
      data <- generate_data(n, out)
      ts_data <- ts(data, frequency=7)
      
      fit <- tryCatch({
        Arima(ts_data,
              order=c(1,0,0),
              seasonal=list(order=c(0,1,1), period=7),
              include.mean=FALSE)
      }, error=function(e) NULL)
      
      phi_hat <- NA
      theta_hat <- NA
      
      if (!is.null(fit)) {
        phi_hat <- fit$coef["ar1"]
        theta_hat <- fit$coef["sma1"]
      }
      
      hasil <- rbind(hasil, data.frame(
        Kategori=kategori,
        N=n,
        Outlier=out,
        Phi=phi_hat,
        Theta=theta_hat
      ))
    }
  }
}

Bagian ini mengatur skenario simulasi yang akan dilakukan. Simulasi dilakukan sebanyak 100 kali untuk setiap kombinasi ukuran sampel dan kondisi data (dengan atau tanpa outlier). Ukuran sampel dibagi menjadi tiga kategori, yaitu kecil, sedang, dan besar.

Selanjutnya dilakukan proses estimasi parameter model SARIMA untuk setiap data hasil simulasi. Proses ini menghasilkan kumpulan data estimasi parameter dari berbagai kondisi simulasi.

Tabel Ringkasan Bias

ringkasan <- hasil %>%
  na.omit() %>%
  group_by(Kategori, N, Outlier) %>%
  summarise(
    Mean_Phi = mean(Phi),
    Mean_Theta = mean(Theta),
    Bias_Phi = Mean_Phi - true_phi,
    Bias_Theta = Mean_Theta - true_theta,
    .groups="drop"
  )

ringkasan

## # A tibble: 6 × 7
##   Kategori     N Outlier Mean_Phi Mean_Theta Bias_Phi Bias_Theta
##   <chr>    <dbl> <lgl>      <dbl>      <dbl>    <dbl>      <dbl>
## 1 Besar      365 FALSE      0.593     -0.503 -0.00671   -0.00272
## 2 Besar      365 TRUE       0.321     -0.679 -0.279     -0.179  
## 3 Kecil       60 FALSE      0.576     -0.541 -0.0242    -0.0410 
## 4 Kecil       60 TRUE       0.281     -0.725 -0.319     -0.225  
## 5 Sedang     180 FALSE      0.597     -0.516 -0.00257   -0.0156 
## 6 Sedang     180 TRUE       0.321     -0.679 -0.279     -0.179

Pada bagian ini dilakukan perhitungan rata-rata estimasi parameter serta biasnya. Bias dihitung sebagai selisih antara nilai estimasi rata-rata dengan parameter asli. Tabel ringkasan ini digunakan untuk mengevaluasi seberapa baik model dalam mengestimasi parameter pada berbagai kondisi.

Output berupa tabel yang berisi nilai rata-rata estimasi dan bias untuk setiap skenario.

Visualisasi Distribusi Estimasi

# Mengubah TRUE/FALSE menjadi teks yang informatif
hasil$Outlier <- ifelse(hasil$Outlier == TRUE, "Ada Training Bursts", "Normal (Tanpa Bursts)")
ringkasan$Outlier <- ifelse(ringkasan$Outlier == TRUE, "Ada Training Bursts", "Normal (Tanpa Bursts)")

ggplot(hasil, aes(x=Kategori, y=Phi, fill=Outlier)) +
  geom_boxplot() +
  geom_hline(yintercept=true_phi, linetype="dashed", color="red") +
  labs(title="Estimasi Parameter Phi",
       y="Nilai Estimasi") +
  theme_bw()

Grafik ini membandingkan hasil estimasi parameter \(\phi\) pada tiga kategori (Besar, Kecil, dan Sedang) di bawah dua kondisi yang berbeda: kondisi “Normal (Tanpa Bursts)” yang diwakili oleh kotak berwarna toska, dan kondisi “Ada Training Bursts” yang diwakili oleh kotak berwarna merah muda. Garis putus-putus merah pada angka 0.6 menunjukkan nilai target atau nilai parameter yang sebenarnya (true value).

Temuan yang paling mencolok dari grafik ini adalah adanya perbedaan akurasi yang sangat drastis antara kedua kondisi tersebut. Pada kondisi “Normal (Tanpa Bursts)”, metode estimasi bekerja dengan sangat baik dan akurat. Hal ini terlihat dari garis tengah (median) pada ketiga kotak toska yang letaknya sangat berimpit dengan garis target 0.6. Sebaliknya, ketika terjadi kondisi “Ada Training Bursts”, performa estimasi menurun drastis dan mengalami bias ke bawah (underestimasi) yang parah. Median dari hasil estimasi pada kondisi ini anjlok ke kisaran 0.25 hingga 0.35, jauh di bawah nilai sebenarnya, apa pun kategori yang digunakan.

Selain masalah akurasi, grafik ini juga menunjukkan informasi mengenai tingkat presisi atau variasi data dari masing-masing “Kategori”. Baik pada kondisi normal maupun saat ada bursts, kategori “Kecil” selalu memiliki rentang kotak dan garis vertikal (whiskers) yang paling panjang, menandakan bahwa varians atau ketidakpastian estimasinya paling tinggi. Seiring dengan perubahan kategori menjadi “Sedang” dan “Besar”, sebaran data cenderung menyempit yang berarti hasil estimasinya menjadi lebih konsisten. Terdapat juga beberapa titik hitam di luar garis batas yang menunjukkan adanya nilai pencilan (outliers).

Sebagai kesimpulan, kehadiran Training Bursts pada data memberikan dampak negatif yang sangat signifikan terhadap proses estimasi parameter \(\phi\), yaitu menyebabkan model gagal menebak nilai asli dan menghasilkan bias ke bawah yang besar. Memperbesar kategori (dari Kecil ke Besar) memang berhasil meningkatkan presisi (mengurangi sebaran data), namun tidak mampu mengoreksi masalah bias yang diakibatkan oleh bursts tersebut.

Visualisasi Bias

ggplot(ringkasan, aes(x=Kategori, y=Bias_Phi, fill=Outlier)) +
  geom_bar(stat="identity", position="dodge") +
  labs(title="Perbandingan Bias Phi",
       y="Bias") +
  theme_minimal()

Grafik batang berjudul “Perbandingan Bias Phi” ini berfungsi untuk mempertegas temuan dari grafik boxplot sebelumnya dengan memvisualisasikan secara langsung nilai bias (selisih antara hasil rata-rata estimasi dengan nilai target yang sebenarnya) pada ketiga kategori. Garis nol (0.0) pada sumbu Y merepresentasikan estimasi yang sempurna tanpa bias. Karena seluruh batang grafik mengarah ke bawah (bernilai negatif), ini secara jelas menunjukkan bahwa seluruh hasil estimasi mengalami underestimasi, yaitu menebak nilai yang lebih rendah dari aslinya.

Pada kondisi “Normal (Tanpa Bursts)” yang ditunjukkan oleh batang berwarna toska, terlihat bahwa nilai biasnya sangat kecil dan nyaris menyentuh angka nol. Meskipun pada kategori “Kecil” terdapat sedikit pelebaran bias negatif dibandingkan “Besar” dan “Sedang”, secara keseluruhan performa estimasi dalam kondisi normal ini terbukti sangat akurat. Hal ini menunjukkan bahwa jika data bersih dari gangguan, metode estimasi tersebut bekerja dengan sangat optimal.

Sebaliknya, kondisi berubah drastis ketika “Ada Training Bursts” (batang berwarna merah muda). Terjadi lonjakan bias negatif yang sangat tajam yang mencapai kisaran -0.28 hingga melebihi -0.30. Dari ketiga kategori tersebut, kategori “Kecil” menanggung dampak bias (penyimpangan) yang paling parah. Meskipun kategori “Besar” dan “Sedang” memiliki tingkat bias yang sedikit lebih baik daripada “Kecil”, nilainya tetap saja sangat jauh dari angka nol.

Secara keseluruhan, grafik ini memberikan kesimpulan yang sangat tegas: gangguan berupa training bursts merusak keakuratan estimasi secara sistematis. Memperbesar ukuran kategori (dari Kecil ke Besar) terbukti tidak efektif untuk mengoreksi atau mengembalikan nilai estimasi ke target yang sebenarnya selama bursts tersebut masih ada di dalam data.