非寿险精算实务部分例题的程序代码.docVIP

非寿险精算实务部分例题的程序代码 (2013-03-09 16:02:56) 转载▼ 分类： 件 下面是2013年3月9-10日非寿险精算实务课程上部分例题的程序代码： expo=c(8000,13600,400,5200,6000,800,2000,2400,1600) loss=c(8400000825000,7020000,7200000,1140000,3600000,4590000,2640000) place=rep(c(A,B,C),each=3) type=as.factor(rep(c(1:3),3)) dat=data.frame(place,type,expo,loss) dat$off=log(expo) dat$pure=loss/expo dat #用线性回归模型 mod=lm(pure~place+type,data=dat,weights=expo) summary(mod) pred=predict(mod,type=response) matrix(pred,3,3) plot(dat$pure,pred,pch=16,cex=1.5);abline(0,1,col=2) #假设loss表示损失次数 m1=glm(loss~place+type,data=dat,family=poisson(link=log),offset=off) summary(m1) pred1=predict(m1,type=response) matrix(pred1,3,3);matrix(loss,3,3) sum(loss);sum(pred1) plot(pred1,loss,pch=16,cex=1.5);abline(0,1,col=2) pred11=pred1/dat$expo plot(pred11,dat$pure);abline(0,1,col=2) matrix(pred11,3,3) #与边际总和法等价 #假设loss表示索赔次数，pure表示索赔频率 m3=glm(pure~place+type,data=dat,family=poisson(link=log),weights=expo) summary(m3) pred3=predict(m3,type=response) pred3 plot(pred3,dat$pure,pch=16,cex=1.5);abline(0,1,col=2) matrix(pred3,3,3) #与边际总和法等价 #假设loss表示总赔款，pure表示索赔强度 m2=glm(pure~place+type,data=dat,family=Gamma(link=log),weights=expo) summary(m2) pred2=predict(m2,type=response) plot(pred2,dat$pure,pch=16,cex=1.5);abline(0,1,col=2) matrix(pred2,3,3);matrix(dat$pure,3,3) sum(dat$pure);sum(pred2) #结果比较 com=data.frame(code=1:9,正态=pred,伽马=pred2,泊松=pred3) com library(ggplot2) library(reshape) com1=melt(com,id=code) com1 qplot(code,value,data=com1,geom=c(line,point),group=variable,colour=variable) #续保率模型 set.seed(111) dat$ratio=runif(9) dat$expo1=round(dat$expo*dat$ratio,0) dat m4=glm(cbind(expo1,expo)~place+type,data=dat,family=binomial(link=logit)) summary(m4) pred4=predict(m4,type=response) pred4 plot(pred4,dat$ratio,pch=16,cex=1.5);abline(0,1,col=2) ============================================================================= 信度模型 library(actuar) dat=data.frame(hachemeister) dat #绘制观察数据 library(reshape) library(ggplot2) dat0=subset(dat,select=1:13)