Modéliser la consommation électrique

A partir de nos prévisions pour la température, on peut tenter de prévoir la consommation électrique. Rappelons que la série de consommation électrique ressemble à ca

plot(electricite[passe,"Load"],type="l")

On peut tenter un modèle assez simple, où la consommation à la date Y_t est fonction d’une tendance linéaire a+bt, de la position dans l’année (sous une forme non-paramétrique), de la température T_t (sous une forme non-linéaire, disons un polynôme de degré 3), de la températeure en t-1, T_{t-1} et de l’indice de production industriel IPI_t (sous une forme linéaire pour faire simple).

model1 = lm(Load~1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

L’idée de la fonction polynomial pour l’impact de la température venait du graphique suivant (où on représente la série de la consommation, une fois enlevée la tendance linéaire)

On peut aussi supposer une forme autorégressive, où Y_{t} sera fonction de Y_{t-1}

model2 = lm(Load~1+Load1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

On peut alors tenter de faire de la prévision. Le soucis avec ce second modèle, c’est la composante autorégressive. Pour prévoir Y_{t+h}, il faut utiliser la prévision faite en t+h-1, \hat Y_{t+h-1}.

base_prevision=data.frame(
Load_prev = NA,  Load1=c(electricite[length(passe),"Load"],rep(NA,110),
Time = electricite[futur,"Time"],
NumWeek = electricite[futur,"NumWeek"],
Temp = MODELES[,"y4"], 
Temp1 = c(electricite[length(passe),"Temp"],
MODELES[1:110,"y4"]),
IPI = electricite[futur,"IPI"])
for(t in 1:110){
 p = predict(model2,newdata=base_prevision)[t]
 base_prevision[t,"Load_prev"] = p
 base_prevision[t+1,"Load1"] = p}

On peut alors comparer la prévision \hat Y_{t} avec l’observation Y_{t}.

plot(1:111,electricite[futur,"Load"],type="l")
lines(1:111,base_prevision[,"Load_prev"],col="red")

On est bon les mois d’été (on capture bien le pic de la première quinzaine du mois d’août) par contre, on sous-estime fortement la consommation l’hiver. Et c’est gênant…

Pour finir, on peut oublier le passage par les variables explicatives, et tenter directement un modèle de série temporelles.

plot(electricite[passe,"Load"],type="l")

Le soucis, c’est le côté hétéroscédastique de la série, avec une pente sur le minimum plus faible que sur le maximum.

n=length(passe)
plot(1:n,electricite[passe,"Load"],type="l")
m=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),min)
M=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),max)
vm=which(electricite[passe,"Load"]%in%m$x)
vM=which(electricite[passe,"Load"]%in%M$x)
points(sort((1:n)[vm]),m$x,pch=19,col="blue")
points(sort((1:n)[vM]),M$x,pch=19,col="blue")
xm=((1:n)[vm])
xM=((1:n)[vM])
regm=lm(m$x~xm,col="blue")
regM=lm(M$x~xM,col="blue")
abline(regm,lty=2)
abline(regM,lty=2)

Une solution classique (et simple) est de passe au logarithme

plot(electricite[passe,"Load"],type="l",log="y")
y=electricite[passe,"Load"]
z = log(y)
plot(z,type="l")
B = data.frame(z=z,t=1:length(z))
abline(lm(z~t,data=B),lty=2,col="red")

Il faut ensuite enlever la tendance linéaire, pour stationnariser la série

z = residuals(lm(z~t,data=B))
Z = ts(z,start=1996,frequency=52)
modelz1 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
modelz1 = arima(Z,order = c(4,0,0), 
seasonal = list(order = c(1,0,1)))

Ce premier modèle est stationnaire, sans racine unité. On peut tenter d’introduire une racine unité saisonnière

modelz2 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(0,1,0)))
modelz2 = arima(Z,order = c(1,0,3), 
seasonal = list(order = c(1,1,0)))

Enfin, le dernier sera un peu plus simple

modelz3 = arima(Z,order = c(1,0,0), 
seasonal = list(order = c(2,0,0)))

On stocke alors toutes les prévisions dans une base

DONNEES_u = data.frame(
  u1 = predict(modelz1,n.ahead = 111)$pred,
  u2 = predict(modelz2,n.ahead = 111)$pred,
  u3 = predict(modelz3,n.ahead = 111)$pred)

On rajoute alors la tendance linéaire à la prévision du résidu

reg = lm(z~t,data=B)
dfutur = data.frame(t=futur)
DONNEES_z = data.frame(
  z1 = DONNEES_u$u1+predict(reg,newdata=dfutur),
  z2 = DONNEES_u$u2+predict(reg,newdata=dfutur),
  z3 = DONNEES_u$u3+predict(reg,newdata=dfutur)

On a besoin de la variance pour la prévision. En effet, classiquement, la prévision sera l’espérance de notre grandeur d’intérêt. Ici, on a fait un modèle linéaire sur \log Y, i.e. \log Y \sim\mathcal{N}(\mu,\sigma^2) et donc \mathbb{E}[Y]=\exp\left(\mu+\sigma^2/2\right)

sigma=summary(reg)$sigma
DONNEES_seu = data.frame(
  seu1=sqrt(predict(modelz1,n.ahead = 111)$se^2+sigma^2),
  seu2=sqrt(predict(modelz2,n.ahead = 111)$se^2+sigma^2),
  seu3=sqrt(predict(modelz3,n.ahead = 111)$se^2+sigma^2))

On fait ici l’hypothèse que les estimateurs des prévisions des deux modèles (la tendance linéaire et celle du modèle autorégressif) sont indépendants et donc on peut sommer les termes de variance. Aussi, les prévisions pour Y sont

DONNEES_y = data.frame(
  y1 = exp(DONNEES_z$z1+1/2*DONNEES_seu$seu1^2),
  y2 = exp(DONNEES_z$z2+1/2*DONNEES_seu$seu2^2),
  y3 = exp(DONNEES_z$z3+1/2*DONNEES_seu$seu3^2))

On peut représenter les prévisions des trois modèles (avec les observations)

plot(futur,electricite[futur,"Load"],type="l")
lines(futur,DONNEES_y$y1,col="red")
lines(futur,DONNEES_y$y2,col="blue")
lines(futur,DONNEES_y$y3,col="purple")

On peut comparer avec la prévision précédente,

lines(futur,base_prevision[,"Load_prev"],
col="orange")

On y perd un peu l’été, par contre, on s’est un peu amélioré l’hiver.

Forecasting Techniques for Demographics

Tomorrow, we will discuss in our cours forecasting tools for demographics. But first, we will see basic static tools. Before playing with longitudinal dataset, let us use “standard” life tables. Some (French) datasets are available on the INED website, but let us use the most popular ones, the TV8890 and TD8890. Those two tables can be downloaded from wikipedia,

url="https://fr.wikipedia.org/wiki/Table_de_mortalit%C3%A9"
download.file(url,"mortalite.html")
library(XML)
tables=readHTMLTable("mortalite.html")

The code to get the dataset is the following (here, we should be careful since there is a space in the number, as a separator for thousands)

TV8890_0=tables[[2]]
a1=as.numeric(as.character(TV8890_0[,1]))
a2=as.numeric(as.character(TV8890_0[,3]))
espace=function(x) gsub("[[:space:]]", "", x)
b1=espace(as.character(TV8890_0[,2]))
b2=espace(as.character(TV8890_0[,4]))
TV8890=data.frame(x=c(a1,a2),lx=as.numeric(c(b1,b2))

One can also read the second life table (note that there is a typo in wikipedia since here, the two tables are exactly the same)

TD8890_0=tables[[1]]
a1=as.numeric(as.character(TD8890_0[,1]))
a2=as.numeric(as.character(TD8890_0[,3]))
b1=espace(as.character(TD8890_0[,2]))
b2=espace(as.character(TD8890_0[,4]))
TD8890=data.frame(x=c(a1,a2),lx=as.numeric(c(b1,b2)))

It is possible to use that survival function to compute some sort to life expectancy at birth

sum(TV8890$lx)/100000-1
[1] 72.01518

One can visualize the survival probability (up to a 100,000 scaling constant)

plot(TV8890,type="l")

or the death probability, i.e. the probability to die at some specific age x, given that you did reach age x, also called the force of mortality

n=nrow(TV8890)
px=(TV8890$lx[1:(n-1)]-TV8890$lx[2:n])/
TV8890$lx[1:(n-1)]
x=TV8890$x[1:(n-1)]
plot(x,px,type="l",xlab="age")

A more popular visualization is obtained with a log scale for the probability

plot(x,px,type="l",log="y")

Finally, we can compute the density of the age at death

pbx=TV8890$lx[1:(n-1)]*px/100000
plot(x,pbx,type="l")

that can also be used to compute life expectancy

sum(x*pbx)
[1] 72.01518

That is for the static case. For longitudinal tables, we can use those from the Human Mortality Database. For instance, for France, we can get information of the number of death at age x during year t, as well as the exposure (number of people alive). For France, there are datasets available here

url="http://freakonometrics.free.fr/FranceDeaths_1x1.txt"
download.file(url,"FRD.txt")
url="http://freakonometrics.free.fr/FranceExposures_1x1.txt"
download.file(url,"FRE.txt")

and for Canada

url="http://freakonometrics.free.fr/CanadaDeaths_1x1.txt"
download.file(url,"CAD.txt")
url="http://freakonometrics.free.fr/CanadaExposures_1x1.txt"
download.file(url,"CAE.txt")

The following code can be use to read those files.

FRD=read.table("FRD.txt",skip = 3,header=TRUE)
tail(FRD)
      Year  Age Female  Male  Total
22195 2015  105 297.96 35.86 333.82
22196 2015  106 182.95 20.39 203.34
22197 2015  107 104.87 10.50 115.37
22198 2015  108  57.27  5.07  62.34
22199 2015  109  31.93  2.59  34.52
22200 2015 110+  33.03  1.61  34.64

Il est trop nul ce prof…

Que celui qui n’a jamais pensé ca une fois dans sa scolarité me jette la première pierre…

Cela dit, il y a ceux qui le pensent, et ceux qui l’affirment haut et fort, genre sur des sites comme rate my professors. J’ai toujours eu du mal avec ce genre de sites (et plus généralement extraire une information utile des évaluations des élèves) car on ne connaît pas les motivations des élèves. Quand on lit “l’examen final est impossible” puis le commentaire d’après “l’examen, c’est juste des exercices vus en classe, c’est vraiment trop facile“, on finit par être perplexe. Oui, assez souvent une classe est composée d’éléments très hétérogènes, et il est impossible de contenter tout le monde. Le plus étonnant quand on enseigne, c’est l’équilibre qui peut en résulter. Le même cours peut être une réussite, et le semestre suivant un fiasco, alors que le contenu et la forme peuvent être identiques. Je me souviens quand j’étais étudiant de ce prof, annoncé par la promo précédante comme le meilleur prof au monde, que je trouvais finalement so-so.

Et puis un jour, certains d’entre nous décident d’aller un cran plus loin, de sauter le pas, de passer d’élève à prof. Sans pour autant que ce soit une vocation, juste un concours de circonstances. Pour ma part, j’avais travaillé en salle de marché, puis comme actuaire à Hong Kong, avant de revenir à la Fédération Francaise des Assurances sur Paris. Et juste avant les vacances de Noël (il y a maintenant une quinzaine d’années) on me propose de donner un cours trois semaines plus tard, au retour des vacances (prof qui lâche son cours au dernier moment, impossible de trouver un remplacant). Je pensais faire de la recherche en entreprise, et j’ai découvert l’enseignement.

Et je n’ai pas arrêté depuis. Avec toutes sortes de formats, toutes sortes de publics… En particulier, j’ai fait beaucoup de formation continue. Et j’en profite pour la minute culturelle du blog : en avril 1792, l’Assemblée législative nationale discute de l’organisation générale de l’instruction publique en France, et Jean Antoine Nicolas Caritat, plus connu sous le nom de marquis de Condorcet, fait un discours très avant-gardiste

Nous avons observé que l’instruction ne devait pas abandonner les individus au moment où il sortent des écoles ; qu’elle devait embrasser tous les âges ; qu’il n’y en avait aucun où il ne fût utile et possible d’apprendre, et que cette seconde instruction est d’autant plus nécessaire, que celle de l’enfance a été resserrée dans des bornes plus étroites. (…) En continuant l’instruction pendant toute la durée de la vie, on empêchera les connaissances acquises dans les écoles de s’effacer promptement de la mémoire, on entretiendra dans les esprits l’activité utile (…), on pourra montrer enfin l’art de s’instruire par soi-même

Or justement, il y a plusieurs mois, le centre de formation Caritat (dédié au monde de l’actuariat de l’assurance et de la protection sociale) lancait un prix pédagogique. L’idée est intéressante : vous pensez qu’il manque de formation sur tel ou tel thème ? et vous pensez pouvoir animer une telle formation ? eh bien, lancez-vous ! on vous monte ce cours sur mesure ! Quand le prix a été lancé, on m’a demandé de faire partie du jury, et j’ai tout de suite accepté. Il s’agissait, pour des personnes sans aucune expérience dans le domaine de la formation (continue) de proposer un plan pour une journée de 7 heures, sur un sujet, avec un fil rouge, de préparer quelques slides. On en a sélectionné certains, et on les a convié pour un oral, en fin d’année… et la remise des prix aura lieu cette semaine. Les vainqueurs pourront se jeter à l’eau dans quelques mois, et monter leur formation, leur cours, pour de vrai. Je trouve l’idée merveilleuse, et j’espère qu’on aura d’autres candidatures intéressantes pour la prochaine édition ! A suivre donc..

Devoir, Méthodes de Prévisions (consommation électrique)

Pour le projet de l’UE de prévision, je demande une prévision de la consommation d’électricité. Pour construire les modèles, je donne accès à un gros fichier qui donne, pour toute la France, la consommation d’électricité, entre 2012 et 2014, heure par heure.

  • comme en cours, la première étape sera de proposer une série de 52 prévisions pour chaque semaine de 2015,
  • la seconde étape sera de proposer deux séries de 31 prévisions, pour chaque jour d’août 2015 et de décembre 2015
  • enfin, la dernière étape sera une prévision de 168 prévisions, pour chaque heure de la semaine du 19 au 25 octobre 2015

Pour aider (?), je mets un lien vers la température en France  dans quelques villes, paris, toulouse, grenoble, nantes, dunkerque (via ftp.ncdc.noaa.gov). Il sera possible d’utiliser la vraie température de 2015 pour faire les prévisions de consommation électrique.  Idéalement, je veux un fichier html (markdown) et 4 fichiers csv avec les 4 séries de prévision.

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…

Using convolutions (S3) vs distributions (S4)

Usually, to illustrate the difference between S3 and S4 classes in R, I mention glm (from base) and vglm (from VGAM) that provide similar outputs, but one is based on S3 codes, while the second one is based on S4 codes. Another way to illustrate is to manipulate distributions.

Consider the case where we want to sum (independent) random variables. For instance two lognormal distribution. Let us try to compute the median of the sum.

The distribution function of the sum of two independent (positive) random variables is F_{S_2}(x)=\int_0^x F_{X_1}(x-y)dF_{X_2}(x)

pSum2 = function(x) integrate(function(y) 
plnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value

Let us visualize that cumulative distribution function

vx=seq(0.1,50,by=.1)
vy=Vectorize(pSum2)(vx)
plot(vx,vy,type="l",ylim=c(0,1))
abline(h=.5,lty=2)

Let us find an upper bound to compute (in a decent time) quantiles

pSum2(350)
[1] 0.99195

and then use the uniroot function to inverse that function

qSum = function(u) uniroot(function(x) 
Vectorize(pSum2)(x)-u, interval=c(0,350))$root
vu=seq(.01,.99,by=.01)
vv=Vectorize(qSum)(vu)

The median is here

qSum(.5)
[1] 14.155

Why not consider the sum of three (independent) distributions ? Its cumulative distribution function can be writen using our previous function F_{S_3}(x)=\int_0^x F_{S_2}(x-y)dF_{X_3}(x)

pSum3 = function(x) integrate(function(y) 
pSum2(x-y)*dlnorm(y,2,2),0,x)$value

If we look at some values we good

pSum3(4)
[1] 0.015624
pSum3(5)
Error in integrate(function(y) plnorm(x - y, 1, 2) * 
dlnorm(y, 2, 1),  : 
  maximum number of subdivisions reached

So obviously, there are computational issues here.

Let us consider the following alternative expression F_{S_3}(x)=\int_0^x F_{X_3}(x-y)dF_{S_2}(x). Of course, it is necessary here to compute the density of the sum of two variables

dSum2 = function(x) integrate(function(y) 
dlnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value
pSum3 = function(x) integrate(function(y) 
dlnorm(x-y,2,2)*dSum2(y),0,x)$value

Again, let us compute some values

pSum3(4)
[1] 0.0090285
pSum3(5)
[1] 0.01186

This one seems to work quite well. But it is just an illusion.

pSum3(9)
Error in integrate(function(y) dlnorm(x - y, 1, 2) *
 dlnorm(y, 2, 1),  : 
  maximum number of subdivisions reached

Clearly, with those S3-type functions, it wlll be complicated to run computations with 3 variables, or more.

Let us consider distributions in the S4-type format of the following package

library(distr)
X1 = Lnorm(mean=1,sd=2)
X2 = Lnorm(mean=2,sd=1)
S2 = X1+X2

To compute the median, we simply have to use

distr::q(S2)(.5)
[1] 14.719

We can also visualize it easily

plot(q(S2))

which looks (very) close to what we got, manually.  But here, it is also possible to work with the sum of 3 (independent) random variables

X3 = Lnorm(mean=2,sd=2)
S3 = X1+X2+X3

To compute the median, use

distr::q(S3)(.5)
[1] 33.208

The function is here

plot(q(S3))

Consommation et production électrique

Une version temporaire des slides est en ligne. Pour le TP, on continuera à utiliser les données suivante

url="http://freakonometrics.free.fr/elec.RData"
download.file(url,"elec.RData")
load("elec.RData")
k=620
n=nrow(electricite)
passe=1:k
futur=(k+1):n
y=electricite$Load[passe]
plot(y,type="l")

La semaine passée, nous avions commencé de modéliser la série de température (cette dernière ayant un impact important comme déterminant de la consommation électrique)

y=electricite$Temp[passe]
plot(y,type="l")
Y=ts(y,start=1996,frequency=52)
library(forecast)

plot(y[(53:k)-52],y[53:k])
abline(lm(y[53:k]~y[(53:k)-52]),col="red")
abline(a=0,b=1,lty=2)

La série temporelle est relativement autocorrélée, à l’ordre 52

acf(y,lag=120)

model1=auto.arima(Y)
acf(residuals(model1),120)

On va garder ce modèle en mémoire, et regarder par la suite ses prévisions. Tentons ici un SARIMA, à la main

model2=arima(Y,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
summary(model2)
model2=arima(Y,order = c(0,0,0), 
seasonal = list(order = c(2,0,0)))
summary(model2)
acf(residuals(model2),120)
model2=arima(Y,order = c(2,0,0), 
seasonal = list(order = c(1,0,0)))            
model2=arima(Y,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
summary(model2)
acf(residuals(model2),120)

Tentons ensuite avec une racine unité saisonnière

Z=diff(Y,52)
model3=arima(Z,order = c(0,0,1), 
seasonal = list(order = c(0,0,1)))
summary(model3)
model3=arima(Z,order = c(0,0,2), 
seasonal = list(order = c(0,0,1)))
summary(model3)
acf(residuals(model3),120)
model3=arima(Y,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
summary(model3)

On peut ensuite tenter un modèle à la Buys-Ballot

reg=lm(Temp~0+as.factor(NumWeek),
data=electricite[passe,])
cycle=predict(reg,newdata=electricite[passe,])
plot(y,type="l")
lines(cycle,col="red")
E=residuals(reg)
plot(E,type="l")
acf(E,120)
model4=arima(E,order = c(1,0,0))
summary(model4)

horizon=n-k
prev4E=predict(model4,n.ahead = horizon)$pred
prev4C=predict(reg,newdata=electricite[futur,])
prev4 =prev4E+prev4C
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(prev4,col="blue")

prev1=predict(model1,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev1,col="purple")

prev2=predict(model2,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev2,col="green")

prev3=predict(model3,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev3,col="orange")

MODELES = data.frame(y1=prev1,
y2=prev2,y3=prev3,y4=prev4, 
y =electricite$Temp[futur])
plot(MODELES$y,type="l")
lines(1:111,MODELES$y1,col="purple")
lines(1:111,MODELES$y2,col="green")
lines(1:111,MODELES$y3,col="orange")
lines(1:111,MODELES$y4,col="blue")

 

Les modèles prédictifs peuvent-il être justes ?

Dans Nosedive (traduit par le titre Chute Libre en France), le premier épisode de la saison 3 de la série télévisée Black Mirror, on découvre la dystopie d’une une société régie par une « cote personnelle », une note, un score allant de 0 à 5. Dans ce monde, chaque personne note les autres, les mieux notés ayant accès à de meilleurs services (priorité dans les services, meilleurs taux, meilleurs prix, etc). Cette tendance à construire des scores dans toutes sortes de domaines (historiquement sur les crédits mais aujourd’hui sur des aspects criminels, voire civiques dans certains pays) ne va-t-elle pas déboucher sur un monde qui serait un concours de popularité sans fin ? Et comment serait-elle conciliable avec une justice sociale, a priori souhaitable ?

Les scores de crédit et les réseaux sociaux

Un score de crédit est, d’un point de vue actuariel, une grandeur proportionnelle à la probabilité de ne pas honorer ses engagements en tant que créditeur. Ça peut être de plus pouvoir payer les échéances trois mois consécutifs, ou juste d’avoir un retard. Dans la vraie vie, comme toujours, c’est un peu plus compliqué. Aux Etats-Unis ou en Grande Bretagne, il n’est pas rare que les étudiants s’endettent sur des dizaines d’années pour avoir l’opportunité de suivre les cours qui les intéressent (même si la motivation est surtout d’obtenir un diplôme en fin de parcours). Mais surtout, dès qu’ils atteignent d’âge de 18 ans, des sociétés de notation de crédit vont surveiller tous leurs déplacements. Souvent à leur insu. Et si un jour, un crédit consommation ou hypothécaire est refusé, les raisons ne sont jamais motivées. Est-ce un retard dans un paiement de loyer ? Des amendes de bibliothèques oubliées ? Une facture d’eau impayée, vieille de plusieurs années ?

Les sociétés de notation de crédit aux États-Unis, mais aussi en Chine, commencent à étudier la possibilité d’utiliser des données provenant de médias sociaux pour améliorer le score de crédit. Compter le nombre de fois où un utilisateur utilise le mot « gaspillé » (« wasted » en anglais) dans ce qu’il poste en ligne ne peut-il pas révéler une information quant au remboursement de dettes ? C’est en tout cas ce que prétend l’analyste de crédit américain FICO: « If you look at how many times a person says ‘wasted’ in their profile, it has some value in predicting whether they’re going to repay their debt (…) It’s not much, but it’s more than zero » (cité dans McLannahan (2015)). En Chine, le prêteur peer-to-peer Jubao a révélé qu’il était plus susceptible de donner des « bonus » aux emprunteurs s’ils étaient des amis Facebook avec des célébrités, tel que le raconte Botsman (2017).

Pour l’instant, les sociétés de notation de crédit utilisent encore les données qu’elles connaissent bien (factures de services publics et cartes de crédit), mais elles imaginent que bien des informations intéressantes doivent être accessibles (d’une manière ou une autre) sur les réseaux sociaux. Mais les données sont encore rares, et difficiles à analyser. Quid de la composante sarcastique ou humoristique dans un tweet utilisant le mot « wasted » ? Comme souvent, la difficulté est que les données réellement pertinentes sont difficiles à obtenir. S’il est possible d’avoir des informations sur le payement du loyer quand un locataire passe par une agence, que faire pour les transactions entre deux particuliers ? Et si c’était possible, comment traiter le cas de colocataires ? Ne pas obtenir de crédit parce qu’un ancien colocataire n’a pas payé dans les temps devient dérangeant. D’autant plus s’il s’agit peut-être d’une facture de téléphone cellulaire réclamée abusivement par la compagnie de téléphonie, alors que l’abonnement avait été résilié.

Mais le gros « malus » dans le score de crédit est bien souvent le fait de ne jamais avoir eu de carte de crédit. On pourrait penser qu’une personne qui n’a pas eu besoin d’une carte de crédit (et se contentait d’une carte de débit, permettant d’acheter chez un commerçant, comme la majorité des cartes bancaires en France) est le propre d’une personne prudente, qui n’a pas besoin de crédit pour des dépenses quotidiennes. Mais pour les établissements de crédit, cette personne n’est pas fiable car on ne la connait pas. Et c’est à elle de prouver qu’elle l’est (on revient à la pratique récurrente d’inversion de la charge de la preuve évoqué dans Charpentier(2016)). C’est étrangement ce qui se passe aujourd’hui quand on veut entrer sur le sol américain sans avoir de page Facebook.

Dans un monde de surveillance généralisée

Et si les établissements de crédit n’étaient pas les seuls intéressés par notre vie ? Que serait un monde si, en plus de savoir si je paye mes factures à temps, certains cherchaient à connaître mes réseaux d’amis, à savoir quels journaux je lis, si je préfère acheter du lait entier ou du lait demi-écrémé ? Quand on visite le musée de la Stasi à Berlin, on découvre ce monde a existé, que 1 personne sur 63 était agent (ou indicateur) de la Stasi (en comptant les indicateurs occasionnels, la proportion peut atteindre une personne sur 6). Le musée décrit un panoptisme total, chacun étant observé en permanence, comme le décrit Foucault (1975). Mais ce cauchemar ne correspond-il pas à monde actuel, de surveillance permanente, plus ou moins consentie. Surveillance via les téléphones cellulaires (géolocalisation pour la fonction la plus courante, mais parfois aussi des enregistrements audio à l’insu de l’utilisateur par certaines applications), via les objets connectés, mais aussi les caméras de surveillance couplées à des algorithmes de reconnaissance faciale de plus en plus performant. Fin 2017, en Chine 170 millions de caméra étaient installées, et le cap des 300 millions devrait être atteint d’ici 2020. Lors d’une expérience tentée par la BBC[1], il a fallu 7 minutes pour retrouver le journaliste John Sudworth qui se promenait dans les rues.

Le danger est que l’on ne sait jamais trop qui contrôle. De plus en plus de compagnies privées de sécurité se sont associées aux gouvernements. Les fournisseurs de messageries électroniques lisent nos messages pour détecter du spam, mais aussi pour revendre certaines informations. Par exemple, dans les règles de confidentialité annexées aux conditions générales d’utilisation de Gmail (Google) on lit « Nos systèmes automatisés analysent vos contenus (y compris les e-mails) afin de vous proposer des fonctionnalités personnalisées sur les produits, telles que (…) des publicités sur mesure ». Les assureurs envisagent de plus de plus l’installation de boitiers GPS dans les voitures, mais en passant par des prestataires externes. Au-delà de la propriété des données (évoqué dans Charpentier & Suire (2016)) on peut s’interroger sur leur revente, et leur utilisation. Savoir que quelqu’un se rend régulièrement dans un centre de transfusion sanguine est une information potentiellement intéressante, surtout couplée à d’autres.

Depuis 2014, le gouvernement chinois travaille sur un système d’évaluation de ses propres citoyens, programmé pour être mis en place en 2020, comme le raconte Trujillo (2017). Ce « système de crédit social » vise à créer un « score citoyen » (pour reprendre l’expression de Galeon & Bergan (2017)), afin de prédire, et prévenir, les dangers potentiels, normalisant les conduites individuelles par des dispositifs panoptiques (par exemple la vidésurveillance), induisant des réflexes d’autodéfense et d’autocontrôle.  Comme le disait Foucault (1975), il s’agit de « faire que la surveillance soit permanente dans ses effets, même si elle est discontinue dans son action; que la perfection du pouvoir tende à rendre inutile l’actualité de son exercice » (même souvent, aujourd’hui, elle est en plus continue dans son action). Certains de ces scores sont utilisés par la police pour savoir où patrouiller pour faire baisser la criminalité, comme PredPol. Mais quand on y regarde de plus près, les prédictions disent, en substance, que les crimes auront lieu (en majorité) dans les zones (historiquement) les plus criminogènes de la ville. La frontière entre la banalité et la tautologie est étroite. Et le réel danger est que bien souvent, les scores transforment les probabilités en quasi-certitudes, et le soupçon devient une preuve, comme le notait Supiot (2015).

Justice prédictive et méthodes actuarielles

En juin 2010, un rapport de l’Académie de médecine préconisait d”« améliorer la pratique des expertises de dangerosité des criminels sexuels en enseignant et en diffusant les méthodes actuarielles ». Ces « méthodes actuarielles » sont tout simplement les techniques de scoring, de « profilage » tel que le défini le règlement européen relatif aux données personnelles du 27 avril 2016 (RGDP). Angèle Christin s’est intéressée aux algorithmes qui estiment la probabilité de récidive dans la justice pénale américaine. Comme elle l’a montré, ces techniques posent nombre de questions, en particulier des biais discriminatoires, l’opacité qui rend difficile les recours, et surtout la difficulté de comprendre ce qui est réellement calculé. Dans l’État de Virginie, un score entre 1 et 10 est renvoyé, convention reprise par Compas (Correctional Offender Management Profiling Alternative Sanctions ) qui offre en plus un code couleur qui prédit le risque de récidive violent. Il s’agit alors d’un outil d’aide à la décision, une machine ne pouvant placer une personne en détention seule (Christin et al. (2015)).

Les conclusions d’un score prédictif dépendent de deux éléments clés : le modèle utilisé, et les données. Dans la majorité des cas aux Etats Unis, les codes des modèles restent opaques (et donc impossible à attaquer), et rares sont ceux qui ont vu les données utilisées pour calibrer ces modèles. Mais on peut se demander si les décisions de justice ne sont pas eux aussi relativement opaques ? Certes, les juges doivent motiver leurs décisions, ce qui les rend critiquables et attaquables, mais si le processus était si transparent, les issues d’un procès (humains) ne devraient-elles pas être alors davantage prévisibles ? Enfin, les différents biais sont assez simples à comprendre. Supposons qu’être riche permet d’avoir un bon avocat, et avoir un bon avocat permet de ne pas avoir certaines condamnations. Dans ce cas, une variable liée à la richesse (le type de véhicule possédé par exemple) sera liée positivement avec le fait de ne pas être coupable (reconnu coupagle), et fera baisser le score de dangerosité. L’autre danger dans les biais de sélection est qu’ils sont parfois complexes à comprendre, voire paradoxaux. Un exemple classique est celui illustré sur la Figure 1. Pendant la seconde guerre mondiale, il a été demandé à des ingénieurs et des statisticiens comment renforcer les bombardiers qui faisaient face au feu ennemi.

Figure1 : Endroits endommagés des avions revenus (source : McGeddon 2016)

Le statisticien Abraham Wald a commencé à collecter des données sur les impacts dans la carlingue, comme le raconte Mangel & Samaniego (1984). A la surprise générale, il a recommandé de blinder les endroits des appareils qui présentaient le moins de dommages. En effet, les avions utilisés dans l’échantillon présentait un biais important : seuls les avions revenus ont été pris en compte. S’ils ont pu revenir avec des trous au bout des ailes, c’est que ces parties sont suffisamment solides. Et comme aucun avion n’est revenu avec des trous au niveau des moteurs des hélices, ce sont ces parties qu’il convenait de renforcer.

Un autre danger est celui où les relations causales sont inversées. Que penser de ce médecin qui prescrit un puissant neuroleptique à un patient mis en examen, de peur que la justice lui reproche de ne pas avoir vu la dangerosité de son patient, et qu’inversement, la justice s’appuie sur cette prescription pour prouver que le patient est dangereux ? Un algorithme mal conçu pourrait comprendre de travers le sens des relations causales.

Mais les modèles prédictifs en matière judiciaire ne sont pas que du côté des juges. Lors d’un accident corporel sur la route, la loi Badinter  (du 5 juillet 1985) prévoit un « droit à indemnisation » pour toute victime d’un accident de la circulation dans lequel est impliqué un véhicule terrestre à moteur. Lorsque la compagnie d’assurance du conducteur propose une indemnité, la victime fait une rapide analyse coût/bénéfice pour savoir s’il va au tribunal. Si elle ne construit pas formellement un modèle prédictif, elle tente de voir, à partir de quelques éléments à sa connaissance, les coûts de demander à un juge de trancher sur le montant de l’indemnité, et ses bénéfices (potentiels).

Autre point important, les juristes appellent ces modèles « prédictifs » des modèles « actuariels ». Or la première fonction des actuaires était d’actualiser, de calculer la valeur du temps. Et le temps judiciaire a des conséquences souvent désastreuses. En quoi une décision humaine, imparfaite, prise au bout de 5 ans de procédure serait « meilleure » qu’une décision automatique prise en 15 jours par une machine ? Nombre de personnes qui ont connu des procédures de plusieurs années, aboutissant à un non-lieu, rêvent de procédures accélérées. Car le « temps perdu » a une valeur, les actuaires le savent bien.

Que penser alors de cette efficacité des modèles algorithmiques ? La justice se doit d’être efficace, mais cette contrainte ne doit pas faire oublier l’objectif central qui est celui de rendre la justice. Que se passe-t-il si l’efficacité devient un objectif, pour ne pas dire le principal objectif ? Car c’est bien la question que posent les modèles prédictifs : quel est l’objectif que l’on cherche à maximiser ? Et comment le formule-t-on de manière simple ?

Aide à la décision, ou justification d’une prise de décision

Aux Etats-Unis, nombre de juges se sont vu reproché de motiver un jugement à l’aide d’outils d’aide à la décision, ce qui laisse planer un doute sur la fonction réelle de ces outils. L’idée était initialement d’apporter une aide. Récemment, plusieurs systèmes mis en place dans les années passées ont été questionnés. Par exemple en Australie, le STMP (Suspect Targeting Management Plan) proposait d’identifier si des préadolescents doivent être surveillés, ou pas. Ce modèle ressemble à s’y méprendre à n’importe quel modèle actuariel, c’est à dire un outil d’évaluation et prédiction des risques, en se concentrant soit sur les récidivistes, soit sur les personnes suspectées de commettre un futur crime. Or un rapport récent montrait que son utilisation n’a eu « aucun impact observable sur la prévention du crime »[2]. Parallèlement aux Etats-Unis, l’outil Compas (Correctional Offender Management Profiling Alternative Sanctions) a été critiqué dans Dressel & Farid (2017) : « Les défenseurs de ces systèmes soutiennent que les données et l’apprentissage automatique avancé rendent ces analyses plus précises et moins biaisées que celles des humains. Cependant, nous montrons que le logiciel d’évaluation des risques Compas, largement utilisé, n’est pas plus précis ou juste que des prédictions faites par des personnes qui ont peu ou pas du tout d’expertise en matière de justice pénale ». En questionnant des personnes recrutées sur Internet, sans compétences en droit, il s’agissait de prévoir si des personnes allaient, ou pas commettre un autre crime dans les deux ans à venir. Compas s’est trompé dans 34.8% des cas, et les internautes dans 33% des cas. Cela dit, on peut se demander ici ce que “se tromper” signifie. En l’occurrence, on ne mesure pas ici la récidive des personnes, mais la condamnation pour récidive des personnes. Et si les modèles (ou les gens) ne s’étaient pas trompés, mais que les juges, en revanche, oui ?

Prévoir, et se tromper

Et si un des soucis ne venait pas dans ce que l’on demande à un outil prédictif ? Prévoir, c’est (fondamentalement) établir une probabilité pour un fait futur. Comme cela avait été rappelé dans un débat sur les sondages et les élections, est-ce qu’on peut dire qu’on se trompe si on annonce qu’un événement peut se produire avec 5% de chances, et qu’il se produit effectivement ? Pour savoir si une technique de prévision est bonne, il faut collecter un ensemble de prévision, et les comparer aux observations. C’est ce que font les météorologues depuis une quinzaine d’années, et qui a été formalisé par Gneiting et al. (2007). Leur idée est qu’un ensemble de distributions prédictives est obtenu par un modèle \{\hat F_t,\hat F_{t+1},..,\hat F_{t+h}\} et il convient de comparer ces distributions aux observations \{y_t,y_{t+1},..,y_{t+h}\} – et non pas \{\hat y_t,\hat y_{t+1},..,\hat y_{t+h}\}. Il faut alors introduire une distance entre les distributions prédictives, et les observations. Dans un système physique, il est possible d’imaginer comprendre les différentes relations causales, et donc de prévoir. Mais dans les relations humaines (et la justice en est un exemple parfait), rien n’est aussi simple, aussi automatique que les lois de mécanique des fluides qui permettent de modéliser des phénomènes météorologiques.

Références

Binet, Jacques-Louis, 2010, La prévention médicale de la récidive chez les délinquants sexuels. Académie de Médecine.

Botsman, Rachel. 2017. Who Can You Trust?: How Technology Brought Us Together – and Why It Could Drive Us Apart. Portfolio Penguin

Charpentier, Arthur & Suire, Raphaël 2016. Données et santé: valeurs, acteurs et santé. Risques, 107

Charpentier, Arthur. 2016. Les dérives du principe de précaution. Risques. 108

Christin, Agnèle, Rosenblat, Alex & Boyd, Danah 2015. Courts and Predictive Algorithms. Datacivilrights

Dressel, Julie & Farid, Hany 2018. The accuracy, fairness, and limits of predicting recidivism. Science Adavances

Foucault, Michel 1975 Surveiller et punir, naissance de la prison. Gallimard

Galeon, Dom & Bergan, Brad 2017. China’s “Social Credit System” Will Rate How Valuable You Are as a Human. Futurism

Gneiting, Balabdaoui & Raftery 2007. Probabilistic forecasts, calibration and sharpness. JRRS-B, 69, 243–268.

Mangel, Marc & Samaniego, Francisco 1984.  Abraham Wald’s work on aircraft survivability , Journal of the American Statistical Association, vol. 79, no 386,‎ 259–267

McLannahan, 2015 Being ‘wasted’ on Facebook may damage your credit score (Octobre 2015, Financial Times)

Supiot, Alain 2015. La gouvernance par les nombres : cours au Collège de France, 2012-2014. Fayard.

Trujillo, Elsa 2017. La Chine met en place un système de notation de ses citoyens pour 2020. Le Figaro, décembre 2017

[1] Décrit dans « In Your Face: China’s all-seeing state » http://www.bbc.com/news/av/world-asia-china-42248056/in-your-face-china-s-all-seeing-state

[2] https://www.numerama.com/politique/300907-un-algorithme-teste-par-la-police-pour-anticiper-les-crimes-des-jeunes-inquiete-laustralie.html

(Advanced) R Crash Course, for Actuaries

The fourth year of the Data Science for Actuaries program started this morning. I will be there for the introduction to R. The slides are available online (created with slidify, the .Rmd file is also available)

A (standard) markdown is also available (as well as the .Rmd file). I have to thank Ewen for his help on slidify (especially for the online quizz, and the integration of leaflet maps or the rgl animated graph….)

Visualizing effects of a categorical explanatory variable in a regression

Recently, I’ve been working on two problems that might be related to semiotic issues in predictive modeling (i.e. instead of a standard regression table, how can we plot coefficient values in a regression model). To be more specific, I have a variable of interest Y that is observed for several individuals i, with explanatory variables \mathbf{x}_i, year t, in a specific region z_i\in\{A,B,C,D,E\}. Suppose that we have a simple (standard) linear model (forget about time here) y_i=\beta_0+\beta_1x_{1,i}+\cdots+\beta_kx_{k,i}+\sum_j \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i

Let us forget the temporal effect to focus on the spatial effect today. And consider some simulated dataset. There will be only one (continuous) explanatory variable. And I will generate correlated covariates, just to be more realistic.

n=1000
library(mnormt)
r=.5
Sigma=matrix(c(1,r,r,1), 2, 2)
set.seed(1)
X=rmnorm(n,c(0,0),Sigma)
X1=cut(X[,1],c(-100,quantile(X[,1],c(.1,.4,.7,.85)),
100),labels=LETTERS[1:5])
X2=X[,2]
Y=5+X[,1]-X[,2]+rnorm(n)/2
db=data.frame(Y,X1,X2)

Here we have y_i=\beta_0+\beta_1x_{1,i}+\sum_{j\in\{A,B,C,D,E\}} \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i The goal here is to get to graph to visualize the vector \hat\alpha=(\hat\alpha_A,\cdots,\hat\alpha_E). Let us run the linear regression

reg1=lm(Y~X1+X2,data=db)
idx=which(substr(names(reg1$coefficients), 1,2)=="X1")
v1=reg1$coefficients[idx]
names(v1)=LETTERS[2:5]
barplot(v1,col=rgb(0,0,1,.4))

Note that it is possible to add some sort of “confidence interval” to discuss significance (or to avoid to spend hours discussing differences in bar heights that are not significantly different)

library(Hmisc)
sv1=summary(reg1)$coefficients[idx,2]
(bp1=barplot(v1,ylim=range(c(0,v1+2*sv1))))
errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

My main concern here is the “reference” that is considered. Should A be the reference? Why not B

db$X1=relevel(db$X1,"B")
reg1=lm(Y~X1+X2,data=db)
idx=which(substr(names(reg1$coefficients),1,2)=="X1")
v1=reg1$coefficients[idx]
names(v1)=LETTERS[c(1,3:5)]
library(Hmisc)
sv1=summary(reg1)$coefficients[idx,2]
(bp1=barplot(v1)
errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

Why not the smallest one? Why not the largest one?… What if there is no simple way to choose. Furthermore, let us get back to the original point, which is that there might be some temporal aspects. More precisely, we can have \hat\alpha^{(t)}=(\hat\alpha_A^{(t)},\cdots,\hat\alpha_E^{(t)}). If we have also \hat\alpha^{(t+1)} and we get another plot, how do we interpret it. If for E the bar is taller, it means that relative to A, the difference has increased. I have the feeling that the interpretation is more complicated because we do not see, on that graph, changes in \hat\alpha^{(t)}_A.

Let us try something else. First, let us get back to the original setting

db$X1=relevel(db$X1,"A")

Consider here the regression without the intercept, so that all values remain

reg1=lm(Y~0+X1+X2,data=db)
idx=which(substr(names(reg1$coefficients),1,2)=="X1")
v1=reg2$coefficients[idx]
names(v1)=LETTERS[1:5]
barplot(v1)

It can be hard to read, especially if Y takes (very) large values, and you think that barplots should start at 0. But still, having those 5 values is nice. Why not rescale that graph?

A natural idea my be to consider the case where no spatial component is considered, and to look at the difference with that reference.

reg1=lm(Y~1+X2,data=db)
reg2=lm(Y~0+X1+X2,data=db)
idx=which(substr(names(reg2$coefficients),1,2)=="X1")
v1=reg2$coefficients[idx]
v2=v1-reg1$coefficients["(Intercept)"]
barplot(v2,col=rgb(0,0,1,.4))
sv2=summary(reg2)$coefficients[idx,2]
(bp2=barplot(v2,ylim=range(c(v2-2*sv2,v2+2*sv2))))
errbar(bp2[,1],v2,v2-2*sv2,v2+2*sv2,add=TRUE)

I like that graph, I should admit it. Now, I still have some remaining questions. For instance, can we insure that when only the intercept is considered, the value of \hat\beta_0 is somewhere between \hat\beta_A,\cdots,\hat\beta_E? Is it possible that \hat\beta_A-\hat\beta_0,\cdots,\hat\beta_E-\hat\beta_0 are all positive? In that case, I would find that hard to interpret.

Actually, if I really want values that can be seen as compared to some average, why not consider a (weighted) average of \hat\beta_A,\cdots,\hat\beta_E? (weights being here proportion in each class, in each region)

w=table(db$X1)
v3=v1-sum(w*v1)/sum(w)
(bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3))))
errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE)

I like that one. But what if, instead of normalizing at the end, we normalize the original dependent variable. By “normalize”, I mean “rescale”, to have a centered variable.

db$Y0=db$Y-mean(db$Y)
reg3=lm(Y0~0+X1+X2,data=db)
sv3=summary(reg3)$coefficients[idx,2]
(bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3))))
errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE)

This one is nice, because it is extremely simple to explain. But what if instead of a linear regression, we add a logistic one (with Y\in\{0,1\})? or a Poisson regression…

So maybe it cannot be the best solution here. Let us try something else… In insurance ratemaking, people like to use “zonier“. It is a two-stage regression. The idea is to run a regression without any spatial components, first. Then, consider the regression of residuals on spatial variables. Here, it would be something like

reg1=lm(Y~1+X2,data=db)
reg4=lm(residuals(reg1)~0+db$X1)

Since we focus on residuals, those are centered, and we have an easy interpretation of respective values

sv4=summary(reg4)$coefficients[idx,2]
v4=reg4$coefficients
(bp4=barplot(v4,names.arg=LETTERS[1:5])))
errbar(bp4[,1],v4,v4-2*sv4,v4+2*sv4,add=TRUE)

I guess that it can also be use in generalized linear models, with Pearson (or deviance) residuals.

Another possible idea can be the following. Again, the goal is not to have the true values, but to visualize on a graph how regions can be different. Here, all of them are significantly different. And in region A, Y is smaller, ceteris paribus (other things equal in the sense that we have taken into account x_1). And in region E it is larger. Here, the graph helps to “see” those differences.

Why not consider a completely different graph. What if we plot vector a instead of \alpha, where a_A can be interpreted as the value of the coefficient if we consider region A against “not region A“. What if we consider 5 regressions where dichotomous versions of Z are considered : Z_j=\mathbf{1}_{Z=j}.

v5=sv5=rep(NA,5)
names(v5)=LETTERS[1:5]
for(k in 1:5){
reg=lm(Y~I(X1==LETTERS[k])+X2,data=db)
v5[k]=reg$coefficients[2]
sv5[k]=summary(reg)$coefficients[2,2]}

We can plot that sequence of values, including some confidence intervals (that would be related to significance with respect to all other regions)

(bp5=barplot(v5,ylim=range(c(v5-2*sv5,v5+2*sv5))))
errbar(bp5[,1],v5,v5-2*sv5,v5+2*sv5,add=TRUE)

Looking at values does not give intuitive results, but I have the feeling that it is easy to explain what we plot (we compare each region to “the rest of the world”), and the ordering of a seems to be consistent with \alpha (but I could not prove it).

Here are some ideas I got. I should be able to provide other graphs, but I would love to discuss with anyone on that topics, to find a proper and nice way to visualize effects of a categorical explanatory variable in a regression model (that can be a logistic one). Comments are open…

Holt-Winters with a Quantile Loss Function

Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord & Snyder (2008)),

when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.

Furthermore, I like it because I think it has nice pedagogical features. Consider simple exponential smoothing, L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1} where \alpha\in(0,1) is the smoothing weight. It is locally constant, in the sense that {}_{t}\hat Y_{t+h} = L_{t}

 library(datasets)
 X=as.numeric(Nile)
 SimpleSmooth = function(a){
  T=length(X)
  L=rep(NA,T)
  L[1]=X[1]
  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
  return(L)
 }
 plot(X,type="b",cex=.6)
 lines(SimpleSmooth(.2),col="red")

When using the standard R function, we get

hw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1])
hw$alpha
[1] 0.2465579

Of course, one can replicate that optimal value

V=function(a){
     T=length(X)
     L=erreur=rep(NA,T)
     erreur[1]=0
     L[1]=X[1]
     for(t in 2:T){
         L[t]=a*X[t]+(1-a)*L[t-1]
         erreur[t]=X[t]-L[t-1] }
     return(sum(erreur^2))
}
optim(.5,V)$par
[1] 0.2464844

Here, the optimal value for \alpha is the one that minimizes the one-step prediction, for the \ell_2 loss function, i.e. \sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2 where here {}_{t-1}\hat Y_t = L_{t-1}. But one can consider another loss function, for instance the quantile loss function, \ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0}). The optimal coefficient is then obtained using

HWtau=function(tau){
loss=function(e) e*(tau-(e<=0)*1)
 V=function(a){
  T=length(X)
  L=erreur=rep(NA,T)
  erreur[1]=0
  L[1]=X[1]
  for(t in 2:T){
  L[t]=a*X[t]+(1-a)*L[t-1]
  erreur[t]=X[t]-L[t-1] }
 return(sum(loss(erreur)))
 }
 optim(.5,V)$par
}

Here is the evolution of \alpha^\star_\tau as a function of \tau (the level of the quantile considered).

T=(1:49)/50
HW=Vectorize(HWtau)(T)
plot(T,HW,type="l")
abline(h= hw$alpha,lty=2,col="red")

Note that the optimal \alpha is decreasing with \tau. I wonder how general this result can be…

Of course, one can consider more general exponential smoothing, for instance the double one, with L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}]andB_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}so that the prediction is now {}_{t}\hat Y_{t+h} = L_{t}+hB_t (it is now locally linear – and no longer constant).

hw=HoltWinters(X,gamma=FALSE,l.start=X[1])
hw$alpha
    alpha 
0.4200241 
hw$beta
      beta 
0.05973389

The code to compute the smoothed series is the following

DoubleSmooth = function(a,b){
  T=length(X)
  L=B=rep(NA,T)
  L[1]=X[1]; B[1]=0
  for(t in 2:T){
  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
 return(L+B)
 }

Here also it is possible to replicate R using the \ell_2 loss function

V=function(A){
     a=A[1]
     b=A[2]
     T=length(X)
     L=B=erreur=rep(NA,T)
     erreur[1]=0
     L[1]=X[1]; B[1]=X[2]-X[1]
     for(t in 2:T){
         L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
         B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] 
         erreur[t]=X[t]-(L[t-1]+B[t-1]) }
     return(sum(erreur^2))
}
optim(c(.5,.05),V)$par
[1] 0.41904510 0.05988304

(up to numerical optimization approximation, I guess). But here also, a quantile loss function can be considered

HWtau=function(tau){
loss=function(e) e*(tau-(e<=0)*1)
 V=function(A){
  a=A[1]
  b=A[2]
  T=length(X)
  L=B=erreur=rep(NA,T)
  erreur[1]=0
  L[1]=X[1]; B[1]=X[2]-X[1]
  for(t in 2:T){
   L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
   B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] 
   erreur[t]=X[t]-(L[t-1]+B[t-1]) }
  return(sum(loss(erreur)))
  }
     optim(c(.5,.05),V)$par
}

and we can plot those values on a graph

T=(1:49)/50
HW=Vectorize(HWtau)(T)
plot(HW[1,],HW[2,],type="l")
abline(v= hw$alpha,lwd=.4,lty=2,col="red")
abline(h= hw$beta,lwd=.4,lty=2,col="red")
points(hw$alpha,hw$beta,pch=19,col="red")

(with \alpha on the x-axis, and \beta on the y-axis). So here, it is extremely simple to change the loss function, but so far, it should be done manually. Of course, one do it also for the seasonal exponential smoothing model.