library(extRemes)
## Loading required package: Lmoments
## Loading required package: distillery
##
## Attaching package: 'extRemes'
## The following objects are masked from 'package:stats':
##
## qqnorm, qqplot
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(POT)
##
## Attaching package: 'POT'
## The following object is masked from 'package:extRemes':
##
## mrlplot
library(readxl)
library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
setwd("C:/Users/LENOVO/OneDrive/Documents/NATIONAL LIFE/Retensi")
data.ajn = read_excel("Data Klaim Meninggal 2020-2023.xlsx")
## New names:
## • `Approved Date` -> `Approved Date...19`
## • `Approved Date` -> `Approved Date...22`
data.ajn = data.ajn[order(data.ajn$`Claim Date`), ]
head(data.ajn, 1500)
## # A tibble: 1,500 × 24
## NO `JENIS PRODUK` `NOMOR KLAIM` NO.Polis `Certificate No`
## <dbl> <chr> <chr> <dbl> <dbl>
## 1 160 PJKN CLD200400000029 1.80e14 180204072190
## 2 2 PJKN CLD200200000006 1.90e14 190295237834
## 3 26 PJKN CLD200200000076 1.90e14 190295075607
## 4 8 PJKN CLD200100000097 1.90e14 190295250485
## 5 6 PJKN CLD200100000093 1.90e14 190295297584
## 6 7 PJKN CLD200100000094 1.90e14 190295275166
## 7 252 PJKN CLD200600000042 1.90e14 190295311582
## 8 319 PJKN CLD200600000176 1.80e14 180204116601
## 9 13 PJKN CLD200200000009 1.90e14 190295233417
## 10 31 PJKN CLD200200000082 1.90e14 190295186133
## # ℹ 1,490 more rows
## # ℹ 19 more variables: `Birth Date` <dttm>, `Start Date` <dttm>,
## # `End Date` <dttm>, `Claim Date` <dttm>, `Sum Insured ND` <dbl>,
## # `Usia Masuk` <dbl>, `Usia Meninggal (thn)` <dbl>, `usia Polis(BLN)` <dbl>,
## # Premium <dbl>, `Submission Claim` <dbl>, `Claim Approved Amount` <dbl>,
## # `Claim Paid Amount` <dbl>, `Submission Date` <dttm>,
## # `Approved Date...19` <dttm>, `Claim NAS` <dbl>, Status <chr>, …
plot(x = data.ajn$`Claim Date`, y = data.ajn$`Sum Insured ND`, type = "h",
xlab = "Years", ylab = "Sum Insured ND")
insured.max = aggregate(data.ajn$`Sum Insured ND`, by = list(as.numeric(substr(data.ajn$`Claim Date`, 1, 4))), max, na.rm = TRUE)
names(insured.max) = c("Year", "Max_Insured_Value")
max_insured_rows = data.frame()
for (i in 1:nrow(insured.max)) {
rows = data.ajn[data.ajn$`Sum Insured ND` == insured.max[i, "Max_Insured_Value"] &
substr(data.ajn$`Claim Date`, 1, 4) == insured.max[i, "Year"],
c("Claim Date", "Sum Insured ND")]
max_insured_rows <- rbind(max_insured_rows, rows)
}
print(max_insured_rows)
## # A tibble: 4 × 2
## `Claim Date` `Sum Insured ND`
## <dttm> <dbl>
## 1 2020-07-24 00:00:00 500000000
## 2 2021-02-13 00:00:00 500000000
## 3 2022-08-08 00:00:00 1500000000
## 4 2023-11-04 00:00:00 726644000
plot(data.ajn$`Claim Date`, data.ajn$`Sum Insured ND`, type = "h",
xlab = "Years", ylab = "Annual maxima loss amount", col = "grey")
points(max_insured_rows, pch = 15, cex = 1, col = "black")
library(evmix)
## Loading required package: splines
## Loading required package: gsl
## 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:POT':
##
## dgpd, mrlplot, pgpd, qgpd, rgpd, tcplot
## The following object is masked from 'package:extRemes':
##
## mrlplot
library(eva)
##
## Attaching package: 'eva'
## The following objects are masked from 'package:evmix':
##
## dgpd, pgpd, qgpd, rgpd
## The following objects are masked from 'package:POT':
##
## dgpd, pgpd, qgpd, rgpd
Amt = data.ajn$`Sum Insured ND`
a = 1/2*mean(Amt)*((mean(Amt))^2/(var(Amt))^2 + 1)
k = 0.5*((mean(Amt))^2/(var(Amt))^2 - 1)
u = a / (k + 1)
u
## [1] 53008246
data = data.ajn$`Sum Insured ND`
plotdist(data, histo = TRUE, demp = TRUE)
summary(data.ajn[,9:10])
## Claim Date Sum Insured ND
## Min. :2020-01-01 00:00:00 Min. :7.500e+05
## 1st Qu.:2020-12-22 00:00:00 1st Qu.:9.988e+06
## Median :2021-08-16 00:00:00 Median :2.800e+07
## Mean :2021-10-18 07:12:00 Mean :5.301e+07
## 3rd Qu.:2022-08-28 00:00:00 3rd Qu.:6.600e+07
## Max. :2023-12-13 00:00:00 Max. :1.500e+09
# Fitting distribusi Normal
fit_norm = fitdist(data, "norm")
summary(fit_norm)
## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 53008246 NA
## sd 72673396 NA
## Loglikelihood: -108924 AIC: 217851.9 BIC: 217865.2
## Correlation matrix:
## [1] NA
# Fitting distribusi Lognormal
fit_lnorm = fitdist(data, "lnorm")
summary(fit_lnorm)
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## meanlog 17.058995 0.01677279
## sdlog 1.252917 0.01186012
## Loglikelihood: -104365 AIC: 208734 BIC: 208747.3
## Correlation matrix:
## meanlog sdlog
## meanlog 1 0
## sdlog 0 1
# Fitting distribusi Weibull
fit_weib = fitdist(data, "weibull")
summary(fit_weib)
## Fitting of the distribution ' weibull ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 8.529559e-01 NA
## scale 5.121703e+07 NA
## Loglikelihood: -104670.8 AIC: 209345.5 BIC: 209358.8
## Correlation matrix:
## [1] NA
# Fitting distribusi Eksponensial
fit_exp = fitdist(data, "exp", method = c("mme"))
summary(fit_exp)
## Fitting of the distribution ' exp ' by matching moments
## Parameters :
## estimate
## rate 1.886499e-08
## Loglikelihood: -104825.6 AIC: 209653.3 BIC: 209659.9
data.frame("Distribution" = c("Normal", "Log-Normal", "Weibull", "Exponential"),
"AIC" = c(fit_norm$aic, fit_lnorm$aic, fit_weib$aic, fit_exp$aic),
"BIC" = c(fit_norm$bic, fit_lnorm$bic, fit_weib$bic, fit_exp$bic),
"Log-Likelihood" = c(fit_norm$loglik, fit_lnorm$loglik, fit_weib$loglik, fit_exp$loglik),
check.names=FALSE)
## Distribution AIC BIC Log-Likelihood
## 1 Normal 217851.9 217865.2 -108924.0
## 2 Log-Normal 208734.0 208747.3 -104365.0
## 3 Weibull 209345.5 209358.8 -104670.8
## 4 Exponential 209653.3 209659.9 -104825.6
#Melihat parameter distribusi data
ks.test(data, "pnorm", mean = fit_norm$estimate[1], sd = fit_norm$estimate[2])
## Warning in ks.test.default(data, "pnorm", mean = fit_norm$estimate[1], sd =
## fit_norm$estimate[2]): ties should not be present for the one-sample
## Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: data
## D = 0.2381, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(data, "plnorm", meanlog = fit_lnorm$estimate[1], sdlog = fit_lnorm$estimate[2])
## Warning in ks.test.default(data, "plnorm", meanlog = fit_lnorm$estimate[1], :
## ties should not be present for the one-sample Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: data
## D = 0.064915, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(data, "pweibull", shape = fit_weib$estimate[1], scale = fit_weib$estimate[2])
## Warning in ks.test.default(data, "pweibull", shape = fit_weib$estimate[1], :
## ties should not be present for the one-sample Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: data
## D = 0.11845, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(data, "pexp", fit_exp$estimate)
## Warning in ks.test.default(data, "pexp", fit_exp$estimate): ties should not be
## present for the one-sample Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: data
## D = 0.15563, p-value < 2.2e-16
## alternative hypothesis: two-sided
Value at Risk
alpha = 0.05 #umumnya 30% untuk premi
p = 0.3#pad (loading factor)
b = exp(fit_lnorm$estimate[1] + ((fit_lnorm$estimate[2])^2)/2)
pstar = 1/(1 + p)
sx = function(x)
{1-pnorm((log(x) - fit_lnorm$estimate[1])/fit_lnorm$estimate[2])}
sxinversalpha = qlnorm(1-alpha, meanlog = fit_lnorm$estimate[1], sdlog = fit_lnorm$estimate[2])
sxinverspstar = qlnorm(1-pstar, meanlog = fit_lnorm$estimate[1], sdlog = fit_lnorm$estimate[2])
retensi = sxinverspstar
sx1 = function(x) 1-pnorm((log(x) - fit_lnorm$estimate[1])/fit_lnorm$estimate[2])
for(a in 10^(7:12)) {
pid = integrate(sx1, a, Inf, stop.on.error = FALSE)
if (pid$message == "OK") break
}
pid = integrate(sx1, retensi, a)$value
deltad=(1+p)*pid
VaR1=retensi+deltad
VaR2=sxinversalpha+deltad
if ( alpha < pstar && pstar < sx(0)
&& sxinversalpha >= (1+p)*b
&& sxinversalpha >= sxinverspstar + deltad)
{ cat("Retensi ada.","\n")
cat("Retensi = ", retensi, "\n")
if ( 0 < retensi && retensi <= sxinversalpha){
cat("Value at Risk 1 = ",VaR1,"\n")
}else
{
cat("Value at Risk 2 = ",VaR2,"\n")
}
}else{cat("Retensi tidak ada.","\n")}
## Retensi ada.
## Retensi = 10185257
## Value at Risk 1 = 71317047
Conditional Tail Method
f=function(x)
{1-pnorm((log(x)-fit_lnorm$estimate[1])/fit_lnorm$estimate[2])}
L = integrate(f,sxinversalpha,retensi)$value
add = L/alpha
CTE1 = retensi+deltad
CTE2 = sxinversalpha+deltad+add
if ( alpha <= pstar && pstar < sx(0))
{ cat("Retensi ada.","\n")
cat("Retensi = ", retensi, "\n")
if ( 0 < retensi && retensi <= sxinversalpha){
cat("Conditional Tail Expectation 1 = ",CTE1,"\n")
}else
{
cat("Conditional Tail Expectation 2 = ",CTE2,"\n")
}
}else {
cat("Syarat tidak terpenuhi. \n")}
## Retensi ada.
## Retensi = 10185257
## Conditional Tail Expectation 1 = 71317047
insured.frame = data.frame("Date" = data.ajn$`Claim Date`, "Insured" = data.ajn$`Sum Insured ND`)
plot(insured.frame, type = "l", col = "black", xlab = "Years", ylab = "Loss exceedance amount")
clip(x1 = min(insured.frame$Date),
x2 = max(insured.frame$Date),
y1 = min(insured.frame$Insured),
y2 = retensi)
lines(insured.frame, lwd = 1, col = "grey")
abline(h = retensi, col="red")
