1 Pendahuluan
Pada kesempatan kali ini, untuk mata kuliah Matematika Aktuaria sudah memasuki pertemuan akhir yang mana setiap mahasiswa mendapat tugas akhir untuk memenuhi syarat kelulusan dari semester 5. Adapun untuk tugas akhirnya diminta untuk menentukan premi asuransi bencana alam yang terjadi di Provinsi Bengkulu sejak tahun 2000 sampai tahun 2020 menggunakan model risiko kolektif dengna premi berdasarkan nilai ekspetasi dan premi berdasarkan nilai standar deviasi. Dengan menggunakan distribusi Poisson akan ditentukan estimasi banyaknya kejadian bencana alam dan menggunakan distribusi Weibull untuk menghitung estimasi besarnya kerugian yang terjadi akibat bencana alam. Dalam melakukan proses perhitungan dan analisis akan menggunakan bantuan bahasa R dengan aplikasi Rstudio.
3 Import Data
# Data kejadian dan kerugian bencana alam
tahun <- c(2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020)
kerugian <- c(13555298411.19, 11029916950.18, 15975856772.64, 15410902354.27, 18985203518.50, 17557268711.70, 10848969106.01, 17395018451.02, 17187708349.74, 19927948832.11, 30431913005.97, 19079543271.32, 18862604188.69, 8834824847.88, 13231795220.00, 23703573276.09, 15743863219.87, 17353216866.92, 22655839742.87, 23433580384.88, 21056987987.80)
kejadian <- c(11, 12, 13, 17, 19, 18, 10, 11, 14, 20, 22, 21, 22, 23, 23, 25, 31, 34, 25, 19, 26)
data <- data.frame(
tahun = tahun,
kerugian = kerugian,
kejadian = kejadian
)
data## tahun kerugian kejadian
## 1 2000 13555298411 11
## 2 2001 11029916950 12
## 3 2002 15975856773 13
## 4 2003 15410902354 17
## 5 2004 18985203518 19
## 6 2005 17557268712 18
## 7 2006 10848969106 10
## 8 2007 17395018451 11
## 9 2008 17187708350 14
## 10 2009 19927948832 20
## 11 2010 30431913006 22
## 12 2011 19079543271 21
## 13 2012 18862604189 22
## 14 2013 8834824848 23
## 15 2014 13231795220 23
## 16 2015 23703573276 25
## 17 2016 15743863220 31
## 18 2017 17353216867 34
## 19 2018 22655839743 25
## 20 2019 23433580385 19
## 21 2020 21056987988 26Data di atas merupakan data bencana alam di provinsi Bengkulu sejak tahun 2000 - 2020 yang mana terdiri dari jumlah kerugian per tahun beserta dengan jumlah kejadiannya dalam tahun tersebut.
# Membuat histogram untuk distribusi kejadian
hist(kejadian, main = "Distribusi Kejadian Bencana Alam", xlab = "Kejadian", ylab = "Frekuensi", col = "skyblue")
# Membuat histogram untuk distribusi kerugian bencana alam
hist(kerugian, main = "Distribusi Kerugian Bencana Alam", xlab = "Kerugian (dalam Rupiah)", ylab = "Frekuensi", col = "salmon")
4 Analisa Deskriptif
Selanjutnya akan dilakukan analisa deskriptif terlebih dahulu untuk mencari tahu gambaran tentang data yang dimiliki, dalam hal ini saya ingin mencari tahu bagaimana menvisualisasikan jumlah kejadian dan kerugian per tahun dalam bentuk yang menarik.
4.1 Jumlah Kejadian Bencana Alam di Indonesia Per Tahun
# Mengatur palet warna untuk setiap bar
palet_warna <- rainbow(length(unique(data$tahun)))
# Membuat plot jumlah kejadian per tahun dengan warna berbeda untuk setiap bar
plot_kejadian <- ggplot(data = data, aes(x = as.factor(tahun), y = kejadian, fill = as.factor(tahun))) +
geom_bar(stat = "identity", color = "black", width = 0.7) +
geom_text(aes(label = kejadian), vjust = -0.5, color = "black", size = 3.5, position = position_dodge(width = 0.7)) +
labs(title = "Jumlah Kejadian Per Tahun",
x = "Tahun",
y = "Jumlah Kejadian") +
scale_fill_manual(values = palet_warna, guide = FALSE) + # Menghilangkan keterangan legend tahun
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Untuk memutar label sumbu x
# Menampilkan plot
plot_kejadian
Berdasarkan grafik tersebut diketahui bahwa kejadian bencana alam di
provinsi Bengkulu paling banyak pada tahun 2017 dengan jumlah sebanyak
34 kejadian. Dan pada tahun 2006 jumlah kejadian bencana paling sedikit
dengan jumlah 10 kejadian. Kemudian, saya ingin mencari tahu secara
deskriptif untuk data kejadian ini:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.00 14.00 20.00 19.81 23.00 34.00## [1] 6.539259Dapat diketahui bahwa nilai rata-rata kejadian bencana alam yang terjadi di provinsi Bengkulu sebanyak 19,81 kejadian per tahunnya. Dengan minimal kejadian 10 dan maksimal kejadian 34, yaitu pada tahun 2017 seperti pada grafik sebelumnya.
4.2 Kerugian Akibat Bencana Alam Per Tahun
# Membuat line chart untuk kerugian per tahun
plot_kerugian <- ggplot(data = data, aes(x = tahun, y = kerugian)) +
geom_line(color = "blue") +
geom_point(color = "indianred", size = 3) +
labs(title = "Kerugian Per Tahun",
x = "Tahun",
y = "Kerugian") +
theme_minimal()
# Menampilkan line chart
print(plot_kerugian)
Berdasarkan grafik di atas, saya dapat menyimpulkan walaupun tahun 2017 jumlah kejadian paling banyak di antara tahun lainnya, tetapi jumlah kerugian tidak terlalu banayk dibandingkan dengan tahun 2010 yang jumlah kejadian tidak sebanyak tahun 2017 tersebut. Menurut saya, hal ini dapat terjadi karena tahun 2017 walaupun kejadiannya banyak, tetapi tidak terlalu berat sehingga kerugian yang dialami pemerintah provinsi tidak begitu besar juga, berbeda dengan tahun 2010 yang kejadiannya lebih sedikit, tetapi kejadinnya berat.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.835e+09 1.541e+10 1.740e+10 1.773e+10 1.993e+10 3.043e+10## [1] 4981967412Kemudian untuk data kerugian dapat diketahui bahwa rata-rata kerugian bencana alam bernilai Rp17.730.000.000 dengan maksimal kerugian sebesar Rp30.431.913.006
5 Distribusi Poisson Untuk Kejadian
rata_kejadian <- mean(data$kejadian)
# Menghitung probabilitas distribusi Poisson
prob_poisson <- ppois(0:max(data$kejadian), lambda = rata_kejadian)
# Membuat histogram distribusi kejadian bencana alam
hist(data$kejadian, breaks = 15, freq = FALSE, main = "Distribusi Kejadian Bencana Alam", xlab = "Kejadian", ylab = "Probabilitas", col = "skyblue")
lines(0:max(data$kejadian), prob_poisson, col = "red") 
6 Distribusi Weibull Untuk Kerugian
## Fitting of the distribution ' weibull ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 3.879656e+00 0.5927046
## scale 1.972176e+10 NaN
## Loglikelihood: -498.5208 AIC: 1001.042 BIC: 1003.131
## Correlation matrix:
## shape scale
## shape 1 NaN
## scale NaN 1params <- coef(fit_weibull)
lambda <- params["shape"]
k <- params["scale"]
expected_value <- k * gamma(1 + 1/lambda)
variance <- k^2 * gamma(1 + 2/lambda) - expected_value^2
expected_value## scale
## 17844991286Didapat nilai E(X(t)) = 17844991286
## scale
## 2.64861e+19Dan nilai Var(X(t)) = 2.64861e+19
# Membuat histogram dari data kerugian
hist(data$kerugian, freq = FALSE, main = "Histogram Data Kerugian dengan Distribusi Weibull")
# Menghasilkan kurva distribusi Weibull
curve(dweibull(x, shape = fit_weibull$estimate[1], scale = fit_weibull$estimate[2]),
col = "blue", lwd = 2, add = TRUE, yaxt = "n")
7 Model Risiko Kolektif
Setelah melakukan uji distribusi Poisson dan Weibull untuk menentukan estimasi dari kejadian dan kerugian, selanjutnya melakukan perhitungan risiko kolektif yang dilakukan dengan proses seperti berikut ini:
7.1 Ekspetasi
## [1] 353500711780Didapat nilai ekspetasinya adalah Rp353.500.711.780
8 Perhitungan Premi Asuransi
Langkah terakhir adalah menghitung besaran premi asuransi dengan 2 metode, yaitu berdasarkan nilai ekspetasi dan berdasarkan nilai standar deviasi, sehingga nantinya dapat membandingkannya.
8.1 Dengan Prinsip Nilai Ekspetasi
# Menginisialisasi nilai a1 hingga a10
a_values <- c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1)
# Menghitung nilai p1 hingga p10 berdasarkan a_values dan Es
p_values <- (1 + a_values) * Es
# Membuat data frame dari nilai a_values dan p_values
data_eks <- data.frame(a = a_values, p_eks = p_values)
# Menampilkan data frame
print(data_eks)## a p_eks
## 1 0.01 357035718898
## 2 0.02 360570726015
## 3 0.03 364105733133
## 4 0.04 367640740251
## 5 0.05 371175747369
## 6 0.06 374710754487
## 7 0.07 378245761604
## 8 0.08 381780768722
## 9 0.09 385315775840
## 10 0.10 3888507829588.2 Dengan Prinsip Standar Deviation
# Menginisialisasi nilai a1 hingga a10
a_values <- c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1)
# Menghitung nilai p1 hingga p10 berdasarkan a_values dan Es
p_values <- Es + a_values * sqrt(Vs)
# Membuat data frame dari nilai a_values dan p_values
data_sd <- data.frame(a = a_values, p_sd = p_values)
# Menampilkan data frame
print(data_sd)## a p_sd
## 1 0.01 354327324993
## 2 0.02 355153938207
## 3 0.03 355980551421
## 4 0.04 356807164634
## 5 0.05 357633777848
## 6 0.06 358460391061
## 7 0.07 359287004275
## 8 0.08 360113617489
## 9 0.09 360940230702
## 10 0.10 361766843916# Nilai a
a <- c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1)
# Nilai p_eks
p_eks <- c(357035718898, 360570726015, 364105733133, 367640740251, 371175747369,
374710754487, 378245761604, 381780768722, 385315775840, 388850782958)
# Nilai p_sd
p_sd <- c(354327324993, 355153938207, 355980551421, 356807164634, 357633777848,
358460391061, 359287004275, 360113617489, 360940230702, 361766843916)
# Membuat data frame dengan menggabungkan nilai a, p_eks, dan p_sd
data_combined <- tibble(a = a, Premi_Ekspetasi = p_eks, Premi_Std_Deviasi = p_sd)
# Menampilkan data frame
print(data_combined)## # A tibble: 10 × 3
## a Premi_Ekspetasi Premi_Std_Deviasi
## <dbl> <dbl> <dbl>
## 1 0.01 357035718898 354327324993
## 2 0.02 360570726015 355153938207
## 3 0.03 364105733133 355980551421
## 4 0.04 367640740251 356807164634
## 5 0.05 371175747369 357633777848
## 6 0.06 374710754487 358460391061
## 7 0.07 378245761604 359287004275
## 8 0.08 381780768722 360113617489
## 9 0.09 385315775840 360940230702
## 10 0.1 388850782958 361766843916Berdasarkan hasil di atas, perhtiungan premi asuransi dengan prinsip nilai ekspetasi memiliki hasil yang lebih besar dibandingkan dengan prinsip standar deviasi. Terlihat bahwa dengan nilai faktor biaya yang awalnya 0.01 atau 1% besasrnya premi terus meningkat juga. Untuk premi terkecil dengan nilai ekspetasi sebesar Rp357.035.718.898 dan untuk nilai standar deviasi sebesar Rp354.327.324.993 dan untuk melihatnya secara visual agar lebih jelas dapat memerhatikan grafik seperti berikut ini:
data_combined <- data.frame(
a = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1),
Premi_Ekspetasi = c(357035718898, 360570726015, 364105733133, 367640740251, 371175747369,
374710754487, 378245761604, 381780768722, 385315775840, 388850782958),
Premi_Std_Deviasi = c(354327324993, 355153938207, 355980551421, 356807164634, 357633777848,
358460391061, 359287004275, 360113617489, 360940230702, 361766843916)
)
# Membuat plot garis dengan titik menggunakan ggplot2
ggplot(data_combined, aes(x = a)) +
geom_line(aes(y = Premi_Ekspetasi, color = "Premi_Ekspetasi"), size = 1) +
geom_point(aes(y = Premi_Ekspetasi, color = "Premi_Ekspetasi"), size = 3) + # Menambahkan titik untuk p_eks
geom_line(aes(y = Premi_Std_Deviasi, color = "Premi_Std_Deviasi"), size = 1) +
geom_point(aes(y = Premi_Std_Deviasi, color = "Premi_Std_Deviasi"), size = 3) + # Menambahkan titik untuk p_sd
labs(title = "Perbandingan Premi Asuransi", x = "Nilai a", y = "Nilai p") +
scale_color_manual(name = "Legend", values = c("Premi_Ekspetasi" = "blue", "Premi_Std_Deviasi" = "red")) +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Berdasarkan data yang disajikan, pola yang terlihat adalah kenaikan nilai premi asuransi bencana alam secara signifikan seiring dengan peningkatan faktor biaya yang diberikan, baik dalam prinsip nilai ekspektasi maupun standar deviasi. Namun, yang menarik adalah perbedaan skala kenaikan premi antara kedua prinsip tersebut. Meskipun kenaikan standar deviasi juga terjadi seiring dengan pertambahan faktor biaya, namun perubahan premi yang dihasilkannya jauh lebih kecil dibandingkan dengan prinsip nilai ekspektasi. Lebih menariknya lagi, pada nilai faktor biaya yang sama, prinsip standar deviasi menghasilkan premi yang lebih rendah dibandingkan dengan prinsip nilai ekspektasi. Dengan demikian, dapat disimpulkan bahwa prinsip standar deviasi bisa menjadi pilihan yang lebih ekonomis dalam perhitungan premi asuransi bencana alam.
9 Alasan Memilih Distribusi Weibull
# Fit distribusi Exp
fit_exp <- fitdist(kerugian, "exp", method = "mme")
# Lakukan uji goodness-of-fit
gof_exp <- gofstat(fit_exp)
# Menampilkan hasil uji goodness-of-fit untuk distribusi Exp
print("Exp distribution:")## [1] "Exp distribution:"## Goodness-of-fit statistics
## 1-mme-exp
## Kolmogorov-Smirnov statistic 0.4101217
## Cramer-von Mises statistic 1.0853520
## Anderson-Darling statistic 5.2387768
##
## Goodness-of-fit criteria
## 1-mme-exp
## Akaike's Information Criterion 1035.130
## Bayesian Information Criterion 1036.175# Fit distribusi Weibull
fit_weibull <- fitdist(kerugian, "weibull", method = "mle")
# Lakukan uji goodness-of-fit
gof_weibull <- gofstat(fit_weibull)
# Menampilkan hasil uji goodness-of-fit untuk distribusi Weibull
print("Weibull distribution:")## [1] "Weibull distribution:"## Goodness-of-fit statistics
## 1-mle-weibull
## Kolmogorov-Smirnov statistic 0.1292899
## Cramer-von Mises statistic 0.0517697
## Anderson-Darling statistic 0.3314834
##
## Goodness-of-fit criteria
## 1-mle-weibull
## Akaike's Information Criterion 1001.042
## Bayesian Information Criterion 1003.131Dari kedua metode tersebut, distribusi Weibull memiliki nilai yang lebih rendah dalam kriteria pengujian Goodness of Fit dengan nilainya lebih rendah dibanding metode sebelumnya, jika nilainya lebih rendah maka dapat diartikan bahwa model yang tercipta lebih sesuai dengan data yang dimiliki. Selain itu, nilai dari statistik Weibull tidak terlalu jauh, sehingga dapat dikatakan konsisten, berbeda dengan metode lain yang cukup jomplang untuk nilai uji statistiknya.