Actuariat IARD

 

Je mets une copie brute des codes utilisés en classe. Désolé de ne pas avoir pris le temps de faire de remise en forme, j’essayerais de faire des billets, à l’occasion, si des questions précises sont posées.

  • Cours 3
baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",
header=TRUE,sep=",")
tail(baseavocat)
baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",")
avocat=(baseavocat$ATTORNEY==1)
sexe=rep("H",length(avocat)); sexe[baseavocat$CLMSEX==2]="F"
marital=rep("M",length(avocat)); marital[baseavocat$MARITAL==2]="C"
marital[baseavocat$MARITAL==3]="V"; marital[baseavocat$MARITAL==4]="D"
ceinture=(baseavocat$SEATBELT==1);
age=baseavocat$CLMAGE
cout=baseavocat$LOSS
base=data.frame(avocat,sexe,marital,ceinture,age,cout)
logit = function(formula, lien="logit", data=NULL) {
glm(formula,family=binomial(link=lien),data) }
exp(-0.87208+0.31670*10)/
(1+exp(-0.87208+0.31670*10))
summary(logit(avocat~cout,data=base))
summary(logit(avocat~sexe,data=base))
# femme
exp(-0.10252+0.33093*0)/
(1+exp(-0.10252+0.33093*0))
# homme
exp(-0.10252+0.33093*1)/
(1+exp(-0.10252+0.33093*1))
table(base$avocat,base$sexe)
log(352/(390+352) / (390/(390+352))
log(352/390)
log(333/265)-log(352/390)
reg=logit(avocat~sexe+cout,data=base)
summary(reg)
newindiv=data.frame(sexe="H",cout=10)
predict(reg,newdata=newindiv)
exp( -1.05051+0.38152*1+0.31994*10)/
(1+exp( -1.05051+0.38152*1+0.31994*10))
predict(reg,newdata=newindiv,type="response")
newindivH=data.frame(sexe="H",cout=seq(0,15))
newindivF=data.frame(sexe="F",cout=seq(0,15))
predH=predict(reg,newdata=newindivH,type="response")
predF=predict(reg,newdata=newindivF,type="response")
plot(seq(0,15),predH,col="blue",type="b",lwd=3)
lines(seq(0,15),predF,col="red",type="b",lwd=3)
reg=logit(avocat~.+cout,data=base)
base=read.table("F://baselogit-tree.txt",sep=";",
header=TRUE)
attach(base)
plot(X1,X2,col=Y+1)
library(tree)
arbre=tree(Y~X1+X2,data=base)
reg=logit(Y~X1+X2,data=base)
predict(reg,newdata=data.frame(X1=20,X2=6),type="response")
plot(arbre)
text(arbre)
arbre=tree(avocat~log(cout)+age,data=base)
  • Cours 4
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
header=TRUE,sep=";")
sinistres=sinistre[sinistre$garantie=="1RC",]
sinistres=sinistres[sinistres$cout>0,]
contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt",
header=TRUE,sep=";")
T=table(sinistres$nocontrat)
T1=as.numeric(names(T))
T2=as.numeric(T)
nombre1 = data.frame(nocontrat=T1,nbre=T2)
I = contrat$nocontrat%in%T1
T1= contrat$nocontrat[I==FALSE]
nombre2 = data.frame(nocontrat=T1,nbre=0)
nombre=rbind(nombre1,nombre2)
baseFREQ = merge(contrat,nombre)
Y=baseFREQ$nbre>0
X1=baseFREQ$ageconducteur 
M=table(Y,X1)
taux=M[2,]/apply(M,2,sum)
plot(as.numeric(colnames(M)),taux,pch=19)
X2=baseFREQ$agevehicule
db=data.frame(Y,X1,X2)
library(tree)
T=tree(Y~X1+X2,split='gini',mincut=5000,data=db)
partition.tree(T)
plot(T)
text(T)
y=predict(T,newdata=data.frame(X1=18:100))
lines(18:100,y,type='s',col='red',lwd=5)
          ################################
Y=baseFREQ$nbre
X1=baseFREQ$carburant
X2=cut(baseFREQ$ageconducteur,c(0,26.5,60.5,100))
E=baseFREQ$exposition
db=data.frame(Y,X1,X2,E)
head(db)
X12=paste(X1,X2)
table(Y,X1)
tapply(Y,X1,sum)
tapply(E,X1,sum)
tapply(Y,X1,sum)/tapply(E,X1,sum)
tapply(Y,X2,sum)
tapply(E,X2,sum)
tapply(Y,X2,sum)/tapply(E,X2,sum)
tapply(Y,X12,sum)
tapply(E,X12,sum)
tapply(Y,X12,sum)/tapply(E,X12,sum)
M=matrix(tapply(Y,X12,sum)/tapply(E,X12,sum),3,2)
colnames(M)=c("D","E")
rownames(M)=c("J","A","S")
  • Cours 5
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
                    header=TRUE,sep=";")
sinistres=sinistre[sinistre$garantie=="1RC",]
sinistres=sinistres[sinistres$cout>0,]
contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt",
                   header=TRUE,sep=";")
T=table(sinistres$nocontrat)
T1=as.numeric(names(T))
T2=as.numeric(T)
nombre1 = data.frame(nocontrat=T1,nbre=T2)
I = contrat$nocontrat%in%T1
T1= contrat$nocontrat[I==FALSE]
nombre2 = data.frame(nocontrat=T1,nbre=0)
nombre=rbind(nombre1,nombre2)
baseFREQ = merge(contrat,nombre)

reg=glm(nbre~s(ageconducteur)+
          offset(log(exposition)),
        family=poisson,data=baseFREQ)

nouveau=
  data.frame(ageconducteur=
               seq(15,120),exposition=1)
lambda=predict(
  reg,newdata=nouveau,
  type="response")

plot(nouveau$age,lambda,type="l",
     col="purple",lwd=3)
abline(h=sum(baseFREQ$nbre)/
         sum(baseFREQ$exposition),col="red")

moyenne=sum(baseFREQ$nbre)/
  sum(baseFREQ$exposition)

plot(nouveau$age,lambda/moyenne*100,type="l",
     col="purple",lwd=3)
abline(h=100,col="red")

reg2=glm(nbre~as.factor(ageconducteur)+
           offset(log(exposition)),
         family=poisson,data=baseFREQ)

nouveau=
  data.frame(ageconducteur=
               seq(18,95),exposition=1)
lambda2=predict(
  reg2,newdata=nouveau,
  type="response")

plot(nouveau$age,lambda2,type="p",
     col="purple",lwd=3,pch=19)
abline(h=sum(baseFREQ$nbre)/
         sum(baseFREQ$exposition),col="red")

lambda=predict(
  reg,newdata=nouveau,
  type="response")
lines(nouveau$age,lambda,col="blue",lwd=3)
summary(reg)

seuils=c(0,25,40,70,120)
reg3=glm(nbre~cut(ageconducteur,seuils)+
           offset(log(exposition)),
         family=poisson,data=baseFREQ)

lambda3=predict(
  reg3,newdata=nouveau,
  type="response",se.fit = TRUE)

lines(nouveau$age,lambda3$fit,col="orange",lwd=2,type="s")
lines(nouveau$age,lambda3$fit+
        2*lambda3$se.fit,col="orange",type="s")
lines(nouveau$age,lambda3$fit-
        2*lambda3$se.fit,col="orange",type="s")
barplot(tapply(baseFREQ$exposition,baseFREQ$ageconducteur,sum))

seuils=c(0,19,23,26,41,65,120)
reg4=glm(nbre~cut(ageconducteur,seuils)+
           offset(log(exposition)),
         family=poisson,data=baseFREQ)

lambda4=predict(
  reg4,newdata=nouveau,
  type="response",se.fit = TRUE)

lines(nouveau$age,lambda4$fit,col="green",lwd=2,type="s")
lines(nouveau$age,lambda4$fit+
        2*lambda4$se.fit,col="green",type="s")
lines(nouveau$age,lambda4$fit-
        2*lambda4$se.fit,col="green",type="s")

library(car)
linearHypothesis(reg4,c( 
"cut(ageconducteur, seuils)(26,41] = cut(ageconducteur, seuils)(41,65]",
"cut(ageconducteur, seuils)(41,65] = cut(ageconducteur, seuils)(65,120]"))

jeune=26
reg5=glm(nbre~ageconducteur+I(ageconducteur<jeune)+
           offset(log(exposition)),
         family=poisson,data=baseFREQ)

lambda5=predict(
  reg5,newdata=nouveau,
  type="response",se.fit = TRUE)

lines(nouveau$age,lambda5$fit,col="red",lwd=6,type="s")
  • Cours 6
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
                    header=TRUE,sep=";")
sinistres=sinistre[sinistre$garantie=="1RC",]
sinistres=sinistres[sinistres$cout>0,]
contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt",
                   header=TRUE,sep=";")
T=table(sinistres$nocontrat)
T1=as.numeric(names(T))
T2=as.numeric(T)
nombre1 = data.frame(nocontrat=T1,nbre=T2)
I = contrat$nocontrat%in%T1
T1= contrat$nocontrat[I==FALSE]
nombre2 = data.frame(nocontrat=T1,nbre=0)
nombre=rbind(nombre1,nombre2)
baseFREQ = merge(contrat,nombre)

reg2=glm(nbre~as.factor(ageconducteur)+
           offset(log(exposition)),
         family=poisson,data=baseFREQ)

nouveau=
  data.frame(ageconducteur=
               seq(18,95),exposition=1)
lambda2=predict(
  reg2,newdata=nouveau,
  type="response")

plot(nouveau$age,lambda2,type="b",
     col="purple",cex=1,pch=19)
abline(h=sum(baseFREQ$nbre)/
         sum(baseFREQ$exposition),col="red")

regpoly=glm(nbre~poly(ageconducteur,degree=1)+
           offset(log(exposition)),
         family=poisson,data=baseFREQ)
lambdapoly=predict(
  regpoly,newdata=nouveau,
  type="response")
lines(nouveau$age,lambdapoly,col="blue")

regpoly=glm(nbre~poly(ageconducteur,degree=20)+
              offset(log(exposition)),
            family=poisson,data=baseFREQ)
lambdapoly=predict(
  regpoly,newdata=nouveau,
  type="response")
lines(nouveau$age,lambdapoly,col="red",lwd=2)

attach(baseFREQ)
library(splines)
B=bs(ageconducteur,degree=1,knots=c(40))

regsp=glm(nbre~bs(ageconducteur,degree=2,
                  knots=c(30))+
              offset(log(exposition)),
            family=poisson,data=baseFREQ)
lambdasp=predict(
  regsp,newdata=nouveau,
  type="response")
lines(nouveau$age,lambdasp,col="blue",lwd=2)

B=bs(ageconducteur,degree=2,
knots=c(30,70))
matplot(ageconducteur,B)

summary(regsp)

prevloc = function(x){
I=abs(baseFREQ$ageconducteur-x)<=5
regloc=glm(nbre~ageconducteur+
            offset(log(exposition)),
            family=poisson,data=baseFREQ[I,])

prevloc = function(x){
  I=dnorm(baseFREQ$ageconducteur,mean=x,sd=3)
  regloc=glm(nbre~ageconducteur+
               offset(log(exposition)),
             family=poisson,data=baseFREQ,
             weights=I)

return(
  predict(regloc,
  newdata=data.frame(exposition=1,
     ageconducteur=x),type="response"))
}

lines(20:80,Vectorize(prevloc)(20:80),lwd=3,col="red")

library(mgcv)

regsp=gam(nbre~s(ageconducteur)+
            offset(log(exposition)),
          family=poisson,data=baseFREQ)

regsp0=gam(nbre~ageconducteur+agevehicule+
             offset(log(exposition)),
           family=poisson,data=baseFREQ)
regsp1=gam(nbre~s(ageconducteur)+s(agevehicule)+
            offset(log(exposition)),
          family=poisson,data=baseFREQ)
regsp2=gam(nbre~s(ageconducteur,agevehicule)+
             offset(log(exposition)),
           family=poisson,data=baseFREQ)

prev0 = function(agec,agev){
predict(regsp0,newdata=data.frame(
  ageconducteur=agec,
  agevehicule=agev,
  exposition=1),type="response")  
}
prev1 = function(agec,agev){
  predict(regsp1,newdata=data.frame(
    ageconducteur=agec,
    agevehicule=agev,
    exposition=1),type="response")  
}
prev2 = function(agec,agev){
  predict(regsp2,newdata=data.frame(
    ageconducteur=agec,
    agevehicule=agev,
    exposition=1),type="response")  
}

x1=seq(18,60,by=2)
y1=seq(0,15)
z0=outer(x1,y1,prev0)
z1=outer(x1,y1,prev1)
z2=outer(x1,y1,prev2)
persp(x1,y1,z1,col="green",
      ticktype="detailed",shade=TRUE,theta=30)
image(x1,y1,z1)
contour(x1,y1,z1,add=TRUE)

persp(x1,y1,z2,col="green",
      ticktype="detailed",shade=TRUE,theta=30)
image(x1,y1,z2)
contour(x1,y1,z2,add=TRUE)
  •  cours 7
y=c(1,2,4,2,6)
x=1:5
db=data.frame(x,y)
plot(db,pch=19,xlim=c(0,6),ylim=c(0,8))
vx=seq(.5,5.5,by=.1)
model1 = glm(y~x,family=gaussian(link="identity"))
y1=predict(model1,newdata=data.frame(x=vx))
lines(vx,y1,col="red")

model2 = glm(y~x,family=gaussian(link="log"))
y2=predict(model2,newdata=data.frame(x=vx),
           type="response")
lines(vx,y2,col="blue")

model3 = glm(y~x,family=gaussian(link="inverse"))
y3=predict(model3,newdata=data.frame(x=vx),
           type="response")
lines(vx,y3,col="purple")

model4 = glm(y~x,family=poisson(link="identity"))
y4=predict(model4,newdata=data.frame(x=vx),
           type="response")
lines(vx,y4,col="blue")
y1=predict(model1,newdata=data.frame(x=vx),se=TRUE)
lines(vx,y1$fit,col="red")
lines(vx,y1$fit+2*y1$se.fit,col="red",lty=2)
lines(vx,y1$fit-2*y1$se.fit,col="red",lty=2)
y4=predict(model4,newdata=data.frame(x=vx),se=TRUE)
lines(vx,y4$fit,col="blue")
lines(vx,y4$fit+2*y4$se.fit,col="blue",lty=2)
lines(vx,y4$fit-2*y4$se.fit,col="blue",lty=2)
y=c(2,2,4,2,6)
x=1:5
db=data.frame(x,y)
plot(db,pch=19,xlim=c(0,6),ylim=c(0,8))
model1 = glm(y~x,family=gaussian(link="identity"))
model4 = glm(y~x,family=poisson(link="identity"))
y1=predict(model1,newdata=data.frame(x=vx),se=TRUE)
lines(vx,y1$fit,col="red")
lines(vx,y1$fit+2*y1$se.fit,col="red",lty=2)
lines(vx,y1$fit-2*y1$se.fit,col="red",lty=2)
y4=predict(model4,newdata=data.frame(x=vx),se=TRUE)
lines(vx,y4$fit,col="blue")
lines(vx,y4$fit+2*y4$se.fit,col="blue",lty=2)
lines(vx,y4$fit-2*y4$se.fit,col="blue",lty=2)
y=c(1,2,4,2,6)
x=1:5
db=data.frame(x,y)
plot(db,pch=19,xlim=c(0,6),ylim=c(0,8))
model1 = glm(y~x,family=gaussian(link="identity"),data=db)
model2 = glm(y~x,family=gaussian(link="log"),data=db)
model5 = glm(y~x,family=poisson(link="log"),data=db)
lines(vx,y1$fit,col="red")
y2=predict(model2,newdata=data.frame(x=vx),se=TRUE,type="response")
lines(vx,y2$fit,col="purple")
lines(vx,y2$fit+2*y2$se.fit,col="purple",lty=2)
lines(vx,y2$fit-2*y2$se.fit,col="purple",lty=2)
y5=predict(model5,newdata=data.frame(x=vx),se=TRUE,type="response")
lines(vx,y5$fit,col="blue")
lines(vx,y5$fit+2*y5$se.fit,col="blue",lty=2)
lines(vx,y5$fit-2*y5$se.fit,col="blue",lty=2)
model0=lm(y~x,data=db)
summary(model0)
library(MASS)
model5 = glm(y~x,family=poisson(link="log"),
             data=db)
model6 = glm.nb(y~x,
             data=db)
  • Cours 8
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
                    header=TRUE,sep=";")
sinistres=sinistre[sinistre$garantie=="1RC",]
sinistres=sinistres[sinistres$cout>0,]
contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt",
                   header=TRUE,sep=";")
T=table(sinistres$nocontrat)
T1=as.numeric(names(T))
T2=as.numeric(T)
nombre1 = data.frame(nocontrat=T1,nbre=T2)
I = contrat$nocontrat%in%T1
T1= contrat$nocontrat[I==FALSE]
nombre2 = data.frame(nocontrat=T1,nbre=0)
nombre=rbind(nombre1,nombre2)
baseFREQ = merge(contrat,nombre)

reg=glm(nbre~carburant+zone+ageconducteur+
          offset(log(exposition)),
        data=baseFREQ,family=poisson(link="log"))

library(MASS)
reg2=glm.nb(nbre~carburant+zone+ageconducteur+
          offset(log(exposition)),
        data=baseFREQ)

library(AER)
dispersiontest(reg)

regq=glm(nbre~carburant+zone+ageconducteur+
          offset(log(exposition)),
        data=baseFREQ,family=quasipoisson(link="log"))

library(pscl)
regzi=zeroinfl(nbre~carburant+zone+ageconducteur+
           offset(log(exposition)) | bs(ageconducteur),
         data=baseFREQ,dist="poisson",link="logit")
summary(regzi)

ages=seq(18,80)
proba=predict(regzi,newdata=data.frame(ageconducteur=ages,
                                       carburant="E",
                                       exposition=1,
                                       zone="A"),
              type="zero")

data=sinistres
head(data)
base=merge(data,contrat)
reg=lm(log(cout)~agevehicule+carburant+zone,data=base)
  • Cours 9
reg1=lm(log(cout)~agevehicule+zone,data=base)
reg1b=glm(log(cout)~agevehicule+zone,data=base,
          family=gaussian(link="identity"))
summary(reg1)
s=summary(reg1)$sigma
age=seq(0,25)
nd=data.frame(agevehicule=age,zone="A")
C.pred=exp(predict(reg1,newdata=nd)+s^2/2)
plot(age,C.pred)
abline(h=mean(base$cout),col="red")

plot(age,predict(reg1,newdata=nd))
abline(h=mean(log(base$cout)),col="red")

plot(age,C.pred,ylim=c(1600,3000))
abline(h=mean(base$cout),col="red")
reg2=glm(cout~agevehicule+zone,data=base,
          family=Gamma(link="log"))
summary(reg2)
C.pred2=predict(reg2,newdata=nd,type="response")
points(age,C.pred2,col="blue")

plot(base$agevehicule,base$cout)
plot(base$agevehicule,log(base$cout))

seuil=200000
# partie 1
Index=base$cout< seuil
reg1.1=lm(log(cout)~agevehicule+zone,data=base[Index,])
reg1.2=glm(cout~agevehicule+zone,data=base[Index,],
         family=Gamma(link="log"))
s.1=summary(reg1.1)$sigma
C1.1=exp(predict(reg1.1,newdata=nd)+s.1^2/2)
C1.2=predict(reg1.2,newdata=nd,type="response")

plot(age,C1.1,ylim=c(1200,2000))
points(age,C1.2,col="blue")
abline(h=mean(base$cout[Index]),col="red")
# partie 2
C2=rep(mean(base$cout[!Index]),nrow(nd))
# partie 3
reg3=glm((cout>seuil)~agevehicule+zone,data=base,
         family=binomial)
P3.1=predict(reg3,newdata=nd,type='response')
plot(age,P3.1)
P3.2=mean(base$cout>seuil)

coutmoyen.1=C1.1*(1-P3.1)+C2*P3.1
coutmoyen.2=C1.1*(1-P3.2)+C2*P3.2
coutmoyen.3=C1.2*(1-P3.2)+C2*P3.2
plot(age,coutmoyen.2,ylim=c(1200,3000))
abline(h=mean(base$cout),col="red")
lines(age,coutmoyen.3,col="blue")

dessin=function(seuil){
# partie 1
Index=base$cout< seuil
reg1.1=lm(log(cout)~agevehicule+zone,data=base[Index,])
reg1.2=glm(cout~agevehicule+zone,data=base[Index,],
           family=Gamma(link="log"))
s.1=summary(reg1.1)$sigma
C1.1=exp(predict(reg1.1,newdata=nd)+s.1^2/2)
C1.2=predict(reg1.2,newdata=nd,type="response")
# partie 2
C2=rep(mean(base$cout[!Index]),nrow(nd))
# partie 3
reg3=glm((cout>seuil)~agevehicule+zone,data=base,
         family=binomial)
P3.1=predict(reg3,newdata=nd,type='response')
P3.2=mean(base$cout>seuil)

coutmoyen.1=C1.1*(1-P3.1)+C2*P3.1
coutmoyen.2=C1.1*(1-P3.2)+C2*P3.2
coutmoyen.3=C1.2*(1-P3.2)+C2*P3.2
plot(age,coutmoyen.2,ylim=c(1500,2300),type="l")
abline(h=mean(base$cout),col="red")
lines(age,coutmoyen.3,col="blue")}
dessin(100000)

library(nnet)
s1=1120
s2=1220
I1=base$cout<s1
I2=base$cout<s2
I3=base$cout<10000
I=I1+I2+I3
base$classes=as.factor(I)
table(base$classes)
library(splines)
regm=multinom(classes~bs(agevehicule)+zone,data=base)
summary(regm)
P=predict(regm,newdata=nd,type="probs")
barplot(t(P),col=c("yellow","red","green","pink"))
  • Cours 10
library(gdata)
  db=read.xls(
 "https://perso.univ-rennes1.fr/arthur.charpentier/SIN_1985_2000-PE.xls",
  sheet=1)
  date=db$DSUR
  D=as.Date(as.character(date),format="%Y%m%d")
  vD=seq(min(D),max(D),by=1)
  sD=table(D)
  d1=as.Date(names(sD))
  d2=vD[-which(vD%in%d1)]
  vecteur.date=c(d1,d2)
  vecteur.cpte=c(as.numeric(sD),rep(0,length(d2)))
  base=data.frame(date=vecteur.date,cpte=vecteur.cpte)
  plot(vecteur.date,vecteur.cpte,type="h",xlim=as.Date(as.character(
  c(19850101,20111231)),format="%Y%m%d"))
 regf=glm(cpte~date,family=poisson,data=base)
 mean(base$cpte)
 abline(h= mean(base$cpte),col="red")
 plot(D,db$COUTSIN,type="h")
 quantile(db$COUTSIN,1-32/2387)
 d=15
 sort(db$COUTSIN,decreasing=TRUE)[1:20]
 c=35
 db$COUTSIN=db$COUTSIN/1e6
 indm=function(x) pmin(35,pmax(x-15,0))
 mean(indm(db$COUTSIN))
 mean(indm(db$COUTSIN))*150
 library(evir)
 seuil=5
 fit=gpd(db$COUTSIN,seuil)
E <- function(yinf,ysup,xi,beta,threshold){
as.numeric(integrate(function(x) 
(x-yinf)*
dgpd(x,xi,mu=threshold,beta),
lower=yinf,upper=ysup)$value+
(1-pgpd(ysup,xi,mu=threshold,beta))*
             (ysup-yinf))
}
 E(15,50,fit$par.ests[1],fit$par.ests[2],5)
 mean(db$COUTSIN>seuil)
 E(100,200,fit$par.ests[1],fit$par.ests[2],5)*
 mean(db$COUTSIN>seuil)*150
  • Cours 11
source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
n=6
L=rep(NA,n-1)
EST.P=PAID
for(j in 1:(n-1)){
L[j]=sum(PAID[1:(n-j),j+1])/sum(PAID[1:(n-j),j])
EST.P[(n-j+1):n,j+1]=L[j]*EST.P[(n-j+1):n,j]
}
C.ULTIME=EST.P[,n]
C.AS.AT=diag(EST.P[,n:1])
R=C.ULTIME-C.AS.AT
sum(R)
PAID2=PAID
PAID2[1,1]=NA
L=rep(NA,n-1)
EST.P=PAID2
for(j in 1:(n-1)){
  L[j]=coefficients(lm(PAID2[,j+1]~0+PAID2[,j],weights=1/PAID2[,j]))
  EST.P[(n-j+1):n,j+1]=L[j]*EST.P[(n-j+1):n,j]
}
EST.P[1,1]=EST.P[1,2]/L[1]
C.ULTIME=EST.P[,n]
C.AS.AT=diag(EST.P[,n:1])
R=C.ULTIME-C.AS.AT
sum(R)
Y=log(L-1)
plot(1:5,Y)
db=data.frame(N=1:5,Y)
reg=lm(Y~N,data=db)
pY=predict(reg,newdata=data.frame(N=6:100))
prod(1+exp(pY))
R=C.ULTIME*prod(1+exp(pY))-C.AS.AT
library(ChainLadder)
MackChainLadder(PAID)
plot(MackChainLadder(PAID))
# ---------------------- REGRESSIONS
INC=PAID
INC[,2:n]=PAID[,2:n]-PAID[,1:(n-1)]
Y=as.vector(INC)
L=rep(2001:2006,n)
C=rep(1:6,each=n)
df=data.frame(Y,L,C)
reg1=lm(log(Y)~as.factor(L)+as.factor(C),data=df)
sigma=summary(reg1)$sigma
mu=predict(reg1,newdata=df)
MY=matrix(exp(mu+sigma^2/2),n,n)
INC.P=EST.P
INC.P[,2:n]=EST.P[,2:n]-EST.P[,1:(n-1)]
df$YP=exp(mu+sigma^2/2)
R=sum(df$YP[is.na(df$Y)])
reg2=glm(Y~as.factor(L)+as.factor(C),data=df,family=poisson)
mu=predict(reg2,newdata=df,type='response')
df$YP2=mu
MY=matrix(mu,n,n)
sum(df$YP2[is.na(df$Y)])
  • Cours 12
plot(cars)
reg=lm(dist~speed,data=cars)
vx=seq(4,26)
yp=predict(reg,newdata=data.frame(speed=vx))
lines(vx,yp,col="red")
vcov(reg)
sqrt(diag(vcov(reg)))
yp=predict(reg,newdata=
data.frame(speed=vx),interval="confidence")
yp
plot(cars)
lines(vx,yp[,1],lwd=2,col="red")
lines(vx,yp[,2],lwd=2,col="blue")
lines(vx,yp[,3],lwd=2,col="blue")
yp=predict(reg,newdata=
 data.frame(speed=vx),interval="prediction")
lines(vx,yp[,2],lwd=2,col="green")
lines(vx,yp[,3],lwd=2,col="green")
bruit=residuals(reg)
X=cars$speed
Xb=sample(X,size=50,replace=TRUE)
Eb=sample(bruit,size=50,replace=TRUE)
Yp=predict(reg,newdata=
             data.frame(speed=Xb))
Yb=Yp+Eb
carsb=data.frame(speed=Xb,dist=Yb)
plot(cars)
points(Xb,Yb,pch=19,col="blue")
abline(reg)
regb=lm(dist~speed,data=carsb)
abline(regb,col="blue")
plot(cars)
abline(reg)
jeu=function(){
Xb=sample(X,size=50,replace=TRUE)
Eb=sample(bruit,size=50,replace=TRUE)
Yp=predict(reg,newdata=
             data.frame(speed=Xb))
Yb=as.vector(Yp)+as.vector(Eb)
carsb=data.frame(speed=Xb,dist=Yb)
regb=lm(dist~speed,data=carsb)
abline(regb,col="blue")}
M=matrix(NA,length(vx),1000)
for(b in 1:1000){
Xb=sample(X,size=50,replace=TRUE)
Eb=sample(bruit,size=50,replace=TRUE)
Yp=predict(reg,newdata=
data.frame(speed=Xb))
Yb=as.vector(Yp)+as.vector(Eb)
carsb=data.frame(speed=Xb,dist=Yb)
regb=lm(dist~speed,data=carsb)
M[,b]=predict(regb,newdata=
data.frame(speed=vx))
}
b025=apply(M,1,function(x)
  quantile(x,.025))
b975=apply(M,1,function(x)
  quantile(x,.975))
plot(cars)
abline(reg)           
lines(vx,b025,col="red")
lines(vx,b975,col="red")
yp=predict(reg,newdata=
           data.frame(speed=vx),
           interval="confidence")
lines(vx,yp[,2],lwd=1,col="blue")
lines(vx,yp[,3],lwd=1,col="blue")
plot(ecdf(bruit))
abline(v=0,col="green")
abline(h=.54,col="green")
u=seq(-50,50,by=.1)
lines(u,pnorm(u,mean(bruit),sd(bruit)),col="red")
# ---------------- reg poisson
plot(cars)
reg=glm(dist~speed,data=cars,family=poisson)
vx=seq(4,26)
yp=predict(reg,
newdata=data.frame(speed=vx),
type="response")
lines(vx,yp,col="red")
ypic=predict(reg,
newdata=data.frame(speed=vx),
type="response",se.fit=TRUE)
ylo=ypic$fit-1.96*ypic$se.fit
yup=ypic$fit+1.96*ypic$se.fit
lines(vx,ylo,col="blue")
lines(vx,yup,col="blue")
bruitpearson=residuals(reg,type="pearson")
jeu=function(){
  Xb=sample(X,size=50,replace=TRUE)
  Eb=sample(bruitpearson,size=50,replace=TRUE)
  Yp=predict(reg,newdata=
             data.frame(speed=Xb),
             type="response")
  Yb=Yp+sqrt(Yp)*Eb
  carsb=data.frame(speed=Xb,dist=Yb)
  regb=glm(dist~speed,data=carsb,family=poisson)
  ypb=predict(regb,
             newdata=data.frame(speed=vx),
             type="response")
  lines(vx,ypb,col="blue")}
M=matrix(NA,length(vx),1000)
for(b in 1:1000){
  Xb=sample(X,size=50,replace=TRUE)
  Eb=sample(bruitpearson,size=50,replace=TRUE)
  Yp=predict(reg,newdata=
               data.frame(speed=Xb),
             type="response")
  Yb=Yp+sqrt(Yp)*Eb
  carsb=data.frame(speed=Xb,dist=Yb)
  regb=glm(dist~speed,data=carsb,family=poisson)
  ypb=predict(regb,
              newdata=data.frame(speed=vx),
              type="response")
  M[,b]=ypb
}
b025=apply(M,1,function(x)
  quantile(x,.025))
b975=apply(M,1,function(x)
  quantile(x,.975))
plot(cars)
lines(vx,ypic$fit)          
lines(vx,b025,col="red")
lines(vx,b975,col="red")
lines(vx,ylo,col="blue")
lines(vx,yup,col="blue")
reg=glm(dist~speed,data=cars,family=quasipoisson)
summary(reg)
ylo2=ypic$fit-1.96*ypic$se.fit*sqrt(5.53)
yup2=ypic$fit+1.96*ypic$se.fit*sqrt(5.53)
lines(vx,ylo2,col="blue",lty=2)
lines(vx,yup2,col="blue",lty=2)
  • Cours 13
source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")

5 thoughts on “Actuariat IARD”

  1. Bonjour,

    Je suis avec l’avis de monsier Charpentier mais R est un peu compliqué j’ai du mal à comprendre les sorties des modèles ou des fonctions !!
    Je vais faire tout les exercices que vous proposez, espérant que ça va m’aider !!
    Merci pour ce site il est très intéressant.

    1. l’actuariat IARD repose sur l’utilisation de données. Il faut un logiciel pour ça. Et je pense que R est le meilleur compromis entre simplicité du langage, richesse des fonctions possibles, facilité pour trouver de l’aide et de la documentation en ligne. En plus, étant universitaire et non pas consultant pour une boite informatique, je me refuse à enseigner un logiciel payant. Or justement, R est gratuit. Et libre….

  2. Je ne sais pas, car je n’utilise plus SAS depuis 12 ans ! je dirais un truc comme

    proc genmod data=base;
    model nbr = X1 | dist = negbin link=log
    type1 type3 offset = logexposition;
    run;

    je ne sais pas si c’est du SAS 9.3. mais voilà en gros comment faire une régression binomiale négative. Maintenant sur le fond, j’ai un soucis à répondre aux questions sur SAS. 1) je n’y connais pas grand chose. 2) SAS coute cher, donc il serait bien qu’à ce prix là, le SAS Institute fournisse un logiciel facile d’accès, ou au moins avec une documentation claire. Si j’aide sur un logiciel aussi cher que SAS, la moindre des choses serait que SAS me rétribue pour faciliter l’utilisation d’un logiciel commercial. Je ne suis pas un représentant de SAS, comme je ne suis pas un représentant de Windows, donc je n’utilise pas Excel pour les méthodes de provisionnement. R est libre, gratuit et communautaire. Le logiciel s’améliore parce qu’il y a des utilisateurs, et tout le monde y a gagne. Passez donc à R, pourquoi continuer à faire du SAS ?

    1. merci pour votre réponse.
      je suis d’accord avec vous sur fond mais malheureusement la plupart des entreprises ont SAS et donc je voulais prendre l’habitude de travailler sur SAS.Et comme je devais par la suite effectuer une simulation de Monte Carlo, chercher les intervalles de confiance….(je sais pas faire tout cela sur R..)
      Sinon sur R, c plus simple?
      Je suis preneur pour tout aide sur R.
      Peut-être que je m’y sentirai plus à l’aise.
      merci.

  3. Bonjour votre blog est très intéressant mais j’aurai voulu savoir comment ajuster des fréquences de sinistres par une binomiale négative sous SAS 9.3
    merci

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.

An Open Lab-Notebook Experiment