{"word_count": 530, "character_count": 5696}
A partir de nos prévisions pour la température, on peut tenter de prévoir la consommation électrique. Rappelons que la série de consommation électrique ressemble à caplot(electricite[passe,"Load"],type="l")On peut tenter un modèle assez simple, où la consommation à la date [latex]Y_t[/latex] est fonction d'une tendance linéaire [latex]a+bt[/latex], de la position dans l'année (sous une forme non-paramétrique), de la température [latex]T_t[/latex] (sous une forme non-linéaire, disons un polynôme de degré 3), de la températeure en [latex]t-1[/latex], [latex]T_{t-1}[/latex] et de l'indice de production industriel [latex]IPI_t[/latex] (sous une forme linéaire pour faire simple).model1 = lm(Load~1+Time+as.factor(NumWeek)+poly(Temp,3)+Temp1+IPI,data=electricite[passe,])L'idée de la fonction polynomial pour l'impact de la température venait du graphique suivant (où on représente la série de la consommation, une fois enlevée la tendance linéaire)On peut aussi supposer une forme autorégressive, où [latex]Y_{t}[/latex] sera fonction de [latex]Y_{t-1}[/latex]model2 = lm(Load~1+Load1+Time+as.factor(NumWeek)+poly(Temp,3)+Temp1+IPI,data=electricite[passe,])On peut alors tenter de faire de la prévision. Le soucis avec ce second modèle, c'est la composante autorégressive. Pour prévoir [latex]Y_{t+h}[/latex], il faut utiliser la prévision faite en [latex]t+h-1[/latex], [latex]\hat Y_{t+h-1}[/latex].base_prevision=data.frame(Load_prev = NA,  Load1=c(electricite[length(passe),"Load"],rep(NA,110),Time = electricite[futur,"Time"],NumWeek = electricite[futur,"NumWeek"],Temp = MODELES[,"y4"], Temp1 = c(electricite[length(passe),"Temp"],MODELES[1:110,"y4"]),IPI = electricite[futur,"IPI"])for(t in 1:110){ p = predict(model2,newdata=base_prevision)[t] base_prevision[t,"Load_prev"] = p base_prevision[t+1,"Load1"] = p}On peut alors comparer la prévision [latex]\hat Y_{t}[/latex] avec l'observation [latex]Y_{t}[/latex].plot(1:111,electricite[futur,"Load"],type="l")lines(1:111,base_prevision[,"Load_prev"],col="red")On est bon les mois d'été (on capture bien le pic de la première quinzaine du mois d'août) par contre, on sous-estime fortement la consommation l'hiver. Et c'est gênant...Pour finir, on peut oublier le passage par les variables explicatives, et tenter directement un modèle de série temporelles.plot(electricite[passe,"Load"],type="l")Le soucis, c'est le côté hétéroscédastique de la série, avec une pente sur le minimum plus faible que sur le maximum.n=length(passe)plot(1:n,electricite[passe,"Load"],type="l")m=aggregate(electricite[passe,"Load"],by=list(as.factor(electricite[passe,"Year"])),min)M=aggregate(electricite[passe,"Load"],by=list(as.factor(electricite[passe,"Year"])),max)vm=which(electricite[passe,"Load"]%in%m$x)vM=which(electricite[passe,"Load"]%in%M$x)points(sort((1:n)[vm]),m$x,pch=19,col="blue")points(sort((1:n)[vM]),M$x,pch=19,col="blue")xm=((1:n)[vm])xM=((1:n)[vM])regm=lm(m$x~xm,col="blue")regM=lm(M$x~xM,col="blue")abline(regm,lty=2)abline(regM,lty=2)Une solution classique (et simple) est de passe au logarithmeplot(electricite[passe,"Load"],type="l",log="y")y=electricite[passe,"Load"]z = log(y)plot(z,type="l")B = data.frame(z=z,t=1:length(z))abline(lm(z~t,data=B),lty=2,col="red")Il faut ensuite enlever la tendance linéaire, pour stationnariser la sériez = residuals(lm(z~t,data=B))Z = ts(z,start=1996,frequency=52)modelz1 = arima(Z,order = c(0,0,0), seasonal = list(order = c(1,0,0)))modelz1 = arima(Z,order = c(4,0,0), seasonal = list(order = c(1,0,1)))Ce premier modèle est stationnaire, sans racine unité. On peut tenter d'introduire une racine unité saisonnièremodelz2 = arima(Z,order = c(0,0,0), seasonal = list(order = c(0,1,0)))modelz2 = arima(Z,order = c(1,0,3), seasonal = list(order = c(1,1,0)))Enfin, le dernier sera un peu plus simplemodelz3 = arima(Z,order = c(1,0,0), seasonal = list(order = c(2,0,0)))On stocke alors toutes les prévisions dans une baseDONNEES_u = data.frame(  u1 = predict(modelz1,n.ahead = 111)$pred,  u2 = predict(modelz2,n.ahead = 111)$pred,  u3 = predict(modelz3,n.ahead = 111)$pred)On rajoute alors la tendance linéaire à la prévision du résidureg = lm(z~t,data=B)dfutur = data.frame(t=futur)DONNEES_z = data.frame(  z1 = DONNEES_u$u1+predict(reg,newdata=dfutur),  z2 = DONNEES_u$u2+predict(reg,newdata=dfutur),  z3 = DONNEES_u$u3+predict(reg,newdata=dfutur)On a besoin de la variance pour la prévision. En effet, classiquement, la prévision sera l'espérance de notre grandeur d'intérêt. Ici, on a fait un modèle linéaire sur [latex]\log Y[/latex], i.e. [latex]\log Y \sim\mathcal{N}(\mu,\sigma^2)[/latex] et donc [latex]\mathbb{E}[Y]=\exp\left(\mu+\sigma^2/2\right)[/latex]sigma=summary(reg)$sigmaDONNEES_seu = data.frame(  seu1=sqrt(predict(modelz1,n.ahead = 111)$se^2+sigma^2),  seu2=sqrt(predict(modelz2,n.ahead = 111)$se^2+sigma^2),  seu3=sqrt(predict(modelz3,n.ahead = 111)$se^2+sigma^2))On fait ici l'hypothèse que les estimateurs des prévisions des deux modèles (la tendance linéaire et celle du modèle autorégressif) sont indépendants et donc on peut sommer les termes de variance. Aussi, les prévisions pour [latex]Y[/latex] sontDONNEES_y = data.frame(  y1 = exp(DONNEES_z$z1+1/2*DONNEES_seu$seu1^2),  y2 = exp(DONNEES_z$z2+1/2*DONNEES_seu$seu2^2),  y3 = exp(DONNEES_z$z3+1/2*DONNEES_seu$seu3^2))On peut représenter les prévisions des trois modèles (avec les observations)plot(futur,electricite[futur,"Load"],type="l")lines(futur,DONNEES_y$y1,col="red")lines(futur,DONNEES_y$y2,col="blue")lines(futur,DONNEES_y$y3,col="purple")On peut comparer avec la prévision précédente,lines(futur,base_prevision[,"Load_prev"],col="orange")On y perd un peu l'été, par contre, on s'est un peu amélioré l'hiver.