## PREDIKSI HARGA MINYAK DUNIA DENGAN METODE AUTOREGRESSIVE FRACTIONALLY INTEGRATED MOVING AVERAGE (ARFIMA)
1 Dimas Kevin Natanael, 2 Diah Safitri, 3 Suparti 1,2,3 Departemen Statistika FSM Universitas Diponegoro Email : dimaskevin25@gmail.com

## ABSTRACT
Autoregressive Fractionally Integrated Moving Average (ARFIMA) model is a development of  the  ARIMA  model.  The  advantage  of  the  ARFIMA  method  is  the  non-integer differentiation value so that it can overcome long memory effect that cannot be solve with the usual ARIMA method. Non-integer differential values can be estimated with a binomial expansion  approach  which  is  an  infinite  weighted  sum  of  past  values  to  solve  the  long memory effect that arises. Some of the advantages of using the ARFIMA model iscapable of modeling high  changes  in the long  term (long  term persistence),  be able to  explain long- term and short-term correlation structures at the same time, to provide models with simple parameters (parsimony) for data with memory long term and short term.Data of world oil price contain long memory effect, then used ARFIMA method to get the best model. The best model obtained is the ARMA([1,7]; 0)  model with the differential value is 0,48937, then the model can be written into ARFIMA ([1,7]; d; 1). The best model chosen has an MSE value of 0,44 and a MAPE value of 3,32%.
Keywords : Sea Passengers, ARIMA Box-Jenkins, Calendar Variation, ARIMAX

## PENDAHULUAN
Indonesia  merupakan  negara  yang berlimpah  akan  berbagai  jenis  Sumber Daya  Alam  seperti  batu  bara,  tembaga, nikel,  pasir  besi,  biji  timah,  dan  lainnya, tak  terkecuali  minyak  mentah  dan  gas bumi.Kementerian  Energi  dan  Sumber Daya Mineral Republik  Produksi minyak mentah Indonesia kini semakin menurun. Kondisi  yang  bertolak  belakang  antara kinerja produksi dan konsumsi minyak ini membuat  Indonesia  mengalami  defisit, yaitu sekitar 800 ribu barel per hari pada tahun 2006 menjadi sekitar 690 ribu barel per  hari  pada  tahun  2015.  Penurunan tersebut  disebabkan  oleh  penurunan produksi sumur – sumur produksi minyak bumi yang umumnya sudah tua sedangkan produksi sumur baru masih relatif terbatas [14].
Dalam  rangka  memenuhi  kebutuhan konsumsi  bahan  bakar  maka  dilakukan kegiatan import. Meskipun  demikian masalah  kembali muncul  kaitannya dengan  harga  minyak import. Harga minyak import sering mengalami fluktuasi yang mengakibatkan harga sering berubah –  ubah,  sehingga  diperlukan  prediksi harga  minyak  yang  akurat.  Prediksi  ini dibutuhkan sebagai bahan informasi yang kaitannya  dengan  perkiraan  harga  bahan bakar  di  masa  yang  akan  datang.  Data harga minyaksendiri merupakan salah satu data  runtun  waktu  yang  dapat  diprediksi nilainya untuk beberapa tahap kedepan.
Adakalanya  suatu  data  runtun  waktu menunjukkan pola memori jangka panjang  ( long  memory ),  ini  terlihat  dari nilai-nilai  autokorelasi  pada  plot  ACF yang  turun  secara  lambat  untuk  jarak waktu (lag) yang semakin meningkat dan hasil  perhitungan  dari  statistik Husrt  (H)
yang  terletak  dalam  interval  0,5  > H  >1 . Identifikasi  ini  mengindikasikan  bahwa nilai  dari  koefisien  pembedabernilai pecahan,  sehingga  model  yang  paling cocok adalah model ARFIMA
( Autoregressive  Fractionally  Integrated Moving Average ) [1].
Granger  (1980)  adalah  orang  yang pertama  kali  memperkenalkan  model
Autoregressive  Fractionally  Integrated Moving Average (ARFIMA) [5] .
Kemudian  oleh      Hosking  (1981) dilakukan pengkajian terhadap sifat – sifat long  memory dari model ARFIMA. M odel ARFIMA dapat  menjelaskan  deretberkala jangka  pendek  ( short  memory )  maupun deret  berkala  jangka  panjang  ( long memory ).  Metode ARFIMA ini  dapat mengatasi kelemahan pada model ARIMA. ARIMA hanya  dapat  menjelaskan  runtun waku jangka pendek, sedangkan ARFIMA dapat  menjelaskan  baik  jangka  pendek maupun jangka panjang [7].
Dalam  penelitian  ini  akan  dilakukan pemodelan  data  harga  minyak  dunia (basket  price) dengan  menggunakan pendekatan  deret  berkala  memori  jangka panjang ARFIMA . Nilai diferensiasi dapat diestimasi  dengan  metode Geweke  and Porter-Hudak (GPH) dan Exact Maximum
Likelihood  (EML) .Dengan  menggunakan metode
Autoregressive Fractionally
Integrated  Moving  Average  (ARFIMA) , harga  minyak  dunia  beberapa  tahap  ke depan  dapat  diprediksi  berdasarkan  data harga minyak dunia dari periode sebelumnya. Hal iniyang melatar belakangi  penulis  dalam  melakukan penelitian  dengan  judul  “Prediksi  Harga Minyak Dunia dengan Metode Autoregressive  Fractionally  Integrated
Moving Average (ARFIMA)”.
Tahapan untuk pemodelan
ARFIMA adalah  pertama  identifikasi model dengan plot ACF dn PACF. Kedua, estimasi parameter pembeda(d) menggunakan  metode Geweke  Porter Hudak (GPH). Ketiga, estimasi parameter menggunakan  metode  Exact Maximum
Likelihood (EML) . Terakhir, uji signfikansi  parameter  dan  verifikasi residual  model  yang  terdiri  dari  uji indepensi residual menggunakan uji Ljung Box ,  dan  uji  homogenitas  residual menggunakan uji Lagrange Multiplier .

## METODE PENELITIAN

## Sumber Data dan Variabel Penelitian
Jenis  data  yang  digunakan  dalam penulisan  skripsi  ini  adalah  jenis  data sekunder  yaitu  data  tentangharga  minyak dunia (basket  price) yang  diperoleh  dari [17]. Data yang diambil adalah data harian dari  bulan  Oktober  tahun  2015  sampai April  tahun  2017.  Data  diolah  dengan menggunakan beberapa sofware seperti R Command, E-views, Minitab , dan Microsoft excel .

## Metode Analisis
Untuk  mencapai  tujuan  penulisan penelitian  ini,  ditempuh  langkah-langkah sebagai berikut.
1. Membuat plot time series data harga minyak  bumi  untuk  mengetahui apakah data tersebut sudah stasioner atau belum.
2. Melakukan  transformasi  jika  ada data  yang  tidak  stasioner  dalam variansi.
3. Membuat plot ACF dan PACF data yang  telah  ditransformasi  untuk mengetahui  adanya  ketergantungan jangka panjang.
4. Melakukan perhitungan Hurst Exponent dengan  metode rescaled range statistics (R/S) untuk melihat efek long memory.
5. Melakukan  estimasi  nilai  d  dengan metode Geweke and Porter-Hudak.
6. Membuat  plot  PACF  dn  ACF dengan  menggunakan  data  yang telah dilakukan diferensi.
7. Melakukan pemodelan dengan metode ARFIMA.
8. Melakukan uji asumsi untuk melihat apakah  residual  memenuhi  asumsi white  noise ,  kesamaan  varian  dan berdistribusi normal.
9. Model yang telah memenuhi semua asumsi akan dibandingkan berdasarkan nilai MSE yang terkecil.
10. Membuat  ramalan  harga  minyak dunia  untuk  10  tahap  ke  depan
dengan menggunakan model ARFIMA terbaik yang diperoleh.
11. Evaluasi  kriteria  model  dilihat  dari nilai Mean  Absolute  Percentage Error (MAPE).

## HASIL PENELITIAN

## Uji Stasioneritas Data
Gambar 1. Box Plot Transformation
Gambar 2. Trend Analysis Plot
Berdasarkan OutputGambar 1 dan
Gambar  2,dapat  disimpulkan  bahwa  data belum  stasioner  dalam  mean  maupun dalam  varian.  Nilai  dari rounded  value adalah 2. Hal ini menunjukan data belum stasioner dalam varian dan perlu
dilakukan  transformasi.  Sesuai  dengan hasil  nilai  rounded  value  yang  diperoleh, transformasi  yang  digunakan  adalah dengan  cara  mencari  nilai  kuadrat  dari data  yang  digunakan.Lalu,  pada  uji Augmented  Dickey  Fuller ,  kedua return saham  tersebut  mempunyai  nilai pvalue sebesar  0,6950  sehingga  belum  stasioner dalam mean secara uji formal.
Pengujian Pola Long Memory
Gambar 3. Plot Autocorrelaion Function
Berdasarkan  Output  Gambar  3  dapat dilihat  bahwa  lag  yang  muncul  turun lambat  secara  hiperbolik.  Secara  visual dapat dikatakan bahwa data mengandung efek  memory  jangka  panjang.  Namun meskipun  secara  visual  terpenuhi,  perlu dilakukan pengujian secara formal.
Selain dilakukan secara visual, identifikasi  efek  memorijangka  panjang juga  perlu  dilakukan  secara  formal.
Identifikasi  model  secara  formal  dapat dilihat  melalui  perhitungan  nilai Hurst Exponent (H) . Digunakan metode rescaled  range  statistics  (R/S) dengan langkah – langkah sebagai berikut:
1. Menghitung  nilai  rata-rata  data (mean)
𝜇 = 1 𝑛 ∑ 𝑋 𝑖 𝑛 𝑖=1 = 217091.31 400 = 42,72
2. Menghitung  rata-rata  tertimbang dari masing-masing data
𝑌 1 = 𝑋 1 − 𝜇 = 44,48 − 42,72 = 1,76 𝑌 2 = 𝑋 2 − 𝜇 = 43,82 − 42,72
= 1,10
⋮ 𝑌 400 = 𝑋 400 − 𝜇 = 51,9 − 42,72 = 9,18
3. Menghitung  simpangan  kumulatif data 𝑍 1 = ∑ 𝑌 𝑖 𝑡 𝑖=1 = 𝑌 1 = 1,76 𝑍 2 = ∑ 𝑌 𝑖 𝑡 𝑖=1 = 𝑌 1 + 𝑌 2 = 1,76 + (1,10) = 2,86 ⋮ 𝑍 400 = ∑ 𝑌 𝑖 = 𝑡 𝑖=1 𝑌 1 + 𝑌 2 + ⋯ + 𝑌 400 = 1,76 + 1,10 + ⋯ + 9,18 = 1,76
4. Menghitung rentangan data
𝑅 1 = max(𝑍 1 ) − min(𝑍 1 ) = 1,76 − 1,76 = 0 𝑅 2 = max(𝑍 1 , 𝑍 2 ) − min(𝑍 1 , 𝑍 2 ) = 2,86 − 1,76 = 1,10 ⋮ 𝑅 400 = max(𝑍 1 , 𝑍 2 , … , 𝑍 400 ) − min(𝑍 1 , 𝑍 2 , … , 𝑍 400 )
= 1171,30
5. Menghitung  standar  deviasi  dari masing-masing data
𝑆 1 = √ 1 𝑡 ∑(𝑋 𝑖 − 𝜇) 2 𝑡 𝑖=1 = 1,76 𝑆 2 = √ 1 𝑡 ∑(𝑋 𝑖 − 𝜇) 2 𝑡 𝑖=1
= 2,02 ⋮
𝑆 400 = √ 1 𝑡 ∑(𝑋 𝑖 − 𝜇) 2 𝑡 𝑖=1 = 0,09 6. Menghitung rescaled range statistics (R/S)
(𝑅/𝑆) 1 = 𝑅 1 𝑆 1 = 0 1,76 = 0 (𝑅/𝑆) 2 = 𝑅 2 𝑆 2 = 1,10 2,02 = 0,54
⋮ (𝑅/𝑆) 400 = 𝑅 400 𝑆 400 =
1171,30
0,09 = 13339,70
7. Menghitung nilai log rescaled range statistics (R/S)
𝑌 1 = log[(𝑅/𝑆) 1 ] = 1 𝑌 2 = log[(𝑅/𝑆) 2 ] = −0,61 ⋮ 𝑌 400 = log[(𝑅/𝑆) 400 ] = 9,45 8. Menghitung nilai log waktu dari data pengamatan
𝑋 1 = log(1) = 0 𝑋 2 = log(2) = 0,7 ⋮
𝑋 400 = log(400) = 5,99 9. Kemudian dapat dilakukan perhitungan nilai Hurst Exponent(H)
𝐻 = ∑ (𝑋 𝑡 − 𝜇 𝑥 𝑡 )(𝑌 𝑡 − 𝜇 𝑌 𝑡 ) 𝑛 𝑡=1 ∑ (𝑋 𝑡 − 𝜇 𝑥 𝑡 ) 2 𝑛 𝑡=1 = 0,8
Estimasi Parameter Pembeda (d) Setelah  sebelumnya  sudah  diketahui bahwa  data  mengandung  efek long memory , maka digunakan metode ARFIMA dalam  proses  peramalannya. metode ARFIMA menggunakan nilai parameter  pembeda  ( d )  dalam  bentuk pecahan.  Maka  dari  itu  perlu  digunakan metode  khusus  dalam  penentuan  nilai parameter  pembedanya.  Metode  yang digunakan  untuk  menentukan  nilai d adalah  metode  Geweke  Porter-Hudak (GPH).
Dalam  proses  perhitungan  nilai d digunakan  bantuan software  R . Package yang  digunakan  dalam software  R ketika mencari  estimasi  nilai d dengan  metode GPH adalah fraccdiff .
Berdasarkan perhitungan dengan software R (lampiran 2)diperoleh  estimasi  nilai d sebesar 1,086094.
Identifikasi Model
Gambar 4. Autocorelaion Function
Gambar 5. Partial Autocorrelation Function
Setelah  melihat  hasil  dari  identifikasi model  pada  Gambar  4  dan  Gambar  5, dapat  diambil  kesimpulan  bahwa  model yang mungkin adalah:
1. ARFIMA (1; d;0)
7. ARFIMA (1; d; 1) 2. ARFIMA ([7]; d; 0)
8. ARFIMA ([7]; d;
1)
3. ARFIMA ([1,7]; d; 0)
9. ARFIMA ([1,7]; d; 1) 4. ARFIMA (0; d; 1) 10. ARFIMA (1; d;[7]) 5. ARFIMA (0; d;[7]) 11. ARFIMA (1; d; [1,7]) 6. ARFIMA (0; d;[1,7]) 12. ARFIMA ([1,7]; d; [1,7])

## Pemilihan Model Terbaik
Dari  beberapa  model  awal  yang  telah dilakukan  pengujian  asumsi,  kita  dapat menentukan  model  mana  yang  terbaik untuk digunakan dalam proses peramalan nantinya.  Hasil  pengujian  dapat  dilihat dalam Tabel 3:
Peramalan Setelah beberapa pengujian sebelumnya,  kemudian  dapat  diperoleh kesimpulan  bahwa  model  terbaik  yaitu model  ARFIMA([1,7];  d;  0)  yang  dapat dituliskan persamaannya,
ɸ(𝐵)(1 − B) 𝑑 𝑍 𝑡 =   𝑎 𝑡 (1 − ɸ 1 𝐵 − ɸ 7 𝐵 7 )(1 − B) 0,48937 𝑍 𝑡 =   𝑎 𝑡 (𝑍 𝑡 − ɸ 1 𝑍 𝑡−1 − ɸ 7 𝑍 𝑡−7 )(1 − B) 0,48937 = 𝑎 𝑡 (𝑍 𝑡 − ɸ 1 𝑍 𝑡−1 − ɸ 7 𝑍 𝑡−7 ) = 𝑎 𝑡 (1 − B) 0,48937 𝑍 𝑡 = ɸ 1 𝑍 𝑡−1 + ɸ 7 𝑍 𝑡−7 + 𝑎 𝑡 (1 − B) 0,48937 𝑍 𝑡 = 0,751𝑍 𝑡−1 + 0,130467𝑍 𝑡−7 + 𝑎 𝑡 (1 − B) 0,48937
dengan (1 − B) 0,48937 dapat didekati dengan ekspansi binomial (1 − B) 0,48937 = 1 − 0,48937𝐵
− 1 2 0,48937(1 − 0,48937)𝐵 2 − 1 6 0,48937(1 − 0,48937)(2 − 0,48937)𝐵 3 − ⋯ = 1 − 0,48937𝐵 − 0,12494𝐵 2 − 0,06291𝐵 3 − ⋯
Setelah  diperoleh  model  terbaik,  maka model  tersebut  dapat  digunakan  untuk melakukan  proses  peramalan.  Peramalan yang  akan  dilakukan  sebanyak  10  tahap kedepan. Data yang di modelkan sebelumnya  telah  dilakukan  transformasi dalam  bentuk  pangkat  dua,  maka  setelah
diperoleh
hasil peramalan harus
dikembalikan kembali ke bentuk asli dari data  dengan  cara  melakukan  akar  pada data hasil peramalan. Data hasil peramalan  yang  telah  dikembalikan disajikan pada Tabel 4.
Tabel 3. Tabel Hasil Uji Asumsi Model Signifikansi Parameter Independensi Residual Kesamaan Varian MSE 1. ARFIMA (1;d ; 0) Parameter Signifikan Independensi Tidak Terpenuhi Varian Sama 0,044 2. ARFIMA (7;d ; 0) Parameter Signifikan Independensi Tidak Terpenuhi Varian Tidak Sama 0,067 3. ARFIMA ([1,7]; d ; 0) Parameter Signifikan Independensi Terpenuhi Varian Sama 0,044 4. ARFIMA (0;d ; 1) Parameter Tidak Signifikan Independensi Tidak Terpenuhi Varian Tidak Sama 0,185 5. ARFIMA (0; d ; 7) Parameter Tidak Signifikan Independensi Tidak Terpenuhi Varian Tidak Sama 0,151 6. ARFIMA (0; d ; [1,7]) Parameter Tidak Signifikan Independensi Tidak Terpenuhi Varian Tidak Sama 2,071 7. ARFIMA (1; d ; 1) ParameterTidak Signifikan IndependensiTidak Terpenuhi Varian Tidak Sama 0,044 8. ARFIMA (7; d ; 1) Parameter Signifikan Independensi Tidak Terpenuhi Varian Sama 0,053 9. ARFIMA ([1,7]; d ; 1) ParameterTidak Signifikan Independensi Terpenuhi Varian Tidak Sama 0,045 10. ARFIMA (1;d ; 7) Parameter Tidak Signifikan Independensi Tidak Terpenuhi Varian Sama 0,045 11. ARFIMA (1; d ; [1,7]) Parameter Tidak Signifikan Independensi Terpenuhi Varian Sama 0,045 12. ARFIMA ([1,7]; d ; [1,7]) Parameter Tidak Signifikan Independensi Terpenuhi Varian Tidak Sama 0,045
Tabel 4. Peramalan dengan model ARFIMA([1,7]; d; 0) Periode Tanggal Nilai Ramalan Interval Kepercayaan Data Aktual Batas atas Batas bawah 401 4/20/2017 52.29 64.20 35.36 50.49 402 4/21/2017 51.47 62.89 32.89 50.01 403 4/24/2017 51.05 62.72 32.56 49.66 404 4/25/2017 50.42 62.44 32.03 49.22 405 4/26/2017 49.87 62.06 31.26 49.65 406 4/27/2017 49.72 62.54 32.19 48.9 407 4/28/2017 49.28 61.73 30.56 49.33 408 5/1/2017 48.82 62.17 31.43 49.19 409 5/2/2017 48.87 61.99 31.08 49.02 410 5/3/2017 48.8 61.86 30.81 48.38
Dalam  proses  peramalan  ini  diperoleh nilai MSE dan MAPE yang  ditampilkan pada Tabel 5.
Tabel 5. Tabel nilai MSE dan MAPE Model MSE MAPE ARFIMA([1,7]; d; 0) 0,044 3,23% Berdasarkan Tabel 5 dapat dilihat bahwa model ARFIMA([1,7]; d; 0) adalah model  dengan  kinerja  peramalan  yang terbaik  karena  mempunyai  nilai MAPE terkecil.  Jika  dibandingkan  dengan  tabel kriteria  pengukuran MAPE, model  ini masuk dalam kriteria kemampuan peramalan sangat baik karena nilai MAPE yang diperoleh kurang dari 10%.

## KESIMPULAN
Berdasarkan  analisis  dan  pembahaan data  harian  harga  minyak  mentah  dunia dari  periode  Oktober  2015  hingga  April 2017  yang  telah  dijelaskan  sebleumnya, dapat  diperoleh  kesimpulan  sebagai berikut:
1. Hasil  pemeriksaan  terhadap  data menunjukan  bahwa  terdapat long memory  effect yang  ditunjukan dengan  plot  ACF  yang  turun lambat secara hiperbolik dan nilai Hurst sebesar 0,8 .
2. Nilai d yang  diperoleh  dari  hasil estimasi  dengan  metode Geweke
Porter-Hudak (GPH) sebesar
1,08609
3. Hanya  ada  satu  model  yang memenuhi  uji  diagnostik  dari  dua belas model yang diujikan. Model tersebut adalah model
ARFIMA([1,7]; d; 0) dengan nilai estimasi d dengan  metode Exact Maximum  Likelihood diperoleh sebesar 0,48937.
4. Nilai MSE yang  diperoleh  dari model ARFIMA([1,7];  d;  0) sebesar 0,044 dan nilai MAPE dari model ARFIMA([1,7];  d;  0) sebesar 3,23%.
5. Model ARFIMA([1,7];  d;  0) dapat dituliskan persamaannya,
𝑍 𝑡 = 0,751𝑍 𝑡−1 + 0,130467𝑍 𝑡−7 + 𝑎 𝑡 (1 − B) 0,48937

## DAFTAR PUSTAKA
[1]  Cahyandari,  R.,  dan  Erviana,  R.  2015. Peramalan Kurs Jual Uang Kertas Mata Uang Singapore Dollar (SGD) terhadap Rupiah Menggunakan ModelARFIMA ( Autoregressive Fractionally Integrated  Moving  Average)  .Kubik Volume 1, Nomor 1
[2] Caraka, R. E., Yasin, H., dan Wawan, S.  2016.  MODEL  LONG  MEMORY DALAM  MEMPREDIKSI  SUHU.
Sekolah
Tinggi Meteorologi Klimatologi dan Geofisika (STMKG)
[3]  Chatfield,  C.  2000. Time  Series Forecasting. Boca  Raton,  Florida: CRC Press
[4]  Cryer,  J.  D.,  dan  Kung,  S.C.  2008. Time Series Analysis With Applications in R 2 nd edition . Springer Science+Business Media, LLC [5]  Granger,  C.  W.  J.  1980.  LONG
MEMORY  RELATIONSHIPS  AND THEAGGREGATION
OF DYNAMIC  MODELS .  Journal  of Econometrics 14.University
Of
California:  North-Hollland  Publishing Company
[6]  Halimi,  R.,  Anggraeni,  W.,  dan Tyasnurita R. 2013. Pembuatan
Aplikasi Peramalan Jumlah Permintaan Produk  dengan  Metode  Time  Series
Exponential Smoothing HOLTS
Winterr  di  PT.  TELEKOMUNIKASI Indonesia Tbk. ITS: JURNAL TEKNIK POMITS Vol. 1, No. 1
[7]  Hosking,J.  R.  M.    1981.  Fractional Differencing . Biometrika , Vol.68 Page 165-176
[8]  Idris,  S.,  Goejantoro,  R.,  dan  Yuki, N.N. 2014. Pemodelan Dan Peramalan Indeks  Harga  Perdagangan  Besar (IHPB) Dengan Menggunakan ARFIMA (Studi Kasus : IHPB Provinsi Kalimantan Timur bulan Januari 2002
–  Desember  2006  dan  Januari  2009  - September 2013). Jurnal EKSPONENSIAL Volume 5, Nomor 2 [9] Montgomery, D. C., Cheryl, L. J., dan Kulahci, M. 2008. Introduction to time series analysis and forecasting . Hoboken, N.J.: Wiley-Interscience.
[10]  Ooms,  M.  and  J.  A.  Doornik.  1999. Inference and forecasting for fractional autoregressive integrated moving average models, with an application to US and UK inflation . Discussion Paper EI  9947/A,  Econometric  Institute, Erasmus University Rotterdam.
[11]  Rosadi,  D.  2012. Ekonometrika  & Analisis Runtun Waktu Terapan dengan Eviews. Andi : Yogyakarta [12]  Soejoeti,  Z.  1987. Materi  Pokok Analisis  Runtun  Waktu .  Jakarta.: Karunika
[13]  Wei,  W.  W.  S.  2006. Time  Series Analysis,  Univariate  and  Multivariate Methods .  Canada.:  Addison  Wesley Publishing Company
[14]  Kementerian  Energi  dan  Sumber Daya  Mineral  Republik  Indonesia. 2017. http://www.esdm.go.id diakses pada  17  Oktober  2017  pukul  20.00 WIB.