This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

# Load necessary libraries
library(readxl)       # For reading Excel files
library(fitdistrplus) # For fitting distributions to data
## Loading required package: MASS
## Loading required package: survival
library(ggplot2)      # For data visualization
library(MASS)         # For statistical functions
library(actuar)       # For actuarial functions
## 
## Attaching package: 'actuar'
## The following objects are masked from 'package:stats':
## 
##     sd, var
## The following object is masked from 'package:grDevices':
## 
##     cm
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Read the data from the Excel file
Claims <- read_excel("C:/Users/regay/Downloads/Telegram Desktop/SIC Claims Summary 2nd Version.xlsx",col_types = c("text", "date", "date","numeric"))
# Step 1: Remove non-positive values
claims_cleaned <- Claims %>%
  filter(`Claims paid` > 0)

# Step 2: Summary statistics
summary(claims_cleaned$`Claims paid`)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     93.1   1588.5   4551.0  14018.2  12000.0 450000.0
# Step 3: Visual inspection with a histogram
hist(claims_cleaned$`Claims paid`, breaks = 30, main = "Claims Paid Distribution", xlab = "Claims Paid")

# Step 4: Boxplot to visually check for outliers
boxplot(claims_cleaned$`Claims paid`, main = "Boxplot of Claims Paid", ylab = "Claims Paid")

# Step 5: Detecting outliers using IQR method
Q1 <- quantile(claims_cleaned$`Claims paid`, 0.25)
Q3 <- quantile(claims_cleaned$`Claims paid`, 0.75)
IQR <- Q3 - Q1

# Define outlier bounds
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR

# Identify outliers
outliers <- claims_cleaned %>%
  filter(`Claims paid` < lower_bound | `Claims paid` > upper_bound)

# Remove outliers (if necessary)
claims_cleaned_no_outliers <- claims_cleaned %>%
  filter(`Claims paid` >= lower_bound & `Claims paid` <= upper_bound)

# Summary statistics after removing outliers
summary(claims_cleaned_no_outliers$`Claims paid`)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    93.06  1385.75  3614.70  5856.58  8000.00 27500.00
# Step 6: Visual inspection after removing outliers
hist(claims_cleaned_no_outliers$`Claims paid`, breaks = 30, main = "Claims Paid Distribution (No Outliers)", xlab = "Claims Paid")

# Step 7: Boxplot after removing outliers
boxplot(claims_cleaned_no_outliers$`Claims paid`, main = "Boxplot of Claims Paid (No Outliers)", ylab = "Claims Paid")

# Step 8: Visual inspection with a Q-Q plot
qqnorm(claims_cleaned_no_outliers$`Claims paid`, main = "Q-Q Plot of Claims Paid (No Outliers)")
qqline(claims_cleaned_no_outliers$`Claims paid`, col = "red")

attach(claims_cleaned_no_outliers)
claims <- as.data.frame(`Claims paid`)  # Extract the claim amounts
# Descriptive statistics
descdist(claims$`Claims paid`)       # Summary statistics of the claims data

## summary statistics
## ------
## min:  93.06   max:  27500 
## median:  3614.705 
## mean:  5856.578 
## estimated sd:  6078.877 
## estimated skewness:  1.535065 
## estimated kurtosis:  4.77776
# Plot histogram of claims data
ggplot(claims, aes(x = `Claims paid`)) +
  geom_histogram(aes(y = after_stat(density)),binwidth = 500, fill = "black", color = "blue", alpha = 0.7) +
  geom_density(color = "red", linewidth = 1) +
  labs(title = "Histogram of Claim Amounts", x = "Claim Amount", y = "Frequency")

# Fit different distributions to the data
claims$`Claims paid`[claims$`Claims paid`<=0] <- min(claims$`Claims paid`[claims$`Claims paid`>0])

fit_norm <- fitdist(claims$`Claims paid`, "norm")       # Fit normal distribution

fit_lognorm <- fitdist(claims$`Claims paid`, "lnorm")   # Fit lognormal distribution

# Method of moments estimates for starting values
mean_claims <- mean(claims$`Claims paid`)
var_claims <- var(claims$`Claims paid`)
shape_start <- mean_claims^2 / var_claims
rate_start <- mean_claims / var_claims

# Fit gamma distribution with method of moments starting values
fit_gamma <- fitdist(claims$`Claims paid`, "gamma", start = list(shape = shape_start, rate = rate_start), method = "mle", lower = 0)

fit_weibull <- fitdist(claims$`Claims paid`, "weibull") # Fit Weibull distribution

fit_exponential <- fitdist(claims$`Claims paid`, "exp",start = list(rate = 1/mean(claims$`Claims paid`)),method = "mme") # Fit exponential distribution

# Parameter estimation for each fitted distribution
fit_norm$estimate       # Estimated parameters for normal distribution
##     mean       sd 
## 5856.578 6076.279
fit_lognorm$estimate    # Estimated parameters for lognormal distribution
##  meanlog    sdlog 
## 8.092515 1.175169
fit_gamma$estimate      # Estimated parameters for gamma distribution
##        shape         rate 
## 0.9911789044 0.0001692133
fit_weibull$estimate    # Estimated parameters for Weibull distribution
##        shape        scale 
##    0.9729841 5784.3892083
fit_exponential$estimate # Estimated parameters for exponential distribution
##         rate 
## 0.0001707482
# Goodness-of-fit tests for each distribution
gofstat(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential),
        fitnames = c("Normal", "Lognormal", "Gamma", "Weibull", "Exponential"))
## Goodness-of-fit statistics
##                                  Normal  Lognormal      Gamma    Weibull
## Kolmogorov-Smirnov statistic  0.1722246 0.05443836 0.07221755 0.06397453
## Cramer-von Mises statistic   12.1409539 0.78023406 1.05202975 0.76870834
## Anderson-Darling statistic   71.3643618 5.11407239 6.77879458 5.53303186
##                              Exponential
## Kolmogorov-Smirnov statistic   0.0740763
## Cramer-von Mises statistic     1.1221339
## Anderson-Darling statistic     7.0985592
## 
## Goodness-of-fit criteria
##                                  Normal Lognormal    Gamma  Weibull Exponential
## Akaike's Information Criterion 23710.74  22638.50 22644.19 22642.78    22642.25
## Bayesian Information Criterion 23720.87  22648.63 22654.32 22652.91    22647.32
# Check summaries
summary(fit_norm)
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##      estimate Std. Error
## mean 5856.578   177.2417
## sd   6076.279   125.3288
## Loglikelihood:  -11853.37   AIC:  23710.74   BIC:  23720.87 
## Correlation matrix:
##      mean sd
## mean    1  0
## sd      0  1
summary(fit_lognorm)
## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters : 
##         estimate Std. Error
## meanlog 8.092515 0.03435637
## sdlog   1.175169 0.02429354
## Loglikelihood:  -11317.25   AIC:  22638.5   BIC:  22648.63 
## Correlation matrix:
##         meanlog sdlog
## meanlog       1     0
## sdlog         0     1
summary(fit_gamma)
## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters : 
##           estimate
## shape 0.9911789044
## rate  0.0001692133
## Loglikelihood:  -11320.1   AIC:  22644.19   BIC:  22654.32
summary(fit_weibull)
## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters : 
##           estimate   Std. Error
## shape    0.9729841   0.02198583
## scale 5784.3892083 184.89480535
## Loglikelihood:  -11319.39   AIC:  22642.78   BIC:  22652.91 
## Correlation matrix:
##           shape     scale
## shape 1.0000000 0.3264904
## scale 0.3264904 1.0000000
summary(fit_exponential)
## Fitting of the distribution ' exp ' by matching moments 
## Parameters : 
##          estimate   Std. Error
## rate 0.0001707482 0.0001707482
## Loglikelihood:  -11320.13   AIC:  22642.25   BIC:  22647.32
# Visualize the fitted distributions
plot.legend <- c("Normal", "Lognormal", "Gamma", "Weibull", "Exponential")

# Density plot
denscomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)

# CDF plot
cdfcomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)

# Q-Q plot
qqcomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)

# P-P plot
ppcomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)

plot.legend <- c("Normal", "Lognormal", "Gamma", "Weibull", "Exponential")
denscomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend) # Density plot

cdfcomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)  # CDF plot

qqcomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)  # Q-Q plot

ppcomp(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential), legendtext = plot.legend)  # P-P plot

# Select the best-fitting distribution based on goodness-of-fit statistics
gof_results <- gofstat(list(fit_norm, fit_lognorm, fit_gamma, fit_weibull, fit_exponential),
                       fitnames = c("Normal", "Lognormal", "Gamma", "Weibull", "Exponential"))
best_fit <- which.min(gof_results$aic) # Find the distribution with the minimum AIC
best_fit_distribution <- c("Normal", "Lognormal", "Gamma", "Weibull", "Exponential")[best_fit] # Get the name of the best-fitting distribution
print(paste("The best fitting distribution is:", best_fit_distribution)) # Print the best-fitting distribution
## [1] "The best fitting distribution is: Lognormal"