1 Pendahuluan

1.1 Latar Belakang

Transportasi udara merupakan salah satu moda transportasi yang memiliki peran strategis dalam mendukung konektivitas antardaerah di Indonesia sebagai negara kepulauan. Jumlah penumpang pesawat domestik menjadi salah satu indikator penting dalam mengukur tingkat mobilitas masyarakat serta kinerja sektor penerbangan nasional.

Berdasarkan data Badan Pusat Statistik (BPS), rata-rata jumlah penumpang pesawat domestik pada periode sebelum pandemi (2017–2019) mencapai sekitar 7,2 juta penumpang per bulan, dengan puncak tertinggi pada Juli 2018 yang mencapai 8,95 juta penumpang. Namun, pandemi COVID-19 yang melanda Indonesia sejak awal 2020 menyebabkan penurunan drastis hingga titik terendah pada Mei 2020 dengan hanya 87.008 penumpang — turun lebih dari 98% dibandingkan kondisi normal. Seiring pelonggaran pembatasan sosial, jumlah penumpang berangsur pulih, meskipun rata-rata periode 2022–2025 masih berada di sekitar 4,97 juta penumpang per bulan.

Fluktuasi yang signifikan ini menunjukkan bahwa data penumpang pesawat domestik memiliki pola musiman yang kuat, terutama dipengaruhi oleh musim liburan, hari raya, dan faktor-faktor eksternal seperti pandemi. Oleh karena itu, diperlukan suatu metode peramalan deret waktu yang mampu menangkap pola musiman sekaligus tren jangka panjang secara akurat.

1.2 Tujuan Analisis

Laporan ini bertujuan untuk:

  1. Mengeksplorasi pola dan karakteristik data deret waktu penumpang pesawat domestik (tren, musiman, stasioneritas).
  2. Mengidentifikasi pola musiman yang terdapat dalam data.
  3. Menentukan parameter model SARIMA yang sesuai melalui analisis ACF dan PACF.
  4. Membandingkan beberapa model kandidat SARIMA berdasarkan kriteria AIC dan BIC.
  5. Mengevaluasi kebaikan model melalui uji diagnostik residual dan akurasi peramalan (MAPE, MSE).
  6. Melakukan peramalan jumlah penumpang pesawat domestik untuk 12 bulan ke depan (2026).

1.3 Sumber Data

Data yang digunakan dalam analisis ini adalah data jumlah penumpang pesawat domestik bulanan yang bersumber dari Badan Pusat Statistik (BPS) Republik Indonesia, mencakup 5 bandar udara utama di Indonesia pada periode Januari 2017 – Desember 2025 (108 observasi). Data ini mencerminkan volume penumpang pada rute domestik dan mencakup periode pandemi COVID-19 (2020–2021) yang menjadi gangguan struktural (structural break) pada deret waktu.

2 Penjelasan Teori

2.1 Deret Waktu (Time Series)

Deret waktu adalah sekumpulan data yang dikumpulkan atau dicatat berdasarkan urutan waktu dengan interval yang sama (misalnya harian, bulanan, atau tahunan). Tujuan analisis deret waktu adalah memahami pola yang terkandung dalam data historis dan menggunakannya untuk membuat peramalan (forecasting) di masa mendatang.

Secara umum, deret waktu dapat mengandung empat komponen utama:

  • Tren (Trend): Kecenderungan data meningkat atau menurun secara jangka panjang.
  • Musiman (Seasonal): Pola yang berulang secara periodik dalam interval waktu tertentu (misalnya setiap 12 bulan untuk data bulanan).
  • Siklus (Cycle): Fluktuasi jangka panjang yang tidak memiliki periode tetap.
  • Residual/Irregular: Variasi acak yang tidak dapat dijelaskan oleh ketiga komponen di atas.

2.2 Model ARIMA

ARIMA (Autoregressive Integrated Moving Average) adalah model deret waktu yang menggabungkan tiga komponen:

  • AR (Autoregressive) orde \(p\): nilai saat ini dipengaruhi oleh \(p\) nilai masa lalu.
  • I (Integrated) orde \(d\): proses differencing sebanyak \(d\) kali untuk mencapai stasioneritas.
  • MA (Moving Average) orde \(q\): nilai saat ini dipengaruhi oleh \(q\) nilai kesalahan (error) masa lalu.

Model ARIMA dinotasikan sebagai ARIMA(p, d, q) dan hanya cocok untuk data yang tidak memiliki pola musiman. Ketika data memiliki komponen musiman yang signifikan — seperti data penumpang pesawat bulanan — diperlukan perluasan model menjadi SARIMA.

2.3 Model SARIMA

SARIMA (Seasonal Autoregressive Integrated Moving Average) merupakan perluasan dari model ARIMA yang mampu menangkap pola musiman dalam data deret waktu. Model SARIMA dinotasikan sebagai:

\[\text{SARIMA}(p, d, q)(P, D, Q)[s]\]

2.3.1 Arti Setiap Parameter

Parameter Jenis Keterangan
\(p\) Non-musiman Orde AR non-musiman: pengaruh \(p\) nilai masa lalu terdekat
\(d\) Non-musiman Jumlah differencing non-musiman untuk mencapai stasioneritas
\(q\) Non-musiman Orde MA non-musiman: pengaruh \(q\) error masa lalu terdekat
\(P\) Musiman Orde AR musiman: pengaruh nilai pada lag musiman ke-\(1\) s.d. \(P\)
\(D\) Musiman Jumlah differencing musiman
\(Q\) Musiman Orde MA musiman: pengaruh error pada lag musiman ke-\(1\) s.d. \(Q\)
\(s\) Periode Panjang satu siklus musiman (\(s = 12\) untuk data bulanan)

Pada analisis ini digunakan \(s = 12\) karena data bersifat bulanan sehingga satu siklus musiman terdiri dari 12 bulan.

2.4 Identifikasi Parameter: Membaca Plot ACF dan PACF

Identifikasi orde model SARIMA dilakukan dengan mengamati pola pada plot ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function) setelah proses differencing. Terdapat dua pola utama yang diperhatikan:

  • Cut off (terpotong): Nilai ACF/PACF langsung tidak signifikan (masuk dalam batas interval kepercayaan) setelah lag tertentu.
  • Dies down (meluruh perlahan): Nilai ACF/PACF menurun secara bertahap dan tidak langsung terpotong.

2.4.1 Aturan Identifikasi Komponen Non-Musiman

Kondisi ACF Kondisi PACF Model yang Diindikasikan
Cut off setelah lag \(q\) Dies down MA(\(q\))
Dies down Cut off setelah lag \(p\) AR(\(p\))
Dies down Dies down ARMA(\(p\), \(q\))

2.4.2 Aturan Identifikasi Komponen Musiman (lag kelipatan \(s\))

Berdasarkan materi perkuliahan (Darmawan, 2024), pola pada lag musiman (lag \(s\), \(2s\), \(3s\), …) diidentifikasi sebagai berikut:

Model Musiman ACF pada Lag Musiman PACF pada Lag Musiman
SAR(\(P\)) Dies down Cut off setelah lag \(Ps\)
SMA(\(Q\)) Cut off setelah lag \(Qs\) Dies down
SARMA(\(P\), \(Q\)) Dies down Dies down

2.4.3 Cara Membaca ACF & PACF pada Analisis Ini

Setelah differencing gabungan (\(d = 1\), \(D = 1\)), pengamatan dilakukan pada:

  1. Lag-lag kecil (lag 1, 2, 3, …): untuk menentukan \(p\) dan \(q\) (komponen non-musiman).
  2. Lag-lag kelipatan 12 (lag 12, 24, 36, …): untuk menentukan \(P\) dan \(Q\) (komponen musiman).

Apabila pada lag-lag kecil, ACF cut off setelah lag 1 → mengindikasikan MA(1) (\(q = 1\)). Apabila PACF cut off setelah lag 1 → mengindikasikan AR(1) (\(p = 1\)). Pada lag musiman (kelipatan 12), apabila ACF cut off setelah lag 12 dan PACF dies down → mengindikasikan SMA(1) (\(Q = 1\)). Sebaliknya jika PACF cut off setelah lag 12 → SAR(1) (\(P = 1\)).

Berdasarkan pola-pola tersebut, beberapa model kandidat SARIMA diajukan dan dibandingkan menggunakan kriteria AIC (Akaike Information Criterion) dan BIC (Bayesian Information Criterion).

2.5 Evaluasi Model

Setelah model terbaik dipilih, dilakukan evaluasi melalui dua tahap:

1. Uji Diagnostik Residual

Model yang baik menghasilkan residual yang bersifat white noise — tidak mengandung autokorelasi tersisa. Pengujian dilakukan menggunakan uji Ljung-Box dengan hipotesis:

  • \(H_0\): Residual merupakan white noise (tidak ada autokorelasi)
  • \(H_1\): Residual mengandung autokorelasi

Jika \(p\text{-value} > 0{,}05\), maka gagal tolak \(H_0\) → residual bersifat white noise → model layak digunakan.

2. Evaluasi Akurasi Peramalan

Akurasi peramalan diukur menggunakan:

  • MAPE (Mean Absolute Percentage Error): rata-rata persentase kesalahan absolut. \[\text{MAPE} = \frac{1}{n}\sum_{t=1}^{n}\left|\frac{Y_t - \hat{Y}_t}{Y_t}\right| \times 100\%\]

    Kriteria MAPE: < 10% = sangat akurat, 10–20% = baik, 20–50% = layak, > 50% = tidak akurat.

  • MSE (Mean Squared Error): rata-rata kuadrat kesalahan, sensitif terhadap kesalahan besar. \[\text{MSE} = \frac{1}{n}\sum_{t=1}^{n}(Y_t - \hat{Y}_t)^2\]

3 Import dan Persiapan Data

data_raw <- read.csv("C:/Users/Fathia S/Downloads/data video adw1.xlsx - data.csv", header = TRUE, sep = ",")
colnames(data_raw) <- trimws(colnames(data_raw))

datatable(data_raw, options = list(pageLength = 6, scrollX = TRUE),
          caption = "Tabel 1. Data Jumlah Penumpang Pesawat Domestik per Bulan")
datats <- ts(data_raw$Total, start = c(2017, 1), frequency = 12)
data   <- datats

n_obs <- length(data)
cat("Jumlah observasi :", n_obs, "\n")
## Jumlah observasi : 108
cat("Periode data     : Januari", start(data)[1], "-",
    month.name[end(data)[2]], end(data)[1], "\n")
## Periode data     : Januari 2017 - December 2025

Interpretasi: Data terdiri dari 108 observasi bulanan dengan frekuensi musiman 12 (data bulanan). Variabel yang dianalisis adalah Total, yaitu jumlah penumpang pesawat domestik per bulan.

4 Eksplorasi Data

4.1 Plot Data Aktual

autoplot(data) +
  labs(title = "Data Aktual Jumlah Penumpang Pesawat Domestik (2017–2025)",
       x = "Tahun", y = "Jumlah Penumpang") +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 12) +
  geom_line(color = "steelblue", linewidth = 0.9)

Interpretasi: Plot menunjukkan adanya tren dan pola musiman yang relatif stabil pada periode sebelum 2020. Pada awal 2020 terjadi penurunan drastis akibat pandemi COVID-19, diikuti oleh proses pemulihan (recovery) secara bertahap hingga mendekati atau melampaui level sebelum pandemi pada tahun-tahun terakhir. Pola ini menunjukkan deret waktu belum stasioner baik dalam rata-rata maupun terdapat unsur musiman.

4.2 Plot Musiman (Seasonal Plot)

ggseasonplot(data, year.labels = TRUE, continuous = TRUE) +
  labs(title = "Seasonal Plot - Jumlah Penumpang Domestik per Bulan",
       y = "Jumlah Penumpang") +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 11)

Interpretasi: Seasonal plot memperlihatkan pola musiman yang berulang setiap tahun, dengan kenaikan pada bulan-bulan tertentu (misalnya menjelang liburan/lebaran/akhir tahun) dan penurunan pada bulan lain. Garis tahun 2020–2021 tampak jauh di bawah garis tahun lainnya, menegaskan dampak pandemi.

4.3 Subseries Plot

ggsubseriesplot(data) +
  labs(title = "Subseries Plot - Rata-rata Penumpang per Bulan",
       y = "Jumlah Penumpang") +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 10)

Interpretasi: Garis horizontal biru pada setiap panel menunjukkan rata-rata jumlah penumpang untuk bulan tersebut sepanjang periode pengamatan. Terlihat variasi rata-rata antarbulan yang mengindikasikan adanya komponen musiman dalam data.

4.4 ACF dan PACF Data Asli

par(mfrow = c(1, 2))
Acf(data,  lag.max = 48, main = "ACF Data Asli")
Pacf(data, lag.max = 48, main = "PACF Data Asli")

par(mfrow = c(1, 1))

Interpretasi: Plot ACF menunjukkan penurunan (decay) yang lambat, mengindikasikan data tidak stasioner dalam rata-rata. Pola gelombang berulang pada lag kelipatan 12 mengindikasikan adanya komponen musiman. Hal ini menjadi dasar perlunya proses differencing.

5 Uji Stasioneritas

5.1 Data Asli

adf_raw <- adf.test(data)
adf_raw
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data
## Dickey-Fuller = -1.7501, Lag order = 4, p-value = 0.6803
## alternative hypothesis: stationary

Interpretasi: Nilai p-value uji ADF sebesar 0.6803 (> 0,05), sehingga gagal tolak H0, artinya data belum stasioner. Karena belum stasioner, dilakukan proses differencing.

5.2 Differencing Non-Musiman (d=1, D=0)

dif_original <- diff(data)

autoplot(dif_original) +
  labs(title = "Differencing Non-Musiman (d=1, D=0)", x = "Tahun", y = "Selisih") +
  theme_minimal(base_size = 12) +
  geom_line(color = "darkorange")

adf_d1 <- adf.test(dif_original)
adf_d1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dif_original
## Dickey-Fuller = -5.6469, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Interpretasi: Setelah differencing satu kali (d=1), p-value ADF menjadi 0.01. Nilai ini ≤ 0,05 sehingga data sudah stasioner dalam rata-rata setelah differencing non-musiman, namun pola musiman pada plot masih dapat terlihat.

5.3 Differencing Musiman (d=0, D=1)

dif_season <- diff(data, lag = 12)

autoplot(dif_season) +
  labs(title = "Differencing Musiman (d=0, D=1)", x = "Tahun", y = "Selisih Musiman") +
  theme_minimal(base_size = 12) +
  geom_line(color = "seagreen")

adf_D1 <- adf.test(dif_season)
adf_D1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dif_season
## Dickey-Fuller = -2.1305, Lag order = 4, p-value = 0.5225
## alternative hypothesis: stationary

Interpretasi: Differencing musiman (lag 12) menghasilkan p-value ADF sebesar 0.5225. Hasil ini menunjukkan differencing musiman saja belum cukup untuk mencapai stasioneritas.

5.4 Differencing Gabungan (d=1, D=1)

dif_oriseason <- diff(dif_original, lag = 12)

autoplot(dif_oriseason) +
  labs(title = "Differencing Gabungan Non-Musiman & Musiman (d=1, D=1)",
       x = "Tahun", y = "Selisih") +
  theme_minimal(base_size = 12) +
  geom_line(color = "firebrick")

adf_dD1 <- adf.test(dif_oriseason)
adf_dD1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dif_oriseason
## Dickey-Fuller = -5.0312, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Interpretasi: Setelah differencing gabungan (d=1, D=1), p-value ADF menjadi 0.01 (≤ 0,05), sehingga H0 ditolak dan data dinyatakan stasioner. Berdasarkan hasil ini, ordo differencing d=1 dan D=1 dipilih sebagai dasar identifikasi model SARIMA.

6 Identifikasi Orde Model (ACF & PACF setelah Differencing d=1, D=1)

conf_bound <- qnorm(0.975) / sqrt(length(dif_oriseason))

par(mfrow = c(2, 1))

acf_res <- acf(dif_oriseason, lag.max = 48, plot = FALSE)
plot(acf_res$acf, type = "h", xaxt = "n",
     xlab = "Lag", ylab = "ACF", main = "ACF Data Setelah Differencing (d=1, D=1)")
axis(1, at = 1:48, labels = 1:48)
abline(h = c(0, conf_bound, -conf_bound), col = c("black", "blue", "blue"), lty = c(1, 2, 2))
abline(v = seq(12, 48, by = 12), col = rgb(1, 0, 0, 0.4), lty = 2)

pacf_res <- pacf(dif_oriseason, lag.max = 48, plot = FALSE)
plot(as.vector(pacf_res$acf), type = "h", xaxt = "n",
     xlab = "Lag", ylab = "PACF", main = "PACF Data Setelah Differencing (d=1, D=1)")
axis(1, at = 1:48, labels = 1:48)
abline(h = c(0, conf_bound, -conf_bound), col = c("black", "blue", "blue"), lty = c(1, 2, 2))
abline(v = seq(12, 48, by = 12), col = rgb(1, 0, 0, 0.4), lty = 2)

par(mfrow = c(1, 1))

Interpretasi: Garis biru putus-putus menunjukkan batas signifikansi (interval kepercayaan 95%). Lag yang melewati batas ini dianggap signifikan.

  • Pada ACF, lag signifikan pada lag-lag awal (non-musiman) mengindikasikan komponen MA non-musiman (q). Lag signifikan pada kelipatan 12 (lag 12, 24, …) mengindikasikan komponen MA musiman (Q).
  • Pada PACF, pola serupa pada lag-lag awal mengindikasikan komponen AR non-musiman (p), sedangkan signifikansi pada kelipatan 12 mengindikasikan komponen AR musiman (P).

Berdasarkan pola decay/cut-off pada kedua plot tersebut, beberapa model kandidat SARIMA(p,d,q)(P,D,Q)[12] dengan d=1 dan D=1 diajukan pada bagian berikut, ditambah hasil pencarian otomatis auto.arima() sebagai pembanding.

7 Estimasi Model Kandidat SARIMA

M1    <- Arima(data, order = c(0,1,1), seasonal = list(order = c(0,1,1), period = 12))
M2    <- Arima(data, order = c(1,1,0), seasonal = list(order = c(1,1,0), period = 12))
M3    <- Arima(data, order = c(0,1,1), seasonal = list(order = c(1,1,0), period = 12))
M4    <- Arima(data, order = c(1,1,0), seasonal = list(order = c(0,1,1), period = 12))
M5    <- Arima(data, order = c(1,1,1), seasonal = list(order = c(1,1,1), period = 12))
Mauto <- auto.arima(data, stepwise = FALSE, approximation = FALSE)

models      <- list(M1, M2, M3, M4, M5, Mauto)
model_names <- c("SARIMA(0,1,1)(0,1,1)[12]",
                 "SARIMA(1,1,0)(1,1,0)[12]",
                 "SARIMA(0,1,1)(1,1,0)[12]",
                 "SARIMA(1,1,0)(0,1,1)[12]",
                 "SARIMA(1,1,1)(1,1,1)[12]",
                 paste0("auto.arima -> ", arimaorder(Mauto)[1], ",",
                        arimaorder(Mauto)[2], ",", arimaorder(Mauto)[3], "x",
                        arimaorder(Mauto)[4], ",", arimaorder(Mauto)[5], ",",
                        arimaorder(Mauto)[6], "[12]"))

aic_df <- data.frame(
  Model = model_names,
  AIC   = round(sapply(models, AIC), 2),
  BIC   = round(sapply(models, BIC), 2)
)

datatable(aic_df, options = list(dom = 't'),
          caption = "Tabel 2. Perbandingan AIC dan BIC Model Kandidat") %>%
  formatStyle(columns = "AIC", backgroundColor = styleEqual(min(aic_df$AIC), "lightgreen"))
best_idx   <- which.min(aic_df$AIC)
best_model <- models[[best_idx]]

cat("Model terbaik berdasarkan AIC terendah:", aic_df$Model[best_idx],
    "\nAIC =", aic_df$AIC[best_idx],
    "| BIC =", aic_df$BIC[best_idx], "\n")
## Model terbaik berdasarkan AIC terendah: SARIMA(1,1,1)(1,1,1)[12] 
## AIC = 2864.06 | BIC = 2876.83

Interpretasi: Tabel di atas membandingkan enam model kandidat SARIMA berdasarkan nilai AIC dan BIC (semakin kecil semakin baik). Model dengan AIC terendah, yaitu SARIMA(1,1,1)(1,1,1)[12] (AIC = 2864.06, BIC = 2876.83), dipilih sebagai model terbaik untuk digunakan pada analisis selanjutnya.

8 Ringkasan Model Terbaik

summary(best_model)
## Series: data 
## ARIMA(1,1,1)(1,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1
##       0.8276  -0.8934  0.2604  -0.9999
## s.e.  0.4644   0.4045  0.1201   0.1486
## 
## sigma^2 = 5.461e+11:  log likelihood = -1427.03
## AIC=2864.06   AICc=2864.73   BIC=2876.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 6788.919 678356.8 417328.6 -11.72763 22.71989 0.3132201 0.0141989

Interpretasi: Output di atas menampilkan estimasi koefisien AR/MA (musiman dan non-musiman) beserta standard error-nya, serta ukuran akurasi in-sample seperti RMSE, MAE, dan MAPE. Koefisien dengan nilai |koefisien| jauh lebih besar dari standard error-nya (rasio > 2) dapat dianggap signifikan secara statistik.

9 Uji Diagnostik Residual

checkresiduals(best_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(1,1,1)[12]
## Q* = 26.64, df = 18, p-value = 0.08601
## 
## Model df: 4.   Total lags used: 22

Interpretasi: Uji Ljung-Box menghasilkan p-value sebesar 0.08601.

  • Jika p-value > 0,05, maka gagal tolak H0, artinya residual bersifat white noise (tidak ada autokorelasi tersisa) — model dianggap layak (adequate).
  • Jika p-value ≤ 0,05, maka residual masih memiliki autokorelasi, sehingga model perlu dievaluasi ulang/diperbaiki.

Pada kasus ini, p-value sebesar 0.08601> 0.05, sehingga residual model memenuhi asumsi white noise (model layak digunakan). Histogram residual pada plot di atas juga menunjukkan apakah residual mendekati distribusi normal.

10 Visualisasi Nilai Aktual vs Fitted

df_fit <- data.frame(
  Time   = as.numeric(time(data)),
  Aktual = as.numeric(data),
  Fitted = as.numeric(fitted(best_model))
)

p_fit <- ggplot(df_fit, aes(x = Time)) +
  geom_line(aes(y = Aktual, color = "Aktual"), linewidth = 0.9) +
  geom_line(aes(y = Fitted, color = "Fitted"), linewidth = 0.9) +
  scale_color_manual(values = c("Aktual" = "steelblue", "Fitted" = "red")) +
  scale_y_continuous(labels = comma) +
  labs(title = "Perbandingan Nilai Aktual vs Fitted",
       x = "Tahun", y = "Jumlah Penumpang", color = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

ggplotly(p_fit)

Interpretasi: Plot interaktif di atas menunjukkan bahwa nilai fitted (merah) secara umum mengikuti pola data aktual (biru), termasuk pola musiman dan penurunan tajam pada periode pandemi. Kedekatan kedua garis mengindikasikan model SARIMA terpilih cukup baik dalam menangkap pola historis data. Arahkan kursor pada grafik untuk melihat nilai pada titik waktu tertentu.

11 Peramalan 12 Bulan ke Depan

fc <- forecast(best_model, h = 12)

fc_df <- data.frame(
  Periode  = format(seq(as.Date(paste0(end(data)[1], "-", end(data)[2], "-01")),
                         by = "month", length.out = 13)[-1], "%b %Y"),
  Forecast = round(as.numeric(fc$mean), 0),
  Lo80     = round(as.numeric(fc$lower[,1]), 0),
  Hi80     = round(as.numeric(fc$upper[,1]), 0),
  Lo95     = round(as.numeric(fc$lower[,2]), 0),
  Hi95     = round(as.numeric(fc$upper[,2]), 0)
)

datatable(fc_df, options = list(pageLength = 12, dom = 't'),
          caption = "Tabel 3. Hasil Peramalan 12 Bulan ke Depan") %>%
  formatCurrency(columns = 2:6, currency = "", interval = 3, mark = ",", digits = 0)
plot(fc,
     main = "Forecast Penumpang Domestik 2026",
     xlab = "Tahun", ylab = "Jumlah Penumpang")

Interpretasi: Berdasarkan tabel dan grafik peramalan, hasil peramalan untuk 12 bulan mendatang (Januari–Desember 2026) menunjukkan prediksi titik (forecast) yang berkisar antara 4.013.081 hingga 5.224.828 penumpang per bulan, beserta interval kepercayaan 80% (Lo80–Hi80) dan 95% (Lo95–Hi95) yang digambarkan sebagai pita biru dan abu-abu pada grafik. Pola musiman pada periode peramalan mengikuti pola historis, dengan nilai yang relatif lebih rendah pada awal hingga pertengahan tahun (Februari–Maret di kisaran 4,0–4,2 juta) dan meningkat pada bulan Juli serta Desember (masing-masing 4.991.996 dan 5.224.828), sesuai karakteristik musiman yang teridentifikasi sebelumnya. Lebar interval kepercayaan yang semakin membesar untuk horizon yang lebih jauh menunjukkan meningkatnya ketidakpastian peramalan jangka panjang.

12 Evaluasi Akurasi Peramalan (Rolling Test: h = 3, 6, 12)

test_lengths <- c(3, 6, 12)
eval_list    <- list()

for (n_test in test_lengths) {
  train <- head(data, -n_test)
  test  <- tail(data, n_test)

  mdl  <- Arima(train,
                order    = best_model$arma[c(1,6,2)],
                seasonal = list(order  = best_model$arma[c(3,7,4)],
                                period = frequency(data)))
  fcast <- forecast(mdl, h = n_test)$mean

  eval_list[[as.character(n_test)]] <- data.frame(
    h    = n_test,
    MAPE = round(mean(abs((test - fcast) / test)) * 100, 4),
    MSE  = round(mean((test - fcast)^2), 2)
  )
}

eval_df <- do.call(rbind, eval_list)

datatable(eval_df, options = list(dom = 't'),
          caption = "Tabel 4. Evaluasi Akurasi Peramalan pada Berbagai Horizon")

Interpretasi: Tabel di atas menunjukkan akurasi model pada data uji (out-of-sample) dengan horizon peramalan h = 3, 6, dan 12 bulan terakhir.

  • MAPE (Mean Absolute Percentage Error) mengukur rata-rata persentase kesalahan absolut. Nilai MAPE di bawah 10% umumnya dikategorikan sebagai peramalan yang sangat baik (highly accurate), 10–20% baik, 20–50% layak, dan di atas 50% tidak akurat.
  • MSE (Mean Squared Error) mengukur rata-rata kuadrat kesalahan; satuannya sama dengan kuadrat satuan data sehingga lebih sensitif terhadap kesalahan besar (outlier).

Pada horizon h=3, MAPE sebesar 3.179%, tergolong sangat akurat. Pada horizon h=12, MAPE sebesar 3.7303%, masih tergolong sangat akurat.

13 Visualisasi Forecast vs Aktual (H = 3)

n          <- length(data)
train      <- head(data, -3)
test_3     <- tail(data, 3)
mdl3       <- Arima(train,
                    order    = best_model$arma[c(1,6,2)],
                    seasonal = list(order  = best_model$arma[c(3,7,4)],
                                    period = frequency(data)))
fc3        <- forecast(mdl3, h = 3)$mean
fitted3    <- fitted(mdl3)

df_plot <- data.frame(
  Time     = 1:n,
  Aktual   = as.numeric(data),
  Fitted   = c(as.numeric(fitted3), rep(NA, 3)),
  Forecast = c(rep(NA, n-3), as.numeric(fc3)),
  IsTest   = c(rep(FALSE, n-3), rep(TRUE, 3))
)

ggplot(df_plot) +
  geom_line(data = df_plot[!df_plot$IsTest, ],
            aes(x=Time, y=Aktual,   color="Data Asli"),     linewidth=1) +
  geom_line(data = df_plot[df_plot$IsTest, ],
            aes(x=Time, y=Aktual,   color="3 Data Testing"), linewidth=1) +
  geom_line(data = df_plot[!df_plot$IsTest, ],
            aes(x=Time, y=Fitted,   color="Fitted"),         linewidth=1) +
  geom_line(data = df_plot[df_plot$IsTest, ],
            aes(x=Time, y=Forecast, color="Forecast"),       linewidth=1) +
  scale_color_manual(values=c("Data Asli"="steelblue",
                              "3 Data Testing"="purple",
                              "Fitted"="red",
                              "Forecast"="green4")) +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size=13) +
  labs(title="SARIMA: Forecast vs Aktual (H=3)",
       x="Time Index", y="Jumlah Penumpang",
       color=NULL) +
  theme(legend.position="bottom")

n      <- length(data)
train  <- head(data, -3)
test_3 <- tail(data, 3)

mdl3    <- Arima(train,
                  order    = best_model$arma[c(1,6,2)],
                  seasonal = list(order  = best_model$arma[c(3,7,4)],
                                  period = frequency(data)))
fc3     <- forecast(mdl3, h = 3)$mean
fitted3 <- fitted(mdl3)

plot_df <- data.frame(
  Time     = 1:n,
  Original = as.numeric(data),
  Fitted   = c(as.numeric(fitted3), rep(NA, 3)),
  Forecast = c(rep(NA, n - 3), as.numeric(fc3)),
  IsLast3  = c(rep(FALSE, n - 3), rep(TRUE, 3))
)

p_h3 <- ggplot() +
  geom_line(data = plot_df[1:(n-3), ],
            aes(x = Time, y = Original, color = "Data Asli"), linewidth = 1) +
  geom_line(data = plot_df[plot_df$IsLast3, ],
            aes(x = Time, y = Original, color = "3 Data Terakhir (Aktual)"), linewidth = 1) +
  geom_line(data = plot_df[1:(n-3), ],
            aes(x = Time, y = Fitted, color = "Fitted Values"), linewidth = 1) +
  geom_line(data = plot_df,
            aes(x = Time, y = Forecast, color = "Forecast"), linewidth = 1) +
  scale_color_manual(values = c("Data Asli" = "steelblue",
                                 "3 Data Terakhir (Aktual)" = "purple",
                                 "Fitted Values" = "red",
                                 "Forecast" = "green4")) +
  scale_y_continuous(labels = comma) +
  coord_cartesian(xlim = c(n - 14, n)) +
  theme_minimal(base_size = 12) +
  labs(title = "SARIMA Forecast vs Aktual (H=3) - 15 Data Terakhir",
       x = "Time Index", y = "Jumlah Penumpang", color = NULL) +
  theme(legend.position = "bottom")

ggplotly(p_h3)

Interpretasi: Grafik di atas memfokuskan pada 15 observasi terakhir, membandingkan:

  • Data Asli (biru): nilai aktual pada data latih.
  • 3 Data Terakhir (ungu): nilai aktual pada periode uji.
  • Fitted Values (merah): nilai dugaan model pada data latih.
  • Forecast (hijau): nilai prediksi model untuk 3 periode uji.

Kedekatan antara garis Forecast (hijau) dan 3 Data Terakhir (ungu) mengindikasikan akurasi peramalan jangka pendek model. Semakin berdekatan kedua garis tersebut, semakin baik kinerja model dalam memprediksi nilai aktual pada horizon dekat.

14 Kesimpulan

Berdasarkan analisis deret waktu terhadap data jumlah penumpang pesawat domestik di Indonesia periode 2017–2025, diperoleh beberapa kesimpulan sebagai berikut:

  1. Karakteristik Data
    Data jumlah penumpang pesawat domestik menunjukkan adanya pola tren dan musiman bulanan yang cukup jelas. Selain itu, terdapat penurunan yang sangat signifikan pada periode pandemi COVID-19 tahun 2020–2021 yang menyebabkan terjadinya gangguan struktural (structural break) pada deret waktu.

  2. Stasioneritas Data
    Hasil uji stasioneritas menggunakan Augmented Dickey-Fuller (ADF) menunjukkan bahwa data asli belum stasioner dengan nilai p-value sebesar 0,6803. Setelah dilakukan differencing non-musiman orde satu (d = 1) dan differencing musiman orde satu (D = 1), data menjadi stasioner dengan p-value sebesar 0,01, sehingga memenuhi asumsi untuk pemodelan SARIMA.

  3. Model Terbaik
    Berdasarkan proses identifikasi dan pemilihan model menggunakan kriteria Akaike Information Criterion (AIC) dan Bayesian Information Criterion (BIC), model terbaik yang diperoleh adalah SARIMA(1,1,1)(1,1,1)[12] dengan nilai AIC = 2864,06 dan BIC = 2876,83.

  4. Diagnostik Residual
    Hasil uji Ljung–Box menghasilkan p-value sebesar 0,08601. Karena nilai tersebut lebih besar dari taraf signifikansi 5%, maka gagal menolak hipotesis nol yang menyatakan bahwa residual tidak mengandung autokorelasi. Dengan demikian, residual model dapat dianggap bersifat white noise sehingga model SARIMA yang diperoleh dinilai memadai (adequate) untuk digunakan dalam peramalan.

  5. Akurasi Peramalan
    Evaluasi akurasi peramalan menunjukkan bahwa model memiliki kinerja yang sangat baik. Nilai Mean Absolute Percentage Error (MAPE) yang diperoleh sebesar 3,179% untuk horizon peramalan 3 periode (h = 3) dan 3,7303% untuk horizon peramalan 12 periode (h = 12). Nilai MAPE yang berada di bawah 10% menunjukkan bahwa hasil peramalan memiliki tingkat akurasi yang sangat baik.

  6. Implikasi Praktis
    Model SARIMA(1,1,1)(1,1,1)[12] dapat digunakan sebagai dasar dalam melakukan peramalan jumlah penumpang pesawat domestik untuk mendukung perencanaan kapasitas transportasi udara, pengambilan keputusan operasional oleh bandara dan maskapai penerbangan, serta penyusunan kebijakan terkait pengembangan sektor transportasi udara di Indonesia.

Secara keseluruhan, model SARIMA(1,1,1)(1,1,1)[12] mampu menangkap pola tren dan musiman pada data jumlah penumpang pesawat domestik Indonesia dengan baik. Hasil diagnostik residual dan evaluasi akurasi menunjukkan bahwa model telah memenuhi asumsi yang diperlukan serta menghasilkan performa peramalan yang baik, sehingga layak digunakan untuk peramalan jangka pendek hingga menengah.

15 Catatan Penggunaan File

File .Rmd ini bersifat reproducible: pastikan file data_pesawat.csv berada pada direktori kerja (working directory) yang sama dengan file .Rmd ini sebelum melakukan Knit. Pustaka R yang dibutuhkan: forecast, tseries, ggplot2, scales, DT, dan plotly.