Freakonometrics

To content | To menu | To search

Wednesday, December 5 2012

Syntaxe pour les SARIMA

La forme générale (enfin, la plus simple) des modèles SARIMA est la suivante

http://latex.codecogs.com/gif.latex?(1-L)^d(1-L^s)^{d_s}\Phi(L)\Phi_s(L^s)%20(X_t-m)%20=%20\Theta(L)\Theta_s(L^s)\varepsilon_t

Afin de comprendre comment les écrire sous R, commençons par un processus simple, autorégressive, de la forme
http://latex.codecogs.com/gif.latex?(1-\phi_1%20L)X_t%20=%20\varepsilon_t

La syntaxe est ici

> arima(X, order = c(p=1, d=0, q=0))

Supposons que l'on souhaite rajouter une composante moyenne mobile, i.e.

http://latex.codecogs.com/gif.latex?(1-\phi_1%20L)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

La syntaxe devient

> arima(X, order = c(p=1, d=0, q=1))

Si on suppose maintenant que pour la composante autorégressive, on a une racine unité et le processus s'écrit

http://latex.codecogs.com/gif.latex?(1-%20L)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

L’estimation de ce modèle se fait avec la commande

> arima(X, order = c(p=0, d=1, q=1))

Si pour finir, on veut calibrer un modèle de la forme

http://latex.codecogs.com/gif.latex?(1-%20L)(1-\phi_1%20L-\phi_2%20L^2)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

on utilise la commande

> arima(X, order = c(p=2, d=1, q=1))

Supposons maintenant que notre série http://latex.codecogs.com/gif.latex?(X_t) ait été obtenu par différenciation saisonnière d'une série http://latex.codecogs.com/gif.latex?(Y_t), au sens où http://latex.codecogs.com/gif.latex?X_t=(1-L^{12})Y_t=\Delta^{12}%20Y_t, alors on devrait écrire

http://latex.codecogs.com/gif.latex?(1-%20L)(1-\phi_1%20L-\phi_2%20L^2)(1-L^{12})Y_t=%20(1-\theta_1%20L)\varepsilon_t

Pour estimer un tel modèle, la syntaxe est alors

> arima(X, order = c(p=2, d=1, q=1),
+ seasonal = list(order = c(0, 1, 0), period = 12))

En particulier, on peut avoir deux modèles assez proches pour modéliser une série cyclique: un bruit blanc avec une intégration saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L^{12})%20X_t%20=%20\varepsilon_t

dont la syntaxe serait

> arima(X, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

et un simple processus autorégressif à l'ordre 12. Mais là encore, soit je veux imposer que seule la dernière composante soit non nulle i.e.

http://latex.codecogs.com/gif.latex?(1-\phi_{12}%20L^{12})%20X_t%20=%20\varepsilon_t

ce qui s'écrit

> arima(X, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(1, 0, 0), period = 12))

soit je prends un processus plus général, de la forme

http://latex.codecogs.com/gif.latex?(1-\phi_1%20L-\cdots%20\phi_{12}%20L^{12})%20X_t%20=%20\varepsilon_t

dont la syntaxe devient

> arima(X, order = c(p=12, d=0, q=0))

En fait, on peut introduire un polynôme autorégressive avec uniquement des retards à l'ordre 12, en plus de la différentiation saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L)(1-L^{12})(1-\phi_1%20L)(1-\phi_{12}%20L^{12})%20X_t%20=%20(1-\theta_1%20L)\varepsilon_t

Ce modèle s'écrire alors sous la forme

> arima(X, order = c(p=1, d=1, q=1),
+ seasonal = list(order = c(1, 1, 0), period = 12))

On peut bien sûr aller plus loin, en autorisant non seulement une composante autorégressive saisonnière, mais pourquoi pas, une composante moyenne mobile saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L)(1-L^{12})(1-\phi_1%20L)(1-\phi_{12}%20L^{12})%20X_t%20=%20(1-\theta_1%20L)\varepsilon_t(1-\theta_{12}L^{12})\varepsilon_t

qui s'écrit, sous R,

> arima(X, order = c(p=1, d=1, q=1),
+ seasonal = list(order = c(1, 1, 1), period = 12))

On a vu la forme la plus générale des SARIMA. Enfin, comme je le disais en préambule, c'est aussi la plus simple. Car si je suppose qu'une autre cycle se superpose au cycle annuel, je peux le faire. Mais on va peut être essayer d'éviter de trop compliquer les notations...

Wednesday, November 14 2012

Série temporelles

Pour la section sur les séries temporelles, je mets des notes de cours en ligne. Les slides seront en ligne d'ici la semaine prochaine. Sinon, quelques séries dans le fichier ci-dessous

> source(
+ "http://freakonometrics.blog.free.fr/public/data/sourcets.R")

Les données sont les suivantes, avec l'utilisation du mot clé épilation sur Google,

> ts.plot(epilation,col="red")

..

ou l'utilisation du mot clé gym,

> ts.plot(gym,col="blue")

On a aussi la série des décès suite à des accidents de la route en Ontario,

> ts.plot(ontario,col="purple")

ou la série du nombre de vente de voitures au Québec

> ts.plot(quebec,col="purple")

La série suivante est la température minimale journalière au parc Montsouris, à Paris,

> ts.plot(temperature,col="green")

Pour finir, deux séries de transport, avec le nombre de voitures qui ont emprunté l'autoroute A7 en France

> ts.plot(autoroute7,col="red")

ou le nombre de voyageur dans un aéroport,

> ts.plot(airline,col="blue")

Tuesday, October 16 2012

R avec un mac

Petit compléments pour les étudiants qui utilisent mac. Le fonctionnement est (beaucoup) plus simple que sous windows. Sur le fonctionnement, je vois essentiellement deux différences. La première est sur la localisation,
> getwd()
[1] "/home/charpentier"
Il suffit d'utiliser des slashs, qui est naturel sous mac (comme sous linux d'ailleurs)
> setwd("/home/charpentier/code-r/")
La seconde est sur l'installation de packages, qui peut se faire en ligne de commandes,
> install.packages("lmtest",dependencies=TRUE)
A part ces deux points, je crois que tout le reste est identique... Sinon, il est possible (sous mac, windows et linux) d'utiliser RStudio. Comme RStudio n'est pas installé dans les salles informatiques, je ne l'utiliserais pas, mais le fonctionnement est relativement simple,

http://freakonometrics.free.fr/rstudio.gif

Monday, October 1 2012

Premier cours de régression

Mercredi aura lieu le premier cours ACT6420 méthodes de prévisions de cette nouvelle session. Le plan de cours est d'ailleurs en ligne. Le nombre d’étudiants inscrits ayant plus que doublé (plus d'une centaines d'étudiants inscrits), la forme changera probablement un peu. En particulier, je pense que je serais obligé d'utiliser des transparents (autant que faire se peu). Je mettrais en ligne les transparents du cours, au fur et à mesure. Sinon je continuerais à illustrer, en R, les différents aspects que nous aborderons.  D'autres informations seront en ligne sur la page ACT6420.

  • données de régression
La première partie du cours portera sur la régression sur des données individuelles.
Parmi les bases utilisées dans les premiers transparents, il y a des données sur la taille père-fils
> library(UsingR)
> data(father.son)
> father.son$fcm=father.son$fheight*2.54
> father.son$scm=father.son$sheight*2.54
> head(father.son)
fheight  sheight      fcm      scm
1 65.04851 59.77827 165.2232 151.8368
2 63.25094 63.21404 160.6574 160.5637
3 64.95532 63.34242 164.9865 160.8897
4 65.75250 62.79238 167.0113 159.4926
5 61.13723 64.28113 155.2886 163.2741
6 63.02254 64.24221 160.0773 163.1752
> attach(father,son)
> plot(fcm,scm,cex=.7)
> boxplot(father.son[,3:4],horizontal=TRUE)
Sinon j'ai aussi des données sur le prestige des métiers, au Canada,
> prestige = read.table("http://freakonometrics.free.fr/
prestige.txt",header=TRUE)
> head(prestige)
education income women prestige census type
GOV.ADMINISTRATORS      13.11  12351 11.16     68.8   1113 prof
GENERAL.MANAGERS        12.26  25879  4.02     69.1   1130 prof
ACCOUNTANTS             12.77   9271 15.70     63.4   1171 prof
PURCHASING.OFFICERS     11.42   8865  9.11     56.8   1175 prof
CHEMISTS                14.62   8403 11.68     73.5   2111 prof
PHYSICISTS              15.64  11030  5.13     77.6   2113 prof
La base que j'utiliserais probablement le plus, pour expliquer la régression simple sera une petite base contenant deux variables, une vitesse de véhicule au moment de freiner, et la distance de freinage.
> data(cars)
> head(cars)
speed dist
1     4    2
2     4   10
3     7    4
4     7   22
5     8   16
6     9   10
Par mes données préférées, il y a aussi cette base avec la taille et le poids d'individus,
> TP = read.table("http://freakonometrics.free.fr/
taillepoids.txt"
,header=TRUE,sep=";") > head(TP) indice sexe poids1 taille1 poids2 taille2 1 1 M 77 182 77 180 2 2 F 58 161 51 159 3 3 F 53 161 54 158 4 4 M 68 177 70 175 5 5 F 59 157 59 155 6 6 M 76 170 76 165
Pour aborder la régression multiple, on aura besoin de plus de variables. En particulier, on dipose de données par quartiers, à Chicago
> chicago=read.table("http://freakonometrics.free.fr/
chicago.txt",header=TRUE,sep=";")
> head(chicago)
Fire   X_1 X_2    X_3
1  6.2 0.604  29 11.744
2  9.5 0.765  44  9.323
3 10.5 0.735  36  9.948
4  7.7 0.669  37 10.656
5  8.6 0.814  53  9.730
6 34.1 0.526  68  8.231

X_1: proportion d'habitations anciennes
X_2: nombre de délits (ou de vols) commis
X_3: revenu médian dans le quartier
Fire
: nombre annuel d'incendies d'appartement (pour 1000 ménages)
Le but sera de modéliser le nombre d'incendies.
Enfin, on va aussi s'intéresser au nombre d'homicides aux États-Unis, à partir de la base
> US=read.table("http://freakonometrics.free.fr/US.txt",
+ header=TRUE,sep=";")

(je renvoie sur le précédant billet pour un descriptif précis). On pourra consulter aussi un ancien billet pour faire de la cartographie, et visualiser les différents états.

  • séries temporelles
Enfin, la seconde partie du cours portera sur la modélisation et la prévision de séries temporelles. Par exemple, on travaillera sur la fréquentation de l'autoroute A7, en France, entre Lyon et Marseille (appelé aussi autoroute du soleil)
> autoroute=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)
> A7
Jan   Feb   Mar   Apr   May   Jun   Jul   Aug
1989
1990 21421 25147 24264 41525 35118 42142 62099 65203
1991 22481 22801 29504 34225 44090 40767 61298 71567
1992 23203 24966 27008 40362 40985 40247 49724 72425
1993 25220 26773 27137 43147 42889 39830 66531 74559
1994 25063 27566 29909 42906 44010 41052 70949 72441
1995 25375 28523 29429 43042 41595 44800 69251 70671
1996 25323 28825 30899 42833 40943 44187 64291 74163
Sep   Oct   Nov   Dec
1989 40236 28874 23874 26545
1990 41830 28462 24957 26415
1991 43894 29799 26534 28489
1992 43086 31670 25571 30551
1993 43979 33113 26378 29959
1994 44447 33277 29176 21328
1995 43874 32676 27334 33004
1996 43275
> plot(A7,col="blue",lwd=2)
On travaillera aussi sur le nombre de ventes de voitures au Québec,

X=read.table( "http://freakonometrics.blog.free.fr/public/
data/car-sales-quebec.csv"
, header=TRUE,sep=";",nrows=108) Xt=ts(X[,2],start=c(1960,1),frequency=12)

Friday, September 21 2012

Transformation logarithmique de séries temporelles

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

plot(X,xlim=c(1,length(X)+20))

Mais on peut avoir l'intuition que modéliser le logarithme de la série pourrait être plus simple,

> X=log(Y)
> plot(X,xlim=c(1,length(X)+20))

On va alors tenter une modélisation par un processus ARMA de cette dernière série,

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

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

> temps=length(X)+1:24
> ciu=(E+2*sqrt(V))
> cil=(E-2*sqrt(V))
> polygon(c(temps,rev(temps)),c(ciu,
+ rev(cil)),col="yellow",border=NA)
> lines(temps,E,col="red",lwd=2)
> lines(temps,ciu,col="red",lty=2)
> lines(temps,cil,col="red",lty=2)

Maintenant, on va devoir remonter. On va utiliser un résultat que l'on a vu sur la transformation logarithmique dans une régression: si après avoir pris le logarithme, on a un modèle simple, Gaussien, c'est que le modèle initial était log-normal. On peut alors utiliser les propriétés de la loi lognormale, dont on connait les moments à partir de ceux de la loi Gaussienne sous-jacente. Pour la prévision, on n'a pas trop le choix,

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

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

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

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

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

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

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> polygon(c(temps,rev(temps)),c(ci1u,
> rev(ci1l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci1u,col="red",lty=2)
> lines(temps,ci1l,col="red",lty=2)

(qui est symétrique et centré sur notre prévision) et dans le second

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)
> polygon(c(temps,rev(temps)),c(ci2u,
> rev(ci2l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci2u,col="red",lty=2)
> lines(temps,ci2l,col="red",lty=2)

Thursday, September 6 2012

Simulation de séries temporelles

Un billet rapide pour reprendre le code tapé en cours, la semaine passée. Considérons  un processus autorégressif d'ordre 1,  où  est un bruit blanc, stationnaire, i.e.  appartient à l'intervalle . Le code pour simuler un tel processus est

n=1000
bruit=rnorm(n)
phi1= .85
X=rep(NA,n)
X[1]=0
for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]}
plot(acf(X),lwd=5,col='blue')
plot(pacf(X),lwd=5,col='blue')

ou avec un autocorrélation au premier ordre négative,

phi1= -0.7

On peut aussi regarder un processus autorégressif au second ordre, 

sur la figure ci-dessous (avec en haut à gauche le triangle de stationnarité du couple de paramètres).

phi1=  0.3
phi2=  0.5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=phi1*X[t-1]+phi2*X[t-2]+bruit[t]}

Histoire de changer un peu, on peut regarder un processus moyenne mobile au premier ordre,  où  est un paramètre dans .

theta1=  .8
X=rep(NA,n)
X[1]=0
for(t in 2:n){
X[t]=bruit[t]+theta1*bruit[t-1]}

ou une moyenne mobile du second ordre,

theta1= -.6
theta2=  .5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=bruit[t]+theta1*bruit[t-1]+
theta2*bruit[t-2]}


Inference and autoregressive processes

Consider a (stationary) autoregressive process, say of order 2,

for some white noise with variance . Here is a code to generate such a process,

> phi1=.5
> phi2=-.4
> sigma=1.5
> set.seed(1)
> n=240
> WN=rnorm(n,sd=sigma)
> Z=rep(NA,n)
> Z[1:2]=rnorm(2,0,1)
> for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]}

Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters.

  • using least square regression

A natural idea is to see here a regression model, and thus, if we consider a matrix formulation,

Here we can run (conditional) ordinary least squares estimation,

> base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)])
> regression=lm(Y~0+X1+X2,data=base)
> summary(regression)
 
Call:
lm(formula = Y ~ 0 + X1 + X2, data = base)
 
Residuals:
Min      1Q  Median      3Q     Max
-4.3491 -0.8890 -0.0762  0.9601  3.6105
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
X1  0.45107    0.05924   7.615 6.34e-13 ***
X2 -0.41454    0.05924  -6.998 2.67e-11 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.449 on 236 degrees of freedom
Multiple R-squared: 0.2561,	Adjusted R-squared: 0.2497
F-statistic: 40.61 on 2 and 236 DF,  p-value: 6.949e-16
 
> regression$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)$sigma
[1] 1.449276
  • using Yule-Walker equations

As we've seen in class, we can easily get the following equations for the autocovariance functions,

which can also be written

So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions

The code is the following

> rho1=cor(Z[1:(n-1)],Z[2:n])
> rho2=cor(Z[1:(n-2)],Z[3:n])
> A=matrix(c(1,rho1,rho1,1),2,2)
> b=matrix(c(rho1,rho2),2,1)
> (PHI=solve(A,b))
[,1]
[1,]  0.4517579
[2,] -0.4155920

Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation)

> estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2)
> sd(estWN)
[1] 1.445706

This estimator is probably not the best one (we can take into account that we've lost two degrees of freedom), but as a starting point, let us consider this one.

  • using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. The standard model is a Gaussian model, i.e.

In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A[1];  phi2=A[2]
+ sigma=A[3]	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
$par
[1]  0.4509685 -0.4144938  1.4430930
 
$value
[1] 425.0164
 
$counts
function gradient
88       NA
 
$convergence
[1] 0
 
$message
NULL
Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient

for the second one,

and finally for the standard deviation of the innovation process

All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).



Friday, June 29 2012

Provisionnement et modélisation individuelle des sinistres

Il y a quelques années, Guillaume Beneteau avait fait son mémoire d'actuariat, en France, sur l'utilisation de données individuelles pour faire du provisionnement (au lieu de travailler sur des donnees aggregées dans des triangles). Dans le même genre d'idées, Ngoc An Dinh et Gilles Chau viennent de soutenir leur mémoire cette semaine à l'ENSAE, encadré par Frédéric Planchet. L'idée était de modéliser les processus de règlements avec deux approches : la première intitulée dynamiques markoviennes qui considère les processus de règlements (avec une certaine dynamique sous-jacente), et une basée sur des modèles linéaires généralisés qui suppose que les règlements appartiennent à une famille de distribution paramétrique et les paramètres dépendent des facteurs liés aux règlements. Le mémoire est en ligne sur le blog, ainsi que le code R.

Saturday, June 23 2012

Actuarial models with R, Meielisalp

I will be giving a short course in Switzerland next week, at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The long version of the slides for the course on Actuarial models with R can be found online with the #Rmetrics tag, and the short version will be uploaded soon. There will be some practicals, based on Swiss mortality table for the part on demography. The datasets can be uploaded using the following code,

DEATH=read.table(
"http://freakonometrics.free.fr/DeathsSwitzerland.txt",
header=TRUE,skip=2)
EXPOSURE=read.table(
"http://freakonometrics.free.fr/ExposuresSwitzerland.txt",
header=TRUE,skip=2)
DEATH$Age=as.numeric(as.character(DEATH$Age))
DEATH=DEATH[-which(is.na(DEATH$Age)),]
EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age))
EXPOSURE=EXPOSURE[-which(is.na(EXPOSURE$Age)),]
  • based on those datasets, plot the log mortality rates for male and female,

  • for those two datasets, plot the log mortality rates in 1900 and 1950, respectively
  • for those two datasets, plot the log mortality rates for the cohort born on 1900 and 1950, respectively
  • on the total dataset (male and female), fit a Lee-Carter model. Plot the age coefficients

  • plot the time coefficients and propose a forecast for that series of estimators.

  • plot the residuals obtained from the regression

  • using those estimates, and the forecasts, project the log-mortality rates

  • extrapolate the survival function of an insured aged 40 in 2000, and compare it with the one obtained on the longitudinal dataset.

  • based on those survival functions, compute actuarial present values for several quantities, e.g. whole life annuities for some insured aged 40, and whole life insurances, and compare those values from 1900 till 2040 (on the graphs below, titles were inverted).

Then, we will briefly mention payment triangles. We will work on the triangle used on http://rworkingparty.wikidot.com/ that can be downloaded below,
OthLiabData=read.csv(
"http://www.casact.org/research/reserve_data/othliab_pos.csv",
header=TRUE, sep=",")
library(plyr)
OL=SumData=ddply(OthLiabData,.(AccidentYear,DevelopmentYear,
DevelopmentLag),summarise,IncurLoss=sum(IncurLoss_H1-BulkLoss_H1),
CumPaidLoss=sum(CumPaidLoss_H1), EarnedPremDIR=
sum(EarnedPremDIR_H1))
LossTri <- as.triangle(OL, origin="AccidentYear",
dev="DevelopmentLag", value="IncurLoss")
Year=as.triangle(OL, origin="AccidentYear",
dev="DevelopmentLag", value="DevelopmentYear")
TRIANGLE=LossTri
TRIANGLE[Year>1997]=NA
Here, the triangle looks like that
> TRIANGLE
dev
origin      1      2      3      4      5      6      7      8      9     10
1988 128747 195938 241180 283447 297402 308815 314126 317027 319135 319559
1989 135147 208767 270979 304488 330066 339871 344742 347800 353245     NA
1990 152400 238665 297495 348826 359413 364865 372436 372163     NA     NA
1991 151812 266245 357430 400405 423172 442329 460713     NA     NA     NA
1992 163737 269170 347469 381251 424810 451221     NA     NA     NA     NA
1993 187756 358573 431410 476674 504667     NA     NA     NA     NA     NA
1994 210590 351270 486947 581599     NA     NA     NA     NA     NA     NA
1995 213141 351363 444272     NA     NA     NA     NA     NA     NA     NA
1996 237162 378987     NA     NA     NA     NA     NA     NA     NA     NA
1997 220509     NA     NA     NA     NA     NA     NA     NA     NA     NA
  • suggest an estimation for the amount of reserves, all years.
  • using a Poisson regression, propose a VaR with level 99.5% for future payments, for all claims that already occurred.

Thursday, June 21 2012

A la recherche des groupes de trolls sur Twitter

Il y a quelques jours avait mis en ligne sur son blog http://olihb.com/ une très jolie carte permettant de visualiser les personnes (ou les comptes) impliquées sur Twitter sur les manifestations récentes au Québec (la carte était basée sur gephi). Beaucoup de monde twittait avec des hashtags, et l'idée était d'utiliser les plus classiques et les plus populaires pour identifier les comptes actifs (#ggi, #manifencours, etc). Mais la semaine passée j'avais l'impression que le débat s’était durci, avec beaucoup plus de trolls dans les discussions. 

Bref, le principe est que les trolls évoluent souvent en petits groupes, donc avec @3wen on a voulu voir si on pouvait identifier facilement les groupes de trolls (sans pour autant faire du sentiment analysis, juste visualiser des personnes qui se suivent beaucoup entre elles, les connivences entre twittos, i.e. du suivi mutuel).

La méthodologie est simple. On a deux bases à notre disposition. Un première contient tous les tweets contenant #ggi entre le 7 juin à 1 heure du matin et le 8 à minuit 45 (soit 14,700 tweets). La seconde contient tous les tweets contenant #ggi entre le 12 juin à 23 heures 30 et le 14 juin à 17 heures (soit 14,784 tweets). Ah oui, les heures sont GMT. Plus précisément, on a la fréquence d'arrivée suivante

Les durées d'observations sont courtes, mais l'idée est que si une personne a tweetté au moins une fois pendant ces périodes, on les a dans notre base, et ça suffit pour les suivre. On ne va pas étudier l'intensité de la mise en ligne de tweets. Je ne rentre pas dans le détail du code, mais c'est basé sur les fonctions présentées dans un précédant billet, et toujours avec la même librairie,

require("RJSONIO")

Sur des deux bases, i.e. 29,484 tweets, on récupère un ensemble de 6,556 twittos uniques (mais répartis dans deux listes: ceux qui ont tweeté le 7 et ceux qui ont tweeté le 14). Pour toutes ces personnes, on est allé voir qui les suivaient. Sur le principe, c'est simple, mais on fonctionne avec des API qui imposent des limites horaires. Bref, il faut bricoler. Pour une personne, on récupère avec la fonction suivante la liste de ses followers,

recuperefollowers=function(id,nbrequetes){
followers=try(scan(paste(
"http://api.twitter.com/1/followers/ids.json?cursor=-1&user_id=",
id,sep=""),what = "character", encoding="latin1"))
if(is.null(attr(followers,"condition"))){
followers=paste(followers[1:length(followers)],collapse=" ")
followers=fromJSON(followers, method = "C")
id_followers=followers$id
next_cursor=followers$next_cursor_str
nbrequetes=nbrequetes+1
plusDe5000=FALSE
while(nbrequetes<150 & next_cursor!="0" &
is.null(attr(followers,"condition"))){
plusDe5000=TRUE
followers=scan(paste(
"http://api.twitter.com/1/followers/ids.json?cursor=",
next_cursor,"&user_id=",id,sep="")
,what = "character", encoding="latin1")
followers=paste(followers[1:length(followers)],collapse=" ")
followers=fromJSON(followers, method = "C")
id_followers=c(id_followers,followers$id)
next_cursor=followers$next_cursor_str
nbrequetes=nbrequetes+1
}
if(plusDe5000 & nbrequetes>=149){
return(c(list(),nbrequetes,TRUE))
}else{
return(c(list(id_followers),nbrequetes,FALSE))
}}else{
nbrequetes=nbrequetes+1
print(paste("Probleme rencontre avec ",id,sep=""))
return(c(list("PROBLEME"),nbrequetes,FALSE))}}

On peut ensuite lancer cette fonction sur tous les identifiants qu'on a récupéré, et on stocke tout le monde dans une (petite) base

temp=try(recuperefollowers(
lesIDAParcourir[unID_index],nbrequetes))
lesFollowers=temp[[1]]
resul=cbind(rep(lesIDAParcourir[unID_index]
,length(lesFollowers)),lesFollowers)
nbrequetes=temp[[2]]
refaire=temp[[3]]
Sys.sleep(runif(1,1,1.5))
compte=compte+1
write.table(resul,paste("auto_resul_",
compte,".txt",sep=""),sep=";")

Si on veut que ça tourne, il faut juste faire patienter un heure après avoir fait 150 requêtes. Bref, dans la boucle, on met une petite fonction qui temporise,

if(nbrequetes>148){
nbrequetes=1
Sys.sleep(60*60+60*trunc(runif(1,2,3))) }

Une fois créées les 6,500 bases de nos 6,500 comptes, contenant la liste de tous les followers, on va les agréger (fort heureusement, on n'avait pas les  gros comptes avec plus d'un million de followers). La liste se récupère avec 

N=list.files("id/petits")

D'ailleurs, pour recoller les 6,500 bases ensemble, Ewen suggérait le code suivant, beaucoup beaucoup plus rapide qu'une gros boucle,

recupere_petits_fichiers=function(x){
temp=read.table(paste("id/petits/",x,sep=""),
header=TRUE,sep=";",comment.char="",
check.names=FALSE,colClasses=
c("character","character"))
colnames(temp)=list("id_twittos","id_follower")
head(temp)

On fait ensuite une récupération des followers pour la première base de tweets, et une autre pour ceux de la seconde. Les deux listes sont obtenues avec le code suivant

present_avant_temp=temp[which(
temp$id_twittos%in%lesID_avant),]
if(nrow(present_avant_temp)>0){
present_avant_temp=
present_avant_temp[which(
present_avant_temp$id_follower%in%lesID_avant),]}
else{
present_avant_temp=NULL}
 
present_apres_temp=
temp[which(temp$id_twittos%in%lesID_apres),]
if(nrow(present_apres_temp)>0){present_apres_temp=
present_apres_temp[
which(present_apres_temp$id_follower%in%lesID_apres),]}
else{present_apres_temp=NULL}
return(c(list(present_avant_temp),
list(present_apres_temp),
list(nbfollowers=nrow(temp))))

Pour constituer notre grosse base, la ruse est d'utiliser non pas une boucle sur tous les petites bases, mais

resul=lapply(N,recupere_petits_fichiers)
liens_avant=data.frame(do.call(
"rbind",lapply(resul,function(x) x[[1]])),
stringsAsFactors=FALSE)
liens_apres=data.frame(do.call(
"rbind",lapply(resul,function(x) x[[2]])),
stringsAsFactors=FALSE)

On a ainsi récupéré une liste avec plusieurs centaines de milliers de comptes Twitter. On a retenu seulement ceux qui avaient twitté avec le hashtag #ggi. On a ainsi pu constituer, pour nos 6,556 twittos une base de personnes qui tweettent, et de personnes qui les suivent (et qui ont twitté pendant cette - courte - période avec le hastag #ggi). Bref, une grosse matrice de correspondance de qui suit qui. On a pu le faire sur les deux bases. C'est précisément la différence entre les deux bases que l'on a cherché à visualiser, et gephi semble permettre de comparer les bases dans le temps. A gauche, on a le nuage de connexions entre les comptes pour la première période, et à droite pour la seconde. 

Bon, le premier est plus compact, mais on voit qu'il y a des différences, en particulier sur les bords, que ce soit en haut

ou à droite du nuage,

Dans les deux cas, on voit que plein de petits comptes sont venus se greffer au nuage, souvent en suivant quelqu'un qui le suit en retour. Et parfois, ce sont les seuls connexions qu'ils ont (qui utilisent le hashtag #ggi). Ce sont ces comptes qui pourraient constituer les nouveaux trolls qu'on espérait identifier. L'avantage est qu'on a les noms de tous ces comptes (on ne les a pas mis sur le graphique), mais dans un second temps, on va aller regarder ce qu'on fait précisément ces comptes. Peut-être interprétons nous mal ces graphiques, qui représentent une arrivée de nouveaux comptes, et que cette dynamique est la même pour tous les nouveaux comptes, qui commencent par de faibles connexions, et ensuite vont s'étoffer au fur et à mesure. Ce qu'il faut que l'on comprenne aussi c'est pourquoi ils sont tous aussi rapprochés. Ils semblent être suivis par des personnes qui sont suivis pas les mêmes personnes. Ce qui laisse à penser que l'on a effectivement identifié une certaine catégorie de comptes.
On peut reprendre les graphs avec une résolution plus fine, avant

et après

ou alors en changeant les options, avant

à comparer avec après

On voit peut être encore plus clairement sur ces deux dernières images les nouveaux comptes se greffer sur le graph existant. Si la qualité ne suffit pas, des graphiques (très) haute résolution sont téléchargeables (15Mo chacun), avant et après (avec les identifiants cette fois). Sinon, en promenant la souris sur le dessin dessous, on voit encore mieux la différence,

Maintenant - pour être complètement honnête - sur la lecture des graphs, j'avoue ne pas avoir tout bien compris. Les représentations semblent claires, mais j'aurais bien aimé être certain d'avoir compris la construction proposée par gephi avant de les mettre en ligne. Car si on regarde sur moins ce compte, on comprend que l'algorithme permet effectivement de bien mettre en avant les réseaux. Parmi les illustrations les plus frappantes, la figure suivante montre des interconnections avec Facebook à gauche, et purement aléatoire à droite (illustration trouvé sur le site de gephi),

La grande difficulté est alors de représenter spatialement les groupes, comme le notent Gastner et Newman dans The Spatial Structure of NetworksWatts et Strogatz dans Collective dynamics of 'small-world' networks ou encore Nisha et Venkatesh dans Small worlds: How and why. Mais pour etre tout à fait complet, je renverrais aussi vers la thèse de doctorat de Jure Leskovec, intitulée Dynamics of large networks ainsi qu'au livre (intégralement téléchargeable) Networks, Crowds, and Markets: Reasoning About a Highly Connected World de David Easley et Jon Kleinberg. Bref, on a découvert un outils probablement très riche, mais il va nous falloir du temps pour comprendre ce qu'il fait vraiment.

Thursday, June 7 2012

La pyramide des âges, au Canada

L'autre jour, Statistique Canada publiait une jolie animation interactive permettant de visualiser la déformation de la "pyramide" des âges (évoquée par exemple sur http://www.lapresse.ca/). Bon, j'ai des données assez proches, tirées de http://www.mortality.org/).

> pop=read.table(
+ "http://freakonometrics.blog.free.fr/PopulationCanada.txt",
+ header=TRUE,skip=2)
> pop$Age=as.numeric(as.character(pop$Age))
> pop$Year=as.numeric(as.character(pop$Year))
> pop=pop[is.na(pop$Age)==FALSE,]

On peut faire assez facilement des pyramides des ages, par année, avec le code suivant,

> library(plotrix)
> seuils=seq(0,110,by=10) > pop$tranche=cut(pop$Age,seuils, right = FALSE) > an=1955 > base=pop[pop$Year==an,] > women=aggregate(base$Female, + list(base$tranche),sum)[,2] > men=aggregate(base$Female, + list(base$tranche),sum)[,2] > nom=as.character(unique(pop$tranche)) > pyramid.plot(men/sum(men)*100, + women/sum(women)*100,labels=nom,gap=2, + lxcol=c("blue","blue","purple","purple","purple", + "purple","red","red","red","red","red"), + rxcol=c("blue","blue","purple","purple","purple", + "purple","red","red","red","red","red"))

http://freakonometrics.blog.free.fr/public/perso6/pyramide-ages-canada2.gif

Au lieu de faire des barres horizontales, on peut cumuler, par tranche d'age, et regarder l'évolution dans le temps des proportions "moins de 20 ans" (bleu), "entre 20 et 60 ans" (mauve) et "plus de 60 ans" (rouge).

> seuils=c(0,20,60,110)
> pop$tranche2=cut(pop$Age,seuils, right = FALSE)
> YEAR=1921:2008
> totaux=matrix(NA,length(YEAR),3)
> for(i in 1:length(YEAR)){
+ base=pop[pop$Year==YEAR[i],]
+ totaux[i,]=aggregate(base$Total,
+ list(base$tranche2),sum)[,2]
+ }
> stotaux=apply(totaux,1,sum)
> X=totaux[,1]/stotaux
> Y=X+totaux[,2]/stotaux

Si on regarde bien la pyramide, on s’aperçoit qu'après guerre, on a une génération importante, qui se déplace ensuite progressivement vers le haut, les baby-boomers. Pour visualiser encore davantage cette génération, on peut l'isoler sur le graphique ci-dessous, en suivant la cohorte née entre 1945 et 1950.

Mais on reviendra un peu plus en détails très bientôt sur cette génération (peut-être plutôt sur des données françaises cette fois afin de limiter les erreurs d'interprétation). A suivre...

Tuesday, May 15 2012

Finding Waldo, a flag on the moon and multiple choice tests, with R

I have to admit, first, that finding Waldo has been a difficult task. And I did not succeed. Neither could I correctly spot his shirt (because actually, it was what I was looking for). You know, that red-and-white striped shirt. I guess it should have been possible to look for Waldo's face (assuming that his face does not change) but I still have problems with size factor (and resolution issues too). The problem is not that simple. At the http://mlsp2009.conwiz.dk/ conference, a price was offered for writing an algorithm in Matlab. And one can even find Mathematica codes online. But most of the those algorithms are based on the idea that we look for similarities with Waldo's face, as described in problem 3 on http://www1.cs.columbia.edu/~blake/'s webpage. You can find papers on that problem, e.g. Friencly & Kwan (2009) (based on statistical techniques, but Waldo is here a pretext to discuss other issues actually), or more recently (but more complex) Garg et al. (2011) on matching people in images of crowds.

What about codes in R ? On http://stackoverflow.com/, some ideas can be found (and thank Robert Hijmans for his help on his package). So let us try here to do something, on our own. Consider the following picture,

With the following code (based on the following file) it is possible to import the picture, and to extract the colors (based on an RGB decomposition),
> library(raster)
> waldo=brick(system.file("DepartmentStoreW.grd",
+ package="raster"))
> waldo
class       : RasterBrick
dimensions  : 768, 1024, 786432, 3 (nrow,ncol,ncell,nlayer)
resolution  : 1, 1  (x, y)
extent      : 0, 1024, 0, 768  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values      : C:\R\win-library\raster\DepartmentStoreW.grd
min values  : 0 0 0
max values  : 255 255 255
My strategy is simple: try to spot areas with white and red stripes (horizontal stripes). Note that here, I ran the code on a Windows machine, the package is not working well on Mac. In order to get a better understanding of what could be done, let us start with something much more simple. Like the picture below, with Waldo (and Waldo only). Here, it is possible to extract the three colors (red, green and blue),
> plot(waldo,useRaster=FALSE)

     

It is possible to extract the red zones (already on the graph above, since red is a primary color), as well as the white ones (green zones on the graphs means a white region on the picture, on the left)
# white component
white = min(waldo[[1]] , waldo[[2]] , waldo[[3]])>220
focalswhite = focal(white, w=3, fun=mean)
plot(focalswhite,useRaster=FALSE)
 
# red component
red = (waldo[[1]]>150)&(max(  waldo[[2]] , waldo[[3]])<90)
focalsred = focal(red, w=3, fun=mean)
plot(focalsred,useRaster=FALSE)
i.e. here we have the graphs below, with the white regions, and the red ones,

From those two parts, it has been possible to extract the red-and-white stripes from the picture, i.e. some regions that were red above, and white below (or the reverse),
# striped component
striped = red; n=length(values(striped)); h=5
values(striped)=0
values(striped)[(h+1):(n-h)]=(values(red)[1:(n-2*h)]==
TRUE)&(values(red)[(2*h+1):n]==TRUE)
focalsstriped = focal(striped, w=3, fun=mean)
plot(focalsstriped,useRaster=FALSE)
So here, we can easily spot Waldo, i.e. the guy with the red-white stripes (with two different sets of thresholds for the RGB decomposition)

 Let us try somthing slightly more complicated, with a zoom on the large picture of the department store (since, to be honest, I know where Waldo is...).

Here again, we can spot the white part (on the left) and the red one (on the right), with some thresholds for the RGB decomposition

Note that we can try to be (much) more selective, playing with threshold. Here, it is not very convincing: I cannot clearly identify the region where Waldo might be (the two graphs below were obtained playing with thresholds)

 And if we look at the overall pictures, it is worst. Here are the white zones, and the red ones,

and again, playing with RGB thresholds, I cannot spot Waldo,

Maybe I was a bit optimistic, or ambitious. Let us try something more simple, like finding a flag on the moon. Consider the picture below on the left, and let us see if we can spot an American flag,

Again, on the left, let us identify white areas, and on the right, red ones

Then as before, let us look for horizontal stripes

Waouh, I did it ! That's one small step for man, one giant leap for R-coders ! Or least for me... So, why might it be interesting to identify areas on pictures ? I mean, I am not Chloe O'Brian, I don't have to spot flags in a crowd, neither Waldo, nor some terrorists (that might wear striped shirts). This might be fun if you want to give grades for your exams automatically. Consider the two following scans, the template, and a filled copy,

A first step is to identify regions where we expect to find some "red" part (I assume here that students have to use a red pencil). Let us start to check on the template and the filled form if we can identify red areas,
exam = stack("C:\\Users\\exam-blank.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)
plot(focalsred,useRaster=FALSE) 
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)
plot(focalsred,useRaster=FALSE)


First, we have to identify areas where students have to fill the blanks. So in the template, identify black boxes, and get the coordinates (here manually)
exam = stack("C:\\Users\\exam-blank.png")
black = max(  exam[[1]] ,exam[[2]] , exam[[3]])<50
focalsblack = focal(black, w=3, fun=mean)
plot(focalsblack,useRaster=FALSE)
correct=locator(20)
coordinates=locator(20)
X1=c(73,115,156,199,239)
X2=c(386,428.9,471,510,554)
Y=c(601,536,470,405,341,276,210,145,79,15)
LISTX=c(rep(X1,each=10),rep(X2,each=10))
LISTY=rep(Y,10)
points(LISTX,LISTY,pch=16,col="blue")

The blue points above are where we look for students' answers. Then, we have to define the vector of correct answers,
CORRECTX=c(X1[c(2,4,1,3,1,1,4,5,2,2)],
X2[c(2,3,4,2,1,1,1,2,5,5)])
CORRECTY=c(Y,Y)
points(CORRECTX, CORRECTY,pch=16,col="red",cex=1.3)
UNCORRECTX=c(X1[rep(1:5,10)[-(c(2,4,1,3,1,1,4,5,2,2)
+seq(0,length=10,by=5))]],
X2[rep(1:5,10)[-(c(2,3,4,2,1,1,1,2,5,5)
+seq(0,length=10,by=5))]])
UNCORRECTY=c(rep(Y,each=4),rep(Y,each=4))
Now, let us get back on red areas in the form filled by the student, identified earlier,
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=5, fun=mean)

Here, we simply have to compare what the student answered with areas where we expect to find some red in,
ind=which(values(focalsred)>.3)
yind=750-trunc(ind/610)
xind=ind-trunc(ind/610)*610
points(xind,yind,pch=19,cex=.4,col="blue")
points(CORRECTX, CORRECTY,pch=1,
col="red",cex=1.5,lwd=1.5)
Crosses on the graph on the right below are the answers identified as correct (here 13),
> icorrect=values(red)[(750-CORRECTY)*
+ 610+(CORRECTX)]
> points(CORRECTX[icorrect], CORRECTY[icorrect],
+ pch=4,col="black",cex=1.5,lwd=1.5)
> sum(icorrect)
[1] 13

In the case there are negative points for non-correct answer, we can count how many incorrect answers we had. Here 4.
> iuncorrect=values(red)[(750-UNCORRECTY)*610+
+ (UNCORRECTX)]
> sum(iuncorrect)
[1] 4
So I have not been able to find Waldo, but I least, that will probably save me hours next time I have to mark exams...

Notes de cours sur les séries temporelles

La session d'hiver n'étant pas terminée, je vais poster mes notes de cours sur la dernière section (sur la modélisation de séries temporelles) pour le cours ACT2040. Il s'agit - comme je l'avais dit en cours - d'une remise au goût du jour de notes tapées il y a une dizaine d'années. J'ai également rajouté du code R, mais il doit resté un certain nombre de coquilles et de fautes de frappe. Je profiterais des jours qui viennent pour réviser cette version.

Tuesday, April 24 2012

Régression médiane et géométrie

Suite à mon exposé d'hier, Pierre me demandait d'argumenter un point que j'avais évoqué oralement sans le justifier: " une régression médiane passe forcément par deux points du nuage ". Par exemple,

x=c(1,2,3)
y=c(3,7,8)
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
library(quantreg)
abline(rq(y~x,tau=.5),col="red")

Essayons de le justifier... de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe "entre" les points (on admettra que la régression médiane sépare l'espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue")

c'est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c'est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue",lwd=2)
for(i in seq(-2,3,by=.5)){
abline(a=i,b=3.2,col="blue",lty=2)}

Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a

d=function(h) sum(abs(y-(h+3.2*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d)(H)
plot(H,D,type="l")

Formellement, l'optimum est ici
> optimize(d,lower=-5,upper=5)
$minimum
[1] -0.2000091
 
$objective
[1] 2.200009

Bref, retenons cette courbe, et notons qu'elle passe par un des points,

Et c'est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de , on passe de à (tant que l'on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de , (là encore tant qu'il n'y a qu'un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s'arrêter dès qu'on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.

Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=optimize(d,lower=-5,upper=5)$minimum,b=3.2,
col="blue",lwd=2) points(x[1],y[1],cex=1.8,lwd=2,col="blue") for(i in seq(-1,5,by=.25)){ abline(a=(y[1]-i*x[1]),b=i,col="blue",lty=2)}

La distance est alors, en fonction de la pente de la droite

d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d2)(H)
plot(H,D,type="l")

Là encore on peut formaliser un peu

> optimize(d2,lower=-5,upper=5)
$minimum
[1] 2.500018
 
$objective
[1] 1.500018

Et si on regarde cette dernière droite, on passe par deux points,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
h=optimize(d2,lower=-5,upper=5)$minimum
abline(a=(y[1]-h*x[1]),
b=h,col="purple",lwd=2)

Pour comprendre pourquoi, comparons les deux cas,

  • en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l'autre, elle diminue

  • en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l'autre, elle diminue, mais en sens inverse par rapport au cas précédent,

Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c'est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu'à l'atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c'est toujours pareil...

Friday, April 13 2012

Ex-æquo aux élections ?

La saison 1 de Mad Men se termine avec l'élection de 1960, Nixon contre Kennedy, qui s'est jouée dans un mouchoir de poche.

L'autre jour, Benoit Rittaud évoquait sur http://images.math.cnrs.fr/ la possibilité que des candidats soient ex-æquo dans une élection, et les aspects légaux (car visiblement la constitution française n'a pas vraiment prévu ce cas). Le plus intéressant aurait été - à mon avis (mais mon avis est fortement biaisé) - de calculer la probabilité qu'un tel évènement survienne. Par exemple, qu'à l'issu du premier tour, le deuxième et le troisième candidat soient ex-æquo, de telle sorte qu'on ne puisse déterminer qui ira au second tour. Justement, un des commentaire affirme "la probabilité qu’un tel événement arrive est environ l’inverse du nombre typique de voix séparant le deuxième homme (ou femme) du troisième". N'ayant réussi à démontrer ce résultat, je me suis dit qu'on pourrait se lancer dans des simulations pour quantifier cette probabilité.

Commençons par supposer qu'il y a quatre candidats (et la possibilité de s'abstenir). On suppose que le taux d'abstention est de 20%, et que les intentions de vote pour chacun des candidats soit de 20%. 

EXE=function(N=100,ns=1000000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("A","B","C","D","E"),size=N,prob=rep(.2,5),
replace=TRUE)
tb=table(M)
Ms=rev(sort(tb[-(names(tb)=="A")]))
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Le soucis est que faire plusieurs centaines de millions de tirages de plus de 40 millions de lois multinomiales pourrait prendre du temps (genre vraiment beaucoup). Une stratégie intuitive peut être de faire des simulations pour une population plus faible, puis d'extrapoler. 

E=function(N) EXE(N)
taille=c(20,50,100,200,500,1000,2000)
proba=Vectorize(E)(taille)
plot(taille,proba,type="b",log="xy")
 
base=data.frame(y=proba,x=taille)
reg=lm(log(y)~I(log(x)),data=base)
s2=summary(reg)$sigma
m=exp(predict(reg,newdata=
data.frame(x=41194689)))*exp(s2^2/2)

L'extrapolation linéaire n'étant peut-être pas parfaite, on peut tenter un ajustement quadratique (comme cela a été fait sur la figure ci-dessus)

reg2=lm(log(y)~I(log(x))+I(log(x)^2),data=base)
s2=summary(reg2)$sigma
m2=exp(predict(reg2,newdata=data.frame(x=41194689)))*exp(s2^2/2)

La différence de prédiction avec les deux modèles est significative comme on le voit sur le graphique. Toutefois, avoir des ex-æquo avec une très très grande population est tellement rare qu'il faut réellement beaucoup de simulations pour avoir un ordre de grandeur valide. L'extrapolation semble être une alternative valide ici (mais je suis ouvert à toute critique, ou discussion).

Maintenant que l'on a une méthodologie (ou presque), utilisons des ordres de grandeurs réalistes. Par exemple en se basant sur l'élection présidentielle de 2002. 41.19 millions d'électeurs, 28.4% d'abstentions, et pour les votes exprimés, un candidat en tête avec 19.9% des votes, un second à 16.8% et le troisième à 16.2% (et les autres plus loin). Les intentions de votes pour les principaux candidats figurent dans la loi multinomiale ci-dessous

EXE=function(N=100,ns=1000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("abs","JC","JMLP","LJ","FB","AL","JPC"
,"NM","autres"),size=N,prob=c(.284,.1988*.716,
.1686*.716,.1618*.716,.0684*.716,.0572*.716,
.0533*.716,.0525*.716,.1714104),replace=TRUE)
tb=table(M);
Ms=rev(sort(tb[-(names(tb)%in%c("abs","autres"))]));
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Si on en croit l'affirmation du commentaire, la probabilité d'avoir eu des ex-æquo (entre le second et le troisième candidat) serait de l'ordre de 

> 1/(41194689*(.1686*.716-.1618*.716))
[1] 4.985823e-06

Je dois avouer que les calculs n'ont rien donné (on verra tout à l'heure ce qui se passe si on change la définition d’ex-æquo). Numériquement, les calculs tournent encore (après 4 jours) sur le serveur... Par contre, si la population ne contenait que 10000 électeurs, la probabilité d'avoir des ex-æquo (entre le 2ème et le 3ème candidat) serait de l'ordre de 3 sur 10000 (avec un intervalle de confiance que je n'ai toujours pas calculé)

> EXE(10000,100000000)
[1] 3.245562e-05

A titre d'illustration, on peut regarder ce qui s'est passé, commune par commune en France (les données sont en ligne sur le site du Ministère de l'Intérieur). L'exercice est biaisé parce qu'il y a de l'hétérogénéité spatiale... mais ça permet toujours d'illustrer... Les nombres de voix obtenues, pour le second et le troisième candidat, par commune est le suivant

> election=read.table(
+ "/Users/UQAM/Documents/Pres2002Tour1.csv",sep=";")
+ election=election[-1,]
> VOIX=election[,seq(14,74,by=4)]
> DEUXIEME=TROISIEME=TOTAL=rep(NA,nrow(VOIX))
> for(i in 1:nrow(VOIX)){
+ u=rev(sort(as.numeric(as.character(
+ election[i,seq(14,74,by=4)]))))
+ TOTAL[i]=sum(u)
+ DEUXIEME[i]=u[2]
+ TROISIEME[i]=u[3]
+ }

Sur 37000 communes,

> nrow(election)
[1] 37513

ou plutôt sur les 24000 communes qui ont eu plus de 5000 votes exprimés

> sum((TOTAL>5000))
[1] 24230

il y a eu 56 cas avec de deuxième et troisième candidats ex-aequo,

> base=data.frame(TOTAL,DEUXIEME,TROISIEME,COMMMUNE=election[,3])
> selection=(TOTAL>5000)
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 56

On a même le détail,

> sousbase[indice,]
TOTAL DEUXIEME TROISIEME                     COMMMUNE
1134   5819      699       699                        Crouy
1507   5386      717       717          Louroux-Bourbonnais
2603   6235     1151      1151                  Saint-Peray
2916   6346      935       935                        Coucy
2970   7529     1413      1413                   Haraucourt
5089   5567      717       717                      Vignats

(etc). On peut alors regarder, en fonction de la taille des communes (ici le nombre de votes exprimés) la probabilité d'avoir un ex-æquo, dans une commune donnée. On va subdivisionner en 50 classes (pour commencer), et compter à chaque fois le nombre d’ex-æquo,

> n=50
> niveau=seq(1/(2*n),1-1/(2*n),by=1/n)
> base$CLASSE=cut(TOTAL,quantile(TOTAL,seq(0,1,by=1/n)),
+ labels=niveau)
> base$EXAEQUO=(base$DEUXIEME==base$TROISIEME)*1
> m=by(base$EXAEQUO, base$CLASSE, sum)

Par exemple pour les 2% des plus grosses villes, on a eu 5 cas de communes avec des ex-aequos 

------------------------------------------------------------
base$CLASSE: 0.99
[1] 5

Si on regarde plus en détails, ce sont les communes avec plus de 7500 votes exprimés,

> quantile(TOTAL,.98)
98%
7427.76
> selection=(TOTAL>quantile(TOTAL,.98))
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 5
> indice=which(sousbase$DEUXIEME==sousbase$TROISIEME)
> sousbase[indice,]
TOTAL DEUXIEME TROISIEME       COMMMUNE CLASSE EXAEQUO
2970   7529     1413      1413     Haraucourt   0.99       1
14359  7496     1392      1392          Renac   0.99       1
25028  8104     1381      1381         Heloup   0.99       1
30014  8087     1340      1340         Barnay   0.99       1
36770  7435      764       764 Pont-sur-Yonne   0.99       1

On peut représenter ces nombres de cas d'ex-æquo et faire une régression de Poisson (car on fait du comptage ici),

> taille=length(TOTAL)/n
> plot(quantile(TOTAL,niveau),m/taille)
> sb=data.frame(niveau,exaequo=as.vector(m))
> library(splines)
> reg=glm(exaequo~bs(niveau),family=poisson,data=sb)
> u=seq(0,1,by=.001)
> lines(quantile(TOTAL,u),predict(reg,
+ newdata=data.frame(niveau=u), + type="response")/taille,lwd=2,col="red")

pour des facilités de lecture, on a en ordonnées des probabilités et en abscisse les tailles moyennes des villes

On peut faire la même chose avec 200 classes,

Bref, pour une commune de taille moyenne (disons moyenne supérieure, avec 5000 votes exprimés), on a 2 chances sur 1000 d'avoir un second et un troisième candidat ex-æquo. Ce qui laisserait entendre que pour 40 millions, la probabilité devrait être faible, voire très faible...

Bon, au delà du problème conceptuel des ex-æquo (ou légal, car je doute qu'il soit possible d'avoir une égalité parfaite, le Ministère de l'Intérieur ferait un recompte et trouverait forcément un résultat différent), il semble que la probabilité semble beaucoup plus faible que celle obtenue par la règle ad-hoc....

Mais peut-être serait-il réaliste de dire que l'élection n'a pas su dégager une majorité (c'est le but d'une élection il me semble) si le nombre de voix séparant deux candidats est trop faible. Par exemple, si moins de 100 voix séparent les candidats, on refait un vote ! Pour revenir sur mon exemple précédant, avec 10000 votants

> EXE100(10000,1000000)
[1] 0.012641

soit une chance sur 100 (je me suis autorisé à faire moins de simulations car l'évènement me semble moins improbable... donc il y a besoin de moins de simulations pour l'observer à l'occasion). Si on va plus loin, et que l'on fait des simulations que l'on projette, comme auparavant, pour estimer la probabilité d'avoir des ex-æquo non pas pour une ville, mais au niveau de la France entière, on est plutôt aux alentours de

> (m=exp(predict(reg,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
1.505422e-12 

pour un modèle linéaire, et

> (m=exp(predict(reg2,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
4.074588e-98  
pour un modèle quadratique (avec la confiance - relativement faible - que l'on peut accorder à une telle projection). Ce qui donne une différence énorme... Si on visualise tout ça,
lines(c(1,100000000),exp(predict(reg,newdata=
data.frame(x=c(1,100000000))))*exp
(s2^2/2),lty=2,col="light blue")

que l'on représente en bleu, en vert un ajustement quadratique

lines(10^(seq(1,8,by=.1)),exp(predict(reg,newdata=
data.frame(x=10^(seq(1,8,by=.1)))))*exp (s2^2/2),lty=2,col="dark green")

alors que le calcul suggéré dans le commentaire donnerait la courbe rouge,

lines(c(1,100000000),1/(c(1,100000000)*
(.1686*.716-.1618*.716)),col="red",lty=2)

Peut-on parler d'évènement improbable ? A titre de comparaison, quand on parle de quelque chose d'improbable, ou qui arrive avec une probabilité infinitésimale, on pense au loto. A l'ancien loto (français), il fallait choisir 6 numéros parmi 49 (sans remise). La probabilité de gagner le gros lot était de 

Pour l'Euromillions, il faut choisir 5 numéros parmi 50 numéros auxquels il faut ajouter les combinaisons possibles avec les Étoiles, soit 2 numéros à choisir parmi 9. La probabilité de gagner le gros lot est alors de

Bref, tomber sur des ex-æquo lors d'une élection nationale, en France, c'est a priori plus rare que gagner le gros lot en jouant une fois au loto... On peut donc dire qu'il ne devrait pas y avoir d’ex-æquo à l'issu du premier tour !

- page 1 of 4