Moyennes de Prédicteurs

La dernière fois, nous avions vu 4 modèles de séries temporelles, pour modéliser la série hebdomadaire de températures. Le premier était obtenu avec auto.arima, les deux étaient ensuite ses modèles SARIMA (le second était une forte autocorélation à l’ordre 52, et le troisième avait une racine unité saisonnière), et le dernier était obtenu par la méthode de Buys-Ballot (cyclique, avec un bruit autocorrélé).

Au lieu de se demander lequel des quatre est le “meilleur”, on peut tenter de trouver une moyenne pondérée, \hat{y}_t = \sum_{i}\omega_i \hat{y}_t^{(i)}\omega=(\omega_i) est un point du simplexe, au sens où \omega_i\geq 0 et \omega_1+...+\omega_k=1. On veut ensuite trouver les “meilleurs” pondérations, par exemple au sens de la norme L^2. On va calibrer nos quatre modèles sur les m premières valeurs, et on compare ensuite la combinaison des prévions sur les 111 prochaines valeurs (et les vraies valeurs), \omega^\star=\text{argmin}\{\sum_{h=1}^{111}\left(y_{t+h}- \sum_{i}\omega_i \hat{y}_{t+h}^{(i)}\right)^2\}

Par exemple, on peut utiliser les 200 premières valeurs.

m=200
Ym=ts(y[1:m],start=1996,frequency=52)

On va ensuite estimer les 4 modèles sur ces 200 valeurs

model1=arima(Ym,order = c(4,0,0), 
seasonal = list(order = c(1,0,0))) 
model2=arima(Ym,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
model3=arima(Ym,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
regL=lm(Temp~0+as.factor(NumWeek),
data=electricite[1:m,])
E=residuals(regL)
model4=arima(ts(E,start=1996,
frequency=52),order = c(1,0,0))

On peut ensuite collecter les prédictions

DONNEES = data.frame(
  y1=predict(model1,n.ahead = 111)$pred,
  y2=predict(model2,n.ahead = 111)$pred,
  y3=predict(model3,n.ahead = 111)$pred,
  y4=predict(model4,n.ahead = 111)$pred+
     predict(regL,newdata=electricite[m+(1:111),]),
  vraie=electricite[m+(1:111),"Temp"])

Pour créer une combinaison linéaire des prévisions, on utilise

a=rep(1/4,4)
y_prev = as.matrix(DONNEES[,1:4]) %*% a

On peut ainsi visualiser les 4 prévisions, leur combinaison linéaire (ici avec des équi-pondérations), et leurs observations

plot(1:111,DONNEES[,5],type="l",lwd=2)
for(i in 1:4) lines(1:111,DONNEES[,i],col="blue")
lines(1:111,yp,col="red")

Pour trouver la valeur “optimale” des pondérations, qui minimisera la somme des carrés des erreurs, on utilise le code suivant (on fera une transformation des données, ici sur le simplexe, pour s’assurer que l’optimisation donnera un élément du simplexe)

erreur=function(a) sum((DONNEES[,1:4] 
%*% a-DONNEES[,5])^2)
erreur2=function(b){
  a=rep(NA,length(b)+1)
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  erreur(a)
}

On obtient l’optimum, et on peut obtenir une prévision

b=optim(par=c(0,0,0),erreur2)$par
a=rep(NA,length(b)+1)
a[1:length(b)]=exp(b)/(1+sum(exp(b)))
a[length(b)+1]=1/(1+sum(exp(b)))
yp=DONNEES[,1:4] %*% a
lines(1:111,yp,col="red",lwd=4)

On peut alors faire changer le nombre d’observations utilisées

Techniquement, ca pose un petit soucis, car il faut s’assurer la convergence de deux algorithmes : celui d’estimation des paramètres des SARIMA, celui de recherche du paramètre de poids.

poids_optimaux=function(m){
  Ym=ts(y[1:m],start=1996,frequency=52)
  TRY=try(model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0))),FALSE)
  if(inherits(TRY, "try-error")) model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0)),method="CSS")
  TRY=try(model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0))),FALSE)
  if(inherits(TRY, "try-error")) model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0)),method="CSS")
  TRY=try(model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1))),FALSE)
  if(inherits(TRY, "try-error")) model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1)),method="CSS")
regL=lm(Temp~0+as.factor(NumWeek),data=electricite[1:m,])
  E=residuals(regL)
model4=arima(ts(E,start=1996,frequency=52),order = c(1,0,0))
  DONNEES = data.frame(
    y1=predict(model1,n.ahead = 111)$pred,
    y2=predict(model2,n.ahead = 111)$pred,
    y3=predict(model3,n.ahead = 111)$pred,
    y4=predict(model4,n.ahead = 111)$pred+
      predict(regL,newdata=electricite[m+(1:111),]),
    vraie=electricite[m+(1:111),"Temp"])
  erreur=function(a) sum((as.matrix(DONNEES[,1:4]) %*% a-DONNEES[,5])^2)
  erreur2=function(b){
    a=rep(NA,length(b)+1)
    a[1:length(b)]=exp(b)/(1+sum(exp(b)))
    a[length(b)+1]=1/(1+sum(exp(b)))
    erreur(a)
  }
  TRY=try(b <- optim(par=c(0,0,0),erreur2)$par, FALSE)
  a=rep(NA,length(b)+1)
  if(!inherits(TRY, "try-error")) {
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  }
  return(a)
}

On peut alors regarder l’évolution des pondérations au cours du temps. On peut vectoriser, mais si l’algorithme aboutit sur une erreur, c’est plus simple d’identifier

VM=seq(200,480,by=10)
VA=matrix(NA,length(VM),4)
for(i in 1:length(VM)) VA[i,]=poids_optimaux(VM[i])
colnames(VA)=c("auto","SARMA","SARIMA","Buys Ballot")

(pour débugger, c’est plus simple les boucles). On peut alors obtenir le graphique suivant, avec en rose Buys-Ballot, et en vert le modèle SARIMA, avec une racine unité saisonnière,

tva=t(VA)
barplot(tva,legend = rownames(counts),
 col=c("yellow","orange","green","pink"))

On se rend compte qu’assez vite, le modèle le plus intéressant est le modèle de Buys Ballot.

Il est possible de changer de fonction de perte, par exemple si on s’intéresse à un quantile à 90%,

tau=.9
loss=function(e) e*(tau-(e<=0)*1)

Dans la fonction, on utilise alors

erreur=function(a) sum(loss(as.matrix(
DONNEES[,1:4]) %*% a-DONNEES[,5]))

C’est, cette fois, beaucoup plus incertain, même si les deux modèles les plus intéressants sont toujours le SARIMA et Buys-Ballot.

Insurance: Risk Pooling and Price Segmentation

Talk this afternoon at the seminar of Telecom ParisTech

Insurance is usually defined as “the contribution of the many to the misfortune of the few”. This idea of pooling risks together using the law of large number legitimates the use of  the expected value as actuarial “fair” premium. In the context of heterogeneous risks, nevertheless, it is possible to legitimate price segmentation based on observable characteristics. But nowadays, intensive segmentation can be observed, with a much wider range of offered premium, on a given portfolio. In this talk, we will briefly get back on statistical approaches of insurance pricing (classical econometric tools vs machine learning). We will then get back on recent experiments (so-called “actuarial pricing game”) organized since 2015, where real actuaries are playing in competitive (artificial) market, that mimic real insurance market. We will get back on conclusions obtained on two editions, the first one, and the most recent one, where a dynamic version of the game was launched.

By the way, there will be soon a fourth version of the “Actuarial Pricing Game”… some information soon, on this blog…