#Constraints Parameter
M = 3
S = 300
L = 6000
R = 60
theta = 0.1
alpha = 0.25
q = 21.41/275.5*(33/92.3)
library(readxl)
setwd("C:/Users/LENOVO/OneDrive/Documents/NATIONAL LIFE/Retensi")
catastrophe = read_excel("BPNB.xlsx", sheet = "Sheet1")
population = read_excel("BPNB.xlsx", sheet = "Sheet2")
head(catastrophe)
## # A tibble: 6 × 4
## Provinsi Bencana Tahun `Jumlah Korban`
## <chr> <chr> <dbl> <dbl>
## 1 Papua Banjir 2019 105
## 2 Sulawesi Barat Tsunami 2021 96
## 3 NTT Banjir 2021 72
## 4 Jawa Barat Tsunami 2022 62
## 5 Sulawesi Selatan Banjir 2019 56
## 6 Jawa Tengah Banjir 2021 50
#Filtering Data
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
catastrophe = catastrophe %>% filter(`Jumlah Korban` >= M)
catastrophe
## # A tibble: 117 × 4
## Provinsi Bencana Tahun `Jumlah Korban`
## <chr> <chr> <dbl> <dbl>
## 1 Papua Banjir 2019 105
## 2 Sulawesi Barat Tsunami 2021 96
## 3 NTT Banjir 2021 72
## 4 Jawa Barat Tsunami 2022 62
## 5 Sulawesi Selatan Banjir 2019 56
## 6 Jawa Tengah Banjir 2021 50
## 7 Riau Tanah Longsor 2023 46
## 8 NTT Banjir 2021 46
## 9 Jawa Barat Tanah Longsor 2021 40
## 10 Jawa Tengah Banjir 2021 40
## # ℹ 107 more rows
#Catastrophe Intensity
years = length(unique(catastrophe$Tahun))
Km = length(catastrophe$Tahun)
lambda.m = Km/years
var.lmbda = function(years) {
total = 0
for (i in years) {
total = total + (aggregate(`Jumlah Korban` ~ `Tahun`, catastrophe, sum)$`Jumlah Korban`[i] - lambda.m)^2
}
varians = 1/(years - 1)*total
return(varians)
}
var.lamda = var.lmbda(years)
library(purrr)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked _by_ '.GlobalEnv':
##
## population
summary_cat = function(catastrophe, population) {
avg_accidents = catastrophe %>%
group_by(Provinsi) %>%
summarise(Average_Accident = length(`Jumlah Korban`)/years)
# Per 100 juta penduduk
summary_df = avg_accidents %>%
left_join(population, by = "Provinsi") %>%
mutate(`Per 100 juta penduduk` = (Average_Accident / Populasi) * 100)
# Variansi Provinsi
sum_prov = catastrophe %>%
group_by(Provinsi, Tahun) %>%
summarise(sum_tahun = length(`Jumlah Korban`))
merged_prov = left_join(sum_prov, avg_accidents, by = "Provinsi")
merged_prov = merged_prov %>%
group_by(Provinsi) %>%
summarise(`total` = sum((sum_tahun - Average_Accident)^2))
summary_df = summary_df %>%
left_join(merged_prov, by = "Provinsi") %>%
mutate(`Varians` = total/(years-1))
# Variansi/Mean
summary_df = summary_df %>%
mutate(`Var/Mean` = Varians/Average_Accident)
# Select and rename columns for final output
summary_df = summary_df %>%
select(Provinsi, `Rerata Kejadian (tahunan)` = Average_Accident, `Per 100 juta penduduk`, `Varians`, `Var/Mean`)
return(summary_df)
}
# Call the function
rangkuman = summary_cat(catastrophe, population)
## `summarise()` has grouped output by 'Provinsi'. You can override using the
## `.groups` argument.
nasional = data.frame(
"Provinsi" = "Indonesia",
"Rerata Kejadian (tahunan)" = lambda.m,
"Per 100 juta penduduk" = lambda.m*100/sum(population$Populasi),
`Varians`= var.lamda,
"Var/Mean" = var.lamda/lambda.m,
check.names = FALSE)
rangkuman = bind_rows(rangkuman, nasional)
rangkuman
## # A tibble: 32 × 5
## Provinsi Rerata Kejadian (tah…¹ Per 100 juta pendudu…² Varians `Var/Mean`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Aceh 0.333 0.00616 0.178 0.533
## 2 Bali 0.167 0.00377 0.139 0.833
## 3 Bangka Beli… 0.333 0.0223 0.178 0.533
## 4 Banten 0.833 0.00680 2.01 2.41
## 5 Bengkulu 0.5 0.0243 0.5 1
## 6 DI Yogyakar… 0.167 0.00443 0.139 0.833
## 7 DKI Jakarta 0.5 0.00468 0.5 1
## 8 Gorontalo 0.167 0.0140 0.139 0.833
## 9 Jawa Barat 2 0.00405 3.6 1.8
## 10 Jawa Tengah 2 0.00540 5.2 2.6
## # ℹ 22 more rows
## # ℹ abbreviated names: ¹`Rerata Kejadian (tahunan)`, ²`Per 100 juta penduduk`
library(eva)
mle_fit = gpdFit(catastrophe$`Jumlah Korban`, threshold = M-1/2, method = "mle")
mle_fit$par.ests
## Scale (Intercept) Shape (Intercept)
## 3.0698290 0.7920318
mle_fit
## Summary of fit:
## Estimate Std. Error z value Pr(>|z|)
## Scale (Intercept) 3.06983 0.53729 5.7135 1.1065e-08 ***
## Shape (Intercept) 0.79203 0.16567 4.7807 1.7470e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '*' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d.xk = log10(round(mean(catastrophe$`Jumlah Korban`), 0))
d.xk
## [1] 1.079181
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following objects are masked from 'package:eva':
##
## dgpd, pgev, pgpd, qgev, qgpd, rgpd
yk = rbetabinom(log10(catastrophe$`Jumlah Korban`)*0.1, catastrophe$`Jumlah Korban`, q, rho = 0)
# Insured Death Distribution
insured_death = function(yk, M){
Yk = numeric(length(yk))
for (i in 1:length(yk)) {
if(yk[i] >= M) {
Yk[i] = yk[i]
}
else {
Yk[i] = 0
}
}
return(Yk)
}
yk = insured_death(yk, M)
library(evmix)
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: gsl
##
## Attaching package: 'gsl'
## The following objects are masked from 'package:VGAM':
##
## erf, erfc, hzeta, zeta
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
##
## backsolve
##
## Attaching package: 'evmix'
## The following objects are masked from 'package:VGAM':
##
## dgpd, pgpd, qgpd, rgpd, rlplot
## The following objects are masked from 'package:eva':
##
## dgpd, pgpd, qgpd, rgpd
library(eva)
qqgpd = function(data, u, tau, xi){
excess = data[data > u]
m = length(excess)
prob = 1:m / (m + 1)
x.hat = u + tau / xi * ((1 - prob)^-xi - 1)
ylim <- xlim <- range(x.hat, excess)
plot(sort(excess), x.hat, xlab = "Sample Quantiles",
ylab = "Fitted Quantiles", xlim = xlim, ylim = ylim, main = "Q-Q Plot")
abline(0, 1, col = "grey", lwd = 2)
}
xk = catastrophe$`Jumlah Korban`
qqgpd(xk, M, mle_fit$par.ests[1], mle_fit$par.ests[2])
# Claim Distribution
zk = yk*R
claim_dist = function(zk, S, L){
Zk = numeric(length(zk))
for (i in 1:length(zk)) {
if(zk[i] < S) {
Zk[i] = 0
}
else {
if (S <= zk[i] && zk[i] < S + L) {
Zk[i] = zk[i] - S
}
else {
Zk[i] = L
}
}
}
return(Zk)
}
zk = claim_dist(zk, S, L)
# Monte Carlo Simulation
library(MonteCarlo)
## Loading required package: snow
N = 10000000
K = rpois(N, lambda.m)
Xk = rgpd(N, M, mle_fit$par.ests[1], mle_fit$par.ests[2])
d.Xk = log10(round(mean(Xk),0))
Y = rbetabinom(log10(round(Xk))*theta, round(Xk), prob = q)
Y = insured_death(Y, M)
Zk = Y*R
# Scenario A
ZK.a = claim_dist(Zk, 0, 10^11)
EC.a = mean(ZK.a)*mean(K)
sd.C.a = sqrt(lambda.m*mean(ZK.a)^2 + var(ZK.a)*lambda.m)
P.a = EC.a + qnorm(1-alpha)*sd.C.a
# Scenario B
ZK.b = claim_dist(Zk, S, L)
EC.b = mean(ZK.b)*mean(K)
sd.C.b = sqrt(lambda.m*mean(ZK.b)^2 + var(ZK.b)*lambda.m)
P.b = EC.b + qnorm(1-alpha)*sd.C.b
# Comparison Table
compare = data.frame(
"Variabel" = c("S", "L", "Expected Claim", "Std. Deviation of Claim", "Cat Reserve"),
"Scenario A" = c("L", "Inf", EC.a, sd.C.a, P.a),
"Scenario B" = c(S, L, EC.b, sd.C.b, P.b),
check.names = FALSE
)
compare
## Variabel Scenario A Scenario B
## 1 S L 300.0000
## 2 L Inf 6000.0000
## 3 Expected Claim 291.855757867997 102.6552
## 4 Std. Deviation of Claim 6579.65404682734 516.6553
## 5 Cat Reserve 4729.76497228921 451.1339
# Conclusion
conclusion = data.frame(
"Catastrophe Risk" = P.a,
"Re Cat Coverage With Premium" = P.b,
"Catastrophe Reserve" = P.a - P.b,
check.names = F
)
conclusion
## Catastrophe Risk Re Cat Coverage With Premium Catastrophe Reserve
## 1 4729.765 451.1339 4278.631
# Comparison Table 2
# Scenario A
ZK.a = claim_dist(Zk, L, 10^10)
EC.a = mean(ZK.a)*mean(K)
sd.C.a = sqrt(lambda.m*mean(ZK.a)^2 + var(ZK.a)*lambda.m)
P.a2 = EC.a + alpha*sd.C.a
# Scenario B
ZK.b = claim_dist(Zk, S, L)
EC.b = mean(ZK.b)*mean(K)
sd.C.b = sqrt(lambda.m*mean(ZK.b)^2 + var(ZK.b)*lambda.m)
P.b2 = EC.b + alpha*sd.C.b
# Comparison Table
compare = data.frame(
"Variabel" = c("S", "L", "Expected Claim", "Std. Deviation of Claim", "Cat Reserve"),
"Scenario A" = c("L", "Inf", EC.a, sd.C.a, P.a2),
"Scenario B" = c(S, L, EC.b, sd.C.b, P.b2),
check.names = FALSE
)
compare
## Variabel Scenario A Scenario B
## 1 S L 300.0000
## 2 L Inf 6000.0000
## 3 Expected Claim 68.6163697268939 102.6552
## 4 Std. Deviation of Claim 6490.15531700024 516.6553
## 5 Cat Reserve 1691.15519897695 231.8190
# Conclusion
conclusion2 = data.frame(
"Catastrophe Risk" = P.a2,
"Re Cat Coverage With Premium" = P.b2,
"Catastrophe Reserve" = P.a2 - P.b2,
check.names = F
)
conclusion2
## Catastrophe Risk Re Cat Coverage With Premium Catastrophe Reserve
## 1 1691.155 231.819 1459.336