Q1
library(readxl);library(ggplot2);library(dplyr)
## Warning: package 'ggplot2' was built under R version 3.5.1
## Warning: package 'dplyr' was built under R version 3.5.1
##
## 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
car.acc0<-read_excel("C:/Users/pj903/Documents/Assignment2/Assignment2-part1.xls",
sheet="Sheet1")
car.acc0$X<-as.numeric(car.acc0$X)
## Warning: 강제형변환에 의해 생성된 NA 입니다
colnames(car.acc0) <- c("X", "N")
car.acc<-car.acc0[1:50,]
##Q1 Histogram of X
ggplot(car.acc, aes(x=X, y=N))+geom_histogram(stat="identity")+
labs(title="자동차 사고 치료비 histogram", x="X(자동차 사고 치료비)", y="N(사고건수)")+
theme_bw()+theme(plot.title=element_text(face='bold',size=17))
## Warning: Ignoring unknown parameters: binwidth, bins, pad
Q2
##Q2 Q-Q plot
car.acc0 <- car.acc0 %>%
mutate(r=cumsum(N),
p=r/(sum(N)+1)) %>%
mutate(normal=qnorm(p),
weibull=log(-log(1-p)),
inverse.weibull=-log(-log(p)),
logistic=log(r/((sum(N)+1)-r)))
library(ggpmisc)
## Warning: package 'ggpmisc' was built under R version 3.5.1
## For news about 'ggpmisc', please, see https://www.r4photobiology.info/
## For on-line documentation see https://docs.r4photobiology.info/ggpmisc/
######
#Pareto_Method1
Pareto_m1<-function(lam){
lm<-lm(I(log(1+X/lam))~I(-log(1-p))+0, data=car.acc0)
summary(lm)$adj.r.squared
}
lambda.pa<-optim(par=5, fn=Pareto_m1, method="BFGS", control=list(fnscale=-1))$par
a.pa<-1/lm(log(1+X/lambda.pa)~I(-log(1-p))+0, data=car.acc0)$coefficients
car.acc0 <- car.acc0 %>%
mutate(pareto=-log(1-p))
#Lognormal Q-Q plot
lognormal.qq<-ggplot(car.acc0[1:50,],aes(normal,log(X)))+geom_point()+
geom_smooth(method="lm",se=F,formula = y ~ x)+
stat_poly_eq(formula = y ~ x,
aes(label = paste(..rr.label.., sep = "~~~")),
parse = TRUE,rr.digits = 4)+
labs(title="Lognormal Q-Q plot",x="inverse of G",y="X(r)")+
theme_bw()+theme(plot.title=element_text(face='bold'))
mu.lognormal<-summary( lm( log(X) ~ normal ,data=car.acc0))$coefficient[1]
sigma.lognormal<-summary( lm( log(X) ~ normal ,data=car.acc0))$coefficient[2]
#Pareto Q-Q plot
pareto.qq<-ggplot(car.acc0[1:50,],aes(pareto,log(1+X/lambda.pa)))+geom_point()+
geom_smooth(method="lm",se=F,formula = y ~ x)+
stat_poly_eq(formula = y ~ x,
aes(label = paste(..rr.label.., sep = "~~~")),
parse = TRUE,rr.digits = 4)+
labs(title="Pareto Q-Q plot",x="inverse of G",y="X(r)")+
theme_bw()+theme(plot.title=element_text(face='bold'))
#Weibull Q-Q plot
weibull.qq<-ggplot(car.acc0[1:50,],aes(weibull,log(X)))+geom_point()+geom_smooth(method="lm",se=F)+
labs(title="Weibull Q-Q plot",x="inverse of G",y="X(r)")+
stat_poly_eq(formula = y ~ x,
aes(label = paste(..rr.label.., sep = "~~~")),
parse = TRUE,rr.digits = 4)+
theme_bw()+theme(plot.title=element_text(face='bold'))
#Inverse Weibull Q-Q plot
inv.weibull.qq<-ggplot(car.acc0[1:50,],aes(inverse.weibull,log(X)))+geom_point()+geom_smooth(method="lm",se=F)+
labs(title="Frechet Q-Q plot",x="inverse of G",y="X(r)")+
stat_poly_eq(formula = y ~ x,
aes(label = paste(..rr.label.., sep = "~~~")),
parse = TRUE,rr.digits = 4)+
theme_bw()+theme(plot.title=element_text(face='bold'))
#Log-logistic Q-Q plot
loglogistic.qq<-ggplot(car.acc0[1:50,],aes(logistic,log(X)))+geom_point()+geom_smooth(method="lm",se=F)+
labs(title="Log-logistic Q-Q plot",x="inverse of G",y="X(r)")+
stat_poly_eq(formula = y ~ x,
aes(label = paste(..rr.label.., sep = "~~~")),
parse = TRUE,rr.digits = 4)+
theme_bw()+theme(plot.title=element_text(face='bold'))
#all of Q-Q plot
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(grid)
grid.arrange(lognormal.qq, pareto.qq,
weibull.qq, inv.weibull.qq,loglogistic.qq, nrow=3)
Q3
##MLE
library(FAdist)
#Pareto_Method1
Pareto_m1<-function(lam){
lm<-lm(I(log(1+X/lam))~I(-log(1-p))+0, data=car.acc0)
summary(lm)$adj.r.squared
}
lam.pa<-optim(par=5, fn=Pareto_m1, method="Nelder-Mead", control=list(fnscale=-1))$par
a.pa<-1/lm(log(1+X/lam.pa)~I(-log(1-p))+0, data=car.acc0)$coefficients
#Pareto
f.pareto<-function(par){
x1<-car.acc$X-1;x2<-car.acc$X;a<-par[1];lam<-par[2]
sum(log(1-(lam/(lam+x2))^a-(1-(lam/(lam+x1))^a))*car.acc$N[1:50])
+log(1-(1-(lam/(lam+50))^a))*(car.acc0$N[51])
}
optim4<-optim(par=c(10,10), fn=f.pareto, method="Nelder-Mead", control=list(fnscale=-1))$par
#frechet
f.frechet<-function(par){
x1<-car.acc$X-1;x2<-car.acc$X;c<-par[1];tau<-par[2]
sum(log(exp(-c/(x2^tau))-exp(-c/(x1^tau)))*car.acc0$N[1:50])+
log(1-exp(-c/(50^tau)))*(car.acc0$N[51])
}
optim5<-optim(par=c(5.41,0.59), fn=f.frechet, method="Nelder-Mead", control=list(fnscale=-1))$par
c.fr<-optim5[1];tau.fr<-optim5[2]
#pareto
xf.pareto<-function(x){
x*a.pa*(lam.pa^a.pa)*(lam.pa+x)^(-a.pa-1)
}
F.pareto<-function(x){
1-(lam.pa/(lam.pa+x))^a.pa
}
#frechet
xf.frechet<-function(x){
x*(c.fr*tau.fr*exp(-c.fr/x^(tau.fr))/x^(tau.fr+1))
}
F.frechet<-function(x){
exp(-c.fr/(x^tau.fr))
}
q3<-data.frame("A"=c(rep(0,5),rep(10,5),rep(20,5)),
"B"=rep(c(50,100,200,500,1000),3),
"보험료(pareto)"=rep(0,15),
"보험료(frechet)"=rep(0,15))
E.N<-sum(car.acc0$N)/271306
for (i in 1:nrow(q3)) {
A<-q3[i,1];B<-q3[i,2]
pareto.E.Y<-integrate(xf.pareto,lower=A,upper=A+B)$value+
B*(1-F.pareto(A+B))-A*(F.pareto(A+B)-F.pareto(A))
frechet.E.Y<-integrate(xf.frechet,lower=A,upper=A+B)$value+
B*(1-F.frechet(A+B))-A*(F.frechet(A+B)-F.frechet(A))
q3[i,3]<-round(pareto.E.Y*E.N*10000,-1)
q3[i,4]<-round(frechet.E.Y*E.N*10000,-1)
}
q3
## A B 보험료.pareto. 보험료.frechet.
## 1 0 50 9540 9500
## 2 0 100 14970 15020
## 3 0 200 22310 22910
## 4 0 500 35520 38310
## 5 0 1000 48770 55090
## 6 10 50 8110 8010
## 7 10 100 13170 13200
## 8 10 200 20220 20830
## 9 10 500 33170 35980
## 10 10 1000 46290 52620
## 11 20 50 7140 7080
## 12 20 100 11890 11990
## 13 20 200 18670 19380
## 14 20 500 31380 34290
## 15 20 1000 44380 50810
Q4
loss<-car.acc[1,1]/2*car.acc[1,2]
for (i in 2:nrow(car.acc)){
loss<-loss+((car.acc[i-1,1]+car.acc[i,1])/2)*car.acc[i,2]
}
total.loss<-loss*10000+500000*3541
q4<-data.frame("총보험료"=c(271306*q3[1,3],271306*q3[1,4]),
"손해율"=c(as.matrix(total.loss/(271306*q3[1,3])*100),
as.matrix(total.loss/(271306*q3[1,4])*100)))
rownames(q4)<-c("Pareto","Frechet")
q4
## 총보험료 손해율
## Pareto 2588259240 100.3373
## Frechet 2577407000 100.7598