Abstrak —Gempa  bumi merupakan peristiwa getaran yang disebabkan dari terjadinya pelepasan akumulasi energi dari lempeng tektonik di dalam bumi. Indonesia merupakan negara yang terletak pada pertemuan antara lempeng tektonik, yang menyebabkan Indonesia rawan terjadi gempa  bumi.  Di Indonesia gempa bumi sebesar 9.1 magnitude pernah terjadi dan menyebabkan tsunami yang memakan korban jiwa lebih dari 160 orang dan 63.977 keluarga kehilangan tempat tinggal. Terjadinya gempa bumi terutama gempa berkekuatan besar dapat menyebabkan besarnya dampak yang ditimbulkan. Risiko kerugian dapat diminimalisir dengan adanya mitigasi bencana sebelum terjadinya peristiwa bencana. Kerugian yang ditimbulkan dari bencana alam dapat diminimalisir dengan mengalihkan risiko kepada perusahaan asuransi. Sebagai lembaga pengalih risiko serta yang mengumpulkan dana premi, perusahaan asuransi harus memperhitungkan risiko yang tidak dapat dimitigasi. Untuk mengoptimalkan manajemen risiko bencana dapat dilakukan analisis mengenai kejadian bencana. Salah satunya dapat menggunakan extreme value theory pada data bencana. Sumber data yang digunakan dalam penelitian ini adalah data gempa bumi di Indonesia tahun 1990-2021 dari katalog gempa bumi National Earthquake Information Center - USGS. Pada penelitian ini digunakan variabel magnitude gempa bumi di Indonesia. Magnitude gempa bumi merupakan kejadian ekstrem. Dalam penelitian ini parameter nilai ekstrem  di estimasi mengikuti distribusi GEV dengan metode Maximum Likelihood Estimation (MLE). Setelah diperoleh parameter dari distribusi GEV akan dilanjutkan estimasi dengan Markov Chain Monte Carlo .  Berdasarkan  pemilihan model terbaik dengan metode DIC, yang diperoleh model terbaik yakni pada data ekstrem maksimum periode 12 bulan. Penelitian dilanjutkan dengan melakukan  perhitungan return level . Hasil dari penelitian ini diharapkan dapat digunakan sebagai referensi dalam perencanaan manajemen risiko bencana alam gempa bumi.
Kata Kunci —Gempa Bumi, Generalized Extreme Value, Magnitude, Markov Chain Monte Carlo, Return Level.

## I. PENDAHULUAN
EMPA bumi merupakan bencana yang termasuk dalam kategori bencana yang disebabkan oleh aktivitas geofisika. Adanya pergerakan lempeng tektonik yang merupakan bagian lapisan kerak bumi keras yang mengapung di atas bagian lapisan astenofer menyebabkan terjadinya gempa bumi. Dengan lempeng utama yang paling berpengaruh pada aktivitas tektonik di Indonesia yakni lempeng Eurasia, lempeng Indo-Australia, lempeng Pasifik, dan lempeng Filipina  [1].  Salah satu gempa bumi besar selama 30 tahun terakhir di Indonesia pernah terjadi di Aceh pada tahun 2004. Gempa bumi berkekuatan sebesar 9.1 magnitudo di Aceh ini menyebabkan tsunami yang memakan lebih dari 160 orang korban jiwa serta 63.977 keluarga kehilangan tempat tinggal. Besar klaim asuransi di Aceh
akibat kejadian tersebut sebesar 3,2 triliun. Risiko kerugian juga dapat diminimalisir dengan adanya mitigasi bencana sebelum terjadinya peristiwa bencana.
Untuk optimalisasi manajemen risiko bencana dapat dilakukan analisis mengenai kejadian bencana. Salah satunya dapat menggunakan extreme value theory pada data bencana. Analisis extreme value bertujuan untuk mengukur perilaku stokastik dari suatu kejadian yang sangat besar atau sangat kecil Parameter yang diestimasi dengan distribusi nilai ekstrem merupakan parameter ditentukan dari data maksimum waktu tertentu, sehingga data yang digunakan terbatas. Untuk mengatasi keterbatasan data yang tersedia digunakan simulasi dengan Bayesian framework untuk estimasi parameter.
Oleh karena itu pada penelitian ini peneliti ingin mengetahui karakteristik data magnitude gempa bumi, hasil estimasi parameter dari data magnitude gempa bumi di Indonesia menggunakan GEV berbasis simulasi MCMC, serta hasil return level magnitude gempa bumi di Indonesia berdasarkan metode GEV berbasis simulasi MCMC. Penelitian ini dilakukan menggunakan data magnitude gempa bumi di Indonesia yang terjadi sejak 1990-2021.
Pemodelan Magnitude Gempa Bumi di Indonesia Menggunakan Generalized Extreme Value (GEV) Berbasis Simulasi
Markov Chain Monte Carlo (MCMC)
Fadhela Anindya, dan Pratnya Paramitha Oktaviana Departemen Aktuaria, Institut Teknologi Sepuluh Nopember (ITS) e-mail : paramita@statistika.its.ac.id
G
Gambar 1. Ilustrasi pengambilan sampel dengan block maxima .
Gambar 2. Histogram magnitude gempa bumi tahun 1990-2021.

## II. TINJAUAN PUSTAKA

## A. Extreme Value Theory
Teori nilai ekstrem  merupakan metode yang dapat digunakan untuk menganalisis perilaku stokastik dari suatu proses yang bernilai sangat besar atau sangat kecil. Metode ini berfokus pada ekor suatu distribusi. Penerapan analisis ekstrem  banyak digunakan untuk penilaian risiko pasar keuangan, pemodelan gelombang laut, aktivitas meteorologi, dan sebagainya.  Dalam analisis nilai ekstrem terdapat pendekatan dengan distribusi GEV dan GPD  [2]. Dalam identifikasi nilai ekstrem terdapat dua pendekatan yakni adalah Block Maxima dan Peak Over Threshold .

## B. Block Maxima
Pada metode Block Maxima nilai ekstrem diperoleh dengan membagi data menjadi beberapa blok periode waktu tertentu.
Di ilustrasikan pada Gambar 1 data dibagi menjadi 5 blok dengan masing-masing blok berisi 5 data. Data yang digunakan untuk dilanjutkan analisis adalah data maksimum yang sudah ditandai hitam. Blok menunjukkan interval waktu yang ditentukan misalnya tahunan, bulanan, dan sebagainya.

## C. Mann Kendal Tren Test
Digunakan untuk mengidentifikasi adanya tren pada data time series , dengan hipotesis:
𝐻𝐻 0 : 𝑆𝑆 = 0 (Tidak terdapat pola tren pada sampel nilai ekstrem) 𝐻𝐻 0 : 𝑆𝑆 ≠ 0 (Terdapat pola tren pada sampel nilai ekstrem) Jika diketahui data 𝑥𝑥 1 , 𝑥𝑥 2 , … 𝑥𝑥 𝑛𝑛 dengan periode waktu 𝑡𝑡 = 1,2, … 𝑛𝑛 , berikut nilai statistik uji Mann Kendal :
𝑍𝑍 = � ( 𝑆𝑆 − 1)/ 𝑉𝑉𝑉𝑉𝑉𝑉 ( 𝑆𝑆 ) 0 , 5 0 ( 𝑆𝑆 + 1)/ 𝑉𝑉𝑉𝑉𝑉𝑉 ( 𝑆𝑆 ) 0 , 5 � untuk 𝑆𝑆 > 0 untuk 𝑆𝑆 = 0 untuk 𝑆𝑆 < 0 (1)
dimana nilai S yakni:
𝑆𝑆 = ∑ ∑ 𝑠𝑠𝑠𝑠𝑛𝑛 ( 𝑥𝑥 𝑗𝑗 − 𝑥𝑥 𝑘𝑘 ) 𝑛𝑛 𝑗𝑗−𝑘𝑘+1 𝑛𝑛−1 𝑘𝑘−1 (2)
dan nilai sgn diperoleh dari persamaan:
𝑠𝑠𝑠𝑠𝑛𝑛�𝑥𝑥 𝑗𝑗 − 𝑥𝑥 𝑘𝑘 � = � 1, untuk 𝑥𝑥 𝑗𝑗 − 𝑥𝑥 𝑘𝑘 > 0 0, untuk 𝑥𝑥 𝑗𝑗 − 𝑥𝑥 𝑘𝑘 = 0 − 1, untuk 𝑥𝑥 𝑗𝑗 − 𝑥𝑥 𝑘𝑘 < 0 (3)
dengan memberikan nilai -1, 0, atau 1 berdasarkan ( 𝑥𝑥 𝑗𝑗 − 𝑥𝑥 𝑘𝑘 ) untuk 𝑗𝑗 > 𝑘𝑘 , dengan 𝑗𝑗 = 𝑡𝑡 + 1 hingga 𝑛𝑛 dan 𝑘𝑘 = 1 hingga 𝑛𝑛 − 1 . Dengan nilai 𝑉𝑉𝑉𝑉𝑉𝑉 ( 𝑆𝑆 ) mengikuti persamaan.
𝑉𝑉𝑉𝑉𝑉𝑉 ( 𝑆𝑆 ) = 1 18 �𝑛𝑛 ( 𝑛𝑛 − 1)(2 𝑛𝑛 + 5) − ∑ 𝑞𝑞 𝑖𝑖 ( 𝑞𝑞 𝑖𝑖 − 1)(2 𝑞𝑞 𝑖𝑖 + 5) 𝑔𝑔 𝑖𝑖=1 � (4)
dimana 𝑠𝑠 adalah jumlah kelompok data, 𝑞𝑞 𝑖𝑖 merupakan jumlah data pada kelompok ke- 𝑖𝑖 . Jika | 𝑍𝑍 | > 𝑍𝑍 𝛼𝛼 2 � atau 𝑝𝑝 − 𝑣𝑣𝑉𝑉𝑣𝑣𝑣𝑣𝑣𝑣 < 𝛼𝛼 maka keputusan tolak H 0 [3].
D. Generalized Extreme Value Terdapat tiga tipe distribusi extreme value yakni tipe I, II, dan III. Ketiga tipe tersebut dikenal dengan Gumbel, Fréchet, dan Weibull.
𝐼𝐼 ∶ 𝐺𝐺 ( 𝑥𝑥 ) = exp �− exp �− � 𝑥𝑥−𝑏𝑏 𝑎𝑎 ��� − ∞ < 𝑥𝑥 < (5)
Gambar 3. Normal probability plot magnitude gempa bumi tahun 1990- 2021 di indonesia.
Gambar 4. Scatter plot pemilihan data nilai ekstrem magnitude gempa bumi dengan block maxima periode 12 bulan.
Gambar 5. DAG data magnitude gempa bumi tahun 1990-2021.
Tabel 1.
Konversi satuan magnitude Skala M agnitude Persamaan Konversi Range Data 𝑀𝑀 𝑠𝑠 𝑀𝑀 𝑤𝑤 = 0,6016 𝑀𝑀 𝑠𝑠 + 2,476 2,8 ≤ M s ≤ 6,1 𝑀𝑀 𝑤𝑤 = 0,69239 𝑀𝑀 𝑠𝑠 + 0,5671 6,2 ≤ M s ≤ 8,7 𝑚𝑚 𝑏𝑏 𝑀𝑀 𝑤𝑤 = 1,0107 𝑚𝑚 𝑏𝑏 + 0,0801 3,7 ≤ m b ≤ 8,2 𝑀𝑀 𝐿𝐿 𝑀𝑀 𝑤𝑤 = 𝑀𝑀 𝐿𝐿 3,0 ≤ M L ≤ 5,8 Tabel 2. Struktur data penelitian Tahun Magnitude 1990 X 1 1990 X 2 ⋮ ⋮ 1990 X 11 1990 1991 X 12 X 13 ⋮ ⋮ 2021 X n Tabel 3.
Statistika deskriptif magnitude gempa bumi tahun 1990-2021 di indonesia Magnitude Gempa Bumi Minimum 2,700 Kuartil 1 4,325 Median 4,628 Rata-rata 4,644 Kuartil 3 4,830 Maksimum 9,100
𝐼𝐼𝐼𝐼 ∶ 𝐺𝐺 ( 𝑥𝑥 ) = � 0 exp �− � 𝑥𝑥−𝑏𝑏 𝑎𝑎 � −𝛼𝛼 � 𝑥𝑥 ≤ 𝑏𝑏 𝑥𝑥 > 𝑏𝑏 (6)
𝐼𝐼𝐼𝐼𝐼𝐼 ∶ 𝐺𝐺 ( 𝑥𝑥 ) = � exp �− �− � 𝑥𝑥−𝑏𝑏 𝑎𝑎 � −𝛼𝛼 �� 1 𝑥𝑥 < 𝑏𝑏 𝑥𝑥 ≥ 𝑏𝑏
(7)
dengan parameter 𝑉𝑉 > 0, 𝑏𝑏 bilangan riil, pada distribusi tipe II dan III 𝛼𝛼 > 0.
Ketiga parameter model tersebut memiliki kelemahan yakni  pada pemilihan tipe yang paling sesuai. Dalam menyelesaikan permasalahan tersebut dikembangkan analisis yang merupakan kombinasi dari ketiganya yang disebut dengan distribusi Generalized Extreme Value (GEV) dengan fungsi distribusi sebagaimana Persamaan (8).
𝐺𝐺 ( 𝑥𝑥 ) = exp �− � 1 + 𝜉𝜉 � 𝑥𝑥−𝜇𝜇 𝜎𝜎 �� −1 𝜉𝜉 � � (8)
Ketiga parameter tersebut yakni location parameter ( 𝜇𝜇 ), scale parameter ( 𝜎𝜎 ),  dan shape parameter ( 𝜉𝜉 ).  Shape parameter pada tiap tipe distribusi yakni, 𝜉𝜉 = 0 untuk distribusi tipe I, 𝜉𝜉 > 0 untuk distribusi tipe II, dan 𝜉𝜉 < 0 untuk distribusi tipe III [2].

## E. Estimasi Parameter
1) Maximum Likelihood Estimation
Merupakan metode untuk mengestimasi parameter model yang mudah beradaptasi dengan  perubahan model  [2]. Estimasi parameter dengan metode MLE dilakukan
Tabel 4.
Uji mann kendall test Periode P-Value 12 bulan 0,3044 6 bulan 0,6091 3 bulan 0,5689
Tabel 5. Hasil estimasi parameter GEV dengan MLE Periode Keterangan 𝜇𝜇̂ 𝜎𝜎� 𝜉𝜉̂ 12 Bulan Parameter 7,2593 0,5144 -0,0860 SE 0,1036 0,0750 0,1419 CI [7,0562; 7,4623] [0,3674; 0,6614] [-0,3641; 0,1921] 6 Bulan Parameter 6,9230 0,4938 -0,0288 SE 0,0704 0,0514 0,1009 CI [6,7850; 7,0610] [0,3930; 0,5945] [-0,2266; 0,1699] 3 Bulan Parameter 6,5301 0,4855 0,0227 SE 0,0490 0,0362 0,0723 CI [6,4340; 6,6261] [0,4145; 0,5565] [-0,1190; 0,1644] Tabel 6. Hasil uji kolmogorov smirnov Periode P-Value 12 bulan 0.7632803 6 bulan 0.6714101 3 bulan 0.8858401 Tabel 7. Hasil estimasi parameter MCMC Periode Parameter Rata- rata Standar Deviasi
MC
Error 12 Bulan 𝜇𝜇̂ 7,2590 0,0301 0,0005 𝜎𝜎� 0,5409 0,0774 0,0006 𝜉𝜉̂ -0,0790 0,0800 0,0010 6 Bulan 𝜇𝜇̂ 6,9210 0,0291 0,0004 𝜎𝜎� 0,5081 0,0502 0,0004 𝜉𝜉̂ -0,0238 0,0673 0,0010 3 Bulan 𝜇𝜇̂ 6,5300 0,0262 0,0002 𝜎𝜎� 0,4932 0,0345 0,0003 𝜉𝜉̂ 0,0254 0,0560 0,0010
Tabel 8.
Pemilihan model terbaik Periode DIC 12 bulan 58,32 6 bulan 113,30 3 bulan 227,00
Tabel 9.
Prediksi Return Level Tahun Return Level 5 8,1192 10 8,5802 20 9,0453 30 9,4325
Gambar  6. Time Series dari history chain data ekstrem magnitude gempa bumi maksimal dalam 12 bulan.
Gambar  7. Time Series dari history  chain data ekstrem magnitude gempa bumi maksimal dalam 6 bulan.
Gambar 8. Time Series dari history chain data ekstrem magnitude gempa bumi maksimal dalam 3 bulan.
denganmemaksimumkan fungsi likelihood , dengan langkah sebagai berikut:
a. Menentukan fungsi likelihood dari distribusi GEV
𝐿𝐿 ( 𝜇𝜇 , 𝜎𝜎 , 𝜉𝜉 | 𝑥𝑥 1 , 𝑥𝑥 2 , … 𝑥𝑥 𝑛𝑛 ) (9)
b. Membuat fungsi ln-likelihood
ln 𝐿𝐿 ( 𝜇𝜇 , 𝜎𝜎 , 𝜉𝜉 | 𝑥𝑥 1 , 𝑥𝑥 2 , … 𝑥𝑥 𝑛𝑛 ) (10)
c. Memaksimumkan fungsi likelihood dengan menurunkan ln- likelihood untuk masing-masing parameter Jika hasil estimasi parameter dengan MLE tidak closed form , maka diperlukan analisis numerik lebih lanjut dengan metode Newton  Raphson [4]. Metode Newton Raphson digunakan turunan pertama dari fungsi ln likelihood fungsi GEV terhadap masing-masing parameternya. Berikut persamaan Newton Raphson :
𝜃𝜃 𝑖𝑖+1 = 𝜃𝜃 𝑖𝑖 − 𝑠𝑠 ( 𝜃𝜃 𝑖𝑖 ) 𝐻𝐻 −1 ( 𝜃𝜃 𝑖𝑖 ) (11)
dimana 𝑠𝑠 ( 𝜃𝜃 𝑖𝑖 ) merupakan vektor gradien berukuran 1 𝑥𝑥𝑝𝑝 dimana 𝑝𝑝 adalah jumlah parameter dan 𝐻𝐻 ( 𝜃𝜃 ) merupakan matriks Hessian berukuran 𝑝𝑝𝑥𝑥𝑝𝑝 yang berisi turunan kedua
terhadap masing-masing parameternya. Iterasi Newton Raphson memiliki nilai awal 𝜃𝜃 0 yang merupakan vector yang elememennya berisi 𝜇𝜇̂ 0 , 𝜎𝜎� 0 ,  dan 𝜉𝜉̂ 0 . Nilai estimasi awal disubstitusikan pada vektor gradien dan matriks Hessian [5].
2) Analisis Bayesian
Pendekatan Bayesian merupakan proses pengambilan keputusan berdasarkan informasi baru dari data yang diamati dan pengetahuan sebelumnya. Diasumsikan 𝜃𝜃 random variabel yang disebut dengan distribusi prior. Terdapat empat macam distribusi prior, c onjugate prior dan non conjugate prior, dimana penentuan prior berdasarkan pola likelihood data, i nformative prior dan non informative prior , penentuan prior yang berdasarkan pada ada atau tidaknya informasi pola distribusi data dari penelitian sebelumnya, p roper prior dan improper prior , penentuan prior berdasarkan pada densitas pada tiap titik, apakah terdistribusi secara uniform atau tidak, p seudo prior , penentuan prior yang pemberian nilainya didasarkan pada hasil elaborasi cara frequentist [6] . Pendekatan Bayesian menggunakan informasi yang diperoleh dari penggabungan dari fungsi likelihood yang
Gambar  11. Statistik gelman-rubin data ekstrem magnitude gempa bumi maksimal dalam 3 bulan.
Gambar  12. Grafik distribusi data ekstrem magnitude gempa bumi maksimal dalam 12 bulan masing-masing parameter.
location chains 1 : 3 start-iteration 4051 5000 6000 7000 scale chains 1 : 3 start-iteration 4051 5000 6000 7000 shape chains 1 : 3 start-iteration 4051 5000 6000 7000 (a) (b) (c) location sample: 21000 location 7.1 7.2 7.3 7.4 scale sample: 21000 scale 0.2 0.4 0.6 0.8 1.0 shape sample: 21000 shape -0.6 -0.4 -0.2 5.55E-17 0.2 (a) (b) (c) Gambar 9. Statistik gelman-rubin data ekstrem magnitude gempa bumi maksimal dalam 12 bulan.
Gambar  10. Statistik gelman-rubin data ekstrem magnitude gempa bumi maksimal dalam 6 bulan. location chains 1 : 3 start-iteration 4051 5000 6000 7000 scale chains 1 : 3 start-iteration 4051 5000 6000 7000 shape chains 1 : 3 start-iteration 4051 5000 6000 7000 (a) (b) (c) location chains 1 : 3 start-iteration 4051 5000 6000 7000 scale chains 1 : 3 start-iteration 4051 5000 6000 7000 shape chains 1 : 3 start-iteration 4051 5000 6000 7000 (a) (b) (c)
diperoleh dari informasi dari data yang diamati 𝑦𝑦 dan informasi dari distribusi prior.
𝑓𝑓 ( 𝜃𝜃 | 𝑦𝑦 ) = 𝑓𝑓 ( 𝜃𝜃 ) 𝑓𝑓 ( 𝑦𝑦 | 𝜃𝜃 ) ∫ Θ 𝑓𝑓 ( 𝜃𝜃 ) 𝑓𝑓� 𝑦𝑦 � 𝜃𝜃 �𝑑𝑑𝜃𝜃 (12)
Keterangan  : 𝑓𝑓 ( 𝑦𝑦 | 𝜃𝜃 ) : Fungsi likelihood 𝑓𝑓 ( 𝜃𝜃 | 𝑦𝑦 ) : Distribusi posterior 𝑓𝑓 ( 𝜃𝜃 ) : Distribusi prior Dalam  model yang lebih kompleks parameter 𝜃𝜃 dapat menyebabkan permasalahan pada integrasi yang rumit. Dapat dilakukan simulasi Markov Chain Monte Carlo [2] . Rantai markov merupakan proses stokastik pada �𝜃𝜃 ( 1 ) , 𝜃𝜃 ( 2 ) , … , 𝜃𝜃 ( 𝑇𝑇 ) � dengan:
𝑓𝑓�𝜃𝜃 ( 𝑡𝑡+1 ) | 𝜃𝜃 ( 𝑡𝑡 ) , … , 𝜃𝜃 ( 1 ) � = 𝑓𝑓�𝜃𝜃 ( 𝑡𝑡+1 ) �𝜃𝜃 ( 𝑡𝑡 ) � (13)
dimana nilai 𝜃𝜃 ( 𝑡𝑡+1 ) ditentukan berdasarkan pada nilai 𝜃𝜃 ( 𝑡𝑡 ) .
Untuk membangkitkan distribusi posterior dimulai dengan menentukan initial value 𝜃𝜃 ( 0 ) , kemudian membangkitkan nilai 𝑇𝑇 hingga konvergen. Monte Carlo merupakan metode simulasi berdasarkan sampel acak dari suatu distribusi. Dengan membangkitkan sampel 𝜃𝜃 ( 1 ) , 𝜃𝜃 ( 2 ) , … , 𝜃𝜃 ( 𝑡𝑡 ) dari distribusi posterior. Dihitung rata-rata sampel dari setiap fungsi dari parameter.
𝐼𝐼̂ = 1 𝑇𝑇 ∑ 𝐺𝐺�𝜃𝜃 ( 𝑡𝑡 ) � 𝑇𝑇 𝑡𝑡=1 (14)
Simulasi digunakan untuk estimasi dan menggambarkan distribusi posterior dari 𝐺𝐺 ( 𝜃𝜃 ) itu sendiri dengan esimasi dari nilai 𝐺𝐺�𝜃𝜃 ( 1 ) � , 𝐺𝐺�𝜃𝜃 ( 2 ) � , … , 𝐺𝐺�𝜃𝜃 ( 𝑇𝑇 ) � [7]. Simulasi MCMC dilakukan dengan membangkitkan variabel acak berdasarkan rantai  Gibbs Sampling yang merupakan salah satu teknik yang digunakan dalam simulasi MCMC untuk menghasilkan variabel acak dari suatu distribusi tertentu tanpa harus mengetahui fungsi densitasnya. Dengan mensimulasikan sampel yang cukup besar, karakteristik dari fungsi densitas dapat diperoleh [8].

## F. Uji Kesesuaian Distribusi
Untuk mengetahui apakah sampel dari sebuah distribusi yang tidak diketahui mengikuti suatu fungsi distribusi tertentu. Berikut hipotesis dari uji Kolmogorov-Smirnov . 𝐻𝐻 0 : Data mengikuti suatu distribusi teoritis tertentu 𝐹𝐹 0 ( 𝑥𝑥 ) 𝐻𝐻 0 : Data tidak mengikuti suatu distribusi teoritis tertentu 𝐹𝐹 0 ( 𝑥𝑥 ) dengan statistik uji:
𝐷𝐷 = max 1≤i≤N �𝐹𝐹 ( 𝑋𝑋 𝑖𝑖 ) − 𝑖𝑖−1 𝑁𝑁 , 𝑖𝑖 𝑁𝑁 − 𝐹𝐹 ( 𝑋𝑋 𝑖𝑖 ) � (15)
dimana 𝑁𝑁 adalah banyak data  dan 𝑖𝑖 adalah urutan data. Keputusan diambil dengan membandingkan 𝐷𝐷 dengan 𝑑𝑑 tabel Kolmogorov-Smirn ov [9].

## G. Deviance Information Criterion
DIC merupakan metode pemilihan model pengembangan dari AIC yang banyak digunakan untuk model yang lebih kompleks.
𝐷𝐷𝐼𝐼𝐷𝐷 = 𝐷𝐷 ( 𝜃𝜃̅ ) + 2 𝑝𝑝 𝐷𝐷 = 𝐷𝐷� + 𝑝𝑝 𝐷𝐷 (16)
dengan,
𝐷𝐷 ( 𝜃𝜃̅ ) : Nilai ekspektasi deviance posterior 𝐷𝐷� : Nilai ekspektasi deviance atas posterior dan nilai 𝑝𝑝 𝐷𝐷 = 𝐷𝐷� − 𝐷𝐷 ( 𝜃𝜃̅ ) . Model dengan nilai DIC terkecil adalah model yang lebih baik.

## H. Return Level
Merupakan nilai maksimum diharapkan dicapai dalam jangka waktu tertentu.  Return level untuk GEV stasioner dinyatakan oleh persamaan.
𝑧𝑧 𝑝𝑝 = � 𝜇𝜇 − 𝜎𝜎 𝜉𝜉 � 1 − { − log (1 − 𝑝𝑝 } −𝜉𝜉 � , 𝜇𝜇 − 𝜎𝜎 log{ − log (1 − 𝑝𝑝 }, 𝜉𝜉 ≠ 0, 1 1 𝜉𝜉 = 0; (17)
Sedangkan persamaan return level untuk GEV non- stasioner dinyatakan dengan persamaan:
𝑧𝑧 𝑝𝑝 = � 𝜇𝜇 ( 𝑡𝑡 ) − 𝜎𝜎 ( 𝑡𝑡 ) 𝜉𝜉 ( 𝑡𝑡 ) � 1 − { − log (1 − 𝑝𝑝 )} −𝜉𝜉 ( 𝑡𝑡 ) � , 𝜇𝜇 ( 𝑡𝑡 ) − 𝜎𝜎 ( 𝑡𝑡 ) log{ − log (1 − 𝑝𝑝 )}, 𝜉𝜉 ( 𝑡𝑡 ) ≠ 0, 1 1 𝜉𝜉 ( 𝑡𝑡 ) = 0; (18)
Gambar  13. Grafik distribusi data ekstrem magnitude gempa bumi maksimal dalam 6 bulan masing-masing parameter.
Gambar  14. Grafik distribusi data ekstrem magnitude gempa bumi maksimal dalam 3 bulan masing-masing parameter. location sample: 21000 location 6.8 6.9 7.0 7.1 scale sample: 21000 scale 0.2 0.4 0.6 0.8 shape sample: 21000 shape -0.4 -0.2 0.0 0.2 0.4 (a) (b) (c) location sample: 21000 location 6.4 6.5 6.6 6.7 scale sample: 21000 scale 0.3 0.4 0.5 0.6 0.7 shape sample: 21000 shape -0.2 0.0 0.2 0.4 (a) (b) (c)
dimana 𝑧𝑧 𝑝𝑝 adalah return level dan 𝑝𝑝 adalah probabilitas tercapainya 𝑧𝑧 𝑝𝑝 [2].

## I. Gempa Bumi di Indonesia
Pertemuaan antar lempeng yang berpotensi menyebabkan terjadinya  aktivitas kegempaan  [10].  Pada mulanya skala gempa bumi dikenal dengan skala richter atau M L . Untuk menjangkau daerah yang lebih luas, dikembangkan skala pengukuran lain yang dapat mendeteksi gempa untuk jarak yang lebih luas dengan menggunakan lebih banyak variabel dalam pengukurannya seperti skala M w , m b ,  dan M S . serta dapat mendeteksi gempa bumi yang memiliki magnitude besar. Dalam penelitian ini, skala magnitude dikonversi dengan persamaan konversi yang merupakan hasil studi hubungan antar magnitude untuk analisis seismisitas di Indonesia. Pada penelitian ini magnitude dikonversi ke skala M w [8] ditunjukkan pada Tabel 1.

## III. METODOLOGI PENELITIAN

## A. Sumber Data dan Variabel Penelitian
Sumber data dalam penelitian ini adalah data gempa bumi di Indonesia tahun 1990-2021  dari katalog gempa bumi National Earthquake Information Center -  United States Geological Survey. Variabel yang digunakan dalam penelitian ini adalah data magnitude gempa bumi. Struktur data pada penelitian ditampilkan pada Tabel 2.

## B. Langkah Analitis
Langkah-langkah analisis pada penelitian ini diantaranya yaitu sebagai berikut:
1. Melakukan pengumpulan data
2. Melakukan data pre processing
3. Melakukan analisis statistika deskriptif
4. Membagi data dalam menjadi beberapa blok. Dalam penelitian ini terdapat 3 periode waktu dengan masing- masing periode yakni 12 bulan, 6 bulan, dan 3 bulan
5. Memilih sampel nilai ekstrem berdasarakan magnitude gempa bumi terbesar tiap blok
6. Melakukan uji stasioneritas data dengan Mann Kendal trend test .
7. Melakukan estimasi parameter GEV stasioner jika data tidak terdapat tren.
8. Uji kesesuaian distribusi dengan uji Kolmogorov Smirnov.
9. Melakukan estimasi parameter GEV menggunakan Bayesian dengan MCMC.
10. Menentukan model terbaik dengan DIC
11. Menentukan return level dengan persamaan return level untuk GEV stasioner 12. Menarik kesimpulan

## IV. HASIL DAN ANALISIS
A. Karakteristik Data Magnitude Gempa Bumi Tahun 1990-2021 di Indonesia
Dalam penelitian ini, satuan dari magnitude dikonversi menjadi moment  magnitude .  Statistika deskriptif dari data magnitude gempa bumi tahun 1990-2021 pada Tabel 3. Nilai kuartil 3 magnitude gempa bumi dengan nilai maksimum magnitude gempa bumi terdapat perbedaan yang
cukup jauh yang menunjukkan adanya indikasi data ekstrem. Identifikasi adanya nilai ekstrem dapat dilakukan secara visual melalui histogram pada Gambar 2.
Berdasarkan histogram Gambar 2  dapat diidentifikasi adanya pola data heavy tail pada data magnitude gempa bumi. Dilakukan pengujian normalitas data dengan normal probability plot dan uji kesesuaian distribusi Kolmogorov Smirnov dengan distribusi normal sebagai distribusi teoritis, Berikut normal probability plot dari data magnitude gempa bumi.
Berdasarkan Gambar 3 menunjukkan bahwa sebaran data yakni titik-titik berwarna hitam tidak mengikuti garis linier, yang menunjukkan bahwa data magnitude gempa bumi tidak berdistribusi normal. Dilakukan uji statistik untuk mengetahui normalitas data dengan uji Kolmogorov-Smirnov dengan hipotesis sebagai berikut:
𝐻𝐻 0 :  Data ekstrem magnitude gempa bumi mengikuti distribusi normal 𝐻𝐻 0 : Data ekstrem magnitude gempa bumi tidak mengikuti distribusi normal Dari hasil uji Kolmogorov-Smirnov diperoleh hasil 𝑝𝑝 − 𝑣𝑣𝑉𝑉𝑣𝑣𝑣𝑣𝑣𝑣 bernilai 0 dan 𝐷𝐷 ℎ𝑖𝑖𝑡𝑡𝑖𝑖𝑛𝑛𝑔𝑔 0,9990. Dengan sampel sebanyak 49.642 maka nilai 𝐷𝐷 𝑡𝑡𝑎𝑎𝑏𝑏𝑡𝑡𝑡𝑡 0,0061. Berdasarkan hasil uji Kolmogorov-Smirnov 𝑝𝑝 − 𝑣𝑣𝑉𝑉𝑣𝑣𝑣𝑣𝑣𝑣 < 0.05 atau 𝐷𝐷 ℎ𝑖𝑖𝑡𝑡𝑖𝑖𝑛𝑛𝑔𝑔 > 𝐷𝐷 0 . 05 maka keputusan yang diperoleh adalah tolak 𝐻𝐻 0 sehingga  data magnitude gempa bumi di Indonesia antara tahun 1990 sampai 2021 tidak berdistribusi normal.
B. Pemilihan Nilai Ekstream dengan Block Maxima
Pada penelitian ini digunakan metode Block Maxima dengan periode 12 bulan, 6 bulan, dan 3 bulan untuk tiap blok. Block Maxima dengan langkah pengambilan sampel data sebagai berikut:
1. Membagi data menjadi beberapa blok. Pada penelitian ini data dibagi menjadi 32 blok untuk model pertama, 64 blok untuk model kedua dan 128 blok untuk model ketiga.
2. Menentukan nilai maksimal pada masing-masing blok. Pada Gambar 4 ditunjukkan salah satu ilustrasi pemilihan data dengan  menggunakan  metode Block Maxima dengan periode 12 bulan.

## C. Uji Stasioneritas Data
Uji stasioneritas pada data ekstrem dapat dilakukan dengan Mann Kendall Trend Test dengan hipotesis sebagai berikut: 𝐻𝐻 0 : 𝑆𝑆 = 0 (Tidak terdapat pola tren pada data ekstrem magnitude gempa bumi)
𝐻𝐻 0 : 𝑆𝑆 ≠ 0 (Terdapat pola tren pada data ekstrem magnitude gempa bumi) Tabel 4 menunjukan hasil uji stasioneritas data magnitude gempa bumi. Dari hasil uji yang diperoleh data ekstrem magnitude gempa bumi  pada masing-masing periode, memiliki nilai p-value < α sehingga tidak terdapat pola tren pada sampel data magnitude gempa bumi ekstrem di Indonesia. Maka digunakan analisis dengan Distribusi GEV untuk data ekstrem yang stasioner atau tidak memiliki pola tren dalam penelitian ini.
D. Estimasi Parameter Generalized Extreme Value Estimasi dilakukan dengan software RStudio menggunakan  metode Maximum Likelihood Estimation
(MLE). Diperoleh hasil estimasi parameter GEV pada Tabel 5.
Dari hasil estimasi dengan MLE fungsi distribusi kumulatif dari distribusi GEV  untuk 32 sampel nilai ekstrem untuk data maksimum pada periode 12 bulan dengan yakni seperti pada persamaan berikut.
𝐺𝐺 ( 𝑥𝑥 ) = exp �− � 1 − 0,0860 � 𝑥𝑥 − 7,2593 0,5144 �� 1 0 , 0860 � �
Berikut  fungsi distribusi kumulatif dari 64 sampel nilai ekstrem. Data berdasarkan pemilihan sampel data maksimum tiap 6 bulan dengan fungsi distribusi kumulatif berikut:
𝐺𝐺 ( 𝑥𝑥 ) = exp �− � 1 − 0,0288 � 𝑥𝑥 − 6,9230 0,4938 �� 1 0 . 088 � �
Selanjutnya merupakan fungsi distribusi kumulatif dari 128 sampel nilai ekstrem magnitude gempa bumi untuk pemilihan data maksimum tiap 3 bulan dengan fungsi distribusi kumulatif.
𝐺𝐺 ( 𝑥𝑥 ) = exp �− � 1 + 0,0227 � 𝑥𝑥 − 6,5301 0,4855 �� −1 0 , 0227 � �

## E. Uji Kesesuaian Distribusi
Dilakukan uji Kolmogorov-Smirnov dengan hipotesis pengujian sebagai berikut:
𝐻𝐻 0 : Data ekstrem magnitude gempa bumi mengikuti suatu distribusi GEV
𝐻𝐻 0 : Data ekstrem magnitude gempa bumi tidak mengikuti suatu distribusi GEV Dari hasil uji diperoleh hasi uji pada Tabel 6. Berdasarkan Tabel 6  diperoleh keputusan Gagal Tolak H 0 untuk ketiga periode. Hal tersebut menunjukkan bahwa data ekstrem magnitude gempa bumi periode maksimum 12 bulan, 6 bulan, dan 3 bulan mengikuti distribusi GEV.
F. Analisis Bayesian Markov Chain Monte Carlo (MCMC) Dengan menggunakan software OpenBugs, pada Gambar 5 ditunjukkan DAG yang memvisualisasikan hubungan antar parameter pada model digunakan Directed Acyclic Graph (DAG). Analisis Bayesian MCMC dapat dilanjutkan jika time series plot dari  distribusi posterior yang dihasilkan adalah stasioner dan parameter yang dihasilkan sudah konvergen. History chain dari  data ekstrem magnitude gempa bumi ditampilkan pada Gambar 6, Gambar 7 dan Gambar 8. Berdasarkan Gambar 6, Gambar 7 dan Gambar 8 masing- masing parameter untuk data ekstrem magnitude gempa bumi masing-masing parameter konvergen, memiliki grafik yang saling tumpang tindih dan menyebar secara acak dan stabil antar rantai.
Statistik Gelman-Rubin digunakan untuk menilai konvergensi dengan membentuk beberapa rantai dari initial value yang over-dispersed dengan membandingkan variabilitas didalam rantai dan diantara rantai. Jumlah rantai yang dihasilkan dilambangkan oleh M dengan panjang setiap rantai 2 T dan lebar rata-rata interval di rantai sejumlah M dilambangkan dengan W . Lebar B diperoleh dari credible interval 100(1-a)% berdasarkan sampel MT . Dengan rasio
R=B/W dari pooled terhadap rata-rata lebar interval harus lebih besar dari 1 jika initial value over-dispersed secara tepat. Nilai tersebut juga cenderung mendekati 1 saat konvergensi didekati. Berikut merupakan output hasil statistik Gelman-Rubin dari penelitian ini untuk data ekstrem magnitude gempa bumi ketiga periode. Berdasarkan Gambar 9, Gambar 10, dan Gambar 11, nilai statistik Gelman-Rubin berkisar pada  angka 1. Hal ini menunjukkan bahwa masing-masing parameter  pada tiap periode konvergen.
Penentuan prior disesuaikan dengan bentuk distribusi teoritis dengan parameter yang diperoleh dari pseudo prior dari estimasi parameter distribusi GEV dengan metode MLE. Berikut simulasi MCMC diperoleh distribusi tiap parameter.
Berdasarkan Gambar 12, Gambar 13, dan Gambar 14, distribusi parameter location dan parameter scale mengikuti distribusi normal dan shape parameter mengikuti distribusi gamma.
Hasil estimasi parameter magnitude gempa bumi memiliki MC error yang  lebih  kecil dari 5% standar deviasi  hasil estimasi parameter yang dapat dilihat pada Tabel 7.

## G. Pemilihan Model Terbaik
Pemilihan model terbaik dilakukan dengan metode DIC. Dari Tabel 8 diperoleh bahwa model terbaik adalah model dengan blok periode 12 bulan dengan DIC sebesar 58,32.

## H. Return Level
Hasil estimasi parameter GEV yang diperoleh pada simulasi MCMC digunakan untuk menentukan nilai return level pada Tabel 9.
Pada  Tabel 9 menunjukkan estimasi ambang batas maksimum nilai magnitude gempa bumi yang terjadi pada periode T dengan peluang kejadian 1/T. Pada periode 5 tahun dengan probabilitas sebesar 1/5 ambang batas maksimum kejadian gempa bumi diestimasi sebesar 8,1109 magnitude . Pada periode 10 tahun dengan probabilitas sebesar 1/10 ambang batas maksimum kejadian gempa bumi diestimasi sebesar 8,5768 magnitude .  Pada periode 20 tahun dengan probabilitas sebesar 1/20 ambang batas maksimum kejadian gempa bumi diestimasi sebesar 9,0482 magnitude . Pada periode 30 tahun dengan probabilitas sebesar 1/30 ambang batas maksimum kejadian gempa bumi diestimasi sebesar 9,3307 magnitude.

## V. KESIMPULAN DAN SARAN
Dari pembahasan pada bab sebelumnya dapat disimpulkan bahwa gempa bumi  yang terjadi di Indonesia antara tahun 1990 sampai 2001 sebesar memiliki rata-rata magnitude sebesar 4,64. Estimasi parameter data magnitude gempa bumi di Indonesia menggunakan GEV berbasis simulasi MCMC menunjukkan hasil yang baik karena nilai MC Error dibawah 5% dari standar deviasi. Dengan metode pemilihan model terbaik menggunakan DIC diperoleh model terbaik merupakan periode 12 bulan. Dengan hasil estimasi parameter 𝜇𝜇̂ adalah 7,25800; parameter 𝜎𝜎� sebesar 0,54080; dan parameter 𝜉𝜉̂ adalah -0,07954. Perhitungan return level dari model terbaik adalah 8,1108 magnitude untuk periode 5 tahun yang berdasarkan data historis pernah terjadi di kawasan pulau Sumatera. Hasil perhitungan return level untuk periode 10 tahun adalah
8,5768 magnitude , untuk periode 20 tahun 9,0481 magnitude , dan 9,3307 magnitude untuk periode 30 tahun. Gempa berkekuataan diatas 8,5768 magnitude berdasarkan data historis pernah terjadi di kawasan pulau Sumatera bagian utara.
Saran untuk penelitian lebih lanjut yakni dapat dilakukan penelitian dengan menambahkan faktor lain seperti, kedalaman, dan menambahkan adanya faktor segmen patahan gempa, atau keadaan gunung api.  Berdasarkan gambaran kejadian gempa dari hasil penelitian, pemerintah, perusahaan asuransi, serta pihak lain yang terkait dapat mempersiapkan perencanaan mitigasi dan manajemen risiko akibat terjadinya gempa bumi.

## DAFTAR PUSTAKA
[1] S. Rohadi, “Studi seismotektonik sebagi indikator potensi gempabumi di wilayah indonesia,” Jurnal Meteorologi dan Geofisika, vol. 10, no. 2, 2009, doi: 10.31172/jmg.v10i2.40.
[2] S. Coles, An Introduction to Statistical Modeling of Extreme Values. Verlag London: Springer Science & Business Media, 2001, ISBN 9781852334598.
[3] D. Rahmayani and Sutikno, “Analisis curah hujan ekstrim non- stasioner dengan pendekatan block maxima di surabaya dan
mojoker,” Jurnal Sains dan Seni ITS, vol. 8, no. 2, pp. 161–168, 2019, doi: 10.12962/j23373520.v8i2.44133.
[4] R. E. Situmorang, D. A. I. Maruddani, and R. Santoso, “Pembentukan portofolio saham dengan metode markowitz dan pengukuran value at risk berdasarkan generalized extreme value (studi kasus: saham perusahaan the idx top ten blue 2017),” Jurnal  Gaussian, vol. 7, no. 2, pp. 212–223, 2018, doi: 10.14710/j.gauss.7.2.212-223.
[5] Y. D. W. Sari and Sutikno, “Estimasi parameter generalized pareto distribution pada kasus identifikasi perubahan iklim di sentra produksi padi jawa timur,” Jurnal Sains dan Seni POMITS, vol. 2, no. 2, pp. 2337–3520, 2013.
[6] M. N. Abdullah, N. Iriawan, and I. Irhamah, “Pendekatan bayesian untuk analisis survival pada kasus demam berdarah dengue pasien RSUD dr. soetomo surabaya,” Jurnal Sains dan Seni ITS, vol. 9, no. 2, pp. 2337–3520, Feb. 2021, doi: 10.12962/j23373520.v9i2.54448.
[7] C. Robert and I. Ntzoufras, “Bayesian modeling using WinBUGS,” CHANCE, vol. 25, no. 2, pp. 60–61, Apr. 2012, doi: 10.1080/09332480.2012.685377.
[8] G. Casella and E. I. George, “Explaining the gibbs sampler,” Am Stat, vol. 46, no. 3, pp. 167–174, 1992, doi: 10.1080/00031305.1992.10475878.
[9] P. P. Oktaviana and Irhamah, “Kolmogorov-smirnov goodness-of-fit test for identifying distribution of the number of earthquakes and the losses due to earthquakes in indonesia,” J Phys Conf Ser, vol. 1821, no. 1, p. 012045, 2021, doi: 10.1088/1742-6596/1821/1/012045.
[10]  Pusat Studi Gempa Nasional, Peta Sumber dan Bahaya Gempa Indonesia Tahun 2017. Bandung: Kementrian Pekerjaan Umum dan Perumahan Rakyat, 2017.