Category Archives: Premium principle

2013, ACT2040, Actuarial science, Automne, Computer, Course, Extreme Value, Premium principle, Reinsurance, Statistics

Pricing reinsurance contracts, another case study

05/11/2013 Arthur Charpentier Leave a comment

A reinsurance case study for tomorrow’s class. The goal will be to price some nonproportional reinsurance contract, for business interruption claims. Consider the following dataset,

> library(gdata)
>  db=read.xls(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/SIN_1985_2000-PE.xls",
+  sheet=1)
Content type 'application/vnd.ms-excel' length 183808 bytes (179 Kb)
open URL
==================================================
downloaded 179 Kb

As for any (standard) insurance contract, there are two parts in the pricing

the expected number of claims
the average cost of individual claims

Here, we do not have covariates (but it might be possible to use some, like the kind of industry, the location, etc).

Let us start with the expected number of claims, per year. Here is the daily frequency,

The data are rather old… but somehow, it is a good thing since after ten years, we can expect that most of the claims have been settled (we’ll discuss claims dynamic starting next week). To plot the graph above, we use

> 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"))

Then, we can get a prediction of the daily number of business interruption claims, e.g. for any day in 2010 (assume that we had to price a reinsurance contract a few years ago), using a (standard) Poisson regression

> regdate=glm(cpte~date,data=base,family=poisson(link="log"))
> nd2010=data.frame(date=seq(as.Date(as.character(20100101),format="%Y%m%d"),
+ as.Date(as.character(20101231),format="%Y%m%d"),by=1))
> pred2010 =predict(regdate,newdata=nd2010,type="response")
> sum(pred2010)
[1] 159.4757

Observe that using old data has drawbacks, since we got much more uncertainty if we use a regression on time (to include some possible trend)

Say we have something like 160 claims over a given year, on average.

> plot(D,db$COUTSIN,type="h")

Let us now focus on the cost of those claims. We have 2,400 claims in our dataset, to fit a model (or at least estimate how much a reinsurance contract might cost us). Assume that we would like to purchase a reinsurance contract for our very large claims. Like the two largest per year. Over 16 years, the decutible should be close to the cost of the 32nd largest claim, which was close to 15 million.

> quantile(db$COUTSIN,1-32/2400)/1e6
98.66667% 
 15.34579 
> abline(h=quantile(db$COUTSIN,1-32/2400),col="blue")

So consider some reinsurance contract with a deductible of 15 million. Unfortunately, we cannot find unlimited covers. So let us assume that a reinsurance company agrees for such a deductible, but with a limited cover of 35 million. The average cost (for the reinsurance company) is $https://latex.codecogs.com/gif.latex?\mathbb{E}(g(X))$ where

$https://latex.codecogs.com/gif.latex?g(x)=\min\{35,\max\{x-15,0\}\}$

A first idea is to look at the first cost, i.e. the empirical average of that indemnity, on our portfolio. The indemnity function is

> indemn=function(x) pmin((x-15)*(x>15),50-15)

we can check on a few losses that it is actually what we wish to compute

> indemn(5)
[1] 0
> indemn(20)
[1] 5
> indemn(50)
[1] 35

Now, if the compute the average repayment by the reinsurance company, over 16 years, we get

> mean(indemn(db$COUTSIN/1e6))
[1] 0.1624292

So, per claim, the reinsurance company will pay, on average 162,430. With 160 claims per year, the pure premium should be close to 26 million

> mean(indemn(db$COUTSIN/1e6))*160
[1] 25.98867

(again, for a 35 million cover, for some claims that should occur, on average, twice a year). As we will see, a standard model in reinsurance is the Pareto distribution (or to be more specific, a Generalized Pareto one),

$F_{(\xi,\mu,\sigma)}(x) = \begin{cases}1 - \left(\ds{1+ \frac{\xi(x-\mu)}{\sigma}\right)^{-1/\xi}} & \text{for }\xi \neq 0, \\1 - \exp \left(-\frac{x-\mu}{\sigma}\right) & \text{for }\xi = 0.\end{cases}$

There are three parameters here

the threshold $https://latex.codecogs.com/gif.latex?\mu$ (that we will consider as fixed, but we will see its impact on reinsurance pricing)
the scale parameter $https://latex.codecogs.com/gif.latex?\sigma$ (called $https://latex.codecogs.com/gif.latex?\beta$ in R)
the tail index $https://latex.codecogs.com/gif.latex?\xi$

The strategy is to consider a threshold below our deductible, e.g. 12 million. Then, given that the loss exceed 12 million, we can fit a Generalized Pareto distribution,

> gpd.PL <- gpd(db$COUTSIN,12e6)$par.ests
> gpd.PL
          xi         beta 
7.004147e-01 4.400115e+06

and compute

$\int_{d}^{d+c}(x-d)\ dF_{(\xi,\mu,\sigma)}(x)$

>  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))
+  }

Here, given that a claim exceeds 12 million, the average repayment is close to 6 million

> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6)
[1] 6058125

Now, we have to take into account the probability to reach 12 million, which is here

> mean(db$COUTSIN>12e6)
[1] 0.02639296

So, if we summarize, we have on average 160 claims per year,

> p
[1] 159.4757

Only 2.6% will exceed 12 million

> mean(db$COUTSIN>12e6)
[1] 0.02639296

So, the yearly frequency of claism larger than 12 million is 4.2 claims

> p*mean(db$COUTSIN>12e6)
[1] 4.209036

And for a claim that exceed 12 million, the average repayment is

> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6)
[1] 6058125

So, the pure premium should be close to

> p*mean(db$COUTSIN>12e6)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6)
[1] 25498867

which (hopefully) is close to the empirical value we got. Actually, it is also possible to look at the impact of the threshold parameter, since it is clearly and intermediate value that could be changed. I mean, why 12 and not 10? Consider

> esp=function(threshold=12e6,p=sum(pred2010)){
+  (gpd.PL <- gpd(db$COUTSIN,threshold)$par.ests)
+  return(p*mean(db$COUTSIN>threshold)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],threshold))
+  }

We can plot the pure premium as a function of that threshold,

> seuils=seq(1e6,15e6,by=1e6)
> plot(seuils,Vectorize(esp)(seuils),type="b",col="red")

which is between 24 and 26 for large thresholds. Again, that is only the first step, and we can price a higher reinsurance layer, like a reinsurance contract with a deductible of 50 million (we have our previous reinsurance contract for claims below that threshold), and a cover of 50 million, for instance. For those high layers, it become interesting to have a parametric model, which should be more robust than the empirical average.

2013, ACT2040, Actuarial science, Automne, Course, GLM, Premium principle, Ratemaking

Régression de Poisson

23/09/2013 Arthur Charpentier Leave a comment

Mercredi, on finira les arbres de classification, et on commencera la modélisation de la fréquence de sinistre. Les transparents sont en ligne.

Comme annoncé lors du premier cours, je suggère de commencer la lecture du Practicionner’s Guide to Generalized Linear Models. Le document correspond au minimum attendu dans ce cours.

2013, ACT2040, Actuarial science, Exams, Hiver, IBNR, Premium principle, Ratemaking

Provisionnement et tarification, examen final

24/04/2013 Arthur Charpentier 4 Comments

L’examen final du cours ACT2040 avait lieu ce matin. L’énoncé est en ligne ainsi que des éléments de correction. En cas d’erreurs, merci de me le faire savoir rapidement, avant que je ne saisisse les notes.

2013, ACT2040, Computer, GLM, GLM, Hiver, Premium principle, Ratemaking

Multiple (smoothed) regression and portfolio exposure

09/03/2013 Arthur Charpentier Leave a comment

Wednesday, in class, we’ve seen how to visualize a multiple regression model (with two continuous explanatory variables). Here, the goal is to predict the average cost of an insurance claim, using some covariates, e.g. the age of the driver, and the age of the car (recall that losses here are liability losses). The prediction obtained from a (standard) generalized linear model, with a log-link

> reg1=glm(cout~ageconducteur+agevehicule,data=base,family=Gamma(link="log"))

The code to visualize the predicted average cost is the following: first, we have to compute predictions for specific values,

> pred=function(x,y){
+ predict(reg,newdata=data.frame(ageconducteur=x,
+ agevehicule=y),type="response")

Then, we use this function to compute values on a grid,

> X=seq(20,80,by=5)
> Y=0:20
> Z=outer(X,Y,p)
> image(X,Y,Z,col=rev(heat.colors(101)))
> contour(X,Y,Z,add=TRUE,
+ levels=c(1400,1800,2000,2200,2400,2600,2800,3000,3200,4000,5000))

If we use factors, and not continuous variates (cut versions of those two variates),

> reg2=glm(cout~cut(ageconducteur,breaks=c(0,22,35,55,80,100))*
+               cut(agevehicule,breaks=c(-1,1,3,5,10,100)),
+ data=base,family=Gamma(link="log"))

(note that we consider the Cartesian product, so values are computed for each product of factors, age of the driver and age of the car) we obtain

Obviously, we’re missing something here: the most expensive class with one model is the cheapeast for the other one! Of course, it might come from our classes (that were chosen a bit randomly), but it might be interesting to use nonlinear functions of the ages. So, let us use splines to smooth those two variables,

> reg3=glm(cout~bs(ageconducteur)+bs(agevehicule),data=base,
+ family=Gamma(link="log"))

With additive smoothed functions, we obtained a symmetric graph (due to the additive property)

while with a bivariate spline

> library(mgcv)
+ reg4=gam(cout~s(ageconducteur,agevehicule),data=base,
+ family=Gamma(link="log"))

(for some odd reasons, I could not use – easily – a bivariate spline in the Generalized Linear Model, but it did work considering a Generalized Additive Model – which is, by no means additive now). We can identify here some regions where the average cost can be extremely expensive… But, as mentioned wednesday, one should keep in mind that some parts of the square above are not reached. More precisely, the distribution of the portfolio, as a function of those two covariates is the following

Thus, the proportion of young drivers driving a brand new car, and the proportion of old drivers driving a very old car is rather small… If the goal is to find niches, one should look at the prediction more carefully, but if the goal is to make that everyone gets an insurance cover, maybe we should allow that some drivers are under-priced (especially when are rare in the portfolio). And one should keep in mind that average costs are extremely sensitive to large losses, as discussed previously http://freakonometrics.hypotheses.org/3490 (and in class)

In the univariate case, I have migrated an old post, we I tried to reproduce (in R and in French) some standard graphs in the insurance industry: it is always interesting to visualize not only the prediction obtained from our models, but also the size of each class in the portfolio,

The post is online here http://freakonometrics.hypotheses.org/1224

2013, ACT2040, GLM, Hiver, Premium principle, Ratemaking

Visualizing overdispersion (with trees)

13/02/2013 Arthur Charpentier Leave a comment

This week, we started to discuss overdispersion when modeling claims frequency. In my previous post, I discussed computations of empirical variances with different exposure. But I did use only one factor to compute classes. Of course, it is possible to use much more factors. For instance, using cartesian products of factors,

> X=as.factor(paste(sinistres$carburant,sinistres$zone,
+ cut(sinistres$ageconducteur,breaks=c(17,24,40,65,101))))
> E=sinistres$exposition
> Y=sinistres$nbre
> vm=vv=ve=rep(NA,length(levels(X)))
>   for(i in 1:length(levels(X))){
+  	   ve[i]=Ei=E[X==levels(X)[i]]
+  	   Yi=Y[X==levels(X)[i]]
+   vm[i]=meani=weighted.mean(Yi/Ei,Ei)    # moyenne 
+   vv[i]=variancei=sum((Yi-meani*Ei)^2)/sum(Ei)    # variance
+  cat("Class ",levels(X)[i],"average =",meani," variance =",variancei,"\n")
+ }
Class D A (17,24]  average = 0.06274415  variance = 0.06174966 
Class D A (24,40]  average = 0.07271905  variance = 0.07675049 
Class D A (40,65]  average = 0.05432262  variance = 0.06556844 
Class D A (65,101] average = 0.03026999  variance = 0.02960885 
Class D B (17,24]  average = 0.2383109   variance = 0.2442396 
Class D B (24,40]  average = 0.06662015  variance = 0.07121064 
Class D B (40,65]  average = 0.05551854  variance = 0.05543831 
Class D B (65,101] average = 0.0556386   variance = 0.0540786 
Class D C (17,24]  average = 0.1524552   variance = 0.1592623 
Class D C (24,40]  average = 0.0795852   variance = 0.09091435 
Class D C (40,65]  average = 0.07554481  variance = 0.08263404 
Class D C (65,101] average = 0.06936605  variance = 0.06684982 
Class D D (17,24]  average = 0.1584052   variance = 0.1552583 
Class D D (24,40]  average = 0.1079038   variance = 0.121747 
Class D D (40,65]  average = 0.06989518  variance = 0.07780811 
Class D D (65,101] average = 0.0470501   variance = 0.04575461 
Class D E (17,24]  average = 0.2007164   variance = 0.2647663 
Class D E (24,40]  average = 0.1121569   variance = 0.1172205 
Class D E (40,65]  average = 0.106563    variance = 0.1068348 
Class D E (65,101] average = 0.1572701   variance = 0.2126338 
Class D F (17,24]  average = 0.2314815   variance = 0.1616788 
Class D F (24,40]  average = 0.1690485   variance = 0.1443094 
Class D F (40,65]  average = 0.08496827  variance = 0.07914423 
Class D F (65,101] average = 0.1547769   variance = 0.1442915 
Class E A (17,24]  average = 0.1275345   variance = 0.1171678 
Class E A (24,40]  average = 0.04523504  variance = 0.04741449 
Class E A (40,65]  average = 0.05402834  variance = 0.05427582 
Class E A (65,101] average = 0.04176129  variance = 0.04539265 
Class E B (17,24]  average = 0.1114712   variance = 0.1059153 
Class E B (24,40]  average = 0.04211314  variance = 0.04068724 
Class E B (40,65]  average = 0.04987117  variance = 0.05096601 
Class E B (65,101] average = 0.03123003  variance = 0.03041192 
Class E C (17,24]  average = 0.1256302   variance = 0.1310862 
Class E C (24,40]  average = 0.05118006  variance = 0.05122782 
Class E C (40,65]  average = 0.05394576  variance = 0.05594004 
Class E C (65,101] average = 0.04570239  variance = 0.04422991 
Class E D (17,24]  average = 0.1777142   variance = 0.1917696 
Class E D (24,40]  average = 0.06293331  variance = 0.06738658 
Class E D (40,65]  average = 0.08532688  variance = 0.2378571 
Class E D (65,101] average = 0.05442916  variance = 0.05724951 
Class E E (17,24]  average = 0.1826558   variance = 0.2085505 
Class E E (24,40]  average = 0.07804062  variance = 0.09637156 
Class E E (40,65]  average = 0.08191469  variance = 0.08791804 
Class E E (65,101] average = 0.1017367   variance = 0.1141004 
Class E F (17,24]  average = 0           variance = 0 
Class E F (24,40]  average = 0.07731177  variance = 0.07415932 
Class E F (40,65]  average = 0.1081142   variance = 0.1074324 
Class E F (65,101] average = 0.09071118  variance = 0.1170159

Again, one can plot the variance against the average,

> plot(vm,vv,cex=sqrt(ve),col="grey",pch=19,
+ xlab="Empirical average",ylab="Empirical variance")
> points(vm,vv,cex=sqrt(ve))
> abline(a=0,b=1,lty=2)

An alternative is to use a tree. The tree can be obtained from another variable (the insured had, or had not, a claim, during the period considered) but it should be rather close to the one we would like to model (the number of claims over the period considered). Here, I did use the whole database (with more that 600,000 lines)

> library(tree)
> T=tree((nombre>0)~as.factor(zone)+as.factor(puissance)+
+ as.factor(marque)+as.factor(carburant)+as.factor(region)+
+ agevehicule+ageconducteur,data=baseFREQ,
+ split =  "gini",minsize =25000)

The tree is the following

> plot(T)
> text(T)

Now, each knot defines a class, and it is possible to use it to define a class. Which is supposed to be homogeneous.

> X=as.factor(T$where)
> E=sinistres$exposition
> Y=sinistres$nbre
> vm=vv=ve=rep(NA,length(levels(X)))
>   for(i in 1:length(levels(X))){
+  	   ve[i]=Ei=E[X==levels(X)[i]]
+  	   Yi=Y[X==levels(X)[i]]
+   vm[i]=meani=weighted.mean(Yi/Ei,Ei)    # moyenne 
+   vv[i]=variancei=sum((Yi-meani*Ei)^2)/sum(Ei)    # variance
+  cat("Class ",levels(X)[i],"average =",meani," variance =",variancei,"\n")
+  }
Class  6 average =   0.04010406  variance = 0.04424163 
Class  8 average =   0.05191127  variance = 0.05948133 
Class  9 average =   0.07442635  variance = 0.08694552 
Class  10 average =  0.4143646   variance = 0.4494002 
Class  11 average =  0.1917445   variance = 0.1744355 
Class  15 average =  0.04754595  variance = 0.05389675 
Class  20 average =  0.08129577  variance = 0.0906322 
Class  22 average =  0.05813419  variance = 0.07089811 
Class  23 average =  0.06123807  variance = 0.07010473 
Class  24 average =  0.06707301  variance = 0.07270995 
Class  25 average =  0.3164557   variance = 0.2026906 
Class  26 average =  0.08705041  variance = 0.108456 
Class  27 average =  0.06705214  variance = 0.07174673 
Class  30 average =  0.05292652  variance = 0.06127301 
Class  31 average =  0.07195285  variance = 0.08620593 
Class  32 average =  0.08133722  variance = 0.08960552 
Class  34 average =  0.1831559   variance = 0.2010849 
Class  39 average =  0.06173885  variance = 0.06573939 
Class  41 average =  0.07089419  variance = 0.07102932 
Class  44 average =  0.09426152  variance = 0.1032255 
Class  47 average =  0.03641669  variance = 0.03869702 
Class  49 average =  0.0506601   variance = 0.05089276 
Class  50 average =  0.06373107  variance = 0.06536792 
Class  51 average =  0.06762947  variance = 0.06926191 
Class  56 average =  0.06771764  variance = 0.07122379 
Class  57 average =  0.04949142  variance = 0.05086885 
Class  58 average =  0.2459016   variance = 0.2451116 
Class  59 average =  0.05996851  variance = 0.0615773 
Class  61 average =  0.07458053  variance = 0.0818608 
Class  63 average =  0.06203737  variance = 0.06249892 
Class  64 average =  0.07321618  variance = 0.07603106 
Class  66 average =  0.07332127  variance = 0.07262425 
Class  68 average =  0.07478147  variance = 0.07884597 
Class  70 average =  0.06566728  variance = 0.06749411 
Class  71 average =  0.09159605  variance = 0.09434413 
Class  75 average =  0.03228927  variance = 0.03403198 
Class  76 average =  0.04630848  variance = 0.04861813 
Class  78 average =  0.05342351  variance = 0.05626653 
Class  79 average =  0.05778622  variance = 0.05987139 
Class  80 average =  0.0374993   variance = 0.0385351 
Class  83 average =  0.06721729  variance = 0.07295168 
Class  86 average =  0.09888492  variance = 0.1131409 
Class  87 average =  0.1019186   variance = 0.2051122 
Class  88 average =  0.05281703  variance = 0.0635244 
Class  91 average =  0.08332136  variance = 0.09067632 
Class  96 average =  0.07682093  variance = 0.08144446 
Class  97 average =  0.0792268   variance = 0.08092019 
Class  99 average =  0.1019089   variance = 0.1072126 
Class  100 average = 0.1018262   variance = 0.1081117 
Class  101 average = 0.1106647   variance = 0.1151819 
Class  103 average = 0.08147644  variance = 0.08411685 
Class  104 average = 0.06456508  variance = 0.06801061 
Class  107 average = 0.1197225   variance = 0.1250056 
Class  108 average = 0.0924619   variance = 0.09845582 
Class  109 average = 0.1198932   variance = 0.1209162

Here, when ploting the empirical variance (per knot) against the empirial average of claims, we get

Here, we can identify classes where remaining heterogeneity.

ACT2040, Computer, GLM, Graphics, Premium principle, Ratemaking

Overdispersion with different exposures

01/02/2013 Arthur Charpentier 1 Comment

In actuarial science, and insurance ratemaking, taking into account the exposure can be a nightmare (in datasets, some clients have been here for a few years – we call that exposure – while others have been here for a few months, or weeks). Somehow, simple results because more complicated to compute just because we have to take into account the fact that exposure is an heterogeneous variable.

The exposure in insurance ratemaking can be seen as a problem of censored data (in my dataset, the exposure is always smaller than 1 since observations are contracts, not policyholders),

the number of claims $https://latex.codecogs.com/gif.latex?N_i$ on the period $https://latex.codecogs.com/gif.latex?[0,1]$ is unobserved
the number of claims $https://latex.codecogs.com/gif.latex?Y_i$ on $https://latex.codecogs.com/gif.latex?[0,E_i]$ is observed (as well as $https://latex.codecogs.com/gif.latex?E_i$ )

And as always, the variable of interest is the unobserved one, because we have to price insurance contract with a cover period of one (full) year. So we have to model the yearly frequency of insurance claims.

In our dataset, we have $https://latex.codecogs.com/gif.latex?(Y_i,E_i)$ ‘s – or more generally also some additional covariates $https://latex.codecogs.com/gif.latex?(Y_i,E_i,\boldsymbol{X}_i)$ ‘s. For ratemaking, we need to estimate $https://latex.codecogs.com/gif.latex?\mathbb{E}(N\vert\boldsymbol{X}=\boldsymbol{x})$ and perhaps also $https://latex.codecogs.com/gif.latex?\text{Var}(N|\boldsymbol{X}=\boldsymbol{x})$ (for instance to test if the Poisson assumption is valid, or not). To estimate the expected value, a natural estimate for $https://latex.codecogs.com/gif.latex?\mathbb{E}(N)$ (forget about covariates as a start) is
$https://latex.codecogs.com/gif.latex?m_N=\frac{\sum_{i=1}^n%20Y_i}{\sum_{i=1}^n%20E_i}$
which is also the weight average of annualized individual counts
$https://latex.codecogs.com/gif.latex?m_N=\sum_{i=1}^n%20\frac{%20E_i}{\sum_{i=1}^n%20E_i}%20\cdot%20\frac{Y_i}{E_i}$
We consider the ratio of the total number of claims to the total exposure-to-
risk. This estimate appears for instance if we consider a Poisson process, so that $https://latex.codecogs.com/gif.latex?N\sim\mathcal{P}(\lambda)$ while $https://latex.codecogs.com/gif.latex?Y\sim\mathcal{P}(\lambda%20\cdot%20E)$ . Then, the likelihood is

$https://latex.codecogs.com/gif.latex?\mathcal{L}(\lambda,\boldsymbol{Y},\boldsymbol{E})=\prod_{i=1}^n%20\frac{e^{-\lambda%20E_i}%20[\lambda%20E_i]^{Y_i}}{Y_i!}$

i.e.

$https://latex.codecogs.com/gif.latex?\log%20\mathcal{L}(\lambda,\boldsymbol{Y},\boldsymbol{E})%20=%20-\lambda%20\sum_{i=1}^n%20E_i%20+\sum_{i=1}^n%20Y_i%20\log[\lambda%20E_i]%20-%20\log\left(\prod_{i=1}^n%20Y_i!\right)$

The first order condition is here

$https://latex.codecogs.com/gif.latex?\frac{\partial}{\partial%20\lambda}\log%20\mathcal{L}(\lambda,\boldsymbol{Y},\boldsymbol{E})%20=%20%20-%20\sum_{i=1}^n%20E_i%20+\frac{1}{\lambda}\sum_{i=1}^n%20Y_i%20=0$

which is satisfied if

$https://latex.codecogs.com/gif.latex?\widehat{\lambda}=\frac{\sum_{i=1}^n%20Y_i}{\sum_{i=1}^n%20E_i}$

So, we do have an estimator for the expected value, and a natural estimator for $https://latex.codecogs.com/gif.latex?\mathbb{E}(N\vert\boldsymbol{X}=\boldsymbol{x})$ is then (if we consider categorical covariates)
$https://latex.codecogs.com/gif.latex?m_{N|\boldsymbol{x}}%20=\frac{\sum_{i,\boldsymbol{X}_i=\boldsymbol{x}}%20Y_i}{\sum_%20{i,\boldsymbol{X}_i=\boldsymbol{x}}%20E_i}$

Now, we need an estimate for the variance, or more precisely the conditional variable. Assume (as a starting point) that all have the same exposure $https://latex.codecogs.com/gif.latex?E$ . For instance, if $https://latex.codecogs.com/gif.latex?E$ is one half, insured were observed only the first six months. Then $https://latex.codecogs.com/gif.latex?N=Y+Y%27$ with $https://latex.codecogs.com/gif.latex?Y\overset{\mathcal%20L}{=}Y%27$ ( $https://latex.codecogs.com/gif.latex?Y$ is the number of claims on the first six months, while $https://latex.codecogs.com/gif.latex?Y%27$ are the number of claims on the last six months), i.e. $https://latex.codecogs.com/gif.latex?\text{Var}(N)=\text{Var}(Y)+%20\text{Var}(Y%27)$ if we assume independent increments. I.e.
$https://latex.codecogs.com/gif.latex?\text{Var}(N)=2\text{Var}(Y)$ , or conversely $https://latex.codecogs.com/gif.latex?E%20\cdot\text{Var}(N)=\text{Var}(Y)$ . More generally, it is reasonable to assume that

$https://latex.codecogs.com/gif.latex?\text{Var}(Y)=E\cdot%20\text{Var}(N)$
for all values of $https://latex.codecogs.com/gif.latex?E$ . And then
$https://latex.codecogs.com/gif.latex?\text{Var}\left(\frac{Y}{E}\right)=\frac{1}{E}\cdot%20\text{Var}(N)$
Thus, it seems legitimate to assume that the empirical variance of $https://latex.codecogs.com/gif.latex?N$ can be written
$https://latex.codecogs.com/gif.latex?S_N^2=E\cdot%20S_{Y/E}^2$
Since the average of $https://latex.codecogs.com/gif.latex?Y_i/E$ is $https://latex.codecogs.com/gif.latex?\overline{N}=m_N$ , then
$https://latex.codecogs.com/gif.latex?S_N^2=E\cdot%20\frac{1}{n}\sum_{i=1}^n%20\left[\frac{Y_i}{E}-\overline{N}\right]^2}%20=%20\frac{1}{n}\sum_{i=1}^n%20E\left[\frac{Y_i}{E}-\overline{N}\right]^2}$
or equivalently
$https://latex.codecogs.com/gif.latex?S_N^2=\frac{1}{n}\sum_{i=1}^n%20\frac{E}{E^2}\left[Y_i-\overline{N}\cdot%20E\right]^2}%20=\frac{1}{n}\sum_{i=1}^n%20\frac{1}{E}[Y_i-\overline{N}\cdot%20E]^2$ i.e.
$https://latex.codecogs.com/gif.latex?S_N^2=\frac{\sum_{i=1}^n%20[Y_i-\overline{N}\cdot%20E]^2%20}{nE}$
Thus, with different $https://latex.codecogs.com/gif.latex?E_i$ ‘s, it would be legitimate (I guess) to consider
$https://latex.codecogs.com/gif.latex?S_N^2=\frac{\sum_{i=1}^n%20[Y_i-\overline{N}\cdot%20E_i]^2%20}{\sum_{i=1}^n%20E_i}$
Thus, an estimator for $https://latex.codecogs.com/gif.latex?\text{Var}(N|\boldsymbol{X}=\boldsymbol{x})$ is
$https://latex.codecogs.com/gif.latex?S_{N|\boldsymbol{x}}^2=\frac{\sum_{i,\boldsymbol{X}_i=\boldsymbol{x}}%20[Y_i-\overline{N}\cdot%20E_i]^2}{\sum_{i,\boldsymbol{X}_i=\boldsymbol{x}%20}%20E_i}$

This can be used to test is the Poisson assumption is valid to model frequency. Consider the following dataset,

>  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)

Here, we do have our two variables of interest, the exposure, per contract,

>  E <- baseFREQ$exposition

and the (observed) number of claims (during that time frame)

>  Y <- baseFREQ$nbre

It is possible to compute without covariates, the average (yearly) number of claims, per contract, and the associated variance

> (mean=weighted.mean(Y/E,E))
[1] 0.07279295
> (variance=sum((Y-mean*E)^2)/sum(E)) 
[1] 0.08778567

It looks like the variance is (slightly) larger than the average (we’ll see in a few weeks how to test it, more formally). It is possible to add covariates, for instance the density of population, in the area where the policyholder lives,

>  X=as.factor(baseFREQ$densite)
>  for(i in 1:length(levels(X))){
+ 	   Ei=E[X==levels(X)[i]]
+ 	   Yi=Y[X==levels(X)[i]]
+  (meani=weighted.mean(Yi/Ei,Ei))    # moyenne 
+  (variancei=sum((Yi-meani*Ei)^2)/sum(Ei))    # variance
+ cat("Density, zone",levels(X)[i],"average =",meani," variance =",variancei,"\n")
+ }
Density, zone 11 average = 0.07962411  variance = 0.08711477 
Density, zone 21 average = 0.05294927  variance = 0.07378567 
Density, zone 22 average = 0.09330982  variance = 0.09582698 
Density, zone 23 average = 0.06918033  variance = 0.07641805 
Density, zone 24 average = 0.06004009  variance = 0.06293811 
Density, zone 25 average = 0.06577788  variance = 0.06726093 
Density, zone 26 average = 0.0688496   variance = 0.07126078 
Density, zone 31 average = 0.07725273  variance = 0.09067 
Density, zone 41 average = 0.03649222  variance = 0.03914317 
Density, zone 42 average = 0.08333333  variance = 0.1004027 
Density, zone 43 average = 0.07304602  variance = 0.07209618 
Density, zone 52 average = 0.06893741  variance = 0.07178091 
Density, zone 53 average = 0.07725661  variance = 0.07811935 
Density, zone 54 average = 0.07816105  variance = 0.08947993 
Density, zone 72 average = 0.08579731  variance = 0.09693305 
Density, zone 73 average = 0.04943033  variance = 0.04835521 
Density, zone 74 average = 0.1188611   variance = 0.1221675 
Density, zone 82 average = 0.09345635  variance = 0.09917425 
Density, zone 83 average = 0.04299708  variance = 0.05259835 
Density, zone 91 average = 0.07468126  variance = 0.3045718 
Density, zone 93 average = 0.08197912  variance = 0.09350102 
Density, zone 94 average = 0.03140971  variance = 0.04672329

Perhaps graphs would be a nice tool to play with, to visualize that information

> plot(meani,variancei,cex=sqrt(Ei),col="grey",pch=19,
+ xlab="Empirical average",ylab="Empirical variance")
> points(meani,variancei,cex=sqrt(Ei))

The size of the circles is related to the size of the group (the area is proportional to the total exposure within the group). The first diagonal corresponds to the Poisson model, i.e. the variance should be equal to the mean. It is also possible to consider other covariates, like the gas type

or the car brand,

It is also possible to consider the age of the driver as a categorical variate

Actually, the age is interesting: we can observe on that dataset a feature that Jean-Philippe Boucher observed also on his own datasets. Let us look more carefully where are the different ages,

On the right, we can observe young (unexperienced) drivers. That was expected. But some classes are below the first diagonal: the expected frequency is large, but not the variance. I.e. we know for sure that young drivers have more car accidents. It is not an heterogeneous class, on the contrary: young drivers can be seen as a relatively homogeneous class, with a high frequency of car accidents.

With the original dataset (here, I use only a subset with 50,000 clients), we do obtain the following graph:

If we do not observe underdispersion for young drivers, observe that those are incredibly homogeneous classes. With a clear impact of experience, since circles are moving downward from age 18 to 25.

Another disturbing story (this was – one more time – suggestion from Jean-Philippe) that it might be possible to consider the exposure as a standard variable, and see if the coefficient is actually equal to 1. Without any covariate,

>  reg=glm(Y~log(E),family=poisson("log"))
>  summary(reg)

Call:
glm(formula = Y ~ log(E), family = poisson("log"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.3988  -0.3388  -0.2786  -0.1981  12.9036  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.83045    0.02822 -100.31   <2e-16 ***
log(E)       0.53950    0.02905   18.57   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 12931  on 49999  degrees of freedom
Residual deviance: 12475  on 49998  degrees of freedom
AIC: 16150

Number of Fisher Scoring iterations: 6

i.e. the parameter is clearly strictly smaller than 1. And it is neither related to significance,

> library(car)
> linearHypothesis(reg,"log(E)",1)
Linear hypothesis test

Hypothesis:
log(E) = 1

Model 1: restricted model
Model 2: Y ~ log(E)

  Res.Df Df  Chisq Pr(>Chisq)    
1  49999                         
2  49998  1 251.19  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

nor to the fact that I did not take into account covariates,

> reg=glm(nbre~log(exposition)+carburant+as.factor(ageconducteur)+as.factor(densite),family=poisson("log"),data=baseFREQ)
>  summary(reg)

Call:
glm(formula = nbre ~ log(exposition) + carburant + as.factor(ageconducteur) + 
    as.factor(densite), family = poisson("log"), data = baseFREQ)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7114  -0.3200  -0.2637  -0.1896  12.7104  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                  -14.07321  181.04892  -0.078 0.938042    
log(exposition)                0.56781    0.03029  18.744  < 2e-16 ***
carburantE                    -0.17979    0.04630  -3.883 0.000103 ***
as.factor(ageconducteur)19    12.18354  181.04915   0.067 0.946348    
as.factor(ageconducteur)20    12.48752  181.04902   0.069 0.945011

(etc). So it might be a too strong assumption to assume that the exposure is an exogenous variate here. But that’s another story !

2011, ACT2040, Computer, Graphics, Premium principle, Ratemaking, Statistics

Visualisation en tarification, avec R

21/09/2011 Arthur Charpentier Leave a comment

Chose promise, chose due: mardi prochain, le 27 septembre, le cours d’Actuariat IARD se déroulera en salle informatique (PK-S1525) à l’horaire habituel. Sinon quelques notes de cours sur la tarification a priori sont en ligne ici.

Quelques références en R

Avant de commencer à programmer en R, quelques références qui peuvent être utiles, en français pour commencer,

“R pour les débutants” d’Emmanuel Paradis, (PDF)
“Introduction à la programmation en S” par Vincent Goulet, (PDF)
“Statistique de l’Assurance” par Arthur Charpentier (PDF) qui sont mes notes de cours de l’an passé, mais qui insistent sur l’utilisation de R, pas sur la programmation en R, qui est supposée un peu connue.

Sinon en anglais, les références sont un peu plus nombreuses,

“R for Beginners” d’Emmanuel Paradis (PDF),
“An Introduction to R” par Longhow Lam (PDF)
“The R language — a short companion” par Marc Vandemeulebroecke (PDF),
“The R Guide” par Jason Owen (PDF),
“Econometrics in R” par Grant Farnsworth (PDF) pour aller plus loin sur les régressions,
“Practical Regression and Anova using R” by Julian Faraway (PDF) sur le meme sujet
“Statistics with R and S-Plus” d’Hugo Quené (PDF)
“Statistical Computing and Graphics Course Notes” par Frank Harrell, (PDF).
“Using R for Data Analysis and Graphics – Introduction, Examples and Commentary” par John Maindonald (PDF).

Le code pour importer les données est le suivant,

> sinistre <- read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
+ header=TRUE,sep=";")
> sinistres=sinistre[sinistre$garantie=="1RC",]
> 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)
> basenb = merge(contrat,nombre)
> head(basenb)
> basesin=merge(sinistres,contrat)
> basesin=basesin[basesin$cout>0,]

Faire des graphiques en R

Dans A practicioner’s guide to Generalized Linear Models (en ligne ici), on peut voir des graphiques comme celui ci-dessous, avec la fréquence de sinistre en fonction de l’age (ou plutôt de classes d’ages). Sur le graphique ci-dessous, on a également une distinction entre les hommes (en bleu) et les femmes (en rouge),

sur ce graphique, les fréquences sont exprimées en logarithme du multiplicateur (i.e. en variation par rapport à la moyenne du portefeuille).
On peut utiliser le code suivant pour générer (automatiquement) des graphiques similaires,

> graphique=function(nom="ageconducteur",
+ niveau=c(17,21,24,29,34,44,64,84,110),
+ continu=TRUE,type=1){
+ if(continu==TRUE){X=cut(basenb[,nom],niveau)}
+ if(continu==FALSE){X=as.factor(basenb[,nom])}
+ E=basenb$exposition
+ Y=basenb$nbre
+ FREQ=levels(X)
+ moyenne=variance=n=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyenne[k] =weighted.mean(Y[X==FREQ[k]]/E[X==FREQ[k]],
+ E[X==FREQ[k]])
+ variance[k]=weighted.mean((Y[X==FREQ[k]]/E[X==FREQ[k]]-
+ moyenne[k])^2,E[X==FREQ[k]])
+ n[k]       =sum(E[X==FREQ[k]])
+}
+ w=barplot(n,names.arg=FREQ,col="light green",axes=FALSE,
+ xlim=c(0,1.2*length(FREQ)+.5))
+ mid=w[,1]
+ axis(2)
+ par(new=TRUE)
+ IC1=moyenne+1.96/sqrt(n)*sqrt(variance)
+ IC2=moyenne-1.96/sqrt(n)*sqrt(variance)
+ moyenneglobale=sum(Y)/sum(E)
+ 
+ if(type==1){
+ plot(mid,moyenne,ylim=range(c(IC1,IC2)),type="b",
+ col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="red")
+ segments(mid-.1,IC1,mid+.1,IC1,col="red")
+ segments(mid-.1,IC2,mid+.1,IC2,col="red")
+ points(mid,moyenne,pch=19,col="red")
+ axis(4)
+ abline(h=moyenneglobale,lty=2,col="red")}
+
+ if(type==2){
+ plot(mid,log(moyenne/moyenneglobale),ylim=
+ range(c(log(IC1/moyenneglobale),log(IC2/moyenneglobale))),
+ type="b",col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,log(IC1/moyenneglobale),mid,
+ log(IC2/moyenneglobale),col="blue")
+ segments(mid-.1,log(IC1/moyenneglobale),mid+.1,
+ log(IC1/moyenneglobale),col="blue")
+ segments(mid-.1,log(IC2/moyenneglobale),mid+.1,
+ log(IC2/moyenneglobale),col="blue")
+ points(mid,log(moyenne/moyenneglobale),pch=19,col="red")
+ axis(4)
+ abline(h=0,lty=2,col="red")}
+
+ mtext("Exposition", 2, line=2, cex=1.2,col="light green")
+ if(type==1){mtext("Fréquence annualisée", 
+     4, line=-2, cex=1.2,col="red")}
+ if(type==2){mtext("Fréquence annualisée (log multiplicateur)", 
+    4, line=-2, cex=1.2,col="red")}
+ }

Par exemple en utilisant un découpage arbitraire par classe d’age (comme cela est fait par défaut dans la fonction),

> graphique()

Mais on pourrait aussi utiliser un découpage assurant d’avoir des classes plus grandes, par exemple en utilisant les quantiles,

> Q=quantile(basenb[,"ageconducteur"],(0:10)/10)
> Q[1]=Q[1]-1
> graphique(nom="ageconducteur",niveau=Q,continu=TRUE)

Ces deux graphiques permettent de visualiser la fréquence empirique de sinistre, par classe d’age, sans modèle paramétrique. Des intervalles de confiance sont également représentés (basés sur une hypothèse de normalité). Notons que l’on peut aussi faire une représentation relative à la fréquence moyenne (en log des multiplicateurs),

> graphique(type=2)

ou encore, si on cherche à analyser la fréquence en fonction de la zone géographique d’habitation

> graphique(nom="zone",continu=FALSE,type=2)

Une autre utilisation peut être faite sur la sévérité des sinistres (coût moyen) et la fréquence, par exemple un assureur,

Il est possible de modifier un peu la fonction pour ajouter au graphique la sévérité des sinistres, e.g.

> graphiquecout=function(nom="ageconducteur",
+ niveau=c(17,21,24,29,34,44,64,84,110),
+ continu=TRUE,type=1){
+ if(continu==TRUE){X=cut(basenb[,nom],niveau)}
+ if(continu==FALSE){X=basenb[,nom]}
+ E=basenb$exposition
+ Y=basenb$nbre
+ FREQ=levels(X)
+ moyennen=variancen=nn=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyennen[k] =weighted.mean(Y[X==FREQ[k]]/E[X==FREQ[k]],
+ E[X==FREQ[k]])
+ variancen[k]=weighted.mean((Y[X==FREQ[k]]/E[X==FREQ[k]]-
+ moyennen[k])^2,E[X==FREQ[k]])
+ nn[k]       =sum(E[X==FREQ[k]])
+ }
+ moyenneglobalen=sum(Y)/sum(E)
+ 
+ if(continu==TRUE){X=cut(basesin[,nom],niveau)}
+ if(continu==FALSE){X=basesin[,nom]}
+ Y=basesin$cout
+ FREQ=levels(X)
+ moyennes=variances=ns=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyennes[k] =mean(Y[X==FREQ[k]])
+ variances[k]=var(Y[X==FREQ[k]])
+ ns[k]=length(Y[X==FREQ[k]])
+ }
+ moyenneglobales=mean(Y)
+ 
+ w=barplot(nn,names.arg=FREQ,col="light green",
+ axes=FALSE,xlim=c(0,1.2*length(FREQ)+.5))
+ mid=w[,1]
+ 
+  par(new=TRUE)
+ IC1=moyennen+1.96/sqrt(nn)*sqrt(variancen)
+ IC2=moyennen-1.96/sqrt(nn)*sqrt(variancen)
+ plot(mid,moyennen,ylim=range(c(IC1,IC2)),type="b",
+ col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="red")
+ segments(mid-.1,IC1,mid+.1,IC1,col="red")
+ segments(mid-.1,IC2,mid+.1,IC2,col="red")
+ points(mid,moyennen,pch=19,col="red")
+ axis(4)
+ abline(h=moyenneglobalen,lty=2,col="red")
+ 
+ par(new=TRUE)
+ IC1=moyennes+1.96/sqrt(ns)*sqrt(variances)
+ IC2=moyennes-1.96/sqrt(ns)*sqrt(variances)
+ plot(mid,moyennes,ylim=range(c(IC1,IC2)),type="b",
+ col="blue",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="blue")
+ segments(mid-.1,IC1,mid+.1,IC1,col="blue")
+ segments(mid-.1,IC2,mid+.1,IC2,col="blue")
+ points(mid,moyennes,pch=19,col="blue")
+ axis(2)
+ abline(h=moyenneglobales,lty=2,col="blue")
+ 
+ mtext("Cout moyen", 2, line=2, cex=1.2,col="blue")
+ mtext("Fréquence annualisée", 4, line=-2, cex=1.2,col="red")
+ }

Si on regarde par classe d’age,

> graphiquecout()

ou encore, en fonction de l’age du véhicule,

> Q=quantile(basenb[,"agevehicule"],(0:10)/10)
> Q[1]=Q[1]-1
> graphiquecout(nom="agevehicule",niveau=Q,continu=TRUE)

ou enfin, en fonction de la zone géographique

> graphiquecout(nom="zone",continu=FALSE)

Actuarial science, Maths, Premium principle, Probability, Publications, Ratemaking

Loi des grands nombres et théorème central limite en assurance

17/06/2011 Arthur Charpentier 4 Comments

Dans le dernier numéro de la revue Risques figure un petit article sur la loi des grands nombres et le théorème central limite en assurance (ici) Le but était de remettre un peu les pendules à l’heure sur la différence entre les implications de ces deux résultats, et surtout d’amorcer une discussion sur l’importance des hypothèses (en particulier celle de risques indépendants).

Actuarial science, Premium principle

Market valuation des garanties des produits d’assurance

08/01/2010 Arthur Charpentier Leave a comment

Pour le dernier cours d’actuariat 2, on parlera de market valuation, i.e. de valorisation financière. Je vais aussi en profiter pour éclairer Bruno (pour qui j’avais déjà commencé à poster des billets sur le contrôle optimal) sur “comment introduit-t-on la mesure risque neutre dans un passif ?” , sous-entendu en assurance-vie (en assurance non-vie c’est à mon avis plus compliqué, mais c’est la question que l’on se pose quand on valoriser des cat-bonds par exemple). Notons que ce problème s’est posé assez tôt dans la littérature économique, je renverrais par exemple au papier de Brennan et Schwartz en 1976.

Valorisation de produits financiers (en marché complet)

Je passerais un peu rapidement (en renvoyant ici ou là pour les grands principes), mais les techniques principales pour valoriser les actifs financiers sont les prix d’Arrow-Debreu, les probabilités risque-neutre, et les déflateurs. Mais auparavant, rappelons que le marché est composé de n actifs financiers (bien distincts) dont les prix et les payoffs sont parfaitement connus, et de n états de la natures, de telle sorte que la matrice des payoffs soit de plein rang. Introduisons le concept de portefeuillle de réplication. Compte tenu du fait que la matrice des payoffs est de plein rang, on a l’existence et l’unicité d’un tel portefeuille. Considérons un ensemble de payoff contingents aux états de la nature

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-17.png$

Il est alors possible de créer un portefeuille (dit de réplication) qui permette de toucher précisément ces payoffs à la date 1. Le prix (à la date 0) d’un tel actif est alors le prix du portefeuille de réplication.
Considérons les actifs rapportant 1 dans un état de la nature, et 0 sinon, à la date 1, i.e.

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-16.png$

appelés actifs d’Arrow-Debreu. On note $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-15.png$ le prix aujourd’hui d’un tel actif. Notons que le prix d’un actif sans risque, rapportant 1 demain dans tous les états du monde est $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-1.png$ . Aussi, si on note $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-14.png$ le taux sans risque, alors

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/correc.png$

Le prix de l’actif rapportant $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-12.png$ à la date 1 est alors $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-10.png$ .
Posons $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-09.png$ tel que $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-08.png$ soit une distribution de probabilité. Alors

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-09.png$

où l’espérance est prix sous la mesure de probabilité $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-pst.png$ appelée probabilité risque-neutre.
Je renvoie ici (par exemple) pour l’introdution de cette probabilité via les arbres binomiaux. Bref, c’est ainsi que l’on arrive au théorème fondamental de la valorisation d’actif en marché complet: le prix d’un actif contingent l’espérance sous probabilité risque neutre, des payoff actualisés.
Aussi, prenons l’exemple d’un put européen, de maturité T et de strike K. Le prix (à la date 0) du put est donné par

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-01.png$

Rappelons d’ailleurs que la formule de Black et Scholes donne un prix pour un tel produit financier (un put européen),

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/BSC01.png$

où

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/BSC02.png$

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/BSC03.png$

Voilà rapidement comment on valorise un produit financier en marché complet, en calculant la valeur d’un portefeuille de réplication.

Valorisation de produits d’assurance

Pour valoriser les produits d’assurance, classiquement, on calcule des primes pures, ou des valeurs actuelles probables. Par exemple en assurance-vie (car c’est essentiellement de ça dont on parle non ?), considérons simplement une garantie plancher, où l’on s’engage à verser 1€ au décès d’un souscripteur. La valeur actuelle probable (que l’on calcule sous la probabilité historique) est alors

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-02.png$

que l’on pourrra réécrire

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-03.png$

(où j’adopte une convention de temps continu pour les taux pour garder une analogie avec le cas financier, mais ça n’est pas vraiment important ici). Bref, par rapport à la valorisation financière on valorise toujours des flux futurs (payoffs) actualisés, mais en prenant l’espérance sous la probabilité historique. Et c’est généralement comme ça que l’on valorise un produit d’assurance: les risques ne sont pas échangés sur des marchés liquides, on n’a donc pas vraimentt de marché complet, et les assureurs espèrent simplement que la mutualisation (via la loi des grands nombres) leur permettra d’avoir égalité entre les flux entrants et les flux sortants.
Remarque: j’avais fait l’an dernier un rapide billet sur la transformée d’Esscher (ici) qui présentait une méthode conciliant les deux approches. Mais je suis toujours à la recherche du papier original (datant de 1932) pour proposer une nouvelle discussion sur ce sujet…

Valorisation des garanties planchers dans les contrats d’assurance vie, une approche heuristique

Là où les choses se compliquent c’est quand on mélange des garanties d’assurance (un décès) avec des garanties financières (sur des taux ou bien un actif financier pour des contats en unité de compte). Considérons pour cela un contrat décès, où l’assureur s’engage à payer, au moment du décès à la date $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-tau.png$ le maximum entre la somme investie en actif risqué à la signature du contrat, $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-07.png$ et le niveau atteint par l’actif lors du décès, $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-06.png$ . Autrement dit, l’assureur s’engage à verser $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-05.png$ . De manière plus général, disons qu’il existe une valeur plancher en dessous de laquelle on ne descendra pas, que l’on notera K.
Notons que l’on peut réécrire le payoff

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-04.png$

Le terme de gauche est la valeur de l’actif, et le second terme peut être interprété comme la valeur de la garantie plancher. Par la suite, on ne s’intéressera qu’à la valeur de cette garantie. Le flux versé (pour cette garantie) à la date t, pour un individu d’âge x à la souscription est alors tout simplement

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-04.png$

(les dates étant supposées discrètes ici). Si l’on cherche à calculer la valeur actuelle probable de l’ensemble des flux, on voit que deux espaces de probabilités intervienne: l’espace des probabilités actuarielles (sous lesquelles on couvre le risque de décès en supposant les risques mutualisables, i.e. l’univers historique $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-05.png$ ) et l’espace des probabilités financières (sous lesquelles on couvre le risque financiers, i.e. l’univers risque neutre $https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-06.png$ ). Autrement dit, la valeur actuelle probable s’écrit

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-07.png$

que l’on pourrrait réécrire

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-08.png$

où le terme à droite est le prix d’un call de strike K de maturité t.
Notons que la commission de contrôle s’était penché sur le sujet il y a quelques années, via le mémoire d’actuariat de Sylvain Merlus et Olivier Pequeux (ici), et c’est ce genre de formule qui avait été proposé, à savoir de faire somme somme de prix de put. Et c’est cette approche que la réglementation a imposé de valoriser cette garantie.
Pour aller plus loin, il avait aussi été précaunisé d’utiliser la formule de Black et Scholes pour valoriser ces puts, ce qui a donné la formule “standard”

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-03.png$

où

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-02.png$

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-01.png$

Voilà comment on valorise une telle garantie. Je n’ai parlé de la couverture, car on voit qu’elle ne serait pas parfaite, loin de là.
Comme toutes les formules standards, elle a l’avantage de la simplicité.
Bon formellement, si on veut faire ça proprement, ça risque d’êetre plus compliqué, et il faudrait parler des mesures risques neutre forward, mais on fera ça une autre fois. Sinon tout ça est repris dans une multitude de documents que l’on peut trouver ici, là, voire ici ou encore là. Et pour aller plus loin, il y a eu récemment pas mal de livres qui abordent ces sujets….

En complément, je pourrais aussi mettre des vieux transparents d’un cours donné il y a plusieurs années maintenant,

Slides titrisation from Arthur Charpentier

Actuarial science, Copulas, Premium principle, Stochastic Processes

Hypothèse d’indépendance en assurance

05/06/2009 Arthur Charpentier Leave a comment

Il y a deux semaines, j’ai reçu un mail dans ma boite à l’ensae (que je consulte rarement) qui posait une vraie question. “Pour le modèle de risque collectif : S=X1+X2+…+XN. On suppose que les Xi et N sont indépendants (montant de sinistres et nombre de sinistres). Question : Dans la pratique (par exemple : assurance automobile), comment peut-on expliquer cette indépendance.“

Sur la forme de la question, on mélange un peu théorie et pratique: c’est la théorie qui impose cette hypothèse, et effectivement, il convient de la confronter à la pratique afin d’utiliser les résultats qui découlent du modèle collectif (en particulier en terme de solvabilité, j’en parlerai à la fin). Mais cela pose toutefois une question très intéressante.

Pour rappels, cette hypothèse a été faite par Filip Lundberg quand il a proposé les premiers résultats en 1903 (voire plus tôt), à l’époque où l’étude des processus était assez peu développée… Charles Spearman n’a introduit sa mesure de corrélation qu’en 1904, par exemple (même si on peut trouver des choses chez Francis Galton en 1888, puis chez Karl Pearson en 1892¹). Autrement dit, le premier modèle reposait sur les hypothèses les plus simples possibles (pour l’époque). Depuis, l’étude des processus ponctuels s’est beaucoup développée, et énormément de résultats ont été obtenus. Autrement dit, il existe des résultats théoriques permettant de s’affranchir de cette indépendance.

Petite remarque en passant, dans les livres avec Michel Denuit, on avait insisté longuement sur cette hypothèse d’indépendance… et c’est précisément sur ce point que Claude Bébéar avait ouvert la préface (ici) et qu’Hans Bühlmann avait clôturé la postface (là). C’est d’ailleurs Hans Bühlmann qui notait en 1963 que “the independence hypothesis is so common to be made that many authors forget to mention it“.
Pour être un peu complet, il y a trois indépendances distinctes dans les hypothèses du modèle collectif:

d’indépendance entre les arrivées des sinistres, on parle plus proprement d’hypothèse d’un processus de Poisson homogène (les durées entre survenances suivent des lois exponentielles indépendantes),
d’indépendance entre les coûts (individuels) des sinistres,
d’indépendance entre les coûts de sinistres, et les arrivées.

C’est surtout la première hypothèse qui a été relâchée. On retrouvera des exemples dans les papiers suivants

Malinovskii, V.K. (1998). Non-Poissonian claims’ arrivals and calculation of the probability of ruin. Insurance: Mathematics and Economics, 22, 2, 123-138.
Gyllenberg, M. & Silvestrov, D.S. (2000). Cramer-Lundberg approximation for nonlinearly perturbed risk processes. Insurance: Mathematics and Economics, 26, 1, 75-90.

Sur les aspects pratiques, on peut retrouver ce phénomène en particulier sur les catastrophes naturelles qui sont très cycliques. L’exemple le plus connu et le plus étudié est les ouragans,

Parisi, F. & Lund, R. (2000). Seasonality and return periods of landfalling Atlantic basin hurricanes. Australian and New Zealand Journal of Statistics, 42, 3, 271-282.

Mais on a la même chose en assurance auto, à cause de la pluie et du brouillard (en France en tous cas²).
Pour la seconde, on peut trouver des choses dans les derniers papiers de Stéphane Loisel (ici).
Mais c’est le dernier point qui était évoqué, et seulement celui là… Avec la même dichotomie que tout à l’heure, on trouve des résultats théoriques, en particulier dans les papiers d’Hansjorg Albrecher

Albrecher, H. & Boxma, O.J. (2004). A ruin model with dependence between claim sizes and claim intervals. Insurance: Mathematics and Economics, 35, 2, 245-254.

Sur les aspects pratiques, l’exemple le plus classique et le plus étudié est celui des tremblements de terre. Pour faire très simple, un tremblement de terre est causé par une accumulation d’énergie: plus on attend, plus le tremblement de terre risque d’être important. Bref, d’un point de vue pratique, cette hypothèse peut être délicate à supposer. On avait utilisé un modèle ACD avec un processus marqué dans le cas des crues de rivières par exemple, dans un papier avec David Sibaï (ici). On modélisait un processus joint prenant en compte les durées des crues, les durées entre crue, et les importances des crues. Là aussi, intensités et durées ne sont pas indépendant.

Mais au delà de l’explication, il faut aussi – peut être – essayer de tester cette indépendance. Si quelqu’un a des données ou des résultats sur ce sujet, je suis preneur.

Pour conclure, je dirais que cette hypothèse d’indépendance est fondamentale. Si on la supprime, on peut perdre beaucoup de résultats classiques, sur la solvabilité par exemple. Pour reprendre une idée commune partagée par beaucoup de monde, dans l’Encyclopaedia of Financial Engineering and Risk Management, (ici) on apprend que “reinsurance is able to offer additional underwriting capacity for cedants, but also to reduce the probability of a direct insurer’s ruin“. En fait, ce résultat devient faux si l’on supprime cette hypothèse d’indépendance: on peut augmenter la probabilité de ruine en se réassurance ! La réassurance diminue alors la solvabilité de la compagnie d’assurance (ou tout du moins peut le faire). C’est dans un papier que j’ai soumis récemment, je garde donc tout ça pour un billet dans les mois à venir… à suivre donc.

¹ Sur la corrélation de Pearson – que j’ai beaucoup critiqué sur ce blog – je renvoie à un très beau papier que je viens (seulement) de découvrir, de Joseph Rodgers et Alan Nicewander, intitulé Thirteen Ways to Look at the Correlation Coefficient (ici)
² Si quelqu’un a des données pour illustrer ce point, je suis preneur !

Freakonometrics