Identification of ARMA processes

Last week (in the MAT8181 course) in order to identify the orders of an ARMA process, we’ve seen the eacf method, and I mentioned the scan method, introduced in Tsay and Tiao (1985). The code below – to produce the output of the scan procedure – has been adapted from an old code by Steve Chen (where I included a visualization of the p-values, with the following colors)

The procedure was described in the course, last Thursday,

arma.scan=function(z,ar.max=15,ma.max=15,alpha=0.01)
{
  ym=function(z,t,m){return(z[t:(t-m)])}
  n=length(z)
  z=z - mean(z)
  cmax=ma.max + 1
  rmax=ar.max + 1
  corref=matrix(0,nrow=rmax,ncol=cmax)
  cmj.table=matrix(0,nrow=rmax,ncol=cmax)
  pv=matrix(0,nrow=rmax,ncol=cmax)
  mark=matrix(rep("X",(rmax)*(cmax)),nrow=rmax,ncol=cmax)
  Rnames=paste("AR",0:(ar.max),sep="-")
  Cnames=paste("MA",0:(ma.max),sep="-")
  rownames(corref)=Rnames
  colnames(corref)=Cnames
  rownames(cmj.table)=Rnames
  colnames(cmj.table)=Cnames
  rownames(pv)=Rnames
  colnames(pv)=Cnames
  rownames(mark)=Rnames
  colnames(mark)=Cnames
  for (m in 0:ar.max)
  {
   m1=m+1
   for (j in 0:ma.max)
   {
   j1=j+1 
   if (m == 0 && j != 0)  
   {
      racf=acf(z,plot=FALSE)$acf[1:(j+1)]    
      lamb=racf[j+1]^2    
      corref[m1,j]=round(lamb,4)
      dmj=1 + 2*sum(racf[1:j]^2)
      cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)
      pvalue =pchisq(cmj,1,lower.tail=FALSE)
      pv[m1,j]=round(pvalue,4)
      cmj.table[m1,j]=round(cmj,4)
      mark[m1,j]=ifelse(pvalue > alpha,"O","X")    
    } 
    else if (m != 0 && j == 0) 
    {
      racf=pacf(z,plot=FALSE)$acf[1:(m+1)]
      lamb=racf[m+1]^2
      corref[m1,j1]=round(lamb,4)
      dmj = 1
      cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)    
      pvalue =pchisq(cmj,1,lower.tail=FALSE)
      pv[m1,j1]=round(pvalue,4)
      cmj.table[m1,j1]=round(cmj,4)    
      mark[m1,j1]=ifelse(pvalue > alpha,"O","X")
    } 
    else
    {        
      mat1=matrix(0,nrow=m1,ncol=m1)
      mat2=matrix(0,nrow=m1,ncol=m1) 
      mat3=matrix(0,nrow=m1,ncol=m1)
      mat4=matrix(0,nrow=m1,ncol=m1)     
      for (t in (j+m+2):n)
      {
         tj1=t-j-1
         ym1=ym(z,tj1,m)
         ym2=ym(z,t,m)    

         mat1=mat1 + as.matrix(ym1)%*%ym1    
         mat2=mat2 + as.matrix(ym1)%*%ym2    
         mat3=mat3 + as.matrix(ym2)%*%ym2    
         mat4=mat4 + as.matrix(ym2)%*%ym1                
      }  
      b1=solve(mat1)%*%mat2
      b2=solve(mat3)%*%mat4
      A=b2%*%b1
      eig <-eigen(A)
      eig.val <-eig$values
      eig.val=Re(eig.val)
      eig.len=length(eig.val)
      eig.vector=eig$vectors
      lamb=min(eig.val)
      eig.vector0=eig.vector[,which.min(eig.val)]
      eig.vector0 = eig.vector0/eig.vector0[1]
      resid=(1:n)*0 
      for (t in (j+m+1):n)
      {
        z0=z[seq(t,t-m,-1)]      
        resid[t]=sum(z0 * eig.vector0)
      } 
      jm1=j + m + 1
      rx=Re(resid[jm1:n])
      racf=acf(rx,plot=FALSE)$acf[1:j]
      dmj=1 + 2*sum(racf^2)
      cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)     
      pvalue =pchisq(cmj,df=1,lower.tail=FALSE)
      corref[m1,j1]=round(lamb,4)     
      pv[m1,j1]=round(pvalue,4)
      cmj.table[m1,j1]=round(cmj,4)    
      mark[m1,j1]=ifelse(pvalue > alpha,"O","X")
    }
   } 
  } 

  cat("\n\nSCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)\n\n")
  cat("Estimates of Squared Canonical Correlation \n\n")
  print(corref)
  cat("\n\nC(m,j)\n\n")
  print(cmj.table)
  cat("\n\nChi-Square(1) Test p-value\n\n")
  print(pv)
  cat("\nSCAN Matrix \n\n")
  print(mark)

plot(0:1,0:1,col="white",xlim=c(0,nrow(pv)-1),ylim=c(0,ncol(pv)-1),axes=FALSE,xlab="AR",ylab="MA")
axis(1); axis(2)
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")[c(1,2,3,5)]
cpv=matrix(as.numeric(cut(as.vector(pv),c(-1,.01,.05,.1,2))),nrow(pv),ncol(pv))
for(i in 1:nrow(pv)){
for(j in 1:ncol(pv)){
 polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
 col=CL[cpv[i,j]])
}}
}

Consider the following simulated time series,

> s=arima.sim(n=200,model=list(ar=c(0,0,0,.4,0,0,0,.5),ma=c(0,0,1))) 
> plot(s,type="l")

The output is here

> arma.scan(s,6,6)

SCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)

Estimates of Squared Canonical Correlation 

       MA-0   MA-1   MA-2   MA-3   MA-4   MA-5   MA-6
AR-0 0.0614 0.0104 0.1862 0.3516 0.0971 0.0128 0.0000
AR-1 0.0302 0.0294 0.1501 0.0943 0.0855 0.0127 0.0385
AR-2 0.3070 0.2781 0.2140 0.0006 0.1589 0.1884 0.2243
AR-3 0.1627 0.0037 0.1927 0.2311 0.1379 0.0207 0.0376
AR-4 0.2087 0.3947 0.3653 0.3075 0.1502 0.1364 0.1013
AR-5 0.1677 0.1219 0.0110 0.0263 0.0332 0.0350 0.0044
AR-6 0.0114 0.0485 0.0561 0.0427 0.0009 0.0089 0.0308

C(m,j)

        MA-0    MA-1    MA-2    MA-3   MA-4   MA-5    MA-6
AR-0  4.1161  0.6585 12.0315 20.6512 4.5388 0.5620  0.0000
AR-1  6.1127  1.9499  9.9356  4.9145 4.7219 0.4642  1.9015
AR-2 72.6011 19.1679 14.3512  0.0337 7.9668 9.6479 11.4573
AR-3 34.9724  0.2386 10.1620 13.4082 6.7875 0.8725  1.4071
AR-4 45.8691 27.5070 19.1422 20.2835 7.3339 5.5374  3.5874
AR-5 35.7981  8.0498  0.6280  1.3543 1.8470 1.7930  0.2338
AR-6  2.2147  3.1466  3.5990  1.9904 0.0511 0.4816  1.6440

Chi-Square(1) Test p-value

       MA-0   MA-1   MA-2   MA-3   MA-4   MA-5   MA-6
AR-0 0.0425 0.4171 0.0005 0.0000 0.0331 0.4534 0.0000
AR-1 0.0134 0.1626 0.0016 0.0266 0.0298 0.4957 0.1679
AR-2 0.0000 0.0000 0.0002 0.8543 0.0048 0.0019 0.0007
AR-3 0.0000 0.6252 0.0014 0.0003 0.0092 0.3503 0.2355
AR-4 0.0000 0.0000 0.0000 0.0000 0.0068 0.0186 0.0582
AR-5 0.0000 0.0046 0.4281 0.2445 0.1741 0.1806 0.6287
AR-6 0.1367 0.0761 0.0578 0.1583 0.8212 0.4877 0.1998

SCAN Matrix 

     MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6
AR-0 "O"  "O"  "X"  "X"  "O"  "O"  "X" 
AR-1 "O"  "O"  "X"  "O"  "O"  "O"  "O" 
AR-2 "X"  "X"  "X"  "O"  "X"  "X"  "X" 
AR-3 "X"  "O"  "X"  "X"  "X"  "O"  "O" 
AR-4 "X"  "X"  "X"  "X"  "X"  "O"  "O" 
AR-5 "X"  "X"  "O"  "O"  "O"  "O"  "O" 
AR-6 "O"  "O"  "O"  "O"  "O"  "O"  "O"

with the following graph

Of course, it is possible to ask for larger values,

> arma.scan(s,12,12)

The graph is now

Voting Twice in France

On the Monkey Cage blog, Baptiste Coulmont (a.k.a. @coulmont) recently uploaded a post entitled “You can vote twice ! The many political appeals of proxy votes in France“, coauthored with Joël Gombin (a.k.a. @joelgombin), and myself. The study was initially written in French as mentioned in a previous post. Baptiste posted additional information on his blog (http://coulmont.com/blog/…) and I also wanted to post some lines of code, to mention a model that was not used in that study (more complex to analyze, but more realistic, and with the same conclusions). The econometric study is based on aggregated voted, with a possible ecological misinterpretation.

  • Regression Model: Possible Explanatory Variables

The first idea was to model proxies using a binomial regression, per pooling station  where  denote the number of proxy vote, per station , and  denotes the number of voters. Proportion  can be a function of possible explanatory variables (on Baptiste’s blog there are additional information about the datasets, obtained from insee.fr and opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

For instance, assume that  is a function of the proportion of owner of the place people live in, denoted  in the neighborhood of the pooling station,

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference]

We can consider a logistic regression

or a logistic regression with splines, if we do not want to assume a linear model

With cubic splines, the code is

> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])

while a standard logistic regression would be

> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

Here it is (the confidence region is for the spline regression) with on blue the first round of the Presidential election, and in red, the second round (in France, it’s a two-round system)

(the legend of the y axis is not correct). We can consider as explanatory variable the rate of H.L.M., low-cost housing or council housing,

If I like the graph, unfortunately, the interpretation of coefficient  might be complicated

> summary(regt1l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset1, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.9549   -1.5722    0.0319    1.6292   13.1303  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.70811    0.01516  -244.6   <2e-16 ***
taux         1.49666    0.04012    37.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12507  on 836  degrees of freedom
Residual deviance: 11065  on 835  degrees of freedom
AIC: 15699

Number of Fisher Scoring iterations: 4

> summary(regt2l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset2, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-15.4872   -1.7817   -0.1615    1.6035   12.5596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24272    0.01230 -263.61   <2e-16 ***
taux         1.45816    0.03266   44.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9424.7  on 836  degrees of freedom
Residual deviance: 7362.3  on 835  degrees of freedom
AIC: 12531

Number of Fisher Scoring iterations: 4

So we did consider a standard linear regression model, for the proxy rate, per station,

(again, either a model with splines, or a standard linear model). The code is

> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

Here is again the evolution as a function of the rate of owner of their homes,

The graph is rather close to the one before, and here, the interpretation of the summary table is more conventional,

> summary(regt1l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.9994 -0.2926  0.0011  0.3173  3.2072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.021268   0.001739   12.23   <2e-16 ***
taux        0.054371   0.004812   11.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared:  0.1326,	Adjusted R-squared:  0.1316 
F-statistic: 127.7 on 1 and 835 DF,  p-value: < 2.2e-16

> summary(regt2l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.9029 -0.4148 -0.0338  0.4029  3.4907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.033909   0.001866   18.17   <2e-16 ***
taux        0.079749   0.005165   15.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared:  0.2221,	Adjusted R-squared:  0.2212 
F-statistic: 238.4 on 1 and 835 DF,  p-value: < 2.2e-16

We have used those codes to produce the graphs mentioned in the post. But before mentioning the residuals of the multiple model we considered, I wanted to share some awesome code that produce maps (I can say that those codes are awesome since Baptiste wrote most of them).

  • Visualization of Residuals on a Map of Paris

To plot the neighborhood of the pooling stations, one more time the post on Baptiste’s blog, explains how the shapefile was obtained from cartelec.net

> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")

To visualize the proxy rate (the average of round one and round two), here is the code

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

To get nice colors, function of the rates, we use

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
>(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

and finally

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

If we consider a model with three explanatory variable, to explain the proxy rate,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

we can plot the residuals using

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
>(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

It might not be a pure random spatial noise… But we could not get better with our small set of covariates.

Bivariate Densities with N(0,1) Margins

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are  distributions.

  • The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins , and with correlation . This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
contour(vx,vy,z,col="blue",add=TRUE)
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green")
lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
dnorm(density(X[,1])$x),col="red")
par(mar=c(3,0,1,1))
barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0, 
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks[1])/
diff(yhist$breaks)[1],col="red")

That was the simple part.

  • The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair , we would like to consider the pair , so that the marginal distributions are . The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z)

Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

  • The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where  and  with independent random variables, given , and that both variables are exponentially distributed, with parameter . As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be . The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that  is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value))
Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy[2],z)*dexp(xy[1],z)*
dgamma(z,a,b),lower=0,upper=Inf)$value
f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y))
vx=seq(-3,3,length=151)
vy=seq(-3,3,length=151)
z=matrix(NA,length(vx),length(vy))
for(i in 1:length(vx)){
for(j in 1:length(vy)){
z[i,j]=f(vx[i],vy[j])}}
set.seed(1)
Theta=rgamma(1500,a,b)
Z=cbind(rexp(1500,Theta),rexp(1500,Theta))
X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2]))

There is a small technical problem, but no big deal.

Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian.

Academic Blogging, a Personal Experience

I wanted to get back on two posts (one on Academic Blogging, and one on Twitter for Academics), and to update them based on the discussion that followed the panel (on Thanksgiving), as well as some more recent discussions. All comments are welcome ! [the post is quite long, a pdf version is available from papers.ssrn.com/2398377]

  • Introduction

To talk in public, to think in solitude, to read and to hear, to inquire and answer inquiries, is the business of a scholar” (from Johnson (1759)). 250 years after, anyone might think that the business has not changed: academics remain in their ivory tower, but from time to time, they have to leave it, to communicate. Recently, Graham (2004) characterized three channels for scholar communications: publication in peer-reviewed journals, conferences and seminars, and a more informal one, that might be called the “new invisible college“, as in  Halavais (2006), that might be related to a “faculty lounge“, using the expression of  Priem, Piwowar, & Hemminger (2012). This third channel is precisely the one that might be related to blogs.

With Internet (emails, blogs, forums, etc), academics now have new mediums to communicate, either within their own community (and launch participative projects), or outside their community. Academic blogs are one medium, among many others.  As explained in Gregg (2006), “blogs have made scholarly work accessible and accountable to a readership outside the academy“. From an insider’s perspective, blogs seem to be extremely popular, because bloggers are active, sharing links, comments, discussions, etc. In 2007, George Siemens was already enthusiastic: “it’s great to see research-focused academics entering the blog space” (see Siemens (2007)). By that time, it was  clearly not as popular as it could be seven years after.

Social medias can be categorized in eight categories, according to Gu & Widén-Wulff (2010) and  Nicolas & Rowlands (2011): blogs, microblogs, RSS, wikis, tagging, social networks, media sharing, and online documents.  For instance, there is also an emerging form of participatory media, with question and answer services, with Quora, and Stackoverflow. But it won’t be mentioned here (some researchers are extremely active there, but it is very difficult to quantify). In this article, we will get back on more than five years of experience, on a scholar blog.  In the first section, we will get back on definitions, on the context, and the origins of the frekonometrics blog, while it started as one of the few official blogs of researchers at Université de Rennes 1, in France. Then, in the second section, we will study the interplay between a blog and other activities within academia. Then, in a third section, there will be a short discussion about microblogging, and the use of Twitter for academics. And finally, the last section will try to explain why, even after almost six years of experience, an academic blogger might still be enthusiastic.

  • The Origins and the Context of the Blog

From a technical point of view, blog is a contracted form of weblog, which is a website made up of ongoing entries, that we will call posts (but might also be called articles). Those posts are published in reverse chronological order: the new material is added to the top of the page, the older materials being automatically archived. So it makes it difficult to follows for infrequent readers. There might be tags and categories, which can be used to distinguish posts related to conference, publications, and teaching, or using dedicated keywords.

  • Hosting Platforms and Identity of the Blog

The first blog was launched because the Université de Rennes 1 (I was just recruited there as a professor, not yet tenured) decided in 2007 to have its own blog platform for its researchers. It was officially Arthur Charpentier’s blog. On this platform, blogs were not anonymous, they were based on the credibility of the researcher, and the name of the university was on the front page, almost as large as the name of the researcher.

Being hosted by the university is a dangerous practice, since it is not clear what is allowed, and what is not. For instance, during the first year, there were a lot of discussions about the status of researchers, in France. Being not tenured put me in a delicate situation. Even if the blog became well referenced on search engines very fast, it was not very popular, by that time. Except perhaps among students (mine, but also students in other programs). Students understood the interest of the blog, as a place to discuss and to interact. The blog became and extension of the class. After the formal lecture, in the room, the blog became a place to share additional documents, datasets, computer codes, etc. After two or three years, the blog started to be popular, within my own community (in economics, econometrics and applied mathematics), not only among students. I started to be recognized in conferences, and I wanted to stop having an eponymous blog (coincidence, or not, it was also by the time I moved to Montréal). Having an eponymous academic blog is a common practice, which makes sense since most bloggers tend to identify their blogs, as both personal and professional. All the more within academia, where boundaries between professional, academic and personal life may be difficult to establish, mainly because all those aspects intertwined constantly in the life of scholars.

Again, having a blog hosted by the university might yield delicate situations. Those institutions would like to use blogs as communication tools, and might not appreciate when the blog is use to publish satire of one’s own department. After a few years, I wanted to transfer the blog somewhere else. Using the term introduced by Charlotte Frost, used in  Mewburn & Thomson (2013), I wanted my blog to be – still – some sort of academic blog, but an “outstitutional” blog, more than an institutional one. I decided to host the blog on my personal internet provider’s platform. After almost three years at Université de Rennes1, the freakonometrics adventure started officially.

With a personal blog, it became necessary to claim that the blog was some sort of academic blog. Marsh (2013) mention that “many people arrive and depart my blog without knowing I’m an academic. I don’t hide the fact. It’s there“.  The name of the researcher appears in the “About page“, or in the resume, which is available on the front page. There is the blogger, freakonometrics and the scholar, Arthur Charpentier. Those two identities overlap, considerably, without being entierely congruent. When one reaches an academic webpage, there is usually no surprise: we expect to find information about the academic background, a list of publications, some information about trained students, lectures notes. But all academic bloggers have their own identities, which is a reflect of its author. Some bloggers like to post pictures taken from places they visit. Others talk about their family, or articles they recently read. Because blogs might have diversified contents, they tend to attract very different audiences (we will get back on that point in the next section). I do not hide myself behind a pseudo (even if I use one) and I still mention my academic background in order to give my blog more credibility. Noel (2013) claims that “a blog does not indicate expertise“.  Traditionally, it is true that a stack of peer-reviewed work automatically indicates some expertise. But it might be time to assess that a well-read blog, with valuable contributions in some field can also indicate some expertise, too.

Several platforms were hosting blogs by that time. Some interesting blogs were hosted by scientific journals (such as The American Scientist or Nature) or societies (in Economics, or Statistics).  But those blogs are somehow institutional. Even if my area of expertise if on the border of social sciences, I decided to migrate to the hypotheses.org platform, a “platform for academic blogs in the humanities and social sciences“. This platform has its own scientific committees, and provides a technical support (not to mention the legal one).

  • Influences and the Blogging Community

Freakonometrics became a blog dedicated to econometrics, economics and applied mathematics. As explained in Quiggin (2011a), “with the arguable exception of law, economics is the academic discipline where blogging has been embraced most enthusiastically“. This might explain why there is such an active – and enthusiastic – community (see also Quiggin, (2011b) for a discussion). On the other hand, Halford (2012) said that “the situation within academic blogging seems to be that we are currently a bunch of islands that are vaguely connected but not really arranged into continents and groups. We are all spread out across the digital world with a fragmented network between us.

So, what can we find in this community? Some researchers like to use their name, like Greg Mankiw’s blog for instance (with subtitle “random observations for students of economics“). But most of them prefer to hide themselves behind a short title, like “Confessions of a Supply-Side Liberal” (by Miles Kimball), “The Conscience of a Liberal” (by Paul Krugman), or longer one, like “Statistical Modeling, Causal Inference, and Social Science” (by Andrew Gelman). But hide is probably too strong since the editor is never hidden: we usually find a short bio, including a picture (most of the time), as well as a link to a webpage hosted by some university. Other examples might be “I’m a bandit“, with subtitle “random topics on optimization, probability and statistics. by Sébastien Bubeck“, or “what’s new“,  with subtitle “updates on my research and expository papers, discussion of open problems, and other maths-related topics. by Terence Tao“. Some blogs use puns (it is a feature that one can find on almost any blog: most of them use humour, just to explain that this is just a blog) like “Hyndsight“, “a blog by Rob Hyndman“.

Most of those popular blogs are blogs of well-know academics. Greg Mankin is not only a popular blogger, he is also the author of a standard undergraduate textbook in macroeconomics, and so is John Cochrane, editor of the “Grumpy Economist” blog, also author of a standard graduate textbook in financial economics. There are many other well-known academics, willing to reach new audiences. But John Quiggin claims that his own blog is much more popular than his research (on Australian policy issues), as well as Tyler Cowen (one of the contributor of the “Marginal Revolution” blog).

Another category might be blogs with several contributors, such as the “Monkey Cage” blog, dedicated to political science research. This blog contains contributions from several researchers, such as John Sides, Erik Voeten, Andrew Gelman, Joshua Tucker and Henry Farrell. But all those researchers have their own blogs besides.

To get back on my personal experience, the name freakonometrics was chosen for two reasons. The first one is that is was a common word among the community of researchers in the Economics department at Ecole Polytechnique. Freakonometrics meant using Econometrics techniques to answer real life questions. That might be opposed to Mathematical Econometrics, as it was taught to students. Since it was exactly the aim of the blog (write posts about simple questions, with a quantitative approach) I liked the idea of using it at a blog’s name. The second reason is related to a major influence for most bloggers in Economics. In 2005, University of Chicago economist Steven Levitt and New York Times journalist Stephen Dubner published a collection of “economic” articles, claiming that economics is, at root, the study of incentives. This is how freakonomics.com/blog/ started (the first post was published in September 2005). My blog was more about Econometrics than Economics. So I decided to use the freakonometrics name. But this was not a very good idea. If all bloggers know freakonomics, there is no doubt that the two blogs are unrelated. But unfortunately, outside this community, there is confusion (as I noticed many time while discussing with journalists for instance).

  • The Structure of the Blog, and the Contents

A blog is usually seen as chaotic, and unstructured. On the opposite, academia, and Science, are places where everything is supposed to be well ordered. Most undergraduate students can not understand where a professor answer “I have no idea what might happen if we relax that assumption“. An open question is a great source of inspiration for a blogger. It might be a start for a bibliographic survey, or on quick simulations. Thus, the blog “offers a unique space to simultaneously achieve efforts related to research, teaching, and service” as mentioned in Grollman (2014). Technically, as recalled in the introduction, a blog is just a series of entries called “posts“, even if, with several thousand words, some posts might also be called “articles“, as suggested in Cohen (2006).

Those posts might be to start a discussion after an open question in the class, to mention in interesting paper recently discovered, to point out an interesting conference, to provide some codes to generate a graph, to upload datasets used in an article about to be published, to criticize an article read in a newspaper, or just to share an experience. Therefore, “academic blog” is a generic term that includes several different genres, from blogs about academic life to pure research blogs (see Walker (2006) for a discussion), even if most of them are hybrid genres. As analyzed in Luzón (2006), making available personal research is a simple way not only to get feedback, but also to increase visibility, to develop “respect and reputation” (as discussed in Gregg (2009)). As claimed by the committee of hypotheses.org, blogs should be seen a lab books, not a substitute to academic journals (a post is not an article, as we will see in the next section) but more a catalogue of personal thoughts (and possibly some findings).

  • Academic Life and Blogging Activities

Traditionnaly, according to Merton (1942), academia was governed by some norms of universalism, communism, disinteresredness, and organized scepticism. More recently,  Ziman (1996) claims, on the contrary, that academia was characterized by inverse norms, of  proprietary, local, local, authoritarian, commissioned and expert norms. Blogging is probbly an answer to this shift. In this section, we will get back on interactions of blogging activity with more conventional academic duties.

  • Academic Hierarchy and Rules

Walker (2006) explains that  “blogging allowed us to circumvent the power structure of academica“. Blogs are somewhere between the conversational style of the conference, and the writen style of articles. But blogs are clearly outside those communities, blogs are much more accessible: you do not have to attend a conference, you don’t even have to sign up, and you are usually allowed to participate in the discussion, without belonging to the community (you should if you want to publish a response to a published article). Furthermore, the bureaucratisation of the universities has had “profound effects on the writing of academics” (as discussed in Brett (1991)). More specifically, it is nearly impossible for academics to provided a “public intellectual function” nowadays because it goes “against the grain of the job“. Blogs can be seen as an answer of academic, to regain some freedom of speech. As mentioned in Gordon (1998), tenure positions were initially created precisely to prevent academics from pressures, to grant them a strong freedom of speech. But in the 2010’s, tenure positions exist to pressure those who have only temporary contracts.

  • Blog Posts and Academic Articles

A blog entry is usually published without any pressure, with no one too impress, just the idea of working through some ideas and concepts, together, the writer and the reader. Academic articles are going through the peer review process. But it is only a peer-to-peer discussion. All those discussions are usually skipped in the final version, and there is no real place for a discussion, except for some controversial papers, where a discussion can be published. As mentioned in Chatterjee & Biswas (2011), some journals now have forums, but there is no place, otherwise to converse about a publish paper. Blogs can be seen as tools for “post-publication review of papers“. This is actually what happened in 2010, when Science published a paper (Sebastiani & al. (2010)), which was then criticized by several science blogs – about methodological fallacies, before being retracted (see Carmichael (2010) and MacArthur (2011) for a full discussion).

Blog posts offer a space not only to summarize a research paper (written in a codified and strict format), but also to discuss the process and the context of the study. As Gregg (2009) suggests, blogs provide a space “distinct from the parent culture of institutions“. A paper might start with a discussion with a colleague, then additional contributors might join to provide additional expertise. The paper is then written in a codified format, and to please potential reviewers. Those reviewers might ask for substantial changes, done because some co-authors really need a publication to apply for a position. Finally, the paper published might have to connexion at all with the initial discussion. Even if bloggers do not want to share with complete strangers the hidden side of a paper, it might be interesting to emphasize dead-ends, attempts that failed, and how the paper started.

  • An Academic Blogger is not a (Science) Journalist

Blogging can revamp the relationship between science and society, simply by making research more comprehensible to others, just like journalists do. There have always been connexions between the two communities: academics, who can understand what exactly is claimed in a research study, and science journalists, who can summarize and explain with less jargon. “Historically, academics have been put in a reactive position, responding to questions from reporters. Blogging places academics in a more proactive position, intervening more effectively in popular debates around the topics they research” as claimed by Jenkins (2008).

Actually, some academics do publish posts in blogs hosted by newspapers, such as “The Conscience of a Liberal” (by Paul Krugman). Somehow, those journals (here the New York Times) host academics the same way newspapers hosted some opinion pages, were academics where invited to give their point of view, a few years back. But, as Quiggin (2011) explains “newspapers are generally reluctant to repost on academic working papers and similar publications unless the conclusions are obviously newsworthly“. Economics in newspapers has to be related to macroeconomics, and science has to be related to medicine or technology (the science page is now a nice advertising page for the most recent smartphone applications). And, as mentioned by Cottrell (2013), a couple of decades ago, experts (not to say scholar) “have functioned as sources for newspaper journalists. Their opinions would emerge often mangled and simplified, always truncated, in articles over which they had no final control.” Now, with blogs, it is possible to read them directly, in a style that is easier to understand, compared with academic (standard) publications. “The general reader has access to expertise that was easily available, a decade ago, only to the insider or the specialist“. Writing a post in an academic blog is neither pretending being a journalist, nor writing an academic article. In the traditional process of research, we discuss with colleagues, possibly in conferences, but only the final publication remains. False starts and heuristics are skipped, because they might appear as un-necessary to get an understanding of the article. And this is exactly where blogging becomes interesting. In newspapers, contributors are required to avoid mathematical equations, and scientific jargon. Blogs are exactly the missing link between a technical research article, “jargonistic and turgid” (as Mewburn & Thomson (2013) wrote), and an article in a newspaper, usually overselling minor contributions.

  • The narrative style of a blog post

There is a stark contrast between the story used in a peer-reviewed paper, and the true, untold, research process. The messy side of a research project, and all the dead-ends might be mentioned, maybe as a quick anecdote, in a conference talk. But for junior researchers, it is not the place to discuss the true story (maybe only the shinny side of it). As mentioned in  Conole (2007), “the blogosphere is offering an alternative style of academic discourse“. A blog post is a place for snippets of the work, more informal, with a completely different style. As climed by Pierre Mounier and Marin Davis (from hypotheses.org) the blog is between the oral and the writen form.

Nevertheless, the tone, or the style of writing, is different from one blog to another one. Post can be writen as a reportage, with a journalistic style (where various sources of information are brought together and synthetised into a story), as a formal or more informal essay (where a story is writen, with opinions), it can be more pedagogical (with a lot of questions) or satirical (see  Mewburn & Thomson (2013) for a detailed comparison, with further examples).

But the time when Karl Sagan was criticized (see e.g. Emanuel (1986)) for mentioning a study in Parade (the  American nationwide Sunday newspaper magazine) before publishing it in Science is now far away. As written in Oreskes & Conway (2011), this has been considered as a”violation of scientific norms“. Similarly, because of its style, some researchers cannot consider a blog post as a serious contribution.

  • Blogging and Microblogging (Twitter)

Blogs are not the only place where academics are now visible. Academics are also extremely active on microblog platforms, such as Twitter. Twitter is different from blogs. Not only because of the lengths of the “posts“, but also because it is extremely codified (while blogs are clearly not). If academics can use blogs to educate, they can use Twitter to inform, to make interesting material available to the public.

  • From Official Rules to Unofficial Uses

Microblogging, on Twitter, is simple. One has to create an account, on twitter.com, and 160 characters can be used for a short bio, with the possible use of an avatar. It is possible to be active, by tweeting (posting online messages with less than 140 characters), or to passive (and then simply follow some discussions). One can use the search window, and then all the tweets containing that word (or that sentence) will appear. It is also possible to follow some hastags (following the # symbol) such as #overlyhonestmethods, or some people (following the @ symbol) such as @freakonometrics. When following that account, all the tweets (posted under the name @freakonometrics) will appear in the so-called TL (or timeline). I will then have a follower (and of course, I might follow back). On can tweet not only text, but also a picture and html links (that will be automatically shortern). Those are the official rules in Twitter.

There are also codes and uses. With a Twitter account, it is possible to write our own tweets, but also to retweet a tweet, which means forwarding another user’s tweet to all of our followers. One can use the retweet button, or use the old fashion RT (for re-tweet) or MT (for modified tweet) if it was necessary to shorten a tweet.  As in Academia, it is important to mention the source of the information when sharing it with some followers (it is called a mention, and it is like using a reference in a research paper). Another way of acknowledging the account that originally shared the content being tweeted is to HT the account (hat tip). It is also possible to discuss with other people. One can reply to a tweet using the reply button (when replying to one of my tweets, the tweet will start with @freakonometrics and only followers also following me will see it in their tweet list). It is also possible to poke another account to make sure that someone reads your tweet. It is possible to favorite a tweet, by clicking the yellow star next to the message. If a lot of the people I follow favorite a tweet, it will be more likeliy to appear in the discover window, even if I do not follow that account. Because of the 140 characters constraint, a lot of strange words can be used on Twitter, such as “OH” which most often means “overheard“. Twitter has also its own habits, such as the #FF follow Friday tradition. On Friday, it is common to share usernames of our favorite twitterers, the accounts you find interesting

  • Are Academics on Twitter ?

Only “one in forty scholars are active on Twitter” as estimated in Priem & al. (2012). And academics are mainly skeptical about Twitter. For most of them, Twitter is either for their kids, or perhaps for graduate students (if they have some ideas about what Twitter is). But not researchers. To get back on my personal experience, I follow (and am followed) by a lot of graduate students, PhD’s or postdocs, and a few more senior researchers. Most of those senior, prominent scientists, share extremely interesting information! But Twitter is time consumming (as blogs), and not everyone is willing to have an account on Twitter (and to be active, by sharing links, or by interacting, more generally).

If not everyone is on Twitter, important people are; at least in Economics. In January, was organized the annual conference of the American Economic Association, in Philadelphia. This is where the job market for PhD students in Economics takes place. It was extremely active on Twitter (updates where obtained by following frequently the #ASSA2014 hashtag). From my perspective (maybe also the people I follow), it looked like everyone there was on Twitter during that (major) event.

  • Twitter as a Bookmark

Like most of my colleagues, I do spend a lot of time online, reading articles, for work, for fun, and sometimes, it can be interesting to keep tracks of those articles, posts on blogs, and articles on arxiv.org or papers.ssrn.com. The first reason I use Twitter is that I need bookmarks. With Twitter, those bookmarks are public. It does not mean that I endorse what I read, it should be understood as I found that interesting, and I want to bookmark it, to find it, someday. As most academics, I do read a lot. And one of our duties is to share information. This was mentioned in Priem & Costello (2010) “the professional impact of Twitter may be particularly prononced for scholars given that sharing information is a central component of their work” (see also Letierce & al. (2010)). They estimated that 30% tweets sent by academics contain a hyperlink to a peer-reviewed ressource (usually a pdf file of a research paper). And it is not necessarily a paper published online 24 hours ago : it can also be a paper rediscovered accidently, or a technical report writen a few decades ago that has just been scanned.

Traditionally, academic visibility is measured using citations, meaning that some work has been accepted by (so called) peers, in the scientific community. Academics need to write to have an impact. But a lot of time is spent on reading. This reading activity is missed by standard citation counts (unless you publish a review, for instance).  On Twitter, you can comment, even (briefly) discuss a publication, some tricks on computer codes, share graphical visualization, etc. It is not like posting an anomynous comment on some scientific blog (actually, several blogs do not have open comments any more). On Twitter, there is some kind of credibility, and not only from the academic resume: people know you and follow you. According to Wagers (2012) “good content is propagated rapidly, bad content is not“. Of course, it is not that simple.

  • Live-Tweet in Conferences

Another popular use among academics is to use Twitter for live-tweets. But as mentioned on a lot of blogs, one should be careful about live-tweeting. Live tweeting is supposed to be fun, but stay polite, and respectful, and to use quotation marks as much as possible. Getting back on the so-called Twittergate (see McMillan Cottom (2012) for a complete summary), Aaron Bady, used that interesting image, about Twitter within the academic community “I conjured up the image of an appropriate cantankerous old professor yelling at a bunch of punk tweeters to get off his lawn, like Clint Eastwood in Gran Twittarino”.

Again, to get back on my personal experience, I usually do not live-tweet, I do not feel comfortable with it (I prefer to take notes in my book, even if I might also write a post on my blog later on, and I always ask the speaker if I can quote what he or she said). “Some worried that having someone tweet their insights before they publish might increase the likelihood that they will be scooped by a colleague — although others regarded that notion as slightly paranoid” as mentioned in Kolowich (2012). “The debate over live tweeting at conferences is, in many ways, about control and access: who controls conference space, presentation content, or access to knowledge?” wrote Risam (2012). Beyond those ethical considerations, there is also a more practical reason: in mathematics, it is quite difficult to live-tweet. It is difficult to write equations in Twitter, and to mention a graph without the formal model is useless (or misleading). There might some interest when there are computational issues, to share a nice visualization for instance.

  •  Twitter as an Online Faculty Lounge

 

It is also possible to start discussions on Twitter, but it has to remain a short discussion (to ask for precisions and references). An important issue in Twitter is to speed up connections between scientists. But there is nothing new here. Traditionally, scientists have always interacted with other scientists, in sort of one-to-one interactions, attending seminars, conferences, and discussion with colleagues. Using the words of Priem & al. (2013), “informal conversations have moved out of the faculty lounge to online social media platforms”, such as Twitter. One of the interests is to join somehow a larger “virtual department” with colleagues that are not next door, but who might be far away and in other areas of research. One can even discuss with real people, outside academia. Since I have interests in risk modeling, finance, climate, computer science, mathematics, I can also discuss with people working on stock markets, in insurance companies, in data visualization startups, even journalists. A important point is that it becomes possible to interact with open minded researchers. As mentioned in Fox (2012) – slightly changing the title – “blogging [and micro blogging] changed how economics share ideas”.

The first step in the scientific process is to find ideas, new ideas or concepts to investigate, datasets to describe. Following interesting people on Twitter can be the first step. The final step is to communicate findings and to disseminate. The time when researchers were studying the table of contents of journals to find interesting articles is behind us. When disseminating on a blog, we can share codes, graphs, datasets, links to additional material. On Twitter, we have to deal with the 140 characters constraint, which makes it hard. One idea can be to use a nice visualization, a graph, or a map.

  • Institutional Accounts and Outstitutional Ones

On Twitter, I mainly follow researchers, only a few institutional accounts. For instance, following @HarvardBiz is interesting to get updates about their blog, and recent articles (like using RSS feeds). But most of them are based on non-explicit rules, and controversial posts and articles are not mentioned. Even if there might be interesting experiences. Recently, academics behind the Nature Chemistry Twitter account went back on 4 years of experience (in Nature Chemistry (2013)). Among lessons learnt, collectively, it was mentioned that with the 140 characters constraint, “clarity is a virtue”, which is indeed one of the things you learn with Twitter. Furthermore, they mention that following important researchers in your scientific community on Twitter is interesting. Not only you can discover interesting information: when following personal accounts, one can also learn more personal information.

  • The Impact of Twitter in Standard Process of Academic Dissemination

Another possible motivation for researchers to be active on Twitter might be for dissemination. Using Twitter can help to reach other researchers, outside the (small) circle of one’s community, as well as journalists, or people working in the industry and for governmental organization, even politics. In a large study, Shuai & al. (2012) proved (following 4,600 papers) that papers mentioned on Twitter are more downloaded, and more cited (see also Eysenbach (2011) for a similar conclusion).

  • Blogging After Almost Six Years

After more than five years as an active blogger, I have the same diagnosis as Kotsko (2006), which mentions that blogging is “especially great for academics who would otherwise be quite isolated from other academics with similar interests“.

  • A Peaceful Island within Academia

To get back on one of the main personal reason why I am still blogging after almost six years, with enthusiasm, is because it is still a lot of fun. And it is a place that I appreciate all the more that academia is currently a very competitive place. I interact with the blogging community because I want to, while (most of the time) I interact with students and colleagues because I have to.  Blogs are probably the last place where we have a complete freedom. Even not to blog, if we do not want to.

There are currently two important issues in academia, when talking about money: tuition fees rise, and the decrease of public funding for research. Those past years, (undergraduate) students became consumers. And as a professor, I am now the seller in the store. Most of our students are no longer interested by the story behind a model; they want some recipes for their future job. They want to know what to use, and when. And since universities evaluate their professors, they ask students to fill evaluation forms. So professors do everything to get good evaluations. That is a simple and extremely rational game. But the pleasure of teaching is lost, sometimes. And similarly, the less money there is to do some research, the more competition. If there are one or two grants in my field of research, in Québec, I have no more colleagues and friends, I have only competitors. We now have to work on different tasks to add lines on our resume. As mentioned previously, it is difficult to find time just to discover something new, and spend some time investigating (without the insurance of getting a publication out of it).

Compared with those conventional academic activities (teaching, doing so research, writing a referee report, filling a form for a grant, etc) blogging is fun. Within the blogosphere, there is no competition, just motivation and stimulation. One can interact with other bloggers, learn from them, and so far, it is still a pleasure to blog. Some bloggers claim that it is a shame that blogging is not recognized (formally) within academia, but actually, it might be seen as a great opportunity. We do blog because we want to: it is not another required task. So it can still be fun…

  • Readers of the Blog

As claimed in Walker (2006), “blogs (…) are inherently social. Whether you have five readers or five hundreds doesn’t really matter, it’s the knowledge that you will be read that is important“. So, I have to confess that I blog mainly for myself, in the sense that I do not want to have a readership waiting for me to post something everyday! Also in the sense that my blog gives me a complete freedom to talk about things I find fun, with a whole person style (and actually I do use my blog to develop my own “writing voice“, as called in Walker (2006)), discussing about personal issues. I remember when I started blogging. At first, I thought that no one will read my blog, and then, friends and colleagues told me that they were. It is actually thrilling (not to say scary) to have 5,000 or 10,000 readers for a blog post, when one think about the number of readers of academic articles. As claimed by O’Connor (2007), “blogs are a more effective medium for intellectual influence than journal articles“. Somehow, it looks like academic journals try to avoid exposure: publishing an article in the Journal of Narrowly Focused Hyper Specialized Field Studies is a great place to hide one’s research.

The process of research is, indeed, a social activity, where we need to keep – and create – interactions with various researchers, read papers, keep our mind open. On the other hand, blogging is definitively a personal activity: I use my blog as a notebook, to keep traces of ideas, or codes, even graphs and interesting readings. From this perspective, blogging is extremely personal. But it is opened, and anyone can access it. So I use my blog to promote my work, and my scholarship. As Greg (2009), I see academic “blogging as conversational scholarship“. Blog are great to encourage conversation, they are coffee houses, in the sense of the XVIIth Century in England. In blog posts, we connect to other blogs, using comments, reactions, and hyperlinks. But actually, Derek de Solla Price (see Price (1986), cited in Horne (2011)) explained that “the prototype of the modern scientific paper is a social device rather than a technique for accumulating quanta of information“. In that sense, having informal discussions is probably the best way to work, as an academic. This is also the idea of Crane (1989), “the growth of scientific knowledge is a kind of diffusion process in which ideas are transmitted from person to person“. Using blogs, we can develop and connect a network of various people, from PhD students to practitioners in the industry, as well as more experienced academics who might share common interests. The blog is read by students, former students, colleagues, probably the dean, and even sometimes the department secretary.

  • A Cost-Benefit Analysis

With a simple cost-benefit utilitarist analysis, one will probably blog if benefits are more important than costs. One component is related to time issue: is blogging costing, or saving time? Academics receive frequently emails, asking for explanations on a technical question (from students, former students, or anyone actually) that will need a detailed email answer. It is possible to consider a “reply to public“, with a blog post. Similarly, while teaching, frequently the same question is asked every a year. From my own experience, the first year, I can write an answer in a blog post, and then, I can integrate it to my notes (blog posts can even be more interesting than lectures notes). I cannot believe that blogging is a waste of time, since I see my blog as a long-term memory. It is also possible to integrate material that cannot be used in notes, such as animations, or videos.

A lot of academics claim that “do not have time to waste blogging” but some of them can write extremely long, detailed (and most of the time, nicely structured) replies. If I write a detailed reply to a specific question (because I found the question interesting), why not share it? Either the answer is wrong (or imprecise), but then someone might post a comment, add a reference, a link, etc, or the answer is correct, and then, why not share it? With blogs, dissemination is immediate, as well as comments and feedback.

A lot of researchers within academia still reject blogs because they’re not serious, and not peer reviewed. Not being serious in the style does not mean that the study is not serious and it can still follow a scientific procedure (this is what we can learn from the Ig Nobels). Furthermore, blogs will be peer reviewed as soon as peers will blog. In blog, comments can even be more constructive than comments one gets from referees in a peer reviewed journal. Blog posts are published on the (open) web, not in some journal so expensive that no one can actually read it. Yes, commenting is not a formal or rigorous as a peer review publication, and having opened comments may contribute to establish quality and credibility of a blog. Having comments from the community is a great benefit.

  • Conclusion

With recent changes in academia, blogs are a perfect place to interact outside institutions (universities, societies and journals). Academics can find there not only a place for free speech, but also a place to interact with scholars outside the circle of hyper-specialized colleagues. Anne-Marie Slaughter recently said that “all the disciplines have become more and more specialized and more and more quantitative, making them less and less accessible to the general public” (quoted in Kristof (2014)). Blogs are the missing link from the jargon-based papers to the oversimplifying journalistic articles, between the conventional academic article, the lab notes, and the talk in a seminar. But one should keep in mind that “academic blogging can be an important medium, when it avoids the meta-narcissistic onanism of blogging about how important academic blogging is”, as claimed by Parr (2012).

  •  References

Adams, R. (2013). Blogging in context: Reviewing the academic library blogosphere. (D. 10.1108/EL-05-2012-0054, Éd.) Electronic Library , 31, 664-667.
Ashlin, A., & Ladle, R. (2006). Science communication: Environmental science adrift in the blogosphere. (D. 10.1126/science.1124197, Éd.) Science , 312.
Batts, S., Anthis, N., & Smith, T. (2008). Advancing science through conversations: Bridging the gap between blogs and the academy. (DOI:10.1371/journal.pbio.0060240, Éd.) PLoS Biology .
Bell, S. (2007, 11 1). To Blog Or Not To Blog – That Is An Academic’s Question. 02 01, 2014, http://acrlog.org/2007/11/01/to-blog-or-not-to-blog-that-is-an-academics-question/
Bik, H., & Goldstein, M. (2013). An Introduction to Social Media for Scientists . (D. 10.1371/journal.pbio.1001535, Éd.) PLoS Biol , 11.
Boyd, D., & Ellison, N. (2007). Social network sites: Definition, history and scholarship. (D. 10.1104/pp.64.6.1070, Éd.) J Comput Mediat Commun , 13.
Brett, J. (1991). The Bureaucratization of Writing: Why so Few Academics are Public Intellectuals. Meanjin , 50, pp. 513-522.
Carmichael, M. (2010, 07 07). The Little Flaw in the Longevity-Gene Study That Could Be a Big Problem. http://www.newsweek.com/little-flaw-longevity-gene-study-could-be-big-problem-74703
Charpentier, A., Coulmont, B., & Gombin, J. (2014, 02 11). Un homme, deux voix : le vote par procuration. 02 11, 2014, La Vie des idées: http://www.laviedesidees.fr/Un-homme-deux-voix-le-vote-par.html
Chatterjee, P., & Biswas, T. (2011). Blogs and Twitter in medical publications – too unreliable to quote, or a change waiting to happen? (D. 10.7196/samj.5213., Éd.) South African Medical Journal , 101, 712-714.
Chemistry, N. (2013). All you can tweet. Nature Chemistry , 5.
Cohen, D. (2006, 08 21). Professors, Start Your Blogs. http://www.dancohen.org/2006/08/21/professors-start-your-blogs/ .
Conole, G. (2007, 10 20). The nature of academic discourse. 02 01, 2014,  http://e4innovation.com/?p=45
Conole, G. (2007, 10 29). The paper vs. blog argument…. 02 01, 2014,  http://e4innovation.com/?p=56
Cottrell, R. (2013, 02 15). Net wisdom.  02 01, 2014, http://www.ft.com/intl/cms/s/2/009050e4-75ea-11e2-9891-00144feabdc0.html
Crane, D. (1989, 02 20). How Scientists Communicate. 02 01, 2014, http://garfield.library.upenn.edu/classics1989/A1989AU43700001.pdf
Emanuel, K. (1986). Nuclear winter: Towards a scientific exercise. Nature , 319.
Esarey, J. (2013, 08 11). Blogs and Academic Tenure. 02 01, 2014,  http://politicalmethodology.wordpress.com/2013/08/11/blogs-and-academic-tenure/
Ewins, R. (2005). Who are you? Weblogs and academic identity. . (D. 10.2304/elea.2005.2.4.368, Éd.) E-Learning , 368–377.
Eysenbach, G. (2011). Can Tweets Predict Citations? Metrics of Social Impact Based on Twitter and Correlation with Traditional Metrics of Scientific Impact. J Med Internet Res , 13.
Fister, B. (2012, 07 12). Serial Scholarship: Blogging as Traditional Academic Practice. 02 01, 2014, http://www.insidehighered.com/blogs/library-babel-fish/serial-scholarship-blogging-traditional-academic-practice
Fitzpatrick, K. (2010). Planned Obsolescence. Information Standards Quarterly , 22, pp. 14-19.
Fox, J. (2012). Can blogging change how ecologists share ideas? In economics, it already has. (D. 10.4033/iee.2012.5b.15.f, Éd.) Ideas in Ecology and Evolution , 5.
Gordon, L. (1998, 09). Tenure on Trial. 02 01, 2014, http://news.stanford.edu/stanfordtoday/ed/9809/9809fea201.shtml
Graham, T. (2004). Scholarly Communication . Serials: The Journal for the Serials Community , 13, 3-11.
Gregg, M. (2009). Blogging from the Ivory Tower Hot-Desk. (D. 10.1177/1354856509342345, Éd.) Convergence , 15, 470-483.
Gregg, M. (2006). Feeling Ordinary: Blogging as conversational scholarship. (DOI: 10.1080/10304310600641604) http://espace.library.uq.edu.au/eserv.php?pid=UQ:7740&dsID=FeelingOrdinary.htm
Grollman, E. (2014, 02 04). Blogging for (a) change. 02 06, 2014, http://conditionallyaccepted.com/2014/02/04/blogging-for-a-change/.
Gu, F., & Widén-Wulff, G. (2010). Scholarly Communication and Possible Changes in the Context of Social Media: a Finnish Case Study. . The Electronic Library , 29, pp. 762-776.
Halavais, A. (2006). Scholarly blogging: Moving towards the visible college. In A. &. J.Jacobs, Uses of Blogs (pp. 117–126). New York: Peter Lang.
Halford, J. (2012, 07 19). No Academic is an Island. http://earlymoderndialogues.wordpress.com/2012/07/19/no-academic-is-an-island/
Halliday, L. (2001). Scholarly communication, scholarly publication and the status of emerging formats. Information Research. , 6.
Herndon, T. (2013, 05 22). The Grad Student Who Took Down Reinhart And Rogoff Explains Why They’re Fundamentally Wrong. 02 01, 2014, Business Insider: http://www.businessinsider.com/herndon-responds-to-reinhart-rogoff-2013-4
Honeycutt, C., & Herring, S. (2009). Beyond Microblogging: Conversation and Collaboration via Twitter. Proceedings of the Forty-Second Hawai’i International Conference on System Sciences .
Horne, D. (2011). Research as a Social Process: Considerations for Academic Libraries. Faculty of Information Quarterly , 3.
Jenkins, H. (2008, 04 08). Why Academics Should Blog….  02 01, 2014,  http://henryjenkins.org/2008/04/why_academics_should_blog.html
Jenkins, R. (2013, 08 08). What’s a Blog Post Worth?  02 01, 2014,  http://chronicle.com/blogs/onhiring/whats-a-blog-post-worth/40591
Johnson, S. (1759). The History of Rasselas, Prince of Abissinia.
Jones, K. (2013, 11 20). Should You Blog Anonymously? 02 01, 2014,  http://www.blogherald.com/2013/11/20/blog-anonymously/
Jurgenson, N. (2012, 01 23). How Academics Can Become Relevant. 02 01, 2014,  http://thesocietypages.org/cyborgology/2012/01/23/how-academics-can-become-relevant/
Kaufman, S. (2007, 11 1). An Enthusiast’s View of Academic Blogs. 02 01, 2014,  http://www.insidehighered.com/views/2007/11/01/kaufman
Kirkup, G. (2010). Academic blogging: Academic practice and academic identity. (D. 10.1080/14748460903557803, Éd.) London Review of Education , 8.
Kolowich, S. (2012, 10 02). The Academic Twitterazzi. 02 01, 2014,  http://www.insidehighered.com/news/2012/10/02/scholars-debate-etiquette-live-tweeting-academic-conferences
Kotsko, A. (2007, 11 1). A Skeptic’s Take on Academic Blogs . 02 01, 2014,  http://www.insidehighered.com/views/2007/11/01/kotsko
Kotsko, A. (2006, 05 23). On Academic Blogging: A Diagnosis. 02 01, 2014,  http://kotsko.blogspot.ca/2006/05/on-academic-blogging-diagnosis.html
Kristof, T. (2014, 02 16). Professors, We Need You!  02 16, 2014,  http://www.nytimes.com/2014/02/16/opinion/sunday/kristof-professors-we-need-you.html
Laden, G. (2008, 12 1). Anonymity & Credibility.  02 01, 2014,  http://scienceblogs.com/gregladen/2008/12/01/anonymity-credibility/
Lariv, M. (2013, 01 29). How Blogging Helped Me Write My Dissertation. 02 01, 2014,  http://chronicle.com/article/How-Blogging-Helped-Me-Write/136893/
Leo, C. (2008). Why the academic world needs blogs. 0 2 01, 2014,  http://christopherleo.com/about/why-the-academic-world-needs-blogs/
Letierce, J., Passant, A., Decker, S., & Breslin, J. (2010). Understanding how twitter is used to spread scientific messages. Web Science Conference,. Raleigh, NC.
Lupton, D. (2012, 05 22). Where are all the sociology blogs? 02 01, 2014,  http://simplysociology.wordpress.com/2012/05/22/where-are-all-the-sociology-blogs
Luzón, M. (2006). Research group blogs: sites for self-presentation and collaboration. 5th AELFE Conference. http://www.unizar.es/aelfe2006/ALEFE06/5.newtechnologies/87.pdf.
Luzón, M. (2009). Scholarly Hyperwriting:The Function of Links in Academic Weblogs. (D. 10.1002/asi.20937, Éd.) Journal of the American Society for Information Science and Technology , 60.
MacArthur, D. (2011, 07 21). Longevity genetics study retracted from Science. 02 01, 2014,  http://www.wired.com/wiredscience/2011/07/longevity-genetics-study-retracted-from-science/
Maitzen, R. (2012). Scholarship 2.0: Blogging and/as academic practice. (D. 10.1080/13555502.2012.689502, Éd.) Journal of Victorian Culture , 17.
Maron, N., & Smith, K. (2006, 11 06). Scholarly Communication.  02 01, 2014,  http://www.arl.org/sc/models/model-pubs/pubstudy/index.shtml
Marsh, A. (2013, 01 28). The boundaries of academic blogging. 02 01, 2014,  http://blogs.lse.ac.uk/impactofsocialsciences/2013/01/28/the-boundaries-of-academic-blogging/
McMillan Cottom, T. (2012, 09 30). An Idea is a Dangerous Thing to Quarantine #twittergate. 02 01, 2014, http://tressiemc.com/2012/09/30/an-idea-is-a-dangerous-thing-to-quarantine-twittergate/
Merton, R. (1942). A Note on Science and Democracy. Journal of Legal and Political Sociology , 1, 115-126.
Mewburn, I., & Thomson, P. (2013). Why do academics blog? An analysis of audiences, purposes and challenges. (D. 10.1080/03075079.2013.835624, Éd.) Studies in Higher Education , 38, 1105-1119 .
Mortensen, T., & Walker, J. (2002). Blogging Thoughts: Personal Publication as an Online Research Tool. In A. Morrison http://possibleworlds.blogs.com/blogsperiment/files/Researching_ICTs_in_context-Ch11-Mortensen-Walker.pdf (Éd.), Research ICTs in Context. University of Olso Press.
Nicolas, D., & Rowlands, I. (2011). Social Media Use in the Research Workflow. Information Sevices and Use , 31, pp. 61-83.
Noel, H. (2013, 08 09). What is tenure for? . 02 01, 2014, http://mischiefsoffaction.blogspot.ca/2013/08/what-is-tenure-for.html
O’Connor, B. (2007, 03 27). Seth Roberts and academic blogging. 02 01, 2014, http://brenocon.com/blog/2007/03/seth-roberts-and-academic-blogging/
Oreskes, N., & Conway, E. (2011). Merchants of Doubt: How a Handful of Scientists Obscured the Truth on Issues from Tobacco Smoke to Global Warming. Bloomsbury Press.
Park, H., & Thelwall, M. (2006). Web-science communication in the age of globalization. (D. 10.1177/1461444806065660, Éd.) New Media & Society , 629-650.
Parr, C. (2012, 11 1). Blog-standard turn-offs for social media neophytes. 02 1, 2014, http://www.timeshighereducation.co.uk/421669.article
Porst, S. (2003, 04 16). Why weblogs are rarely used to document research (2). 02 01, 2014, mathemagenic: http://blog.mathemagenic.com/2003/04/16.html
Price, D. (1986). Little Science, Big Science… and Beyond. Columbia University Press.
Priem, J., & Costello, K. (2010, 10 22). How and why scholars cite on Twitter.02 01, 2014, http://jasonpriem.org/self-archived/Priem_Costello_Twitter.pdf
Priem, J., Piwowar, H., & Hemminger, B. (2012, March 20). Altmetrics in the wild: Using social media to explore scholarly impact. arXiv:1203.4745.http://arxiv.org/abs/1203.4745.
Quiggin, J. (2011a). Economic Blogs. (D. 10.1111/j.1759-3441.2011.00148.x, Éd.) Economic Papers , 30, 437-440.
Quiggin, J. (2011b.) Economic Blogs and Blog Economics. In A. B. Jacobs, Uses of Blogs.
Reid, A. (2011). On the value of academic blogging . http://alex-reid.net/2011/03/on-the-value-of-academic-blogging.html .
Riesch, H., & Mendel, J. (2014). Science Blogging: Networks, Boundaries and à Limitations. (DOI:10.1080/09505431.2013.801420, Éd.) Science as Culture , 23 (5), pp. 51–72.
Risam, R. (2012, 09 30). Conference Live Tweets: Twitter Good or #Twittergate? 02 01, 2014, http://roopikarisam.com/2012/09/30/conference-live-tweets-twitter-good-or-twittergate/
Sebastiani, P., & al., e. (2010, 07 01). Genetic signatures of exceptional longevity in humans. Science .
Shema, H., Bar-Ilan, J., & Thelwall. (2012, 05 11). Research Blogs and the Discussion of Scholarly Information . (D. 10.1371/journal.pone.0035869, Éd.) PLOS One .
Shuai, X., Pepe, A., & Bollen, J. (2012). How the Scientific Community Reacts to Newly Submitted Preprints: Article Downloads, Twitter Mentions, and Citations. DOI: 10.1371/journal.pone.0047523.
Siemens, G. (2007, 10 06). Academic bloggers. 02 01, 2014, on http://www.elearnspace.org/blog/2007/10/06/academic-bloggers/
Silva, L. (2013, 05 12). So You Want to Blog (Academic Edition) . From http://www.insidehighered.com/blogs/university-venus/so-you-want-blog-academic-edition
Snow, C. (1963). The Two Cultures: a second look. . (C. U. Press., Éd.)
Tenopir, C., Volentine, R., & King, D. (2013). Social media and scholarly reading. (D. 10.1108/OIR-04-2012-0062, Éd.) Online Information Review , 193-216.
Veletsianos, G. (2012). Higher education scholars’ participation and practices on Twitter. (D. 10.1111/j.1365-2729.2011.00449.x, Éd.) Journal of Computer Assisted Learning , 28, 336-349 .
Wagers, S. (2012, 12 07). Don’t Have Time to Tweet-bollocks! Twitter can even save you time as a scientist. 02 01, 2014, http://mendelspod.com/blog/dont-have-time-to-tweet-bollocks
Walker, J. (2006). Blogging from inside the ivory tower. In A. B. Jacobs, Uses of blogs.
Walker, J. (2009). Weblog. Dans M. K.-L. David Herman, & T. &. Francis (Éd.), Routledge Encyclopedia of Narrative Theory.
Weller, K., & Puschmann, C. (2011). Twitter for scientific communication: How can citations/references be identified and measured? http://journal.webscience.org/500/1/153_​paper.pdf. (Éd.), Web Science Conference, 2011. Germany.
Wilcox, C. (2012). It’s time to e-volve: Taking responsibility for science communication in a digital age. (http://www.biolbull.org/content/222/2/85.full, Éd.) Biological Bulletin , 222, 85-87.
Wilkins, J. (2008). The roles, reasons and restrictions of science blogs. (D. 10.1016/j.tree.2008.05.004, Éd.) Trends Ecol Evol , 23.
Wilkinson, D., Harries, G., Thelwall, M., & Price, E. (2003). Motivations for academic web site interlinking: evidence for the Web as a novel source of information on informal scholarly communication . (D. 10.1177/016555150302900105, Éd.) Journal of Information Science , 29.
Ylijoki, O.-H. (2013). Boundary work between work and life in the high-speed university. (D. 10.1080/03075079.2011.577524, Éd.) Studies in Higher Education , 38, 242-255.
Ziman, J. (1996). Post-Academic Science: Constructing Knowledge With Netwroks and Norms. Science Studies , 9, 67-80.

Somewhere else, part 116

Some writings worth reading,

This fear of being unmasked as the incompetent you “really” are is so common that it actually has a clinical name: impostor syndrome. A shocking number of successful people (particularly women), believe that they haven’t really earned their spots, and are at risk of being unmasked as frauds at any moment. Many people deliberately seek out easy tests where they can shine, rather than tackling harder material that isn’t as comfortable. If they’re forced into a challenge they don’t feel prepared for, they may even engage in what psychologists call “self-handicapping”: deliberately doing things that will hamper their performance in order to give themselves an excuse for not doing well. Self-handicapping can be fairly spectacular: in one study, men deliberately chose performance-inhibiting drugs when facing a task they didn’t expect to do well on. “Instead of studying,” writes the psychologist Edward Hirt, “a student goes to a movie the night before an exam. If he performs poorly, he can attribute his failure to a lack of studying rather than to a lack of ability or intelligence. On the other hand, if he does well on the exam, he may conclude that he has exceptional ability, because he was able to perform well without studying.”

 If there were such a thing, it would probably be rule No. 1 in the teaching manual for instructors of aspiring suicide bombers: Don’t give lessons with live explosives. In what represented a cautionary tale for terrorist teachers, and a cause of dark humor for ordinary Iraqis, a commander at a secluded terrorist training camp north of Baghdad unwittingly used a belt packed with explosives while conducting a demonstration early Monday for a group of militants, killing himself and 21 other members of the Islamic State of Iraq and Syria, army and police officials said. Iraqi citizens have long been accustomed to daily attacks on public markets, mosques, funerals and even children’s soccer games, so they saw the story of the fumbling militants as a dark — and delicious — kind of poetic justice, especially coming amid a protracted surge of violence led by the terrorist group, including a rise in suicide bombings. Just last week a suicide bomber struck a popular falafel shop near the Ministry ofForeign Affairs here, killing several people. On Monday evening Raad Hashim, working the counter at a liquor store near the site of the attack, burst out laughing when he heard the news.

P values have always had critics. In their almost nine decades of existence, they have been likened to mosquitoes (annoying and impossible to swat away), the emperor’s new clothes (fraught with obvious problems that everyone ignores) and the tool of a “sterile intellectual rake” who ravishes science but leaves it with no progeny. One researcher suggested rechristening the methodology “statistical hypothesis inference testing”, presumably for the acronym it would yield. The irony is that when UK statistician Ronald Fisher introduced the P value in the 1920s, he did not mean it to be a definitive test. He intended it simply as an informal way to judge whether evidence was significant in the old-fashioned sense: worthy of a second look. The idea was to run an experiment, then see if the results were consistent with what random chance might produce. Researchers would first set up a ‘null hypothesis’ that they wanted to disprove, such as there being no correlation or no difference between two groups. Next, they would play the devil’s advocate and, assuming that this null hypothesis was in fact true, calculate the chances of getting results at least as extreme as what was actually observed. This probability was the P value. The smaller it was, suggested Fisher, the greater the likelihood that the straw-man null hypothesis was false.

  • and some thoughts about comsumption

avec un peu de lecture en français,

Did I miss something?

Ἀγεωμέτρητος μηδεὶς εἰσίτω

comme disait Platon (si on en croit la légende). J’aurais pu aussi mettre comme titre “retour sur un malentendu“. Le débat de jeudi soir, sur le Big Data, était quelque peu déstabilisant. Pour revenir un peu sur l’histoire complète (j’avais promis de faire un billet sur tout ce que j’aurais pu apprendre pendant la soirée), j’avais été contacté pour participer à un débat sur le Big Data il y a quelques semaines. J’avais alors accepté, et suggéré d’inviter aussi Yves-Alexandre, que je n’avais jamais rencontré, mais dont j’avais vu passer des études, et que je suivais sur Twitter. Après discussion avec Sophie, on s’était dit qu’il serait pertinent d’avoir des intervenants complémentaires, connaissant en particulier les réseaux, le droit des réseaux, ou pour parler de data visualization, ou de data journalism. Bref, Sophie a contacté Vincent Gautrais et Jean-Hugues Roy, ce qui permettais d’élargir le panel, et ne pas avoir uniquement le point de vue de statisticiens. J’ai découvert le titre en voyant l’affiche, et je me disais que d’un point de vue communication, associer Big Data et Big Brother était une stratégie qui devait pouvoir se justifier. Et effectivement, en quelques heures toutes les places avaient été réservées. Depuis, j’avoue avoir beaucoup travaillé, car j’étais anxieux à l’idée de débattre sur un sujet aussi technique, dans un pavillon du campus des sciences. Je voulais éviter des imprécisions sur les données utilisées en climatologie et en météorologie. Je voulais éviter d’embrouiller si on parlait d’aspects techniques de parallélisation des algorithmes. Je voulais éviter d’être confus si on parlait inférence causale (même si sur ce point, Yves-Alexandre aurait pu éclairer le débat). Mais je n’avais pas imaginer que le débat tourne autour de Big Brother, du sentiment de violation de la vie privée, et des angoisses que cela causait. Après coup, je me rends compte que j’étais probablement un usurpateur (je n’ai pas été le seul à m’en rendre compte, compte tenu de quelques réactions que nous avons eu par la suite) dans le “débat” car il aurait mieux valu avoir un psychologue, ou un sociologue, pour parler de psychose collective face aux technologies, ou du sentiment de viol (de sa vie privée). Je pense aussi qu’un spécialiste en bureautique, pour expliquer comment désactiver les options de publicité ciblée sous Google, aurait été utile puisque plusieurs remarques portait sur ce point. Cela dit, je mettais de guillemets à “débat” car beaucoup d’interventions de la salle étaient plus des commentaires que des questions destinées aux quatre “experts” sollicités (mais peut-être est-ce une manière classique de débattre dans certains lieux). J’avoue avoir été surpris par le fait que, dans un débat tenu dans un cadre universitaire, le public applaudisse un discours partisan, une forme de témoignage. Comme je l’avais dit, je n’ai pas l’habitude de participer à ce genre de débat public, et j’ai trouvé cela déconcertant à plusieurs reprises. Pour l’anecdote, j’avais l’impression en discutant avec Sophie, qui animait le débat, qu’une opposition serait faite entre “sciences dures” et “sciences molles” (mathématiciens opposés au journaliste et au juriste), mais finalement, on a surtout ressenti une animosité de la part de certaines personnes dans la salle (dont plusieurs ont pris la parole). Afin d’illustrer la malentendu évoqué au début, je me contenterais de m’auto-troller, afin d’offrir un témoignage d’un des intervenants, justement, en mettant en commentaire un message qu’il s’est permis de nous envoyer à tous les quatre vendredi après-midi.

Maintenant, il y a eu quelques commentaires (voire quelques questions pertinentes) sur lesquelles je voulais revenir. Et mentionner des points qui ont été évoqué lorsque la discussion s’est poursuivie à quatre (avec les autres intervenants, dans un vrai bar cette fois, pas “des Sciences”). Histoire d’introduire le débat, j’avais voulu évoquer un exemple qui a beaucoup fait parler ces derniers jours (je pense à Who owns real-time sports data?), et qui me semblait parfaitement introduire le débat. Un match de basket, c’est ça, par exemple, avec plusieurs joueurs qui bougent vite, des joueurs qui attaquent, d’autre qui défendent, et un ballon.

Maintenant, il faut savoir que toute l’action se retrouve sous forme de data, via des caméras qui permettent de tout numériser. Bref, on passe d’observations en direct, en trois dimensions, qui évoluent dans le temps, à des “données”. Ce que j’appellerais des données complexes. Avec en plus, une reconnaissance des joueurs, qui peuvent être suivis, même s’ils bougent (très) vite.

Si on fait un petit dessin, on peut faire une analyse spatio-temporelle de l’action de jeux (et analyser les prises de décision).

Si on en croit wikipedia, le Big Data, c’est les 3V: volume (avec des gros volumes de données – ici on part d’images vidéos, analysées, dont on extrait des informations, au dixième de seconde près), vélocité (le but est de faire des analyses en temps réel, de comprendre les décisions des joueurs, comme tirer, ou passer la balle à un coéquipier) et variété (on a toute sorte de données, effectivement, avec des localisations, des noms de joueurs, des propriétés sur les joueurs, car on a accès aux statistiques des joueurs, on sait qui est bon pour marquer à trois points, qui défend bien). Bref, j’avais envie d’illustrer le Big Data avec cet exemple qui me semble représentatif (même si un peu idyllique, car peu de monde a accès à des données aussi incroyables). Cela dit, ces données sont les observations d’individus consentants, qui semblent éloignées des préoccupations de certaines personnes qui souhaitaient que le débat porte sur la vie privée.

Parmi les points que j’ai noté, Vincent a évoqué une intervention de Michel Serres qui permet de réfléchir, effectivement,

Par la suite, on a pu parler longuement de ce qu’est l’anonymat des données (qui a été peu abordé pendant le débat proprement dit, malheureusement) en revenant sur Unique in the Crowd: The privacy bounds of human mobility, étude à laquelle a participé Yves-Alexandre, sur les métadonnées, et le fait que “human mobility traces are highly unique. In fact, in a dataset where the location of an individual is specified hourly, and with a spatial resolution equal to that given by the carrier’s antennas, four spatio-temporal points are enough to uniquely identify 95% of the individuals” (on pourra aussi relire MIT News à ce sujet)

La discussion sur les métadonnées était passionnante: les métadonnées sont suffisantes pour faire des études simples (comme l’étude sur la grippe de Google, évoquée dans un précédant billet, que j’avais écrit en vue de préparer le débat, mais on peut penser aux cartes en temps réel indiquant l’état du trafic sur les routes). Mais cette étude remet beaucoup de chose en cause, sur l’anonymat des métadonnées. Sujet passionnant, s’il en est !

Je n’ai pas noté grand chose, malheureusement, donc si des personnes ont pris des notes, avec des exemples ou des commentaires intéressants, je serais ravi de les relayer.

Data, Datum

Ce soir, je quitte ma cabane d’ermite pour intervenir dans la conférence Big Data: Big Brother organisé au Cœur des Sciences, où un peu plus de trois cent personnes se sont inscrites en quelques heures.

On va parler de data. Et depuis, je me prépare autant que possible. Comme le rappelle l’article Data is data, or are they?,

Data emerged in 1646 as the plural of the Latin datum, which according to the Oxford English Dictionary was the past participle of dare (“give”) and meant “a thing given or granted; a thing known or assumed as a fact, and made the basis of reasoning or calculation; a fixed starting point for a series of measurements etc.” Datum remains standard and retains the general meaning of “a unit of information”, though it tends to appear mostly in academic and specialist disciplines such as philosophy, surveying, geodesy, topography, technical drawing, and cartography. The meaning of the derived plural data has changed somewhat over the centuries. The OED definition from the late 19th century (“Facts, esp. numerical facts, collected together for reference or information”) seems to testify to the broadening influence of the hard sciences. In the 20th century, the rapidly expanding fields of information technology incorporated the word into a huge variety of IT- and computer-related compound nouns, such as database, data entry, data flow, data mining, data processing, data protection, and data stream. The plural data is used in many scientific, technical, academic and other formal contexts, though different practices prevail in different places. In computing jargon, social sciences, and everyday use, data is often treated as an abstract mass noun, like information. It has the general meaning “mass of information”

Ce n’est pas dans mes habitudes de parler devant autant de monde, pour une discussion ouverte selon les termes de Sophie. J’ai quand même préparer deux images pour illustrer ce qu’est le big data, en ppt. La suite ce soir (ou dans un prochain billet si j’arrive à prendre suffisamment de notes au fur et à mesure des interventions).

Temperatures Series as Random Walks

Last year, I did mention in a post that unit-root tests are dangerous, because they might lead us to strange models. For instance, in a post, I did obtain that the temperature observed in January 2013, in Montréal, might be considered as a random walk process (or at leat an integrated process). The code to extract the data has changed (since the website has been updated), so here, we use

library(RCurl)
library(XML)
options(RCurlOptions = list(useragent = "R"))
HEURE=0:23
extracttemp=function(Y,M,D){
url=paste(
"http://climate.weather.gc.ca/climateData/hourlydata_e.html?timeframe=1&Prov=QC&StationID=5415&Year=",Y,"&Month=",
M,"&Day=",D,sep="")
wp <- getURLContent(url)
doc <- htmlParse(wp, asText = TRUE) 
docName(doc) <- url
tmp <- readHTMLTable(doc)
basejour=data.frame(Year=Y,Month=M,Day=D,
Hour=HEURE,Temp=as.numeric(as.character(data.frame(tmp[2])[,2]))[2:25])
return(basejour)}
B=NULL
for(y in 1955:2013){
for(d in 1:31){
B=rbind(B,extracttemp(y,1,d))}}

Here are all the temperatures observed, and 2013,

plot(B$X,B$Temp,cex=.5,col="light blue",xlab="January, in Montreal",ylab="Temperature (Celsius)")
I=which(B$Year==2013)
lines(B$X[I],B$Temp[I],col="red")

In the previous post, one test only was used, and one year was considered. I was wondering if this behavior was observed only with temperature of 2013 (or not), and how the other tests (mentioned in a previous post too) were performing.

I might need a function, because those tests cannot be used if there is a missing value, even only one. So I did use the value observed one hour before (just to make sure that the tests can be done)

correcty=function(Y){
I=which(is.na(Y))
	if(length(I)==0){Yc=Y}
	if(length(I)>0){Yc=Y;for(i in I) Yc[i]=Yc[i-1]}	
return(Yc)
}

Now, we can compute the p-values, for all the years, and the three different three (keeping in mind that two test if the series is non-stationary, and one if the series is stationary)

DF=matrix(NA,2013-1954,3)
library(urca)
for(y in 1955:2013){
Z=B$Temp[which(B$Year==y)]
	Zc=correcty(Z)
	DF[y-1954,2]=as.numeric(pp.test(Zc)$p.value)	
	DF[y-1954,1]=as.numeric(kpss.test(Zc)$p.value)	
	DF[y-1954,3]=as.numeric(adf.test(Zc)$p.value)
}

Visually, if red means stationary, and blue means non-stationary, we get

DFP=DF
DFP[,1]=DF[,1]<.05
DFP[,2:3]=DF[,2:3]>.05
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")
plot(0:1,0:1,xlim=c(1950,2015),ylim=c(0,3),axes=FALSE,xlab="",ylab="")
axis(1)
text(1952,.5,"KPSS")
text(1952,1.5,"PP")
text(1952,2.5,"ADF")
for(y in 1955:2013){
for(i in 1:3){
polygon(y+c(-1,-1,1,1)/2.2,i-.5+c(-1,1,1,-1)/2.2,col=CL[1+(DFP[y-1954,i]==1)*5],border=NA)}}

Quite frequently, we conclude that the temperature is a random walk. Which does not make sense (from a physical point of view). But again, it might come from the fact that temperature are stationary, but with some fractional behavior (as suggested in the previous post).

Unit Root Tests

This week, in the MAT8181 Time Series course, we’ve discussed unit root tests. According to Wold’s theorem, if https://latex.codecogs.com/gif.latex?(Y_t) is  (weakly) stationnary then

https://latex.codecogs.com/gif.latex?%20%20%20%20%20Y_{t}=\sum%20_{{j=0}}^{\infty%20}\psi_{j}\varepsilon%20_{{t-j}}+\xi%20_{t}

where https://latex.codecogs.com/gif.latex?%20(\varepsilon%20_{{t}}) is the innovation process, and where https://latex.codecogs.com/gif.latex?%20(\xi%20_{{t}}) is some deterministic series (just to get a result as general as possible). Observe that

https://latex.codecogs.com/gif.latex?\sum%20_{{j=0}}^{{\infty%20}}|\psi_{{j}}|^{2}%20%3C%20\infty

as discussed in a previous post. To go one step further, there is also the Beveridge-Nelson decomposition : an integrated of order one process, defined as

https://latex.codecogs.com/gif.latex?%20%20%20%20%20\Delta%20Y_{t}=(1-L)%20Y_t=\sum%20_{{j=0}}^{\infty%20}\psi_{j}\varepsilon%20_{{t-j}}+\xi=\Psi(L)\varepsilon%20_{{t}}+\xican be represented as

a linear trend https://latex.codecogs.com/gif.latex?+ a random walk https://latex.codecogs.com/gif.latex?+ a stationary remaining term

i.e.

https://latex.codecogs.com/gif.latex?%20Y_{t}=\underbrace{Y_0%20+%20\xi%20t%20%C2%A0}+\underbrace{\Psi(1)\sum_{j=1}^t\varepsilon%20_{{i}}}+\underbrace{\tilde\Psi(L)\varepilon_0-\tilde\Psi(L)\varepsilon_t}

where https://latex.codecogs.com/gif.latex?%20\tilde\Psi(\cdot) is the polynomial with terms https://latex.codecogs.com/gif.latex?%20\tilde\psi_j, where

https://latex.codecogs.com/gif.latex?%20\tilde\psi_j%20=\sum_{i=j+1}^\infty\psi_i

For unit-root tests, we will use various representation of the process. In order to illustrate the implementation of those tests, consider the following series

> E=rnorm(240)
> X=cumsum(E)
> plot(X,type="l")
  • Dickey Fuller (standard)

Here, for the simple version of the Dickey-Fuller test, we assume that

https://latex.codecogs.com/gif.latex?%20Y_t=\alpha+\beta%20t+\varphi%20Y_{t-1}+\varepsilon_t

and we would like to test if https://latex.codecogs.com/gif.latex?%20\varphi=1 (or not). We can write the previous representation as

https://latex.codecogs.com/gif.latex?%20\Delta%20Y_t=\alpha+\beta%20t+[\varphi-1]%20Y_{t-1}+\varepsilon_t

so we simply have to test if the regression coefficient in the linear regression is – or not – null. Which can be done with Student’s test. If we consider the previous model without the linear drift, we have to consider the following regression

> lags=0
> z=diff(X)
> n=length(z)
> z.diff=embed(z, lags+1)[,1]
> z.lag.1=X[(lags+1):n]
> summary(lm(z.diff~0+z.lag.1 ))

Call:
lm(formula = z.diff ~ 0 + z.lag.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.84466 -0.55723 -0.00494  0.63816  2.54352 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.005609   0.007319  -0.766    0.444

Residual standard error: 0.963 on 238 degrees of freedom
Multiple R-squared:  0.002461,	Adjusted R-squared:  -0.00173 
F-statistic: 0.5873 on 1 and 238 DF,  p-value: 0.4442

Our testing procedure will be based on the Student’s t value,

> summary(lm(z.diff~0+z.lag.1 ))$coefficients[1,3]
[1] -0.7663308

which is exactly the value computed using

> library(urca)
> df=ur.df(X,type="none",lags=0)
> df

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: -0.7663

The interpretation of this value can be done using critical values (99%, 95%, 90%)

> qnorm(c(.01,.05,.1)/2)
[1] -2.575829 -1.959964 -1.644854

If the statistics exceeds those values, then the series is not stationnary, since we cannot reject the assumption that https://latex.codecogs.com/gif.latex?%20\varphi-1=0. So we might conclude that there is a unit root. Actually, those critical values are obtained using

> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 

Call:
lm(formula = z.diff ~ z.lag.1 - 1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.84466 -0.55723 -0.00494  0.63816  2.54352 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.005609   0.007319  -0.766    0.444

Residual standard error: 0.963 on 238 degrees of freedom
Multiple R-squared:  0.002461,	Adjusted R-squared:  -0.00173 
F-statistic: 0.5873 on 1 and 238 DF,  p-value: 0.4442

Value of test-statistic is: -0.7663 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

The problem with R is that there are several packages that can be used for unit root tests. Just to mention another one,

> library(tseries)
> adf.test(X,k=0)

	Augmented Dickey-Fuller Test

data:  X
Dickey-Fuller = -2.0433, Lag order = 0, p-value = 0.5576
alternative hypothesis: stationary

We do have here also a test where the null hypothesis is that there is a unit root. But the p-value is quite different. What is odd is that we have

> 1-adf.test(X,k=0)$p.value
[1] 0.4423705
> df@testreg$coefficients[4]
[1] 0.4442389

(but I think it is a coincidence).

  • Augmented Dickey Fuller

It is possible to had some lags in the regression. For instance, we can consider

https://latex.codecogs.com/gif.latex?%20\Delta%20Y_t=\alpha+\beta%20t+[\varphi-1]%20Y_{t-1}+\psi%20\Delta%20Y_{t-1}+\varepsilon_t

Again, we have to check if one coefficient is null, or not. And this can be done using Student’s t test.

> lags=1
> z=diff(X)
> n=length(z)
> z.diff=embed(z, lags+1)[,1]
> z.lag.1=X[(lags+1):n]
> k=lags+1
> z.diff.lag = embed(z, lags+1)[, 2:k]
> summary(lm(z.diff~0+z.lag.1+z.diff.lag ))

Call:
lm(formula = z.diff ~ 0 + z.lag.1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87492 -0.53977 -0.00688  0.64481  2.47556 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657

Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared:  0.003292,	Adjusted R-squared:  -0.005155 
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777

> summary(lm(z.diff~0+z.lag.1+z.diff.lag ))$coefficients[1,3]
[1] -0.7328138

This value is the one obtained using

> df=ur.df(X,type="none",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 

Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87492 -0.53977 -0.00688  0.64481  2.47556 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657

Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared:  0.003292,	Adjusted R-squared:  -0.005155 
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777

Value of test-statistic is: -0.7328 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

And again, other pckages can be used:

> adf.test(X,k=1)

	Augmented Dickey-Fuller Test

data:  X
Dickey-Fuller = -1.9828, Lag order = 1, p-value = 0.5831
alternative hypothesis: stationary

Hopefully, the conclusion is the same (we should reject the assumption that the series is stationary, but I am not sure about the computation of the p-value).

  • Augmented Dickey Fuller with trend and drift

So far, we have not included the drift in our model. But this is simple to do (this will be called the augmented version of the previous procedure): we just have to include a constant in the regression,

> summary(lm(z.diff~1+z.lag.1+z.diff.lag ))

Call:
lm(formula = z.diff ~ 1 + z.lag.1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.91930 -0.56731 -0.00548  0.62932  2.45178 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared:  0.02313,	Adjusted R-squared:  0.01482 
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393

The statistics of interest are obtained here considering some analysis of variance outputs, where this model is compared with the one without the integrated part, and the drift,

> summary(lm(z.diff~1+z.lag.1+z.diff.lag ))$coefficients[2,3]
[1] -2.303948
> anova(lm(z.diff ~ z.lag.1 + 1 + z.diff.lag),lm(z.diff ~ 0 + z.diff.lag))$F[2]
[1] 2.732912

Those two values are the ones obtained also with

> df=ur.df(X,type="drift",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.91930 -0.56731 -0.00548  0.62932  2.45178 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared:  0.02313,	Adjusted R-squared:  0.01482 
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393

Value of test-statistic is: -2.3039 2.7329 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

And we can also include a linear trend,

> temps=(lags+1):n
> summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))

Call:
lm(formula = z.diff ~ 1 + temps + z.lag.1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87727 -0.58802 -0.00175  0.60359  2.47789 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
temps       -0.0004194  0.0009767  -0.429   0.6680  
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared:  0.0239,	Adjusted R-squared:  0.01139 
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287

> summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))$coefficients[3,3]
[1] -1.98282
> anova(lm(z.diff ~ z.lag.1 + 1 + temps+ z.diff.lag),lm(z.diff ~ 1+ z.diff.lag))$F[2]
[1] 2.737086

while R function returns

> df=ur.df(X,type="trend",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87727 -0.58802 -0.00175  0.60359  2.47789 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
tt          -0.0004194  0.0009767  -0.429   0.6680  
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared:  0.0239,	Adjusted R-squared:  0.01139 
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287

Value of test-statistic is: -1.9828 1.8771 2.7371 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47
  • KPSS test

Here, in the KPSS testing procedure, two models can be considerd : with a drift, or with a linear trend. Here, the null hypothesis is that the series is stationnary.

With a drift, the code is

> summary(ur.kpss(X,type="mu"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 4 lags. 

Value of test-statistic is: 0.972 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.73

while it will be, in the case there is a trend

> summary(ur.kpss(X,type="tau"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: tau with 4 lags. 

Value of test-statistic is: 0.5057 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

One more time, it is possible to use another package to get the same test (but again, a different output)

> kpss.test(X,"Level")

	KPSS Test for Level Stationarity

data:  X
KPSS Level = 1.1997, Truncation lag parameter = 3, p-value = 0.01

> kpss.test(X,"Trend")

	KPSS Test for Trend Stationarity

data:  X
KPSS Trend = 0.6234, Truncation lag parameter = 3, p-value = 0.01

At least, there is some kind of consistency, since we keep rejecting the stationnary assumption, for that series.

  • Philipps-Perron test

The Philipps-Perron test is based on the ADF procedure. The code is here

> PP.test(X)

	Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller = -2.0116, Truncation lag parameter = 4, p-value = 0.571

with again, a possible alternative with the other package

> pp.test(X)

	Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller Z(alpha) = -7.7345, Truncation lag parameter = 4, p-value
= 0.6757
alternative hypothesis: stationary
  •  Comparison

I will not spend more time comparing the different codes, in R, to run those tests. Let us spend some additional time on a quick comparison of those three procedure. Let us generate some autoregressive processes, with more or less autocorrelation, as well as some random walk, and let us see how those tests perform :

> n=100
> AR=seq(1,.7,by=-.01)
> P=matrix(NA,3,31)
> M1=matrix(NA,1000,length(AR))
> M2=matrix(NA,1000,length(AR))
> M3=matrix(NA,1000,length(AR))

> for(i in 1:(length(AR)+1)){
+ for(s in 1:1000){
+ if(i==1) X=cumsum(rnorm(n))
+ if(i!=1) X=arima.sim(n=n,list(ar=AR[i]))
+ library(urca)
+ M2[s,i]=as.numeric(pp.test(X)$p.value)
+ M1[s,i]=as.numeric(kpss.test(X)$p.value)
+ M3[s,i]=as.numeric(adf.test(X)$p.value)
+ }}

Here, we would like to count how many times the p-value of our tests exceed 5%,

> prop05=function(x) mean(x>.05)
+ P[1,]=1-apply(M1,2,prop05)
+ P[2,]=apply(M2,2,prop05)
+ P[3,]=apply(M3,2,prop05)
+ }
> plot(AR,P[1,],type="l",col="red",ylim=c(0,1),ylab="proportion of non-stationnary 
+ series",xlab="autocorrelation coefficient")
> lines(AR,P[2,],type="l",col="blue")
> lines(AR,P[3,],type="l",col="green")
> legend(.7,1,c("ADF","KPSS","PP"),col=c("green","red","blue"),lty=1,lwd=1)

 

We can see here how poorly Dickey-Fuller test behave, since a 50% (at least) of our autoregressive processes are considered as non-stationnary.

Le Vote par Procuration en France

La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“.

Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition.

  • Modèle de régression, et recherche de variables explications

Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que  où  est le nombre de procurations dans le bureau de vote , et où  est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que , la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

Si on suppose que  est fonction  le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote),

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference]

il est légitime de tenter une régression logistique,

voire un lissage par splines, si on pense que le lien peut ne pas être linéaire,

Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques)

> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])

et, pour la régression logistique standard (linéaire)

> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration).

On peut faire la même chose sur le taux de logement HLM, dans le quartier,

Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de  laisse à désirer (toute suggestion est la bienvenue !).

> summary(regt1l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset1, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.9549   -1.5722    0.0319    1.6292   13.1303  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.70811    0.01516  -244.6   <2e-16 ***
taux         1.49666    0.04012    37.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12507  on 836  degrees of freedom
Residual deviance: 11065  on 835  degrees of freedom
AIC: 15699

Number of Fisher Scoring iterations: 4

> summary(regt2l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset2, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-15.4872   -1.7817   -0.1615    1.6035   12.5596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24272    0.01230 -263.61   <2e-16 ***
taux         1.45816    0.03266   44.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9424.7  on 836  degrees of freedom
Residual deviance: 7362.3  on 835  degrees of freedom
AIC: 12531

Number of Fisher Scoring iterations: 4

 

On a alors voulu comparer avec un modèle qui me semble moins juste, mais qui est plus simple à interpréter, où on suppose que le taux de procuration (par bureau) est expliqué par un modèle linéaire

(que l’on peut aussi lisser, pour vérifier que le lien est effectivement linéaire). Le code est ici

> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre)

Cette fois, on a une sortie de régression plus classique,

> summary(regt1l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.9994 -0.2926  0.0011  0.3173  3.2072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.021268   0.001739   12.23   <2e-16 ***
taux        0.054371   0.004812   11.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared:  0.1326,	Adjusted R-squared:  0.1316 
F-statistic: 127.7 on 1 and 835 DF,  p-value: < 2.2e-16

> summary(regt2l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.9029 -0.4148 -0.0338  0.4029  3.4907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.033909   0.001866   18.17   <2e-16 ***
taux        0.079749   0.005165   15.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared:  0.2221,	Adjusted R-squared:  0.2212 
F-statistic: 238.4 on 1 and 835 DF,  p-value: < 2.2e-16

On note plusieurs choses de ces graphiques. (i) les deux types de régression donnent des modèles pour les taux de procuration très très proches. Donc autant prendre le plus simple à interpréter. (ii) le lissage n’apporte rien, et le modèle linéaire semble pertinent. On a ainsi regardé plusieurs variables, et on en a retenu un certain nombre, pour tenter un modèle multiple. Avant de parler des résidus de notre modèle, je devrais peut être prendre quelques lignes pour parler un peu de cartographie (une nouvelle fois, Baptiste m’a fait découvrir de belles fonctions).

  • Visualiser les bureaux de vote à Paris

Pour récupérer le fond de carte, avec les bureaux de vote, on utilise (là encore, je renvoie au blog de Baptiste, qui explique l’utilisation de données de cartelec.net)

> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")

Si on veut visualiser, par exemple, le taux de procuration (disons la moyenne entre les deux tours), on utilise les données suivantes

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur,

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

Reste à conclure, en faisant une visualisation graphique de nos données

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

Dans ce cas, la visualisation des résidus donne

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux.

Somewhere else, part 115

Some writings worth readings (with a lot of posts on big data, and data science)

“The theme of computer science and its applications is focused on the various aspects of computer science and its applications for advances in computer science and its applications and provides an opportunity for academic and industry professionals to discuss the latest issues and progress in the area of computer science and its applications. Therefore this book will be include the various theories and practical applications in computer science and its applications.”

with temperature overlay http://earth.nullschool.net/#current/wind/…, HT @RoseGWhite

et un peu de lecture en français,


Did I miss something?

La Loi des Petits Nombres

Comme nous l’avions vu dans Charpentier (2010), la loi des grands nombres est souvent évoquée pour justifier la mutualisation des risques indépendants : plus la mutualité sera grande, plus petite sera la variabilité. A condition que les risques ne soient pas trop grands. Un cas problématique sont les risques catastrophiques : ils sont rares, et parfois tellement (potentiellement) coûteux que l’hypothèse d’existence de la variance doit être remise en cause. La modélisation de ces événements rares repose sur la loi des petits nombres, pour reprendre le nom de l’ouvrage de Ladislaus Bortkiewicz. Comme nous allons le voir, la loi de Poisson est à la loi des petits nombres ce que la loi normale est à la loi des grands nombres. Nous verrons donc en détails l’importance de la loi de Poisson pour modéliser les événements rares. Et nous verrons pourquoi la probabilité qu’un événement ne survienne pas est toujours de 37%. Ou presque.

  • La loi de Siméon-Denis Poisson

Siméon-Denis Poisson a travaillé sur les calculs des probabilités pendant presque 20 ans, de 1820 à 1840. Comme pour beaucoup de ses contemporains, ses premiers travaux portèrent sur des problèmes de jeux, avec une communication à l’Académie des Sciences sur l’avantage du banquier au jeu de trente et quarante. Mais son premier travail conséquent a porté sur un problème longuement étudié par Laplace (qui a été un de ses professeurs) sur la proportion des naissances des filles et des garçons (problème classique de statistique) publié en 1830. Il y présente en particulier une démonstration de la loi des grands nombres pour une loi de Bernoulli, qu’il modifiera par la suite. Dans ce mémoire, il y présente une loi, qui sera appelée plus tard la loi de Poisson : « la probabilité qu’une événement dont la chance à chaque épreuve est la fraction très petite x/mu n’arrive pas plus de  fois dans un très grand nombre d’épreuve mu d’épreuve (pour reprendre la terminologie utilisée en 1837, ce qui correspond à la fonction de répartition avec une terminologie plus contemporaine) est »

https://latex.codecogs.com/gif.latex?P=\left(1+x+\frac{x^2}{1\cdot%202}+\frac{x^3}{1\cdot%202\cdot%203}+\cdots+\frac{x^n}{1\cdot%202\cdots%20n}\right)e^{-x%20}

On verra réapparaître cette loi dans son fameux traité, paru en 1837, Recherches sur la probabilité des jugements. En particulier, dans le chapitre 8, il obtient sa loi comme limite de la loi binomiale https://latex.codecogs.com/gif.latex?B(T,\lambda%20/T) lorsque https://latex.codecogs.com/gif.latex?T devient grand, mais passe assez rapidement à d’autres considérations. Il n’étudie pas cette loi limite, et ne propose pas vraiment de l’utiliser. Il faudra attendre les travaux de Ladislaus Bortkiewicz, presque un siècle plus tard, pour voir des applications, dans Das Gesetz der kleinen Zahlen (la loi des petits nombres). Et comme toujours en mathématiques, lui attribuer tout le crédit est un peu excessif, puisqu’un siècle auparavant, en 1718, de Moivre avait obtenu la loi même loi, toujours comme limite de la loi Binomiale. Si le nom de Poisson est resté à la postérité, c’est essentiellement parce que Boltzmann le cite, en 1868, ainsi que Seidel en 1876, et surtout Tchebychev. Cela dit, Poisson était loin d’être un inconnu, en tant que scientifique. Ses travaux l’ont amené à travailler sur des problèmes d’électrostatique (la fameuse équation de Poisson), à l’équation de la chaleur avec Fourier, il s’est opposé à Fresnel sur des problèmes d’optique, et il a présidé à deux reprises l’Académie des Sciences.

La distribution d’un nombre d’événements obtenu comme somme de variables de Bernoulli, dans une grande population, peut s’approcher par la loi de Poisson. Ou lorsque la probabilité de survenance est faible, par rapport à la taille de l’échantillon. Formellement, si https://latex.codecogs.com/gif.latex?N suit une loi binomiale https://latex.codecogs.com/gif.latex?B(n,p), avec , alors

https://latex.codecogs.com/gif.latex?P[N=k]=\binom{n}{k}p^k(1-p)^{n-k}\sim%20e^{-\lambda}\frac{\lambda^k}{k!}

On parlera de loi des petits nombres car on compte ici des évènements rares (la probabilité de survenance d’un évènement étant inversement proportionnelle à ).

Ce résultat se traduit de la manière suivante: si on se donne un échiquier, si on lance 100 pièces et si on compte le nombre de pièces par case, la distribution suivra une loi de Poisson de moyenne 1 (le ratio entre le nombre de pièces et de cases). Une illustration est évoquée sur la Figure suivante

Nombre de pièces par case

Fréquence

Loi de Poisson

0

36

36,78

1

39

36,78

2

16

18,39

3

7

6,13

4

2

1,53

5 et plus

0

0,37

  • L’utilisation de la loi de Poisson, et les mélanges Poissonniens

A la fin du XIXème siècle, à l’université de Göttingen, Wilhelm Lexis (connu en démographie pour les diagrammes qui portent aujourd’hui son nom) avait envie d’étudier les applications de cette loi. Alors qu’il étudiait l’utilisation de la loi de Poisson dans un contexte démographique, un de ces étudiants, Ladislaus Bortkiewicz, proposa de l’utiliser pour modéliser des nombres d’accidents. Dans l’exemple fameux de Bortkiewicz, datant de 1898, il avait étudié le nombre de soldats de cavaliers morts par ruade de cheval, entre 1875 et 1894, dans 10 corps (soit 200 corps annuels). Il avait obtenu la distribution donnée dans la table suivante. La colonne de droite donne le nombre de décès donné par une loi de Poisson de même moyenne. On note que l’ajustement est très bon, et on comprend vite la fascination de Bortkiewicz pour cette loi.

Nombre de décès, par corps

Fréquence

Loi de Poisson

0

109

108,67

1

65

66,21

2

22

20,22

3

3

4,11

4

1

0,63

5 et plus

0

0,08

Maintenant, pour être tout à fait honnête, la loi que nous avons énoncée comme une loi des petits nombres est très différente de celle proposée par Bortkiewicz (celle qu’il a énoncé dans son ouvrage lui a valu un cinglant commentaire de Corrado Gini, qui écrivait en 1907, de manière provocatrice que « the law of small numbers does not exist »). Mais l’utilisation de la loi de Poisson pour modéliser les accidents, et plus généralement les petits nombres venait de débuter.

Dans les années 1930, Filip Lunderg avait noté l’intérêt (théorique) du processus de Poisson pour modéliser l’arrivée des sinistres, ainsi que plusieurs actuaires de l’école scandinave (Esscher en 1932, Sgerdahl en 1939, Lüders en 1934, etc). Depuis, tous les actuaires ont utilisé cette loi pour modéliser toute sorte d’évènements “rares“, comme la survenance annuelle d’ouragans, aux États-Unis (Figure ci-dessous).

Nombre d’ouragans

Fréquence

Loi de Poisson

0

30

27,16

1

48

47,99

2

37

42,41

3

29

24,98

4

8

11,03

5

3

3,90

6

3

1,15

7

1

0,29

8 et plus

0

0,08

  • De la loi de Poisson à la période de retour

Un concept fondamental en gestion des risques extrêmes a été introduit par Emil Gumbel en 1958, liant le temps qui s’écoule (en années) entre deux événements consécutifs et la probabilité (annuelle) de survenance. Pour des événements se produisant avec une probabilité annuelle , indépendamment les uns des autres, le temps moyen d’attente entre deux événements est , appelé période de retour, et la probabilité qu’aucun événement ne survienne pendant années (consécutives) est alors

On peut résumer ceci dans le tableau ci-dessous, où on étudie probabilité qu’un événement ne survienne pas pendant n années (en ligne) en fonction de la période de retour . Pour un événement centenaire, il y a 36,60% chance pour qu’il ne survienne pas en cent ans.

Période de retour (en années)

Nombre d’années
sans catastrophe

10

20

50

100

200

10

34,86%

59,87%

90,43%

90,43%

95,11%

20

12,15%

35,84%

66,76%

81,79%

90,46%

50

0,51%

7,69%

36,41%

60,50%

77,83%

100

0,00%

0,59%

13,26%

36,60%

60,57%

200

0,00%

0,00%

1,75%

13,39%

36,69%

  • La probabilité qu’un événement ne survienne pas est… 37%

Comme on le voit dans le tableau ci-dessus, sur la diagonale, la probabilité qu’un événement ne survienne pas pendant années quand sa période de retour est est de l’ordre de 37%. D’ailleurs, si on revient un instant sur notre échiquier évoqué auparavant, la probabilité de n’avoir aucune pièce sur une case était de… 37%

Reprenons ici un exemple qui avait fait polémique il y a quelques années, sur le risque nucléaire. Dans un article (Laponche et Dessus, 2011), on apprenait que la probabilité d’avoir un incident majeur sur réacteur nucléaire, pour une année, était de l’ordre de 0,0003 (3 chances sur 10,000). Or comme il y a 143 réacteurs nucléaires, sur 30, la probabilité d’avoir un incident majeur est (selon le calcul des auteurs)

D’où la conclusion savoureuse, « la probabilité d’occurrence d’un accident majeur sur ces parcs serait donc […] de plus de 100% pour l’Union européenne ». Comme on l’a noté, la probabilité d’avoir un incident majeur sur une période lorsque la probabilité annuelle est https://latex.codecogs.com/gif.latex?p s’écrit

https://latex.codecogs.com/gif.latex?P[N\leq%20n]=1-(1-p)^n%20\sim%20np

en utilisant un développement limité non justifié ici ! Il convient ici d’utiliser le modèle de Poisson. La probabilité d’avoir au moins un incident majeur est

https://latex.codecogs.com/gif.latex?P[N\leq%20n]=1-(1-p)^n%20\sim%201-e^{-np}

soit ici 72,47%. Si la probabilité était traduite en terme de durée

https://latex.codecogs.com/gif.latex?p=\frac{1}{T}

avec

https://latex.codecogs.com/gif.latex?T=\frac{1}{0,0003\times143}\sim%2023,31

on voit qu’ici, le temps moyen d’attente entre deux incidents majeurs en Europe est de 23 ans. Ce qui fait que sur 23 ans, la probabilité de n’avoir aucun incident majeur est de l’ordre de 37%.

  • Conclusion

La loi de Poisson est présente partout en assurance, car c’est la loi centrale pour modéliser le comptage d’événements rares. Le nombre de décès dans un portefeuille d’assurance-vie suit une loi de Poisson, tout comme le nombre d’accidents par unité de temps en assurance automobile. On retrouve même l’estimateur Chain Ladder du montant de provisions pour sinistres a payer quand on utilise une régression de Poisson sur les incréments de paiements.