# I Fought the (distribution) Law (and the Law did not win)

A few days ago, I was asked if we should spend a lot of time to choose the distribution we use, in GLMs, for (actuarial) ratemaking. On that topic, I usually claim that the family is not the most important parameter in the regression model. Consider the following dataset

> db <- data.frame(x=c(1,2,3,4,5),y=c(1,2,4,2,6))
> plot(db,xlim=c(0,6),ylim=c(-1,8),pch=19) To visualize a regression model, use the following code

> nd=data.frame(x=seq(0,6,by=.1))
+ prd1=predict(reg,newdata=nd,se.fit = TRUE,type="response")
+ y1=prd1$fit + y1_upp=prd1$fit+prd1$residual.scale*1.96* prd1$se.fit
+ y1_low=prd1$fit-prd1$residual.scale*1.96*
prd1$se.fit + polygon(c(nd$x,rev(nd$x)),c(y1_upp, rev(y1_low)),col="light green",angle=90, density=40,border=NA) + lines(nd$x,y1,col="red",lwd=2)
+ }

For instance, with a Poisson regression (with a log link function) we get

> plot(db)
+ data=db) while, with a Gaussian regresion (but still with a log link function), we get

> plot(db)
+ data=db)
> add_predict(reg2) If we just care about the expected value of our prediction, the output is more or less the same

> plot(db)
> lines(nd$x,predict(reg1,newdata=nd, + type="response"),col="red",lwd=1.5) > lines(nd$x,predict(reg2,newdata=nd,
+ type="response"),col="blue",lwd=1.5) So, indeed, forget about the (distribution) law when running a GLM. Not convinced? Consider – on the same dataset – a Poisson regression (with an identity link function this time)

> plot(db)
+ data=db) while, with a Gaussian regresion (but still with an identity link function), we get

> plot(db)
+ data=db)
> add_predict(reg2) Again, if we just plot the expected value of our prediction, the output is more or less the same

> plot(db)
> lines(nd$x,predict(reg1,newdata=nd, + type="response"),col="red",lwd=1.5) > lines(nd$x,predict(reg2,newdata=nd,
+ type="response"),col="blue",lwd=1.5) So clearly, the simplistic message you should not care too much about the (distribution) law seems to be valid…

# Log-transform kernel density estimation of income distribution

Our paper Log-transform kernel density estimationof income distribution, written with Emmanuel Flachaire is now available on http://papers.ssrn.com/id=2514882,

Standard kernel density estimation methods are very often used in practice to estimate density function. It works well in numerous cases. However, it is known not to work so well with skewed, multimodal and heavy-tailed distributions. Such features are usual with income distributions, defined over the positive support. We first show that a preliminary logarithmic transformation of the data, combined with standard kernel density estimation methods, can provide a much better fit of the overall density estimation. Then, we show that the fit of the bottom of the distribution may not be satisfactory, even if a better fit of the upper tail can be obtained in general.

# Exposure with binomial responses

Last week, we’ve seen how to take into account the exposure to compute nonparametric estimators of several quantities (empirical means, and empirical variances) incorporating exposure. Let us see what can be done if we want to model a binomial response. The model here is the following: ,

• the number of claims  on the period  is unobserved
• the number of claims  on  is observed (as well as )

that can be visualize below Consider the case where the variable of interest is not the number of claims, but simply the indicator of the occurrence of a claim. Then we wish to model the event versus , interpreted as non-occurrence and occurrence. Given the fact that we can only observe  versus . Having an inclusion is not enough to derive a model. Actually, with a Poisson process model, we can get easily that

With words, it means that the probability of not having a claim in the first six months of the year is the square root of not have a claim over a year. Which makes sense. Assume that the probability of not having a claim can be explained by some covariates, denoted , through some link function (using the GLM terminology),

Now, since we do observe  – and not  – we have

The dataset we will use is always the same

> 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) > sinistres = merge(contrat,nombre) > sinistres$nonsin = (sinistres$nbre==0) The first model we can consider is based on the standard logistic approach, i.e. That’s nice, but difficult to handle with standard functions. Nevertheless, it is always possible to compute numerically the maximum likelihood estimator of given . > Y=sinistres$nonsin
> X=cbind(1,sinistres$ageconducteur) > E=sinistres$exposition
> logL = function(beta){
+ 	pi=(exp(X%*%beta)/(1+exp(X%*%beta)))^E
+ 	-sum(log(dbinom(Y,size=1,prob=pi)))
+ }
> optim(fn=logL,par=c(-0.0001,-.001),
+ method="BFGS")
$par  2.14420560 0.01040707$value
 7604.073
$counts function gradient 42 10$convergence
 0
$message NULL > parametres=optim(fn=logL,par=c(-0.0001,-.001), + method="BFGS")$par

Now, let us look at alternatives, based on standard regression models. For instance a binomial-log model. Because the exposure appears as a power of the annual probability, everything would be fine if  was the exponential function (or  was the log link function), since

Now, if we try to code it, it starts quickly to be problematic,

> reg=glm(nonsin~ageconducteur+offset(exposition),
Error: no valid set of coefficients has been found: please supply starting values

I tried (almost) everything I could, but I could not get rid of that error message,

> startglm=c(0,-.001)
> names(startglm)=c("(Intercept)","ageconducteur")
> etaglm=rep(-.01,nrow(sinistresI))
> etaglm[sinistresI$nonsin==0]=-10 > muglm=exp(etaglm) > reg=glm(nonsin~ageconducteur+offset(exposition), + data=sinistresI,family=binomial(link="log"), + control = glm.control(epsilon=1e-5,trace=TRUE,maxit=50), + start=startglm, + etastart=etaglm,mustart=muglm) Deviance = NaN Iterations - 1 Error: no valid set of coefficients has been found: please supply starting values So I decided to give up. Almost. Actually, the problem comes from the fact that is closed to 1. I guess everything would be nicer if we could work with probability close to 0. Which is possible, since where is close to 0. So we can use Taylor’s expansion, Here, the exposure does no longer appears as a power of the probability, but appears multiplicatively. Of course, there are higher order terms. But let us forget them (so far). If – one more time – we consider a log link function, then we can incorporate the exposure, or to be more specific, the logarithm of the exposure. > regopp=glm((1-nonsin)~ageconducteur+offset(log(exposition)), + data=sinistresI,family=binomial(link="log")) which now works perfectly. Now, to see a final model, perhaps we should get back to our Poisson regression model since we do have a model for the probability that . > regpois=glm(nbre~ageconducteur+offset(log(exposition)), + data=sinistres,family=poisson(link="log")) We can now compare those three models. Perhaps, we should also include the prediction without any explanatory variable. For the second model (actually, it does run without any explanatory variable), we run > regreff=glm((1-nonsin)~1+offset(log(exposition)), + data=sinistres,family=binomial(link="log")) so that the prediction is here > exp(coefficients(regreff)) (Intercept) 0.06776376 This value is comparable with the logistic regression, > logL2 = function(beta){ + pi=(exp(beta)/(1+exp(beta)))^E + -sum(log(dbinom(Y,size=1,prob=pi)))} > param=optim(fn=logL2,par=.01,method="BFGS")$par
> 1-exp(param)/(1+exp(param))
 0.06747777

But is quite different from the Poisson model,

> exp(coefficients(glm(nbre~1+offset(log(exposition)),
(Intercept)
0.07279295

Let us produce a graph, to compare those models,

> age=18:100
> yml1=exp(parametres+parametres*age)/(1+exp(parametres+parametres*age))
> plot(age,1-yml1,type="l",col="purple")
> yp=predict(regpois,newdata=data.frame(ageconducteur=age,
+ exposition=1),type="response")
> yp1=1-exp(-yp)
> ydl=predict(regopp,newdata=data.frame(ageconducteur=age,
+ exposition=1),type="response")
> plot(age,ydl,type="l",col="red")
> lines(age,yp1,type="l",col="blue")
> lines(age,1-yml1,type="l",col="purple")
> abline(h=exp(coefficients(regreff)),lty=2) Observe here that the three models are quite different. Actually, with two models, it is possible to run more complex regression, e.g. with splines, to visualize the impact of the age on the probability of having – or not – a car accident. If we compare the Poisson regression (still in red) and the log-binomial model, with Taylor’s expansion, we get The next step is to see how to incorporate the exposure in a tree. But that’s another story…

# Modélisation et prévision, cas d’école Quelques lignes de code que l’on reprendra au prochain cours, avec une transformation en log, et une tendance linéaire. Considérons la recherche du mot clé headphones, au Canada, la base est en ligne sur l’ancien blog, à l’adresse freakonometrics.blog.free.fr/…

> report=read.table(
> source("http://freakonometrics.blog.free.fr/public/code/H2M.R")
> plot(headphones) Mais le modèle linéaire ne devrait pas convenir, car la série explose,

> n=length(headphones)
> X1=seq(12,n,by=12)
> X2=seq(6,n,by=12)
> points(time(headphones)[X2],Y2,pch=19,col="blue") Il est alors naturel de prendre le logarithme de la série,

> plot(headphones,log="y") C’est cette série que l’on va modéliser (mais c’est bien entendu la première série, au final, qu’il faudra prévoir). On commence par ôter la tendance (ici linéaire)

> X=as.numeric(headphones)
> Y=log(X)
> n=length(Y)
> T=1:n
> B=data.frame(Y,T)
> reg=lm(Y~T,data=B)
> plot(T,Y,type="l")
> lines(T,predict(reg),col="purple",lwd=2) On travaille alors sur la série résiduelle.

> Z=Y-predict(reg)
> acf(Z,lag=36,lwd=6)
> pacf(Z,lag=36,lwd=6) On peut tenter de différencier de manière saisonnière,

> DZ=diff(Z,12)
> acf(DZ,lag=36,lwd=6)
> pacf(DZ,lag=36,lwd=6) On ajuste alors un processus ARIMA, sur la série différenciée,

> mod=arima(DZ,order=c(1,0,0),
+ seasonal=list(order=c(1,0,0),period=12))
> mod

Coefficients:
ar1     sar1  intercept
0.7937  -0.3696     0.0032
s.e.  0.0626   0.1072     0.0245

sigma^2 estimated as 0.0046:  log likelihood = 119.47

Mais comme c’est la série de base qui nous intéresse, on utilise une écriture SARIMA,

> mod=arima(Z,order=c(1,0,0),
+ seasonal=list(order=c(1,1,0),period=12))

On fait alors la prévision de cette série.

> modpred=predict(mod,24)
> Zm=modpred$pred > Zse=modpred$se

On utilise aussi le prolongement de la tendance linéaire,

> tendance=predict(reg,newdata=data.frame(T=n+(1:24)))

Pour revenir enfin à notre série initiale, on utilise les propriétés de la loi lognormales, et plus particulièrement la forme de la moyenne, pour prédire la valeur de la série,

> Ym=exp(Zm+tendance+Zse^2/2)

Graphiquement, on a

> plot(1:n,X,xlim=c(1,n+24),type="l",ylim=c(10,90))
> lines(n+(1:24),Ym,lwd=2,col="blue") Pour les intervalles de confiance, on peut utiliser les quantiles de la loi lognormale,

> Ysup975=qlnorm(.975,meanlog=Zm+tendance,sdlog=Zse)
> Yinf025=qlnorm(.025,meanlog=Zm+tendance,sdlog=Zse)
> Ysup9=qlnorm(.9,meanlog=Zm+tendance,sdlog=Zse)
> Yinf1=qlnorm(.1,meanlog=Zm+tendance,sdlog=Zse)
> polygon(c(n+(1:24),rev(n+(1:24))),
+ c(Ysup975,rev(Yinf025)),col="orange",border=NA)
> polygon(c(n+(1:24),rev(n+(1:24))),
+ c(Ysup9,rev(Yinf1)),col="yellow",border=NA) # Le passage au log dans les modèles linéaires

Un billet rapide pour compléter et illustrer le passage au log dans un modèle linéaire (que l’on abordera cette semaine en cours). Le point de départ est le modèle linéaire, où on suppose que, conditionnellement à  suit une loi normale. Pour rappel, si on a une loi normale, , alors  et . Les intervalles de confiance à 90% et 95% sont symétriques par rapport à la moyenne (qui est aussi la médiane, soit dit en passant), Dans un modèle Gaussien avec homoscédasiticité,  i.e.  alors que . On a alors les bandes de confiance suivantes, pour un modèle de régression linéaire, Bon, maintenant, que se passe-t-il si on prend l’exponentiel ? Pour la loi normale, rappelons que l’on obtient une loi lognormale, i.e. , les deux paramètres étant liés à la loi normale sous jacente, car désormais

alors que

Graphiquement, on a la loi suivante, avec les intervalles de confiance à 90% et 95% représentés ci-dessous. Le point noir est   alors que le point bleu est l’espérance de la loi lognormale. On notera que le quantile de la loi log-normale est l’exponentiel du quantile de la loi normale. En effet, si alors . En particulier,  n’est pas la moyenne de , mais la médiane (puisque  était la médiane de ).

Mais il n’est pas rare de voir utilisé un intervalle de confiance de la forme

qui est la forme classique de l’intervalle de confiance Gaussien (symétrique autour de la moyenne). Ici, on aurait les niveaux suivants Notons qu’il n’y a aucune raison ici d’avoir une probabilité d’être dans l’intervalle de confiance obtenu avec les quantiles de la loi normale.

Maintenant, si on prend l’exponentiel d’un modèle linéaire (i.e. le logarithme de la variable d’intérêt est modélisé par un modèle linéaire) on a

avec une variance (conditionnelle) qui dépend de la variable explicative

Là encore, le plus naturel est d’utiliser comme bornes de l’intervalle de confiance des quantiles associés à la loi lognormale, mais il n’est pas rare de voir utilisé des intervalles de type Gaussiens, On perd là encore en interprétation car les bornes n’ont plus rien à voir avec les quantiles.

# Transformation logarithmique de séries temporelles

Pour poursuivre une discussion amorcée en fin de cours, dans certains cas, on peut avoir l’impression que modéliser une série pourrait être compliqué,

plot(X,xlim=c(1,length(X)+20)) Mais on peut avoir l’intuition que modéliser le logarithme de la série pourrait être plus simple,

> X=log(Y)
> plot(X,xlim=c(1,length(X)+20)) On va alors tenter une modélisation par un processus ARMA de cette dernière série,

> md=arima(X,c(12,0,1))
> P=predict(md,24)
> E=P$pred > V=P$se^2

On peut alors faire une prévision sur cette série plus simple à modéliser, et visualiser cette prévision.

> temps=length(X)+1:24
> ciu=(E+2*sqrt(V))
> cil=(E-2*sqrt(V))
> polygon(c(temps,rev(temps)),c(ciu,
+ rev(cil)),col="yellow",border=NA)
> lines(temps,E,col="red",lwd=2)
> lines(temps,ciu,col="red",lty=2)
> lines(temps,cil,col="red",lty=2) Maintenant, on va devoir remonter. On va utiliser un résultat que l’on a vu sur la transformation logarithmique dans une régression: si après avoir pris le logarithme, on a un modèle simple, Gaussien, c’est que le modèle initial était log-normal. On peut alors utiliser les propriétés de la loi lognormale, dont on connait les moments à partir de ceux de la loi Gaussienne sous-jacente. Pour la prévision, on n’a pas trop le choix,

> mu=exp(E+.5*V)

Par contre, pour construire un intervalle de confiance, soit on utilise la variance de notre loi lognormale pour avoir la variance de notre processus, et on oublie cette histoire de loi lognormale pour construire un intervalle Gaussien,

> sig2=(exp(V)-1)*exp(2*E+V)
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)

ou alors on utilise le fait que comme la transformation est monotone, l’intervalle de confiance peut etre vu comme une transformation du précédant intervalle de confiance,

> ci2u=exp(E+2*sqrt(V))
> ci2l=exp(E-2*sqrt(V))

Si on compare visuellement les deux, on a dans le premier cas,

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> polygon(c(temps,rev(temps)),c(ci1u,
> rev(ci1l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci1u,col="red",lty=2)
> lines(temps,ci1l,col="red",lty=2) (qui est symétrique et centré sur notre prévision) et dans le second > plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)
> polygon(c(temps,rev(temps)),c(ci2u,
> rev(ci2l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci2u,col="red",lty=2)
> lines(temps,ci2l,col="red",lty=2)

# Des calculs en base 12 à la magie des log-rendements Un des gros avantages dans la culture anglo-saxonne, c’est que culturellement, les anglais ont été de grands amateurs du système duodécimal (en base 12). Par exemple, il y a 12 pouces dans un pied, ou le fait que jusqu’aux années 70, il y avait 12 pences dans un schilling… Et ce que j’apprécie, c’est que ce système pousse à manipuler les fractions (2, 3, 4, 6 étant des diviseurs de 12) alors que le système décimal pousse à utiliser les nombres à virgule. Et compte tenu de la proximité avec les États-Unis, même si ici le système métrique a été adopté, il reste des fonds de ces calculs en base 12.
Par exemple, à l’école à Montréal, les tables d’addition et de multiplication vont jusqu’à 12. Aussi, l’autre soir, je faisais réviser les tables (de multiplication) de 11 et 12 à mon fils, et on a passé du temps sur et . Car sa première réponse à combien font huit fois douze ? a été ben…. cent ?. J’ai bien entendu protesté (c’est un euphémisme) pour la forme, tout en sachant qu’il est naturel de croire que , ou que . N’est ce pas ce qu’on pense quand on apprend qu’un prix a augmenté de 20% avant de baisser de 20% le lendemain ? Car c’est la même erreur que l’on commet. Et mine de rien, raisonner avec des rendements, c’est pénible…. (tous ceux qui ont discuté taux d’intérêt avec un banquier le savent).
Pour rappel, on définie le rendement comme la variation relative, i.e. L’intérêt de cette formule est évident Et si on regarde ce qui se passe sur deux jours Aussi, si un jour on a un rendement de , et le lendemain un rendement de  Autrement dit, si on subit une baisse de 10% puis une hausse de 10% (ou le contraire d’ailleurs), on ne revient pas à la position initiale, car (ou ce que je voyais avec mon fils quand on notait que ). En fait, on perd de l’argent… Damned, c’est pénible ça, non ? D’autant plus que si on répète cette alternance (entre une baisse de 10% puis une hausse de 10%) sur un an, au final, on a perdu 85% de valeur… Autrement dit, avec des rendements d’espérance nulle, et symétriques, la volatilité à elle seule fait perdre de l’argent…
Si on s’intéresse maintenant aux log-rendement, cette fois, on note que (on retrouve très explicitement la croissance exponentielle des prix). En faisant des développements limités, i.e. si les rendements sont faibles, alors on observe que . En effet Aussi, si vaut 0.05, on parlera (abusivement) de hausse de 5%. Si on regarde une évolution sur deux jours, Autrement dit, si on subit une baisse de 10% (en log-rendement) puis une hausse de 10% (ou le contraire), on revient à la position initiale (et donc on a une variation nulle) car .
C’est probablement pour ça que les financiers utilisent des log-rendements: ils n’ont pas (en toute rigueur) la même interprétation que la variation relative, mais au moins on peut faire les calculs de manière très intuitive.

# Encore une année record…

Comme l’a noté la blogosphère ces dernières heures, il semble que des records aient (encore) été battus en terme de catastrophes naturelles: “catastrophes naturelles : une année 2010 exceptionnelle” (ici) ou “2010, année record des catastrophes naturelles” (). Maintenant on peut remonter un peu dans le temps aussi: “catastrophes naturelles: 2008 bat des records” (ici) ou “2008 a battu des records de catastrophes naturelles” (), “en 2007, les catastrophes naturelles ont atteint un nombre record” (ici), “catastrophes naturelles : 2005 bat des records en matière d’intensité et de coûts” (ici) ou “les catastrophes ont entraîné une facture record en 2005” (), “déjà record, le bilan 2004 des catastrophes naturelles s’alourdit” (ici), “année 2000 : record de catastrophes naturelles” (ici) et “catastrophes naturelles : encore une année record” (en 2000, ) etc, etc. Bon, en tant que statisticien, les rares résultats que je connaisse sur les records sont des théorèmes en log (i.e. en variation lente). Si est une suite de variables i.i.d. et que l’on définie les durées entre records,  , avec . Alfred Rényi a montré que alors que Marcel Neuts a montré que En particulier, on note que Autrement dit, plus le temps avance, plus les records doivent devenir rares.
Par exemple, regardons la simulation suivante. Considérons un processus Poissonnien homogène, , avec les dates des catastrophes, et leur coût. On note le nombre annuel de catastrophes (éventuellement, on ne considère que les catastrophes dont le cout dépasse un seuil prédéfini, comme cela semble généralement être le cas). On peut définir le processus de records, qui ressemble à ça Bon, maintenant, si on a des records tous les ans (ou presque) c’est probablement que l’hypothèse i.i.d. n’était pas valide.
Supposons par exemple qu’il y ait de l’inflation, i.e. les coût sont . Dans ce cas, le processus de nombre d’événements annuels est et celui des records Si on répète l’expérience un grand nombre de fois, le nombre moyen d’événements par année me cesse de croître (en rouge ci-dessous) et le processus de records moyen (avec le cas i.i.d. en bleu), Autrement dit, si l’inflation est trop importante, il est possible d’avoir régulièrement de nouveaux records.
Une autre possibilité est de supposer qu’il n’y a pas d’inflation, mais que le nombre de sinistres ne cesse de croître (ici un processus de Poisson à intensité linéaire), c’est à dire une fréquence plus importante des événements climatiques extrêmes Le nombre annuel de catastrophes est alors soit, si l’on considère le processus de records Là encore, si l’on répète l’expérience un grand nombre de fois, le nombre moyen de catastrophes ne cesse de croître alors que le processus de records a la moyenne suivante Bref, il est possible d’avoir effectivement des records presque tous les ans. Cela signifie probablement que les sinistres coûtent de plus en plus cher (hausse de la valeur à risque), mais aussi probablement que les catastrophes sont de plus en plus fréquentes. On retrouve finalement ce que je disais il y a quelques mois maintenant (ici).