# Prédire le nombre de morts, suite

Il y a quelques jours, un tweet de Baptiste Coulmont m’inspirait un billet sur l’estimation de la surmortalité, en France. Ce soir, c’est un tweet d’Alexandre Blanchet, qui parlait de surmortalité, au Québec. J’étais surpris, parce que la prévision est bien au dessus des années passés contrairement à ce qu’on pouvait observer en France. En repartant des codes d’Alexandre (data_prep.R puis predictions.R) j’ai remis un peu en forme le graphique, pour avoir la vision de long-terme, pour comprendre la tendance linéaire utilisée par Alexandre (qui utilise un modèle à la Buys-Ballot) Je me disais qu’on pouvait tenter des modèles de lissages exponentiel, à la Holt-Winters,

library(forecast) x=ts(df.l_training$deces,start = 2010, frequency = 12) x.hw &lt;- ets(x, model = "AAA") autoplot(forecast(x.hw,24)) en version additive, ou multiplicative x.hw2 &lt;- ets(x, model = "MAM") autoplot(forecast(x.hw2,24)) Je laisse à d’autre le soin de comparer aux nombres de morts observés depuis 2 mois, je reste fasciné pour ma part par la tendance à la hausse de la mortalité, depuis 10 ans… Et j’essayerais de comparer ces trois modèles dès que j’aurais finis mes corrections de copies de la session d’hiver. # Prédire l’espérance de vie (à la naissance) Pour finir les exemples de modèles de prévision, on a parlé de calculs d’espérance de vie (à la naissance). Le point de départ des modèles démographiques, ce sont les tables de mortalité, dont je parlais dans le précédant billet. Mais cette vision est très biaisée, car elle suppose qu’il n’y aura pas d’amélioration des conditions de vie. Pour faire les choses proprement, on utilise des données plus complètes, avec des nombres de décès en fonction de l’âge $x$, mais aussi de la date $t$. url="http://freakonometrics.free.fr/FranceDeaths_1x1.txt" download.file(url,"FRD.txt") url="http://freakonometrics.free.fr/FranceExposures_1x1.txt" download.file(url,"FRE.txt") DECES=read.table("FRD.txt",skip = 3,header=TRUE) EXPOS=read.table("FRE.txt",skip = 3,header=TRUE) On notera $D_{x,t}$ le nombre de décès observé, et $E_{x,t}$ l’exposition. De telle sorte que la probabilité de décéder dans l’année, pour une personne d’âge $x$ à la date $t$ est $q_{x,t}=D_{x,t}/E_{x,t}$. Classiquement, ces données sont stockées dans une matrice pour la visualisation, et dans une base de données pour les régressions annees = 1816:2015 ages = 0:110 MDF=matrix(DECES$Female,length(ages),length(annees)) MDH=matrix(DECES$Male,length(ages),length(annees)) MEF=matrix(EXPOS$Female,length(ages),length(annees)) MEH=matrix(EXPOS$Male,length(ages),length(annees)) MEF[MEF&lt;=0]=1 MEH[MEH&lt;=0]=1 QF=MDF/MEF QH=MDH/MEH QF[QF&lt;=0]=NA QH[QH&lt;=0]=NA (il faut bricoler un peu pour éviter les problèmes de valeurs nulles car (i) on fait des ratios (ii) on passera ensuite au logarithme). On peut visualiser cette surface, en fonction de $x$ et $t$. persp(log(QF)) ou (plus visuel) library(rgl) persp3d(ages,annees,log(QH),col="light blue") Pour modèliser l’évolution de $q_{x,t}$ on peut s’inspirer du modèle de Lee & Carter (1992) qui suppose que $\log(q_{x,t})=A_x+B_x\cdot K_t$. Les $\mathbf{A}=(A_0,A_1,\cdots,A_{110})$ c’est en quelques sortes l’évolution moyenne de $\log(q_{x,t})$. $\mathbf{K}=(K_{1816},K_{1817},\cdots,K_{2015})$ va permettre de comprendre les gains de conditions de vie, la baisse de la probabilité de décéder dans l’année. Ces gains ne sont pas uniformes, et on utilise $\mathbf{B}=(B_0,B_1,\cdots,B_{110})$ pour faire dépendre les gains de l’âge. Pour estimer les paramères, $\mathbf{A}$, $\mathbf{B}$ et $\mathbf{K}$, on peut être tentés d’utiliser un modèle binomial, au sens où $D_{x,t}\sim \mathcal{B}(E_{x,t},q_{x,t})$, qui est le modèle de base en assurance-vie. Autrement dit, ici $D_{x,t}\sim \mathcal{B}(E_{x,t},\exp[A_x+B_x\cdot K_t])$. Une autre piste est d’utiliser la loi des petits nombres, qui dit que si la probabilité est faible (et c’est le cas quand on s’intéresse à la probabilité de mourir dans l’année), la loi binomiale peut être approchée par une loi de Poisson. Autrement dit, $D_{x,t}\sim \mathcal{P}(E_{x,t}\cdot q_{x,t})$, i.e. $D_{x,t}\sim \mathcal{P}(\cdot\exp[\log(E_{x,t})+A_x+B_x\cdot K_t])$. On reconnait ici une régression de Poisson, avec comme variables explicatives l’âge $x$, l’année $t$, et l’exposition apparait comme variable offset. Le seul soucis est que ce n’est pas une régression linéaire. On a ici un modèle non linéaire, puisque $\mathbb{E}[D_{x,t}]=(\cdot\exp[\log(E_{x,t})+A_x+B_x\cdot K_t])$. On a à la fois une somme et un produit. Bref, on va partir vers un modèle nonlinéaire généralisé. base=data.frame(DH=DECES$Male,DF=DECES$Female, EH=EXPOS$Male,EF=EXPOS$Female, age=rep(ages,length(annees)), annee=rep(annees,each=length(ages))) idx=which((base$DH&lt;=0)|(base$DF&lt;=0)|(base$EH&lt;=0)|(base$EF&lt;=0)) base=base[-idx,] library(gnm) reg=gnm( DH ~ offset(log(EH)) + as.factor(age) + Mult(as.factor(age),as.factor(annee)), family = poisson(link="log"), data=base) On a nos estimateurs $\hat{\mathbf{A}}$, $\hat{\mathbf{B}}$ et $\hat{\mathbf{K}}$. Ax=reg$coefficients[2:111] Bx=reg$coefficients[112:222] Kt=reg$coefficients[223:length(reg$coefficients)] On va passer sous silence la non identifiabilité du modèle. On peut représenter nos trois jeux de coefficients. Tout d’abord les $\hat{\mathbf{A}}$ qui nous donnent l’évolution moyenne, plot(ages[-1],Ax) On peut aussi tracer $\hat{\mathbf{K}}$ qui donne l’évolution temporelle. plot(annees,-Kt) Encore une fois, le modèle n’est pas identifiable. Si $(\hat{\mathbf{A}},\hat{\mathbf{B}},\hat{\mathbf{K}})$ est solution, il en est de même pour $(\hat{\mathbf{A}},-\hat{\mathbf{B}},-\hat{\mathbf{K}})$. Bref, le sens d’évolution n’a pas de sens. On peut représenter $-\hat{\mathbf{K}}$ qui présente l’avantage d’être décroissant (et de décrire l’amélioration des conditions de vie). Enfin, traçons aussi $-\hat{\mathbf{B}}$ plot(ages,-Bx) La difficulté est que pour prévoir l’espérance de vie, on a besoin de calculer $q_{t,x}$ pour des valeurs importantes de $t$ (non encore observée). Par exemple, on pourrait avoir envie de connaître $q_{50,2020}$ (pour quelqu’un né en 1970). On a envie d’utiliser $q_{50,2020}=\exp(\hat{A}_{50}+\hat{B}_{50}\hat{K}_{2020})$. Le soucis c’est que $\hat{K}_{2020}$ ne fait pas partie des grandeurs estimées $\hat{\mathbf{K}}$. L’idée est de construire un prédicteur $\tilde{K}_{2020}$ à partir de $\{\hat{K}_{1816},\cdots,\hat{K}_{2015}\}$. Classiquement (et c’est l’idée originale de Lee & Carter (1992)) on peut tenter un modèle $ARIMA$. Ici, la série est un peu trop régulière. On peut tenter soit un modèle exponentiel, soit un modèle linéaire (sur la série originale $\hat{\mathbf{K}}$, au delà de 1950) idx=which(annees&gt;=1950) regk=lm(log(Kt[idx])~annees[idx]) futur=2016:2125 prevKt=exp(regk$coefficients+regk$coefficients*futur) plot(annees,Kt,xlim=c(1815,2125),ylim=range(c(Kt,prevKt))) points(futur,prevKt,col="red") regk2=lm(Kt[idx]~annees[idx]) prevKt2=(regk2$coefficients+regk2$coefficients*futur) points(futur,prevKt2,col="blue") On peut alors construire la série des prédictions sur les données passées, $\hat{q}_{x,t}=\exp[\hat{A}_{x}+\hat{B}_{x}\hat{K}_{t}]$, et sur les données futures $\tilde{q}_{x,t}=\exp[\hat{A}_{x}+\hat{B}_{x}\tilde{K}_{t}]$. Si on reste sur les données passées, c’est assez simple, avec ici les probabilités de décès pour 1880 BASE=data.frame( x=base$age, t=base$annee, indice_x=1+base$age, indice_t=base$annee-1815) BASE$Ax=Ax[BASE$indice_x] BASE$Bx=Bx[BASE$indice_x] BASE$Kt=Kt[BASE$indice_t] cste=reg$coefficients BASE$Qpred=exp(cste+BASE$Ax+BASE$Bx*BASE$Kt) plot(BASE$x[BASE$t==1880],BASE$Qpred[BASE$t==1880], log="y") On peut faire la même chose pour le futur (ici pour l’année 2050) avec les deux modèles

BASE2=data.frame( x=rep(ages,length(futur)), t=rep(futur,each=length(ages)), indice_x=1+rep(ages,length(futur)), indice_t=rep(futur,each=length(ages))-2015) BASE2$Ax=Ax[BASE2$indice_x] BASE2$Bx=Bx[BASE2$indice_x] BASE2$Kt1=prevKt[BASE2$indice_t] BASE2$Kt2=prevKt2[BASE2$indice_t] cste=reg$coefficients BASE2$Qpred1=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt1) BASE2$Qpred2=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt2) head(BASE2) plot(BASE2$x[BASE2$t==2050],BASE2$Qpred1[BASE2$t== 2050],log="y") pour la prévision exponentielle, et pour la prévision linéaire On va ensuite bricoler pour créer le vecteur $\mathbf{q}$ qui contiendra soit $\hat{q}_{h,t+h}$, soit $\tilde{q}_{h,t+h}$. En effet, pour faire un calcul d’espérance de vie, on ne regarde pas à un instant donnée, mais on suit une “cohorte” dans le temps. Par exemple, pour une personne née en 1968, on a la probabilité de décès (dans l’année) suivante annee = 1968 vt=annee+ages vx=ages vq=rep(NA,length(vx)) sbase=data.frame(x=vx,t=vt) for(i in 1:nrow(sbase)){ if(sbase$t[i] &lt;= 2015) {vq[i]=BASE[(BASE$x==sbase$x[i])&amp;(BASE$t==sbase$t[i]),"Qpred"]} if(sbase$t[i] &gt; 2015) {vq[i]=BASE2[(BASE2$x==sbase$x[i])&amp;(BASE2$t==sbase$t[i]),"Qpred2"]} } plot(ages,vq,log="y") abline(v=2015-1968+.5,lty=2) avec à gauche les valeurs estimées par notre modèle (on aurait pu utiliser les valeurs brutes), et à droite, les valeurs prédites (puisqu’on parle de probabilités de décès en 2030, entre autres). Pour calculer l’espérance de vie à la naissance, on utilise le code suivant sum(cumprod(exp(-vq[1:110])))  77.62047 On peut alors faire une petite fonction, afin de visualiser l’évolution de cette espérance de vie Ex=function(annee = 1968,past="obs",fut="2"){ vt=annee+ages vx=ages sbase=data.frame(x=vx,t=vt,vq=NA) for(i in 1:nrow(sbase)){ if(sbase$t[i] &lt;= 2015) {sbase$vq[i]=BASE[(BASE$x==sbase$x[i])&amp;(BASE$t==sbase$t[i]),paste0("Q",past)]} if(sbase$t[i] &gt; 2015) {sbase$vq[i]=BASE2[(BASE2$x==sbase$x[i])&amp;(BASE2$t==sbase$t[i]),paste0("Qpred",fut)]} } vP = cumprod(exp(-(sbase$vq[1:110]))) sum(vP)}

Aussi

ANNEES=1930:2010 E2=Vectorize(Ex)(ANNEES) plot(ANNEES,E2) Si on regarde en variation, on voit qu’il y a (en gros) une variation de 0.25 d’année en année (d’où le fameux “on gagne un trimestre d’espérance de vie tous les ans”)

plot(ANNEES[-1],diff(E2),ylim=c(0,1)) Comme on plafonne ici les âges (à 110 ans) on comprend que ce gain ne peut que dimunuer, ce qu’on observe ici.

Par contre, si on prend une prévision qui conserve la variation exponentielle de $K_t$, on obtient

E1=Vectorize(function(a) Ex(a,"obs","1"))(ANNEES) plot(ANNEES,E1) ce qui n’est pas réaliste. Encore moins en considérant la variation de cette courbe (beaucoup trop optimiste)

plot(ANNEES[-1],diff(E1),ylim=c(0,1)) # Prédire ou Prévoir, réflexions (probablement) personnelles

Il y a quelques années, j’avais fait un billet où je m’énervais contre une phrase que j’avais entendu (“l’homme qui avait prédit la crise). Le point de départ était surtout un énervement lié au fait que beaucoup de journalistes vivent avec des oeillères (par exemple en associant “économiste” à “type à la télé qui parle d’économie“), et ne connaissent pas la diversité de la communauté des économistes (surtout que la pluralité n’est pas vraiment une vertue présente dans les grands médias, la manière dont sont présentés les différents candidats à la future élection présidentielle peut en être une illustration). Mais je m’écarte…. en fait, un de mes points portait sur l’utilisation du terme “prédit“.

Avec un peu de recul, je ne sais pas si la langue française à changé depuis, mais j’aurais envie de revenir sur ma distinction entre “prévoir” et “prédire“. Prévoir et prédire ont deux racines latines distinctes, prévoir venant du latin praevidere, et prédire de praedicere. Beaucoup notent que “prédire” a longtemps eu le sens de “prophétiser“, où l’annonce (de ce qui doit arriver) est souvent d’inspiration divine. En anglais, c’est presque plus compliqué (car il y a plus de synonymes) mais la même confusion existe. On peut retrouver “foretell” construit sur “fore” et “tell” (qui est à relier à une notion de prémonition), “forespell” construit sur “fore” et “spell” (on retroue la notion de sort, proche du néerlandais “voorspellen”), “foresee” construit sur “fore” et “see”, ou encore “forecast“, dont on reparlera très bientôt. A côté de ces mots, construits sur le préfixe “fore”, il existe le mot “predict“, tiré du latin praedicere, construit sur “prae” ‎(avant, que l’on retrouvé dans “pré” en français) et “dicere” ‎(dire). Comme le note le dictionnaire, “predicted” a un sens proche de “foretold” mais “without the mystical flavour”. Bref, “prédit” a longtemps eu une connotation mystique, divine.

Pourtant, les deux semblent parfois être interchangeables. Au point que fin 2006, Le Devoir écrivait « parce que cette administration n’a pas su prévoir ce qui était prévisible, nous en sommes rendus à siéger à 6 h un lundi matin! ». “prévoir ce qui était prévisible” vaste programme, non ?

Si on reste sur l’anglais, on peut tenter une distinction. Par exemple “forecast“, c’est à dire “prévoir“, est associé aux ouvrages suivants, sur amazon.com Autrement dit, “forecast” (ou “prévision“) correspond à un aspect temporel. On cherche à prévoir, à un instant donné $t$ la valeur d’une variable $x$ à une date ultérieur, disons $t+h$. Classiquement, on cherche (et ça sera notre “prévision”)

$$\mathbb{E}_t[X_{t+h}]=\mathbb{E}[X_{t+h}|x_t,x_{t-1},\cdots]$$

C’est pour cela que bon nombre de cours de séries temporelles, ou séries chronologiques sont des cours de “modèles de prévision“. Mais on commence à voir arriver des ouvrages de “méthodes prédictives“, “predictive modeling“, proposant des “prédictions Si on regarde de plus près, ces techniques sont utilisées non pas dans un contexte de séries temporelles, mais de données individuelles. Ici, on cherche à prévoir une variable d’intérêt $y$ à l’aide de covariables (variables dites explicatives $x_1,\cdots,x_k$. La prédiction sera classiquement

$$\mathbb{E}_{\boldsymbol{x}}[Y]=\mathbb{E}[Y|x_1,\cdots,x_k]$$ # Forecasting Natural Catastrophes (is rather difficult)

Following my previous post, I wanted to spend more time, on the time series with “global weather-related disaster losses as a proportion of global GDP” over the time period 1990-2016 that Roger Pilke sent me last night.

db=data.frame(year=1990:2016, ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))

In my previous post, I spend some time explaining that we should provide some sort of ‘confidence interval’ when we try to predict a pattern. That was what we call ‘model uncertainty’. But there are two (important) issues that I did not mention. (1) it is a time series, so why not use techniques dedicated to time series objects ? (2) we do not really care actually about ‘model uncertainty’ (unless we want to assess if a decreasing trend is significant, or not), and we care more about real prediction uncertainty: in the next ten years, what could be the range for the this ratio, with some given probability (say 95%)? Could we say that with 95% chance the global weather-related disaster losses as a proportion of global GDP should be (each year) within 0 and 0.35 or 0 and 0.7?

A first idea might be to use exponential smoothing techniques (without a seasonal component here).

ratio=ts(db$ratio,start=1990,frequency=1) plot(ratio,xlim=c(1990,2030)) hw=HoltWinters(ratio,gamma=FALSE) phw=predict(hw,n.ahead=15,prediction.interval = TRUE) plot(hw,phw,xlim=c(1990,2030)) polygon(c(2017:2031,rev(2017:2031)), c(phw[,2],rev(phw[,3])),border=NA,col=rgb(0,0,1,.2)) The decreasing trend is coming from the fact that exponential smoothing is here a linear regression, with weight exponentially decaying with time (the older, the smaller the weight). But we cannot use that prediction, since the ratio cannot (obviously) be negative. So why not consider, here, the logarithm of the ratio plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") hw=HoltWinters(log(ratio),gamma=FALSE) phw=predict(hw,n.ahead=15,prediction.interval = TRUE) abline(v=2016,lty=2,col="grey") lines(2017:2031,exp(phw[,2]),col="blue") lines(2017:2031,exp(phw[,3]),col="blue") lines(c(1992:2016,2017:2031),c(exp(hw$fitted[,1]),exp(phw[,1])),col="red") polygon(c(2017:2031,rev(2017:2031)),exp(c(phw[,2],rev(phw[,3]))),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) The confidence band is huge, here. What if we consider some ARIMA model here?

fit=auto.arima(ratio) farma=forecast(fit,15) farma=cbind(as.numeric(farma$fitted)[1:15],as.numeric(farma$lower[,1]),as.numeric(farma$upper[,1]),as.numeric(farma$lower[,2]),as.numeric(farma$upper[,2])) plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") lines(2017:2031,farma[,4],col="blue") lines(2017:2031,farma[,5],col="blue") lines(2017:2031,farma[,1],col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,4],rev(farma[,5])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Here, there is an intercept, but no dynamics for the time series (which is considered, here, as a pure white noise). We get exactly the same if we consider the average value of the series fit=lm(ratio~1,data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Here, we get back to my previous post, if we want to consider a possible trend (and not only an intercept)

fit=lm(ratio~year,data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Again, the confidence region is not based on inference related error, but on model uncertainty: we try to visualize where future observations might be with (say) 95% chance. Note we can also consider (why not?) a quadratic regression fit=lm(ratio~poly(year,2),data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) I am usually not a huge fan of those polynomial regression, but recently, I’ve seen that a lot in economic papers (like “if it’s not linear, add a squared version of the explanatory variable”, which is a rather odd strategy, I’ll publish some posts on that issue this year).

Here again, it might be more clever to consider a logarithmic transformation of the ratio, to insure that the ratio remains positive

fit=lm(log(ratio)~year,data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(exp(pf+s^2/2),exp(pf-1.96*s),exp(pf+1.96*s)) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(exp(predict(fit)+s^2/2),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Observe that future trend is mainly driven by the three latest observations, that were rather low (compared with older observations). What if we remove them? dbna=db db$ratio[25:27]=NA fit=lm(ratio~1,data=dbna) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016-3,lty=2,col="grey") ndb=data.frame(year=2014:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2014:2031,farma[,2],col="blue") lines(2014:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red") polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) More funny, if we consider a quadratic regression, we obtain an increasing trend for the future fit=lm(ratio~poly(year,2),data=dbna) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016-3,lty=2,col="grey") ndb=data.frame(year=2014:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2014:2031,farma[,2],col="blue") lines(2014:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red") polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) As we can see, it is rather difficult to get relevant prediction for the future, based on 25 observations…. If anyone has a suggestion, comments are open…

# Seasonal, or periodic, time series

Monday, in our MAT8181 class, we’ve discussed seasonal unit roots from a practical perspective (the theory will be briefly mentioned in a few weeks, once we’ve seen multivariate models). Consider some time series , for instance traffic on French roads,

> autoroute=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
> X=autoroute$a100 > T=1:length(X) > plot(T,X,type="l",xlim=c(0,120)) > reg=lm(X~T) > abline(reg,col="red") As discussed in a previous post, if there is a trend, we should remove it, and work on the residual > Y=residuals(reg) > acf(Y,lag=36,lwd=3) We can observe that there is some seasonality, here. A first strategy might be to assume that there is a seasonal unit root, so we consider , and we try to find some ARMA process. Consider the empirical autocorrelation function of that time series, > Z=diff(Y,12) > acf(Z,lag=36,lwd=3) or the partial autocorrelation function > pacf(Z,lag=36,lwd=3) The first graph might suggest a MA(1) structure, while the second graph might suggest an AR(1) time series. Let us try both. > model1=arima(Z,order=c(0,0,1)) > model1 Call: arima(x = Z, order = c(0, 0, 1)) Coefficients: ma1 intercept -0.2367 -583.7761 s.e. 0.0916 254.8805 sigma^2 estimated as 8071255: log likelihood = -684.1, aic = 1374.2 > E1=residuals(model1) > acf(E1,lag=36,lwd=3) which can be considered as a white noise (if you’re not convinced, try either Box-Pierce or Ljung-Box test). Similarly, > model2=arima(Z,order=c(1,0,0)) > model2 Call: arima(x = Z, order = c(1, 0, 0)) Coefficients: ar1 intercept -0.3214 -583.0943 s.e. 0.1112 248.8735 sigma^2 estimated as 7842043: log likelihood = -683.07, aic = 1372.15 > E2=residuals(model2) > acf(E2,lag=36,lwd=3) which can be also considered as a white noise. So what we have, so far is for some white noise . This suggest the following SARIMA structure on , > model2b=arima(Y,order=c(1,0,0), + seasonal = list(order = c(0, 1, 0), + period=12)) > model2b Call: arima(x = Y, order = c(1, 0, 0), seasonal = list(order = c(0, 1, 0), period = 12)) Coefficients: ar1 -0.2715 s.e. 0.1130 sigma^2 estimated as 8412999: log likelihood = -685.62, aic = 1375.25 So far, so good. Now, what if we consider that we do not have a seasonal unit root, but simply a large autoregressive coefficient in some AR structure. Let us try something like where a natural guess is that this coefficient should – probably – be close to one. Let us try this one, > model3c=arima(Y,order=c(1,0,0), + seasonal = list(order = c(1, 0, 0), + period = 12)) > model3c Call: arima(x = Y, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12)) Coefficients: ar1 sar1 intercept -0.1629 0.9741 -684.9455 s.e. 0.1170 0.0115 3064.4040 sigma^2 estimated as 8406080: log likelihood = -816.11, aic = 1640.21 which is comparable with what we got previously (somehow), so we might assume that this model can be considered as an interesting one. We will discuss further the fact that the first coefficient might be considered as non-significant. What is the difference from those two models? With a short term horizon, the two models are comparable. Clearly > library(forecast) > previ=function(model,h=36,b=40000){ + prev=forecast(model,h) + T=1:85 + Tfutur=86:(85+h) + plot(T,Y,type="l",xlim=c(0,85+h),ylim=c(-b,b)) + polygon(c(Tfutur,rev(Tfutur)),c(prev$lower[,2],rev(prev$upper[,2])),col="orange",border=NA) + polygon(c(Tfutur,rev(Tfutur)),c(prev$lower[,1],rev(prev$upper[,1])),col="yellow",border=NA) + lines(prev$mean,col="blue")
+ lines(Tfutur,prev$lower[,2],col="red") + lines(Tfutur,prev$upper[,2],col="red")
+ }

Now, on a (very) long term perspective, the models are quite different: one is stationnary, so the forecast will tend to the average value (here 0, since the trend was removed), while the other one is (seasonaly) integrated, so the confidence interval will increase. For the non stationry, we get

> previ(model2b,600,b=60000) and for the stationary one

> previ(model3c,600,b=60000) But as mentioned in the introduction of this course, forecasts with those models are relevent only for short-term horizon (say not too large). And in that case, the prediction is almost the same here,

> previ(model2b,36,b=60000) > previ(model3c,36,b=60000) Now, if we come back on our second model, we did mention previously that the autoregressive coefficient might be considered as non-significant. What if we remove it?

> model3d=arima(Y,order=c(0,0,0),
+               seasonal = list(order = c(1, 0, 0),
+               period = 12))
> (model3d)

Call:
arima(x = Y, order = c(0, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12))

Coefficients:
sar1  intercept
0.9662  -696.5661
s.e.  0.0134  3182.3017

sigma^2 estimated as 8918630:  log likelihood = -817.03,  aic = 1640.07

If we look at a (short-term) forecast, we get

> previ(model3d,36,b=32000) Do you see any difference? To be honest, I don’t… If we look at the figures, we get

> cbind(forecast(model2b,12)$mean,forecast(model3c,12)$mean,forecast(model3d,12)$mean) Time Series: Start = 86 End = 97 Frequency = 1 1 -4908.4920 -5092.8999 -5520.8780 2 -10012.7837 -9640.8103 -9493.0339 3 -3880.2202 -3841.1960 -3828.2611 4 -18102.5211 -17638.4086 -17499.1828 5 -20602.7346 -20090.9117 -19934.1066 6 -10463.2212 -10209.0139 -10132.0439 7 2458.1538 2376.4897 2351.2377 8 -1680.3342 -1654.4844 -1647.0057 9 876.6837 836.2342 823.4934 10 18046.5642 17561.6520 17413.1463 11 21531.4820 20956.3451 20780.2836 12 -3217.6103 -3152.0446 -3132.4112 Figures are different, but not significantly (keep in mind the size of the confidence interval). This might explain why, in R, when we ask for an autoregressive process or order , then we get a model with parameters to estimate, and even if some are not significant, we usually keep them for the forecast. Most of the time, from forecasting point of view, it’s no big deal. # Prévision des vente de voitures au Québec Après discussion, la date limite pour me renvoyer le second devoir est fixée à vendredi 5 octobre midi, par courriel. Je veux un fichier pdf par équipe, avec une page de garde indiquant le nom des membres du groupe, et le mot clé retenu. Maintenant, je vais en profiter pour mettre en ligne le code tapé de matin, en cours. On avait travaillé sur la modélisation de la série des ventesde voiture, au Québec. La série est en ligne, X=read.table( "http://freakonometrics.blog.free.fr/public/ data/car-sales-quebec.csv", header=TRUE,sep=";",nrows=108) Xt=ts(X[,2],start=c(1960,1),frequency=12) On peut regarder l’évolution de cette série temporelles (car les dessins, et la visualisation, c’est important), plot(Xt) On note une tendance linéaire (ou qu’on pourrait supposer linéaire), que l’on va estimer, X=as.numeric(Xt) temps=1:length(X) base=data.frame(temps,X) reg=lm(X~temps) et que l’on pourra représenter plot(temps,X,type="l",xlim=c(0,120)) abline(reg,col="red") Si on veut une série stationnaire (et c’est effectivement ce que l’on cherche), on va retrancher cette tendance à notre série, et travailler sur la série résiduelle. La série résiduelle est ici Y=X-predict(reg) plot(temps,Y,type="l") que l’on devrait pouvoir supposer stationnaire. Classiquement, on regarde les autocorrélations, acf(Y,36,lwd=4) où on repère un beau cycle annuel, mais on va surtout utiliser les autocorrélations partielles pour identifier l’ordre de la composante autorégressive. pacf(Y,36,lwd=4) La 12ème est non-nulle, et tous les autres ensuite sont significativement nulles (ou presque). On va donc tenter un . > fit.ar12=arima(Y,order=c(12,0,0)) > fit.ar12 Series: Y ARIMA(12,0,0) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.1975 0.0832 -0.1062 -0.1212 0.1437 -0.1051 0.0319 s.e. 0.0838 0.0809 0.0826 0.0843 0.0850 0.0833 0.0854 ar9 ar10 ar11 ar12 intercept -0.0332 -0.0616 0.2635 0.4913 -148.3180 s.e. 0.0853 0.0840 0.0840 0.0841 384.5095 sigma^2 estimated as 2177974: log likelihood=-946.75 AIC=1921.51 AICc=1926.03 BIC=1959.06 La douzième composante semble significative, si on repense au test de Student, u=seq(-6,6,by=.1) plot(u,dnorm(u),type="l") points(0.4913/0.0841,0,pch=19,col="red") En revanche, la huitième ne l’est pas points(-0.1018/0.0847,0,pch=19,col="blue") (voilà pour le petit retour sur la lecture des tests). Vérifions maintenant que le bruit résiduel est bien blanc… Si on regarde les autocorrélations acf(residuals(fit.ar12),36,lwd=4) Mais ce n’est pas un test de bruit blanc. En particulier, si une autocorrélation semblerait significative, on pourrait accepter que – globalement – les autocorrélations soient très proches de 0. On va alors faire des tests de Box-Pierce > Box.test(residuals(fit.ar12),lag=12, + type='Box-Pierce') Box-Pierce test data: residuals(fit.ar12) X-squared = 7.7883, df = 12, p-value = 0.8014 Pour l’interpétation du test, c’est toujours pareil: on a une statistique de l’ordre de 7, et la loi sous-jacente (la somme des carrés des 12 premières autocorrélations du bruit) est une loi du chi-deux, à 12 degrés de liberté, u=seq(0,30,by=.1) plot(u,dchisq(u,df=12),type="l") points(7.7883,0,pch=19, col="red") On accepte l’hypothèse que les 12 premières autocorrélations soient nulles. Si on va un peu plus loin, > Box.test(residuals(fit.ar12),lag=18, + type='Box-Pierce') Box-Pierce test data: residuals(fit.ar12) X-squared = 20.3861, df = 18, p-value = 0.3115 on a la même conclusion, Pour rappels, on peut retrouver à la main la p-value, > 1-pchisq(20.3861,df=18)  0.3115071 Faisons notre petite dessin de toutes les p-value, pour s’assurer que le bruit est bien un bruit blanc,  BP=function(h) Box.test(residuals(fit.ar12),lag=h, type='Box-Pierce')$p.value
plot(1:24,Vectorize(BP)(1:24),type='b',
ylim=c(0,1),col="red")
abline(h=.05,lty=2,col="blue") Cette fois c’est bon, on tient notre premier modèle, un . Si on veut aller plus loin, on peut regarder ce que donnerait de la sélection automatique,

> library(caschrono)
> armaselect(Y,nbmod=5)
p q      sbc
[1,] 14 1 1635.214
[2,] 12 1 1635.645
[3,] 15 1 1638.178
[4,] 12 3 1638.297
[5,] 12 4 1639.232

On peut être tenté d’aller voir ce qui se passerait en rajoutant une composante moyenne mobile à notre modèle précédant,

> fit.arma12.1=arima(Y,order=c(12,0,1))
> fit.arma12.1
Series: Y
ARIMA(12,0,1) with non-zero mean

Coefficients:
ar1     ar2      ar3      ar4     ar5      ar6     ar7
0.0301  0.1558  -0.0941  -0.1461  0.1063  -0.0688  -0.002
s.e.  0.1235  0.0854   0.0757   0.0784  0.0807   0.0774   0.080
ar9     ar10    ar11    ar12     ma1  intercept
-0.0646  -0.0798  0.2538  0.5786  0.2231  -131.3495
s.e.   0.0802   0.0766  0.0751  0.0861  0.1393   368.8156

sigma^2 estimated as 2127759:  log likelihood=-945.65
AIC=1921.31   AICc=1926.52   BIC=1961.54

Le nouveau coefficient,  est à peine significatif. Éventuellement avec un seuil à 10%… Pourquoi pas? Si on regarde les résidus, sans grande surprise, on a toujours un bruit blanc, encore plus blanc qu’auparavant (en violet sur le dessin ci-dessous)

BP=function(h) Box.test(
residuals(fit.arma12.1),lag=h,
type='Box-Pierce')$p.value plot(1:24,Vectorize(BP)(1:24),type='b', ylim=c(0,1),col="red") abline(h=.05,lty=2,col="blue") BP=function(h) Box.test(residuals(fit.ar12),lag=h, type='Ljung-Box')$p.value
lines(1:24,Vectorize(BP)(1:24),col="purple",
type="b") On peut aussi aller voir parmi les modèles proposés, en particulier le modèle .

> fit.ar14=
+ arima(Y,order=c(14,0,0),method="CSS")
> fit.ar14
Series: Y
ARIMA(14,0,0) with non-zero mean

Coefficients:
ar1     ar2      ar3      ar4     ar5      ar6     ar7
0.2495  0.2105  -0.0584  -0.1569  0.1282  -0.1152  0.0268
s.e.  0.0956  0.0972   0.0854   0.0830  0.0838   0.0840  0.0847
ar9     ar10    ar11    ar12     ar13     ar14  intercept
-0.0327  -0.1116  0.2649  0.5887  -0.1575  -0.1572    80.5
s.e.   0.0855   0.0851  0.0853  0.0886   0.1031   0.0999   338.9

sigma^2 estimated as 2218612:  part log likelihood=-942.31

C’est un peu tiré par les cheveux, mais on pourrait accepter l’hypothèse que soit significativement non-nuls. Mais on est encore limite…  Allez, on l’accepte aussi dans notre gang de modèle.

On a finalement trois modèles. Si on fait un peu de backtesting, sur les 12 derniers mois,

T=length(Y)
backtest=12
subY=Y[-((T-backtest+1):T)]
subtemps=1:(T-backtest)
plot(temps,Y,type="l")
lines(subtemps,subY,lwd=4)
fit.ar12.s=arima(subY,
order=c(p=12,d=0,q=0),method="CSS")
fit.arma12.1.s=arima(subY,
order=c(p=12,d=0,q=1),method="CSS")
fit.ar14.s=arima(subY,
order=c(p=14,d=0,q=0),method="CSS")
p.ar12=predict(fit.ar12.s,12)
pred.ar12=as.numeric(p.ar12$pred) p.arma12.1=predict(fit.arma12.1.s,12) pred.arma12.1=as.numeric(p.arma12.1$pred)
p.ar14=predict(fit.ar14.s,12)
pred.ar14=as.numeric(p.ar14$pred) on obtient les prévisions suivantes (avec les valeurs observées dans la première colonne) > (M=cbind(observé=Y[((T-backtest+1):T)], + modèle1=pred.ar12, + modèle2=pred.arma12.1, + modèle3=pred.ar14)) observé modèle1 modèle2 modèle3 97 -4836.2174 -5689.3331 -5885.4486 -6364.2471 98 -3876.4199 -4274.0391 -4287.2193 -4773.8116 99 1930.3776 1817.8411 2127.9915 2290.1460 100 3435.1751 4089.3598 3736.1110 4039.4150 101 7727.9726 6998.9829 7391.6694 7281.4797 102 2631.7701 3456.8819 3397.5478 4230.5324 103 -509.4324 -2128.6315 -2268.9672 -2258.7216 104 -1892.6349 -3877.7609 -3694.9409 -3620.4798 105 -4310.8374 -3384.0905 -3430.4090 -2881.4942 106 2564.9600 -504.6883 -242.5018 183.2891 107 -1678.2425 -1540.9904 -1607.5996 -855.7677 108 -4362.4450 -3927.4772 -3928.0626 -3718.3922 > sum((M[,1]-M[,2])^2)  19590931 > sum((M[,1]-M[,3])^2)  17293716 > sum((M[,1]-M[,4])^2)  21242230 I.e. on aurait envie de retenir le second modèle, . On va maintenant l’utiliser pour faire un peu de prévision, library(forecast) fit.arma12.1= arima(Y,order=c(12,0,1)) fit.arma12.1 PREDARMA=forecast(fit.arma12.1,12) plot(Y,xlim=c(1,120),type="l") temps=T+1:12 polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2],
rev(PREDARMA$upper[,2])),col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1],
rev(PREDARMA$upper[,1])),col="orange",border=NA) lines(temps,PREDARMA$mean,col="red") Cela dit, on peut aussi aller beaucoup plus loin dans la prévision,

PREDARMA=forecast(fit.arma12.1,120)
plot(Y,xlim=c(1,210),type="l")
temps=T+1:120
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2], rev(PREDARMA$upper[,2])),col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1], rev(PREDARMA$upper[,1])),col="orange",border=NA)
lines(temps,PREDARMA$mean,col="red") Bon, on y est presque, car on a modélisé , la série obtenue en enlevant la tendance linéaire. Pour remonter sur , on va rajouter la tendance à la prévision faite auparavant, X=as.numeric(Xt) temps=1:length(X) plot(temps,X,type="l",xlim=c(0,210), ylim=c(5000,30000)) base=data.frame(temps,X) reg=lm(X~temps) abline(reg,col="red") PREDTENDANCE=predict(reg,newdata= data.frame(temps=T+1:120)) temps=T+1:120 polygon(c(temps,rev(temps)),c(PREDTENDANCE+ PREDARMA$lower[,2],rev(PREDTENDANCE+PREDARMA$upper[,2])), col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDTENDANCE+ PREDARMA$lower[,1],rev(PREDTENDANCE+PREDARMA$upper[,1])), col="orange",border=NA) lines(temps,PREDTENDANCE+PREDARMA$mean,col="red") # Unit root, or not ? is it a big deal ?

Consider a time series, generated using

set.seed(1)
E=rnorm(240)
X=rep(NA,240)
rho=0.8
X=0
for(t in 2:240){X[t]=rho*X[t-1]+E[t]}

The idea is to assume that an autoregressive model can be considered, but we don’t know the value of the parameter. More precisely, we can’t choose if the parameter is either one (and the series is integrated), or some value strictly smaller than 1 (and the series is stationary). Based on past observations, the higher the autocorrelation, the lower the variance of the noise.

rhoest=0.9; H=260
u=241:(240+H)
P=X*rhoest^(1:H)
s=sqrt(1/(sum((rhoest^(2*(1:300))))))*sd(X)

Now that we have a model, consider the following forecast, including a confidence interval,


plot(1:240,X,xlab="",xlim=c(0,240+H),
ylim=c(-9.25,9),ylab="",type="l")
V=cumsum(rhoest^(2*(1:H)))*s
polygon(c(u,rev(u)),c(P+1.96*sqrt(V),
rev(P-1.96*sqrt(V))),col="yellow",border=NA)
polygon(c(u,rev(u)),c(P+1.64*sqrt(V),
rev(P-1.64*sqrt(V))),col="orange",border=NA)
lines(u,P,col="red")
Here, forecasts can be derived, with any kind of possible autoregressive coefficient, from 0.7 to 1. I.e. we can chose to model the time series either with a stationary, or an integrated series, As we can see above, assuming either that the series is stationary (parameter lower – strictly – than 1) or integrated (parameter equal to 1), the shape of the prediction can be quite different. So yes, assuming an integrated model is a big deal, since it has a strong impact on predictions.

# Sondages et prévisions de séries temporelles

Ce soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant

sondage=read.table("http://freakonometrics.blog.free.fr/
public/data/sondage1ertour2012.csv",
sondage[36,1]="11/03/12"
sondeur=unique(sondage$SONDEUR) sondage$date=as.Date(as.character(sondage$DATE), "%d/%m/%y") (je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance, CL=brewer.pal(6, "RdBu") couleur=c(CL,CL,"orange","grey") datefin=as.Date("2012/04/12") k=1 plot(sondage$date,sondage[,k+2],col=couleur[k],
pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40),
xlim=c(min(sondage$date),datefin)) rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
for(k in 2:4){
points(sondage$date,sondage[,k+2], col=couleur[k],pch=19,cex=.7) rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
} Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste

rl=lowess(sondage$date,sondage[,k+2]) Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser), D=data.frame(date=min(sondage$date)+0:148)
library(splines)
k=1
plot(sondage$date,sondage[,k+2], col=couleur[k],pch=19,cex=.7, xlab="",ylab="",ylim=c(0,40),xlim= c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL,border=NA)
lines(D$date,prl,lwd=2,col=CL) points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,)
for(k in 2:4){
points(sondage$date,sondage[,k+2], col=couleur[k],pch=19,cex=.7) } avec, là encore, la commande suivante pour lisser rs=lm(sondage[,k+2]~bs(date),data=sondage) Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets), On a ici la courbe suivante avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy, k=2 plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin)) rs=lm(sondage[,k+2]~bs(date),data=sondage) prl=predict(rs,newdata=D) prlse=sqrt(prl/100*(1-prl/100)/1000)*100 polygon(c(D$date,rev(D$date)),c(prl+2*prlse, rev(prl-2*prlse)),col=CL,border=NA) lines(D$date,prl,lwd=2,col=CL)
points(sondage$date,sondage[,k+2],col= couleur[k],pch=19,cex=.7,) for(k in c(1,3:4)){ points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7)
} Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête, Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“.

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait…

Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via

D=data.frame(date=min(sondage$date)+0:148) prl=predict(rs,newdata=D) On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections) On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante library(mnormt) simulation=function(S){ proba=c(S,100-sum(S))/100 variance= -proba%*%t(proba) diag(variance)=proba*(1-proba) variance=variance/1000 return(rmnorm(1,proba[1:4],variance[1:4,1:4]))} simulsondages=function(M){ Z=rep(NA,ncol(M)) for(i in 1:nrow(M)){ Z=rbind(Z,simulation(M[i,])) } return(Z[-1,])} prediction4=matrix(NA,nrow(D),4) for(k in 1:4){ rs=lm(sondage[,k+2]~bs(date),data=sondage) prediction4[,k]=predict(rs,newdata=D) } Par exemple, si on la fait tourner une fois, on obtient le graphique suivant set.seed(1) S4=simulsondages(prediction4) S100=100*S4 k=1 plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin)) rs=lm(sondage[,k+2]~bs(date),data=sondage) lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k], pch=19,cex=.7) rs=lm(sondage[,k+2]~bs(date),data=sondage) lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
} (les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

library(forecast)
k=1
X=S100[,k]
ETS=ets(X)
F=forecast(ETS,h=60)
fdate=max(D$date)+1:60 k=1 plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin)) rs=lm(sondage[,k+2]~bs(date),data=sondage) lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k], pch=19,cex=.7) rs=lm(sondage[,k+2]~bs(date),data=sondage) lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]), rev(as.numeric(F$upper[,2]))),col=CL,border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL,border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL) Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes ETS=ets(X) F=forecast(ETS,h=60) Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante On peut bien entendu faire une projection similaire pour Nicolas Sarkozy, k=2 X=S100[,k] ETS = ets(X) F=forecast(ETS,h=60) k=1 plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin)) rs=lm(sondage[,k+2]~bs(date),data=sondage) lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k], pch=19,cex=.7) rs=lm(sondage[,k+2]~bs(date),data=sondage) lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]), rev(as.numeric(F$upper[,2]))),col=CL,border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL,border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL) On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet. On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus, set.seed(1) for(sim in 1:20){ Ssim=simulsondages(prediction4) F1=forecast(ets(100*Ssim[,1]),h=60) F2=forecast(ets(100*Ssim[,2]),h=60) lines(fdate,as.numeric(F1$mean),col=CL)
lines(fdate,as.numeric(F2$mean),col=CL)} Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour set.seed(1) VICTOIRE=rep(NA,1000) for(sim in 1:1000){ Ssim=simulsondages(prediction4) F1=forecast(ets(100*Ssim[,1]),h=60) F2=forecast(ets(100*Ssim[,2]),h=60) VICTOIRE[sim]=(F1$mean>F2$mean)} mean(VICTOIRE) Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes ! # Combien de temps profite-t-on de ses grands parents ? Ce Hier matin, je suis tombé un peu par hasard sur deux graphiques de l’INSEE (en France) avec l’age moyen des mères à l’accouchement, en fonction de rang de naissance de l’enfant, avec tout d’abord 1905-1965, puis 1960-2000 Ces graphiques sont passionnants en soi – comme en ont témoigné pas mal de followers sur Twitter – mais ils m’ont fait m’interroger. En particulier, sur la croissance observée depuis 30 ans, qui me faisait penser à la tendance croissante observée sur les durées de vie. On n’a – malheureusement – pas accès au données complètes sur le site, mais on peut trouver d’autres donnéesintéressantes (en l’occurrence l’age moyen à la naissance). > agenaissance=read.table("http://freakonometrics.blog.free.fr/ public/data/agenaissance.csv",header=TRUE,sep=",") > agenaissance$Age=as.character(agenaissance$AGE) > agenaissance$AGE=as.numeric(substr(agenaissance$Age,1,2))+ + as.numeric(substr(agenaissance$Age,4,4))/10
> plot(agenaissance$ANNEE+.5,agenaissance$AGE,
+ type="l",lwd=2,col="blue")

Visuellement, on retrouve la courbe en bleu foncée sur les graphiques ci-dessus, On peut alors aller en cran plus loin, en se demandant non pas quel était l’âge moyen de la mère, mais de la grand-mère (au sens la mère de la mère)

> agenaissance$NAIS.MERE=(agenaissance$ANNEE+.5)-
+ agenaissance$AGE > w=(trunc(agenaissance$NAIS.MERE-.5))
> rownames(agenaissance)=agenaissance$ANNEE > a1=agenaissance[as.character(w),]$NAIS.MERE
> a2=agenaissance[as.character(w+1),]$NAIS.MERE > p=agenaissance$NAIS.MERE-(w+.5)
> agenaissance$NAIS.GRD.MERE=(1-p)*a1+p*a2 > agenaissance$age.GRD.MERE=agenaissance$ANNEE+.5- + agenaissance$NAIS.GRD.MERE
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE
2000  2000 30.3 30,3     1970.2       1942.87        57.63
2001  2001 30.4 30,4     1971.1       1943.80        57.70
2002  2002 30.4 30,4     1972.1       1944.92        57.58
2003  2003 30.5 30,5     1973.0       1945.95        57.55
2004  2004 30.5 30,5     1974.0       1947.05        57.45
2005  2005 30.6 30,6     1974.9       1948.04        57.46
> plot(agenaissance$ANNEE+.5,agenaissance$age.GRD.MERE,
+ type="l",lwd=2,col="red")

Là encore, on peut visualiser l’âge de la grand-mère maternelle à la naissance A partir de là, on peut se demander combien de temps on profite de ses grands-parents (ou tout du moins ici de sa grand mère maternelle), en se basant sur les calculs d’espérance de vie résiduelle. En utilisant le modèle de Lee-Carter pour modéliser les taux de décès annuel, et en extrapolant sur le siècle en cours, on peut extrapoler les espérances de vie résiduelles.

> Deces <- read.table("http://freakonometrics.free.fr/
> Deces$Age <- as.numeric(as.character(Deces$Age))
> Deces$Age[is.na(Deces$Age)] <- 110
> Expo$Age <- as.numeric(as.character(Expo$Age))
> Expo$Age[is.na(Expo$Age)] <- 110
>  library(forecast)
>  library(demography)
>  YEAR <- unique(Deces$Year);nC=length(YEAR) > AGE <- unique(Deces$Age);nL=length(AGE)
>  MUF  <- matrix(Deces$Female/Expo$Female,nL,nC)
>  POPF <- matrix(Expo$Female,nL,nC) > BASEF <- demogdata(data=MUF, pop=POPF,ages=AGE, + years=YEAR, type="mortality", + label="France", name="Femmes", lambda=1) > LCF <- lca(BASEF) > LCFf<-forecast(LCF,h=100) > A <- LCF$ax
> B <- LCF$bx > K1 <- LCF$kt
> K2 <- K1[length(K1)]+LCFf$kt.f$mean
> K <- c(K1,K2)
> MU <- matrix(NA,length(A),length(K))
> for(i in 1:length(A)){
+ for(j in 1:length(K)){
+ MU[i,j] <- exp(A[i]+B[i]*K[j]) }}
> esp.vie = function(xentier,T){
+ s <- seq(0,99-xentier-1)
+ MUd <- MU[xentier+1+s,T+s-1898]
+ Pxt <- cumprod(exp(-diag(MUd)))
+ ext <- sum(Pxt)
+ return(ext) }
> EVIE = function(x,T){
+ x1 <- trunc(x)
+ x2 <- x1+1
+ return((1-(x-x1))*esp.vie(x1,T)+(x-x1)*esp.vie(x2,T)) }
> agenaissance$EV=NA > for(i in 1:100){ + t <- 2006-i + agenaissance$EV[agenaissance$ANNEE==t]= + EVIE(x=agenaissance$age.GRD.MERE[
+ agenaissance$ANNEE==t],t) } > tail(agenaissance) ANNEE AGE Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE EV 2000 30.3 30,3 1970.2 1942.87 57.63 29.13876 2001 30.4 30,4 1971.1 1943.80 57.70 29.17047 2002 30.4 30,4 1972.1 1944.92 57.58 29.39027 2003 30.5 30,5 1973.0 1945.95 57.55 29.52041 2004 30.5 30,5 1974.0 1947.05 57.45 29.72511 2005 30.6 30,6 1974.9 1948.04 57.46 29.80398 Autrement dit, sur la dernière ligne, l’espérance de vie (résiduelle) pour une femme de 57.46 ans en 2005 était d’environ 29.80 ans. On peut alors visualiser non seulement l’âge moyen de sa grand-mère à la naissance, mais son espérance de vie résiduelle, > plot(agenaissance$ANNEE+.5,agenaissance$EV, + type="l",lwd=2,col="purple") On note que depuis 30 ans, en France, la durée (moyenne) pendant laquelle les petits-enfants vont profiter de leur grands parents s’est stabilisé à une trentaine d’années. On peut aussi continuer, et remonter d’un cran (en refaisant tourner le code avec quelques modifications): on a alors l’âge (moyen) de son arrière grand mère à la naissance, et la durée de vie (résiduelle) des arrière grand mères On manque ici de données, mais il semble que l’on profite – en moyenne – environ 5 ans de son arrière grand mère. Maintenant on peut aussi s’interroger sur les limites de cette étude rapide. En particulier, de même qu’il existe une corrélation forte entre les durées de vie de conjoints (e.g. broken heart syndrom de Jagger & Sutton (1991)), on peut se demander si la naissance d’enfants et de petits-enfants a un impact sur la durée de vie résiduelle d’une personne (ou si on peut supposer l’indépendance comme on l’a fait ici). # Time horizon in forecasting, and rules of thumb I recently received an email about forecasting and rules of thumb. “Dans la profession […] se transmet une règle empirique qui voudrait que l’on prenne un historique du double de l’horizon de prévision : 20 ans de données pour une prévision à 10 ans, etc… Je souhaite savoir si cette règle n’aurait pas, par hasard, un fondement théorique quitte à ce que le rapport ne soit pas de 2 pour 1, mais de 3 pour 1 ou de 1 pour 1 par exemple.” To summarize briefly, the rule is to consider a 2-1 ratio for the period of observation vs. forecast horizon. And the interesting question is if there are justifications for such a rule… At first, I remembered a rules of thumb, from the book by Box and Jenkins, which states that it is meaningless to look at autocorrelations when lags exceed the sample size over 6. So with 12 years of data, autocorrelations with a lag higher than two years are useless. But it is not what is mentioned here. So I looked at some dataset, and some standard time series models. • It depends on the series It might obvious… but if it is the case, it means that it will be difficult to have a general rule of thumb. Consider e.g. the number of airline passengers, library(forecast) X = AirPassengers ETS = ets(X) plot(forecast(ETS,h=length(X)/2)) or some sales in a big store,  or car casualties in France, or the temperature in Nottingham Castle,  or the water level at Lack Hurron, or the flow of the Nile river,  or see also here for forecasting techniques in demography. Actually, in the case of life insurance, actuaries have to forecast future demography, i.e. try to assess death rates of those who currently purchase retirement contracts, who might be 20 years old. So they have to forecast death rate until 2100, say. One the one hand, it sounds difficult to make forecast over a century (it is already difficult for climate, I guess it is even more complex for human life). On the other hand, a 2-1 ratio means that we have to use data from 1800… Here again, it is difficult to justify that mortality in the 1850 could be interesting to say anything about mortality in 2050. So I guess it will be difficult to justify the use of general rules of thumb…. • It depends on the model Consider the following (simulated) series. Several models can be fitted. And the shape on the forecast (and the forecast error) will depend on the model considered. The benchmark can be the model without any dynamics, i.e. we assume that observations are i.i.d. Or more classically, assume that it is simple a white noise, i.e. an i.i.d centered process. Then the forecast is the following, With that kind of assumption, we see that the 2-1 ratio is useless since we can get forecasts up to any horizon…. But that does not seem very robust. For instance, if we consider exponential smoothing techniques, we can obtain Which is rather different. And with the 2-1 ratio, obviously, there is a lot of uncertainty at the end ! It would be even worst if we assume that we look at a random walk. Because actually a dozen models – at least – can be considered, from ARIMA, seasonal ARIMA, Holt Winters, Exponential Smoothing, etc… So I do not see any theoretical justification of that rule of thumb. Obviously, the maximum horizon can not be extremely far away if the series is non-stationary, with a very irregular pattern, and with a lot of noise… So we’re back at the beginning. If anyone is willing to share his or her experience, comments are open. # I really need to find hot (and sexy) topics 50 days ago (here), I was supposed to be very optimistic about the probability that I could reach a million viewed pages on that blog (over a bit more than two years). Unfortunately, the wind has changed and today, the probability is quite low…  base=read.table("millionb.csv",sep=";",header=TRUE) X1=cumsum(base$nombre)
X2=cumsum(base$nombre) X=X1+X2 D=as.Date(as.character(base$date),"%m/%d/%Y")
kt=which(D==as.Date("01/06/2010","%d/%m/%Y"))
D0=as.Date("08/11/2008","%d/%m/%Y")
D=D0+1:length(X1)
P=rep(NA,(length(X)-kt)+1)
for(h in 0:(length(X)-kt)){
model  <- arima(X[1:(kt+h)],c(7 1,7),method="CSS")
forecast <- predict(model,200)
u=max(D[1:kt+h])+1:300
k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
(P[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))
}
plot( D[length(D)-length(P)]+1:220,c(P,rep(NA,220-length(P))),
ylab="Probability to reach 1,000,000",xlab="",
type="l",col="red",ylim=c(0,1)) So, I guess my posts on multiple internal rates of return, or Young’s inequality will have to wait next year… I really need to find some more sexy post to attract readers.. Challenge accepted ! 