1 PENDAHULUAN
1.1 Latar Belakang
Statistika ialah ilmu yang berhubungan dengan pengumpulan, pengolahan, dan penarikan kesimpulan juga pengambilan keputusan dari suatu data pada berbagai bidang. Suatu peramalan terkadang diperlukan untuk menentukansebuah kebijakan untuk masa mendatang pada berbagai bidang kehidupan seperti ekonomi, sosial, teknologi, bahkan kesehatan. Penerapan ilmu statistika dapat digunakan dalam melakukan peramalan. Salah satu metode analisis statistika yang dapat diterapkan dalam peramalan adalah analisis deret waktu.
Data berdasarkan waktu pengambilan dapat dibedakan menjadi dua jenis, yakni data cross section atau data yang dikumpulkan dalam satu waktu tertentu, dan data time series atau data yang diambil selama kurun waktu tertentu (Murray dan Larry, 2007). Dalam penerapan analisis deret waktu jenis data yang digunakan yakni data time series. Analisis deret waktu adalah analisis untuk memodelkan data time series serta meramalkan nilai data pada masa mendatang menggunakan data masa lampau.
Salah satu pemodelan yang dapat digunakan untuk memodelkan data deret waktu (time series) adalah model Autoregressive Integrated Moving Averange (ARIMA). Model ARIMA adalah model yang mengabaikan independensi variabel dalam membuat peramalan dan mengansumsikan bahwa data stationer. Peramalan dapat digunakan sebagai pertimbangan dalam membuat keputusan atau kebijakan karena efektif atau tidaknya suatu keputusan bergantung pada beberapa faktor yang tidak dapat dilihat pada waktu keputusan diambil (Aswi dan Sukmana, 2006).
Kereta Api merupakan salah satu pilihan akomodasi saat ingin bepergian. PT. Kereta Api Indonesia (Persero) adalah Badan Usaha Milik Negara Indonesia yang menyelenggarakan jasa angkutan kereta api. Layanan PT. Kereta Api Indonesia (Persero) meliputi angkutan penumpang dan barang.PT Kereta Api Indonesia telah melakukan berbagai peningkatan pelayanan seperti pelarangan penumpang berdiri, pelarangan penjual asongan, dan penambahan AC pada setiap kereta kelas ekonomi. Oleh karena itu perlunya dilakukan peramalan jumlah penumpang sehingga membantu PT KAI mengantisipasi peningkatan penumpang kereta api.
Berdasarkan penjelasan tersebut, maka penelitian ini akan memodelkan jumlah penumpang kereta api di Pulau Jawa menggunakan metode Autoregressive Interated Moving Averange (ARIMA). Data yang digunakan adalah jumlah penumpang kereta api di Pulau Jawa. Hasil penelitian akan mengetahui model ARIMA terbaik dan diharapkan mampu memberi informasi kepada PT.KAI tentang jumlah penumpang kereta api di Pulau Jawa. Sehingga dari pihak PT KAI dapat mengantisipasi peningkatan penumpang kereta api dengan metode time series, khususnya metode ARIMA.
1.2 Statistika Deskriptif
Statisika merupakan kumpulan teknik atau metode yang berguna untuk membuat keputusan tentang sebuah proses atau populasi berdasarkan analisis informasi yang terkandung dalam contoh dari populasi itu (Montgomery, 2009).Metode statistika dibagi menjadi dua, yaitu statistika deskriptif dan statistika inferensial (statistik induksi). Statistika deskriptif adalah metode-metode yang berkaitan dengan pengumpulan dan penyajian suatu data sehingga memberikan informasi yang berguna (Walpole, 1992). Statistika deskriptif hanya memberikan informasi mengenai data dan tidak menarik kesimpulan mengenai apapun tentang data induknya yang lebih besar. Penyajian data biasanya dalam bentuk tabel, grafik, maupun diagram. Dalam statistika deskriptif juga menunjukkan ukuran pemusatan dan penyebaran data. Ukuran pemusatan data terdiri atas rata-rata, median, dan modus, sedangkan ukuran penyebaran data terdiri atas ragam, simpangan baku, jangkauan dan lain sebagainya.
1.3 Model ARIMA
Model deret waktu yang popular adalah model Autoregresive Integrated Moving Averange (ARIMA) yang dikembangkan oleh George E. P. Box dan Gwilyam M. Jenkins. Identifikasi model dapat dilihat dari hasil ACF dan PACF deret waktu. ARIMA sangat baik dalam ketepatan peramalan jangka pendek, tetapi kurang tepat dalam peramalan jangka panjang. ARIMA adalah model yang mengabaikan variabel independen dalam membuat sebuah peramalan dan model yang mengansumsikan data harus stationer (Wei, 1990).
Tujuan model ARIMA adalah menentukan hubungan statistik yang baik antar variabel yang diramalkan dengan nilai masa lalu variabel tersebut sehingga peramalan dapat dilakukan dengan model tersebut. Model ARIMA merupakan gabungan dari model Autoregressive (AR) dan model Moving Average (MA).
1.3.1 Model Autoregressive (AR)
Model Autoregressive (AR) adalah model yang menunjukkan \(Y_t\) sebagai fungsi linier dari \(Y_t\) pada periode sebelumnya. Model AR yang memiliki ordo 𝑝dinotasikan dengan AR (𝑝), dengan bentuk umum dinyatakan dalam persamaan (1).
\[ Y_{t} = \phi_1Y_{t-1} + \phi_2Y_{t-2}+...+\phi_pY_{t-p}+a_t \tag{1} \]
1.3.2 Model Moving Average (MA)
Model Moving Average (MA) adalah model deret waktu dengan karateristik data periode sekarang merupakan kombinasi linier dari white noise atau nilai kesalahan periode sebelumnya dengan bobot tertentu. Model MA yang memiliki ordo 𝑞 dinotasikan dengan MA (𝑞), dengan bentuk umum dinyatakan dalam persamaan (2).
\[ Y_{t} = a_t - \theta_1Y_{t-1} - \theta_2Y_{t-2} -...- \theta_pY_{t-p} \tag{2} \]
1.3.3 Model Autoregressive Moving Averange (ARMA)
Model Autoregresive Moving Averange (ARMA) merupakan gabungan dari model AR dengan ordo 𝑝 dan MA dengan ordo 𝑞. Sehingga bentuk umum model ARMA (𝑝, 𝑞) dinyatakan dalam persamaan (3).
\[ Y_{t} = \phi_1Y_{t-1} + \phi_2Y_{t-2}+...+\phi_pY_{t-p}+a_t - \theta_1Y_{t-1} - \theta_2Y_{t-2} -...- \theta_pY_{t-p} \tag{3} \]
1.3.4 Model Autoregressive Interated Moving Averange (ARIMA)
Jika terdapat selisih pada deret waktu dan dilakukan analisis dengan model ARMA (𝑝, 𝑞) pada data tersebut, maka dapat dinyatakan deret waktu yang asli adalah ARIMA (𝑝, 𝑑, 𝑞). Artinya model ARIMA digunakan apabila data tidak stationer, maka perlu dilakukan pembedaan (differencing/ d) karena kata Integrated mengacu pada pembedaan. Bentuk umum model ARIMA (𝑝, 𝑑, 𝑞) dinyatakan dalam persamaan (4).
\[ \Delta Y_{t} = \phi_1\Delta Y_{t-1} + \phi_2\Delta Y_{t-2}+...+\phi_p\Delta Y_{t-p}+ a_t - \theta_1\Delta Y_{t-1} - \theta_2\Delta Y_{t-2} -...- \theta_p\Delta Y_{t-p} \tag{4} \]
1.3.5 Model Seasonal Autoregressive Interated Moving Averange (SARIMA)
Musiman adalah kecenderungan mengulangi pola tingkah gerak dalam periode musiman, umumnya satu tahun pada data bulanan. Model ARIMA Musiman merupakan model ARIMA yang digunakan untuk menyelesaikan data time series musiman yang terdiri dari dua bagian, yaitu bagian tidak musiman (non-musiman) dan bagian musiman. Bagian non-musiman dari metode ini adalah model ARIMA. Secara umum bentuk model ARIMA musiman atau ARIMA (p,d,q)(P,D,Q)S sebagai berikut.
\[ \phi_p(B)\Phi_P(B^S)(1-B)^d(1-B^S)^DY_t = \theta_q(B)\Theta_Q(B^S) a_t \tag{5} \]
1.4 Stationeritas
Stationeritas adalah asumsi yang sangat penting dalam deret waktu. Suatu deret pengamatan dikatakan stationer apabila proses tidak berubah dengan berjalannya waktu. Artinya rata-rata deret pengamatan konstan di sepanjang waktu.Dalam sebuah data nyata memungkinkan data tidak stationeritas karena nilai rata-rata maupun nilai ragam yang tidak konstan. Jika data yang dimiliki tidak stationer dapat dimodifikasi agar data menjadi stationer.
Kita dapat memodifikasi varian agar stationer dengan mengubah data yakni:
- Jika simpangan baku suatu deret sebanding dengan level, maka mengambil logaritma natural agar deret waktu konstan dalam varian.
- Jika varian deret asli sebanding dengan level, maka mengambil the square root akan memodifikasi varian menjadi konstan.
Apabila data tidak stationer dalam varian, maka dilakukan transformasi menggunakan Transformasi Box-Cox agar varian menjadi konstan (Myers, 1990). Data yang telah stationer dapat diliat dari nilai 𝜆 = 1 (Wei, 1990). Stationer atau tidak suatu data deret waktu dalam mean dapat dilihat dari plot ACF dan PACF. Dilihat plot ACF, jika nilai-nilai dari autokorelasi cenderung turun lambat menuju nilai nol, ketidakstationeran suatu data dalam mean dapat diatasi dengan proses yang bernama Differencing. Untuk menguji kestationeritasan terhadap rata-rata dapat pula menggunakan Uji Akar Unit atau Unit Root Test yaitu dengan melihat nilai Augmented Dickey-Fuller (ADF).
1.5 Diagnostik Model
Diagnostik model digunakan untuk menguji apakah model telah terpilih dengan benar. Uji kesesuaian model dilakukan dengan uji white noise pada sisaan. Asumsi white noise pada sisaan menunjukkan bahwa tidak ada autokorelasi sehingga sisaan bersifat acak atau tidak memiliki pola tertentu (rata-rata nol dan ragam konstan). Jika sisaan bersifat acak, maka dapat dikatakan bahwa model yang terpilih adalah model yang terbaik. Model dapat diuji dengan uji Ljung-Box (Wei, 1990).
1.6 Uji Asumsi Normalitas Sisaan
Uji asumsi normalitas sisaan digunakan untuk mengetahui apakah sisaan menyebar secara normal atau tidak. Pengujian asumsi normalitas dapat dilakukan dengan beberapa metode, salah satunya yaitu Lilliefors test. Dalam pengujian menggunakan Lilliefors test harus memenuhi beberapa syarat, yaitu (Sudjana, 2005)
- Pengamatan \(y_1,y_2,...,y_n\) dijadikan bilangan baku \(z_1,z_2,...,z_n\)
- Untuk setiap bilangan baku menggunakan distribusi normal baku, kemudian dilakukan perhitungan peluang \(𝐹(z_i) =𝑃(𝑦 ≤ z_i)\).
- Selanjutnya dihitung proporsi \(z_1,z_2,...,z_n\) yang lebih kecil atau sama dengan \(z_i\). Jika proporsi dinyatakan dalam \(S(z_i)\).
- Hitung selisih dari \(F(z_i) − 𝑆(z_i)\) kemudian ditentukan nilai mutlaknya.
- Ambil harga yang paling besar diantara harga-harga multak tersebut.
- Apabila \(L_{hitung} < L_{tabel}\), maka dapat dinyatakan bahwa sisaan berdistribusi normal.
1.7 Data
Data yang digunakan dalam penelitian ini adalah data sekunder yang bersumber dari Badan Pusat Statistik Indonesia. Data yang digunakan adalah data jumlah penumpang kereta api di Pulau Jawa pada bulan Januari 2016 sampai Februari 2020.
2 SOURCE CODE
2.1 Library yang Dibutuhkan
> library(TSA)
> library(tseries)
> library(lmtest)
> library(forecast)
> library(tidyverse)
> library(rmarkdown)
> library(readxl)2.2 Memanggil Data
> Penumpang_Kereta <- read_excel("C:/Users/Asus/Downloads/Penumpang Kereta.xlsx")
> Data<-as.numeric(Penumpang_Kereta$Jumlah)
> data<-ts(Data, start=c(2016,1), frequency = 12)
> paged_table(as.data.frame(Penumpang_Kereta))Function read_excel() digunakan untuk membaca dokumen excel yang berisi data yang akan dianalisis. Argumen yang diisikan dalam function adalah letak dokumen excel disimpan. Function as.numeric() adalah function yang digunakan untuk mengubah type data menjadi numeric. Argumen yang diisikan dalam fuction adalah data yang ingin diubah. Function paged_table() digunakan untuk membuat sebuah tabel dalam RStudio.
2.3 Plot Deret Waktu
> plot(data, xlab="Periode", ylab="Jumlah", main="Plot Time Series") 2.4 Uji Stationeritas
> FitAR::BoxCox.ts(data)> adf.test(data)
Augmented Dickey-Fuller Test
data: data
Dickey-Fuller = -1.1846, Lag order = 3, p-value = 0.9005
alternative hypothesis: stationary> acf(data, lag=100)> pacf(data, lag=100)> adf.test(diff(data))
Augmented Dickey-Fuller Test
data: diff(data)
Dickey-Fuller = -6.3315, Lag order = 3, p-value = 0.01
alternative hypothesis: stationaryFunction BoxCox.ts() digunakan untuk menguji kestationeran data terhadap ragam. Dan function adf.test() digunakan untuk menguji kestationeran data terhadap rata-rata. Function diff() digunakan untuk operasi differencing pada data deret waktu yang belum stationer dalam rata-rata.
2.5 Identifikasi Model
> acf(diff(data))> pacf(diff(data))2.6 Estimasi Parameter
> m1=arima(data, order=c(1,1,0), seasonal=list(order=c(0,1,0), period=12), method="ML")
> lmtest::coeftest(m1)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.45283 0.14978 -3.0233 0.002501 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> m2=arima(data, order=c(0,1,1), seasonal=list(order=c(0,1,0), period=12), method="ML")
> lmtest::coeftest(m2)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ma1 -0.54288 0.11596 -4.6817 2.845e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> m3=arima(data, order=c(1,1,1), seasonal=list(order=c(0,1,0), period=12), method="ML")
> lmtest::coeftest(m3)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.14132 0.24033 -0.5880 0.55653
ma1 -0.46282 0.18845 -2.4559 0.01405 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Function arima() digunakan untuk melakukan estimasi parameter pada model ARIMA. Argumen yang diisikan dalam function adalah variabel x atau data deret waktu yang akan dimodelkan, dan order=c(p,d,q) yang merupakan dugaan orde ARIMA berdasarkan identifikasi model. Function coeftest() digunakan untuk menguji tingkat signifikansi parameter berdasarkan hasil estimasi.
2.7 Diagnostik Model
> aic.model<-data.frame(Nama=c("m1","m2","m3"), Model=c("ARIMA(1,1,0)","ARIMA(0,1,1)","ARIMA(1,1,1)"), AIC = c(m1$aic, m2$aic, m3$aic))
> aic.model
Nama Model AIC
1 m1 ARIMA(1,1,0) 628.8802
2 m2 ARIMA(0,1,1) 625.6377
3 m3 ARIMA(1,1,1) 627.3161> checkresiduals(m2)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(0,1,0)[12]
Q* = 11.397, df = 9, p-value = 0.2494
Model df: 1. Total lags used: 10
Fumction checkresiduals() digunakan melihat apakah model yang terbentuk atau terpilih telah sesuai dan layak digunakan untuk tahap peramalan data.
2.8 Forecast
> ramalan<-forecast::forecast(data, model=m2, h=12)
> ramalan
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Mar 2020 35032.76 33621.90 36443.61 32875.04 37190.48
Apr 2020 35070.76 33519.48 36622.03 32698.29 37443.23
May 2020 34478.76 32798.76 36158.76 31909.42 37048.09
Jun 2020 34225.76 32426.22 36025.30 31473.60 36977.91
Jul 2020 38267.76 36356.14 40179.37 35344.19 41191.32
Aug 2020 34506.76 32489.28 36524.24 31421.29 37592.22
Sep 2020 34579.76 32461.70 36697.81 31340.47 37819.04
Oct 2020 35778.76 33564.69 37992.83 32392.63 39164.88
Nov 2020 35192.76 32886.67 37498.85 31665.90 38719.62
Dec 2020 36674.76 34280.19 39069.33 33012.57 40336.94
Jan 2021 33436.76 30956.85 35916.66 29644.07 37229.44
Feb 2021 31643.76 29081.36 34206.15 27724.91 35562.60
> plot(ramalan)Function forecast() digunakan untuk melakukan peramalan atau prediksi beberapa waktu ke depan. Argumen yang diisikan dalam function adalah data, model dan h. Model yang diisikan adalah model yang telah terpilih dari diagnostik model, dan h digunakan untuk seberapa banyak prediksi yang ingin dilakukan.
3 HASIL DAN PEMBAHASAN
3.1 Statistika Deskriptif
> mean(data)
[1] 32631.48Berdasarkan Data jumlah penumpang kereta api perhari, diperoleh nilai rataan sebesar 32631 penumpang perhari. Dapat dilihat pada plot data time series terdapat indikasi sebuah pola tren yang cenderung naik setiap periodenya, di mana mengindikasikan pula bahwa data deret waktu tersebut tidak stationer. Dapat pula dilihat pada plot terdapat pola musiman yang terjadi setiap tahunnya. Hal ini dikarenakan pada Hari Raya Lebaran penumpang kereta api akan melonjak kemudian stabil dan kembali terulang disetiap tahunnya.
3.2 Uji Stationeritas
Uji yang pertama dilakukan adalah pengujian stationeritas data deret waktu terhadap ragam. Dengan menggunakan function BoxCox() didapatkan hasil \(\lambda=1\), yang artinya data deret waktu jumlah penumpang kereta api telah memenuhi kstationeritasan terhadap ragam dan tidak perlu dilakukan transformasi.
Selanjutnya dilakukan uji stationeritas data deret waktu terhadap rata-rata menggunakan uji Augmented Dickey-Fuller (ADF).
Dengan Hipotesis:
\(H_0\): data tidak stationer
\(H_1\): data stationer
Dengan menggunakan function adf.test() didapatkan hasil nilai p-value sebesar 0.9005, nilai ini \(>\alpha\) yang artinya Terima \(H_0\), sehingga data jumlah penumpang kereta api tidak memenuhi asumsi kestationeritasan terhadap rata-rata dan perlu dilakukan differencing. Ketidak-stationeritasan ini juga dapat dilihat pada plot ACF data yakni terlihat pola menurun secara eksponensial yang berarti data tidak stationer terhadap rata-rata. Dan indikasikan adanya musiman karena terlihat pada plot ACF dan PACF yang memiliki pola tertentu.
Setelah dilakukan differencing sebanyak satu kali, data kembali di uji dengan uji yang sama seperti sebelumnya. Dan didapatkan nilai p-value sebesar 0.01 nilai ini \(<\alpha\) yang artinya Tolah \(H_0\), sehingga dapat disimpulkan data jumlah penumpang kereta api telah memenuhi asumsi kestationeritasan terhadap rata-rata. Dan dapat dilanjutkan untuk tahap berikutnya yakni identifikasi model.
3.3 Identifikasi Model
Identifikasi model ARIMA dapat dilakukan dengan melihat plot ACF dan PACF dari data yang telah memenuhi asumsi Stationeritas. Berdasarkan Plot ACF dan PACF data stationer, dapat terlihat bahwa lag yang signifikan pada grafik ACF adalah lag 11 sehingga model MA maksimal berordo 𝑞 = 7, dan lag pada grafik PACF adalah lag 10 sehingga model AR maksimal berordo 𝑝 = 10. Karena data dilakukan differencing sebanyak satu kali maka ordo d = 1.
Pada identifikasi kali ini akan menggunakan model tentative ARIMA (1,1,0), ARIMA (0,1,1) dan ARIMA (1,1,1) karena kedua plot ACF dan PACF terjadi lag signifikan pada lag pertama juga. Pada Identifikasi awal terdapat pola musiman, sehingga digunakan model musiman ARIMA \((0,1,0)_{12}\).
3.4 Estimasi Parameter
Dilakukan estimasi parameter pada 3 model tentative ARIMA, lalu dilanjutkan dengan pengujian signifikansi parameter untuk menentukan model yang terbaik. Model yang terbaik adalah model yang di mana semua parameter dalam model terbukti signifikan. Dari ketiga model tentative terdapat 2 model yang memenuhi signifikansi parameter, yakni model ARIMA \((1,1,0)(0,1,0)_{12}\) dan model ARIMA \((0,1,1)(0,1,0)_{12}\). Selanjutnya akan dipilih dari kedua model tersebut yang terbaik.
3.5 Diagnostik Model
Pemilihan model terbaik dapat menggunakan nilai AIC dari model yang paling kecil. Dari ketiga model tentative nilai AIC yang terkecil adalah model ARIMA \((0,1,1)(0,1,0)_{12}\) dengan nilai AIC sebesar 625.6377.
Selanjutnya dilakukan uji Ljung-Box untuk mengetahui apakah model yang terpilih memenuhi asumsi white noise dari sisaan modelsekaligus menguji kelayakan model untuk tahap selanjutnya.
Dengan Hipotesis:
\(H_0\): sisaan white noise
\(H_1\): sisaan tidak white noise
Dengan menggunakan function checkreduals() didapatkan hasil nilai p-value sebesar 0.2494, nilai ini \(>\alpha\) yang artinya Terima \(H_0\), dapat disimpulkan bahwa sisaan model memenuhi asumsi white noise dan model layak digunakan untuk tahap peramalan data.
Selanjutnya dapat dilihat pada beberapa grafik dari sisaan untuk melihat normalitas sisaan model. Pada plot sisaan terlihat bahwa data telah menyebar sekitar nilai 0 dan pada plot ACF sisaan dapat dilihat bahwa tidak terdapat lag yang signifikan. Hal tersebut telah membuktikan bahwa sisaan model telah memenuhi uji normalitas sisaan.
3.6 Forecasting
Berdasarkan beberapa tahap yang telah dilakukan didapatkan model ARIMA yang terpilih untuk melakukan peramalan adalah model ARIMA \((0,1,1)(0,1,0)_{12}\). Dengam persamaan sebegai berikut.
\[ (1-B)^1(1-B^{12})^1Y_t = (1-(-0.45283))a_t \\ (1-B^{12}-B+B^{13})Y_t=a_t+0.45283a_t \\ Y_t-BY_t-B^{12}Y_t+B^{13}Y_t=a_t+0.45283a_t \\ Y_t = Y_{t-1}+Y_{t-12}-Y_{t-13}+a_t+0.45283a_t \]
Selanjutnya dilakukan peramalan atau prediksi untuk 12 bulan kedepan pada data jumlah penumpang kereta api di Pulau Jawa. Dapat dilihat pada Plot Forecast di mana dari bulan Februari 2020 terjadi kenaikan pada bulan Maret 2020, dan kenaikan yang signifikan terjadi pada bulan Juli 2020 hal ini dikarenakan banyaknya masyarakat yang melakukan mudik pada Hari Raya Lebaran. Pola prediksi 12 bulan kedepan mirip dengan pola-pola yang terjadi sebelumnya.
4 DAFTAR PUSTAKA
Aswi dan Sukmana. 2006. Analisis Deret Waktu. Makassar: Andira Publisher.
Murray R. Spiegel and Larry J. Stephens. 2007. Schaum’s Outlines Teori Dan Soal-Soal Statistik. Jakarta: Erlangga.
Montgomery, Douglas. 2009. Introduction to Statistical Quality Control. New York: John Wiley & Sons, Inc.
Myers, R. H. 1990. Clasican and Modern Regression with Application. 2nd edn. Psw-Kent Publishing Com.
Sudjana. 2005. Metode Statistika. Jakarta: Tristo.
Walpole, R. E. 1992. Pengantar Statistika. 3rd edn. Jakarta: PT Gramedia Pustaka Utama.
Wei, W. W. (1990). Time Series Analysis: Univariate and Multivariate Methods. Addison Wesley Publishing Company, Inc.