# Fairness and discrimination, PhD Course, #2 Insurance and risk classes

For the second course, we will get back a little bit on insurance pricing in a context of heterogeneous portfolio, and risk classification (slides are still online on the github repository). The starting point will be the pure premium.

See our online textbook, with Michel Denuit, Non Life Insurance Mathematics, for additional motivation. If we have some risk related variables $\boldsymbol{x}=(x_1,\cdots,x_k)$, the pure premium will be the conditional expectation,

Here also, we have some law of numbers, for the conditional expected value,

This relationship, which defines the conditional expected value using the limiting value of a conditional frequency cannot be used to define properly $\mathbb{P}[Y|\boldsymbol{X}=\boldsymbol{x}]$ and $\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}]$. One can consider a limit,$$\mathbb{P}\big(Y\in \mathcal{A}\big\vert X = x\big)=\lim_{\epsilon\to0}\frac{\mathbb{P}(\{Y\in \mathcal{A}\}\cap\{|X -x|\leq \epsilon\})}{\mathbb{P}(\{|X -x|\leq \epsilon\})}$$or$$\mathbb{P}\big(Y\in \mathcal{A}\big\vert X = x\big)=\lim_{\epsilon\to0}\mathbb{P}\big(Y\in \mathcal{A}\big\vert |X -x|\leq \epsilon\big)$$as in the law of the unconscious statistician or as Proschan and Presnell (1998) wrote it

statisticians make liberal use of conditioning arguments to shorten what would otherwise be long proofs

We can now compute conditional frequency, given some risk characteristics, for some quantity of interest $y$, such as the age of death, in life insurance contracts.

#### Demographic risk and heterogeneity

First, we will see some gender-based life tables, starting with the one obtained by Nicolaas Struyck (see e.g. Alberts et al. (2014))

More recently, in France, some wealth based life tables were obtained, with various quantiles

And finally, we will see some life tables obtained 50 years ago in the US, with racial distinction

#### Mean and variance decomposition

About pure premiums, an important property is the law of total expectations, and a desirable property, that we will name “balance property”

We will also mention variance and variance decomposition, depending if we take into heterogeneity, or not. With homogenous pricing, we have

If we use the “true” underlying risk factor, $\Theta$, we have the standard variance decomposition, also called law of total variance

i.e.

And finally, if we do not observe $\Theta$, but we have a collection of covariates, $\boldsymbol{X}=(X_1,\cdots,X_k)$,

#### Some historical perspectives

In the textbook, Insurance: Biases, Discrimination and Fairness, I have several paragraph about an historical perspective, starting with insurance as clubs, without segmentation. Then segmentation started, with risk classes and groups. For example, according to Issues And Needed Improvements In State Regulation Of The Insurance Business, by Harry Havens, in 1979,

The price which a person pays for automobile insurance depends on age, sex, marital status, place of residence and other factors. This risk classification system produces widely differing prices for the same coverage for different people. Questions have been raised about the fairness of this system, and especially about its reliability as a predictor of risk for a particular individual. While we have not tried to judge the propriety of these groupings, and the resulting price differences, we believe that the questions about them warrant careful consideration by the State insurance departments. In most States the authority to examine classification plans is based on the requirement that insurance rates are neither inadequate, excessive, nor unfairly discriminatory. The only criterion for approving classifications in most States is that the classifications be statistically justified — that is, that they reasonably reflect loss experience. Relative rates with respect to age, sex, and marital status are based on the analysis of national data. A youthful male driver, for example, is charged twice as much as an older driver all over the country} (…) t has also been claimed that insurance companies engage in redlining – the arbitrary denial of insurance to everyone living in a particular neighborhood. Community groups and others have complained that State regulators have not been diligent in preventing redlining and other forms of improper discrimination that make insurance unavailable in certain areas. In addition to outright refusals to insure, geographic discrimination can include such practices as: selective placement of agents to reduce business in some areas, terminating agents and not renewing their book of business, pricing insurance at un-affordable levels, and instructing agents to avoid certain areas. We reviewed what the State insurance departments were doing in response to these problem. To determine if redlining exists, it is necessary to collect data on a geographic oasis. Such data should include current insurance policies, new policies being written, cancellations, and non-renewals. It is also important to examine data on losses by neighborhoods within existing rating territories because marked discrepancies within territories would cast doubt on the validity of territorial boundaries. Yet, not even a fifth of the States collect anything other than loss data, and that data is gathered on a territory-wide basis.

According to The Role of Risk Classification in Property and Casualty Insurance: A Study of the Risk Assessment Process : Final Report, by Barbara Casey, Jacques Pezier and Carl Spetzler, in 1976,

On the other hand, the opinion that distinctions based on sex, or any other group variable, necessarily violate individual rights reflects ignorance of the basic rules of logical inference in that it would arbitrarily forbid the use of relevant information. It would be equally fallacious to reject a classification system based on socially acceptable variables because the results appear discriminatory. For example, a classification system may be built on use of car, mileage, merit rating, and other variables, excluding sex. However, when verifying the average rates according to sex one may discover significant differences between males and females. Refusing to allow such differences would be attempting to distort reality by choosing to be selectively blind. The use of rating territories is a case in point. Geographical divisions, however designed, are often correlated with socio-demographic factors such as income level and race because of natural aggregation or forced segregation according to these factors. Again we conclude that insurance companies should be free to delineate territories and assess territorial differences as well as they can. At the same time, insurance companies should recognize that it is in their best interest to be objective and use clearly relevant factors to define territories lest they be accused of invidious discrimination by the public. (…) One possible standard does exist for exception to the counsel that particular rating variables should not be proscribed. What we have called equal treatment’ standard of fairness may precipitate a societal decision that the process of differentiating among individuals on the basis of certain variables is discriminatory and intolerable. This type of decision should be made on a specific, statutory basis. Once taken, it must be adhered to in private and public transactions alike and enforced by the insurance regulator. This is, in effect, a standard for conduct that by design transcends and preempts economic considerations. Because it is not applied without economic cost, however, insurance regulators and the industry should participate in and inform legislative deliberations that would ban the, use of particular rating variables as discriminatory.

And then, more recently, we started to talk about personalization, as in Barry and Charpentier (2020). And next week, we will start talking about predictive modeling, and machine learning.

# INF7100, statistiques

La seconde partie de mon intervention sur la science des données, dans le cadre du cours INF7100 portera sur les statistiques, univariées et multivarirées. Le plan sera le suivant

• 201: De la Statistique aux Sciences de Données pdf video (14:24)
• 211: Fonctions Usuelles en Statistique (fonction de réparition, densité, histogramme) pdf video (28:37)
• 221: Indicateurs Statistique: Valeur Centrale (moyenne) pdf video (32:56)
• 222: Indicateurs Statistique: Dispersion (variance, inégalités) pdf video (22:21)
• 223: Indicateurs Statistique: Approximations (approximation normale) pdf video (18:42)
• 224: Indicateurs Statistique: Quantiles pdf video (24:54)
• 231: Inférence (statistique bayésienne) pdf video (39:33)
• 241: Tests Statistiques (1) (tests, significativité, p-value) pdf video (43:41)
• 242: Tests Statistiques (2) (erreurs) pdf video (16:51)
• 261: Statistiques Bivariées pdf video (25:16)
• 271: Statistiques Multivariées: Projections pdf video (29:06)
• 272: Statistiques Multivariées: Clusters pdf video (32:21)
• 281: Réseaux et Graphs pdf video (32:40)
• 291: Données Chronologiques pdf video (29:01)

Usually, when I give a course on GLMs, I try to insist on the fact that the link function is probably more important than the distribution. In order to illustrate, consider the following dataset, with 5 observations

x = c(1,2,3,4,5) y = c(1,2,4,2,6) base = data.frame(x,y)

Then consider several model, with various distributions, and either an identity link (and in that case $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbf{x}^T\mathbf{\beta}$) or a log link function (so that $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=e^{\mathbf{x}^T\mathbf{\beta}}$)

regNId = glm(y~x,family=gaussian(link="identity"),data=base) regNlog = glm(y~x,family=gaussian(link="log"),data=base) regPId = glm(y~x,family=poisson(link="identity"),data=base) regPlog = glm(y~x,family=poisson(link="log"),data=base) regGId = glm(y~x,family=Gamma(link="identity"),data=base) regGlog = glm(y~x,family=Gamma(link="log"),data=base) regIGId = glm(y~x,family=inverse.gaussian(link="identity"),data=base) regIGlog = glm(y~x,family=inverse.gaussian(link="log"),data=base

One can also consider some Tweedie distribution, to be even more general

library(statmod) regTwId = glm(y~x,family=tweedie(var.power=1.5,link.power=1),data=base) regTwlog = glm(y~x,family=tweedie(var.power=1.5,link.power=0),data=base)

Consider the prediction obtained in the first case, with the linear link function

library(RColorBrewer) darkcols = brewer.pal(8, "Dark2") plot(x,y,pch=19) abline(regNId,col=darkcols[1]) abline(regPId,col=darkcols[2]) abline(regGId,col=darkcols[3]) abline(regIGId,col=darkcols[4]) abline(regTwId,lty=2)

The predictions are very very close, aren’t they ? In the case of the exponential prediction, we obtain

plot(x,y,pch=19) u=seq(.8,5.2,by=.01) lines(u,predict(regNlog,newdata=data.frame(x=u),type="response"),col=darkcols[1]) lines(u,predict(regPlog,newdata=data.frame(x=u),type="response"),col=darkcols[2]) lines(u,predict(regGlog,newdata=data.frame(x=u),type="response"),col=darkcols[3]) lines(u,predict(regIGlog,newdata=data.frame(x=u),type="response"),col=darkcols[4]) lines(u,predict(regTwlog,newdata=data.frame(x=u),type="response"),lty=2)

We can actually look closer. For instance, in the linear case, consider the slope obtained with a Tweedie model (that will include all the parametric familes mentioned here, actually)

pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base))$coefficients[2,1:2] Vgamma = seq(-.5,3.5,by=.05) Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3,ylim=c(.965,1.03),xlab="power",ylab="slope") The slope here is always very very close to one ! Even more if we add a confidence interval plot(Vgamma,Vpente[1,]) lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2) Heuristically, for the Gamma regression, or the Inverse Gaussian one, because the variance is a power of the prediction, if the prediction is small (here on the left), the variance should be small. So, on the left of the graph, the error should be small with a higher power for the variance function. And that’s indeed what we observe here erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(-.1,.04),xlab="power",ylab="error") abline(h=0,lty=2) Of course, we can do the same with the exponential models pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base))$coefficients[2,1:2] Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3)

or, if we add the confidence bands, we obtain

plot(Vgamma,Vpente[1,],ylim=c(0,.8),type="l",lwd=3,xlab="power",ylab="slope") lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2)

So here also, the “slope” is rather similar… And if we look at the error we make on the left part of the graph, we obtain

erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(.001,.32),xlab="power",ylab="error")

So my point is that the distribution is usually not the most important point on GLMs, even if chapters of books on GLMs are distribution based… But as mentioned in an another post, if you consider a nonlinear transformation, like we have with GAMs, the story is more complicated…

# Unbiased Estimators vs. Minimizing a Quadratic Loss Function

Unbiased estimators are important in statistics. I guess because of Cramér Rao bound, for the variance. In the sense that if  $\mathbb{E}[\widehat{\theta}]=\theta$, then  $\text{Var}[\widehat{\theta}]\geq I_\theta^{-1}$, where  $I_\theta$ denotes Fisher information (the proof was writen in an old post).

But what could we be the variance if  $\widehat{\theta}$ is not unbiased ?

Consider the following simple case, with a Gaussian i.id. sample  $\{X_1,\cdots,X_n\}$ from a  $\mathcal{N}(\theta,\sigma^2)$. We know that the estimator of the Method of Moments is the same as the Maximum Likelihood estimator, i.e.  $\widehat{\theta}=\overline{X}$. And this estimator is efficient, in the sense that its variance is equal to Cramér-Rao lower bound,  $\text{Var}[\widehat{\theta}]=n^{-1}\sigma^2$.

But what if we consider another estimator? For instance  $\tilde{\theta}=\alpha\cdot\overline{X}$, with  $\alpha$ not necessarily equal to 1. In that case, this estimator is (usually) biased since

$\mathbb{E}[\tilde{\theta}]=\alpha\theta$

And its variance is

$\text{Var}[\tilde{\theta}]=\alpha^2\text{Var}[\overline{X}]=\frac{\alpha^2\sigma^2}{n}$

We can visualise those two functions (the bias and the variance) using

 n=10 alpha=seq(0,2,by=.01) b=1-alpha v=alpha^2/n plot(alpha,b,xlab="alpha",col="red",type="l") par(new=TRUE) plot(alpha,v,col="blue",type="l",axes=FALSE axis(4,) mtext("bias",side=2,line=-1,col="red") mtext("variance",side=4,line=-1,col="blue") 

Observe that if $\alpha$ is small (smaller than 1), the variance is smaller than the Cramér-Rao lower bound. And here, the mean squared error, defined as

$\text{mse}[\tilde{\theta}]=\mathbb{E}[(\tilde{\theta}-\theta)^2]=\text{bias}[\tilde{\theta}]^2+\text{Var}[\tilde{\theta}]$

which is, here,

$\text{mse}[\tilde{\theta}]=[(\alpha-1)\theta]^2+\frac{\alpha^2\sigma^2}{n}$

The optimal value is obtained when the first order condition is satisfied

$\frac{\partial\text{mse}[\tilde{\theta}]}{\partial\alpha}=-2(\alpha-1)\theta^2-\frac{2\alpha\sigma^2}{n}=0$

i.e.

$\alpha=\frac{n\theta^2}{n\theta^2+\sigma^2}$

So biased estimators can be more interesting than unbiased estimators, if the goal is the minimize the mean square error.

# Heuristics on bias and variance for kernel density estimators

Consider the simple case of a moving histogram (which is a very simple kernel). The idea is to recall that

$f(x)=\lim_{h\downarrow 0} f_h(x)$

where

$f_h(x)=\frac{1}{h}\left[F\left(x+\frac{h}{2}\right)-F\left(x-\frac{h}{2}\right)\right]$

is the slope close to point $x$.

Then we use the empirical cumulative density to approximate the slope, i.e.

$\hat f_h(x)=\frac{1}{h}\left[\hat F_n\left(x+\frac{h}{2}\right)-\hat F_n\left(x-\frac{h}{2}\right)\right]$

which can also be writen

$\hat f_h(x)=\frac{1}{nh}\sum_{i=1}^n \boldsymbol{1}\left(x_i\in\left[x-\frac{h}{2},x+\frac{h}{2}\right]\right)$

Consider now the density seen as a random variable

$\hat f_h(x)=\frac{1}{nh}\sum_{i=1}^n \underbrace{\boldsymbol{1}\left(X_i\in\left[x-\frac{h}{2},x+\frac{h}{2}\right]\right)}_{Y_i}$

where the$Y_i$‘s are i.i.d. where $Y_i\sim\mathcal{B}(p_x)$, with

$p_x=\mathbb{P}\left(X_i\in\left[x-\frac{h}{2},x+\frac{h}{2}\right]\right)=h\cdot f_h(x)$

Thus, observe that $\mathbb{E}(\hat f_h(x))=f_h(x)$, but that’s not what we’re looking for… From Taylor’s expansion,

$f_h(x)=f(x)+\frac{h^2}{24}f''(x)+o(h^2)$

thus

$\mathbb{E}(\hat f_h(x))\sim f(x)+\frac{h^2}{24}f''(x)$

where the bias comes from the approximation of the density by some string. About the variance,

$\text{Var}(\hat f_h(x))=\frac{1}{n^2h^2}\cdot n p_x[1-p_x]$

thus, since $p_x=h\cdot f_h(x)=o(h)$,

$\text{Var}(\hat f_h(x))\sim\frac{nh}{n^2h^2}\cdot f_h(x)$

i.e.

$\text{Var}(\hat f_h(x))\sim\frac{1}{nh}\cdot f(x)$

We can observe that

$\text{bias}(\hat f_h(x))=\mathbb{E}(\hat f_h(x))-f(x)\sim\frac{h^2}{24}f''(x)$

is decreasing as $h \downarrow 0$, while the variance is increasing as $h \downarrow 0$. This is the standard bias-variance tradeoff in statistics.

# Variance of the Average of a Sequence

In the case where $\{Y_1,\cdots,Y_n\}$ are i.i.d. random variables, then

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\text{Var}(Y_t)}{n}$

Now, what if $\{Y_1,\cdots,Y_n\}$ are identically distributed, but no longer independent. What if we have an autoregressive process? Assume that

$Y_t=\phi Y_{t-1}+\varepsilon_t$

Then

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right) = \frac{1}{n^2}\left(\sum_{t=1}^n \text{Var}(Y_t)+\sum_{s\neq t}\text{Cov}(Y_s,Y_t) \right)$

can be written

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right) = \frac{1}{n^2}\left(n\gamma(0) + \sum_{h=1}^{n-1} 2 (n-h) \gamma(h)\right)$

Here, we will express the variance as a function of $\gamma(0)$ and $\phi$, but it is possible to use also $\sigma^2$, since, in the context of an $AR(1)$,

$\gamma(0)=\frac{\sigma^2}{1-\phi^2}$

Now, since $\gamma(h)=\phi^h \gamma(0)$ we get

$\text{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right) = \frac{\gamma(0)}{n^2}\left(n\gamma(0) + \sum_{h=1}^{n-1} 2 (n-h) \phi^h\right)$

which can be simplified, since

$\sum_{h=1}^{n-1} 2(n-h) \phi^h=2\phi^{n-1} \sum_{h=1}^{n-1} (n-h) \left(\frac 1\phi\right)^{n-h-1}=2\phi^{n-1} \frac {\partial}{\partial x}\left. \sum_{i=1}^{n} x^{n-i} \right\vert_{x=1/\phi}$

i.e.

$\sum_{h=1}^{n-1} 2(n-h) \phi^h=2\phi^{n-1} \frac{(n-1)\phi^{-(n+1)} - n\phi^{-n} + \phi^{-1}}{\phi^{-1}(\phi^{-1}-1)^2}=2 \frac{(n-1)\phi^{-1} - n + \phi^{n-1}}{(\phi^{-1}-1)^2}$

So, the variance of the mean can be writen as

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\gamma(0)}{n^2}\left[n + 2 \frac{(n-1)\phi^{-1} - n + \phi^{n-1}}{(\phi^{-1}-1)^2}\right]$

Observe that if $n$ is large enough,

$V=\frac{\gamma(0)}{n^2}\left[n + 2 \frac{(n-1)\phi^{-1} - n + \phi^{n-1}}{(\phi^{-1}-1)^2}\right]\sim \frac{\gamma(0)}{n}\frac{1+\phi}{1-\phi}$

This asymptotic relationship is well known actually. A simple way to get it is the following. One can can write

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\gamma(0)}{n}\left[ \sum_{h=-n+1}^{n-1}\left(1-\frac{\vert h\vert}{n}\right)\rho(h) \right]$

or equivalently

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)=\frac{\gamma(0)}{n}\left[ 1+2\sum_{h=1}^{n-1}\left(1-\frac{h}{n}\right)\rho(h) \right]$

But actually, the first relationship is probably more interesting to get an asymptotic approximation,

$V=\mathrm{Var}\left(\frac{1}{n}\sum_{t=1}^n Y_t\right)\sim\frac{\gamma(0)}{n}\left[ \sum_{h=-\infty}^{\infty}\rho(h) \right]$

In the context of an $AR(1)$ process, this can be writen

$\sum_{h=-\infty}^{\infty}\rho(h) = \frac{1+\phi}{1-\phi}$

Thus, we get the following well-known relationship

$V=\frac{\gamma(0)}{n}\cdot \frac{1+\phi}{1-\phi}$

In the case where $\{Y_1,\cdots,Y_n\}$ is an i.i.d. sequence, i.e. $\phi=0$, then we get the relationship mentioned initially. And in the case of a random walk… unfortunately, we cannot use that relationship. But observe that

$V=\frac{1}{n}\text{Var}\left(\sum_{t=1}^n\sum_{h=1}^t\varepsilon_h\right)$

i.e.

$V=\frac{1}{n}\text{Var}\left(n\varepsilon_1+(n-1)\varepsilon_2+\cdots+\varepsilon_n\right)$which can be written

$V=\frac{\text{Var}(\varepsilon)}{n}\sum_{h=1}^n h^2=\frac{(2n+1)(n+1)}{6n}\text{Var}(\varepsilon)$

If we compare the true value and the approximation, we get the following graph,

> V=function(phi,s2=1,n=100){
+ g0=s2/(1-phi^2)
+ if(phi<1){
+ if(phi==0){v1=g0/n}
+ if(phi>0){v1=g0/n^2*(n+2*((n-1)*
+ phi^(-1)-n+phi^(n-1))/(phi^(-1)-1)^2)}
+ v2=g0/n*(1+phi)/(1-phi)
+ }
+ if(phi==1){
+ v1=(2*n+1)*(n+1)*s2/(6*n)
+ v2=NA
+ }
+ return(c(v1,v2))}
>
> Vphi=function(phi) V(phi,1,100)
> x=seq(.01,1,by=.02)
> M=matrix(unlist(lapply(x,V)),nrow=2)
> plot(x,M[1,],type="l",col="red",log="y",
+ ylab="Variance of the average (log scale)",
+ xlab="Autoregressive coefficient")
> lines(x,M[2,],col="blue")

# Visualizing overdispersion (with trees)

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.

# Overdispersion with different exposures

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 $N_i$ on the period $[0,1]$ is unobserved
• the number of claims $Y_i$ on $[0,E_i]$ is observed (as well as $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 $(Y_i,E_i)$‘s – or more generally also some additional covariates $(Y_i,E_i,\boldsymbol{X}_i)$‘s. For ratemaking, we need to estimate $\mathbb{E}(N\vert\boldsymbol{X}=\boldsymbol{x})$ and perhaps also $\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 $\mathbb{E}(N)$ (forget about covariates as a start) is
$m_N=\frac{\sum_{i=1}^n Y_i}{\sum_{i=1}^n E_i}$
which is also the weight average of annualized individual counts
$m_N=\sum_{i=1}^n \frac{ E_i}{\sum_{i=1}^n E_i} \cdot \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 $N\sim\mathcal{P}(\lambda)$ while $Y\sim\mathcal{P}(\lambda \cdot E)$. Then, the likelihood is

$\mathcal{L}(\lambda,\boldsymbol{Y},\boldsymbol{E})=\prod_{i=1}^n \frac{e^{-\lambda E_i} [\lambda E_i]^{Y_i}}{Y_i!}$

i.e.

$\log \mathcal{L}(\lambda,\boldsymbol{Y},\boldsymbol{E}) = -\lambda \sum_{i=1}^n E_i +\sum_{i=1}^n Y_i \log[\lambda E_i] - \log\left(\prod_{i=1}^n Y_i!\right)$

The first order condition is here

$\frac{\partial}{\partial \lambda}\log \mathcal{L}(\lambda,\boldsymbol{Y},\boldsymbol{E}) = - \sum_{i=1}^n E_i +\frac{1}{\lambda}\sum_{i=1}^n Y_i =0$

which is satisfied if

$\widehat{\lambda}=\frac{\sum_{i=1}^n Y_i}{\sum_{i=1}^n E_i}$

So, we do have an estimator for the expected value, and a natural estimator for $\mathbb{E}(N\vert\boldsymbol{X}=\boldsymbol{x})$ is then (if we consider categorical covariates)
$m_{N|\boldsymbol{x}} =\frac{\sum_{i,\boldsymbol{X}_i=\boldsymbol{x}} Y_i}{\sum_ {i,\boldsymbol{X}_i=\boldsymbol{x}} E_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 $E$. For instance, if $E$ is one half, insured were observed only the first six months. Then $N=Y+Y'$ with $Y\overset{\mathcal L}{=}Y'$ ($Y$ is the number of claims on the first six months, while $Y'$ are the number of claims on the last six months), i.e. $\text{Var}(N)=\text{Var}(Y)+ \text{Var}(Y')$ if we assume independent increments. I.e.
$\text{Var}(N)=2\text{Var}(Y)$, or conversely $E \cdot\text{Var}(N)=\text{Var}(Y)$. More generally, it is reasonable to assume that

$\text{Var}(Y)=E\cdot \text{Var}(N)$
for all values of $E$. And then
$\text{Var}\left(\frac{Y}{E}\right)=\frac{1}{E}\cdot \text{Var}(N)$
Thus, it seems legitimate to assume that the empirical variance of $N$ can be written
$S_N^2=E\cdot S_{Y/E}^2$
Since the average of $Y_i/E$ is $\overline{N}=m_N$, then
$S_N^2=E\cdot \frac{1}{n}\sum_{i=1}^n \left[\frac{Y_i}{E}-\overline{N}\right]^2} = \frac{1}{n}\sum_{i=1}^n E\left[\frac{Y_i}{E}-\overline{N}\right]^2}$
or equivalently
$S_N^2=\frac{1}{n}\sum_{i=1}^n \frac{E}{E^2}\left[Y_i-\overline{N}\cdot E\right]^2} =\frac{1}{n}\sum_{i=1}^n \frac{1}{E}[Y_i-\overline{N}\cdot E]^2$i.e.
$S_N^2=\frac{\sum_{i=1}^n [Y_i-\overline{N}\cdot E]^2 }{nE}$
Thus, with different $E_i$‘s, it would be legitimate (I guess) to consider
$S_N^2=\frac{\sum_{i=1}^n [Y_i-\overline{N}\cdot E_i]^2 }{\sum_{i=1}^n E_i}$
Thus, an estimator for $\text{Var}(N|\boldsymbol{X}=\boldsymbol{x})$ is
$S_{N|\boldsymbol{x}}^2=\frac{\sum_{i,\boldsymbol{X}_i=\boldsymbol{x}} [Y_i-\overline{N}\cdot E_i]^2}{\sum_{i,\boldsymbol{X}_i=\boldsymbol{x} } E_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",
>  sinistres=sinistre[sinistre$garantie=="1RC",] > sinistres=sinistres[sinistres$cout>0,]
>  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 !

# De la difficulté de faire parler les chiffres…

Parution d’un court article intituléde la difficulté de faire parler des chiffres pour analyser la gravité des accidents de la route” dans le dernier numéro de Variance. Le numéro complet est en ligne sur http://ensae.org/…. Sinon, tous les articles de vulgarisation sont en ligne sur http://freakonometrics.hypotheses.org/….

Le code pour le premier graphique (sur les tuées) est

base=read.table(
"http://freakonometrics.free.fr/base-graph-accidents-graves.txt",
base$date=as.Date(base$date)
base$dateavant=as.Date(base$dateavant)
base$dateapres=as.Date(base$dateapres)
plot(base$date,base$compte,main="Blessés graves sur route entre 2002 et 2009",
xlab="Date",
ylab="Nombre de blessés sur la route, par jour",col="white")
points(base$dateavant,base$compteavant,col="light green")
lines(base$dateavant,base$tendanceavant,col="red",lty=2)
lines(base$dateavant,base$splinesavant,lwd=3,col="red")
points(base$dateapres,base$compteapres,col="light blue")
lines(base$dateapres,base$tendanceapres,col="red",lty=2)
lines(base$dateapres,base$splinesapres,lwd=3,col="red")

alors que pour les seconds (sur les blessés)

base=read.table(
"http://freakonometrics.free.fr/base-graph-accidents-deces.txt",
base$date=as.Date(base$date)
base$dateavant=as.Date(base$dateavant)
base$dateapres=as.Date(base$dateapres)
plot(base$date,base$compte,main="Blessés graves sur route entre 2002 et 2009",
xlab="Date",
ylab="Nombre de blessés sur la route, par jour",col="white")
points(base$dateavant,base$compteavant,col="light green")
lines(base$dateavant,base$tendanceavant,col="red",lty=2)
lines(base$dateavant,base$splinesavant,lwd=3,col="red")
points(base$dateapres,base$compteapres,col="light blue")
lines(base$dateapres,base$tendanceapres,col="red",lty=2)
lines(base$dateapres,base$splinesapres,lwd=3,col="red")`

# Estimer une variance… ou un écart-type ?

Lors du cours de demain, nous allons parler (longuement) de l’estimation de la variance (comme cela avait commencé en TD vendredi dernier). Mais en pratique, c’est souvent l’écart-type qui nous intéresse. On va passer par mal de temps à chercher un estimateur sans biais de la variance, et montrer que

devrait convenir, mais cela ne nous donne pas vraiment d’information sur le biais de l’estimateur naturel de l’écart-type, à savoir

car .
Bref, il y a sûrement un biais, et ça pourrait être intéressant de le quantifier. Classiquement, on suppose disposer d’un échantillon  i.i.d. On va faire une hypothèse forte, à savoir que leur loi est une loi normale , où  est inconnu.
Naturellement, on voudrait considérer

On va poser, afin de simplifier un peu,

et on va chercher à calculer

où  est la densité de cette variable. Or on sait que  suit une loi du chi-deux à  degrés de liberté, dont la densité est la loi gamma suivante

En notant que l’on peut réécrire

on a que cette espérance s’écrit

soit en faisant un peu de calculs,

ou encore

Ce qui va se simplifier en utilisant la fonction gamma

On a ainsi