{"word_count": 561, "character_count": 5900}
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 [latex]SARIMA[/latex] (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, [latex]\hat{y}_t = \sum_{i}\omega_i \hat{y}_t^{(i)}[/latex] où [latex]\omega=(\omega_i)[/latex] est un point du simplexe, au sens où [latex]\omega_i\geq 0[/latex] et [latex]\omega_1+...+\omega_k=1[/latex]. On veut ensuite trouver les "meilleurs" pondérations, par exemple au sens de la norme [latex]L^2[/latex]. On va calibrer nos quatre modèles sur les [latex]m[/latex] premières valeurs, et on compare ensuite la combinaison des prévions sur les 111 prochaines valeurs (et les vraies valeurs), [latex display="true"]\omega^\star=\text{argmin}\{\sum_{h=1}^{111}\left(y_{t+h}- \sum_{i}\omega_i \hat{y}_{t+h}^{(i)}\right)^2\}[/latex]Par exemple, on peut utiliser les 200 premières valeurs.m=200Ym=ts(y[1:m],start=1996,frequency=52)On va ensuite estimer les 4 modèles sur ces 200 valeursmodel1=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édictionsDONNEES = 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 utilisea=rep(1/4,4)y_prev = as.matrix(DONNEES[,1:4]) %*% aOn peut ainsi visualiser les 4 prévisions, leur combinaison linéaire (ici avec des équi-pondérations), et leurs observationsplot(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évisionb=optim(par=c(0,0,0),erreur2)$para=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] %*% alines(1:111,yp,col="red",lwd=4)On peut alors faire changer le nombre d'observations utiliséesTechniquement, ca pose un petit soucis, car il faut s'assurer la convergence de deux algorithmes : celui d'estimation des paramètres des [latex]SARIMA[/latex], celui de recherche du paramètre de poids.poids_optimaux=function(m){  Ym=ts(y[1:m],start=1996,frequency=52)  TRY=try(model1 &#60;- arima(Ym,order = c(4,0,0),     seasonal = list(order = c(1,0,0))),FALSE)  if(inherits(TRY, "try-error")) model1 &#60;- arima(Ym,order = c(4,0,0),     seasonal = list(order = c(1,0,0)),method="CSS")  TRY=try(model2 &#60;- arima(Ym,order = c(1,0,0),     seasonal = list(order = c(1,1,0))),FALSE)  if(inherits(TRY, "try-error")) model2 &#60;- arima(Ym,order = c(1,0,0),     seasonal = list(order = c(1,1,0)),method="CSS")  TRY=try(model3 &#60;- arima(Ym,order = c(0,0,2),     seasonal = list(order = c(0,1,1))),FALSE)  if(inherits(TRY, "try-error")) model3 &#60;- 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 &#60;- 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'identifierVM=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 [latex]SARIMA[/latex], 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=.9loss=function(e) e*(tau-(e&#60;=0)*1)Dans la fonction, on utilise alorserreur=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 [latex]SARIMA[/latex] et Buys-Ballot.