Tag Archives: prévision

Modèles de Prévision, suite

La seconde partie de transparents pour le cours ACT6420 est en ligne. Il s’agit d’un document qui permettra à tout le monde garder des traces, mais je n’entends pas forcément les suivre… Par exemple, ce soir, je propose de jouer un peu avec les données d’incendies à Chicago

> chicago=read.table("http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)

ou une base sur les états américains (cf un ancien billet)

> US=read.table("http://freakonometrics.free.fr/US.txt", + header=TRUE,sep=";")

Sinon, la semaine passée, on avait évoqué la matrice chapeau (dont on reparlera davantage lorsque l’on évoquera la phase de diagnostique). L’idée était de noter que

https://latex.codecogs.com/gif.latex?%20%20%20%20%20\widehat{\boldsymbol{Y}}%20=%20H%20\boldsymbol{Y}

avec

https://latex.codecogs.com/gif.latex?%20%20%20%20%20H%20=%20\boldsymbol{X}%20\left(\boldsymbol{X}^\top%20\boldsymbol{X}%20\right)^{-1}%20\boldsymbol{X}^\top

On peut visualiser cette matrice avec le code suivant

> X=cbind(1,cars$speed)
> M=X %*% solve( t(X) %*% X) %*% t(X)
> Q=matrix(as.numeric(cut(M,c(-10,quantile(M,(1:5)/6),10),labels=1:6)),50,50)
> plot(1:50,1:50,col="white",axes=FALSE,xlab="",ylab="",ylim=c(50,0))
> library(RColorBrewer) 
> couleur=rev(brewer.pal(6,"RdBu")) 
> for(i in 1:50){for(j in 1:50){ points(i,j,col=couleur[Q[i,j]],pch=15,cex=1)}}
> axis(1)
> axis(2)

avec des cases rouges pour les valeurs importantes (et positive) et des case bleues pour des valeurs importantes (mais négatives)

Modèles de Prévision

Cet session d’été, je redonne le cours ACT6420, modèles de prévision. Le plan de cours est maintenant en ligne (comme souvent, les salles de cours ne sont pas encore attribuées définitivement mais elles seront disponibles sur le site de l’UQAM). Les transparents pour la première partie de cours sont en ligne (il s’agit de rappels de statistique, et d’une introduction aux modèles de régression).

Des exercices et des codes informatiques seront mis en ligne sur /courses/modeles-de-prevision.

Durant le cours, j’essayerais d’illustrer avec le maximum d’exemples sur des données réelles, à partir de codes R. Quelques séances de démonstrations se feront en salle informatique, afin de découvrir les rudiments de R. Pour aller plus loin, je recommande

  • “R pour les débutants”d’Emmanuel Paradis (PDF)
  • “Brise Glace-R” d’Andrew Robinson er Arnaud Schloesing (PDF)
  • “Introduction à la programmation en R” de Vincent Goulet (PDF)

ainsi que quelques tutoriels en ligne, par exemple, dans un ancien billet, mais aussi sur https://youtube.com/.

Elections Sans Bipartisme

Hier, sur Twitter, @JF_Godbout partageait un joli graphique relatif aux élections québécoises, avec les nombres de votes obtenus (ici en pourcentage des votes totaux) et le pourcentage de sièges que cela donne,

Il faut dire qu’hier, c’était jour d’élection. Ce sont des élections à un tour, avec plusieurs partis (disons 4 grands partis si on se limite aux élections les plus récentes). Jusqu’à lundi, le PQ (Parti Québécois) avait une majorité (relative), et depuis lundi, le PLQ (Parti Libéral) est revenu au pouvoir. Et il y a à côté de ces deux grands partis d’autres partis, que l’on va regrouper. @JF_Godbout m’ayant envoyé ses données, j’ai pu reproduire les graphiques.

La grande différence ici, entre nos approches est – je pense – la méthode de lissage: alors que  @JF_Godbout faisait des régressions locales, j’ai opté pour des splines.

reg=glm(Seat~bs(Election),data=d,family=binomial)
yp=predict(reg,newdata=data.frame(Election=1970:2015),type="response")

J’utilise des splines parce que ça se généralise facilement à d’autres lois, et qu’on peut aussi les utiliser quand on régresse sur plusieurs variables continues. Si on rajoute les résultats de l’élection de lundi, on obtient le graphique suivant

Ce qui me gêne dans cette modélisation, c’est qu’on considère nos trois dernières régressions comme indépendantes…. alors que ce n’est pas le cas: on a trois modalités possibles, et la somme doit forcément faire 1. D’ailleurs, quand on somme les prédictions, on dépasse 1.

Avec deux partis, notre régression logistique nous garantie d’avoir une somme à un, mais plus ici. A moins d’utiliser un modèle multinomial. Si on oublie le lissage, l’idée est d’écrire

qui est tout simplement l’extension du modèle logistique. La difficulté technique est que, sous R, quand on fait une régression logistique, la variable d’intérêt doit être une variable prenant trois facteurs, et non pas un vecteur de probabilités. Au lieu d’utiliser des taux, ou des pourcentages, je devrais utiliser des nombres d’observations de chaque modalités. On va donc crééer une fausse base

> vote=data.frame(
+ plq=plq$Vote,
+ pq=pq$Vote,
+ autre=autre$Vote,
+ Election=pq$Election)
> vote
      plq     pq  autre Election
1  0.5465 0.3022 0.1513     1972
2  0.3378 0.4137 0.2485     1976
3  0.4608 0.4926 0.0466     1981
4  0.5599 0.3869 0.0532     1985
5  0.4995 0.4016 0.0989     1989
6  0.4440 0.4475 0.1085     1994
7  0.4355 0.4287 0.1358     1998
8  0.4599 0.3324 0.2077     2003
9  0.3308 0.2835 0.3857     2007
10 0.4208 0.3517 0.2275     2008
11 0.3120 0.3195 0.3685     2012
12 0.4140 0.2538 0.3312     2014
> M=t(round(vote[,1:3]*1000))
> M
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
plq    546  338  461  560  500  444  436  460  331   421   312   414
pq     302  414  493  387  402  448  429  332  284   352   320   254
autre  151  249   47   53   99  108  136  208  386   227   368   331

On réplique alors autant de fois le nom du parti qu’on a de (pseudo) votes

> Y=X=NULL
> for(i in 1:3){
+ for(j in 1:ncol(M)){
+ Y=c(Y,rep(c("PLQ","PQ","Autre")[i],M[i,j]))
+ X=c(X,rep(pq$Election[j],M[i,j]))
+ }}
> basev=data.frame(Y=as.factor(Y),X=X)

C’est un peu artificiel, j’en conviens, mais ça marche… On peut alors faire la régression

> library(nnet)
> regv=multinom(Y~bs(X),data=basev)
# weights:  15 (8 variable)
initial  value 13186.643301 
iter  10 value 12241.931805
final  value 12221.591172 
converged

et des prévisions

> ypv=predict(regv,newdata=
+ data.frame(X=1973:2014),type="probs")
> ypv
        Autre       PLQ        PQ
1  0.18426281 0.4705415 0.3451957
2  0.16051060 0.4767408 0.3627486
3  0.14149153 0.4805514 0.3779571
4  0.12631826 0.4826100 0.3910718
5  0.11426236 0.4834165 0.4023212
6  0.10474046 0.4833543 0.4119053
7  0.09729304 0.4827117 0.4199953
8  0.09156246 0.4817016 0.4267360
9  0.08727323 0.4804776 0.4322491
10 0.08421535 0.4791477 0.4366370
11 0.08223095 0.4777835 0.4399855
12 0.08120364 0.4764292 0.4423672
13 0.08105030 0.4751066 0.4438431
14 0.08171457 0.4738203 0.4444651

Et on notera que cette fois, nos prévisions somment à un.

Cela dit, visuellement, on obtient des choses proches de ce qu’on avait auparavant,

mais intellectuellement, je trouve plus logique d’utiliser autre chose que des régressions logistiques quand on est dans une élection où plusieurs partis sont en compétition.

Estimation et prévision pour des séries temporelles

Pour la fin de cours de modèles de prévision , quelques transparents sur l’identification et l’estimation des modèles SARIMA, quelques compléments sur les tests (racines unités et non-stationnarité, ainsi que saisonnalité), et enfin, quelques pistes pour construire des prédictions (avec une quantification de l’incertitude), avec des codes R. Les transparents sont en ligne ici (même si la page de garde est identique aux autres, il s’agit de nouveau matériel). Maintenant que les transparents sont finis (et en ligne) les prochains billets seront orientés autour de la modélisation et des aspects computationnels.

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 https://latex.codecogs.com/gif.latex?AR(12).

> 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 https://latex.codecogs.com/gif.latex?\phi_{12} 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 https://latex.codecogs.com/gif.latex?\phi_{8} 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)
[1] 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 https://latex.codecogs.com/gif.latex?AR(12). 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, https://latex.codecogs.com/gif.latex?\theta_1 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 https://latex.codecogs.com/gif.latex?AR(14).

> 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 https://latex.codecogs.com/gif.latex?\phi_{14}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)
[1] 19590931
>  sum((M[,1]-M[,3])^2)
[1] 17293716
>  sum((M[,1]-M[,4])^2)
[1] 21242230

I.e. on aurait envie de retenir le second modèle, https://latex.codecogs.com/gif.latex?ARMA(12,1). 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é https://latex.codecogs.com/gif.latex?(Y_t), la série obtenue en enlevant la tendance linéaire. Pour remonter sur https://latex.codecogs.com/gif.latex?(X_t), 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")