Pour reprendre une question posée au téléphone par François, et par Eric dans un commentaire (ici), il est possible de comparer l’évolution du nombre de catastrophes et de l’évolution de la population (mondiale). Pour les données de population on peut en trouver sur wikipedia (ici), sur le site des nations unies (là) après 1950, ou encore sur le site du census.gov (là). Ensuite on peut moyenner les doublons.
POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt", sep=";",header=TRUE) population=tapply(POP$POP,POP$YEAR,mean) years=as.numeric(names(population)) reg=lm(population~bs(years,3)) lines(years,predict(reg),col="red") plot(years,population) baseline=predict(reg,newdata=data.frame(years=1900:2010))
Revenons à nos catastrophes “naturelles“. Par exemple les tremblements de terre. La première piste (évoquée là) était de lisser cette série pour l’utiliser en baseline (multiplicatif). On retrouve ici une série qui oscille autour de 1.
earthquake=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet.csv", sep=",",header=TRUE) flood=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_2_.csv", sep=",",header=TRUE) storm=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_3_.csv", sep=",",header=TRUE) library(splines) n=nrow(earthquake) earthquake=earthquake[-((n-3):n),] DATE=1900:2010 EARTHQUAKE=FLOOD=STORM=rep(0,length(DATE)) EARTHQUAKE[as.numeric(as.character(earthquake$Year))-1899]= as.numeric(as.character(earthquake$Number)) plot(DATE,EARTHQUAKE,type="b") reg=lm(EARTHQUAKE~bs(DATE,6)) lines(DATE,predict(reg),col="red") plot(DATE,EARTHQUAKE/predict(reg),type="b") abline(h=1,lty=2,col="blue")
Si on utilise cette fois la population mondiale comme baseline, on obtient
POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt", sep=";",header=TRUE) population=tapply(POP$POP,POP$YEAR,mean) years=as.numeric(names(population)) reg=lm(population~bs(years,3)) lines(years,predict(reg),col="red") plot(years,population) baseline=predict(reg,newdata= data.frame(years=1900:2010))/1000000 plot(DATE,EARTHQUAKE/baseline,type="b")
La série de points présente une jolie rupture. Afin de trouver la date optimale de rupture, utilisons un test basé sur le critère d’Akaike.
Y=EARTHQUAKE/baseline X=DATE DUMMY=1910:2000 VAIC=rep(NA,length(DUMMY)) for(x in DUMMY){ IX=(X>x) reg=lm(Y~1+IX) VAIC[x-1909]=AIC(reg) } plot(DUMMY,VAIC) k=which.min(VAIC) DUMMY[k]
Manifestement, il y a une rupture en 1975 (que je ne sais trop comment interpréter).
Si on regarde une autre série, comme celle des inondations.
n=nrow(flood) flood=flood[-((n-3):n),] DATE=1900:2010 FLOOD[as.numeric(as.character(flood$Year))-1899]= as.numeric(as.character(flood$Number)) plot(DATE,FLOOD,type="b") reg=lm(FLOOD~bs(DATE,6)) lines(DATE,predict(reg),col="red") population=tapply(POP$POP,POP$YEAR, mean) years=as.numeric(names(population)) reg=lm(population~bs(years,3)) baseline=predict(reg,newdata=data.frame(years=1900:2010)) plot(DATE,FLOOD/baseline,type="b")
On retrouve ici clairement le fait que la croissance de la population ne peut expliquer la croissance du nombre d’inondations.