PEMODELAN SARIMA UNTUK PERAMALAN JUMLAH ANGKUTAN BARANG KERETA API DI PULAU JAWA PADA MASA PANDEMI COVID-19

1 PENDAHULUAN

1.1 Latar Belakang

      Tersedianya jaringan prasarana kereta api mendukung sistem distribusi logistik nasional serta menjadi angkutan penumpang massal perkotaan dan antarkota (Prasidi & Rifni, 2020). Besarnya potensi pasar pada angkutan kereta api terkait keunggulannya dibandingkan moda transportasi lain menurut Utomo (2009) adalah jangkauan pelayanan dengan kapasitas angkut besar, ketepatan waktu karena menggunakan jalur tersendiri, penggunaan energi berdampak efisiensi harga, polusi kecil, kecepatan dapat divariasikan dan aksesbilitas baik. Pada masa pandemi saat ini, tidak dapat dipungkiri bahwa segala jenis sektor kehidupan terkena dampaknya, Muhyiddin (2020) menyebutkan bahwa COVID-19 membuat perekonomian Indonesia melemah. Pemerintah telah melakukan berbagai cara untuk dapat memutus rantai penyebaran virus penyebab COVID-19, mulai dari penerapan protokol kesehatan hingga pembatasan mobilitas pada sektor transportasi. Dengan demikian, moda transportasi kereta api juga terkena dampak dari pandemi karena berkuranganya jumlah penumpang, sehingga upaya pemulihan kereta api harus dilakukan dengan menyesuaikan keadaan pada masa pandemi.
      Sejarah mencatat bahwa jasa pengangkutan barang sebagai tahap awal perkembangan perkeretaapian di tanah air. Oleh karena itu, Walaupun pandemi ini berdampak pada pengurangan jumlah penumpang, pelayanan terhadap berbagai komoditas angkutan barang tetap dijalankan oleh PT Kereta Api Indonesia (Persero). Data jumlah barang melalui transportasi kereta api di Jawa setiap bulan selalu dilakukan pembaharuan oleh Badan Pusat Statistik (BPS). Data tersebut dapat dikatakan sebagai data runtun waktu. Data runtun waktu adalah data dari satu atau kelompok unit yang diobservasi dalam beberapa periode berurutan. Supriadi (2020) juga menyebutkan bahwa data runtun waktu adalah jenis data yang dikumpulkan menurut urutan waktu dalam suatu rentang waktu tertentu. Dari data runtun waktu tersebut dapat dilakukan analisis runtun waktu yang dapat digunakan untuk meramalkan keputusan di masa yang akan datang (Lailani, 2019). Menurut Prasetya & Lukiastuti (2009), esensi peramalan adalah perkiraan peristiwa-peristiwa di waktu yang akan datang atas dasar pola-pola di waktu yang lalu, dan penggunaan kebijakan terhadap proyeksi-proyeksi dengan pola-pola di waktu yang lalu.
      Salah satu metode yang digunakan untuk peramalan data runtun waktu adalah Autoregressive Integrated Moving Average (ARIMA). Syarifuddin & Pratomo (2013) menyebutkan bahwa metode ARIMA adalah metode peramalan yang tidak menggunakan teori atau pengaruh antar variabel seperti model regresi, sehingga metode ini tidak memerlukan penjelasan terkait variabel bebas dan variabel terikat. Selain itu, terdapat pula metode Seasonal Autoregressive Integrated Moving Average (SARIMA) yang merupakan pengembangan dari model ARIMA pada data runtun waktu yang memiliki pola musiman atau memiliki jangka waktu panjang (Durrah et al, 2018). Penelitian terdahulu dengan topik yang serupa telah dilakukan oleh Utomo & Fanani (2020) dengan metode SARIMA dan menggunakan data jumlah penumpang kereta api di Indonesia. Didapatkan model terbaik untuk peramalan adalah ARIMA(1,1,2)(0,1,1)₁₂ dengan nilai MAPE sebesar 6.26%, nilai ini memiliki range dibawah 10% yang artinya model peramalan sangat baik. Sedangkan penelitian ini memfokuskan pada variabel yang berbeda yaitu jumlah barang angkutan melalui kereta api di Jawa. Dengan demikian penelitian pada topik ini akan memberikan peramalan jumlah barang angkutan agar PT Kereta Api Indonesia (Persero) dapat mengoptimalkan layanan angkutan barang yang secara tidak langsung dapat membantu pemulihan perekonomian di Indonesia pada masa pandemi. Mengingat kereta api dikatakan sebagai instrumen vital bagi negara dalam meraih kemajuan perekonomian.

1.2 Analisis Deret Waktu

      Analisis deret waktu merupakan metode peramalan di masa mendatang atau pengamatan suatu variabel dari waktu ke waktu dengan selang waktu yang tetap. Kemudian setiap pengamatan tersebut akan dinyatakan sebagai variabel Y_t dengan indeks waktu (t) tertentu sebagai urutan waktu pengamatan. Sehingga analisis deret waktu pada dasarnya digunakan untuk menganalisis data dengan mempertimbangkan waktu. Data deret waktu sendiri biasanya disajikan dalam kurun waktu harian, mingguan, bulanan, maupun tahunan.
      Berdasarkan analisis deret waktu, data akan digambarkan dengan suatu grafik untuk mengetahui perkembangan atau pola yang sedang diamati. Tahapan penting dalam memilih metode deret waktu yang tepat adalah dengan mempertimbangkan jenis pola data (Nasir, 2015). Pola data deret waktu dapat dibedakan menjadi empat sebagai berikut (Firdaus, 2020).

1. Pola trend, yaitu menunjukkan peningkatan atau penurunan dalam jangka panjang selama periode waktu yang diamati.
2. Pola musiman, yaitu adanya fluktuasi data yang berulang setiap beberapa ari, minggu, atau bulan karena afaktor cuaca, hari raya, dan lain-lain.
3. Pola siklis, yaitu terjadinya pola musiman dalam jangka waktu lebih panjang dan berulang biasanya setiap lima hingga sepuluh tahun.
4. Pola acak, ketika terjadinya fluktuasi dalam data yang disebabkan oleh variasi selain ketiganya di atas, yaitu faktor-faktor yang tidak diantisipasi.

1.2.1 Autocorrelation Function (ACF)

      Autocorrelation Function (ACF) merupakan salah satu konsep yang berkaitan dengan analisis deret waktu. Koefisien korelasi merupakan statistik kunci dalam analisis deret waktu, yaitu menyatakan ukuran korelasi deret waktu dengan dirinya sendiri dengan selisih waktu (lag) 0, 1, 2 periode atau lebih (Nasir, 2015). Berdasarkan Wei (2006), nilai autokorelasi antara Y_t dan Y_t+k dapat diformulasikan sebagai berikut.

\[ \rho_k =\frac {cov(Y_t,Y_{t+k})} {\sqrt{var(Y_t)} \sqrt{var(Y_{t+k})}} =\frac {\gamma_k} {\gamma_0} \tag{1} \] Dimana \(Var(Y_t )=Var(Y_{t+k} )=\gamma_0\). Sebagai fungsi dari k, \(\gamma_k\) merupakan fungsi autokovarian dan \(\rho_k\) disebut dengan fungsi autokorelasi atau ACF. Berikut sifat-sifat penting bagi \(\gamma_k\) dan \(\rho_k\).

1. \(\gamma_0=var(Y_t)\) dan \(\rho_0=1\)
2. \(|\gamma_k| \le \gamma_0\) dan \(\rho_k \le 1\)
3. \(\gamma_k = \gamma_{-k}\) dan \(\rho_k = \rho_{-k}\)
4. \(\sum^{n}_{i=1}{\sum^{n}_{j=1}}{\alpha_i \alpha_j \gamma_{|t_i-t_j|}} \ge 0\) dan \(\sum^{n}_{i=1}{\sum^{n}_{j=1}}{\alpha_i \alpha_j \rho_{|t_i-t_j|}} \ge 0\)

1.2.2 Partial Autocorrelation Function (PACF)

      Partial Autocorrelation Function (PACF) juga merupakan ukuran korelasi pada analisis deret waktu. Fungsi autokorelasi parsial adalah suatu fungsi yang menunjukkan besarnya korelasi antara pengamatan ke t yaitu Y_t denga pengamatan waktu-waktu sebelumnya yaitu 𝑍_𝑡−1, 𝑍_𝑡−2, … , 𝑍_𝑡−𝑘. Berdasarkan Wei (2006), nilai parsial autokorelasi pada lag ke-k adalah sebagai berikut.

\[ \phi_k= \frac {\left| \begin{array}{cc} 1&\rho_1&\rho_2&\ldots&\rho_{k-2}&\rho_1\\ \rho_1&1&\rho_1&\ldots&\rho_{k-3}&\rho_2\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ \rho_{k-1}&\rho_{k-2}&\rho_{k-3}&...&\rho_1&\rho_k \end{array} \right| } {\left| \begin{array}{cc} 1&\rho_1&\rho_2&...&\rho_{k-2}&\rho_{k-1}\\ \rho_1&1&\rho_1&...&\rho_{k-3}&\rho_{k-2}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ \rho_{k-1}&\rho_{k-2}&\rho_{k-3}&...&\rho_1&1 \end{array} \right|} \tag{2} \]

1.3 Model Seasonal Autoregressive Integrated Moving Average (SARIMA)

      Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan metode peramalan time series untuk model data fluktuatif dengan pola data musiman (Fahrudin et all, 2020). Model ini juga merupakan pengembangan dari model Autoregressive Integrated Moving Average (ARIMA). Sehingga model Seasonal ARIMA atau SARIMA dapat dinotasikan dengan ARIMA(p, d, q)(P, D, Q)^s, dimana,

• (p, d, q) : ordo AR, differencing, dan MA non-seasonal
  
• (P, D, Q) : ordo AR, differencing, dan MA seasonal
  
• s : jumlah periode setiap musim

Model SARIMA jika dinyatakan dalam sebuah persamaan akan berbentuk sebagai berikut. \[ \phi_p(B^s)\emptyset_p(B)(1-B)^d(1-B^s)^DZ_t=\theta_q(B)\Theta_Q(B^S)\alpha_t \tag{3} \] Keterangan:

• \(\emptyset_p(B):\) AR non-seasonal
• \(\phi_p(B^s):\) AR seasonal
• \((1-B)^d:\) differencing non-seasonal
• \((1-B)^D:\) differencing seasonal
• \(\theta_q(B):\) MA non-seasonal
• \(\Theta_Q(B^S):\) MA seasonal

1.3.1 Identifikasi Model

      Pengidentifikasian karakteristik data penting untuk dilakukan sebelum menganalisis deret waktu. Salah satunya dengan melihat pola data yang terbentuk. Jika data deret waktu memiliki pola yang berulang setiap s periode waktu, maka data tersebut memiliki unsur musiman yang kemudian dimodelkan dengan metode SARIMA. Berdasarkan plot data tersebut sekaligus dapat mengetahui kestasioneran data. Kestasioneran suatu data terdiri dari dua hal yaitu stasioner dalam ragam dan rata-rata.
      Stasioneritas data terjadi ketika pola data konstan disekitar rata-rata dan variansi atau tidak terjadi kenaikan dan penurunan selama periode waktu tertentu. Stasioneritas ragam dapat diperiksa dengan melihat nilai lamda (λ) pada plot uji Box-cox. Jika nilai λ menunjukkan sama dengan satu, maka data sudah dapat dikatakan stasioner terhadap ragam. Jika kondisi tersebut tidak terpenuhi, maka harus dilakukan transformasi Box-cox, dengan bentuk transformasi sebagai berikut (Nasir, 2015).

\[ T(Y_t)=\frac {Y_t^\lambda-1} {\lambda} \tag{4} \] Dimana 𝜆 dinamakan parameter transformasi. Di bawah ini adalah beberapa nilai 𝜆 beserta transformasinya.

Tabel 1. Nilai 𝜆 beserta Transformasinya

Nilai Lamda (λ)	Transformasi
\(-1\)	\(\frac{1} {Y_t}\)
\(-0.5\)	\(\frac{1} {\sqrt{Y_t}}\)
\(0\)	\(Ln(Y_t)\)
\(0.5\)	\(\sqrt{Y_t}\)
\(1\)	\(Y_t\)

      Stasioneritas rata-rata dapat dideteksi dengan melihat plot deret waktu, plot fungsi autokorelasi atau ACF, dan menggunakan uji Augmented Dickey-Fuller (ADF). Hipotesis yang akan diuji jika menggunakan uji ADF dan kriteria uji adalah sebagai berikut.

H₀: Data tidak stasioner terhadap rata-rata
H₁: Data stasioner terhadap rata-rata

Dengan kriteria pengujian,

1. Jika p-value < nilai alpha (0.05), maka H₀ ditolak. H₀ ditolak berarti data stasioner terhadap rata-rata.
2. Jika p-value > nilai alpha (0.05), maka H₀ diterima. H₀ diterinma berarti data tidak stasioner terhadap rata-rata.

Berdasarkan Munawaroh (2010), apabila data tidak stasioner pada rata-rata, maka dapat dikonversikan menjadi deret stasioner melalui differensiasi. Bentuk differensiasi pertama (d = 1) adalah sebagai berikut. \[ \nabla Y_t=Y_t-Y_{t-1} \tag{5} \] Sedangkan bentuk differensiasi kedua (d = 2) adalah sebagai berikut. \[ \nabla^2 Y_t=\nabla Y_t-\nabla Y_{t-1} \tag{6} \]

      Selanjutnya pada pengidentifikasian model atau penentuan orde SARIMA dapat dibantu dengan mengamati plot fungsi autokorelasi (ACF) dan fungsi parsial autokorelasi (PACF). Ordo p ditentukan berdasarkan lag non-seasonal yang keluar atau signifikan pada plot PACF dan orde P ditentukan berdasarkan lag seasonal yang keluar atau signifikan pada plot PACF. Sedangkan ordo q ditentukan berdasarkan lag non-seasonal yang keluar atau signifikan pada plot ACF dan orde Q ditentukan berdasarkan lag seasonal yang keluar atau signifikan pada plot ACF.

1.3.2 Estimasi dan Pengujian Parameter

      Parameter dalam model SARIMA perlu dilakukan estimasi. Metode yang digunakan untuk estimasi parameter dalam model SARIMA adalah metode Maximum Likelihood. Kemudian masing-masing estimasi parameter yang didapatkan akan diuji signifikansinya dalam model. Signifikansi parameter tersebut bertujuan untuk mengetahui kelayakan model sehingga dapat digunakan. Secara umum, misalkan \(\phi\) adalah parameter pada model SARIMA dan \(\hat\phi\) adalah nilai estimasi parameter tersebut, serta \(s.e(\hat\phi)\) adalah standard error dari nilai estimasi tersebut.
      Hipotesis yang akan diuji pada signifikansi parameter adalah sebagai berikut.

H₀: \(\phi=0\) (Parameter pada model signifikan)
H₁: \(\phi\ne0\) (Parameter pada model tidak signifikan)

Dengan kriteria pengujian,

1. Jika p-value < nilai alpha (0.05), maka H₀ ditolak. H₀ ditolak berarti parameter pada model signifikan.
2. Jika p-value > nilai alpha (0.05), maka H₀ diterima. H₀ diterima berarti parameter pada model tidak signifikan.

Statistik uji yang digunakan adalah sebagai berikut. \[ t_{hitung}=\frac {\hat\phi}{s.e(\hat\phi)} \tag{7} \]

1.3.3 Diagnostik Model

      Diagnostik model dilakukan pada setiap model tentatif yang didapatkan. Dimana residual model tersebut yang kemudian akan diperiksa apakah telah memenuhi asumsi atau tidak. Diagnostik model terdiri dari uji asumsi kenormalan dan white noise.
      Pengujian kenormalan digunakan untuk mengetahui apakah data telah memenuhi asumsi kenormalan atau belum. Salah satu uji yang dapat digunakan untuk mengambil keputusan tersebut adalah dengan uji Jarque Bera (JB). Berikut hipotesis beserta kriteria uji JB.

H₀: residual berdistribusi normal
H₁: residual tidak berdistribusi normal

Dengan kriteria pengujian,

1. Jika p-value < nilai alpha (0.05), maka H₀ ditolak. H₀ ditolak berarti residual tidak berdistribusi normal.
2. Jika p-value > nilai alpha (0.05), maka H₀ diterima. H₀ diterinma berarti residual berdistribusi normal.

Statistik uji pada Jarque Bera (JB) adalah (Gujarati, 2006), \[ JB=\frac {n}{6}(S^2+\frac {(K-3)^2}{4}) \tag{8} \] Dimana n merupakan ukuran sampel, S menyatakan kemencengan, dan K menyatakan peruncingan. Statistik JB di atas mengikuti distribusi chi-square dengan d.k 2 secara asimtotis (dalam sampel besar) sebagai berikut. \[ JB_{asy} \sim \chi^2_{(2)} \tag{9} \]

      Pengujian asumsi white noise dilakukan dengan menggunakan uji Ljung-Box. Berikut hipotesis beserta kriteria uji Ljung-Box berdasarkan (Sitepu & Sinaga, 2018).

H₀: residual white noise
H₁: residual tidak white noise

Dengan kriteria pengujian,

1. Jika \(Q>\chi^2\) (0.05), maka H₀ ditolak. H₀ ditolak berarti residual tidak white noise.
2. Jika \(Q<\chi^2\), maka H₀ diterima. H₀ diterinma berarti residual white noise.

Statistik uji Ljung-Box adalah, \[ Q=n(n+2) \sum^K_{k=1} \frac{\hat\rho^2_k}{(n-k)}, n>K \tag{10} \] Keterangan:

• \(\hat\rho^2_k:\) autokorelasi residual pada lag ke-k
• \(n:\) banyaknya pengamatan
• \(K:\) jumlah dari lag waktu yang dimasukkan dalam uji

1.3.4 Pemilihan Model Terbaik

      Dalam menentukan model terbaik dari beberapa model tentatif yang telah melalui proses pemeriksaan asumsi, perlu dilakukan perhitungan model residual yang sesuai berdasarkan kesalahan peramalan. Beberapa kriteria yang digunakan dalam pemilihan model terbaik yaitu (Putri & Aghsilni, 2019),

1. Akaike’s Information Criterion (AIC)
   Akaike’s Information Criterion (AIC) pertama kali diperkenalkan oleh Akaike untuk mengidentifikasi model dari suatu kumpulan data. Model dikatakan baik jika memiliki nilai AIC terkecil. Berikut persamaan AIC dalam pemilihan model terbaik. \[ AIC=\log\hat\sigma^2+\frac {2k}{n} \tag{11} \] Keterangan:

• \(\log\hat\sigma^2:\) ukuran likelihood
• \(k:\) jumlah parameter
• \(n:\) banyaknya pengamatan

2. Bayesin Information Criterion (BIC)
   Bayesin Information Criterion (BIC) merupakan suatu tipe metode pemilihan model dengan pendekatan Penalized Maximum Likelihood. Pendekatan tersebut pertama kali diperkenalkan oleh Schwartz. Model dikatakan baik jika memiliki nilai BIC terkecil. Berikut persamaan BIC dalam pemilihan model terbaik. \[ BIC=\log\hat\sigma^2+\frac {k \log n}{n} \tag{12} \] Keterangan:

• \(\log\hat\sigma^2:\) ukuran likelihood
• \(k:\) jumlah parameter
• \(n:\) banyaknya pengamatan

1.4 Kereta Api

      PT Kereta Api Indonesia (Persero) merupakan Badan Usaha Milik Negara (BUMN) kategori perusahaan yang mendapat pelimpahan tugas menyelenggarakan pelayanan publik dalam bidang angkutan kereta api (Samawati, 2019). Kereta api merupakan transportasi darat yang bergerak di atas rel. Berdasarkan PP No.72/2009 tentang Lalu Lintas dan Angkutan Kereta Api pasal 2 ayat 2 menyebutkan jaringan pelayanan perkeretaapian terdiri dari dua jaringan pelayanan yaitu jaringan pelayanan perkeretaapian antar kota dan perkotaan (Biomantara et al, 2018). Sehingga layanan kereta api secara umum cukup menjadi andalan masyarakat sebagai alternatif transportasi umum untuk perjalanan luar kota ataupun sekedar dalam kota. Hal tersebut dikarenakan transportasi kereta api lebih efektif dan efisien jika dibandingkan dengan moda transportasi darat lainnya.
      Selain sebagai transportasi untuk penumpang, kereta api juga menyediakan layanan angkutan barang. PT Kereta Api Indonesia (Persero) saat ini sedang memaksimalkan layanan angkutan barang di tengah pandemi COVID-19. Mengingat angkutan barang juga merupakan salah satu bagian utama yang dijalankan oleh PT Kereta Api Indonesia (Persero) sekaligus sebagai awal mula perkembangan perkeretaapian di Indonesia. Pada kereta api angkutan barang yang sebagian besar terkait dengan para pengusaha angkutan dan pengiriman barang corporate yaitu konteiner atau peti kemas, semen, pupuk, dan lain-lain (Biomantara et al, 2018). Sihombing (2011) menyebutkan bahwa perkeretaapian merupakan salah satu instrumen yang berfungsi sebagai tulang punggung sistem logistik dan distribusi nasional di dalam perekonomian Indonesia ke depan.

1.5 Data

      Dalam penelitian ini menggunakan data sekunder yang diperoleh dari situs resmi Badan Pusat Statistik (bps.go.id). Data yang digunakan adalah data jumlah barang angkutan melalui transportasi kereta api di Pulau Jawa selama 70 bulan dari bulan Januari 2016 hingga Oktober 2021. Berikut ringkasan data tersebut yang digunakan dalam penelitian.

2 SOURCE CODE

2.1 Library yang Dibutuhkan

> library(readxl)
> library(FitAR)
> library(tseries)
> library(forecast)
> library(lmtest)
> library(normtest)

Source code di atas digunakan untuk mengaktifkan package agar fungsi yang terdapat pada package dapat digunakan. Package “readxl” digunakan untuk membaca file excel yang diimport ke R. Fungsi dari package “FitAR” biasa digunakan untuk plot uji Box-Cox dalam pemeriksaan stasioneritas ragam. Package “tseries” berisi fungsi-fungsi dalam analisis time series seperti uji ADF dengan fungsi adf.test, garch, dan lain sebagainya. Sedangkan package “lmtest” biasa digunakan untuk menguji model regresi linier. Namun, dalam hal ini memanfaatkan fungsi coeftest dari “lmtest” untuk membantu uji signifikansi parameter model deret waktu. Package terakhir adalah “normtest” yang berguna untu uji normalitas.

2.2 Import Data

> data_barang <- read_excel("E:/#Univ Brawijaya/SEMESTER 5/MPPI/UAS/Data Jumlah Barang.xlsx", 
+                           col_types = c("text", "numeric"))
> head(data_barang)
# A tibble: 6 x 2
  Waktu         `Jumlah Barang`
  <chr>                   <dbl>
1 Januari 2016              927
2 Februari 2016             734
3 Maret 2016                785
4 April 2016                967
5 Mei 2016                  873
6 Juni 2016                 945

Data yang digunakan dalam analisis berasal dari file excel dan diimport dengan memanfaatkan fungsi read_excel. Kemudian data disimpan dalam obyek bernama “data_barang”. Fungsi head digunakan untuk menampilkan enam baris pengamatan pertama dari data set “data_barang”.

2.3 Plot Data

> plot(data_barang$`Jumlah Barang`, ylab='Jumlah Angkutan Barang KAI', xlab='Waktu', main='Gambar 1. Plot Data Jumlah Angkutan Barang Kereta Api di Pulau Jawa', type='o')

Fungsi plot digunakan untuk membuat dan menampilkan plot data “Jumlah Barang” dari data set “data_barang” dengan judul “Gambar 1. Plot Data Jumlah Angkutan Barang Kereta Api di Pulau Jawa”. Simbol titik yang akan ditampilkan pada plot time series adalah lingkaran.

2.4 Stasioneritas Data

2.4.1 Stasioneritas Ragam

> BoxCox.ts(data_barang$`Jumlah Barang`)

Analisis Box-Cox untuk deret waktu variabel “Jumlah Barang” dari data set “data_barang” menggunakan fungsi BoxCox.ts. Fungsi ini digunakan untuk memeriksa stasioneritas ragam dengan melihat nilai lamda (λ) yang dihasilkan.

2.4.2 Stasioneritas Rata-rata

> #PLOT ACF PACF AWAL
> acf(data_barang[,2],lag.max = 50)

> pacf(data_barang[,2],lag.max=50)

Source code di atas digunakan untuk menampilkan plot ACF dan PACF dari data pada kolom kedua data set “data_barang” yaitu “Jumlah Barang” dengan maksimum lag 50. Fungsi acf untuk menghitung estimasi fungsi autokorelasi, sedangkan pacf untuk autokorelasi parsial.

> #STASIONERITAS RATA-RATA NON MUSIMAN
> adf.test(data_barang$`Jumlah Barang`)

    Augmented Dickey-Fuller Test

data:  data_barang$`Jumlah Barang`
Dickey-Fuller = -1.9542, Lag order = 4, p-value = 0.5938
alternative hypothesis: stationary

Fungsi adf.test digunakan untuk menghitung uji Augmented Dickey-Fuller. Perintah tersebut berguna untuk menguji stasioneritas data “Jumlah Barang” terhadap rata-rata.

> #DIFFERENCING NON MUSIMAN
> diffnon=diff(data_barang$`Jumlah Barang`)
> adf.test(diffnon)

    Augmented Dickey-Fuller Test

data:  diffnon
Dickey-Fuller = -3.642, Lag order = 4, p-value = 0.03634
alternative hypothesis: stationary
> acf(diffnon, lag.max = 50)

Fungsi diff berguna untuk mendifferensiasi data deret waktu yang tidak stasioner terhadap rata-rata yaitu “Jumlah Barang” yang berada pada data set “data_barang”. Kemudian data hasil differensiasi disimpan ke dalam variabel bernama “diffnon”. Data “diffnon” kembali diuji stasioneritas rata-ratanya dengan fungsi adf.test dan membentuk plot fungsi autokorelasi dengan maksimum lag 50, yang memanfaatkan fungsi acf.

> #DIFFERENCING MUSIMAN
> diffmus1=diff(diffnon,lag=12)
> acf(diffmus1,lag.max = 50)

> diffmus2=diff(diff(diffnon,lag=12),lag=12)

Membentuk variabel baru bernama “diffmus1” dari hasil differensiasi data “diffnon” dengan bantuan fungsi diff. lag = 12 pada fungsi diff menandakan bahwa differensiasi yang dilakukan terhadap musiman yang berulang setiap 12 periode. Dari data “diffmus1” kembali dibentuk plot fungsi autokorelasi dengan maksimum lag 50, yang memanfaatkan fungsi acf. Kemudian membentuk variabel bernama “diffmus2” dari hasil differensiasi musiman kedua data “diffnon”.

2.5 Identifikasi Model

> acf(diffmus2, lag.max = 50)

> pacf(diffmus2, lag.max = 50)

Pembentukan plot fungsi autokorelasi dengan fungsi acf dan plot autokorelasi parsial menggunakan fungsi pacf dengan lag maksimum yang ditampilkan hingga 50. Objek deret waktu yang digunakan untuk pembentukan plot ACF dan PACF adalah “diffmus2”. Kedua plot tersebut berguna untuk mengidentifikasi model, dimana plot PACF membantu sebagai penentu ordo AR dan plot ACF sebagai penentu ordo MA.

2.6 Estimasi Parameter Model SARIMA

> #ESTIMASI PARAMETER SARIMA (1,1,1)(0,2,0)
> model1=Arima(data_barang$`Jumlah Barang`, order = c(1,1,1),
+                                  seasonal = list(order = c(0,2,0),period=12))
> summary(model1)
Series: data_barang$`Jumlah Barang` 
ARIMA(1,1,1)(0,2,0)[12] 

Coefficients:
          ar1      ma1
      -0.0728  -0.5347
s.e.   0.2575   0.2203

sigma^2 estimated as 29359:  log likelihood=-298.29
AIC=602.58   AICc=603.17   BIC=608

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE     MASE
Training set 2.525088 134.2932 87.96716 0.6970509 8.961646 0.773116
                     ACF1
Training set 0.0003740404
> coeftest(model1)

z test of coefficients:

     Estimate Std. Error z value Pr(>|z|)  
ar1 -0.072786   0.257507 -0.2827  0.77744  
ma1 -0.534668   0.220329 -2.4267  0.01524 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perintah Arima digunakan untuk analisis model ARIMA dari data time series “Jumlah Barang” yang berasal dari data set “data_barang”. Ordo non-seasonal model ARIMA yang digunakan adalah ordo AR 1, ordo differencing 1. dan ordo MA 1 (1,1,1), sedangkan ordo seasonal dalam 12 periode dengan urutan yang sama adalah (0,2,0). Hasil analisis disimpan dalam obyek bernama “model1”. Fungsi summary digunakan untuk menampilkan ringkasan hasil dari analisis deret waktu seasonal ARIMA atau SARIMA yang telah tersimpan pada obyek “model1”. Selain itu, terdapat fungsi coeftest yang berguna dalam menampilkan nilai z atau hasil uji wald dari estimasi parameter “model1” untuk uji signifikansi parameter model SARIMA yang terbentuk.

> #ESTIMASI PARAMETER SARIMA (1,1,0)(0,2,0)
> model2=Arima(data_barang$`Jumlah Barang`, order = c(1,1,0),
+              seasonal = list(order = c(0,2,0),period=12))
> summary(model2)
Series: data_barang$`Jumlah Barang` 
ARIMA(1,1,0)(0,2,0)[12] 

Coefficients:
          ar1
      -0.4858
s.e.   0.1299

sigma^2 estimated as 29986:  log likelihood=-299.21
AIC=602.43   AICc=602.71   BIC=606.04

Training set error measures:
                  ME   RMSE      MAE       MPE    MAPE      MASE        ACF1
Training set 2.15584 137.29 90.72997 0.5039281 9.26722 0.7973975 -0.05406871
> coeftest(model2)

z test of coefficients:

    Estimate Std. Error z value  Pr(>|z|)    
ar1 -0.48583    0.12994 -3.7388 0.0001849 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perintah Arima digunakan untuk analisis model ARIMA dari data time series “Jumlah Barang” yang berasal dari data set “data_barang”. Ordo non-seasonal model ARIMA yang digunakan adalah ordo AR 1, ordo differencing 1. dan ordo MA 0 (1,1,0), sedangkan ordo seasonal dalam 12 periode dengan urutan yang sama adalah (0,2,0). Hasil analisis disimpan dalam obyek bernama “model2”. Fungsi summary digunakan untuk menampilkan ringkasan hasil dari analisis deret waktu seasonal ARIMA atau SARIMA yang telah tersimpan pada obyek “model2”. Selain itu, terdapat fungsi coeftest yang berguna dalam menampilkan nilai z atau hasil uji wald dari estimasi parameter “model2” untuk uji signifikansi parameter model SARIMA yang terbentuk.

> #ESTIMASI PARAMETER SARIMA (0,1,1)(0,2,0)
> model3=Arima(data_barang$`Jumlah Barang`, order = c(0,1,1),
+              seasonal = list(order = c(0,2,0),period=12))
> summary(model3)
Series: data_barang$`Jumlah Barang` 
ARIMA(0,1,1)(0,2,0)[12] 

Coefficients:
          ma1
      -0.5816
s.e.   0.1152

sigma^2 estimated as 28748:  log likelihood=-298.33
AIC=600.67   AICc=600.95   BIC=604.28

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE
Training set 2.570215 134.4249 88.37225 0.7275877 8.999508 0.7766763
                    ACF1
Training set -0.03030203
> coeftest(model3)

z test of coefficients:

    Estimate Std. Error z value  Pr(>|z|)    
ma1 -0.58159    0.11521  -5.048 4.466e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perintah Arima digunakan untuk analisis model ARIMA dari data time series “Jumlah Barang” yang berasal dari data set “data_barang”. Ordo non-seasonal model ARIMA yang digunakan adalah ordo AR 0, ordo differencing 1. dan ordo MA 1 (0,1,1), sedangkan ordo seasonal dalam 12 periode dengan urutan yang sama adalah (0,2,0). Hasil analisis disimpan dalam obyek bernama “model3”. Fungsi summary digunakan untuk menampilkan ringkasan hasil dari analisis deret waktu seasonal ARIMA atau SARIMA yang telah tersimpan pada obyek “model3”. Selain itu, terdapat fungsi coeftest yang berguna dalam menampilkan nilai z atau hasil uji wald dari estimasi parameter “model3” untuk uji signifikansi parameter model SARIMA yang terbentuk.

2.7 Diagnostik Model

> #PERIKSA RESIDUAL MODEL SARIMA (1,1,0)(0,2,0)
> checkresiduals(model2)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(0,2,0)[12]
Q* = 10.369, df = 9, p-value = 0.3215

Model df: 1.   Total lags used: 10
> jb.norm.test(model2$residuals)

    Jarque-Bera test for normality

data:  model2$residuals
JB = 16.197, p-value = 0.004

Perintah untuk memeriksa residual dari model deret waktu apakah white noise atau tidak, menggunakan fungsi checkresiduals terhadap model seasonal ARIMA yang telah terbentuk yaitu “model2”. Fungsi tersebut salah satunya akan menghasilkan output uji Ljung-Box yang digunakan untuk uji asumsi white noise. Selanjutnya terdapat fungsi jb.norm.test yang berguna untuk uji normalitas residual dari “model2”.

> #PERIKSA RESIDUAL MODEL SARIMA (0,1,1)(0,2,0)
> checkresiduals(model3)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(0,2,0)[12]
Q* = 10.311, df = 9, p-value = 0.3259

Model df: 1.   Total lags used: 10
> jb.norm.test(model3$residuals)

    Jarque-Bera test for normality

data:  model3$residuals
JB = 41.147, p-value = 0.001

Perintah untuk memeriksa residual dari model deret waktu apakah white noise atau tidak, menggunakan fungsi checkresiduals terhadap model seasonal ARIMA yang telah terbentuk yaitu “model3”. Fungsi tersebut salah satunya akan menghasilkan output uji Ljung-Box yang digunakan untuk uji asumsi white noise. Selanjutnya terdapat fungsi jb.norm.test yang berguna untuk uji normalitas residual dari “model3”.

2.8 Model Terbaik

> AIC(model2)
[1] 602.4252
> AIC(model3)
[1] 600.6674

Fungsi AIC digunakan untuk menghitung Akaike’s Information Criterion (AIC) dari model seasonal ARIMA yang telah terbentuk yaitu “model2” dan “model3”. Nilai AIC yang diperoleh akan digunakan untuk menentukan model deret waktu terbaik berdasarkan nilai yang terkecil.

2.9 Peramalan

> fmodel3 = forecast(model3, h=3)
> fmodel3
   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
71       840.8195 623.5301 1058.109 508.5041 1173.135
72       937.8195 702.2769 1173.362 577.5882 1298.051
73      1067.8195 815.3400 1320.299 681.6854 1453.954
> #PLOT PERAMALAN
> plot(fmodel3)

Penggunaan fungsi forecast untuk peramalan model deret waktu yang telah terbentuk yaitu “model3”, dengan jumlah periode peramalan sebanyak tiga. Data hasil peramalan disimpan dalam obyek bernama “fmodel3”. Kemudian hasil peramalan “fmodel3” digambarkan melalui plot time series dengan fungsi plot

3 HASIL DAN PEMBAHASAN

3.1 Hasil Penelitian

3.1.1 Identifikasi Plot Time Series

      Pengidentifikasian plot time series berguna untuk mengetahui bentuk atau tren dari data penelitian yang digunakan. Berdasarkan Gambar 1 pada sub-bab 2.3, terlihat data jumlah angkutan barang kereta api di Pulau Jawa memiliki pola musiman dimana data akan naik dan turun dalam jangka waktu yang tetap dan berulang. Terlihat pula pola data yang tidak konstan disekitar rata-rata sehingga ada indikasi data belum stasioner.

3.1.2 Identifikasi Stasioneritas Data

      Berdasarkan plot data pada Gambar 1 yang mengindikasikan data belum stasioner, maka akan dilakukan identifikasi stasioneritas data dahulu sebelum melanjutkan analisis. Hal pertama yang perlu dilakukan adalah mengidentifikasi stasioneritas data terhadap ragam menggunakan plot uji Box-Cox pada Gambar 2 di sub-bab 2.4.1. Nilai lambda yang dihasilkan adalah satu, maka data sudah stasioner terhadap ragam. Oleh karena itu, data tidak perlu dilakukan transformasi Box-cox. Selanjutnya dilakukan identifikasi kestasioneran data terhadap rata-rata dengan terlebih dahulu melihat plot Autocorrelation Function (ACF) terhadap data aktual.
      Berdasarkan Plot ACF pada Gambar 3 di sub-bab 2.4.2, terlihat bahwa terdapat lebih dari tiga lag pertama yang signifikan sehingga mengindikasikan data belum stasioner terhadap rata-rata. Sedangkan secara value dapat dilakukan uji Augmented Dickey Fuller (ADF) untuk mengetahui secara pasti apakah data belum stasioner terhadap rata-rata. Hipotesis yang digunakan pada uji ADF adalah sebagai berikut.

H₀: data tidak stasioner terhadap rata-rata
H₁: data stasioner terhadap rata-rata

      Hasil uji ADF pada sub-bab 2.4.2, diperoleh p-value sebesar 0.5938 lebih besar dari alpha (0.05) maka terima H₀. Dengan demikian, dapat disimpulkan bahwa data tidak stasioner terhadap rata-rata sehingga perlu dilakukan penanganan yaitu dengan differencing.
      Selanjutnya dilakukan kestasioneran data dalam rata-rata dengan melakukan differencing pertama pada data aktual. Hasil uji ADF setelah differencing non musiman satu kali dengan hipotesis uji yang sama adalah p-value (0.03634) lebih kecil dari alpha (0.05) maka tolak H₀. Dengan demikian, dapat disimpulkan bahwa data non musiman telah stasioner terhadap rata-rata.
      Kemudian membentuk Plot ACF setelah data di differencing non musiman pada Gambar 4 di sub-bab 2.4.2. Berdasarkan Gambar tersebut, terlihat bahwa hanya lag pertama yang signifikan sehingga data non musiman yang telah diuji ADF sebelumnya benar telah stasioner terhadap rata-rata. Namun, masih terdapat tiga yang signifikan pada lag berkelipatan 12 sehingga efek musiman belum stasioner. Oleh karena itu, perlu dilakukan differencing musiman dimana lag order yang digunakan sebesar 12.
      Setelah dilakukan differencing musiman sebanyak dua kali, diilakukan pembentukan Plot ACF kembali pada Gambar 5 di sub-bab 2.5. Berdasarkan Gambar tersebut, terlihat bahwa sudah tidak ada tiga lag musiman atau lebih yang signifikan. Dengan demikian, data sudah stasioner terhadap rata-rata baik non musiman maupun musiman, sebab telah dilakukan differencing non musiman sebanyak satu kali dan differencing musiman sebanyak dua kali. Tahap selanjutnya adalah mengidentifikasi model tentatif.

3.1.3 Identifikasi Model

      Apabila data telah stasioner terhadap ragam dan rata-rata, maka akan diperoleh model tentatif. Identifikasi model tentatif atau penentuan orde SARIMA berdasarkan lag yang signifikan pada plot ACF dan PACF yang telah stasioner pada sub-bab 2.5. Berdasarkan Gambar 5 terlihat bahwa lag 1 pada plot ACF signifikan sehingga dapat ditentukan MA(1) untuk ordo non musiman yaitu q = 1. Sedangkan pada plot PACF yang terlihat pada Gambar 6, lag 1 juga signifikan sehingga dapat ditentukan AR(1) untuk ordo non musiman yaitu p = 1. Diketahui bahwa data stasioner setelah dilakukan differencing non musiman sebanyak satu kali dan differencing musiman sebanyak dua kali sehingga ordo differencing non musiman yang diperoleh adalah d = 1, sedangkan ordo differencing musiman adalah D = 1. Berdasarkan penentuan orde tersebut, diperoleh model tentatif sebanyak tiga yaitu SARIMA(1,1,1)(0,2,0)₁₂, SARIMA(1,1,0)(0,2,0)₁₂, dan SARIMA(0,1,1)(0,2,0)₁₂.

3.1.4 Estimasi dan Uji Signifikansi Parameter Model

      Hasil estimasi parameter model tentatif dan uji signifikansinya terdapat pada sub-bab 2.6 dengan hipotesis uji yang digunakan adalah sebagai berikut.

H₀: parameter model tidak signifikan
H₁: parameter model signifikan
Berdasarkan hasil estimasi parameter pada sub-bab 2.6, model SARIMA(1,1,1)(0,2,0)₁₂ memiliki dua parameter yaitu AR(1) dengan p-value (0.7774) lebih besar dari alpha (0.05) maka terima H₀ dan parameter MA (1) dengan p-value (0.0152) lebih kecil dari alpha (0.05) maka tolak H₀. Dengan demikian, parameter AR(1) tidak signifikan sedangkan parameter MA(1) signifikan sehingga dapat disimpulkan model SARIMA(1,1,1)(0,2,0)₁₂ tidak lolos uji signifikansi.

      Sedangkan hasil estimasi parameter model SARIMA(1,1,0)(0,2,0)₁₂ memiliki parameter AR(1) dengan p-value (0.00018) lebih kecil dari alpha (0.05) maka tolak H₀ sehingga parameter tersebut signifikan. Begitupun model SARIMA(0,1,1)(0,2,0)₁₂ juga memiliki satu parameter yaitu MA(1) dengan p-value (4.466 x 10^-7) lebih kecil dari alpha (0.05) maka tolak H₀ sehingga parameter tersebut signifikan. Dengan demikian, model SARIMA(1,1,0)(0,2,0)₁₂ dan SARIMA(0,1,1)(0,2,0)₁₂ lolos uji signifikansi untuk selanjutnya dilakukan uji diagnostik sisaan model.

3.1.5 Diagnostik Model

      Diagnostik sisaan model meliputi, uji asumsi white noise dengan statistik uji Ljung-Box dan uji asumsi normalitas dengan statistik uji Jarque Bera. Pada model yang telah memiliki parameter signifikan akan dilakukan uji asumsi white noise menggunakan uji Ljung-Box dengan hipotesis sebagai berikut.

H₀: residual white noise
H₁: residual tidak white noise

      Hasil uji Ljung-Box pada sub-bab 2.7 menunjukkan bahwa model SARIMA(1,1,0)(0,2,0)₁₂ dan SARIMA(0,1,1)(0,2,0)₁₂ memiliki p-value berturut-turut adalah 0.3215 dan 0.3259 yang lebih besar dari alpha (0.05) maka terima H₀. Dengan demikian, dapat disimpulkan bahwa model SARIMA(1,1,0)(0,2,0)₁₂ dan SARIMA(0,1,1)(0,2,0)₁₂ memiliki residual yang white noise. Karena uji asumsi white noise telah terpenuhi, maka selanjutnya dilakukan pengujian asumsi normalitas residual dengan hipotesis sebagai berikut.

H₀: residual berdistribusi normal
H₁: residual tidak berdistribusi normal

      Berdasarkan hasil uji Jarque Berra pada sub-bab 2.7, p-value kedua model berturut-turut adalah 0.008 dan 0.0005 yang lebih kecil dari alpha (0.05) maka tolak H₀. Dengan demikian, residual model SARIMA(1,1,0)(0,2,0)₁₂ dan SARIMA(0,1,1)(0,2,0)₁₂ tidak berdistribusi normal. Namun, menurut Dalil Limit Pusat yang menyatakan bahwa kurva distribusi sampling (untuk ukuran sampel 30 atau lebih) akan berpusat pada nilai parameter dan akan memiliki semua sifat-sifat distribusi normal. Dari teorema tersebut dapat ditarik kesimpulan bahwa data berdistribusi normal karena jumlah observasi lebih dari 30.

3.1.6 Pemilihan Model Terbaik

      Setelah lolos melewati berbagai uji, maka dalam memilih model terbaik untuk meramalkan periode selanjutnya adalah dengan cara menentukan nilai Akaike’s Information Criterion (AIC) setiap model. Model dengan nilai AIC terkecil adalah model terbaik yang dipilih. Berdasarkan hasil analisis pada sub-bab 2.8 Model Terbaik, model SARIMA(0,1,1)(0,2,0)₁₂ memiliki nilai AIC sebesar 600.6674, nilai tersebut lebih kecil dibandingkan nilai AIC model SARIMA(1,1,0)(0,2,0)₁₂ yaitu 602.4252. Dengan demikian, dapat disimpulkan model terbaik yang terpilih adalah SARIMA(0,1,1)(0,2,0)₁₂.

3.1.7 Peramalan

      Peramalan jumlah barang angkutan kereta api di Pulau Jawa untuk 12 periode kedepan menggunakan model terbaik yang telah terpilih yaitu SARIMA(0,1,1)(0,2,0)₁₂. Nilai koefisien parameter yang diperoleh pada sub-bab 2.6 yaitu, MA(1) atau \(\hatθ_1\) adalah -0.5816. Bentuk persamaan model SARIMA(0,1,1)(0,2,0)₁₂ yang terbentuk adalah sebagai berikut.

\[ Y_t=Y_{t-1}+2(Y_{t-12}-Y_{t-13})-Y_{t-24}+Y_{t-25}+0.5816e_{t-1}+e_t \]

      Hasil peramalan untuk 3 periode kedepan menggunakan model SARIMA(0,1,1)(0,2,0)₁₂ adalah sebagai berikut.

1. November 2021: 840.8195
2. Desember 2021: 937.8195
3. Januari 2022 : 1067.8195

Untuk melihat tren hasil prediksi tersebut dapat dilihat pada plot peramalan di sub-bab 2.9. Pada bulan November 2021, jumlah prediksi barang angkutan kereta api sebanyak 841 barang, jumlah tersebut menurun dari bulan sebelumnya yaitu pada bulan Oktober 2021 yang berjumlah 1078 barang. Namun, jumlah barang angkutan kereta api akan cenderung meningkat kembali hingga awal tahun 2022 yaitu di bulan Januari sebanyak 1068 barang.

3.2 Pembahasan

      Data runtun waktu jumlah barang angkutan melalui transportasi kereta api di Pulau Jawa selama 70 bulan dari bulan Januari 2016 hingga Oktober 2021 diprediksikan akan mengalami penurunan pada bulan November 2021 dan akan kembali meningkat pada periode selanjutnya. Model yang digunakan untuk peramalan adalah SARIMA(0,1,1)(0,2,0)₁₂. Model tersebut dapat dikatakan layak untuk peramalan karena telah memiliki parameter yang signifikan dan memenuhi asumsi residual yang white noise. Selain itu, model tersebut dinyatakan sebagai model terbaik karena memiliki nilai AIC terkecil diantara model SARIMA lain yang telah memenuhi uji diagnostik sisaan model.
      Selama pandemi COVID-19 hingga ditetapkannya PPKM darurat, PT Kereta Api Indonesia mencoba tetap stabil walapun pandemi berdampak pada pengurangan jumlah penumpang, namun pelayanan terhadap berbagai komoditas angkutan barang tetap dijalankan. Berdasarkan hasil peramalan yang diperoleh yaitu November 2021 hingga Januari 2022 jumlah angkutan barang kereta api di Pulau Jawa akan cenderung kembali meningkat meskipun diawali denan penurunan dari bulan Oktober 2021. Meningkatnya jumlah angkutan barang tersebut kemungkinan disebabkan oleh semakin berkurangnya kasus COVID-19 yang menandakan semakin meredanya pandemi sehingga permintaan pada layanan angkutan barang diprediksikan akan meningkat. Dengan demikian, pihak PT Kereta Api Indonesia dapat mempertimbangkan hasil peramalan yang diperoleh untuk dapat mengoptimalkan layanan angkutan barang yang secara tidak langsung dapat membantu pemulihan perekonomian di Indonesia pasca pandemi.

4 DAFTAR PUSTAKA

Biomantara, K., & Herdiansyah, H. (2018). Peran Kereta Api Indonesia (KAI) sebagai Infrastruktur Transportasi Wilayah Perkotaan. Jurnal Humaniora Bina Sarana Informatika, 4-7.

Durrah, F. I., Yulia, & Parhusip, T. P. (2018). Peramalan Jumlah Penumpang Pesawat Di Bandara Sultan Iskandar Muda Dengan Metode SARIMA (Seasonal Autoregressive Integrated Moving Average). Journal of Data Analysis, 3.

Fahrudin, R., & Sumitra, I. D. (2020). Peramalan Inflasi Menggunaan Metode SARIMA dan Single Exponential Smoothing (Studi Kasus: Kota Bandung). 112.

Firdaus, M. (2020). 8Aplikasi Ekonometrika dengan E-Views, Stata, dan R. Bogor: PT Penerbit IPB Press.

Gujarati, D. N. (2006). Dasar-dasar Ekonometrika Edisi Ketiga. Jakarta: Erlangga.

Lailani, P. (2019). Analisis Runtun Waktu Peminat Program Studi Pendidikan Matematika FKIP Universitas Jember Menggunakan Metode Exponential Smoothing. 2.

Muhyiddin. (2020). COVID-19, New Normal dan Perencanaan Pembangunan di Indonesia. The Indonesian Journal of Development Planning, 245.

Munawaroh, S. (2010). Analisis Model ARIMA Box-Jenkis pada Data Fluktuasi Harga Emas. 27.

Nasir, W. Y. (2015). Peramalan Jumlah Penumpang dari Pelayaran Dalam Negeri di Pelabuhan Kota Makassar Menggunakan Metode Seasonal Autoregressive Integrated Moving Average (SARIMA). 16-27.

Prasetya, H., & Lukiastuti, F. (2009). Manajemen Operasi. Yogyakarta: MedPress.

Prasidi, A., & Rifni, M. (2020). Kapasitas Infrastruktur dan Fasilitas pada Kereta Api Angkutan Barang dan Logistik. Jurnal Logistik Indonesia, 32.

Putri, D. M., & Aghsilni. (2019). Estimasi Model Terbaik Untuk Peramalan Harga Saham PT. Polychem Indonesia Tbk. dengan ARIMA. MAP Journal, 6-7.

Samawati, P. (2019). Demonopolisasi PT KAI (Persero) Dan PT Pelindo (Persero) Penguatan Sistem Ekonomi Demokrasi. Mmbar Hukum, 316.

Sihombing, Y. F. (2011). Penerapan Kebijakan Revitalisasi Perkeretaapian dan Implikasinya terhadap Perekonomian Indonesia: Pendekatan Pengganda Social Accounting Matrix (Periode 2005-2010). 4.

Sitepu, R. K.-K., & Sinaga, B. (2018). Aplikasi Model Ekonometrika Estimasi, Simulasi, dan Peramalan Menggunakan SAS 9.2. Bogor: PT Penerbit IPB Press.

Supriadi, I. (2020). Metode Riset Akuntansi. Yogyakarta: Deepublish.

Syarifuddin, & Pratomo, W. A. (2013). Efektivitas Penggunaan ARIMA dan VAR dalam Memproyeksi Permintaan Kredit di Indonesia. 3.

Utomo, P., & Fanani, A. (2020). Peramalan Jumlah Penumpang Kereta Api di Indonesia Menggunakan Metode Seasonal Autoregressive Integrated Moving Average (SARIMA). Jurnal Mahasiswa Matematika ALGEBRA, 169-178.

Utomo, S. H. (2009). Jalan Rel. Yogyakarta: Beta Offset Yogyakarta.

Wei, W. W. (2006). Time Series Analysis Univariate and Multivariate Methods. United States of America : Greg Tobin.