La dernière fois, nous avions vu 4 modèles de séries temporelles, pour modéliser la série hebdomadaire de températures. Le premier était obtenu avec auto.arima, les deux étaient ensuite ses modèles SARIMA (le second était une forte autocorélation à l’ordre 52, et le troisième avait une racine unité saisonnière), et le dernier était obtenu par la méthode de Buys-Ballot (cyclique, avec un bruit autocorrélé).
Au lieu de se demander lequel des quatre est le “meilleur”, on peut tenter de trouver une moyenne pondérée, \hat{y}_t = \sum_{i}\omega_i \hat{y}_t^{(i)} où \omega=(\omega_i) est un point du simplexe, au sens où \omega_i\geq 0 et \omega_1+...+\omega_k=1. On veut ensuite trouver les “meilleurs” pondérations, par exemple au sens de la norme L^2. On va calibrer nos quatre modèles sur les m premières valeurs, et on compare ensuite la combinaison des prévions sur les 111 prochaines valeurs (et les vraies valeurs), \omega^\star=\text{argmin}\{\sum_{h=1}^{111}\left(y_{t+h}- \sum_{i}\omega_i \hat{y}_{t+h}^{(i)}\right)^2\}
Par exemple, on peut utiliser les 200 premières valeurs.
m=200 Ym=ts(y[1:m],start=1996,frequency=52) |
On va ensuite estimer les 4 modèles sur ces 200 valeurs
model1=arima(Ym,order = c(4,0,0), seasonal = list(order = c(1,0,0))) model2=arima(Ym,order = c(1,0,0), seasonal = list(order = c(1,1,0))) model3=arima(Ym,order = c(0,0,2), seasonal = list(order = c(0,1,1))) regL=lm(Temp~0+as.factor(NumWeek), data=electricite[1:m,]) E=residuals(regL) model4=arima(ts(E,start=1996, frequency=52),order = c(1,0,0)) |
On peut ensuite collecter les prédictions
DONNEES = data.frame( y1=predict(model1,n.ahead = 111)$pred, y2=predict(model2,n.ahead = 111)$pred, y3=predict(model3,n.ahead = 111)$pred, y4=predict(model4,n.ahead = 111)$pred+ predict(regL,newdata=electricite[m+(1:111),]), vraie=electricite[m+(1:111),"Temp"]) |
Pour créer une combinaison linéaire des prévisions, on utilise
a=rep(1/4,4) y_prev = as.matrix(DONNEES[,1:4]) %*% a |
On peut ainsi visualiser les 4 prévisions, leur combinaison linéaire (ici avec des équi-pondérations), et leurs observations
plot(1:111,DONNEES[,5],type="l",lwd=2) for(i in 1:4) lines(1:111,DONNEES[,i],col="blue") lines(1:111,yp,col="red") |
Pour trouver la valeur “optimale” des pondérations, qui minimisera la somme des carrés des erreurs, on utilise le code suivant (on fera une transformation des données, ici sur le simplexe, pour s’assurer que l’optimisation donnera un élément du simplexe)
erreur=function(a) sum((DONNEES[,1:4] %*% a-DONNEES[,5])^2) erreur2=function(b){ a=rep(NA,length(b)+1) a[1:length(b)]=exp(b)/(1+sum(exp(b))) a[length(b)+1]=1/(1+sum(exp(b))) erreur(a) } |
On obtient l’optimum, et on peut obtenir une prévision
b=optim(par=c(0,0,0),erreur2)$par a=rep(NA,length(b)+1) a[1:length(b)]=exp(b)/(1+sum(exp(b))) a[length(b)+1]=1/(1+sum(exp(b))) yp=DONNEES[,1:4] %*% a lines(1:111,yp,col="red",lwd=4) |
On peut alors faire changer le nombre d’observations utilisées
Techniquement, ca pose un petit soucis, car il faut s’assurer la convergence de deux algorithmes : celui d’estimation des paramètres des SARIMA, celui de recherche du paramètre de poids.
poids_optimaux=function(m){ Ym=ts(y[1:m],start=1996,frequency=52) TRY=try(model1 <- arima(Ym,order = c(4,0,0), seasonal = list(order = c(1,0,0))),FALSE) if(inherits(TRY, "try-error")) model1 <- arima(Ym,order = c(4,0,0), seasonal = list(order = c(1,0,0)),method="CSS") TRY=try(model2 <- arima(Ym,order = c(1,0,0), seasonal = list(order = c(1,1,0))),FALSE) if(inherits(TRY, "try-error")) model2 <- arima(Ym,order = c(1,0,0), seasonal = list(order = c(1,1,0)),method="CSS") TRY=try(model3 <- arima(Ym,order = c(0,0,2), seasonal = list(order = c(0,1,1))),FALSE) if(inherits(TRY, "try-error")) model3 <- arima(Ym,order = c(0,0,2), seasonal = list(order = c(0,1,1)),method="CSS") regL=lm(Temp~0+as.factor(NumWeek),data=electricite[1:m,]) E=residuals(regL) model4=arima(ts(E,start=1996,frequency=52),order = c(1,0,0)) DONNEES = data.frame( y1=predict(model1,n.ahead = 111)$pred, y2=predict(model2,n.ahead = 111)$pred, y3=predict(model3,n.ahead = 111)$pred, y4=predict(model4,n.ahead = 111)$pred+ predict(regL,newdata=electricite[m+(1:111),]), vraie=electricite[m+(1:111),"Temp"]) erreur=function(a) sum((as.matrix(DONNEES[,1:4]) %*% a-DONNEES[,5])^2) erreur2=function(b){ a=rep(NA,length(b)+1) a[1:length(b)]=exp(b)/(1+sum(exp(b))) a[length(b)+1]=1/(1+sum(exp(b))) erreur(a) } TRY=try(b <- optim(par=c(0,0,0),erreur2)$par, FALSE) a=rep(NA,length(b)+1) if(!inherits(TRY, "try-error")) { a[1:length(b)]=exp(b)/(1+sum(exp(b))) a[length(b)+1]=1/(1+sum(exp(b))) } return(a) } |
On peut alors regarder l’évolution des pondérations au cours du temps. On peut vectoriser, mais si l’algorithme aboutit sur une erreur, c’est plus simple d’identifier
VM=seq(200,480,by=10) VA=matrix(NA,length(VM),4) for(i in 1:length(VM)) VA[i,]=poids_optimaux(VM[i]) colnames(VA)=c("auto","SARMA","SARIMA","Buys Ballot") |
(pour débugger, c’est plus simple les boucles). On peut alors obtenir le graphique suivant, avec en rose Buys-Ballot, et en vert le modèle SARIMA, avec une racine unité saisonnière,
tva=t(VA) barplot(tva,legend = rownames(counts), col=c("yellow","orange","green","pink")) |
On se rend compte qu’assez vite, le modèle le plus intéressant est le modèle de Buys Ballot.
Il est possible de changer de fonction de perte, par exemple si on s’intéresse à un quantile à 90%,
tau=.9 loss=function(e) e*(tau-(e<=0)*1) |
Dans la fonction, on utilise alors
erreur=function(a) sum(loss(as.matrix( DONNEES[,1:4]) %*% a-DONNEES[,5])) |
C’est, cette fois, beaucoup plus incertain, même si les deux modèles les plus intéressants sont toujours le SARIMA et Buys-Ballot.
