R Notebook - PROJET 9 - PEYRONNE David

—————————————————————————————-

Rappel du projet

Prédire la consommation électrique

—————————————————————————————-

Données

Analyses

Afin d’introduire votre analyse, effectuez une brève description des données (analyses univariées et bivariées).

Bilan

Questions / Réponses

Données

—————————

1 Importation des données

—————————

1.1 Sources des données

• -> Jeu de données fichier xlsx

1.2 Installation des packages

• -> L’installation des packages ne se fait qu’une seule fois, je les ai commentées pour éviter de réisntaller à chaque chargement du script *

#install.packages(c("Factoshiny"))
#install.packages("ggpubr")
#install.packages("rJava")
#install.packages("xlsx") 
#install.packages("tseries")
#install.packages("dataseries")
#install.packages("tidyverse")
#install.packages("caschrono")

1.3 Installation des libraries

• -> L’activation des librairies se fait à chaque démarrage du script*
  
• -> corrplot doit être appelé après ggplot2 pour éviter un warning*

1.4 Configuration des variables pour chaque fichier

• -> définition d’une variable répertoire pour faciliter le changement d’OS *

repertoire_sources <- "/home/user/Documents/formations/oc/Formation Analyst/P9/sources/"
repertoire_images <- "/home/user/Documents/formations/oc/Formation Analyst/P9/images/"
fichier <- paste (repertoire_sources,"tab.csv",sep = "")

1.5 Création des Data Frame par lecture du fichier températures

• -> les critères d’importation, séparateurs, virgules décimales, entête etc, ont étés définis après observation des fichiers sources à l’aide d’une éditeur de texte ou d’un tableur*

data <- read.table(fichier, header=TRUE, sep=',', fileEncoding = "ISO-8859-15")
head(data)

valeurs maxi mini de la période

data_max=max(data$Mois)
data_min=min(data$Mois)
cat("La série de données est entre : " , data_min,"et : " , data_max)

On garde les consommations de la France pour les données et

On vérifie qu’il n’y ait pas de valeurs nulles pour la consommation

df_fr <- subset(data, Territoire == "France")
df_fr <- df_fr[,c("Mois","Consommation.totale")]
sapply(df_fr, function(x) sum(is.na(x)))

On converti la variable Mois en Date

df_fr$Mois<-as.Date(as.yearmon(df_fr$Mois, "%Y-%m"))

affichage de la série temporelle

df_fr %>%
  ggplot(aes(x = Mois, y = `Consommation.totale`)) +
  geom_line()
  ggsave("img/01 - consommationTotale.png", width = 7.29, height = 4.5)

• -> On peut voir une saisonnalité de la consommation.
• -> Un pic de consommation important en période hivernale lié au chauffage. Et un plus petit pic en plein été lié à l’utilisation de la climatisation. On peut supposer qu’il s’agit d’une série temporelle additive.

1.6 Création des Data Frame DJU par lecture du fichier

fichier_dju <- paste (repertoire_sources,"dju.csv",sep = "")
dju <- read.csv2(fichier_dju,header = TRUE,sep = ",")

conversion du champ Mois au format date

dju$Mois<-as.Date(as.yearmon(dju$Mois, "%Y-%m"))

affichage de la série temporelle

dju %>%
  ggplot(aes(x = Mois, y = `DJU`)) +
  geom_line()
  ggsave("img/02 - DJU.png", width = 7.29, height = 4.5)

• -> On constate qu’il y a la même saisonnarité que la consommation

1.7 Combiner les deux tableaux de données pour les dates communes

on observe sur le graphique précédent que les dernières valeurs sont à zéro

tail(dju,10)

on vérifie pour les consommations :

tail(df,10)

On restrain le jeu de données entre janvier 2012 et mai 2020 inclus

On combine les deux tableaux de données

On renomme les colonnes

df_fr <- subset(df_fr, Mois < "2020-06-01")
dju <- subset(dju, Mois < "2020-06-01")
df <- cbind( df_fr, dju$DJU) 
names(df)[names(df) == 'Consommation.totale'] <- 'conso'
names(df)[names(df) == 'dju$DJU'] <- 'dju'
names(df)[names(df) == 'Mois'] <- 'mois'
dju %>%
  ggplot(aes(x = Mois, y = `DJU`)) +
  geom_line()

1.8 Données “centré réduit” pour superposer les courbes

df_scaled <- data.frame(df, scale(df[,2:3]))
ggplot(df_scaled, aes(x = mois)) +
  geom_line(aes(y = conso.1, color = "conso")) +
  geom_line(aes(y = dju.1, color = "dju")) 
  ggsave("img/03 - consommation et DJU.png", width = 7.29, height = 4.5)

• -> les courbes se superposent très bien à quelques exceptions près tels que les maximums, et l’usage de la climatisation car le dju donne l’écart de température par rapport au froid ce qui laisse supposer les besoins en chauffage alors que les écarts en été sont dus à la climatisation

—————————

2 Régression linéaire

—————————

plot des données

plot(conso~dju, data = df)

• -> on y voit un diagonale se dessiner on peut donc envisager une régression linéaire

2.1 Plot avec la régression linéaire

model <- lm(conso~dju, data = df)
plot(conso~dju, data = df) + abline(model, col = "red")

Image sauvegarde du graphique

png("img/04 - regressionLineaire.png", width = 7.25, height = 4.5, units = 'in', res = 300)
plot(conso~dju, data = df) + abline(model, col = "red")

Summary du modèle de la régression

summary(model)

ANALYSE du modèle

• -> les *** ( trois étoiles) montrent les indices les plus pertinents
• -> R² à 0.95 indique qu’il y a une forte corrélation entre la consommation et la DJU
• -> la P-Value < 5% indique que les coefficients sont significatifs
• -> nous récupérons le cofficient de DJU pour corriger la consommation par rapport au DJU

model$coefficients[2]

2.2 Consommation corrigée par la regression linéaire

df$corrige<-df$conso-model$coefficients[2]*df$dju
ggplot(df, aes(x = mois)) +
  geom_line(aes(y = conso, color = "réelle")) +
  geom_line(aes(y = corrige, color = "corrigée"))
 ggsave("img/05 - consoReelleetCorrigee.png", width = 7.29, height = 4.5)

• -> on observe que malgré la correction on peut y voir une saisonnalité.
• -> la régression ne retire pas la saisonnalité de nos données
• -> elle a permis de mettre en évidence les résidus

2.3 Résidus de la régression linéaire

res<-resid(model)
plot(df$mois,res,main="Résidus") + abline (v = as.Date("2012-08-01") , col = "blue") + abline (v = as.Date("2013-08-01") , col = "blue") + abline (v = as.Date("2014-08-01") , col = "blue")  + abline (v = as.Date("2015-08-01") , col = "blue") - abline(h=0,col="red")

• -> les lignes bleues sont les mois d’août
• -> Les résidus positifs les plus importants apparaissent en été, avec la différence de consommation dûe au climatiseur qui n’est pas prise en compte avec les calculs du chauffage.

—————————

3 Saisonnalité

—————————

• -> pour désaisonnaliser la série, on va la décomposer

consommation <- ts(df$conso, start=c(2012, 1), end=c(2020, 4), frequency = 12)
consommation_dec <- decompose(consommation) #,'additive') je ne vois pas la différence avec ou cette valeur.
plot(consommation_dec)

png("img/05b - boxPlotConsommationsMensuelles.png", width = 7.25, height = 4.5, units = 'in', res = 300)
boxplot(consommation~ cycle(consommation))

boxplot(consommation~ cycle(consommation))

• -> ce boxplot montre sur l’ensemble des années la consommation pour chaque mois .
• -> on y voit une forte consommation entre les mois 11 et 3 et faibles en été avec un petit rebond en juillet ( climatisation )
• -> certains mois ont des écarts importants janvier / février, d’autres sont très réguliers tels que Juin / Juillet

corrigee <- ts(df$corrige, start=c(2012, 1), end=c(2020, 4), frequency = 12)
corrigee_dec <- decompose(corrigee) #,'additive') je ne vois pas la différence avec ou cette valeur.

png("img/05C - decompositionConsoCorrigee.png", width = 7.25, height = 4.5, units = 'in', res = 300)
plot(corrigee_dec)

plot(corrigee_dec)

• -> la tendance à la baisse qui n’était pas perceptible sur le graphe des consommations
• -> la saisonnalité est confirmée
• -> la valeur résiduelle est représentée

Désaisonnalisation - consommation corrigée des variations saisonnières

• -> La décomposition permet d’obtenir et de retirer l’aspect saisonnier d’une série temporelle ## 3.1 Saisonnalité sur la consommation

consommation_cvs = consommation-consommation_dec$seasonal

png("img/06 - consoDessainalisee.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(consommation, consommation_cvs, xlab="t", ylab="consommations",col=c(1,2),lwd=c(1,2))
title("Consommation réelle sans les variations saisonnières")

ts.plot(consommation, consommation_cvs, xlab="t", ylab="consommations",col=c(1,2),lwd=c(1,2))
title("Consommation réelle sans les variations saisonnières")

3.2 saisonnalité sur la consommation corrigée ( par les DJU )

corrigee_cvs = corrigee-corrigee_dec$seasonal

png("img/07 - consoCorrigeeDessainalisee.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(corrigee, corrigee_cvs, xlab="temps", ylab="consommation",col=c(1,2),lwd=c(1,2))
legend("bottomleft",legend=c("corrigée","désaisonnalisée"),col=c(1,2),lwd=c(1,2))
title("Consommation corrigée des variations saisonnières")

ts.plot(corrigee, corrigee_cvs, xlab="temps", ylab="consommation",col=c(1,2),lwd=c(1,2))
legend("bottomleft",legend=c("corrigée","désaisonnalisée"),col=c(1,2),lwd=c(1,2))
title("Consommation corrigée des variations saisonnières")

le test de Ljung-Box nous permet de confirmer la blancheur des résidus. Ils sont donc indépendants du temps si les p-value de chaque lag est >5%

Box.test.2(corrigee_dec$random,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> le test de blancheur n’est pas rejeté donc le modèle que l’on vient de trouver est correct

test de la normalité des résidus avec la méthode Shapiro

shapiro.test(corrigee_dec$random)

• -> p-value > 5% donc les résidus respectent une loi normale

—————————

4 PREDICTION - HOLT WINTERS

—————————

prédiction méthode de Holt Winters (lissage exponentiel) sur la série corrigée modèle AAA pour prendre en compte l’erreur, la tendance et la saisonnalité

4.1 Holts-Winters - application de la méthode

howi=ets(corrigee,model="AAA")
summary(howi)

• -> les smooting parameters aident à déterminer si la prédiction se fait par rapport aux dernières valeurs ou à des valeurs plus anciennes. Plus on est proche de 1 plus les dernières valeurs comptent pour prédire les nouvelles. Plus on s’approche de zéro plus il faudra se baser sur des valeurs anciennes.

4.2 Analyse a posteriori Holt Winters pour vérifier l’efficacité de ce modèle sur un tronc connu

On tronque la série de la période de janvier 2019 à janvier 2020, et on la prévoit sur les données précédentes

x_tronc=window(corrigee,end=c(2019,4))
x_a_prevoir=window(corrigee,start=c(2019,5))

On applique la méthode sur cette série tronquée

modelhowi_tronc=ets(x_tronc,model="AAA")
summary(modelhowi_tronc)

Estimation de la prédiction et de l’intervalle de 95% en plus et en moins

pred_modelhowi_tronc=forecast(modelhowi_tronc,h=12,level=95)
predhowi_tronc=pred_modelhowi_tronc$mean
predhowi_l_tronc=ts(pred_modelhowi_tronc$lower,start=c(2019,5),frequency=12)
predhowi_u_tronc=ts(pred_modelhowi_tronc$upper,start=c(2019,5),frequency=12)

affichage du résultat

png("img/08 - analysePosterioriHW.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(x_a_prevoir,predhowi_tronc,predhowi_l_tronc,predhowi_u_tronc,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(3,3,2,2))
legend("topleft",legend=c("X","X_prev"),col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Analyse a posteriori : Holt Winters")

ts.plot(x_a_prevoir,predhowi_tronc,predhowi_l_tronc,predhowi_u_tronc,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(3,3,2,2))
legend("topleft",legend=c("X","X_prev"),col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Analyse a posteriori : Holt Winters")

• -> La prédiction est bien contenue dans l’intervealle des 95%
• -> L’erreur de prédiction finale est due à un événément exceptionnel, température élevée pour un mois de janvier donc arrêt prématuré des chauffages electriques

4.3 Les résidus

Box.test.2(pred_modelhowi_tronc$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

les rédisus sont indépendants, il s’agit d’un bruit blanc, le modèle que l’on vient de trouver est correct

Normalité des résidus - méthode Shapiro

shapiro.test(pred_modelhowi_tronc$residuals)

Les résidus suivent une loi normale car p-value > 5%

4.4 Précision : RMSE et MAPE

rmse_howi=sqrt(mean((x_a_prevoir-predhowi_tronc)^2))
mape_howi=mean(abs(1-predhowi_tronc/x_a_prevoir))*100
cat('RMSE Holt Winters : ',rmse_howi,'\n')
cat('MAPE Hot Winters : ',mape_howi,'\n')

4.5 Prévision à l’aide du modèle Holt Winters sur un an

pred_model_howi=forecast(howi,h=12,level=95) # prédiction modèle HW sur séries de 12 avec un niveau de 5%
pred_howi=pred_model_howi$mean # moyenne 
pred_howi_l=ts(pred_model_howi$lower,start=c(2020,5),frequency=12) # estimation basse de l'intervalle
pred_howi_u=ts(pred_model_howi$upper,start=c(2020,5),frequency=12) # estimation haute

sauvegarde image

png("img/09 - previsionHW.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(corrigee,pred_howi,pred_howi_l,pred_howi_u,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : Holt Winters")

ts.plot(corrigee,pred_howi,pred_howi_l,pred_howi_u,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : Holt Winters")

Enregistrement en image

png("img/10 - previsionHW_ZOOM.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(window(corrigee,start=c(2019,5)),pred_howi,pred_howi_l,pred_howi_u,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : Holt Winters - ZOOM ")

ts.plot(window(corrigee,start=c(2019,5)),pred_howi,pred_howi_l,pred_howi_u,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : Holt Winters - ZOOM ")

—————————

5 Prédiction - SARIMA

—————————

méthode SARIMA - Arima sur série avec Saisonnarité : Auto Regressive Integration Mean Average donne les valeurs p / d / q P/D/Q

Cette methode s’applique sur les séries stationnaires, si on peut le supposer visuellement, on peut utiliser une fonction qui va donner un résultat ( pour pouvoir l’intégrer dans des algorythmes d’analyse par exemple )

5.1 STATIONARITÉ

Pour confirmer que la série n’est pas stationnaire on peut utiliser ADF Augmented Dickey-Fuller

# la série est-elle stationnaire ?^
stationnaire_test<-adf.test(corrigee,k=12)
stationnaire_test

Analyse : l’hypothèse H0 est que la série est stationnaire.

or la p-value est supérieure à 5% donc on rejette H0 et on retient l’Hypothèse alternative Ha, la série n’est pas stationnaire.

il faut en premier stationnariser la série temporelle par différenciation

5.2 Stationnarisation

Autocorrélogramme Simple

acf=acf(corrigee,lag.max=36,ylim=c(-1,1))
acf$lag=acf$lag*12

sauvegarde de l’image

png("img/11 - ACF01 .png", width = 7.25, height = 4.5, units = 'in', res = 300)
acf(corrigee,lag.max=36,ylim=c(-1,1))

la décroissance est très lente, ce qui empêche cette sortie ACF d’être un autocorrélogramme simple

les traits bleus pointillés sont les limites de la significativité à 95%

5.2.1 Différenciation avec une saisonnalité de 12

y_dif_1_12=diff(corrigee,lag = 12,differences = 1)
acf=acf(y_dif_1_12,lag.max=36,ylim=c(-1,1), plot=FALSE)
acf$lag=acf$lag*12
plot(acf)

sauvegarde de l’image

png("img/12 - ACF02-différenciée-Lag12 .png", width = 7.25, height = 4.5, units = 'in', res = 300)
acf(y_dif_1_12,lag.max=36,ylim=c(-1,1))

###€ La sortiE ACF de la série différenciée semble pouvoir être interprétée comme un autocorrélogramme simple empirique. On ne voit plus les décroissances lentes précédentes. La série est donc potentiellement stationnaire. On identifiera donc un modèle SARIMA . ## 5.2.2 ACF - La série est-elle stationnaire ?

stationnaire_test<-adf.test(y_dif_1_12)
stationnaire_test

p-value de l’adf < 5% donc série potentiellement stationnaire

5.2.3 PACF ( partialy auto correlation function)

pacf=pacf(y_dif_1_12,lag.max=36,ylim=c(-1,1), plot=FALSE)
pacf$lag=pacf$lag*12
plot(pacf,ylim = c(-1,1))

IMAGE sauvegarde de l’image

png("img/13 - PACF01-différenciée-Lag12 .png", width = 7.25, height = 4.5, units = 'in', res = 300)
pacf(y_dif_1_12,lag.max=36,ylim=c(-1,1))

La série temporelle étant potentiellement stationnaire, on passe à l’identification de modèle

5.3 RECHERCHE D’UN MODÈLE SARIMA POUR LA PRÉDICTION

On va rechercher de façon empirique quelle sera la meilleure combinaison des coefficients pour le modèle

SARIMA (p,d,q)(P,D,Q)[12] pour une saisonnalité de 12

p pour AR, P pour SAR, d pour la tendance, D pour la saisonnalité, q pour MA et Q pour SMA

deux fonctions nous aident à différencier et choisir les coefficients du modèle pour d et D

d = ndiffs(corrigee)
D = nsdiffs(corrigee)
cat('d = ',d,' - différentiation en tendances / D = ', D,'différenciations en saisonnalité')

on va travailler sur la série corrigée et on y applique les paramètres que nous venons de trouver et nous appliquerons la methode CSS-ML ( maximum de vraisemblance conditionnelle )

5.3.1 SARIMA (000)(111)

————————–

5.3.1.1 Modélisation

model1=Arima(corrigee,order=c(0,0,0),list(order=c(1,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model1)

5.3.1.2 Coefficients

t_stat(model1)

• P-valeurs > 5% : Les coefficients ne sont pas significatifs

5.3.1.3 Test de blancheur - Ljung-Box

Box.test.2(model1$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

5.3.1.4 Normalité des résidus - méthode Shapiro

shapiro.test(model1$residuals)

5.3.2 SARIMA (001)(111)

————————-

5.3.2.1 Modélisation

model2=Arima(corrigee,order=c(0,0,1),list(order=c(1,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model2)

5.3.2.2 Coefficients

t_stat(model2)

• -> P-valeurs > 5% : Les coefficients ne sont pas significatifs sauf sma1

5.3.2.3 Test de blancheur - Ljung-Box

Box.test.2(model2$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

5.3.2.4 Normalité des résidus - Shapiro

shapiro.test(model2$residuals)

5.3.3 SARIMA (101)(111)

—————–

5.3.3.1 Modélisation

model3=Arima(corrigee,order=c(1,0,1),list(order=c(1,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model3)

5.3.3.2 Coefficients

t_stat(model3)

• -> P-valeurs > 5% : Les coefficients ne sont pas significatifs sauf sma1

5.3.3.3 Test de blancheur - Ljung-Box

Box.test.2(model3$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

5.3.3.4 Normalité des résidus - Shapiro

shapiro.test(model3$residuals)

• -> la p-valeur étant inférieure à 5% l’hypothèse de normalité est rejetée.

5.3.4 SARIMA (011)(011)

—————–

5.3.4.1 Modélisation

model4=Arima(corrigee,order=c(0,1,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model4)

5.3.4.2 Coefficients

t_stat(model4)

• -> P-valeurs > 5% : Les coefficients sont significatifs

5.3.4.3 Test de blancheur - Ljung-Box

Box.test.2(model4$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

5.3.4.4 Normalité des résidus - Shapiro

shapiro.test(model4$residuals)

• -> la p-valeur étant supérieure à 5% l’hypothèse de normalité n’est pas rejetée, donc les résidus sont potentiellement gaussiens

5.3.5 SARIMA (101)(011)

—————–

5.3.5.1 Modélisation

model5=Arima(corrigee,order=c(1,0,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model5)

5.3.5.2 Coefficients

t_stat(model5)

• -> P-valeurs > 5% : Les coefficients sont significatifs

5.3.5.3 Test de blancheur - Ljung-Box

Box.test.2(model5$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

5.3.5.4 Normalité des résidus - Shapiro

shapiro.test(model5$residuals)

—————————

6 SARIMA AUTO

—————————

6.1 Modélisation

model_auto <- auto.arima(corrigee, stepwise=TRUE, approximation=TRUE, allowdrift=FALSE)

6.2 Coefficients

t_stat(model_auto)

• -> P-valeurs > 5% : Les coefficients ne sont pas significatifs sauf sma1

6.3 Test de blancheur - Ljung-Box

Box.test.2(model_auto$residuals,nlag = c(6,12,18,24,30,36), type = "Ljung-Box",decim = 5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

6.4 Normalité des résidus - Shapiro

shapiro.test(model_auto$residuals)

—————————

7 Comparaison des modèles SARIMA 4 et 5 sur la période tronquée

—————————

7.1 MODÈLE 5

7.1.1 Modélisation

model5tronc=Arima(x_tronc,order=c(1,0,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model5tronc)

7.1.2 Coefficients

t_stat(model5tronc)

seul SMA1 est significatif

7.1.3 Test de blancheur - Ljung-Box

Box.test.2(model5tronc$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

7.1.4 Normalité des résidus - Shapiro

shapiro.test(model5tronc$residuals)

• -> la p-valeur étant supérieure à 5% l’hypothèse de normalité n’est pas rejetée, donc les résidus sont potentiellement gaussiens

7.2 MODÈLE 4

7.2.1 Modélisation

model4tronc=Arima(x_tronc,order=c(0,1,1),list(order=c(0,1,1),period=12),include.mean=FALSE,method="CSS-ML")
summary(model3tronc)

7.2.2 Coefficients

t_stat(model4tronc)

les coefficients sont significatifs

7.2.3 Test de blancheur - Ljung-Box

Box.test.2(model4tronc$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)

• -> P-valeurs > 5% : les résidus sont un bruit blanc

7.2.4 Normalité des résidus - Shapiro

shapiro.test(model4tronc$residuals)

• -> la p-valeur étant supérieure à 5% l’hypothèse de normalité n’est pas rejetée, donc les résidus sont potentiellement gaussiens

7.3 CHOIX DU MODÈLE SARIMA : modèle 4

le modèle retenu est le modèle 4 SARIMA (0,1,1 ) (0,1,1) [12]

• pour lequel les coefficients sont significatifs
• les résidus sont des bruits blancs

7.4 PRECISION DU MODÈLE SUR INTERVALLE CONNU

pred_model4tronc=forecast(model4tronc,h=12,level=95)
pred4_tronc=pred_model4tronc$mean
pred4_l_tronc=ts(pred_model4tronc$lower,start=c(2019,5),frequency=12)
pred4_u_tronc=ts(pred_model4tronc$upper,start=c(2019,5),frequency=12)
ts.plot(window(x_a_prevoir),pred4_tronc,pred4_l_tronc,pred4_u_tronc,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(3,3,2,2))
legend("topleft",legend=c("X","X_prev"),col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Analyse a posteriori : SARIMA ZOOM modèle 4 (0,1,1)(0,1,1)[12]")

IMAGE - sauvegarde du graphique

png("img/14 - Analyse a posteriori SARIMA ZOOM modèle 4.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(window(x_a_prevoir),pred4_tronc,pred4_l_tronc,pred4_u_tronc,xlab="t",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(3,3,2,2))
legend("topleft",legend=c("X","X_prev"),col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Analyse a posteriori : SARIMA ZOOM modèle 4 (0,1,1)(0,1,1)[12]")

Précision du modèle 4 avec RMSE et MAPE

rmse_sarima=sqrt(mean((x_a_prevoir-pred4_tronc)^2))
mape_sarima=mean(abs(1-pred4_tronc/x_a_prevoir))*100
cat('RMSE SARIMA : ',rmse_sarima,'\n')
cat('MAPE SARIMA : ',mape_sarima,'%\n')

7.5 PRÉVISION SUR UN AN

pred_model4=forecast(model4,h=12,level=95) # prédiction modèl3 sur séries de 12 avec un niveau de 5%
pred4=pred_model4$mean
pred4_l=ts(pred_model4$lower,start=c(2020,5),frequency=12)
pred4_u=ts(pred_model4$upper,start=c(2020,5),frequency=12) # estimation haute 
ts.plot(corrigee,pred4,pred4_l,pred4_u,xlab="temps",ylab="Consommation corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : SARIMA")

ts.plot(window(corrigee,start=c(2018,5)),pred4,pred4_l,pred4_u,xlab="temps",ylab="Conso corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : SARIMA - Zoom ")

png("img/15 - Prédiction SARIMA ZOOM avec modèle 4.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(window(corrigee,start=c(2018,5)),pred4,pred4_l,pred4_u,xlab="temps",ylab="Conso corrigée",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
legend("topleft",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
title("Prévision : SARIMA - Zoom ")

—————————

8 BILAN - Comparaison des modèles les plus performants holt-winters et SARIMA

—————————

8.1 Analyses a postériori superposées

ts.plot(window(x_a_prevoir),pred4_tronc,predhowi_tronc,xlab="t",ylab="Consommation corrigée",col=c(1,2,3),lty=c(1,1,1),lwd=c(3,3,3))
legend("topleft",legend=c("Réelle","SARIMA"),col=c(1,2),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("Holt Winters"),col=c(3),lty=c(1),lwd=c(3))
title("Analyse a posteriori : SARIMA + Holt Winters ")

png("img/16 - Analyse Posteriori SARIMA et Holt-Winters.png", width = 7.25, height = 4.5, units = 'in', res = 300)
ts.plot(window(x_a_prevoir),pred4_tronc,predhowi_tronc,xlab="t",ylab="Consommation corrigée",col=c(1,2,3),lty=c(1,1,1),lwd=c(3,3,3))
legend("topleft",legend=c("Réelle","SARIMA"),col=c(1,2),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("Holt Winters"),col=c(3),lty=c(1),lwd=c(3))
title("Analyse a posteriori : SARIMA + Holt Winters ")

8.2 Précision de la prédiction avec RMSE et MAPE

RMSE est dépendant de l’échelle, MAPE indice en pourcentage

cat('RMSE Holt Winters : ',rmse_howi,'\n')
cat('RMSE SARIMA : ',rmse_sarima,'\n')
cat('MAPE Hot Winters : ',mape_howi,'\n')
cat('MAPE SARIMA : ',mape_sarima,'\n')

• -> Ces indicateurs mettent en évidence les écarts par rapport à la réalité, plus ils sont faibles, plus ils sont proches des vrais valeurs et donc ils sont meilleurs.
• -> Ainsi, la méthode la meilleure est, pour ce cas, SARIMA

—————————

MODÈLE RETENU SARIMA (0.1.1)(0.1.1)[12]

—————————