2.1 Frekuensi Klaim

summary(data$Frekuensi)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    34.0   195.2   387.0   806.8  1355.2  3277.0

sd(data$Frekuensi)

## [1] 865.4049

var(data$Frekuensi)

## [1] 748925.6

var(data$Frekuensi)/mean(data$Frekuensi)

## [1] 928.2284

2.2 Besaran Klaim

summary(data$Besaran)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 8.939e+08 4.970e+09 1.037e+10 1.833e+10 2.303e+10 1.615e+11

sd(data$Besaran)

## [1] 27923285243

var(data$Besaran)

## [1] 7.797099e+20

2.3 Visualisasi Data

par(mfrow=c(2,2))
plot(data$t, data$Frekuensi, type = "l", col = "blue",
     main = "Frekuensi Klaim Bulanan",
     xlab = "Waktu (Bulan)", ylab = "Frekuensi")
points(data$t, data$Frekuensi, pch = 19, cex = 0.6,col = "blue")

plot(data$t, data$Besaran, type = "l", col = "red",
     main = "Total Besaran Klaim Bulanan",
     xlab = "Waktu (Bulan)", ylab = "Besaran (Rp)")
points(data$t, data$Besaran, pch = 19, cex = 0.6, col = "red")

hist(data$Frekuensi, breaks = 20, col = "lightblue", probability = TRUE,
     main = "Distribusi Frekuensi Klaim",
     xlab = "Frekuensi")

hist(data$Besaran, breaks = 20, col = "lightcoral", probability = TRUE,
     main = "Distribusi Besaran Klaim",
     xlab = "Besaran (Rp)")

par(mfrow=c(1,1))

4.1 Regresi Poisson Non-Homogen

model_poisson <- glm(Frekuensi ~ t + dummy_1 + dummy_2 + dummy_3,
              family = poisson(link = "log"), data = data)
summary(model_poisson)

## 
## Call:
## glm(formula = Frekuensi ~ t + dummy_1 + dummy_2 + dummy_3, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.713515   0.064639  88.391   <2e-16 ***
## t            0.001580   0.001384   1.142    0.254    
## dummy_1      1.919350   0.057833  33.188   <2e-16 ***
## dummy_2     -0.046285   0.038652  -1.197    0.231    
## dummy_3      1.319292   0.024724  53.360   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 38545  on 47  degrees of freedom
## Residual deviance: 10716  on 43  degrees of freedom
## AIC: 11107
## 
## Number of Fisher Scoring iterations: 5

4.2 Uji Overdispersi

dispersiontest(model_poisson, trafo = 2)

## 
##  Overdispersion test
## 
## data:  model_poisson
## z = 3.8761, p-value = 5.307e-05
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
##     alpha 
## 0.1774416

deviance(model_poisson) / df.residual(model_poisson)

## [1] 249.2074

sum(residuals(model_poisson, type = "pearson")^2) / df.residual(model_poisson)

## [1] 235.3227

4.3 Regresi Binomial Negatif

nb_model <- glm.nb(Frekuensi ~ t + dummy_1 + dummy_2 + dummy_3, data = data)
summary(nb_model)

## 
## Call:
## glm.nb(formula = Frekuensi ~ t + dummy_1 + dummy_2 + dummy_3, 
##     data = data, init.theta = 2.415612934, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 5.343968   0.872563   6.124  9.1e-10 ***
## t           0.009773   0.018925   0.516 0.605548    
## dummy_1     2.241702   0.798955   2.806 0.005019 ** 
## dummy_2     0.138803   0.501538   0.277 0.781969    
## dummy_3     1.398393   0.362681   3.856 0.000115 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.4156) family taken to be 1)
## 
##     Null deviance: 137.286  on 47  degrees of freedom
## Residual deviance:  51.273  on 43  degrees of freedom
## AIC: 698.3
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.416 
##           Std. Err.:  0.467 
## 
##  2 x log-likelihood:  -686.300

4.4 Pemilihan Model Terbaik (AIC/BIC)

data.frame(
  Model = c("Poisson", "Negative Binomial"),
  AIC = c(AIC(model_poisson), AIC(nb_model)),
  BIC = c(BIC(model_poisson), BIC(nb_model))
)

##               Model        AIC        BIC
## 1           Poisson 11106.8578 11116.2138
## 2 Negative Binomial   698.2997   709.5269

lr_stat    <- 2 * (as.numeric(logLik(nb_model)) - as.numeric(logLik(model_poisson)))
p_value_lr <- pchisq(lr_stat, df = 1, lower.tail = FALSE) / 2
lr_stat

## [1] 10410.56

p_value_lr

## [1] 0

4.5 Diagnostik Model

par(mfrow = c(2, 2))

plot(fitted(nb_model), residuals(nb_model, type = "deviance"),
     xlab = "Fitted Values", ylab = "Deviance Residuals",
     main = "Residual vs Fitted - NB")
abline(h = 0, col = "red", lty = 2)

qqnorm(residuals(nb_model, type = "deviance"))
qqline(residuals(nb_model, type = "deviance"), col = "red")

plot(fitted(nb_model), sqrt(abs(residuals(nb_model, type = "deviance"))),
     xlab = "Fitted Values", ylab = "√|Deviance Residuals|",
     main = "Scale-Location - NB")

plot(cooks.distance(nb_model), type = "h",
     main = "Cook's Distance - NB", xlab = "Observation")
abline(h = 4/nrow(data), col = "red", lty = 2)

par(mfrow = c(1, 1))

pakai_nb <- AIC(nb_model) < AIC(model_poisson)
model_frekuensi <- if (pakai_nb) nb_model else model_poisson
data$fitted_freq <- fitted(model_frekuensi)
 
ggplot(data, aes(x = t)) +
  geom_point(aes(y = Frekuensi)) +
  geom_line(aes(y = fitted_freq), color = "red") +
  scale_x_continuous(breaks = data$t[seq(1, nrow(data), by = 3)],
                      labels = data$Bulan[seq(1, nrow(data), by = 3)]) +
  labs(title = "Frekuensi Klaim: Aktual vs Model Terpilih",
       x = "Bulan", y = "Jumlah Klaim") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

5.1 Uji Kruskal-Wallis

kw_test <- kruskal.test(Besaran ~ tahun, data = data)
kw_test

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Besaran by tahun
## Kruskal-Wallis chi-squared = 24.44, df = 3, p-value = 2.022e-05

5.2 Validasi Data Ekstrem

data %>%
  dplyr::filter(Besaran > quantile(Besaran, 0.90)) %>%
  dplyr::select(Bulan, Frekuensi, Besaran)

## # A tibble: 5 × 3
##   Bulan   Frekuensi       Besaran
##   <chr>       <int>         <dbl>
## 1 06-2020      2654  30208893540 
## 2 07-2020      3277  35165702210 
## 3 02-2023       234  42786417273.
## 4 10-2023       317 120323720454 
## 5 11-2023       426 161517828319

5.3 Pemilihan Threshold

persentil_uji <- c(0.60, 0.70, 0.80, 0.90)
hasil_threshold <- data.frame()

for (p in persentil_uji) {
  u_uji <- as.numeric(quantile(data$Besaran, p)) / 1e9
  fit_uji <- tryCatch(
    fitgpd(data$Besaran / 1e9, threshold = u_uji, est = "mle"),
    error = function(e) NULL,
    warning = function(w) NULL
  )
  if (!is.null(fit_uji)) {
    hasil_threshold <- rbind(hasil_threshold, data.frame(
      persentil = p,
      threshold_miliar = round(u_uji, 2),
      n_exceed  = fit_uji$nat,
      xi        = round(fit_uji$fitted.values["shape"], 3),
      se_xi     = round(fit_uji$std.err["shape"], 3)
    ))
  } else {
    hasil_threshold <- rbind(hasil_threshold, data.frame(
      persentil = p, threshold_miliar = round(u_uji, 2),
      n_exceed = NA, xi = NA, se_xi = NA
    ))
  }
}
print(hasil_threshold, row.names = FALSE)

##  persentil threshold_miliar n_exceed    xi se_xi
##        0.6            12.73       19 0.447 0.292
##        0.7            21.06       15 0.932 0.448
##        0.8            23.74       10 1.719 1.014
##        0.9            28.28       NA    NA    NA

mrlplot(data$Besaran)

u <- as.numeric(quantile(data$Besaran, 0.60))
u

## [1] 12728394866

5.4 Estimasi Parameter GPD

fit_gpd <- fitgpd(data$Besaran/1e9, threshold = u/1e9, est = "mle")
summary(fit_gpd)

## List of 23
##  $ fitted.values : Named num [1:2] 13.189 0.447
##   ..- attr(*, "names")= chr [1:2] "scale" "shape"
##  $ std.err       : Named num [1:2] 4.706 0.292
##   ..- attr(*, "names")= chr [1:2] "scale" "shape"
##  $ std.err.type  : chr "observed"
##  $ var.cov       : num [1:2, 1:2] 22.1493 -0.6753 -0.6753 0.0852
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "scale" "shape"
##   .. ..$ : chr [1:2] "scale" "shape"
##  $ fixed         : NULL
##  $ param         : Named num [1:2] 13.189 0.447
##   ..- attr(*, "names")= chr [1:2] "scale" "shape"
##  $ deviance      : num 153
##  $ corr          : NULL
##  $ convergence   : chr "successful"
##  $ counts        : Named int [1:2] 36 21
##   ..- attr(*, "names")= chr [1:2] "function" "gradient"
##  $ message       : NULL
##  $ threshold     : num 12.7
##  $ nat           : int 19
##  $ pat           : num 0.396
##  $ data          : num [1:48] 23.7 23.8 23.5 14.2 15.3 ...
##  $ exceed        : num [1:19] 23.7 23.8 23.5 14.2 15.3 ...
##  $ scale         : Named num 13.2
##   ..- attr(*, "names")= chr "scale"
##  $ var.thresh    : logi FALSE
##  $ est           : chr "MLE"
##  $ logLik        : num -76.5
##  $ opt.value     : num 76.5
##  $ hessian       : num [1:2, 1:2] 0.0595 0.4716 0.4716 15.468
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "scale" "shape"
##   .. ..$ : chr [1:2] "scale" "shape"
##  $ threshold.call: chr "12.728394866"
##  - attr(*, "class")= chr [1:2] "uvpot" "pot"

fit_gpd$nat

## [1] 19

fit_gpd$pat

## [1] 0.3958333

fit_gpd$fitted.values

##      scale      shape 
## 13.1894112  0.4467473

fit_gpd$std.err

##     scale     shape 
## 4.7063089 0.2919568

sigma <- as.numeric(fit_gpd$fitted.values["scale"]) * 1e9
sigma

## [1] 13189411202

xi    <- as.numeric(fit_gpd$fitted.values["shape"])
xi

## [1] 0.4467473

5.5 Diagnostik Model GPD

plot(fit_gpd, npy = 12)

5.6 Return Level GPD

return_level_gpd <- function(sigma, xi, u, pat, npy, T_years) {
  m <- T_years * npy
  u + (sigma / xi) * ((m * pat)^xi - 1)
}
sapply(c(1, 2, 5), function(T) return_level_gpd(sigma/1e9, xi, u/1e9, fit_gpd$pat, 12, T))

## [1]  42.42603  63.92105 104.75030

6.1 Skenario Normal

data_baru <- data.frame(t = nrow(data) + 1, dummy_1 = 0, dummy_2 = 0, dummy_3 = 0)
EN <- as.numeric(predict(model_frekuensi, newdata = data_baru, type = "response"))
 
VarN <- if (identical(model_frekuensi, nb_model)) {
  theta <- model_frekuensi$theta
  EN + (EN^2) / theta
} else {
  EN
}
EN

## [1] 337.9399

VarN

## [1] 47615.11

average_klaim <- mean(data$average)
ES <- EN * average_klaim
ES

## [1] 13735403599

6.2 Skenario Ekstrem

var_gpd <- function(p, u, sigma, xi, pat) {
  u + (sigma / xi) * (((1 - p) / pat)^(-xi) - 1)
}
 
VaR_95 <- var_gpd(0.95, u, sigma, xi, fit_gpd$pat)
VaR_95

## [1] 57607232681

VaR_99 <- var_gpd(0.99, u, sigma, xi, fit_gpd$pat)
VaR_99

## [1] 135908278330

Frekuensi Klaim

ringkasan_frekuensi <- data.frame(
  Komponen = c("Model terpilih", "AIC", "Rasio varians/mean (indikasi awal overdispersi)",
               "Proyeksi E(N) bulan berikutnya", "Proyeksi Var(N) bulan berikutnya"),
  Nilai = c(
    ifelse(pakai_nb, "Binomial Negatif", "Poisson Non-Homogen"),
    round(AIC(model_frekuensi), 2),
    round(var(data$Frekuensi) / mean(data$Frekuensi), 2),
    round(EN, 2),
    round(VarN, 2)
  )
)
ringkasan_frekuensi

##                                          Komponen            Nilai
## 1                                  Model terpilih Binomial Negatif
## 2                                             AIC            698.3
## 3 Rasio varians/mean (indikasi awal overdispersi)           928.23
## 4                  Proyeksi E(N) bulan berikutnya           337.94
## 5                Proyeksi Var(N) bulan berikutnya         47615.11

Besaran Klaim

ringkasan_besaran <- data.frame(
  Komponen = c("Threshold (u)", "Observasi di atas threshold",
               "Shape (xi)", "Scale (sigma)",
               "p-value Kruskal-Wallis (kestabilan antar tahun)"),
  Nilai = c(
    format(round(u, 0), big.mark = ","),
    paste0(fit_gpd$nat, " dari ", nrow(data)),
    round(xi, 4),
    format(round(sigma, 0), big.mark = ","),
    round(kw_test$p.value, 4)
  )
)
ringkasan_besaran

##                                          Komponen          Nilai
## 1                                   Threshold (u) 12,728,394,866
## 2                     Observasi di atas threshold     19 dari 48
## 3                                      Shape (xi)         0.4467
## 4                                   Scale (sigma) 13,189,411,202
## 5 p-value Kruskal-Wallis (kestabilan antar tahun)              0

Total Klaim

ringkasan_total <- data.frame(
  Ukuran = c("Estimasi skenario normal", "VaR 95% skenario ekstrem", "VaR 99% skenario ekstrem"),
  Nilai_Rupiah = c(
    format(round(ES, 0), big.mark = ","),
    format(round(VaR_95, 0), big.mark = ","),
    format(round(VaR_99, 0), big.mark = ",")
  )
)
ringkasan_total

##                     Ukuran    Nilai_Rupiah
## 1 Estimasi skenario normal  13,735,403,599
## 2 VaR 95% skenario ekstrem  57,607,232,681
## 3 VaR 99% skenario ekstrem 135,908,278,330