PERAMALAN KECEPATAN ANGIN RATA-RATA HARIAN DI SURABAYA MENGGUNAKAN METODE BAYESIAN MODEL AVERAGING DENGAN PENDEKATAN EXPECTATION MAXIMIZATION

PERAMALAN KECEPATAN ANGIN RATA-RATA HARIAN DI SURABAYA MENGGUNAKAN METODE BAYESIAN MODEL AVERAGING DENGAN PENDEKATAN EXPECTATION MAXIMIZATION Nama NRP Jurusan Dosen Pembimbing

: Diah Kusumawati : 1307 100 049 : Statistika FMIPA-ITS : Dr. Rer. Pol. Heri Kuswanto, S.Si, M.Si

Abstrak Perubahan iklim yang terjadi dewasa ini menjadikan cuaca sangat sulit untuk diprediksi. Salah satu variabel cuaca yang terkena dampak dari perubahan iklim ini adalah kecepatan angin. Kecepatan angin sendiri selain memiliki beberapa fungsi penting, namun jika mencapai kisaran tertentu juga dapat memberikan dampak merugikan bagi kehidupan manusia. Dalam penelitian ini akan diperkenalkan suatu metode baru untuk mengkalibrasi permalan kecepatan angin, yaitu metode BMA-EM. Data peramalannya sendiri merupakan data ensembel tiruan yang dibangkitkan dari model ARIMA. Sebelum dikalibrasi, diketahui bahwa peramalan ensemble hanya dapat menangkap sebagian kecil observasi (17% untuk Lead ke-1 dan 12% untuk Lead ke-7). Sedangkan jika dibandingkan dengan hasil kalibrasi dengan metode BMA-EM, presentase ini mengalami peningkatan yang cukup besar, yaitu mencapai 92% untuk Lead ke-1 dan 76% untuk Lead ke-7. Berdasarkan nilai CRPS terkecil, dapat diketahui bahwa metode BMA-EM ini dapat menghasilkan peramalan yang lebih akurat dan reliabel. Kata kunci : Kecepatan angin, peramalan ensemble, BMA-EM. 1.

PENDAHULUAN Perubahan iklim global yang terjadi dewasa ini cukup memberikan dampak yang mengkhawatirkan bagi kelangsungan makhluk hidup di seluruh dunia, termasuk juga bagi Indonesia. LAPAN (2009) menyatakan bahwa dampak ekstrem dari perubahan iklim terutama adalah terjadinya kenaikan temperatur serta pergeseran musim. Selain itu, perubahan iklim yang terjadi ini dapat memberikan dampak pada berubahnya suhu udara, curah hujan, tekanan udara, dan kecepatan angin. Dari keempat faktor ini, kecepatan angin merupakan faktor yang memberikan pengaruh penting dalam kehidupan sehari-hari manusia. Antara lain sebagai input pembangkit listrik tenaga angin, faktor penentu penerbangan dan pelayaran, pengaruh kecepatan air dalam pengairan tanaman, serta merupakan faktor penting dalam mempengaruhi prakiraan cuaca. Sebaliknya, kecepatan angin dalam kisaran tertentu (melebihi batas maksimum kondisi aman) juga dapat menimbulkan kerugian bagi manusia. Pentingnya meramalkan kecepatan angin untuk menghindari timbulnya kerugian bagi kehidupan manusia inilah yang menjadi latar belakang dari adanya penelitian ini, terutama bagi kota besar seperti kota Surabaya. Beberapa penelitian mengenai peramalan kecepatan angin telah dilakukan sebelumnya, antara lain oleh Irhamah, dkk (2010), dan Faulina (2010). Perbedaan yang mendasar pada penelitian sebelumnya dengan penelitian ini adalah pada metode taksiran yang digunakan. Pada penelitian sebelumnya, peramalan yang dihasilkan berupa taksiran titik. Di banyak pusat prediksi di Negara-negara maju seperti NAEFS, telah dikembangkan model peramalan dengan ensemble (Ensemble Prediction System, EPS) yang akan memberikan luaran berupa taksiran interval. Zhu (2005) juga menyatakan, kelebihan dari penggunaan model ensemble dibandingkan peramalan yang menghasilkan taksiran titik adalah dapat menangkap adanya unsur / faktor ketidakpastian. Dalam prakteknya metode ensemble yang digunakan pada EPS dibangkitkan dari suatu proses numerik yang sangat rumit dan dengan tingkat komputasi yang tinggi, selain itu keterbatasan data juga merupakan faktor penghambat dalam peramalan dengan menggunakan ensemble ini. Oleh karena itu, dalam penelitian ini pembangkitan data ensemble dilakukan dengan alternatif lain, yaitu dengan memanfaatkan empat model time series berbeda. Raftery, Gneiting, Balabdaoui, dan Polakowski (2005) menyatakan salah satu karakteristik dari peramalan ensemble adalah sering kali bersifat underdispersive, yaitu hasil peramalan yang dihasilkan cenderung terpusat pada suatu titik tertentu dengan nilai varians yang relatif rendah. Oleh karena itu, dalam penelitian ini akan dilakukan suatu proses kalibrasi dengan menggunakan metode BMA-EM, yang dimaksudkan untuk menggeser atau melakukan penyesuaian terhadap nilai mean dan varian dari hasil peramalan sedemikian rupa sehingga dapat mendekati nilai observasi. Dengan penggunaan metode ini diharapkan dapat menghasilkan peramalan yang lebih akurat dan reliabel bila dibandingkan dengan menggunakan metode peramalan ensemble tiruan dengan 1

memanfaatkan beberapa model time series pada umumnya. Sehingga dapat menjadi alternatif ataupun metode baru dalam melakukan peramalan cuaca terutama untuk kecepatan angin. 2.

TINJAUAN PUSTAKA

Model Time series Deret waktu (time series) adalah rangkaian pengamatan yang diamati secara berurutan menurut urutan waktu. Beberapa model time series yang umum digunakan antara lain Model AR (q) dengan persamaan ∅ = + ; Model MA (q) dengan persamaan = + ; Model ARMA (p,q) dengan persamaan ∅ = + ; serta Model ARIMA ARIMA (p,d,q) dengan persamaan ∅ 1 − = + . Adapun tahapan-tahapan prosedur Box-Jenkins dalam Wei (2006) untuk pembentuk model time series antara lain identifikasi model, estimasi parameter dan pengujian signifikansi parameter model time series, diagnostic checking, kriteria pemilihan model, serta peramalan. Namun karena dalam penelitian ini data peramalan ensembel tiruan dibangkitkan dengan 4 model yang berbeda, maka kriteria pemilihan model tidak perlu dilakukan. Selain itu, tahapan peramalan yang akan dilakukan dari 4 macam model time series dalam penelitian ini bukanlah peramalan sekaligus dalam banyak periode, melainkan dilakukan dengan memperbarui model peramalan dari hari ke hari. Peramalan ini umumnya disebut dengan peramalan ensembel tiruan. Bayesian Model Averaging (BMA) Bayesian Model Averaging pertama kali diperkenalkan oleh Raftery, dkk pada tahun 2005, dimana metode ini merupakan salah satu metode kalibrasi yang akan menghasilkan fungsi densitas (pdf) hasil prediksi yang lebih akurat. Kalibrasi ini perlu dilakukan dikarenakan peramalan ensembel umumnya memiliki sifat underdispersive (terpusat pada suatu titik tertentu dengan nilai varians yang relatif rendah) ataupun overdispersive (varians tinggi). Prosedur kerja dari BMA pada dasarnya adalah memberikan bobot yang sesuai terhadap peramalan ensembel, dimana bobot yang diberikan didasarkan pada kemampuan prediksi dari masing-masing model ensembel yang digunakan. Misalkan = , , … , merupakan prediksi ensembel yang diperoleh dari K model yang berbeda, dan y merupakan peramalan ensembel yang terkalibrasi. Dimana dalam BMA, setiap peramalan ensembel dari , (k = 1, 2, …, K) terkait dengan suatu pdf bersyarat | yang dapat diartikan sebagai pdf bersyarat dari y pada , dimana merupakan peramalan terbaik dalam ensembel. Model prediksi BMA untuk peramalan ensembel dinyatakan sebagai | , , … , = ∑ | Dimana menunjukkan probabilitas posterior bahwa ramalan ke k adalah ramalan terbaik yang selalu bernilai positif dan jika dijumlahkan akan sama dengan satu. dapat dipandang sebagai bobot yang mencerminkan kontribusi ataupun kemampuan prediksi dari masing-masing model selama periode training (Vrugt, dkk. 2008). Berdasarkan penelitian sebelumnya (oleh Sloughter, dkk pada 2010), dapat diketahui bahwa data peramalan ensemble tiruan dari kecepatan angin mengikuti distribusi Gamma dengan parameter bentuk α dan parameter skala β. Sehingga jika pada prosedur BMA merupakan pdf bersyarat dari y pada , dimana merupakan peramalan terbaik dalam ensembel, maka bentuk pdf persamaan adalah sebagai berikut ) | = "! $!% exp − !

#$!

!

Parameter dari distribusi gamma tergantung pada ensembel ramalan anggota fk, Melalui hubungan, * = + + + dan , = - + - dimana * = . / adalah mean dari distribusi dan , = 0. / adalah standar deviasi (Sloughter, dkk. 2010). Expectation Maximization (EM) Algoritma EM merupakan suatu prosedur iteratif dari estimasi maksimum likelihood. Untuk mengimplementasikan algoritma EM pada metode BMA, digunakan suatu kuantitas ℎ yang tidak teramati, dimana akan mempunyai nilai 1 jika anggota ensembel k merupakan prediksi terbaik pada waktu ke t dan bernilai 0 jika sebaliknya. Pertama kali perlu ditentukan nilai awal dari bobot dan varians untuk kemudian algoritma EM akan berjalan secara iteratif antara expectation dan maximization sampai diperoleh kondisi yang konvergen. Pada step expectation, nilai ℎ diestimasi kembali dengan diberikan nilai bobot dan varian terbaru mengikuti persamaan berikut 2

3 ℎ2 =

4! 5)6 |7!6 ,89:;

A:; ∑! |?@> ,σ

Dimana j merupakan jumlah iterasi dan | , , 3% adalah pdf bersyarat dari anggota k yang berpusat pada . Pada step maximization nilai dari bobot dan varian diperbarui menggunakan nilai 3 estimasi dari ℎ yaitu ℎ2 seperti pada persamaan berikut 3 3 2 3 − = ∑C ℎ2 ; , 3 = ∑C ∑ ℎ

C

C

Dengan iterasi antara expectation dan maximization, algoritma EM memperbarui nilai wk dan σ2. Iterasi berhenti sampai nilainya konvergen dengan toleransi yang sangat kecil (Vrugt, dkk. 2008).

Evaluasi Kebaikan Model dengan Continous Ranked Probability Score (CRPS) Hasil peramalan dari kalibrasi ensembel berupa sebuah taksiran interval dalam bentuk pdf. Oleh karena itu, evaluasi kebaikan model tidak bisa dilakukan dengan prosedur standar seperti MSE atau MAPE. Sehingga dalam penelitian ini akan digunakan evaluasi kebaikan model dengan prosedur CRPS. CRPS berkaitan dengan probabilitas rank / pangkat, yang membandingkan distribusi hasil ramalan dengan pengamatan, dimana keduanya direpresentasikan sebagai fungsi distribusi kumulatif (cdf).

S∞

7

DEFG = CHIJK ∑CHIJK O LS%∞ MNO P − NO PQ RP

7 NO P : cdf dari hasil peramalan (ke – i) NO P : cdf dari pengamatan sebenarnya (ke – i) T- UV : jumlah peramalan Hasil dari CRPS ini berupa suatu nilai, dimana peramalan akan dikatakan bagus bila CRPS yang dihasilkan sangat kecil atau mendekati nol.

Kecepatan Angin Angin adalah udara yang bergerak yang diakibatkan oleh rotasi bumi dan juga karena adanya perbedaan tekanan udara di sekitarnya. Angin bergerak dari tempat bertekanan udara tinggi ke tempat yang bertekanan udara rendah. Sifat angin meliputi kekuatan angin, arah angin, dan kecepatan angin. Kekuatan dan arah angin dapat diketahui dengan bermacam cara, antara lain dengan bendera angin. Sedangkan pengukuran kecepatan angin dapat dilakukan dengan alat anemometer, yang diletakkan pada ketinggian seratus meter diatas permukaan air. Wilayah di Indonesia sendiri memiliki karakteristik kecepatan angin rata-rata yang berada pada skala rendah, yaitu antara 3 m/s hingga 6 m/s (Nehemia, 2009). Kecepatan angin juga merupakan parameter meteorologi yang memiliki fungsi penting, diantaranya untuk mendeteksi polusi udara terutama di daerah perkotaan, prakiraan cuaca, kecepatan aliran air untuk pengairan tanaman, take off dan landing pesawat, pelayaran, dan lain-lain 3.

METODOLOGI PENELITIAN Data yang digunakan dalam penelitian ini adalah data kecepatan angin rata-rata perhari (knot) tahun 2008 hingga tahun 2009 pada stasiun pengamatan Juanda Surabaya. Data diperoleh dari Badan Meteorologi Klimatologi dan Geofisika Indonesia. Data kecepatan angin yang digunakan dalam penelitian ini terlebih dahulu akan dibagi menjadi dua, yaitu data in-sample dan data out-sample. Untuk pembentukan model kalibrasi dan peramalan (in-sample), data yang digunakan adalah data pada bulan Januari 2008 sampai dengan September 2009. Sedangkan untuk validasi hasil peramalan atau pemilihan peramalan terbaik (out-sample) digunakan data 3 bulan berikutnya, yaitu bulan Oktober 2009 sampai Desember 2009. Langkah-langkah analisis yang dilakukan pada penelitian ini sebagai berikut. 1. Membagi data menjadi 2 bagian, yaitu data in-sample (Januari 2008 – September 2009) dan data outsample (Oktober 2009 – Desember 2009). 2. Mengidentifikasi data series melalui plot Time series, ACF, dan PACF-nya (terhadap data in-sample). Jika data telah stasioner maka dapat menuju ke proses selanjutnya. Sebaliknya, jika data belum stasioner terhadap mean dilakukan differencing dan jika data belum stasioner terhadap varians dilakukan transformasi. 3. Mengidentifikasi model time series berdasarkan plot ACF dan PACF. Dari identifikasi ini akan dicari empat macam model time series yang sesuai untuk data kecepatan angin di 4. Membangkitkan data ensembel tiruan dari keempat model time series. Proses ini dilakukan dengan memperbarui model peramalan dari hari ke hari, sehingga bukan peramalan sekaligus dalam banyak periode seperti pada umumnya, melainkan peramalan satu tahap kedepan. 3

5. Kalibrasi dan peramalan dengan BMA menggunakan pendekatan EM. Melalui langkah kalibrasi ini akan didapatkan parameter-parameter kalibrasi yang akan digunakan dalam mendapatkan taksiran interval dengan BMA menggunakan pendekatan EM. 6. Evaluasi ketepatan peramalan antara kedua metode yang digunakan (ensembel tiruan dan BMA – EM) dengan menggunakan CRPS. 4.

ANALISIS DAN PEMBAHASAN Kecepatan angin rata-rata harian di Surabaya untuk data in-sample memiliki kecepatan rata-rata 8,594 km/jam, dengan kecepatan angin terendah 0,5 km/jam yang terjadi pada 22 April 2009 dan kecepatan angin tertinggi 26,6 km/jam terjadi pada 22 Pebruari 2008. Nilai skewness yang positif (1,13) menunjukkan bahwa pola penyebaran/ distribusi data kecepatan angin yang digunakan tidak simetris serta sebagian besar dari data yang digunakan cenderung berada di bawah nilai rata-rata. Selain itu, nilai kurtosis yang positif (2,43) menunjukkan distribusi data yang memuncak/ runcing. 30

Lower CL

Upper CL Lambda

8

(using 95.0% confidence)

25

Estimate

7 6 StDev

wind

20

15

4

5

3 2 1

64

128

192

256

320 Index

384

448

512

0.35 0.69

Rounded Value

0.50

5

10

0

0.52

Lower CL Upper CL

576

Limit -2

-1

0

1 2 Lambda

3

4

5

Gambar 1 Plot Time series dan Box-Cox Data In-sample Kecepatan Angin Rata-rata Harian di Surabaya

Berdasarkan Gambar 1(a), secara visual dapat diketahui bahwa variansi antara titik-titik pengamatan cukup besar, hal ini ditunjukkan pada titik-titik tertentu terdapat kecepatan angin yang cenderung terletak jauh dari nilai rata-rata (8,594). Sedangkan melalui Gambar 1(b), dapat diketahui bahwa dengan menggunakan tingkat kepercayaan 95% diperoleh nilai optimal lambda sebesar 0,5. Nilai lambda yang kurang dari 1 ini menunjukkan bahwa data yang digunakan belum stasioner terhadap varian sehingga perlu dilakukan transformasi akar, dan data hasil dari transformasi akar inilah yang akan digunakan dalam analisis selanjutnya (with 5% significance limits for the partial autocorrelations) 1.0

0.8

0.8

0.6

0.6

Partial Autocorrelation

Autocorrelation

(with 5% significance limits for the autocorrelations) 1.0

0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8

0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8

-1.0

-1.0 1

5

10

15

20

25

30

35 40 Lag

45

50

55

60

65

70

1

5

10

15

20

25

30

35 40 Lag

45

50

55

60

65

70

Gambar 2 Plot ACF dan PACF Data Kecepatan Angin (Hasil Transformasi)

Setelah diperoleh data pengamatan baru (data hasil transformasi akar) yang telah stasioner dalam varian, langkah selanjutnya adalah mengidentifikasi stasioneritas dalam mean terhadap data hasil transformasi tersebut. Plot ACF (Gambar 2(a)) dari data hasil transformasi menunjukkan pola turun secara eksponensial/ cepat pada setiap lag-lagnya. Sedangkan pada plot PACF (Gambar 2(b)), lag pertama hingga lag ketiga keluar dari batas signifikansi. Hal ini mengindikasikan bahwa data yang digunakan telah stasioner terhadap mean.

4

Estimasi Parameter

Model Constant P-value ϕ1 P-value ϕ2 P-value ϕ3 P-value θ3 P-value θ7 P-value θ12 P-value θ19 P-value White Noise Residual Normal

Tabel 1. Evaluasi 3 Model ARIMA Yang Digunakan (3,0,(19)) (3,0,(12,19)) (3,0,(12)) (3,0,(7,19))

AR (3) 2.8926 0,000 0.43878 0,000 0.17791 0,000 0.14491 0,000 -

2.8899 0,000 0.44336 0,000 0.17868 0,000 0.14889 0,000 -

2.88957 0,000 0.43331 0,000 0.18063 0,000 0.14579 0,000 -

-

-

-

-

-

-

0.09878 0.0154 Tidak Memenuhi

-0.08957 0.0289 0.09844 0.0157 Memenuhi Memenuhi

Tidak Memenuhi

(3,0,(3,7,19))

(3,0,(3,7))

0.07597 0.0635* -

2.89659 0,000 0.4306 0,000 0.16088 0,000 0.27491 0,000 0.15859 0.0112 0.10526 0.014 -

2.9 0,000 0.42706 0,000 0.16123 0,000 0.25965 0,000 0.14402 0.0241 0.09415 0.0277 -

0.10371 0.0111 Memenuhi Memenuhi

0.10755 0.0076 Memenuhi Memenuhi

2.89286 0,000 0.42964 0,000 0.17994 0,000 0.14272 0,000 -

2.88982 0,000 0.44725 0,000 0.18654 0,000 0.15823 0,000 -

-0.08829 0.0316 Memenuhi Memenuhi

Memenuhi Memenuhi

Dalam penelitian ini tidak hanya digunakan satu model terbaik saja, namun akan digunakan empat model ARIMA yang berbeda untuk membangkitkan data ensemble. Tabel 1 menunjukkan hasil evaluasi keempat model ARIMA yang akan digunakan dalam penelitian (yaitu ARIMA (3,0,(12,19)), ARIMA (3,0,(12)), ARIMA (3,0,(3,7,19)) dan ARIMA (3,0,(3,7))), dimana dapat diketahui bahwa estimasi perameter untuk setiap model telah signifikan, serta syarat residual white noise dan berdistribusi normal telah terpenuhi. Keempat model ARIMA yang telah teridentifikasi ini selanjutnya akan digunakan sebagai model referensi untuk membangkitkan data ensemble tiruan pada periode-periode berikutnya. • Peramalan Kecepatan Angin Dengan Menggunakan Ensemble Tiruan Peramalan ensemble tiruan yang dilakukan dalam penelitian ini bukanlah peramalan sekaligus dalam banyak periode, melainkan dengan memperbarui model peramalan dari hari ke hari. Peramalan ensemble tiruan dari keempat model berbeda ini akan dilakukan dengan dua lead time, yaitu peramalan untuk satu tahap/ hari kedepan (pada lead ke-1) serta peramalan untuk tujuh tahap/ hari kedepan (pada lead ke-7). Peramalan Dengan Ensemble Tiruan Untuk Lead ke-1 Untuk peramalan 1 Oktober 2009 (Lead ke-1), diperoleh dengan memodelkan data in-sample (1 Januari 2008 – 30 September 2009) mengikuti model ARIMA (3,0,(12,19)). Dengan memperbarui data in-sample (mengurangi satu data terlama dan menambahkan satu data terbaru) menjadi 2 Januari 2008 – 1 Oktober 2009 dan dengan model ARIMA (3,0,(12,19)), akan diperoleh hasil peramalan untuk 2 Oktober 2009. Pembaruan data in-sample ini dilakukan terus-menerus hingga didapatkan peramalan untuk 31 Desember 2009. Selain itu, langkah ini juga dilakukan terhadap 3 model lainnya (ARIMA (3,0,(12)), ARIMA (3,0,(3,7,19)) dan ARIMA (3,0,(3,7))). Dari hasil peramalan ini kemudian dikembalikan ke bentuk awal dari transformasi (dikuadratkan). Sehingga hasil untuk peramalan ensemble tiruan (lead ke1) tersaji pada Gambar 3 berikut ini. 2 2 .5

V a r iab le O bs M .1 M .2 M .3 M .4

2 0 .0 1 7 .5 1 5 .0 1 2 .5 1 0 .0 7 .5 5 .0

1

9

18

27

36

45

54

63

72

81

Gambar 3 Hasil Peramalan Ensemble Tiruan (Lead ke-1)

5

90

Setelah didapatkan 4 buah peramalan ensemble untuk masing-masing model, dapat diperoleh nilai µ dan σ2 (perhari) dari keempat model peramalan ensemble tiruan tersebut. Berdasarkan penelitian sebelumnya oleh Sloughter dkk (2010), data dari masing-masing ensemble mengikuti distribusi gamma. Sehingga parameter gamma (α dan β) akan diestimasi dari nilai µ dan σ2 yang telah ditemukan sebelumnya untuk kemudian dibentuk menjadi sebuah pdf, seperti pada Gambar 4 berikut. Gamma, Shape=1253.38, Scale=0.006745

Gamma, Shape=351.278, Scale=0.024769


8.7

10 0.8

1.4

0.7

1.2

0.6

1.0 0.8

0.4 0.3

0.4

0.2

0.0

0.025

0.0

8.928

0.025

0.025 7.815

10/01/09

0.0

9.634

9.814 11/01/09

Gamma, Shape=308.25, Scale=0.0235197 8.1


7.4

1.6

5.5 1.8

1.0

1.6

1.4 0.8

1.2

1.4 1.2

0.8 0.6

0.6

Density

Density

1.0 Density

0.025 9.077

10/31/09


0.4

1.0 0.8 0.6

0.4

0.4

0.2 0.2 0.025 0.0

1.0

0.5

0.1

0.025 7.992

1.5

0.5

0.6

0.025

2.0

Density

1.6

0.2

12.4

0.9

Density

Density

1.8

0.025

0.025 5.906

0.0

6.952 11/30/09

0.025 6.463

0.2

0.025

0.0

8.081

0.025 5.437

6.418

12/01/09

12/31/09

Gambar 4 PDF Kecepatan Angin untuk Lead ke-1

Berdasarkan Gambar 4 dapat diketahui bahwa sebagian besar nilai observasi (garis vertikal) berada diluar selang kepercayaan 95% dari hasil peramalan ensemble tiruan. Selain itu, dapat diketahui bahwa dari total peramalan sebanyak 92 hari (1 Oktober 2009 – 31 Desember 2009), hanya terdapat 16 hari dimana nilai observasi masuk ke dalam interval peramalan. Hal ini mengindikasikan bahwa hasil peramalan ensemble tiruan umunya bersifat underdispersive, yaitu memiliki varians yang relatif kecil sehingga akan sulit untuk menangkap nilai observasi. Peramalan Dengan Ensemble Tiruan Untuk Lead ke-7 Langkah yang harus dilakukan dalam peramalan ensemble tiruan untuk 7 hari kedepan ini hampir sama dengan peramalan ensemble tiruan untuk 1 tahap kedepan. Yang menjadi perbedaan adalah data insample digunakan untuk meramalkan tujuh hari kedepan, namun hanya mengambil hasil peramalan pada lead ke-7. Sehingga untuk data in-sample awal (1 Januari 2008 – 30 September 2009) digunakan untuk meramalkan tanggal 7 Oktober 2009. Dan data in-sampel 2 Januari 2008 – 1 Oktober 2009 digunakan untuk meramalan 8 Oktober 2009. Langkah ini dilakukan terus menerus hingga didapatkan hasil peramalan ensemble untuk 31 Desember 2009, untuk masing-masing model yang digunakan seperti pada Gambar 5. 2 2 .5

V a r ia b le O bs M .1 M .2 M .3 M .4

2 0 .0 1 7 .5 1 5 .0 1 2 .5 1 0 .0 7 .5 5 .0

1

9

18

27

36

45

54

63

72

81

Gambar 5 Hasil Peramalan Ensemble Tiruan (Lead ke-7)

Berdasarkan Gambar 5, terlihat bahwa hasil peramalan ensemble tiruan untuk lead ke-7 terletak jauh dari nilai observasi perharinya. Hal ini membuktikan bahwa lead time peramalan yang panjang akan menghasilkan peramalan dengan tingkat kesalahan dan ketidakpastian yang lebih tinggi (jauh berbeda dengan hasil peramalan ensemble untuk lead ke-1). Seperti pada peramalan untuk lead ke-1, setelah mendapatkan hasil peramalan ensemble serta nilai µ dan σ2 dari peramalan perharinya, langkah berikutnya 6

adalah membentuk fungsi distribusi probabilitas mengikuti distribusi gamma dengan parameter α dan β, seperti pada Gambar 6 berikut ini. Gamma, Shape=2830.71, Scale=0.00292444



8.9

10

12.4

0.8 4

2.5 0.7 2.0

0.6

3

1.0

Density

Density

Density

0.5 1.5

0.4

2

0.3 0.2

1

0.5 0.1 0.025 0.0

0.025

0.025 7.976

0.0

8.586

0.025

0.025 7.626

10/07/09

0

9.659

Gamma, Shape=905.312, Scale=0.0086703 8.1


7.4

5.5

2.0

1.8 1.6

1.4 1.2

1.4

1.5

1.2

0.8

Density

Density

1.0 Density

0.025 9.574 11/01/09


1.6

1.0

0.6

1.0 0.8 0.6

0.4

0.5

0.4

0.2 0.025 0.0

9.192

10/31/09

0.025 7.346

8.369 11/30/09

0.025 0.0

0.025 7.644

8.491

0.2 0.0

12/01/09

0.025 6.795 12/31/09

0.025 7.771

Gambar 6 PDF Kecepatan Angin untuk Lead ke-7

Gambar 6, menunjukkan hasil peramalan untuk lead ke-7 ini tidak jauh berbeda dengan hasil peramalan ensemble dengan menggunakan lead ke-1 (Gambar 4), dimana sebagian besar nilai observasi berada diluar selang kepercayaan 95% dari hasil peramalan ensemble. Untuk hasil peramalan ensemble dengan mengunakan lead ke-7 ini, dari total 86 hari peramalan (7 Oktober 2009 – 31 Desember 2009) hanya terdapat 10 interval peramalan ensemble yang dapat memuat nilai observasi. Hasil peramalan ini menunjukkan bahwa peramalan ensemble memiliki sifat underdispersive sehingga sulit menangkap nilai observasi. Jika dibandingkan dengan hasil peramalan ensemble dengan menggunakan lead ke-1, peramalan dengan lead ke-7 (peramalan jangka panjang) akan menghasilkan kesalahan peramalan yang lebih besar daripada menggunakan lead ke-1. • Kalibrasi Peramalan Kecepatan Angin Menggunakan Metode Bayesian Model Averaging Karena peramalan untuk 1 Oktober – 31 Desember 2009 dengan menggunakan metode ensemble tiruan (baik untuk lead ke-1 dan lead ke-7) bersifat underdispersive dan masih memberikan hasil yang kurang baik, maka hasil peramalan ensemble tiruan perlu dikalibrasi. Metode kalibrasi yang digunakan dalam penelitian ini adalah metode Bayesian Model Averaging dimana parameter-parameternya (varians dan bobot) akan diestimasi dengan pendekatan Expectation Maximization. Selain itu, dalam penelitian ini akan digunakan beragam nilai m (window training), antara lain 10, 15, 20, dan 25. Peramalan Dengan BMA-EM Untuk Lead ke-1 Penggunaan nilai m yang berbeda-beda akan menyebabkan perbedaan hari/ tanggal pertama yang dapat diramalkan. Untuk peramalan menggunakan Lead ke-1, hari/ tanggal pertama yang dapat diramalkan dimulai dari satu hari setelah m. Selain itu, penggunaan nilai m yang berbeda-beda juga akan menyebabkan hasil peramalan yang berbeda pula (meskipun dalam tanggal yang sama). Setelah menentukan ukuran m yang digunakan, langkah berikutnya adalah mendapatkan nilai +

dan + , σ2 dan weight untuk setiap model dalam peramalannya. Pada dasarnya nilai + dan + digunakan untuk menghilangkan bias dari peramalan ensemble, dimana proses ini bertujuan untuk mengeser nilai µ peramalan ensemble tiruan agar dapat mendekati nilai observasi melalui prosedur regresi. Nilai + dan + dapat diperoleh dengan cara meregresikan peramalan ensemble terhadap observasi dengan jumlah data yang digunakan sebanyak m data sebelum tanggal peramalan. Sebagai ilustrasi, misalkan akan dikalibrasi peramalan ensemble tanggal 31 Desember 2009 dengan m = 10, maka data yang digunakan adalah data mulai 21 Desember 2009 - 30 Desember 2009. Melalui proses regresi ini akan diperoleh koefisien regresi + = 16,31995 dan + = -0,32153. Koreksi bias didapatkan melalui rumusan * = + + + , dimana merupakan peramalan ensemble untuk tiap model pada 31 Desember 2009, maka akan diperoleh * = 7,028284; * = 6,78842; *W = 7,698124; dan *X = 7,076041. Nilai varian dan weight diestimasi dengan pendekatan EM. Dari pendekatan EM ini, akan diperoleh koefisien - = 6,94E-10 dan - = 1,5E-01, dengan standar deviasi , = - + - . Selanjutnya didapatkan , = 0,8974; , = 0,9205; ,W = 0,8327; dan ,X = 0,8927. Serta nilai bobot untuk masing7

masing model adalah = 0,284662; = 0,014391; W = 0,678097; dan X = 0,02285. Semakin besar bobot dari model itu, semakin tinggi tingkat kontribusi peramalan ensemble terhadap hasil peramalan terkalibrasi. Karena peramalan ensemble untuk kecepatan angin mengikuti distribusi gamma, maka parameter gamma diestimasi dari nilai * dan , yang telah ditemukan sebelumnya. Dan diperoleh . = 61.34391, / = 0.114572; . = 54.38414, / = 0.1248424; .W = 85.47348, /W = 0.090064; dan .X = 62.82457, /X = 0.112632. Sehingga dapat diperoleh model prediksi BMA untuk peramalan ensemble sebagai berikut. | , , W , X = ∑X | = | + | + W W |W + X X |X = . , / + . , / + W W .W , /W + X X .X , /X 0.45

0.4

a

0.4 0.35

data1 data2 data3 data4 data5

b

0.35 0.3 Probability density

Probability density

0.3 0.25 0.2 0.15

0.2 0.15 0.1

0.1

0.05

0.05 0

0.25

3

4

5

6

7 8 Kecepatan Angin

9

10

11

0

12

0.35

1

3

4

5 6 7 Kecepatan Angin

4

6 8 Kecepatan Angin

8

9

10

11

d 0.3

0.25

0.25 Probability density

Probability density

c 0.3

0.2

0.15

0.2

0.15

0.1

0.1

0.05

0.05

0

2

0.35

0

2

4

6 Kecepatan Angin

8

10

12

0

0

2

10

12

14

Gambar 7 PDF Kecepatan Angin Dengan Metode BMA-EM untuk 31 Desember 2009 (Lead ke-1), dimana (a) menggunakan m = 10, (b) m = 15, (c) m = 20, dan (d) m = 25.

Gambar 7 menunjukkan hasil peramalan metode BMA-EM untuk 31 Desember 2009 dengan menggunakan 4 macam nilai m yang berbeda. Kurva biru dalam gambar menunjukkan pdf predictive untuk Model 1, kurva ungu untuk Model 2, kurva hijau untuk Model 3, kurva merah untuk peramalan ensemble Model 4 dan kurva hitam adalah hasil peramalan terkalibrasi atau predictive pdf gabungan (mixture). Sedangkan untuk garis vertikal berwarna biru merupakan data observasi dan garis vertikal berwarna merah merupakan selang kepercayaan 95% dari hasil peramalan terkalibrasi. Melalui Gambar 7 dapat diketahui bahwa dari empat ukuran m berbeda, hanya pada m=10 yang hasil peramalan terkalibrasinya (dengan selang kepercayaan 95%) tidak dapat menangkap nilai observasi yaitu 5,5. Pada Gambar 7(b) dapat dilihat bahwa hasil peramalan terkalibrasi berhimpit dengan predictive pdf Model 2. Hal ini menunjukkan bahwa peramalan ensemble tiruan untuk Model 2 adalah peramalan yang berkontribusi paling besar terhadap hasil peramalan terkalibrasi dengan kata lain peramalan ensemble Model 2 memiliki weight yang paling besar (Tabel 2). Tabel 2. Parameter Distribusi untuk 31 Desember 2009 (Lead ke-1) m=10 m=15 m=20 m=25 M1 0.28466 9.31E-05 1.58E-03 0.1051 M2 0.01439 0.9999068 9.98E-01 0.0039 w M3 0.67810 1.94E-11 8.77E-17 0.891 M4 0.02285 1.11E-07 1.16E-15 7E-05 µ -terkalibrasi 7.4820 6.1565 6.2743 6.1931 σ2 -terkalibrasi 0.8304 1.3000 1.8099 2.2814 batas bawah 5.6959 3.9218 3.6375 3.2327 batas atas 9.2681 8.3913 8.9111 9.1536

Berdasarkan Gambar 7 dan Tabel 2, dapat diketahui bahwa hasil peramalan terkalibrasi terbaik untuk 31 Desember 2009 adalah menggunakan m=15, dikarenakan hasil peramalan terkalibrasi memiliki bias paling kecil, sehingga nilai meannya mendekati observasi (5,5) serta menghasilkan interval 8

peramalannya yang tidak terlalu lebar. Berikut ini akan disajikan hasil peramalan metode BMA (optimum) untuk hari lainnya. data1 data2 data3 data4 data5 Probability density

Probability density

0.3 0.25 0.2 0.15 0.1

0.18

0.16

0.16

0.14

0.14

0.12

0.12

0.1 0.08 0.06 0.04

0.05 0

0.18

Probability density

0.4 0.35

6

7

8


12

13

14

15

0

0.06 0.04

0.02

5

0.1 0.08

0.02

0

2

4

6

10/31/09 (m=15)


14

16

18

0

20

0

5

10 15 Kecepatan Angin

11/01/09 (m=25) 0.35

0.25

20

25

11/15/09 (m=25) 0.35

0.3

0.3

0.25

0.25

0.15

0.1

Probability density

Probability density

Probability density

0.2

0.2

0.15

0.2

0.15

0.1

0.1

0.05

0.05

0.05

0 -2

0

2

4

6 8 Kecepatan Angin

10

12

14

16

0

0

2

4

6 8 Kecepatan Angin

11/30/09 (m=20)

10

12

14

0

0

2

4

12/01/09 (m=15)

6 8 Kecepatan Angin

10

12

14

12/15/09 (m=20)

Gambar 8 PDF Kecepatan Angin Dengan Metode BMA-EM Untuk Beberapa Hari Lainnya (Lead ke-1)

Berdasarkan Gambar 8, dapat diketahui bahwa hanya terdapat sebagian kecil nilai observasi yang berada diluar batas selang kepercayaan 95% dari hasil peramalan ensemble terkalibrasi (lead ke-1). Sehingga bila dibandingkan dengan hasil peramalan ensemble tiruan (untuk lead ke-1), dapat diketahui bahwa hasil peramalan dengan menggunakan metode BMA-EM akan lebih akurat dan reliabel karena lebih bisa menangkap nilai observasi. Peramalan Dengan BMA-EM Untuk Lead ke-7 Untuk peramalan menggunakan Lead ke-7 ini, hari/ tanggal pertama yang dapat diramalkan adalah tujuh hari setelah m-hari. Langkah-langkah yang harus dilakukan untuk peramalan ensemble terkalibrasi untuk lead ke-7 ini sama dengan langkah untuk peramalan ensemble tiruan lead ke-1, yaitu menentukan ukuran m yang digunakan, mendapatkan koefisien b0 dan b1 melalui proses regresi dan mengestimasi nilai µk, serta mengestimasi σk2 dan wk dengan pendekatan EM. 0.45

0.45

a

0.4


Probability density

0.3 0.25 0.2 0.15

0.25 0.2 0.15

0.1

0.1

0.05

0.05

0

1

2

3

4

5 6 Kecepatan Angin

7

8

9

0

10

0.4

0

c

2

3

4 5 Kecepatan Angin

6

7

8

9

10

11

d

0.35 0.3 Probability density


1

0.4

0.35

0.25 0.2 0.15

0.25 0.2 0.15

0.1

0.1

0.05

0.05

0


b

0.4 0.35

0.35

1

2

3

4


8

9

10

11

0

1

2

3

4


8

9

Gambar 9 PDF Kecepatan Angin Dengan Metode BMA-EM untuk 31 Desember 2009 (Lead ke-7)

Pada Gambar 9 dapat diketahui bahwa hasil peramalan terkalibrasi (menggunakan setiap nilai m) dengan selang kepercayaan 95% telah dapat menangkap nilai observasi, yaitu 5,5. Dari Gambar 9(a), 9(c), dan 9(d) juga dapat dilihat bahwa weight untuk masing-masing model berbeda. Hal ini menunjukkan bahwa setiap peramalan ensemble memiliki kontribusi yang berbeda-beda terhadap hasil peralamalan 9

terkalibrasi. Sebaliknya, pada Gambar 9(b) hasil peramalan ensemble terkalibrasi berhimpit dengan Model 2, yang menunjukkan bahwa peramalan ensemble untuk Model 2 memberikan kontribusi yang paling besar terhadap hasil peramalan terkalibrasi dibanding model lainnya. Selain itu, berdasarkan Gambar 9 dapat diketahui bahwa hasil peramalan terbaik untuk 31 Desember 2009 adalah dengan menggunakan m=10, dilihat dari sisi bias dan interval peramalannya. Berikut adalah hasil peramalan metode BMA untuk hari lainnya. 0.2

0.35

0.14


0.18 0.3 0.16

0.2

0.15

0.12



Probability density

0.25

0.12 0.1 0.08 0.06

0.1

0.08

0.06

0.04

0.04 0.02

0.05 0.02 0

4

6

8

10 Kecepatan Angin

12

14

0

16

0

2

4

10/31/09 (m=15)

6


14

16

18

0 -5

20

0

5

11/01/09 (m=25)

0.25

0.25

0.2

0.2

10 15 Kecepatan Angin

20

25

30

11/15/09 (m=25) 0.35

0.3

0.15

0.1

Probability density

Probability density

Probability density

0.25 0.15

0.1

0.2

0.15

0.1

0.05

0.05

0.05

0 -2

0

2

4

6 8 Kecepatan Angin

10

12

14

16

0 -2

11/30/09 (m=20)

0 0

2

4

6 8 Kecepatan Angin

10

12

14

16

0

2

12/01/09 (m=15)

4

6 8 Kecepatan Angin

10

12

14

12/15/09 (m=10)

Gambar 10 PDF Kecepatan Angin Dengan Metode BMA-EM Untuk Beberapa Hari Lainnya (Lead ke-7)

Sama halnya dengan peramalan ensemble terkalibrasi menggunakan lead ke-1, melalui Gambar 10 dapat diketahui bahwa hanya terdapat sebagian kecil nilai observasi yang berada diluar batas selang kepercayaan 95% dari hasil peramalan ensemble terkalibrasi (lead ke-7). Sehingga bila dibandingkan dengan hasil peramalan ensemble tiruan (untuk lead ke-7), dapat diketahui bahwa hasil peramalan dengan menggunakan metode BMA-EM akan menghasilkan peramalan yang lebih baik (akurat dan reliabel) karena lebih bisa menangkap nilai observasi. Namun jika hasil peramalan ensemble terkalibrasi pada lead ke-7 ini dibandingkan dengan peramalan ensemble terkalibrasi menggunakan lead ke-1, maka hasil terbaik adalah dengan menggunakan peramalan untuk lead ke-1. Dikarenakan, semakin panjang lead time yang digunakan dalam peramalan maka unsur ketidakpastian semakin besar.

• Perbandingan Antara Kedua Metode Yang Digunakan Evaluasi hasil peramalan antara kedua metode yang digunakan ini dilakukan dengan menggunakan Continous Ranked Probability Score (CRPS), dimana prosedur dasarnya adalah membandingkan distribusi hasil ramalan dengan pengamatan, yang keduanya direpresentasikan sebagai fungsi distribusi kumulatif (cdf). Tabel 3 Evaluasi Hasil Peramalan Lead Ke-1 Ke-7

CRPS Ensemble Tiruan BMA – EM Ensemble Tiruan BMA – EM

m = 10 1.83889 1.44442 2.22787 1.51923

m = 15 1.67675 1.27134 2.25877 1.52070

m = 20 1.66256 1.32036 2.32267 1.73238

m = 25 1.67322 1.29918 2.44868 1.72472

Berdasarkan Tabel 3 dapat diketahui bahwa hasil peramalan dengan menggunakan metode BMAEM akan lebih mendekati nilai observasi daripada hasil peramalan dengan metode ensemble tiruan, selain itu hasil interval peramalannya juga padat. Hal ini ditunjukkan oleh nilai CRPS metode kalibrasi BMAEM yang lebih kecil daripada ensemble tiruan. Sedangkan untuk penggunaan m, dapat diketahui bahwa peramalan terbaik untuk Lead ke-1 adalah dengan menggunakan nilai m sebesar 15, dan untuk Lead ke-7 sebesar 10, dikarenakan akan menghasilkan nilai CRPS yang optimum/ terkecil dibandingkan dengan menggunakan nilai m lainnya. 10

Lead ke-7 (m=10) Lead ke-1 (m=15)

Lead ke-1

16 (17%) 76 (83%) Dapat menangkap obs.

Lead ke-7

Tidak dapat

10 (12%) 76 (88%)

0

20

40

60

80

71 (92%) 6 (8%) Dapat menangkap obs. Tidak dapat

58 (76%) 18 (24%)

0

20

40

60

80

Gambar 11. Evaluasi Hasil Peramalan

Gambar 11 menunjukkan bahwa penggunaan kalibrasi dengan metode BMA-EM akan meningkatkan keakuratan dari hasil peramalan. Hal ini dibuktikan dengan banyaknya observasi yang semula tidak dapat ditangkap ensemble tiruan, namun setelah dilakukan kalibrasi dengan metode BMAEM, sebagian besar observasi telah dapat ditangkap oleh hasil peramalan. 5.

KESIMPULAN Kesimpulan yang dapat diambil berdasarkan hasil analisis dan pembahasan adalah sebagai berikut. 1. Dengan menggunakan metode peramalan ensemble tiruan dapat diketahui bahwa hanya 17% (untuk Lead ke-1) dan 12% (untuk Lead ke-7) dari interval peramalan ensemble yang dapat memuat nilai observasi. Hal ini menunjukkan bahwa peramalan dengan menggunakan ensemble tiruan masih memberikan hasil yang kurang baik. 2. Penggunaan kalibrasi metode BMA-EM terhadap peramalan ensemble dengan training window optimum (m=15 untuk lead ke-1 dan m=10 untuk lead ke-7) menunjukkan bahwa terjadi peningkatan presentase peramalan yang dapat memuat nilai observasi, menjadi 92% untuk Lead ke-1 dan 76% untuk Lead ke-7. Hal ini membuktikan bahwa metode kalibrasi BMA-EM akan meningkatkan keakuratan dan reliabiltas dari hasil peramalan ensembel. 3. Melalui perhitungan CRPS (Continuous Ranked Probability Score), dapat diketahui bahwa metode kalibrasi dengan BMA-EM dapat meghasilkan peramalan yang lebih baik dari ensembel tiruan yang belum terkalibrasi. Karena metode BMA-EM akan menghasilkan nilai CRPS yang lebih kecil Pada penelitian ini hanya digunakan empat model ARIMA yang berbeda. Sehingga untuk penelitian selanjutnya disarankan pembangkitkan data ensemblenya dapat menggunakan jumlah model yang lebih banyak, mengingat semakin banyak model yang digunakan maka semakin besar faktor ketidakpastian yang dapat ditangkap oleh model tersebut. DAFTAR PUSTAKA Faulina, Ria. (2010). Adaptive Neuro-Fuzzy Inference System Untuk Peramalan Kecepatan Angin RataRata Harian Di Sumenep. Tugas Akhir. Mahasiswa Jurusan Statistika FMIPA ITS, Surabaya. Irhamah, Prasetyo, D.D., dan Fithriasari, K. (2010). Pengembangan Model Ramalan Kecepatan Angin Menggunakan Hybrid Time series dan Algoritma Genetika. Penelitian Produktif ITS – 2010. Surabaya : ITS. LAPAN. (2009). (http://iklim.dirgantaraDampak Perubahan Iklim. lapan.or.id/index.php?option=com_content&view=article&id=60&Itemid=37, diakses pada 2 Februari lapan.or.id/index.php?option=com_content&view=article&id=60&Itemid=37, 2011). Nehemia, Ronald. (2009). Optimalisasi Ekstraksi Energi Angin Kecepatan Rendah Di Indonesia Dengan Aplikasi Konverter Boost. (http://konversi.wordpress.com/2009/01/24/optimalisasi-ekstraksi-energiangin-kecepatan-rendah-di-indonesia-dengan-aplikasi-konverter-boost/, diakses pada 6 Februari 2011). Raftery, A. F., Gneiting, T., Balabdaoui, F., dan Polakowski, M. (2005). Using Using Bayesian Model to Calibrate Forecast Ensembles. Monthly Weather Review, Vol 133. Hal 1155-1173. Sloughter, M.J., Gneiting, T., dan Raftery, A.E. (2010). Probabilistic Wind Speed Forecasting Using Ensembles and Bayesian Model Averaging. Journal of the American Statistical Association, Vol 105, No 489. Hal 25-34. 11

Vrugt, J.A., Diks, C.G., dan Clark, M.P. (2008). Ensemble Bayesian Model Averaging Using Markov Chain Monte Carlo Sampling. Environmental Fluid Mechanics, Vol 8, No 5-6. Hal 579-595. Wei, W.W.S. (2006). Time series Analysis, Univariate and Multivariate Methods. Canada : Addison Wesley Publishing Company. Zhu, Yuejian. (2005). Ensemble Forecast: A New Approach To Uncertainty And Predictability. Advance In Atmosphere Science, Vol 22, No 6. Hal 781-788.

12

PERAMALAN KECEPATAN ANGIN RATA-RATA HARIAN DI SURABAYA MENGGUNAKAN METODE BAYESIAN MODEL AVERAGING DENGAN PENDEKATAN EXPECTATION MAXIMIZATION

Recommend Documents