1. Pendahuluan

Laporan ini menganalisis data konsentrasi karbon dioksida (CO2) atmosfer bulanan yang diukur di Observatorium Mauna Loa, Hawaii, dari tahun 1959 hingga 1997. Data ini tersedia sebagai dataset bawaan R (co2) dan merupakan salah satu deret waktu klasik yang menunjukkan pola trend naik jangka panjang beserta pola musiman (seasonal) yang stabil, sehingga cocok untuk dianalisis menggunakan model ARIMA (Autoregressive Integrated Moving Average).

Tujuan analisis:

  1. Mendeskripsikan karakteristik data deret waktu CO2.
  2. Menguji stasioneritas data dan melakukan differencing bila diperlukan.
  3. Mengidentifikasi dan membentuk model ARIMA yang sesuai.
  4. Melakukan diagnostik residual dan mengevaluasi akurasi forecast.
  5. Menginterpretasikan model akhir dan menarik kesimpulan.

2. Dasar Teori

2.1 Deret Waktu dan Stasioneritas

Deret waktu (time series) adalah sekumpulan observasi yang dikumpulkan secara berurutan berdasarkan waktu. Salah satu syarat penting sebelum memodelkan deret waktu dengan ARIMA adalah stasioneritas, yaitu kondisi di mana rata-rata, varians, dan autokovariansi data tidak berubah terhadap waktu. Data yang memiliki trend atau pola musiman biasanya bersifat non-stasioner dan perlu ditransformasi (misalnya melalui differencing) sebelum dimodelkan.

Uji stasioneritas yang umum digunakan adalah Augmented Dickey-Fuller (ADF) test, dengan hipotesis:

  • H0: data memiliki unit root (non-stasioner)
  • H1: data stasioner

2.2 Model ARIMA

Model ARIMA(p, d, q) menggabungkan tiga komponen:

  • AR(p): komponen autoregressive, menyatakan nilai saat ini dipengaruhi oleh p nilai sebelumnya.
  • I(d): differencing sebanyak d kali untuk mencapai stasioneritas.
  • MA(q): komponen moving average, menyatakan nilai saat ini dipengaruhi oleh q error masa lalu.

Karena data CO2 memiliki pola musiman tahunan (siklus 12 bulan), model yang lebih sesuai adalah SARIMA (Seasonal ARIMA), dengan notasi ARIMA(p,d,q)(P,D,Q)[s], di mana komponen (P,D,Q) menangkap pola musiman dengan periode s (di sini s = 12).

Identifikasi orde model umumnya dibantu dengan plot ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function), serta pemilihan otomatis berbasis AIC melalui fungsi auto.arima() dari paket forecast.

3. Deskripsi Data

3.1 Struktur dan Ringkasan Data

data(co2)
summary_co2 <- as.data.frame(t(as.matrix(summary(co2))))
knitr::kable(summary_co2, caption = "Ringkasan statistik deskriptif CO2 (ppm)", digits = 2)
Ringkasan statistik deskriptif CO2 (ppm)
Min. 1st Qu. Median Mean 3rd Qu. Max.
313.18 323.53 335.17 337.05 350.26 366.84

Berdasarkan statistik deskriptif, konsentrasi CO2 atmosfer selama periode januari 1959 hingga Desember 1997 memiliki nilai minimum sebesar 313.18 ppm dan maksimum sebesar 366.84 ppm. Nilai rata-rata konsentrasi CO2 adalah 337.05 ppm, sedangkan nilai mediannya sebesar 335.17 ppm. Selisih antara nilai minimum dan maksimum menunjukkan adanya variasi konsentrasi CO2 yang cukup besar selama periode pengamatan. Selain itu, nilai kuartil pertama sebesar 323.53 ppm menunjukkan bahwa sekitar 25% observasi memiliki konsentrasi CO2 di bawah atau sama dengan nilai tersebut, sedangkan kuartil ketiga sebesar 350.26 ppm menunjukkanbahwa sekitar 75% observasi berada di bawah atau sama dengan nilai tersebut. Secara keseluruhan, statistik deskriptif menunjukkan adanya peruahan konsentrasi CO2 selama periode pengamatan yang selanjutnya perlu dikaji melalui visualisasi deret waktu untuk melihat pola tren dan musiman.

Dataset co2 merupakan objek ts dengan frekuensi 12 (bulanan), dimulai dari Januari 1959 hingga Desember 1997, sehingga terdapat 468 observasi.

co2_df <- data.frame(
  Tahun = floor(time(co2))[1:12],
  Bulan = month.abb,
  CO2 = as.numeric(window(co2, start = c(1959,1), end = c(1959,12)))
)
knitr::kable(co2_df, caption = "Cuplikan data CO2 tahun 1959 (12 bulan pertama)")
Cuplikan data CO2 tahun 1959 (12 bulan pertama)
Tahun Bulan CO2
1959 Jan 315.42
1959 Feb 316.31
1959 Mar 316.50
1959 Apr 317.56
1959 May 318.13
1959 Jun 318.00
1959 Jul 316.39
1959 Aug 314.65
1959 Sep 313.68
1959 Oct 313.18
1959 Nov 314.66
1959 Dec 315.43

Pada Tabel menunjukkan 12 observasi konsentrasi CO2 pada tahun 1959. Konsentrasi CO2 meningkat dari 315.42 ppm pada Januari hingga mencapai 318.13 ppm pada Mei, kemudian mengalami penurunan hingga 313.18 ppm pada Oktober sebelum kembali meingkat pada November dan Desember. Pola kenaikan dan penurunan dalam satu tahun tersebut memberikan indikasi awal adanya pola musiman tahunan pada konsentrasi CO2, yang selanjutnya akan dianalisis melalui visualisasi runtun waktu dan dekomposisi.

3.2 Visualisasi Data

plot(co2, main = "Konsentrasi CO2 Atmosfer Bulanan (Mauna Loa, 1959-1997)",
     ylab = "CO2 (ppm)", xlab = "Tahun", col = "darkgreen")

Berdasarkan plot deret waktu, konsentrasi CO₂ menunjukkan tren peningkatan dalam jangka panjang selama periode 1959–1997. Konsentrasi CO₂ pada awal periode pengamatan berada pada kisaran sekitar 315 ppm, kemudian meningkat hingga mendekati 367 ppm pada akhir periode. Selain tren jangka panjang, terlihat adanya pola fluktuasi yang berulang setiap tahun, yang menunjukkan adanya komponen musiman dengan periode 12 bulan. Dengan adanya tren dan pola musiman tersebut, data pada level awal diduga belum stasioner. Oleh karena itu, diperlukan pengujian stasioneritas secara formal menggunakan Augmented Dickey-Fuller (ADF) test sebelum dilakukan pemodelan.

plot(decompose(co2))

Hasil dekomposisi menunjukkan bahwa deret waktu CO2 terdiri atas beberapa komponen utama, yaitu trend, seasonal, dan random. Komponen trend menunjukkan peningkatan konsentrasi CO2 secara bertahap selama periode pengamatan. Komponen seasonal menunjukkanadanya pola berulang periodek setiap tahun, sedangkan komponen random menunjukkan variasi yang tidak dijelaskan oleh komponen trend dan seasonal. Keberadaan pola musiman dengan periode 12 bulan mendukung penggunaan model SARIMA, karena model tersebut mengakomodasi ketergantungan pada pola musiman selain pola non-musiman.

4. Tahapan Analisis

4.1 Uji Stasioneritas (ADF Test)

Hipotesis

\[ H_0: \text{Data memiliki unit root (tidak stasioner)} \qquad H_1: \text{Data stasioner} \] Taraf signifikansi: \(\alpha = 5\%\)

adf_awal <- adf.test(co2)
tabel_adf_awal <- data.frame(
  Statistik = c("Dickey-Fuller", "Lag Order", "p-value"),
  Nilai = c(round(unname(adf_awal$statistic), 4),
            unname(adf_awal$parameter),
            round(adf_awal$p.value, 4))
)
knitr::kable(tabel_adf_awal, caption = "Hasil ADF Test pada Data CO2 (sebelum differencing)")
Hasil ADF Test pada Data CO2 (sebelum differencing)
Statistik Nilai
Dickey-Fuller -2.8299
Lag Order 7.0000
p-value 0.2269

Keputusan dan kesimpulan: Diperoleh statistik Dickey-Fuller sebesar -2.8299 dengan p-value sebesar 0.2269. Karena p-value > 0,05, maka \(H_0\) tidak ditolak. Dengan demikian, data CO2 pada level belum stasioner, sehingga diperlukan proses differencing sebelum pemodelan SARIMA.

4.2 Differencing

Differencing non-musiman orde satu dituliskan sebagai:

\[ \nabla Y_t = Y_t - Y_{t-1} \] Differencing dilakukan berurutan: non-musiman terlebih dahulu, kemudian musiman.

co2_diff <- diff(co2, differences = 1)
co2_diff_seasonal <- diff(co2_diff, lag = 12)
plot(co2_diff_seasonal, main = "Data CO2 Setelah Differencing (d=1, D=1, s=12)",
     ylab = "Selisih CO2", col = "steelblue")

4. Uji Stasioneritas Setelah Differencing

Hipotesis \[ H_0: \text{Data memiliki unit root (tidak stasioner)} \qquad H_1: \text{Data stasioner} \]

Taraf signifikansi: \(\alpha = 5\%\)

adf_diff <- adf.test(co2_diff_seasonal)
tabel_adf_diff <- data.frame(
  Statistik = c("Dickey-Fuller", "Lag Order", "p-value"),
  Nilai = c(round(unname(adf_diff$statistic), 4),
            unname(adf_diff$parameter),
            round(adf_diff$p.value, 4))
)
knitr::kable(tabel_adf_diff, caption = "Hasil ADF Test Setelah Differencing")
Hasil ADF Test Setelah Differencing
Statistik Nilai
Dickey-Fuller -8.9846
Lag Order 7.0000
p-value 0.0100

Keputusan dan kesimpulan: Setelah differencing non-musiman (\(d=1\)) dan musiman (\(D=1\), \(s=12\)), diperoleh statistik Dickey-Fuller sebesar -8.9846 dengan p-value sebesar 0.01. Karena p-value < 0,05, maka \(H_0\) ditolak, sehingga data hasil differencing dinyatakan stasioner dan dapat digunakan untuk tahap identifikasi dan pembentukan model SARIMA.

4.4 Identifikasi Model (ACF/PACF)

par(mfrow = c(1,2))
acf(co2_diff_seasonal, main = "ACF")
pacf(co2_diff_seasonal, main = "PACF")

par(mfrow = c(1,1))

ACF digunakan untuk melihat hubungan antara suatu observasi dengan observasi pada lag-lag sebelumnya, sedangkan PACF melihat hubungan parsial setelah pengaruh lag di antaranya diperhitungkan. Identifikasi orde \((p,q)\) dan \((P,Q)\) secara manual dari pola cut-off/tailing-off pada plot di atas dilengkapi dengan pemilihan otomatis berbasis AIC menggunakan auto.arima() pada bagian berikut. Interpretasi lag spesifik mana yang cut-off sebaiknya dibaca langsung dari bentuk plot hasil knit, bukan diasumsikan tanpa melihatnya.

4.4 Pembentukan Model dengan auto.arima()

model_auto <- auto.arima(co2, seasonal = TRUE, stepwise = FALSE, approximation = FALSE)
ord <- arimaorder(model_auto)
koef <- coef(model_auto)
se <- sqrt(diag(vcov(model_auto)))
tabel_koef <- data.frame(
  Parameter = names(koef),
  Estimate = round(unname(koef), 4),
  `Std. Error` = round(unname(se), 4),
  check.names = FALSE
)
knitr::kable(tabel_koef, caption = paste0("Koefisien Model ", arimaorder(model_auto) |> (\(o) paste0("SARIMA(", o[1], ",", o[2], ",", o[3], ")(", o[4], ",", o[5], ",", o[6], ")[", o[7], "]"))()))
Koefisien Model SARIMA(0,1,3)(0,1,1)[12]
Parameter Estimate Std. Error
ma1 -0.3394 0.0475
ma2 -0.0180 0.0497
ma3 -0.0973 0.0467
sma1 -0.8538 0.0256 
tabel_fit <- data.frame(
  Ukuran = c("AIC", "AICc", "BIC", "Log-Likelihood", "Sigma^2"),
  Nilai = c(round(model_auto$aic, 2), round(model_auto$aicc, 2),
            round(model_auto$bic, 2), round(model_auto$loglik, 2),
            round(model_auto$sigma2, 4))
)
knitr::kable(tabel_fit, caption = "Ukuran Kebaikan Model (Goodness of Fit)")
Ukuran Kebaikan Model (Goodness of Fit)
Ukuran Nilai
AIC 176.8600
AICc 177.0000
BIC 197.4700
Log-Likelihood -83.4300
Sigma^2 0.0852

Interpretasi: Model terpilih adalah SARIMA(\(0,1,3\))(\(0,1,1\))\(_{12}\). Komponen non-musiman terdiri atas MA(\(3\)) dan komponen musiman SMA(\(1\)), dengan differencing non-musiman sebanyak \(d=1\) kali dan musiman \(D=1\) kali pada periode 12 bulan. Koefisien MA menggambarkan hubungan nilai saat ini dengan error pada beberapa periode sebelumnya setelah differencing, sedangkan koefisien SMA menggambarkan hubungan dengan error pada periode musiman (12 bulan) sebelumnya. Model menghasilkan AIC sebesar 176.86, AICc sebesar 177, dan BIC sebesar 197.47 — nilai yang lebih kecil menunjukkan model lebih baik dibanding kandidat lain yang dievaluasi auto.arima().

4.6 Diagnostik Residual

Hipotesis

\[ H_0: \rho_1 = \rho_2 = \cdots = \rho_h = 0 \qquad H_1: \text{setidaknya terdapat satu } \rho_k \neq 0 \] dengan \(\rho_k\) merupakan autokorelasi residual pada lag ke-\(k\). Statistik Ljung-Box:

\[ Q = n(n+2) \sum_{k=1}^{h} \frac{\hat\rho_k^2}{n-k} \] Taraf signifikansi: \(\alpha = 5\%\)

checkresiduals(model_auto, plot = TRUE)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,3)(0,1,1)[12]
## Q* = 17.337, df = 20, p-value = 0.631
## 
## Model df: 4.   Total lags used: 24
lb_test <- Box.test(residuals(model_auto), lag = 24, fitdf = length(coef(model_auto)), type = "Ljung-Box")
tabel_lb <- data.frame(
  Statistik = c("X-squared", "df", "p-value"),
  Nilai = c(round(unname(lb_test$statistic), 4),
            unname(lb_test$parameter),
            round(lb_test$p.value, 4))
)
knitr::kable(tabel_lb, caption = "Hasil Uji Ljung-Box pada Residual Model")
Hasil Uji Ljung-Box pada Residual Model
Statistik Nilai
X-squared 17.3369
df 20.0000
p-value 0.6310

Keputusan dan kesimpulan: Diperoleh statistik \(Q^*\) sebesar 17.3369 dengan derajat bebas 20 dan p-value sebesar 0.631 pada 24 lag. Karena p-value > 0,05, maka \(H_0\) tidak ditolak. Dengan demikian, residual model `r ifelse(lb_test$p.value > 0.05, “dapat dianggap white noise, sehingga model SARIMA yang terbentuk telah memenuhi kriteria diagnostik residual.”, “masih menunjukkan indikasi autokorelasi, sehingga model perlu ditinjau ulang.”)

4.7 Evaluasi Akurasi (Data Latih vs Uji)

Data dibagi menjadi data latih (Januari 1959–Desember 1994) dan data uji (Januari 1995–Desember 1997).

co2_train <- window(co2, end = c(1994, 12))
co2_test  <- window(co2, start = c(1995, 1))

model_train <- auto.arima(co2_train, seasonal = TRUE)
forecast_test <- forecast(model_train, h = length(co2_test))
tabel_akurasi <- round(as.data.frame(accuracy(forecast_test, co2_test)), 4)
knitr::kable(tabel_akurasi, caption = "Ukuran Akurasi Forecast: Data Latih vs Data Uji")
Ukuran Akurasi Forecast: Data Latih vs Data Uji
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.0169 0.2832 0.2256 0.0050 0.0675 0.1829 0.0053 NA
Test set 0.4133 0.5226 0.4470 0.1139 0.1233 0.3624 0.2939 0.3963 
plot(forecast_test, main = "Forecast vs Data Aktual (1995-1997)")
lines(co2_test, col = "red")
legend("topleft", legend = c("Forecast", "Aktual"), col = c("blue", "red"), lty = 1)

Interpretasi: Pada data uji periode 1995-1997, model menghasilkan RMSE sebesar 0.5226 ppm, MAE sebesar 0.447 ppm, dan MAPE sebesar 0.1233%. RMSE dan MAE mengukur rata-rata besar kesalahan prediksi dalam satuan ppm, sedangkan MAPE menyatakannya dalam persentase relatif terhadap nilai aktual. Nilai-nilai tersebut perlu dibandingkan dengan skala data (CO2 berkisar sekitar 300-370 ppm) untuk menilai apakah kesalahan tergolong kecil atau besar secara substantif — bukan dinilai hanya dari besar-kecilnya angka semata.

4.8 Forecast ke Depan

forecast_depan <- forecast(model_auto, h = 24)
plot(forecast_depan, main = "Forecast Konsentrasi CO2 24 Bulan ke Depan")

Interpretasi: Model SARIMA(\(0,1,3\))(\(0,1,1\))\(_{12}\) digunakan untuk memproyeksikan konsentrasi CO2 24 bulan ke depan berdasarkan pola historis (trend dan musiman) yang telah dipelajari model. Interval prediksi melebar seiring bertambahnya horizon peramalan, menandakan ketidakpastian yang meningkat. Nilai forecast spesifik per bulan tidak dituliskan di sini karena tidak ditampilkan sebagai tabel di laporan ini — bila diperlukan, dapat ditambahkan lewat as.data.frame(forecast_depan).

5. Model Akhir dan Interpretasi

Berdasarkan seluruh tahapan analisis, model yang dipilih adalah SARIMA(\(0,1,3\))(\(0,1,1\))\(_{12}\), menggunakan differencing non-musiman sebanyak \(d=1\) kali dan musiman \(D=1\) kali dengan periode 12 bulan.

Diagnostik residual melalui uji Ljung-Box menghasilkan \(Q^*=17.3369\), derajat bebas 20, dan p-value 0.631. Karena p-value > 0,05, tidak terdapat bukti cukup adanya autokorelasi residual, sehingga residual dapat dianggap white noise.

Pada evaluasi data uji periode 1995-1997, model menghasilkan RMSE sebesar 0.5226 ppm, MAE sebesar 0.447 ppm, dan MAPE sebesar 0.1233%. Model kemudian digunakan untuk menghasilkan forecast konsentrasi CO2 selama 24 bulan ke depan

6. Kesimpulan

Data konsentrasi CO2 atmosfer Mauna Loa periode 1959-1997 terdiri atas 468 observasi bulanan dan menunjukkan trend naik jangka panjang serta pola musiman tahunan. Uji ADF pada data awal menghasilkan statistik Dickey-Fuller -2.8299 dengan p-value 0.2269, sehingga \(H_0\) tidak ditolak dan data pada level dinyatakan belum stasioner.

Setelah differencing non-musiman (\(d=1\)) dan musiman (\(D=1\), lag 12), uji ADF menghasilkan statistik -8.9846 dengan p-value 0.01, sehingga \(H_0\) ditolak, dan data hasil differencing dinyatakan stasioner.

Model SARIMA(\(0,1,3\))(\(0,1,1\))\(_{12}\) terpilih melalui auto.arima(), dengan uji Ljung-Box menghasilkan \(Q^*=17.3369\) dan p-value 0.631, sehingga residual dapat dianggap white noise. Pada data uji 1995-1997, model menghasilkan RMSE 0.5226 ppm, MAE 0.447 ppm, dan MAPE 0.1233%. Model selanjutnya digunakan untuk forecast konsentrasi CO2 24 bulan ke depan berdasarkan pola historis data, dengan catatan ketidakpastian yang membesar seiring horizon proyeksi dan tanpa memperhitungkan potensi perubahan struktural (misalnya kebijakan iklim) di luar pola historis.