LAPORAN PKL

LIBRARY

library(dplyr)
library(ggplot2)
library(strucchange)
library(readxl)
library(AER)
library(MASS)
library(POT)

1. LOAD DATASET

data <- read_excel("DATASET BULANAN.xlsx")

data$Frekuensi <- as.integer(data$Frekuensi)
data$Besaran <- as.numeric(data$Besaran)
data$t <- 1:nrow(data)
data$tahun <- substr(data$Bulan, 4, 7)
data$average <- data$Besaran / data$Frekuensi

str(data)

## tibble [48 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Bulan    : chr [1:48] "01-2020" "02-2020" "03-2020" "04-2020" ...
##  $ Frekuensi: int [1:48] 2082 1874 1953 1185 1277 2654 3277 2310 2383 1835 ...
##  $ Besaran  : num [1:48] 2.37e+10 2.38e+10 2.35e+10 1.42e+10 1.53e+10 ...
##  $ t        : int [1:48] 1 2 3 4 5 6 7 8 9 10 ...
##  $ tahun    : chr [1:48] "2020" "2020" "2020" "2020" ...
##  $ average  : num [1:48] 11368970 12696007 12013940 11995723 11982443 ...

summary(data)

##     Bulan             Frekuensi         Besaran                t        
##  Length:48          Min.   :  34.0   Min.   :8.939e+08   Min.   : 1.00  
##  Class :character   1st Qu.: 195.2   1st Qu.:4.970e+09   1st Qu.:12.75  
##  Mode  :character   Median : 387.0   Median :1.037e+10   Median :24.50  
##                     Mean   : 806.8   Mean   :1.833e+10   Mean   :24.50  
##                     3rd Qu.:1355.2   3rd Qu.:2.303e+10   3rd Qu.:36.25  
##                     Max.   :3277.0   Max.   :1.615e+11   Max.   :48.00  
##     tahun              average         
##  Length:48          Min.   : 10649693  
##  Class :character   1st Qu.: 12705849  
##  Mode  :character   Median : 20612120  
##                     Mean   : 40644520  
##                     3rd Qu.: 28997660  
##                     Max.   :379570096

2. STATISTIK DESKRIPTIF

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))

3. UJI PERUBAHAN STRUKTURAL

bp <- breakpoints(Frekuensi ~ t, data = data)
summary(bp)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = Frekuensi ~ t, data = data)
## 
## Breakpoints at observation number:
##                       
## m = 1   10            
## m = 2   10       33   
## m = 3   10       33 40
## m = 4   7  14    33 40
## m = 5   7  14 25 33 40
## 
## Corresponding to breakdates:
##                                                                     
## m = 1   0.208333333333333                                           
## m = 2   0.208333333333333                                     0.6875
## m = 3   0.208333333333333                                     0.6875
## m = 4   0.145833333333333 0.291666666666667                   0.6875
## m = 5   0.145833333333333 0.291666666666667 0.520833333333333 0.6875
##                          
## m = 1                    
## m = 2                    
## m = 3   0.833333333333333
## m = 4   0.833333333333333
## m = 5   0.833333333333333
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS 2.814e+07 1.382e+07 8.399e+06 5.393e+06 4.582e+06 4.329e+06
## BIC 7.853e+02 7.628e+02 7.505e+02 7.409e+02 7.447e+02 7.536e+02

plot(bp)

breakdates(bp)

## [1] 0.2083333 0.6875000 0.8333333

data$dummy_1 <- ifelse(data$t <= 10, 1, 0)
data$dummy_2 <- ifelse(data$t > 10 & data$t <= 33, 1, 0)
data$dummy_3 <- ifelse(data$t > 33 & data$t <= 40, 1, 0)

4. PEMODELAN FREKUENSI KLAIM

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. PEMODELAN BESARAN KLAIM

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

mrlplot(data$Besaran)

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

## [1] 24570461014

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] 6.49 1.33
##   ..- attr(*, "names")= chr [1:2] "scale" "shape"
##  $ std.err       : Named num [1:2] 5.579 0.915
##   ..- attr(*, "names")= chr [1:2] "scale" "shape"
##  $ std.err.type  : chr "observed"
##  $ var.cov       : num [1:2, 1:2] 31.124 -2.905 -2.905 0.837
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "scale" "shape"
##   .. ..$ : chr [1:2] "scale" "shape"
##  $ fixed         : NULL
##  $ param         : Named num [1:2] 6.49 1.33
##   ..- attr(*, "names")= chr [1:2] "scale" "shape"
##  $ deviance      : num 67.2
##  $ corr          : NULL
##  $ convergence   : chr "successful"
##  $ counts        : Named int [1:2] 39 31
##   ..- attr(*, "names")= chr [1:2] "function" "gradient"
##  $ message       : NULL
##  $ threshold     : num 24.6
##  $ nat           : int 8
##  $ pat           : num 0.167
##  $ data          : num [1:48] 23.7 23.8 23.5 14.2 15.3 ...
##  $ exceed        : num [1:8] 30.2 35.2 24.6 27.5 26.3 ...
##  $ scale         : Named num 6.49
##   ..- attr(*, "names")= chr "scale"
##  $ var.thresh    : logi FALSE
##  $ est           : chr "MLE"
##  $ logLik        : num -33.6
##  $ opt.value     : num 33.6
##  $ hessian       : num [1:2, 1:2] 0.0475 0.165 0.165 1.7672
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "scale" "shape"
##   .. ..$ : chr [1:2] "scale" "shape"
##  $ threshold.call: chr "24.570461014"
##  - attr(*, "class")= chr [1:2] "uvpot" "pot"

fit_gpd$nat

## [1] 8

fit_gpd$pat

## [1] 0.1666667

fit_gpd$fitted.values

##   scale   shape 
## 6.48764 1.33307

fit_gpd$std.err

##     scale     shape 
## 5.5789076 0.9149306

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

## [1] 6487639746

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

## [1] 1.33307

5.5 Diagnostik Model GPD

plot(fit_gpd, npy = 12)

5.6 Return Level GPD

retlev(fit_gpd, period = c(1, 2, 5), npy = 12)

## Warning in plot.window(...): "period" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "period" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "period" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "period" is not a
## graphical parameter

## Warning in box(...): "period" is not a graphical parameter

## Warning in title(...): "period" is not a graphical parameter

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]  31.96482  50.59405 124.48985

6. TOTAL KLAIM

6.1 Hasil Model Frekuensi

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

6.2 Hasil Model Besaran

EX    <- u + sigma / (1 - xi)
VarX  <- (sigma^2) / (((1 - xi)^2) * (1 - 2*xi))
EX

## [1] 5092140303

VarX

## [1] -2.27715e+20

6.3 Estimasi Total Klaim

ES   <- EN * EX
VarS <- EN * VarX + VarN * (EX^2)
ES

## [1] 1.720837e+12

VarS

## [1] 1.157701e+24

sqrt(VarS)

## [1] 1.075965e+12

6.4 Simulasi Monte Carlo

set.seed(123)
n_sim <- 10000
 
N_sim <- if (identical(model_frekuensi, nb_model)) {
  rnbinom(n_sim, size = model_frekuensi$theta, mu = EN)
} else {
  rpois(n_sim, lambda = EN)
}
 
S_sim <- sapply(N_sim, function(n) {
  if (n == 0) return(0)
  sum(rgpd(n, loc = u, scale = sigma, shape = xi))
})
quantile(S_sim, 0.95)

##          95% 
## 6.825331e+14

quantile(S_sim, 0.99)

##          99% 
## 4.973398e+15

hist(S_sim, breaks = 50, main = "Distribusi Simulasi Total Klaim Bulanan",
     xlab = "Total Klaim (Rp)", col = "skyblue")
abline(v = quantile(S_sim, c(0.95, 0.99)), col = c("orange", "red"), lwd = 2)
legend("topright", legend = c("VaR 95%", "VaR 99%"),
       col = c("orange", "red"), lwd = 2)

HASIL AKHIR

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) 24,570,461,014
## 2                     Observasi di atas threshold      8 dari 48
## 3                                      Shape (xi)         1.3331
## 4                                   Scale (sigma)  6,487,639,746
## 5 p-value Kruskal-Wallis (kestabilan antar tahun)              0

Total Klaim

ringkasan_total <- data.frame(
  Ukuran = c("E(S) - analitik", "SD(S) - analitik",
             "VaR 95% - simulasi Monte Carlo", "VaR 99% - simulasi Monte Carlo"),
  Nilai_Rupiah = c(
    format(round(ES, 0), big.mark = ","),
    format(round(sqrt(VarS), 0), big.mark = ","),
    format(round(quantile(S_sim, 0.95), 0), big.mark = ","),
    format(round(quantile(S_sim, 0.99), 0), big.mark = ",")
  )
)
ringkasan_total

##                           Ukuran Nilai_Rupiah
## 1                E(S) - analitik 1.720837e+12
## 2               SD(S) - analitik 1.075965e+12
## 3 VaR 95% - simulasi Monte Carlo 6.825331e+14
## 4 VaR 99% - simulasi Monte Carlo 4.973398e+15