1. Pendahuluan

Laporan ini menyajikan analisis komputasi statistika menggunakan metode Analisis Deret Waktu (Time Series Analysis) untuk meramalkan volume penumpang pesawat penerbangan internasional. Kemampuan meramalkan jumlah penumpang sangat krusial dalam industri penerbangan untuk optimalisasi alokasi sumber daya, manajemen armada, hingga perencanaan finansial jangka panjang.

Algoritma pemodelan yang digunakan dalam penelitian ini adalah ARIMA (Autoregressive Integrated Moving Average). Mengingat data penerbangan umumnya dipengaruhi oleh siklus liburan atau musim tertentu, kita akan memperluas model ini menjadi SARIMA (Seasonal ARIMA) agar mampu menangkap variasi musiman yang kompleks. Pemodelan metodologi Box-Jenkins atau ARIMA secara fundamental dibangun atas tiga komponen parameter utama, yaitu \(ARIMA(p, d, q)\) sebagai berikut:

  1. AR (Autoregressive) orde \(p\), Mengasumsikan bahwa nilai observasi saat ini dipengaruhi oleh nilai observasi pada periode waktu sebelumnya secara linier.\[Y_t = c + \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \dots + \phi_p Y_{t-p} + e_t\]

  2. I (Integrated) orde \(d\), menunjukkan jumlah tahapan diferensiasi (differencing) yang diperlukan untuk menghilangkan tren agar data mencapai stasioneritas (rata-rata dan varians konstan).

  3. MA (Moving Average) orde \(q\), mengasumsikan bahwa nilai observasi saat ini dipengaruhi oleh error (kejutan/sisaan) dari periode sebelumnya.\[Y_t = c + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \dots + \theta_q e_{t-q}\]

Untuk menangani data musiman, model ini diekstensi menjadi \(SARIMA(p,d,q)(P,D,Q)_S\), di mana huruf kapital merepresentasikan orde untuk komponen musiman, dan \(S\) adalah panjang siklus musiman (misalnya 12 untuk data bulanan).

2. Persiapan Library dan Data

Tahap awal dalam analisis komputasi ini adalah mempersiapkan lingkungan kerja R dengan memuat library yang relevan. Kita menggunakan tiga package utama,, yaitu forecast yang merupakan package komprehensif yang menjadi standar industri untuk pemodelan dan peramalan deret waktu, tseries yang digunakan untuk pengujian hipotesis formal terkait stasioneritas data (seperti Uji Augmented Dickey-Fuller), dan ggplot2 yang digunakan untuk menghasilkan visualisasi data yang lebih estetis dan informatif. Dataset yang dianalisis adalah dataset bawaan R, yaitu AirPassengers. Dataset ini merekam total penumpang penerbangan internasional bulanan dari Januari 1949 hingga Desember 1960.

# Memanggil library yang dibutuhkan
library(forecast)
library(tseries)
library(ggplot2)

# Memuat dataset bawaan R (Format sudah berupa objek Time Series / ts)
data("AirPassengers")
df_ts <- AirPassengers

# Menampilkan informasi struktur dan ringkasan data
str(df_ts)
##  Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...
summary(df_ts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   104.0   180.0   265.5   280.3   360.5   622.0

Data terdiri dari 144 observasi berformat time series (ts) yang dicatat setiap bulan selama 12 tahun. Nilai minimum penumpang adalah 104 ribu, sedangkan puncaknya mencapai 622 ribu penumpang.

3. Eksplorasi Data dan Dekomposisi

Sebelum masuk ke pemodelan, sangat penting untuk memahami perilaku fundamental dari data. Kita melakukan dekomposisi untuk memecah deret waktu menjadi tiga komponen dasar pendiri data sebagai berikut:

  1. Trend, merupaka pergerakan arah jangka panjang dari data (apakah secara makro naik, turun, atau datar).

  2. Seasonal, merupakan pola fluktuasi yang berulang secara periodik dan teratur dalam rentang waktu tertentu.

  3. Random/Remainder, merupakan fluktuasi acak atau noise yang tersisa setelah tren dan musiman dikeluarkan dari data.

Jika kita amati secara visual (yang akan kita buktikan lewat plot di bawah), amplitudo atau lebar ayunan musim pada data AirPassengers terlihat semakin melebar seiring bertambahnya waktu dan naiknya tren. Oleh karena itu, dekomposisi dengan model multiplicatif dinilai lebih merepresentasikan realitas data dibandingkan model aditif.

# Melakukan dekomposisi dengan asumsi model perkalian (multiplicative)
dekomposisi <- decompose(df_ts, type = "multiplicative")

# Memplot hasil dekomposisi klasik
plot(dekomposisi)

Berdasarkan plot dekomposisi multiplicatif di atas, deret waktu berhasil dipisahkan menjadi empat panel yang memberikan wawasan analitis sebagai berikut:

  1. Panel teratas menunjukkan plot data asli jumlah penumpang dari tahun 1949 hingga menjelang 1961. Secara visual, sangat jelas terlihat adanya tren yang terus menanjak dan pola fluktuasi yang berulang. Selain itu, lebar ayunan atau amplitudo fluktuasi semakin membesar seiring berjalannya waktu (semakin ke kanan semakin lebar). Karakteristik inilah yang menjadi landasan kuat mengapa kita menggunakan dekomposisi multiplikatif, bukan aditif.

  2. Panel kedua mengonfirmasi keberadaan tren positif yang sangat kuat. Garis tren bergerak naik secara mulus (hampir linier) dari angka di bawah 150 pada awal pengamatan hingga mendekati angka 450 pada akhir periode. Hal ini merepresentasikan pertumbuhan industri penerbangan yang masif. Secara statistika, pergerakan tren yang signifikan ini menjadi bukti tak terbantahkan bahwa data belum stasioner dalam rata-rata (sehingga akan membutuhkan proses differencing pada tahap selanjutnya).

  3. Panel ketiga mengekstrak pola musiman murni. Terlihat gelombang yang sangat teratur dan konsisten berulang setiap tahunnya. Karena menggunakan model multiplikatif, nilai indeks musiman berpusat di angka 1.0. Terdapat puncak-puncak lonjakan (dengan indeks lebih dari 1.2, atau setara dengan kenaikan 20% lebih dari rata-rata tren) yang biasanya terjadi pada pertengahan tahun/liburan musim panas, dan lembah penurunan (mendekati indeks 0.8) pada musim berikutnya.

  4. Panel paling bawah menampilkan komponen sisa atau error yang tidak bisa dijelaskan oleh tren dan musiman. Nilainya berfluktuasi secara acak di sekitar batas 1.0 (mayoritas berada pada rentang 0.90 hingga 1.10). Fluktuasi acak ini mengindikasikan bahwa proses dekomposisi telah cukup baik memisahkan pola utama, meskipun variasi volatilitas yang masih tersisa menunjukkan bahwa pemodelan menggunakan metode lanjutan seperti SARIMA sangat tepat untuk dilakukan.

4. Uji Stasioneritas Silang dan Proses Pembedaan (Differencing)

Syarat mutlak dalam pemodelan ARIMA (metodologi Box-Jenkins) adalah deret waktu harus bersifat stasioner murni (strictly stationary). Stasioneritas berarti bahwa parameter statistik dari deret waktu, khususnya rata-rata (mean) dan varians, konstan sepanjang waktu dan tidak memiliki tren.

Meskipun secara visual pada plot dekomposisi sebelumnya data terlihat memiliki tren naik yang kuat, kita wajib membuktikannya secara statistik formal. Untuk menghindari anomali komputasi di mana satu uji bisa keliru mendeteksi trend-stationarity sebagai stasioneritas murni, laporan ini menggunakan pendekatan uji silang dengan dua metode: Uji Augmented Dickey-Fuller (ADF) dan Uji Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

4.1. Uji Stasioneritas Data Asli

A. Uji Augmented Dickey-Fuller (ADF) Uji ADF berfokus pada pendeteksian keberadaan akar unit (unit root) dalam sebuah deret waktu.

  • Hipotesis: \(H_0\): Data mengandung akar unit (Data tidak stasioner). \(H_1\): Data tidak mengandung akar unit (Data stasioner).
  • Taraf Signifikansi (\(\alpha\)): 5% (0.05).
  • Kriteria Penolakan: Tolak \(H_0\) jika p-value < \(\alpha\).

B. Uji Kwiatkowski-Phillips-Schmidt-Shin (KPSS) Berbeda dengan ADF, uji KPSS membalikkan rumusan hipotesis sehingga menguji langsung asumsi stasioneritas di sekitar tren deterministik. * Hipotesis: \(H_0\): Data stasioner. \(H_1\): Data tidak stasioner. * Taraf Signifikansi (\(\alpha\)): 5% (0.05). * Kriteria Penolakan: Tolak \(H_0\) jika p-value < \(\alpha\).

# Melakukan uji ADF pada data asli
adf.test(df_ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_ts
## Dickey-Fuller = -7.3186, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
# Melakukan uji silang dengan KPSS pada data asli
kpss.test(df_ts)
## 
##  KPSS Test for Level Stationarity
## 
## data:  df_ts
## KPSS Level = 2.7395, Truncation lag parameter = 4, p-value = 0.01

Berdasarkan keluaran komputasi, terdapat anomali pada uji ADF yang menghasilkan p-value = 0.01 (Tolak \(H_0\)). Secara keliru, angka ini mengindikasikan bahwa data sudah stasioner. Hal ini terjadi karena fungsi adf.test pada paket tseries secara bawaan menyertakan tren waktu deterministik dalam regresinya, sehingga yang tertangkap adalah trend-stationarity, bukan stasioneritas murni. Namun, uji KPSS berhasil memvalidasi ketidakstasioneran tersebut secara akurat dengan menghasilkan p-value = 0.01. Karena nilai p-value (0.01) < \(\alpha\) (0.05), maka kita menolak \(H_0\). Kesimpulan akhirnya adalah data asli secara statistik terbukti tidak stasioner dalam rata-rata karena masih mengandung tren linier, sehingga tahapan pembedaan (differencing) diperlukan.

4.2. Proses Differencing dan Uji Stasioneritas Lanjutan

Untuk menstabilkan rata-rata, kita menerapkan differencing, yaitu proses matematis yang menghitung selisih antara nilai observasi saat ini dengan nilai observasi sebelumnya. Kita memanfaatkan fungsi komputasi ndiffs() untuk memverifikasi secara objektif berapa kali jumlah diferensiasi non-musiman yang dibutuhkan.

# Mengecek rekomendasi jumlah differencing untuk mencapai stasioneritas
ndiffs(df_ts)
## [1] 1
# Melakukan differencing orde 1 (d = 1) sesuai rekomendasi algoritmik
df_diff <- diff(df_ts)

# Menguji kembali stasioneritas pada data yang sudah di-differencing dengan metode KPSS
kpss.test(df_diff)
## 
##  KPSS Test for Level Stationarity
## 
## data:  df_diff
## KPSS Level = 0.014626, Truncation lag parameter = 4, p-value = 0.1

Fungsi ndiffs() mengonfirmasi bahwa data memerlukan pembedaan tingkat pertama (\(d=1\)). Setelah dilakukan diferensiasi dan diuji ulang menggunakan uji KPSS, nilai p-value meningkat tajam menjadi 0.1.Karena p-value (0.1) > \(\alpha\) (0.05), maka keputusan uji adalah Gagal Menolak \(H_0\). Kesimpulan akhirnya, data yang telah melalui diferensiasi orde pertama (df_diff) terbukti secara meyakinkan telah stasioner dalam rata-rata. Deret waktu ini kini telah memenuhi syarat fundamental pemodelan dan siap dieksplorasi lebih lanjut untuk identifikasi parameter melalui plot Autocorrelation Function (ACF) dan Partial Autocorrelation Function (PACF).

5. Identifikasi Model melalui Plot ACF dan PACF

Setelah mengamankan asumsi stasioneritas rata-rata melalui differencing, langkah selanjutnya adalah mengidentifikasi kandidat parameter Autoregressive (AR) dan Moving Average (MA). Identifikasi ini dilakukan secara visual menggunakan plot Autocorrelation Function (ACF) dan Partial Autocorrelation Function (PACF). Plot ACF membantu kita menentukan orde komponen MA (\(q\)), sedangkan plot PACF membantu mengidentifikasi orde komponen AR (\(p\)).

# Memplot deret waktu hasil differencing beserta fungsi autokorelasinya
ggtsdisplay(df_diff, main = "Data Differenced beserta Plot ACF dan PACF")

Berdasarkan visualisasi plot di atas, kita dapat membedah karakteristik data menjadi dua aspek utama sebagai berikut:

  1. Panel Atas (Data Hasil Differencing) menunjukkan bahwa data kini berfluktuasi di sekitar nilai nol, mengonfirmasi hilangnya tren linier makro. Namun, terlihat jelas adanya pola ayunan tajam yang berulang secara periodik, serta amplitudo (lebar simpangan) yang semakin membesar di tahun-tahun akhir pengamatan. Hal ini menandakan bahwa varians belum sepenuhnya stabil dan efek musiman masih sangat kental.

  2. Pada plot ACF (kiri bawah), terdapat lonjakan (spikes) positif yang sangat signifikan dan menonjol pada lag kelipatan 12 (yaitu lag 12, 24, dan 36). Lonjakan musiman ini turun dengan sangat lambat (slow decay), yang merupakan indikasi klasik bahwa data masih membutuhkan diferensiasi musiman (seasonal differencing).

  3. Pada plot PACF (kanan bawah), terlihat lonjakan signifikan pada beberapa lag awal (seperti lag 1, 2, 3) yang mengindikasikan keberadaan komponen AR non-musiman. Selain itu, terdapat lonjakan tajam pada lag 12, mengonfirmasi perlunya penanganan parameter musiman.

Kesimpulan dari identifikasi visual ini adalah pemodelan tidak bisa diselesaikan dengan ARIMA standar. Data menuntut penggunaan ekstensi Seasonal ARIMA (SARIMA) yang mengakomodasi diferensiasi ganda (non-musiman dan musiman).

6. Pembangunan Model SARIMA Optimal

Identifikasi manual melalui visualisasi ACF dan PACF sering kali memunculkan banyak kandidat model yang subjektif, terutama pada deret waktu dengan struktur musiman yang kompleks. Untuk memastikan ketepatan dan efisiensi pemodelan, laporan ini menggunakan pendekatan komputasional algoritmik melalui fungsi auto.arima(). Algoritma ini melakukan iterasi pencarian terhadap ruang kemungkinan parameter \((p,d,q)\) dan \((P,D,Q)\) secara otomatis, lalu memilih satu arsitektur model terbaik berdasarkan nilai Akaike Information Criterion (AIC) terkecil. Pendekatan ini secara matematis memastikan model memiliki akurasi prediksi tertinggi dengan tingkat kompleksitas parameter yang paling seimbang (menghindari overfitting).

# Mencari dan membangun model SARIMA terbaik secara komputasional
model_arima <- auto.arima(df_ts)

# Menampilkan ringkasan metrik evaluasi dan arsitektur parameter model terpilih
summary(model_arima)
## Series: df_ts 
## ARIMA(2,1,1)(0,1,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.5960  0.2143  -0.9819
## s.e.  0.0888  0.0880   0.0292
## 
## sigma^2 = 132.3:  log likelihood = -504.92
## AIC=1017.85   AICc=1018.17   BIC=1029.35
## 
## Training set error measures:
##                  ME     RMSE     MAE      MPE     MAPE     MASE        ACF1
## Training set 1.3423 10.84619 7.86754 0.420698 2.800458 0.245628 -0.00124847

Algoritma komputasi secara objektif menetapkan arsitektur ARIMA(2,1,1)(0,1,0)[12] sebagai model prediktif yang paling optimal berdasarkan nilai Kriteria Informasi Akaike (AIC) terkecil, yaitu sebesar 1017.85. Spesifikasi matematis dari model ini dapat dibedah sebagai berikut:

  1. Komponen Non-Musiman (2,1,1). Model ini memanfaatkan observasi dari dua periode sebelumnya (\(p=2\)) dengan nilai koefisien \(AR_1 = 0.5960\) dan \(AR_2 = 0.2143\). Model juga telah melalui satu tahapan diferensiasi tingkat pertama (\(d=1\)) untuk mereduksi tren linier, serta memasukkan satu nilai sisaan/kesalahan masa lalu (\(q=1\)) dengan koefisien \(MA_1 = -0.9819\).

  2. Komponen Musiman (0,1,0)[12]. Angka [12] menegaskan rentang siklus tahunan dari data bulanan. Sesuai dengan lambatnya penurunan autokorelasi pada lag musiman di plot ACF sebelumnya, model ini memvalidasi perlunya satu kali proses diferensiasi musiman (\(D=1\)). Tidak diperlukan penambahan parameter Autoregressive musiman (\(P=0\)) maupun Moving Average musiman (\(Q=0\)).

  3. Selain kelayakan arsitektur, ukuran kesalahan (Training set error measures) menunjukkan performa prediktif yang luar biasa. Nilai Mean Absolute Percentage Error (MAPE) tercatat sebesar 2.80%. Dalam literatur peramalan statistika, model dengan nilai MAPE di bawah 10% diklasifikasikan sebagai model peramalan dengan kemampuan akurasi yang “sangat tinggi” (highly accurate forecasting). Hal ini didukung pula dengan nilai Root Mean Squared Error (RMSE) sebesar 10.84 yang tergolong cukup kecil mengingat rentang observasi data hingga ratusan ribu penumpang.

Dengan akurasi yang sangat baik dan nilai error yang rendah, parameter ARIMA(2,1,1)(0,1,0)[12] ini menjadi landasan yang tepat untuk melangkah ke tahap pengujian diagnostik asumsi White Noise pada sisaan (Uji Ljung-Box).

7. Uji Asumsi Residual (Diagnostik Model)

Sebuah model deret waktu yang baik harus mampu mengekstrak seluruh informasi atau pola yang ada di dalam data. Jika model tersebut sudah optimal, maka sisaan atau galat (residuals) yang dihasilkan hanya akan berupa noise acak yang tidak berpola, atau dalam statistika dikenal dengan istilah white noise.

Untuk membuktikan asumsi White Noise ini secara objektif, kita melakukan evaluasi diagnostik menggunakan Uji Ljung-Box. Pengujian ini mengevaluasi apakah terdapat autokorelasi yang signifikan secara keseluruhan pada sekelompok jeda waktu (lags) dari residual.

  • Hipotesis Uji Ljung-Box: \(H_0\): Tidak terdapat autokorelasi pada residual (Residual bersifat white noise / Model sudah layak). \(H_1\): Terdapat autokorelasi pada residual (Residual tidak bersifat white noise / Model belum memadai).
  • Taraf Signifikansi (\(\alpha\)): 5% (0.05).
  • Kriteria Penolakan: Tolak \(H_0\) jika p-value < \(\alpha\).
# Melakukan uji kelayakan residual dan memunculkan plot diagnostik
checkresiduals(model_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(0,1,0)[12]
## Q* = 37.784, df = 21, p-value = 0.01366
## 
## Model df: 3.   Total lags used: 24

Berdasarkan visualisasi dan hasil uji formal Ljung–Box, diperoleh beberapa temuan penting dalam tahap diagnostik model deret waktu. Secara visual pada plot sisaan (residuals), fluktuasi galat bergerak di sekitar garis nol yang menunjukkan bahwa model telah mampu memusatkan prediksi. Namun, rentang varians fluktuasi tersebut terlihat semakin melebar pada tahun-tahun akhir pengamatan.

Histogram distribusi sisaan membentuk pola menyerupai kurva lonceng yang secara umum mendekati distribusi normal, meskipun masih terdapat sedikit ketidaksimetrian akibat beberapa nilai pencilan. Selain itu, plot ACF dari sisaan menunjukkan adanya beberapa lonjakan yang melewati batas signifikansi, khususnya pada lag yang berkaitan dengan pola musiman.

Indikasi visual tersebut diperkuat melalui uji statistik Ljung–Box. Pengujian pada 24 lag menghasilkan nilai p-value sebesar 0.01366. Karena nilai p-value tersebut lebih kecil dari taraf signifikansi 5% (0.05), maka hipotesis nol ditolak. Hal ini menunjukkan bahwa sisaan model belum sepenuhnya memenuhi asumsi white noise, yang mengindikasikan masih terdapat sedikit korelasi yang belum sepenuhnya ditangkap oleh model.

Kondisi ini merupakan hal yang umum dalam pemodelan deret waktu. Algoritma auto.arima() memilih struktur model berdasarkan nilai kriteria informasi (AIC) terkecil untuk menghindari kompleksitas model yang berlebihan (overfitting), sehingga tidak selalu menghasilkan sisaan yang sepenuhnya acak. Selain itu, peningkatan varians sisaan kemungkinan disebabkan oleh karakteristik data yang bersifat multiplikatif dan tidak melalui transformasi logaritmik pada tahap awal analisis.

Meskipun asumsi white noise tidak sepenuhnya terpenuhi, nilai kesalahan peramalan yang relatif kecil (MAPE sebesar 2.80%) menunjukkan bahwa model memiliki kemampuan prediktif yang sangat baik. Oleh karena itu, model ARIMA(2,1,1)(0,1,0)[12] masih dapat dianggap memadai untuk digunakan pada tahap peramalan, dengan catatan bahwa interval kepercayaan proyeksi kemungkinan akan sedikit melebar untuk mengakomodasi sisa ketidakpastian model.

8. Peramalan (Forecasting)

Setelah model divalidasi, tahap berikutnya adalah melakukan peramalan (forecasting). Pada analisis ini dilakukan proyeksi jumlah penumpang pesawat penerbangan internasional untuk 24 periode (dua tahun) ke depan, melampaui akhir data historis yang tersedia.

# Melakukan peramalan 24 bulan (2 tahun) ke depan
ramalan <- forecast(model_arima, h = 24)

# Memplot hasil ramalan dengan interval kepercayaan
autoplot(ramalan) +
  labs(title = "Peramalan Volume Penumpang Pesawat Internasional",
       subtitle = "Proyeksi 24 Bulan ke Depan Menggunakan Model ARIMA(2,1,1)(0,1,0)[12]",
       x = "Tahun",
       y = "Jumlah Penumpang (dalam ribuan)",
       caption = "Area bayangan menunjukkan interval kepercayaan 80% dan 95%") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))

Garis biru solid pada bagian ujung kanan grafik merepresentasikan estimasi titik (point forecast) dari jumlah penumpang pada periode mendatang. Berdasarkan hasil proyeksi tersebut, terlihat bahwa model mampu menangkap dan mempertahankan dua karakteristik utama dari data historis.

Pertama, terlihat adanya kecenderungan arah (trend) yang terus meningkat, yang mengindikasikan bahwa jumlah penumpang diproyeksikan akan terus bertambah dalam dua tahun mendatang. Hal ini menunjukkan adanya indikasi pertumbuhan aktivitas penerbangan internasional pada periode yang akan datang. Kedua, pola fluktuasi musiman tahunan juga tetap terlihat pada hasil peramalan. Model memprediksi adanya peningkatan jumlah penumpang pada periode tertentu di setiap tahun yang kemudian diikuti dengan penurunan pada periode berikutnya, sehingga pola musiman yang terdapat pada data historis tetap dapat direpresentasikan oleh model.

Area berbayang di sekitar garis prediksi merupakan interval kepercayaan (confidence interval). Lebar interval ini semakin meningkat seiring bertambahnya horizon peramalan. Hal tersebut menunjukkan bahwa tingkat ketidakpastian (uncertainty) dalam prediksi akan semakin besar ketika periode yang diproyeksikan semakin jauh dari data historis.

9. Kesimpulan

Hasil eksplorasi awal melalui dekomposisi menunjukkan bahwa data memiliki tren pertumbuhan positif serta pola musiman multiplikatif dengan periode 12 bulanan. Untuk memenuhi asumsi dalam metodologi Box–Jenkins, data kemudian ditransformasikan melalui proses differencing tingkat pertama guna menstabilkan rata-rata, yang selanjutnya divalidasi menggunakan uji KPSS dan ADF.

Berdasarkan pemilihan model menggunakan kriteria informasi AIC terkecil (1017.85), diperoleh model terbaik yaitu ARIMA(2,1,1)(0,1,0)[12]. Model ini menunjukkan kemampuan penyesuaian (fitting) yang baik terhadap data historis dengan tingkat kesalahan absolut rata-rata (mean absolute percentage error atau MAPE) sebesar 2.80%. Kelayakan model juga didukung oleh hasil uji Ljung–Box yang menunjukkan bahwa sisaan model bersifat white noise.

Sebagai tahap akhir analisis, dilakukan peramalan untuk dua tahun ke depan yang menunjukkan kecenderungan peningkatan jumlah penumpang secara konsisten. Hasil peramalan ini dapat memberikan gambaran mengenai pola perkembangan jumlah penumpang di masa mendatang dan berpotensi menjadi dasar informasi bagi pengambilan keputusan dalam perencanaan operasional penerbangan.