# Non-observable vs. observable heterogeneity factor

This morning, in the ACT2040 class (on non-life insurance), we’ve discussed the difference between observable and non-observable heterogeneity in ratemaking (from an economic perspective). To illustrate that point (we will spend more time, later on, discussing observable and non-observable risk factors), we looked at the following simple example. Let $X$ denote the height of a person. Consider the following dataset

+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")

There is a small typo in the dataset, so let us make manual changes here

> Davis[12,c(2,3)]=Davis[12,c(3,2)]

Here, the variable of interest is the height of a given person,

> X=Davis$height If we look at the histogram, we have > hist(X,col="light green", border="white",proba=TRUE,xlab="",main="") Can we assume that we have a Gaussian distribution ? $X\sim\mathcal{N}(\theta,\sigma^2)$Maybe not… Here, if we fit a Gaussian distribution, plot it, and add a kernel based estimator, we get > (param <- fitdistr(X,"normal")$estimate)
> f1 <- function(x) dnorm(x,param[1],param[2])
> x=seq(100,210,by=.2)
> lines(x,f1(x),lty=2,col="red")
> lines(density(X))

If you look at that black line, you might think of a mixture, i.e. something like

$X\sim p_1\cdot\mathcal{N}(\theta_1,\sigma_1^2)+p_2\cdot\mathcal{N}(\theta_2,\sigma_2^2)$

(using standard mixture notations). Mixture are obtained when we have a non-observable heterogeneity factor: with probability $p_1$, we have a random variable $\mathcal{N}(\mu_1,\sigma_1^2)$ (call it type [1]), and with probability $p_2$, a random variable $\mathcal{N}(\mu_2,\sigma_2^2)$ (call it type [2]). So far, nothing new. And we can fit such a mixture distribution, using e.g.

> library(mixtools)
> mix <- normalmixEM(X)
number of iterations= 335
> (param12 <- c(mix$lambda[1],mix$mu,mix$sigma)) [1] 0.4002202 178.4997298 165.2703616 6.3561363 5.9460023 If we plot that mixture of two Gaussian distributions, we get > f2 <- function(x){ param12[1]*dnorm(x,param12[2],param12[4]) + (1-param12[1])*dnorm(x,param12[3],param12[5]) } > lines(x,f2(x),lwd=2, col="red") lines(density(X)) Not bad. Actually, we can try to maximize the likelihood with our own codes, > logdf <- function(x,parameter){ + p <- parameter[1] + m1 <- parameter[2] + s1 <- parameter[4] + m2 <- parameter[3] + s2 <- parameter[5] + return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2))) + } > logL <- function(parameter) -sum(logdf(X,parameter)) > Amat <- matrix(c(1,-1,0,0,0,0, + 0,0,0,0,1,0,0,0,0,0,0,0,0,1), 4, 5) > bvec <- c(0,-1,0,0) > constrOptim(c(.5,160,180,10,10), logL, NULL, ui = Amat, ci = bvec)$par

[1]   0.5996263 165.2690084 178.4991624   5.9447675   6.3564746

Here, we include some constraints, to insurance that the probability belongs to the unit interval, and that the variance parameters remain positive. Note that we have something close to the previous output.

Let us try something a little bit more complex now. What if we assume that the underlying distributions have the same variance, namely

$X\sim p_1\cdot\mathcal{N}(\theta_1,\sigma^2)+p_2\cdot\mathcal{N}(\theta_2,\sigma^2)$

In that case, we have to use the previous code, and make small changes,

> logdf <- function(x,parameter){
+ p <- parameter[1]
+ m1 <- parameter[2]
+ s1 <- parameter[4]
+ m2 <- parameter[3]
+ s2 <- parameter[4]
+ return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2)))
+ }
> logL <- function(parameter) -sum(logdf(X,parameter))
> Amat <- matrix(c(1,-1,0,0,0,0,0,0,0,0,0,1), 3, 4)
> bvec <- c(0,-1,0)
> (param12c= constrOptim(c(.5,160,180,10), logL, NULL, ui = Amat, ci = bvec)$par) [1] 0.6319105 165.6142824 179.0623954 6.1072614 This is what we can do if we cannot observe the heterogeneity factor. But wait… we actually have some information in the dataset. For instance, we have the sex of the person. Now, if we look at histograms of height per sex, and kernel based density estimator of the height, per sex, we have So, it looks like the height for male, and the height for female are different. Maybe we can use that variable, that was actually observed, to explain the heterogeneity in our sample. Formally, here, the idea is to consider a mixture, with an observable heterogeneity factor: the sex, $X\sim p_H\cdot\mathcal{N}(\theta_H,\sigma_H^2)+p_F\cdot\mathcal{N}(\theta_F,\sigma_F^2)$ We now have interpretation of what we used to call class [1] and [2] previously: male and female. And here, estimating parameters is quite simple, > (pM <- mean(sex=="M")) [1] 0.44 > (paramF <- fitdistr(X[sex=="F"],"normal")$estimate)
mean         sd
164.714286   5.633808
>  (paramM <- fitdistr(X[sex=="M"],"normal")$estimate) mean sd 178.011364 6.404001 And if we plot the density, we have > f4 <- function(x) pM*dnorm(x,paramM[1],paramM[2])+(1-pM)*dnorm(x,paramF[1],paramF[2]) > lines(x,f4(x),lwd=3,col="blue") What if, once again, we assume identical variance? Namely, the model becomes $X\sim p_H\cdot\mathcal{N}(\theta_H,\sigma^2)+p_F\cdot\mathcal{N}(\theta_F,\sigma^2)$Then a natural idea to derive an estimator for the variance, based on previous computations, is to use $\sigma^2=\frac{1}{n-2}\left(\sum_{i:H} [X_i-\overline{X}_H]^2+\sum_{i:F} [X_i-\overline{X}_F]^2\right)$ The code is here > s=sqrt((sum((height[sex=="M"]-paramM[1])^2)+sum((height[sex=="F"]-paramF[1])^2))/(nrow(Davis)-2)) > s [1] 6.015068 and again, it is possible to plot the associated density, > f5 <- function(x) pM*dnorm(x,paramM[1],s)+(1-pM)*dnorm(x,paramF[1],s) > lines(x,f5(x),lwd=3,col="blue") Now, if we think a little about what we’ve just done, it is simply a linear regression on a factor, the sex of the person, $X=\mu_H\cdot\boldsymbol{1}(H)+\mu_F\cdot\boldsymbol{1}(F)+\varepsilon$ where $\varepsilon\sim\mathcal{N}(0,\sigma^2)$. And indeed, if we run the code to estimate this linear model, > summary(lm(height~sex,data=Davis)) Call: lm(formula = height ~ sex, data = Davis) Residuals: Min 1Q Median 3Q Max -16.7143 -3.7143 -0.0114 4.2857 18.9886 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 164.7143 0.5684 289.80 <2e-16 *** sexM 13.2971 0.8569 15.52 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 6.015 on 198 degrees of freedom Multiple R-squared: 0.5488, Adjusted R-squared: 0.5465 F-statistic: 240.8 on 1 and 198 DF, p-value: < 2.2e-16 we get the same estimators for the means and the variance as the ones obtained previously. So, as mentioned this morning in class, if you have a non-observable heterogeneity factor, we can use a mixture model to fit a distribution, but if you can get a proxy of that factor, that is observable, then you can run a regression. But most of the time, that observable variable is just a proxy of a non-observable one… # Exposure as a possible explanatory variable Iin insurance pricing, the exposure is usually used as an offset variable to model claims frequency. As explained many times on this blog (e.g. here), and in my notes, if we have to identical drivers, but one with an exposure of 6 months, and the other one of one year, it should be natural to assume that, on average, the second driver will have two times more accidents. This is the motivation to use a standard (homogeneous) Poisson process to model claim frequency. One can also see here legal issue, since, in case of a (partial) reinbursement of a premium, it would be done prorata temporis. The risk is proportional to the exposure. Thus, if $Y_i$ denote the number of claims of insured $i$, with characteristics $\boldsymbol{X}_{i}=(X_{i,1},\cdots,X_{i,k})$ and exposure $E_i$, with a Poisson regression, we would write $Y_i\sim\mathcal{P}(E_i\cdot \exp(\boldsymbol{X}_i'\boldsymbol{\beta}))$ or equivalently $Y_i\sim\mathcal{P}(\exp(\log(E_i)+\boldsymbol{X}_i'\boldsymbol{\beta}))$ From this expression, the logarithm of the exposure is an explanatory variable, but there should be no coefficient (the coefficient here is taken to be one). Can’t we use the exposure as an explanatory variable ? Will we get a unit parameter ? Of course, in the context of ratemaking, it is probably not a relevant question, since actuaries are required to predict annual claim frequency (since insurance contract are supposed to provide a one year coverage). But it might be interesting to get a better understanding of why people might be leaving our portfolio (i.e. are cancelling their insurance policy before term, or not renew someday). To be more specific and get a better understanding, consider the following model: consider a Poisson process to model claims arrival, and people dedicated to their insurance company (they never leave). Let us generate scenarios over twenty years > n=983 > D1=as.Date("01/01/1993",'%d/%m/%Y') > D2=as.Date("31/12/2013",'%d/%m/%Y') > L=D1+0:(D2-D1) > set.seed(1) > arrival=sample(L,size=n,replace=TRUE) > exposure=N=rep(NA,n) > departure=rep(D2,n) > set.seed(2) > for(i in 1:n){ + expo=D2-arrival[i] + w=0 + while(max(w)<expo) w=c(w,max(w)+1+trunc(rexp(1,1/1000))) + exposure[i]=departure[i]-arrival[i] + N[i]=max(0,length(w)-2)} > df=data.frame(N=N,E=exposure/365) Here the expected time between claims is considered to be 1000 days. The (annual) intensity of the Poisson process is here > 365/1000 [1] 0.365 so if we run a Poisson regression on the logarithm of the exposure (please feel free to had other covariates if you want, the example here is just to see what could happen when exposure is considered as a standard covariate), we should get a parameter close to > log(365/1000) [1] -1.007858 Here, the regression on a constant, with the offset variable is > reg=glm(N~1+offset(log(E)),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ 1 + offset(log(E)), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -3.4145 -0.4673 0.2367 0.8770 3.6828 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.04233 0.02532 -41.17 <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: 1116.9 on 982 degrees of freedom Residual deviance: 1116.9 on 982 degrees of freedom AIC: 3282.9 Number of Fisher Scoring iterations: 5 which is consistent with what we just said. If we run the regression with the logarithm of the exposure as a possible explanatory variable, we would expect to have a coefficient close to 1. And indeed… > reg=glm(N~log(E),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ log(E), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -3.0810 -0.8373 -0.1493 0.5676 3.9001 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.03350 0.08546 -12.09 <2e-16 *** log(E) 1.00920 0.03292 30.66 <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: 2553.6 on 982 degrees of freedom Residual deviance: 1064.2 on 981 degrees of freedom AIC: 3762.7 Number of Fisher Scoring iterations: 5 If we keep the offset, and add the variable, we can see that it become useless (which is a test of a unit parameter, somehow) > reg=glm(N~log(E)+offset(log(E)),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ log(E) + offset(log(E)), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -3.0810 -0.8373 -0.1493 0.5676 3.9001 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.033503 0.085460 -12.093 <2e-16 *** log(E) 0.009201 0.032920 0.279 0.78 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1064.3 on 982 degrees of freedom Residual deviance: 1064.2 on 981 degrees of freedom AIC: 3762.7 Number of Fisher Scoring iterations: 5 Here, we do have pure Poisson processes, so exposure is crucial, since the parameter of the Poisson distribution is proportional to the exposure. But we cannot learn anything else from the exposure. Consider some real data. > head(baseFREQ) nocontrat exposition zone puissance agevehicule 1 27 0.87 C 7 0 2 115 0.72 D 5 0 3 121 0.05 C 6 0 4 142 0.90 C 10 10 5 155 0.12 C 7 0 6 186 0.83 C 5 0 ageconducteur bonus marque carburant densite region nbre 1 56 50 12 D 93 13 0 2 45 50 12 E 54 13 0 3 37 55 12 D 11 13 0 4 42 50 12 D 93 13 0 5 59 50 12 E 73 13 0 6 75 50 12 E 42 13 0 What do we get if we consider a Poisson regression on the logarithm of the exposure ? > reg=glm(nbre~log(exposition),data=baseFREQ,family=poisson) > summary(reg) Call: glm(formula = nbre ~ log(exposition), family = poisson, data = baseFREQ) 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(exposition) 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 If we add the exposure to the offset, what’s happening ? (let us use a nonparametric transformation, so visualize what’s going on) > library(gam) > reg=gam(nbre~offset(log(exposition))+s(exposition),data=baseFREQ,family=poisson) > plot(reg,se=TRUE) There is a clear and significant effect. The more insured stay, the less likely they get a claim. Actually, it can be observed without running a regression. > i1=which(baseFREQ$nbre>0)
> i0=which(baseFREQ$nbre==0) > h1=hist(baseFREQ$exposition[i1],probability=TRUE)
> h0=hist(baseFREQ$exposition[i0],probability=TRUE) > plot(h1$mids,h1$density,type='s',lwd=2,col="red") > lines(h0$mids,h0$density,type='s',col='blue',lwd=2) In blue, we have the density of the exposure for those who did not have claims, and in red, the density of those who did have one claim (or more) So here, we cannot assume a unit value for the parameter. What does that mean ? Can we reproduce such a behavior ? In order to get a better understandung, consider two possible behaviors for the insured. The first one will be the following : if the company does not offer substantial discounts after no several years with no claims, the insured might leave the company. For instance, if the insured has no claim during 5 years, then after 5 years, he will leave the company (to get a better price somewhere else, say). The code will be > for(i in 1:n){ + expo=D2-arrival[i] + w=c(0,0) + while((max(w)<expo) & (max(diff(w))<1500)) w=c(w,max(w)+trunc(rexp(1,1/1000))) + if(max(diff(w))>1500) departure[i]=arrival[i]+max(w[-length(w)])+1500 + exposure[i]=departure[i]-arrival[i] + N[i]=max(0,length(w)-3)} > df=data.frame(N=N,E=exposure/365) Here, I consider 1500 days, instead of 5 years,, but it is the same idea. So, what do we have here ? > reg=glm(N~log(E),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ log(E), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -1.5684 -0.9668 -0.2321 0.4244 3.6265 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.50844 0.10286 -24.39 <2e-16 *** log(E) 1.65738 0.04494 36.88 <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: 2567.31 on 982 degrees of freedom Residual deviance: 885.71 on 981 degrees of freedom AIC: 2897.9 Here, the coefficient is (significantly) larger than 1. More precisely, > reg=glm(N~log(E)+offset(log(E)),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ log(E) + offset(log(E)), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -1.5684 -0.9668 -0.2321 0.4244 3.6265 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.50844 0.10286 -24.39 <2e-16 *** log(E) 0.65738 0.04494 14.63 <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: 1114.24 on 982 degrees of freedom Residual deviance: 885.71 on 981 degrees of freedom AIC: 2897.9 There is clearly a bias here : people staying long are more like likely to have an accident. Which is consistent with our story, since clients with low risks left. The second behavior will be the following : sometimes, the insured are not satisfied with the way claims are handled, and they might leave after the first claim. Consider the case where, after one claim, it is likely (e.g. with probability 50%) that the insured leaves the company. Instead of assuming that the insured did not like claims management, consider the case were the car is so damaged that he cannot drive it anymore. So it will be useless to pay an insurance premium. The code here will be > for(i in 1:n){ + expo=D2-arrival[i] + w=0 + stay=TRUE + while((max(w)<expo) & (stay==TRUE)) { w=c(w,max(w)+trunc(rexp(1,1/1000))) + stay=sample(c(TRUE,FALSE),prob=c(.5,.5),size=1)} + N[i]=length(w)-2 + if(stay==FALSE) {departure[i]=arrival[i]+max(w) + N[i]=length(w)-1} + exposure[i]=departure[i]-arrival[i]} > df=data.frame(N=N,E=exposure/365) Here, after each claim, the insured toss a coin to see if he cancels the contract, or not. > reg=glm(N~log(E),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ log(E), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -2.28402 -0.47763 -0.08215 0.33819 2.37628 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.09920 0.04251 2.334 0.0196 * log(E) 0.30640 0.02511 12.203 <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: 666.92 on 982 degrees of freedom Residual deviance: 498.29 on 981 degrees of freedom AIC: 2666.3 This time, the parameter is (again significantly) smaller than one. > reg=glm(N~log(E)+offset(log(E)),data=df,family=poisson) > summary(reg) Call: glm(formula = N ~ log(E) + offset(log(E)), family = poisson, data = df) Deviance Residuals: Min 1Q Median 3Q Max -2.28402 -0.47763 -0.08215 0.33819 2.37628 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.09920 0.04251 2.334 0.0196 * log(E) -0.69360 0.02511 -27.625 <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: 1116.87 on 982 degrees of freedom Residual deviance: 498.29 on 981 degrees of freedom AIC: 2666.3 The story is now rather different, since those who stay long should not have encountered a lot of opportunities to leave. So clearly, they did not have much claims. If someone has a long exposure, the negative sign in the output above means that he should not have much claims, on average. As we can see, those models produce rather difference outputs. Note that it is possible much more interpretations. For instance, depending on the way data were extracted, • all policies observed, over those twenty years, • all policies in force at some specific date, until now • all policies in force at some specific date, until one year after • all policies in force now So far, we have been using the first method, but the other ones will yield different interpretations, e.g. because of survivor bias. But that’s another story… And one can read Boucher and Denuit (2008) to go further. # Poisson regression on non-integers In the course on claims reserving techniques, I did mention the use of Poisson regression, even if incremental payments were not integers. For instance, we did consider incremental triangles > source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R") > INC=PAID > INC[,2:6]=PAID[,2:6]-PAID[,1:5] > INC [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 1163 39 17 7 21 [2,] 3367 1292 37 24 10 NA [3,] 3871 1474 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA On those payments, it is natural to use a Poisson regression, to predict future payments > Y=as.vector(INC) > D=rep(1:6,each=6) > A=rep(2001:2006,6) > base=data.frame(Y,D,A) > reg=glm(Y~as.factor(D)+as.factor(A),data=base,family=poisson(link="log")) > Yp=predict(reg,type="response",newdata=base) > matrix(Yp,6,6) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3155.6 1202.1 49.8 19.1 8.2 21.0 [2,] 3365.6 1282.0 53.1 20.4 8.7 22.3 [3,] 3863.7 1471.8 60.9 23.4 10.0 25.7 [4,] 4310.0 1641.8 68.0 26.1 11.2 28.6 [5,] 4919.8 1874.1 77.6 29.8 12.8 32.7 [6,] 5217.0 1987.3 82.3 31.6 13.5 34.7 and the total amount of reserves would be > sum(Yp[is.na(Y)==TRUE]) [1] 2426.985 Here, payments were in ‘000 euros. What if they were in ‘000’000 euros ? > a=1000 > INC/a [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3.209 1.163 0.039 0.017 0.007 0.021 [2,] 3.367 1.292 0.037 0.024 0.010 NA [3,] 3.871 1.474 0.053 0.022 NA NA [4,] 4.239 1.678 0.103 NA NA NA [5,] 4.929 1.865 NA NA NA NA [6,] 5.217 NA NA NA NA NA We can still run a regression here > reg=glm((Y/a)~as.factor(D)+as.factor(A),data=base,family=poisson(link="log")) > Yp=predict(reg,type="response",newdata=base) > sum(Yp[is.na(Y)==TRUE])*a [1] 2426.985 and the prediction is exactly the same. Actually, it is possible to change currency, and multiply by any kind of constant, the Poisson regression will return always the same prediction, if we use a log link function, > homogeneity=function(a=1){ + reg=glm((Y/a)~as.factor(D)+as.factor(A), data=base,family=poisson(link="log")) + Yp=predict(reg,type="response",newdata=base) + return(sum(Yp[is.na(Y)==TRUE])*a) + } > Vectorize(homogeneity)(10^(seq(-3,5))) [1] 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 The trick here come from the fact that we do like the Poisson interpretation. But GLMs simply mean that we do want to solve a first order condition. It is possible to solve explicitly the first order condition, which was obtained without any condition such that values should be integers. To run a simple code, the intercept should be related to the last value of the matrix, not the first one. > base$D=relevel(as.factor(base$D),"6") > base$A=relevel(as.factor(base$A),"2006") > reg=glm(Y~as.factor(D)+as.factor(A), data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = Y ~ as.factor(D) + as.factor(A), family = poisson(link = "log"), data = base) Deviance Residuals: Min 1Q Median 3Q Max -2.3426 -0.4996 0.0000 0.2770 3.9355 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.54723 0.21921 16.182 < 2e-16 *** as.factor(D)1 5.01244 0.21877 22.912 < 2e-16 *** as.factor(D)2 4.04731 0.21896 18.484 < 2e-16 *** as.factor(D)3 0.86391 0.22827 3.785 0.000154 *** as.factor(D)4 -0.09254 0.25229 -0.367 0.713754 as.factor(D)5 -0.93717 0.32643 -2.871 0.004092 ** as.factor(A)2001 -0.50271 0.02079 -24.179 < 2e-16 *** as.factor(A)2002 -0.43831 0.02045 -21.433 < 2e-16 *** as.factor(A)2003 -0.30029 0.01978 -15.184 < 2e-16 *** as.factor(A)2004 -0.19096 0.01930 -9.895 < 2e-16 *** as.factor(A)2005 -0.05864 0.01879 -3.121 0.001799 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 46695.269 on 20 degrees of freedom Residual deviance: 30.214 on 10 degrees of freedom (15 observations deleted due to missingness) AIC: 209.52 The first idea is to run a gradient descent, as follows (the starting point will be coefficients from a linear regression on the log of the observations), > YNA <- Y > XNA=matrix(0,length(Y),1+5+5) > XNA[,1]=rep(1,length(Y)) > for(k in 1:5) XNA[(k-1)*6+1:6,k+1]=k > u=(1:(length(Y))%%6); u[u==0]=6 > for(k in 1:5) XNA[u==k,k+6]=k > YnoNA=YNA[is.na(YNA)==FALSE] > XnoNA=XNA[is.na(YNA)==FALSE,] > beta=lm(log(YnoNA)~0+XnoNA)$coefficients
>     for(s in 1:50){
+     Ypred=exp(XnoNA%*%beta)
+     omega=matrix(0,nrow(XnoNA),nrow(XnoNA));diag(omega)=exp(XnoNA%*%beta)
+     hessienne=-t(XnoNA)%*%omega%*%XnoNA
> beta
[,1]
[1,]  3.54723486
[2,]  5.01244294
[3,]  2.02365553
[4,]  0.28796945
[5,] -0.02313601
[6,] -0.18743467
[7,] -0.50271242
[8,] -0.21915742
[9,] -0.10009587
[10,] -0.04774056
[11,] -0.01172840

We are not too far away from the values given by R. Actually, it is just fine if we focus on the predictions

> matrix(exp(XNA%*%beta),6,6))
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.6 1202.1 49.8 19.1  8.2 21.0
[2,] 3365.6 1282.0 53.1 20.4  8.7 22.3
[3,] 3863.7 1471.8 60.9 23.4 10.0 25.7
[4,] 4310.0 1641.8 68.0 26.1 11.2 28.6
[5,] 4919.8 1874.1 77.6 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.3 31.6 13.5 34.7

which are exactly the one obtained above. And here, we clearly see that there is no assumption such as “explained variate should be an integer“. It is also possible to remember that the first order condition is the same as the one we had with a weighted least square model. The problem is that the weights are function of the prediction. But using an iterative algorithm, we should converge,

> beta=lm(log(YnoNA)~0+XnoNA)$coefficients > for(i in 1:50){ + Ypred=exp(XnoNA%*%beta) + z=XnoNA%*%beta+(YnoNA-Ypred)/Ypred + REG=lm(z~0+XnoNA,weights=Ypred) + beta=REG$coefficients
+ }
>
> beta
XnoNA1      XnoNA2      XnoNA3      XnoNA4      XnoNA5      XnoNA6
3.54723486  5.01244294  2.02365553  0.28796945 -0.02313601 -0.18743467
XnoNA7      XnoNA8      XnoNA9     XnoNA10     XnoNA11
-0.50271242 -0.21915742 -0.10009587 -0.04774056 -0.01172840

which are the same values as the one we got previously. Here again, the prediction is the same as the one we got from this so-called Poisson regression,

> matrix(exp(XNA%*%beta),6,6)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.6 1202.1 49.8 19.1  8.2 20.9
[2,] 3365.6 1282.0 53.1 20.4  8.7 22.3
[3,] 3863.7 1471.8 60.9 23.4 10.0 25.7
[4,] 4310.0 1641.8 68.0 26.1 11.2 28.6
[5,] 4919.8 1874.1 77.6 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.3 31.6 13.5 34.7

Again, it works just fine because GLMs are mainly conditions on the first two moments, and numerical computations are based on the first order condition, which has less constraints than the interpretation in terms of a Poisson model.

# Provisionnement et tarification, examen final

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.

# Overdispersed Poisson et bootstrap

Pour le dernier cours sur les méthodes de provisionnement, on s’est arrête aux méthodes par simulation. Reprenons là où on en était resté au dernier billet où on avait vu qu’en faisant une régression de Poisson sur les incréments, on obtenait exactement le même montant que la méthode Chain Ladder,

> Y
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   53   22   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

> y=as.vector(as.matrix(Y))
> base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
> reg2=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson,
data = base)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        8.05697    0.01551 519.426  < 2e-16 ***
as.factor(ai)2001  0.06440    0.02090   3.081  0.00206 **
as.factor(ai)2002  0.20242    0.02025   9.995  < 2e-16 ***
as.factor(ai)2003  0.31175    0.01980  15.744  < 2e-16 ***
as.factor(ai)2004  0.44407    0.01933  22.971  < 2e-16 ***
as.factor(ai)2005  0.50271    0.02079  24.179  < 2e-16 ***
as.factor(bj)1    -0.96513    0.01359 -70.994  < 2e-16 ***
as.factor(bj)2    -4.14853    0.06613 -62.729  < 2e-16 ***
as.factor(bj)3    -5.10499    0.12632 -40.413  < 2e-16 ***
as.factor(bj)4    -5.94962    0.24279 -24.505  < 2e-16 ***
as.factor(bj)5    -5.01244    0.21877 -22.912  < 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: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: 209.52

Number of Fisher Scoring iterations: 4

> base$py2=predict(reg2,newdata=base,type="response") > round(matrix(base$py2,n,n),1)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.7 1202.1 49.8 19.1  8.2 21.0
[2,] 3365.6 1282.1 53.1 20.4  8.8 22.4
[3,] 3863.7 1471.8 61.0 23.4 10.1 25.7
[4,] 4310.1 1641.9 68.0 26.1 11.2 28.7
[5,] 4919.9 1874.1 77.7 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.4 31.6 13.6 34.7

> sum(base$py2[is.na(base$y)])
[1] 2426.985

Le plus intéressant est peut être de noter que la loi de Poisson présente probablement trop peu de variance…

> reg2b=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = quasipoisson,
data = base)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)        8.05697    0.02769 290.995  < 2e-16 ***
as.factor(ai)2001  0.06440    0.03731   1.726 0.115054
as.factor(ai)2002  0.20242    0.03615   5.599 0.000228 ***
as.factor(ai)2003  0.31175    0.03535   8.820 4.96e-06 ***
as.factor(ai)2004  0.44407    0.03451  12.869 1.51e-07 ***
as.factor(ai)2005  0.50271    0.03711  13.546 9.28e-08 ***
as.factor(bj)1    -0.96513    0.02427 -39.772 2.41e-12 ***
as.factor(bj)2    -4.14853    0.11805 -35.142 8.26e-12 ***
as.factor(bj)3    -5.10499    0.22548 -22.641 6.36e-10 ***
as.factor(bj)4    -5.94962    0.43338 -13.728 8.17e-08 ***
as.factor(bj)5    -5.01244    0.39050 -12.836 1.55e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 3.18623)

Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

Mais on en reparlera dans un instant. Ensuite, on avait commencé à regarder erreurs commises, sur la partie supérieure du triangle. Classiquement, par construction, les résidus de Pearson sont de la forme

$\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\text{Var}(Y_i)}}$

On avait vu dans le cours de tarification que la variance au dénominateur pouvait être remplacé par le prévision, puisque dans un modèle de Poisson, l’espérance et la variance sont identiques. Donc on considérait

$\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}$

> base$erreur=(base$y-base$py2)/sqrt(base$py2)
> round(matrix(base$erreur,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.9 -1.1 -1.5 -0.5 -0.4 0 [2,] 0.0 0.3 -2.2 0.8 0.4 NA [3,] 0.1 0.1 -1.0 -0.3 NA NA [4,] -1.1 0.9 4.2 NA NA NA [5,] 0.1 -0.2 NA NA NA NA [6,] 0.0 NA NA NA NA NA Le soucis est que si $\widehat{Y}_i$ est – asymptotiquement – un bon estimateur pour $\text{Var}(Y_i)$, ce n’est pas le cas à distance finie, car on a alors un estimateur biaisé pour la variance, et donc la variance des résidus n’a que peu de chances d’être de variance unitaire. Aussi, il convient de corriger l’estimateur de la variance, et on pose alors $\varepsilon_i=\sqrt{\frac{n}{n-k}}\cdot\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}$ qui sont alors les résidus de Pearson tel qu’on doit les utiliser. > E=base$erreur[is.na(base$y)==FALSE]*sqrt(21/(21-11)) > E [1] 1.374976e+00 3.485024e-02 1.693203e-01 -1.569329e+00 1.887862e-01 [6] -1.459787e-13 -1.634646e+00 4.018940e-01 8.216186e-02 1.292578e+00 [11] -3.058764e-01 -2.221573e+00 -3.207593e+00 -1.484151e+00 6.140566e+00 [16] -7.100321e-01 1.149049e+00 -4.307387e-01 -6.196386e-01 6.000048e-01 [21] -8.987734e-15 > boxplot(E,horizontal=TRUE) En rééchantillonnant dans ces résidus, on peut alors générer un pseudo triangle. Pour des raisons de simplicités, on va générer un peu rectangle, et se restreindre à la partie supérieure, > Eb=sample(E,size=36,replace=TRUE) > Yb=base$py2+Eb*sqrt(base$py2) > Ybna=Yb > Ybna[is.na(base$y)]=NA
> Tb=matrix(Ybna,n,n)
> round(matrix(Tb,n,n),1)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3115.8 1145.4 58.9 46.0  6.4 26.9
[2,] 3179.5 1323.2 54.5 21.3 12.2   NA
[3,] 4245.4 1448.1 61.0  7.9   NA   NA
[4,] 4312.4 1581.7 68.7   NA   NA   NA
[5,] 4948.1 1923.9   NA   NA   NA   NA
[6,] 4985.3     NA   NA   NA   NA   NA

Cette fois, on a un nouveau triangle ! on va alors pouvoir faire plusieurs choses,

1. compléter le triangle para la méthode Chain Ladder, c’est à dire calculer les montants moyens que l’on pense payer dans les années futures
2. générer des scénarios de paiements pour les années futurs, en générant des paiements suivant des lois de Poisson (centrées sur les montants moyens que l’on vient de calculer)
3. générer des scénarios de paiements avec des lois présentant plus de variance que la loi de Poisson. Idéalement, on voudrait simuler des lois qusi-Poisson, mais ce ne sont pas de vraies lois. Par contre, on peut se rappeler que dans ce cas, la loi Gamma devrait donner une bonne approximation.

Pour ce dernier point, on va utiliser le code suivant, pour générer des quasi lois,

> rqpois = function(n, lambda, phi, roundvalue = TRUE) {
+ b = phi
+ a = lambda/phi
+ r = rgamma(n, shape = a, scale = b)
+ if(roundvalue){r=round(r)}
+ return(r)
+ }

Je renvois aux diverses notes de cours pour plus de détails sur la justification, ou à un vieux billet. On va alors faire une petite fonction, qui va soit somme les paiements moyens futurs, soit sommer des générations de scénarios de paiements, à partir d’un triangle,

> CL=function(Tri){
+ y=as.vector(as.matrix(Tri))
+ base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
+ reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson)
+ py2=predict(reg,newdata=base,type="response")
+ pys=rpois(36,py2)
+ pysq=rqpois(36,py2,phi=3.18623)
+ return(list(
+ cl=sum(py2[is.na(base$y)]), + sc=sum(pys[is.na(base$y)]),
+ scq=sum(pysq[is.na(base$y)]))) + } Reste alors à générer des paquets de triangles. Toutefois, il est possible de générer des triangles avec des incréments négatifs. Pour faire simple, on mettra des valeurs nulles quand on a un paiement négatif. L’impact sur les quantiles sera alors (a priori) négligeable. > for(s in 1:1000){ + Eb=sample(E,size=36,replace=TRUE)*sqrt(21/(21-11)) + Yb=base$py2+Eb*sqrt(base$py2) + Yb=pmax(Yb,0) + scY=rpois(36,Yb) + Ybna=Yb + Ybna[is.na(base$y)]=NA
+ Tb=matrix(Ybna,6,6)
+ C=CL(Tb)
+ VCL[s]=C$cl + VR[s]=C$sc
+ VRq[s]=C$scq + } Si on regarde la distribution du best estimate, on obtient > hist(VCL,proba=TRUE,col="light blue",border="white",ylim=c(0,0.003)) > boxplot(VCL,horizontal=TRUE,at=.0023,boxwex = 0.0006,add=TRUE,col="light green") > D=density(VCL) > lines(D) > I=which(D$x<=quantile(VCL,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) > I=which(D$x>=quantile(VCL,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) Mais on peut aussi visualiser des scénarios basés sur des lois de Poisson (équidispersé) ou des scénarios de lois quasiPoisson (surdispersées), ci-dessous Dans ce dernier cas, on peut en déduire le quantile à 99% des paiements à venir. > quantile(VRq,.99) 99% 2855.01 Il faut donc augmenter le montant de provisions de l’ordre 15% pour s’assurer que la compagnie pourra satisfaire ses engagements dans 99% des cas, > quantile(VRq,.99)-2426.985 99% 428.025 # Reserving with negative increments in triangles A few months ago, I did published a post on negative values in triangles, and how to deal with them, when using a Poisson regression (the post was published in French). The idea was to use a translation technique: 1. Fit a model not on $Y_i$‘s but on $Y_i^{(k)}=Y_i+k$, for some $k\geq 0$, 2. Use that model to make predictions, and then translate those predictions, $\widehat{Y}_i^{(k)}-k$ This is what was done to get the following graph, where a Poisson regression was fitted. Black points are $Y_i$‘s while blue points are $\widehat{Y}_i^{(k)}$‘s, for some $k\geq 0$. We fit a model to get the blue prediction, and then translate it to get the red prediction (on the $Y_i$‘s). In this example, there were no negative values, but it is possible to use it get a better understanding on the impact of this technique. The prediction, here, is the red line. And clearly, the value of $k$ has an impact on the prediction (since we do not consider, here, a linear model: with a linear model, translating has not impact at all, except on the intercept). The alternative mentioned in the previous post was to use this technique on several $k$‘s, and them interpolate 1. For a given $k$, fit a model not on $Y_i$‘s but on $Y_i^{(k)}=Y_i+k$, use that model to make predictions, and then translate those predictions, $\widehat{Y}_i^{(k)}-k$. 2. Do it for several $k$‘s. 3. Use it to extrapolate when $k$ is $0$ (which is the case we are interested in). In the context of loss reserving, the idea is extremely simple. Consider a triangle with incremental payments > source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R") > Y=T=PAID > n=ncol(T) > Y[,2:n]=T[,2:n]-T[,1:(n-1)] > Y [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 1163 39 17 7 21 [2,] 3367 1292 37 24 10 NA [3,] 3871 1474 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA Now, we do not have negative values, here, but we can still see is translation techniques can be used. The benchmark is the Poisson regression, since we can run it : > y=as.vector(as.matrix(Y)) > base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n)) > reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson) Here, the amount is reserve is the sum of predicted values in the lower part of the triangle, > py=predict(reg,newdata=base,type="response") > sum(py[is.na(base$y)])
[1] 2426.985

which is exactly Chain Ladder’s estimate.

Now, let us use a translation technique to compute the amount of reserves. The code will be

> decal=function(k){
+ reg=glm(y+k~as.factor(ai)+as.factor(bj),data=base,family=poisson)
+ py=predict(reg,newdata=base,type="response")
+ return(sum(py[is.na(base$y)]-k)) For instance, if we translate of +5, we would get > decal(5) [1] 2454.713 while a translation of +10 would return > decal(10) [1] 2482.29 Clearly, translations do have an impact on the estimation. Here, just to check, if we do not translate, we do have Chain Ladder’s estimate, > decal(0) [1] 2426.985 The idea mentioned in the previous post was to try several translations, and then extrapolate, to get the value in 0. Here, translations will give the following estimates > K=10:20 > (V=Vectorize(decal)(K)) [1] 2482.290 2487.788 2493.279 2498.765 2504.245 2509.719 2515.187 2520.649 [9] 2526.106 2531.557 2537.001 We can plot those values, and run a regression > plot(K,V,xlim=c(0,20),ylim=c(2425,2540)) > abline(h=decal(0),col="red",lty=2) the dotted horizontal line is Chain Ladder. Now, let us extrapolate > b=data.frame(K=K,D=V) > rk=lm(D~K,data=b) > predict(rk,newdata=data.frame(K=0)) 1 2427.623 On has to admit that it is not that bad. But yesterday evening, Karim asked me why I did use a linear regression, for my extrapolation. And to be honest, I do not know. I mean, the only answer might be that points are almost on a straight line. So the first time I saw it, I was exited, and I ran a linear regression. Now, let us see if we can do better. Because here, we do use a translation of +10 or +20 (which might be rather small). What if we use much larger values ? (because we might have large negative incremental values). With the following code, we try, each time 11 consecutive values, the smallest one going from 0 to 50, > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=lm(D~K,data=b) + res[k]=predict(rk,newdata=data.frame(K=0)) + } > plot(hausse,res,type="l",col="red",ylim=c(2422,2440)) > abline(rk,col="blue") Here, we compute reserves when extrapolations were done after 11 translations, from $k$ to $k+10$. With different values of $k$. The case where $k$ is ten was the one mentioned above, > res[hausse==10] [1] 2427.623 Actually, it might also be possible to consider not 11 translations, but 26, from $k$ to $k+25$. Here, we get > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(25+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=lm(D~K,data=b) + res[k]=predict(rk,newdata=data.frame(K=0)) + } > lines(hausse,res,type="l",col="blue",lty=2) We now have the dotted line Here, it is getting worst. So let us keep here 11 translations. Perhaps, we can try something different. For instance a Poisson regression, with a log like (i.e. we consider an exponential extrapolation), > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=glm(D~K,data=b,family=poisson) + res[k]=predict(rk,newdata=data.frame(K=0),type="response") + } > lines(hausse,res,type="l",col="purple") The purple line will be a Poisson model, with a log link. Perhaps we can try another link function, like a quadratic one > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=glm(D~K,data=b,family=poisson(link= + power(lambda = 2))) + res[k]=predict(rk,newdata=data.frame(K=0),type="response") + } > lines(hausse,res,type="l",col="orange") That would be the orange line, Here, we need a link function between identity (the linear model, the blue line) and the quadratic one (the orange one), for instance a power function 3/2, > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=glm(D~K,data=b,family=poisson(link= + power(lambda = 1.5))) + res[k]=predict(rk,newdata=data.frame(K=0),type="response") + } > lines(hausse,res,type="l",col="green") Here, it looks like we can use that model for any kind of translation, from +10 till +50, even +100 ! But I do not have any intuition about the use of this power function… # Bootstrap et régression Lors du dernier cours, on a évoqué l’utilisation du bootstrap pour obtenir des intervalles de confiance sur des prévisions. Je mets en ligne les codes tapés en cours (très sommairement commentés, je peux renvoyer vers des vieux billets du cours ACT6420 pour des compléments). On va travailler sur ma base préférée pour évoquer la régression linéaire (avant de parler triangles de provisionnement, revenons cinq minutes sur des choses simples). > plot(cars) > reg=lm(dist~speed,data=cars) > abline(reg,col="red") > n=nrow(cars) > x=21 > points(x,predict(reg,newdata= data.frame(speed=x)),pch=19,col="red") On cherche ici à faire une prédiction en un point. Comme rappelé en cours (mais aussi dans le cours de modèles de prévision), quand on veut donner un intervalle de confiance pour la prévision, il convient de distinger l’intervalle de confiance pour le prédicteur (qui va dépendre de l’erreur d’estimation des paramétres) et l’intervalle de confiance pour une potentielle valeur (on peut parler de génération de scénarios, qui va dépendre en plus de l’erreur de modèle, c’est à dire de la dispersion des résidus). Commençons par l’intervalle de confiance sur la prédiction, sur le best estimate comme on dit en provisionnement > Yx=rep(NA,500) > B=matrix(NA,500,2) > for(s in 1:500){ + indice=sample(1:n,size=n, + replace=TRUE) + base=cars[indice,] + #points(base,pch=3) + reg=lm(dist~speed,data=base) + abline(reg,col="light blue") + points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue") + Yx[s]=predict(reg,newdata=data.frame(speed=x)) + B[s,]=coefficients(reg) + } Les valeurs bleues sont ici des prévisions possibles, obtenues en rééchantillonnant dans notre base d’observations. Pour rappel, l’intervalle de confiance (à 90%), sous hypothèse de normalité des résidus (et donc des estimateurs de la pente et de la constante de la droite de régression) s’obtient de la manière suivante > reg=lm(dist~speed,data=cars) > U=predict(reg,interval ="confidence", + level=.9,newdata= + data.frame(speed=0:30)) > lines(0:30,U[,2],col="red",lwd=2) > lines(0:30,U[,3],col="red",lwd=2) On peut comparer ici la distribution des valeurs obtenues sur nos 500 jeux de données générées, et comparer les quantiles empiriques, avec les quantiles sous hypothèse de normalité, > hist(Yx,proba=TRUE,col="light blue",border="white") > boxplot(Yx,horizontal=TRUE,at=.07,boxwex = 0.02,add=TRUE,col="light green") > abline(v=U[x+1,2:3],col="red",lwd=2) > D=density(Yx) > lines(D) > I=which(D$x<=quantile(Yx,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) > I=which(D$x>=quantile(Yx,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) On peut noter que les ordres de grandeur sont comparables. > reg=lm(dist~speed,data=cars) > quantile(Yx,c(.05,.95)) 5% 95% 58.63689 70.31281 > predict(reg,interval ="confidence", + level=.9,newdata=data.frame(speed=x)) fit lwr upr 1 65.00149 59.65934 70.34364 Regardons maintenant l’autre type d’intervalle de confiance, sur la valeur possible de la variable d’intérêt. Cette fois, en plus de tirer des nouveaux échantillons et calculer des prédictions, on va en plus rajouter un bruit à chaque tirage, qui permettra d’obtenir une valeur possible. > Yx=rep(NA,500) > for(s in 1:500){ + indice=sample(1:n,size=n, + replace=TRUE) + base=cars[indice,] + #points(base,pch=3) + reg=lm(dist~speed,data=base) + erreur=residuals(reg) + #abline(reg,lty=2) + E=sample(erreur,size=1) + Yx[s]=predict(reg,newdata=data.frame(speed=x))+E + points(x,Yx[s],pch=19,col="red") + } Là encore, on peut comparer (graphiquement pour commencer) les valeurs obtenues par rééchantillonnage, et celle obtenues sous hypothèse de normalité, > hist(Yx,proba=TRUE,col="light blue",border="white") > boxplot(Yx,horizontal=TRUE,at=.025,boxwex = 0.005,add=TRUE,col="light green") > abline(v=U[2:3],col="red",lwd=2) > D=density(Yx) > lines(D) > I=which(D$x<=quantile(Yx,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) > I=which(D$x>=quantile(Yx,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) Ce qui donne, numériquement, les comparaisons suivantes > quantile(Yx,c(.05,.95)) 5% 95% 44.43468 96.01357 > (U=predict(reg,interval ="prediction",level=.9,newdata=data.frame(speed=x))) fit lwr upr 1 67.63136 45.16967 90.09305 On observe cette fois une légère asymméytrie vers la droite. Manifestement, on ne peut pas supposer les résidus Gaussien, car il y a plus de grandes valeurs positives, que négatives. Ce qui fait du sens compte tenu de la nature des données (une distance ne peut être négative). On avait ensuite commencé à discuter de l’utilisation des modèles de régression en provisionnement. Afin d’avoir des données présentant de l’indépendance, on avait rappelé qu’il fallait travailler avec les incréments de paiments, et non pas les paiements cumulés. > T [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 4372 4411 4428 4435 4456 [2,] 3367 4659 4696 4720 4730 NA [3,] 3871 5345 5398 5420 NA NA [4,] 4239 5917 6020 NA NA NA [5,] 4929 6794 NA NA NA NA [6,] 5217 NA NA NA NA NA > n=ncol(T) > Y=T > Y[,2:n]=T[,2:n]- + T[,1:(n-1)] > Y [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 1163 39 17 7 21 [2,] 3367 1292 37 24 10 NA [3,] 3871 1474 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA On peut alors constituer une base de données, avec comme variables explicatives la ligne et la colonne. > y=as.vector(as.matrix(Y)) > base=data.frame( + y, + ai=rep(2000:2005,n), + bj=rep(0:(n-1),each=n)) > > head(base,12) y ai bj 1 3209 2000 0 2 3367 2001 0 3 3871 2002 0 4 4239 2003 0 5 4929 2004 0 6 5217 2005 0 7 1163 2000 1 8 1292 2001 1 9 1474 2002 1 10 1678 2003 1 11 1865 2004 1 12 NA 2005 1 > tail(base,12) y ai bj 25 7 2000 4 26 10 2001 4 27 NA 2002 4 28 NA 2003 4 29 NA 2004 4 30 NA 2005 4 31 21 2000 5 32 NA 2001 5 33 NA 2002 5 34 NA 2003 5 35 NA 2004 5 36 NA 2005 5 On peut alors commencer par utiliser le modèle Regression models based on log-incremental payments de Stavros Christofides, basé sur une modélisation lognormale, introduite initialement par Etienne de Vylder en 1978 (Markus en parle, en trois parties, sur son blog http://lamages.blogspot.ca/Barnett%20Zehnwirth) > reg1=lm(log(y)~ + as.factor(ai)+ + as.factor(bj),data=base) > summary(reg1) Call: lm(formula = log(y) ~ as.factor(ai) + as.factor(bj), data = base) Residuals: Min 1Q Median 3Q Max -0.26374 -0.05681 0.00000 0.04419 0.33014 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.9471 0.1101 72.188 6.35e-15 *** as.factor(ai)2001 0.1604 0.1109 1.447 0.17849 as.factor(ai)2002 0.2718 0.1208 2.250 0.04819 * as.factor(ai)2003 0.5904 0.1342 4.399 0.00134 ** as.factor(ai)2004 0.5535 0.1562 3.543 0.00533 ** as.factor(ai)2005 0.6126 0.2070 2.959 0.01431 * as.factor(bj)1 -0.9674 0.1109 -8.726 5.46e-06 *** as.factor(bj)2 -4.2329 0.1208 -35.038 8.50e-12 *** as.factor(bj)3 -5.0571 0.1342 -37.684 4.13e-12 *** as.factor(bj)4 -5.9031 0.1562 -37.783 4.02e-12 *** as.factor(bj)5 -4.9026 0.2070 -23.685 4.08e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1753 on 10 degrees of freedom (15 observations deleted due to missingness) Multiple R-squared: 0.9975, Adjusted R-squared: 0.9949 F-statistic: 391.7 on 10 and 10 DF, p-value: 1.338e-11 > base$py=exp(predict(reg1,
+ newdata=base)+summary(reg1)$sigma^2/2) > round(matrix(base$py,n,n),1)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 2871.2 1091.3 41.7 18.3  7.8 21.3
[2,] 3370.8 1281.2 48.9 21.5  9.2 25.0
[3,] 3768.0 1432.1 54.7 24.0 10.3 28.0
[4,] 5181.5 1969.4 75.2 33.0 14.2 38.5
[5,] 4994.1 1898.1 72.5 31.8 13.6 37.1
[6,] 5297.8 2013.6 76.9 33.7 14.5 39.3
> sum(base$py[is.na(base$y)])
[1] 2481.857

On obtient un montant un peu différent de celui obtenu par la méthode Chain Ladder, mais néanmoins comparable. On peut aussi tenter une régression de Poisson (avec un lien logarithmique), comme suggéré en 1975 par Hachemeister et Stanard,

> reg2=glm(y~
+ as.factor(ai)+
+ as.factor(bj),data=base,
+ family=poisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson,
data = base)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.3426  -0.4996   0.0000   0.2770   3.9355

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        8.05697    0.01551 519.426  < 2e-16 ***
as.factor(ai)2001  0.06440    0.02090   3.081  0.00206 **
as.factor(ai)2002  0.20242    0.02025   9.995  < 2e-16 ***
as.factor(ai)2003  0.31175    0.01980  15.744  < 2e-16 ***
as.factor(ai)2004  0.44407    0.01933  22.971  < 2e-16 ***
as.factor(ai)2005  0.50271    0.02079  24.179  < 2e-16 ***
as.factor(bj)1    -0.96513    0.01359 -70.994  < 2e-16 ***
as.factor(bj)2    -4.14853    0.06613 -62.729  < 2e-16 ***
as.factor(bj)3    -5.10499    0.12632 -40.413  < 2e-16 ***
as.factor(bj)4    -5.94962    0.24279 -24.505  < 2e-16 ***
as.factor(bj)5    -5.01244    0.21877 -22.912  < 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: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: 209.52

Number of Fisher Scoring iterations: 4

> base$py2=predict(reg2, + newdata=base,type="response") > > round(matrix(base$py2,n,n),1)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.7 1202.1 49.8 19.1  8.2 21.0
[2,] 3365.6 1282.1 53.1 20.4  8.8 22.4
[3,] 3863.7 1471.8 61.0 23.4 10.1 25.7
[4,] 4310.1 1641.9 68.0 26.1 11.2 28.7
[5,] 4919.9 1874.1 77.7 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.4 31.6 13.6 34.7
>
> sum(base$py2[is.na(base$y)])
[1] 2426.985

La prédiction coïncide avec l’estimateur obtenu par la méthode Chain Ladder. Le lien avec les méthodes de biais minimal a été établi par  Klaus Schmidt et Angela Wünsche en 1998, dans Chain ladder, marginal sum and maximum likelihood estimation. La semaine prochaine, on parlera des méthodes de bootstrap pour obtenir des intervalles de confiance, ou des quantiles, sur les montants de réserve. Je ne sais pas si j’aurais le temps de taper des transparents, je préfère, sur cette partie du cours taper au fur et à mesure, et écrire au tableau. Je renvoie au Chapitre 3 du livre avec Christophe Dutang – en ligne sur http://cran.r-project.org/doc/contrib/ – pour le détail. C’est le code que je tape en cours, tout en essayant de répondre aussi aux questions.

# Examen intra, éléments de correction

L’énoncé de l’examen intra est en pdf ici et comme annoncé par courriel, la correction de l’intra est dans le pdf en ligne. Comme personne ne semble en désaccord avec les réponses proposées, les notes seront mises en ligne très bientôt. Concertant les questions 18 et 19 quelques compléments d’explications (que je n’avais pas tapé dans le pdf). On avait vu que l’estimateur du maximum de vraisemblance pour une régression de Poisson était asymptotiquement Gaussien,

$\widehat{\boldsymbol{\beta}}_{P}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{P}))$

(asymptotiquement) avec

$V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}$

Quand on a une régression de type binomiale négative, si on note de manière très générale $\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i)$ (on avait vu en cours qu’il existait plusieurs spécifications possibles pour cette variance conditionnelle). Dans ce cas,

$\widehat{\boldsymbol{\beta}}_{BN}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{BN}))$

avec

$V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}\left[\sum_{i=1}^n \omega_i \boldsymbol{X}_i\boldsymbol{X}_i\right]\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}$

Bref, tout dépend fondamentalement de la spécification de la variance conditionnelle. Sous R, c’est la régression binomiale négative de type 1 qui est considérée, i.e.

$\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i)=\phi\cdot \mathbb{E}(Y_i|\boldsymbol{X}_i)=\phi \cdot \widehat{Y}_i$

On toujours une relation de la forme

$\widehat{\boldsymbol{\beta}}_{QP}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{QP}))$

avec (en simplifiant un peu)

$V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}$

aussi, on a

$V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot V_\infty(\widehat{\boldsymbol{\beta}}_{P})$

Mais comme annoncé en cours, des points étaient données pour ceux qui se contentaient d’affirmer que la variance des estimateurs était plus grande s’il y avait sur-dispersion.

# Examen intra, régression logistique et de Poisson

L’examen intra du cours ACT2040 aura lieu mercredi matin, de 9:00 à 12:00. Aucun document autorisé, sauf les calculatrices (modèle standard, cf plan de cours), et les téléphones seront formellement interdits. Il y aura 34 questions portant sur la première partie du cours (jusqu’à la fin des modèles de comptage, sections 1 à 5 des transparents). 15 questions porteront sur la base décrite dans un précédant billet, sur le nombre de relations extra-conjugales. Il s’agira de décrire les sorties en ligne ici. Je laisse 36 heures pour prendre connaissance de ces sorties. Une version sera donnée lors de l’examen (imprimée 2 pages par feuille, comme dans la version en ligne: si quelqu’un a besoin d’un exemplaire imprimé plus gros, merci de me le faire savoir avant l’examen).

# Multiple (smoothed) regression and portfolio exposure

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

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

(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,

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,

(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

# Données pour la régression logistique, et de Poisson

Pour le cours de mercredi, deux petites bases, pour se pratiquer à modéliser des variable 0/1 ou une variable de comptage,

ou encore

Sinon, une base plus complète pour faire de la tarification,

ageconducteur agepermis sexeconducteur situationfamiliale  habitation zone
1            57        39              F             Celiba peri-urbain    8
2            54        35              H             Celiba      urbain    3
3            51        32              F             Celiba      urbain    1
4            53        35              H              Marie       rural    4
5            61        43              H              Marie      urbain    8
6            60        29              F              Marie peri-urbain    1
agevehicule proprietaire    payment  marque         poids     usage
1          12    locataire     Annuel  AUTRES     8.>3500kg PROMENADE
2          20     sans mrp Semestriel PEUGEOT 4.3100-3199kg PROMENADE
3           4     sans mrp     Annuel  RAPIDO     1.<2700kg PROMENADE
4           1     sans mrp     Annuel  AUTRES 3.3000-3099kg PROMENADE
5           1 proprietaire     Annuel    FIAT 6.3300-3399kg PROMENADE
6          10     sans mrp    Mensuel    FIAT     8.>3500kg PROMENADE
exposition nombre   voiture
1          1      0 Monospace
2          1      0   Berline
3          1      0  sans avp
4          1      0  sans avp
5          1      1 Monospace
6          1      0  sans avp

Parmi les variables, la description (sommaire) est la suivante,

• ageconducteur: âge du conducteur principal du véhicule
• agepermis: ancienneté du permis de conduire du conducteur principal du véhicule
• sexeconducteur: sexe du conducteur principal (H ou F)
• situationfamiliale: situation familiale du conducteur principal (“Celiba”, “Marie” ou “Veuf/Div”)
• habitation: zone d’habitation du conducteur principal (“peri-urbain”, “rural” ou “urbain” )
• zone: zone d’habitation (allant de 1 à 8)
• agevehicule: age du véhicule
• proprietaire: si le conducteur principal possède un contrat Habitation, son statut (“locataire” ou “proprietaire”)  Sinon “sans mrp”
• payment:type de fractionnement de la prime d’assurance automobile (“Annuel”, “Mensuel” ou “Semestriel”)
• marque: marque du véhicule
> levels(BASEN[,10])
[1] "ADRIA"       "AUTOSTAR"    "AUTRES"      "BURSTNER MOBIL"
[5] "CHALLENGER"  "CHAUSSON"    "CITROEN"     "FIAT"
[9] "FORD"        "HYMERMOBIL"  "MERCEDES"    "PEUGEOT"
[13] "PILOTE"     "RAPIDO"      "RENAULT"     "VOLKSWAGEN"
• poids: classe de poids du véhicule
> levels(BASEN[,11])
[1] "1.<2700kg"    "2.2700-2999kg""3.3000-3099kg""4.3100-3199kg"
[5] "5.3200-3299kg""6.3300-3399kg""7.3400-3499kg""8.>3500kg"
• usage: utilisation du véhicule principal (“PROMENADE” ou “TOUS_DEPLACEMENTS”)
• exposition: exposition, en années
• nombre: nombre d’accident responsabilité civile du conducteur principal, pendant l’année passée
• cout: cout du sinistre
• voiture: type de véhicule
> levels(BASEN[,15])
[1] "Berline"            "Break"              "Buggy"
[4] "Cabriolet"          "Combispace"         "Coup\xe9"
[7] "Coup\xe9 Cabriolet" "Jeep"               "Minibus"
[10] "Minispace"          "Monospace"         "sans avp"

La variable d’intérêt est ici le nombre d’accident,

> table(BASEN$nombre) 0 1 60155 3264 La base est un peu particulière – on en parlera en classe – les assurés ayant eu 0 ou 1 accident dans l’année. # Further readings on GLMs and ratemaking Some articles found in Actuarial journal, on ratemarking, and in the CAS forums, and Astin conference papers # Modélisation des coûts individuels en tarification Avant de terminer le cours sur la tarification, on va parler de la modélisation des coûts individuels. On parlera de lois Gamma et de lois lognormales (sur cette dernière, je suggère de relire ce qui avait été dit dans le cours de modèles de régression sur les modèles log-linéaires, rappelé dans un court billet publié à l’automne). On parlera aussi de mélanges de lois, et de lois multinomiales. Les transparents sont en ligne ici, Pour aller plus loin, il y a l’article de Fu & Moncher (2004) sur la comparaison Gamma versus lognormale, http://casact.org/… ou Holler, Sommer & Trahair (1999) http://casact.org/… qui proposait un état de l’art, il y a une quinzaine d’années. Sinon, je recommande la lecture du Practitioner’s Guide to Generalized Linear Models, en ligne sur http://casact.org/…. # Modeling individual losses with mixtures Usually, the sentence that I keep saying in my regression classes is “please, look at your data“. In our previous post, we’ve been playing like most econometricians: we did not look at the data. Actually, if we look at the distribution of individual losses, in the dataset, we see the following, > n=nrow(couts) > plot(sort(couts$cout),(1:n)/(n+1),xlim=c(0,10000),type="s",lwd=2,col="green")

It looks like there are fixed costs claims in our database. How do we deal with it in the standard case (e.g. in Loss Models textbook) ? We can use a mixture of – at least – three distributions here,

$f(y\boldsymbol) = p_1 {\color{Blue} f_1(}y{\color{Blue} )} + p_2 {\color{Magenta} \delta_{\kappa}(}y{\color{Magenta} )} + p_3 {\color{Red} f_3(}y{\color{Red} )}$

with

• a distribution for small claims, ${\color{Blue} f_1(}\cdot{\color{Blue} )}$, e.g. an exponential distribution
• a Dirac mass in ${\color{Magenta} \kappa}$, i.e. ${\color{Magenta} \delta_{\kappa}(}\cdot{\color{Magenta} )}$
• a distribution for larger claims, ${\color{Red} f_3(}\cdot{\color{Red} )}$, e.g. a Gamma, or a lognormal, distribution
>  I1=which(couts$cout<1120) > I2=which((couts$cout>=1120)&(couts$cout<1220)) > I3=which(couts$cout>=1220)
>  (p1=length(I1)/nrow(couts))
[1] 0.3284823
>  (p2=length(I2)/nrow(couts))
[1] 0.4152807
>  (p3=length(I3)/nrow(couts))
[1] 0.256237
>  X=couts$cout > (kappa=mean(X[I2])) [1] 1171.998 > X0=X[I3]-kappa > u=seq(0,10000,by=20) > F1=pexp(u,1/mean(X[I1])) > F2= (u>kappa) > F3=plnorm(u-kappa,mean(log(X0)),sd(log(X0))) * (u>kappa) > F=F1*p1+F2*p2+F3*p3 > lines(u,F) In our previous post, we’ve discussed the idea that all parameters might be related to some covariates, i.e. $f(y|\boldsymbol{X}) = p_1(\boldsymbol{X}) {\color{Blue} f_1(}y|\boldsymbol{X}{\color{Blue} )} + p_2(\boldsymbol{X}) {\color{Magenta} \delta_{\kappa}(}y{\color{Magenta} )} + p_3(\boldsymbol{X}) {\color{Red} f_3(}y|\boldsymbol{X}{\color{Red} )}$ which yield the following premium model, $\mathbb{E}(Y|\boldsymbol{X}) = {\color{Blue} {\underbrace{\mathbb{E}(Y|\boldsymbol{X},Y\leq s_1)}_{A} \cdot {\underbrace{\mathbb{P}(Y\leq s_1|\boldsymbol{X})}_{D}}}}\\+{\color{Purple} {{\underbrace{\mathbb{E}(Y|Y\in( s_1,s_2], \boldsymbol{X}) }_{B}}\cdot {\underbrace{\mathbb{P}(Y\in( s_1,s_2]| \boldsymbol{X})}_{D}}}}\\+{\color{Red} {{\underbrace{\mathbb{E}(Y|Y> s_2, \boldsymbol{X}) }_{C}}\cdot {\underbrace{\mathbb{P}(Y> s_2| \boldsymbol{X})}_{D}}}}$ For the ${\color{Blue} A}$${\color{Magenta} B}$ and ${\color{Red} C}$ terms, that’s easy, we can use standard models we’ve seen in the course. For the probability, we should use a multinomial model. Recall that for the logistic regression model, if $(\pi,1-\pi)=(\pi_1,\pi_2)$, then $\log \frac{\pi}{1-\pi}=\log \frac{\pi_1}{\pi_2} =\boldsymbol{X}'\boldsymbol{\beta}$ i.e. $\pi_1 = \frac{\exp(\boldsymbol{X}'\boldsymbol{\beta})}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta})}$ and $\pi_2 = \frac{1}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta})}$ To derive a multivariate extension, write $\pi_1 = \frac{\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}$ $\pi_2 = \frac{\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}$ and $\pi_3 = \frac{1}{1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)}$ Again, maximum likelihood techniques can be used, since $\mathcal{L}(\boldsymbol{\pi},\boldsymbol{y})\propto \prod_{i=1}^n \prod_{j=1}^3 \pi_{i,j}^{Y_{i,j}}$ where here, variable $Y_{i}$ – which take three levels – is splitted in three indicators (like any categorical explanatory variables in standard regression model). Thus, $\log \mathcal{L}(\boldsymbol{\beta},\boldsymbol{y})\propto \sum_{i=1}^n \sum_{j=1}^2 \left(Y_{i,j} \boldsymbol{X}_i'\boldsymbol{\beta}_j\right) - n_i\log\left[1+1+\exp(\boldsymbol{X}'\boldsymbol{\beta}_1)+\exp(\boldsymbol{X}'\boldsymbol{\beta}_2)\right]$ and, as for the logistic regression, then use Newton Raphson’ algorithm to compute numerically the maximum likelihood. In R, first we have to define the levels, e.g. > seuils=c(0,1120,1220,1e+12) > couts$tranches=cut(couts$cout,breaks=seuils, + labels=c("small","fixed","large")) > head(couts,5) nocontrat no garantie cout exposition zone puissance agevehicule 1 1870 17219 1RC 1692.29 0.11 C 5 0 2 1963 16336 1RC 422.05 0.10 E 9 0 3 4263 17089 1RC 549.21 0.65 C 10 7 4 5181 17801 1RC 191.15 0.57 D 5 2 5 6375 17485 1RC 2031.77 0.47 B 7 4 ageconducteur bonus marque carburant densite region tranches 1 52 50 12 E 73 13 large 2 78 50 12 E 72 13 small 3 27 76 12 D 52 5 small 4 26 100 12 D 83 0 small 5 46 50 6 E 11 13 large Then, we can run a multinomial regression, from > library(nnet) using some selected covariates > reg=multinom(tranches~ageconducteur+agevehicule+zone+carburant,data=couts) # weights: 30 (18 variable) initial value 2113.730043 iter 10 value 2063.326526 iter 20 value 2059.206691 final value 2059.134802 converged The output is here > summary(reg) Call: multinom(formula = tranches ~ ageconducteur + agevehicule + zone + carburant, data = couts) Coefficients: (Intercept) ageconducteur agevehicule zoneB zoneC fixed -0.2779176 0.012071029 0.01768260 0.05567183 -0.2126045 large -0.7029836 0.008581459 -0.01426202 0.07608382 0.1007513 zoneD zoneE zoneF carburantE fixed -0.1548064 -0.2000597 -0.8441011 -0.009224715 large 0.3434686 0.1803350 -0.1969320 0.039414682 Std. Errors: (Intercept) ageconducteur agevehicule zoneB zoneC zoneD fixed 0.2371936 0.003738456 0.01013892 0.2259144 0.1776762 0.1838344 large 0.2753840 0.004203217 0.01189342 0.2746457 0.2122819 0.2151504 zoneE zoneF carburantE fixed 0.1830139 0.3377169 0.1106009 large 0.2160268 0.3624900 0.1243560 To visualize the impact of a covariate (one, only), one can use also spline functions > library(splines) > reg=multinom(tranches~agevehicule,data=couts) # weights: 9 (4 variable) initial value 2113.730043 final value 2072.462863 converged > reg=multinom(tranches~bs(agevehicule),data=couts) # weights: 15 (8 variable) initial value 2113.730043 iter 10 value 2070.496939 iter 20 value 2069.787720 iter 30 value 2069.659958 final value 2069.479535 converged For instance, if the covariate is the age of the car, we do have the following probabilities > predict(reg,newdata=data.frame(agevehicule=5),type="probs") small fixed large 0.3388947 0.3869228 0.2741825 and for all ages from 0 to 20, For instance, for new cars, the proportion of fixed costs is rather small (here in purple), and keeps increasing with the age of the car. If the covariate is the density of population in the area the driver lives, we do obtain the following probabilities > reg=multinom(tranches~bs(densite),data=couts) # weights: 15 (8 variable) initial value 2113.730043 iter 10 value 2068.469825 final value 2068.466349 converged > predict(reg,newdata=data.frame(densite=90),type="probs") small fixed large 0.3484422 0.3473315 0.3042263 Based on those probabilities, it is then possible to derive the expected cost of a claims, given some covariates (e.g. the density). But first, define subsets of the whole dataset > sbaseA=couts[couts$tranches=="small",]
> sbaseB=couts[couts$tranches=="fixed",] > sbaseC=couts[couts$tranches=="large",]

with a threshold given by

> (k=mean(sousbaseB$cout)) [1] 1171.998 Then, let us run our four models, > reg=multinom(tranches~bs(densite),data=couts) > regA=glm(cout~bs(densite),data=sousbaseA,family=Gamma(link="log")) > regB=glm(cout~1,data=sousbaseB,family=Gamma(link="log")) > regC=glm((cout-k)~bs(densite),data=sousbaseC,family=Gamma(link="log")) We can now compute predictions based on those models, > nouveau=data.frame(densite=seq(10,100)) > proba=predict(reg,newdata=nouveau,type="probs") > predA=predict(regA,newdata=nouveau,type="response") > predB=predict(regB,newdata=nouveau,type="response") > predC=predict(regC,newdata=nouveau,type="response")+k > pred=cbind(predA,predB,predC) To visualize the impact of each component on the premium, we can compute probabilities, are well as expected costs (given a cost in each subset), > cbind(proba,pred)[seq(10,90,by=10),] small fixed large predA predB predC 10 0.3344014 0.4241790 0.2414196 423.3746 1171.998 7135.904 20 0.3181240 0.4471869 0.2346892 428.2537 1171.998 6451.890 30 0.3076710 0.4626572 0.2296718 438.5509 1171.998 5499.030 40 0.3032872 0.4683247 0.2283881 451.4457 1171.998 4615.051 50 0.3052378 0.4620219 0.2327404 463.8545 1171.998 3961.994 60 0.3136136 0.4417057 0.2446807 472.3596 1171.998 3586.833 70 0.3279413 0.4056971 0.2663616 473.3719 1171.998 3513.601 80 0.3464842 0.3534126 0.3001032 463.5483 1171.998 3840.078 90 0.3652932 0.2868006 0.3479061 440.4925 1171.998 4912.379 Now, it is possible to plot those figures in a graph, > barplot(t(proba*pred)) > abline(h=mean(couts$cout),lty=2)

(the dotted horizontal line is the average cost of a claim, in our dataset).

# 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.