Projet 9 - Prédire la demande d'électricité
Table des matières
- Paramètres
- 1. Analyses des données
- 2. Correction de l'effet température sur la consommation
- 3. Désaisonnalisation de la consommation corrigée
- 4. Prévision de la consommation corrigée passée
Paramètres¶
In [1]:
import warnings
def fxn():
warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
fxn()
# importations de modules "local"
import sys
sys.path.append('functions/')
import base_notebook
import css_base_dark
import fonctions_perso
In [2]:
# import des librairies et des paramètres
from base_notebook import *
from css_regression import *
from css_base_dark import *
# creation dossier grahs
path_graph = "data/exports/graphs/"
os.makedirs(path_graph, exist_ok=True)
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 100
# repère de début de traitement pour calculer le temps total d'exécution du notebook
start_time_m4 = time.time()
In [3]:
# Import des données de production et de consommation d'électricité
df = pd.read_csv("data/inputs/eCO2mix_RTE_energie_M.csv", index_col=0, parse_dates=True, sep='\t', engine='python', encoding='ANSI')
In [4]:
# affichage aléatoires de 3 observations
df.sample(3)
# affichage de la taille du df
df.shape
# affichage des colonnes
df.columns
# affichage de l'index
df.index
Out[4]:
Out[4]:
Out[4]:
Out[4]:
Analyse du daraframe
- on va s'intéresser à la demande en electricité au niveau national :
- on limite les données au
TerritoireFrance - on limite les colonnes à la seule colonne "Consommation Totale" (que l'on renomme en 'conso' - l'unité est le GWh)
- on limite les données au
- on renomme l'index en 'date'
- on va s'intéresser à la demande en electricité au niveau national :
In [5]:
df_fr = df.loc[df['Territoire'] == 'France', ['Consommation totale']]
df_fr.columns = ['conso']
df_fr=df_fr.rename_axis('date')
df_fr
Out[5]:
In [6]:
# vérification de présence de valeurs nulles
df_fr.isna().sum()
Out[6]:
In [7]:
# Visualisons la série temporelle de la consommation d'électricité en France
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso'])
plt.title("Consommation totale d'électricité en France (en GWh)", y=1.01)
plt.ylabel("Consommation (en GWh)")
plt.tight_layout()
plt.savefig(path_graph+"01-serie_temporelle_conso_France_brute.png", dpi=200)
plt.show();
Ainsi, on constate que :
- la tendance semble stable au cours du temps
- la saisonnalité est annuelle
In [8]:
# calculons le minimum et le maximum de la consommation pour chaque année (en indiquant également le mois)
# années de la periode
# on limite les données aux années complètes (on ne compte pas 2021)
years_list = np.arange(2012, 2021, 1).tolist()
# création d'un df df_fr_year avec pour index les années
df_fr_year=pd.DataFrame(
index=years_list,
columns=["avg", "min", "min_m", "max", "max_m"])
# création d'un df groupby par année
df_fr2=df_fr.copy()
df_fr2=df_fr2[df_fr2.index < "2021-01.-01"]
df_fr2['year']=pd.DatetimeIndex(df_fr2.index).year
df_fr2["month"] = pd.DatetimeIndex(df_fr2.index).month
df_fr3=df_fr2.groupby(['year'])
# initialisation des listes avec valeurs et numéro du mois
avg_list_value=[]
sum_list_value=[]
min_list_value=[]
min_list_month=[]
max_list_value=[]
max_list_month=[]
# on boucle sur chaque année et on récupère les données dans les listes
for name, group in df_fr3:
avg_list_value.append(group['conso'].mean())
sum_list_value.append(group['conso'].sum())
min_list_value.append(group['conso'].min())
min_list_month.append(group.loc[group['conso'] == group['conso'].min(), "month"].values[0])
max_list_value.append(group['conso'].max())
max_list_month.append(group.loc[group['conso'] == group['conso'].max(), "month"].values[0])
# on remplit df_fr_year avec les listes précédentes
df_fr_year['avg'] = avg_list_value
df_fr_year['avg'] = df_fr_year['avg'].apply(lambda x: round(x))
df_fr_year['conso_annee'] = sum_list_value
df_fr_year['conso_annee'] = df_fr_year['conso_annee'].apply(lambda x: "{:,}".format(x).replace(',', ' '))
df_fr_year['min'] = min_list_value
df_fr_year['min_m'] = min_list_month
df_fr_year['max'] = max_list_value
df_fr_year['max_m'] = max_list_month
df_fr_year
Out[8]:
- Représentation graphique de la saisonnalité
In [9]:
# prepa des données
years = df_fr_year.index
df_fr['year']=pd.DatetimeIndex(df_fr.index).year
df_fr["month"] = pd.DatetimeIndex(df_fr.index).month
months_chiffre = np.arange(0,12,1).tolist()
months=["Janv", "Fev", "Mars", "Avril", "Mai", "Juin", "Juil", "Aout", "Sept", "Oct", "Nov", "Dec"]
df_fr["month"]=df_fr["month"].apply(lambda x: months[x-1])
df_fr.dtypes
# prepa des couleurs
np.random.seed(100)
mycolors = np.random.choice(list(mpl.colors.XKCD_COLORS.keys()), len(years), replace=False)
# Dessiner le graphique
plt.figure(figsize=(12,6), dpi= 100)
for i, y in enumerate(years) :
if i > 0 :
plt.plot('month', 'conso', data=df_fr.loc[df_fr.year==y, :], color=mycolors[i], label=y)
plt.text(df_fr.loc[df_fr.year==y, :].shape[0]-.9, df_fr.loc[df_fr.year==y, 'conso'][-1 :].values[0], y, fontsize=8, color=mycolors[i])
#Décoration
plt.gca().set(xlim=(-0.3, 12), ylim=(25000, 60000))
plt.ylabel("Consommation (en GWh)")
plt.yticks(fontsize=11, alpha=.7)
plt.title("Saisonnalité de la consommation brute d'électricité en France")
plt.tight_layout()
plt.savefig(path_graph+"02-saisonnalite_consommation_electricite_brute.png", dpi=200)
plt.show();
- On constate que la consommation brute d'électricité suit approximativement le même schéma chaque année : baisse jusqu'en juillet-août, puis augmentation jusqu'à la fin de l'année
-> d'où une saisonnalité de 12 mois - De plus, on constate que les courbes annuelles se superposent : les niveaux de consommation sont équivalents chaque année
-> d'où une tendance stable dans le temps (on remarque l'année particulière de 2020 du fait de la crise sanitaire et qui se manifeste par une consommation moindre d'électricité)
In [10]:
with plt.style.context("seaborn-ticks"):
fig, axes = plt.subplots(1, 2, figsize=(15,5))
config_moy={"marker":"o",
"markerfacecolor":"white",
"markeredgecolor":"black",
"markeredgewidth":"0.5",
"markersize":"6"}
config_moy_little={"marker":"o",
"markerfacecolor":"white",
"markeredgecolor":"black",
"markeredgewidth":"0.5",
"markersize":"4"}
q_colors = sns.color_palette("Set1", 12, .7)
_=sns.boxplot(x='year', y='conso', data=df_fr[df_fr['year'] != 2021], showmeans=True, meanprops=config_moy, palette=(q_colors[i] for i,j in enumerate(df_fr['year'][:-1])), ax=axes[0]);
qualitative_colors = sns.color_palette("Set3", 18);
_=sns.boxplot(x='month', y='conso', data=df_fr.loc[~df_fr.year.isin([2012, 2020]), :], palette=(qualitative_colors[i] for i,j in enumerate(df_fr['month'])), meanprops=config_moy_little, showmeans=True);
_=axes[0].set_title('BoxPlot par année\n(la Tendance)', color="darkred", fontsize=14);
_=axes[0].set_ylabel("Consommation brute (en GWh)");
_=axes[0].set_xlabel("");
_=axes[1].set_title('BoxPlot par mois\n(La saisonnalité)', color="darkred", fontsize=14);
_=axes[1].set_ylabel("Consommation brute (en GWh)");
_=axes[1].set_xlabel("");
plt.tight_layout(w_pad=3)
plt.savefig(path_graph+"03-boxplot_annee_mois_conso_brute.png", dpi=200)
plt.show();
Boxplot par année
- cf l'évolution de la consommation mensuelle annualisée
- Confirmation du constat précédent : l'évolution de la consommation d'électricité est comparable d'une année sur l'autre : les niveaux sont relativement constants.
- On constate que c'est en 2017 que le mois avec la consommation la plus importante sur la période étudiée a eu lieu.
- A l'inverse, c'est 2020 qui connait le mois durant lequel la consommation a été la plus faible (et les moyennes mensuelles les plus faibles) : cf. l'effet covid
Boxplot par mois
- on peut classer les mois en 2 catégories :
- de novembre à mars : niveaux et dipersion élevés de la consommation
- d'avril à octobre : faibles niveaux et dispersion de la consommation
=> rend compte de l'importance des faibles températures sur le niveau et la dispersion de la consommation d'électricité
- on peut classer les mois en 2 catégories :
Intégration des données météo
Pour corriger les données de l'effet température, on va intégrer les Degrés Jour Unifié :- Le degré jour est une valeur représentative de l’écart entre la température d’une journée donnée et une température de référence, ici 18 °C
- On utilisera le DJU de chauffe, c'est-à-dire celui correspondant à la situation où la température moyenne de la journée est inférieure à la température de référence
--> donc si la température extérieure, pour une journée donnée, est de 9 °C, alors dju = 18 - 9 = 9 (si température supérieure à 18°, DJU de chauffe = 0)
---> puis les données sont agrégées par mois
Remarques
- Etant donné le nombre élevé de stations différentes (1 ou 2 par département), nous avons décidé de nous limiter à une seule station, celle de Paris
- Concernant la méthode de calcul, nous avons retenu la méthode simplifiée "Météo"
In [11]:
# import des dju chauffe
dju = pd.read_csv("data/inputs/dju_paris_meteo_chauffe.csv", index_col=[0], sep=";", decimal=",")
dju
Out[11]:
In [12]:
# nos données de consommation commencent le 1/1/2012 -> on ne garde que les dju depuis cette date
dju=dju[dju.index>2011]
# suppression de la colonne "total"
dju=dju.drop(columns="Total")
# on ordonnes les données par année (ordre croissant)
dju=dju.sort_index()
dju
Out[12]:
On va intégrer les dju dans le df de la consommation
On va au préalable limiter les données de df_fr aux périodes pour lesquelles on dispose des valeurs de dju :
--> jusqu'en mai 2020 inclus
In [13]:
# limitation des données de df_fr
df_fr = df_fr[df_fr.index < '2020-06-01']
# initialisation liste qui contiendra les valeurs de dju
li=[]
# pour chaque index
for i in dju.index:
# on ajoute les valeurs des colonnes
li.extend([dju.loc[i,col] for col in dju.columns])
# on ne garde que les 5 premiers mois de l'années 2020 -> on ne garde donc pas les 7 dernières valeurs de la liste
li=li[:-7]
# on verifie que les tailles de df_fr et de li sont identiques
if len(li) == df_fr.shape[0]:
df_fr['dju']=li
df_fr
Out[13]:
In [14]:
df_fr.dtypes
Out[14]:
In [15]:
plt.rc('grid', alpha=0.3, linestyle='--')
plt.rc('figure', figsize=(10, 6), dpi=100, titlesize=14)
plt.rc('axes', titlecolor="darkred", titlesize=14, labelcolor="gray", labelsize=12)
plt.rc('xtick', labelsize=11)
plt.rc('ytick', labelsize=11)
In [16]:
fig=plt.figure(figsize=(12,4))
plt.rc('grid', alpha=0.3, linestyle='--')
ax = fig.add_subplot(211)
ax.plot(df_fr['conso'], label="conso", color="teal")
plt.ylabel("Consommation (en GWh)", color="gray", fontsize=10)
plt.legend(fontsize=9)
ax = fig.add_subplot(212)
ax.plot(df_fr['dju'], label="DJU", color="red")
plt.ylabel("DJU (en °C)", fontsize=10)
plt.suptitle('Comparaison graphique "consommation électrique" et graphique "DJU"', y=0.93, color="darkred")
plt.legend(fontsize=9)
plt.tight_layout()
plt.savefig(path_graph+"04-serie_temporelle_conso_France_brute_et_dju_2_graphs.png", dpi=200)
plt.show();
In [17]:
fig, ax1 = plt.subplots(figsize=(15,5))
color = 'teal'
ax1.plot(df_fr["conso"], label="Consommation (en GWh)", color=color)
plt.ylabel("Conso (en GWh)", color=color, labelpad=10)
plt.legend(loc="upper left")
ax2 = plt.gca().twinx()
color2='red'
ax2.plot(df_fr["dju"], label="DJU (en d °C)", color=color2)
plt.ylabel("DJU (en °C)", color=color2, labelpad=10)
plt.legend(loc="upper right")
plt.suptitle('Superposition graphique "consommation électrique" et graphique "DJU"', y=0.95)
plt.tight_layout()
plt.savefig(path_graph+"05-serie_temporelle_conso_France_brute_et_dju_1_graph.png", dpi=200)
plt.show();
In [18]:
#export des data
df_fr.to_csv("data/exports/df_fr.csv", header=True, index=True)
df_fr_year.to_csv("data/exports/df_fr_year.csv", header=True, index=True)
- Remarque :
- on constate la présence, chaque année, de pics de consommation en juillet, qui ne semblent pas être pris en compte par les DJU.
-> or nous n'avons intégré que les DJU de chauffe : on peut donc penser que ces pics correspondent à l'utilisation des climatisations.
(mais nous nous limiterons aux seuls DJU de chauffe pour notre analyse)
- on constate la présence, chaque année, de pics de consommation en juillet, qui ne semblent pas être pris en compte par les DJU.
In [19]:
from fonction_regression_lineaire import *
# modèle de régression linéaire
reg = ols('conso~dju', data=df_fr).fit()
In [20]:
# verification des hypothèses d'application du modèle linéaire
res_hyp = hypotheses_regression_lineaire('conso', ['dju'], df_fr, 'conso_dju_chauffe')
In [21]:
# analyse du modèle
res_analyse = analyse_reg(reg, droite=True, dju="dju", df=df_fr, version_name="conso_dju_chauffe")
- Correction des données de consommation mensuelles de l'effet de température
-> on va donc supprimer de la série temporelle la variabilité de la consommation causée par les faibles températures (et donc par l'utilisation de chauffage)
In [22]:
# on récupère les coefficients
coef_dju = res_analyse['coef']['dju']
# serie conso corrigée de effet de chauffe = serie conso brute - saisonnalité DJU de chauffe
df_fr['conso_corr_dju'] = df_fr['conso'] - df_fr['dju']*coef_dju
In [23]:
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso'], color="teal", label="Consommation brute")
plt.plot(df_fr['conso_corr_dju'], color="red", label="Consommation corrigée")
plt.title("Consommation d'électricité brute et corrigée de l'effet des températures en France (en GWh)")
plt.ylabel("Consommation (en GWh)")
plt.legend(fontsize=8)
plt.tight_layout()
plt.savefig(path_graph+"06-serie_temporelle_conso_France_brute_et_corrigee_effet_temperature.png", dpi=200)
plt.show();
In [24]:
df_fr.sample(1)
Out[24]:
In [25]:
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso_corr_dju'], color="red", label="Consommation corrigée", zorder=3)
for i in range(9):
datte = df_fr.index[i*12]
plt.axvline(x=datte, color='lightgray', linestyle='--', zorder=2)
plt.ylabel("Consommation (en GWh)")
plt.title("Consommation d'électricité corrigée de l'effet des températures en France (en GWh)")
plt.legend(loc="best", fontsize=10)
plt.tight_layout()
plt.savefig(path_graph+"07-serie_temporelle_conso_France_corrigee_effet_temperature.png", dpi=200)
plt.show();
Une série temporelle est caractérisée par les éléments suivants :
une tendance ($T_t$) :
- cf. l'évolution à long terme de la série
- dans notre cas, la tendance est constante (cf. graph 01 "Consommation totale d'électricité" et graph 02 "Saisonnalité de la consommation d'électricité")
- plusieurs méthodes possibles pour modéliser la tendance :
- moyennes mobiles,
- différenciation ($∆y_t = y_t − y_{t−1}$)
une saisonnalité ($S_t$) :
- cf. l'existence d'un modèle distinct, doté d'un pattern propre et répété à intervalles réguliers en raison de facteurs saisonniers (l'heure, le jour, le mois, le trimestre...)
- toutes les séries n'ont pas nécessairement une saisonnalité
- dans notre cas, nous avons précédemment traité de la partie liée aux températures ; nous traiterons de la saisonnalité totale dans la prochaine partie en appliquant une moyenne mobile
- Les traitements s’appliquant à la tendance s’appliquent également à la saisonnalité (moyennes mobiles, différenciation...)
- un cycle ($C_t$)
- cf. l'existence d'un phénomène répétitif, de fréquence inconnue ou changeante (non périodiques)
- toutes les séries n'ont pas nécessairement un cycle
- notre série n'a pas a priori de cycle particulier
- un résidu ($\epsilon_t$)
- cf. l'erreur (différence entre valeur observée et prédite)
-> elle correspond à la partie aléatoire de la série
-> on cherchera, si ce n'est déjà le cas, à la rendre stationaire - dans notre cas, nous obtiendrons ces résidus par différence :
- pour l'instant, nous avons la conso corrigée (CR) de l'effet température ;
- dans la prochaine partie, nous appliquerons une moyenne mobile de période 12 sur CR : on obtiendra ainsi une primo estimation de la tendance
-> il restera alors une double composante "saisonnalité + erreur"
-> on appliquera alors une régression linéaire sur cette double composante afin d'obtenir la version corrigées des variations saisonnières, c'est-à-dire le résidus
-> on décomposera alors la série temporelle en ses différentes composantes grâce à la fonctionseasonal_decomposede 'statsmodels :- on vérifiera les valeurs obtenues par cette fonction avec celles calculées :
- valeurs de la tendance (obtenues via les MM)
- valeurs du résidus et de la saisonnalité (obtenues par régression linéaire)
- on vérifiera les valeurs obtenues par cette fonction avec celles calculées :
- cf. l'erreur (différence entre valeur observée et prédite)
Séries temporelles additives et multiplicatives :
En fonction de la nature de la tendance et de la saisonnalité, la série temporelle peut être modélisée :
- soit comme une série additive : chaque observation peut s'exprimer comme la somme des 3 composantes :
- Valeur = Tendance + Saisonnalité + Erreur ==> $X_t = T_t + S_t + {\large\epsilon}_t$
- Ainsi, une série est additive si l'importance de l'effet saisonnier sur les données ne dépend pas de la taille des données
-> l'importance de l'effet saisonnier ne varie pas, que la série augmente ou diminue : il est constant dans le temps
- soit comme une série multiplicative : chaque observation peut s'exprimer comme le produit des 3 composantes :
- Valeur = Tendance Saisonnalité Erreur ==> $X_t = T_t \times S_t \times {\large\epsilon}_t$
- Ainsi, une série est multiplicative si l'importance de l'effet saisonnier sur les données dépend de la taille des données
-> l'importance de l'effet saisonnier augmente au fur et à mesure que les valeurs de la série augmentent, et diminue lorsque les valeurs de la série diminuent.
- soit comme une série additive : chaque observation peut s'exprimer comme la somme des 3 composantes :
Concernant notre série temporelle, le modèle additif semble le mieux convenir au vu des graphiques précédents
3.1. Calcul des différentes composantes¶
3.1.1. Recherche de la tendance de la série de consommation corrigée¶
On va maintenant enlever la saisonnalité de la consommation corrigée de l'effet température
Pour rappel, la saisonnalité est ici de 12 mois
Méthode utilisée :
- On va lisser la série temporelle de la consommation corrigée de l'effet température en appliquant une moyenne mobile simple sur 12 mois : M12
-> cela va permettre d'éliminer les fluctuations les moins significatives
- On va lisser la série temporelle de la consommation corrigée de l'effet température en appliquant une moyenne mobile simple sur 12 mois : M12
3 objectifs dans la recherche de la MM optimale :
-> on cherche la MM :- qui laisse la tendance invariante
- qui absorbe la saisonnalité
- qui réduit la variance du résidu
In [26]:
# calcul des moyennes mobiles
def moy_mob(serie, p):
nb_elemt = len(serie)
k = p // 2
mb = []
for i in range(k):
mb.append(np.nan)
if p%2 != 0:
for i in range(nb_elemt-p+1):
s=serie[i:i+p]
mb.append(sum(s)/p)
else:
for i in range(nb_elemt-p):
s=serie[i:i+p+1]
mb.append((0.5*s[0]+0.5*s[-1]+sum(s[1:-1]))/p)
for i in range(k):
mb.append(np.nan)
mbr = [round(x,2) for x in mb]
return mbr
On calcule M12 sur la série corrigée
conso_corr_dju
-> on obtient ainsi la tendance de la série corrigéeOn ajoute cette tendance dans notre df, dans la colonne
conso_corr_lissee
In [27]:
x=df_fr['conso_corr_dju'].tolist()
y=moy_mob(x,12)
df_fr['conso_corr_lissee']=y
df_fr.head(10)
Out[27]:
In [28]:
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso_corr_dju'], color="red", label="Consommation corrigée", zorder=3)
plt.plot(df_fr['conso_corr_lissee'], color="green", label="Consommation corrigée lissée", zorder=4)
for i in range(9):
datte = df_fr.index[i*12]
plt.axvline(x=datte, color='lightgray', linestyle='--', zorder=2)
plt.ylabel("Consommation (en GWh)")
plt.title("Consommation d'électricité corrigée et corrigée lissée via M12")
plt.legend(loc="best", fontsize=9)
plt.tight_layout()
plt.savefig(path_graph+"08-serie_temporelle_conso_France_corrigee_et_corrigee_et_lissee_M12.png", dpi=200)
plt.show();
In [29]:
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso'], color="lightseagreen", label="Consommation brute")
plt.plot(df_fr['conso_corr_dju'], color="red", label="Consommation corrigée")
plt.plot(df_fr['conso_corr_lissee'], color="green", label="Consommation corrigée lissée", zorder=4)
plt.title("Consommation d'électricité brute, corrigée, et corrigée lissée (en GWh)")
plt.ylabel("Consommation (en GWh)")
plt.legend(loc="best", fontsize=8)
plt.tight_layout()
plt.savefig(path_graph+"09-serie_temporelle_conso_France_brute_et_corrigee_et_lissee.png", dpi=200)
plt.show();
In [30]:
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso_corr_lissee'], color="green", label="Consommation corrigée lissée", zorder=3)
plt.ylabel("Consommation (en GWh)")
plt.title("Consommation d'électricité corrigée et lissée via M12 (en GWh) : la tendance")
plt.legend(loc="best", fontsize=8)
plt.tight_layout()
plt.savefig(path_graph+"10-serie_temporelle_conso_France_corrigee_lissee_M12.png", dpi=200)
plt.show();
3.1.2. Désaisonnaliser la double composante "saisonnalité - résidus"¶
- La différence (car modèle additif) entre la série corrigée de l'effet température (
conso_corr_dju) et la tendance (conso_corr_lissee) est la double composante "saisonnalté et erreur"
-> on l'ajoute dans notre df dans une colonnedouble_compo
In [31]:
df_fr['double_compo'] = df_fr['conso_corr_dju'] - df_fr['conso_corr_lissee']
df_fr.tail(10)
Out[31]:
In [32]:
# copie du df
dft = df_fr.copy()
# limitations des données aux valeurs non NaN et années entières (multiple de 12)
dft = dft[(dft.index>"2012-06-01") & (dft.index<"2019.07.01")]
dft.shape
dft.head(2)
dft.tail(2)
Out[32]:
Out[32]:
Out[32]:
On va appliquer le modèle linéaire suivant : $Y_t = a + bt + \sum_{i=1}^{12} (c_i1)$
Mais problème de colinéarité : la somme des indicatrices est égale à 1 ; or "1" est déjà une variable
=> Solution : pour chaque $Y_t$ on enlève la valeur ajustée moyenne des indicatrices $\hat{c_i}$
In [33]:
# on isole la série à traiter
y=dft[['double_compo']]
# base tendancielle
t=np.arange(1,85,1)
# base saisonnière
# 12 indicatrices : Si vaudra "1" pour le mois i, quelque soit l'année considérée
for i in range(12):
su = np.repeat(0, repeats=12)
su[i] = 1
s = np.tile(su, 84 // len(su) + 1)[:84]
vars()['s' + str(i+1)] = s
In [34]:
from sklearn import linear_model
# on applique la réfression linéaire
reg = linear_model.LinearRegression(fit_intercept=False)
t=reg.fit(np.array([t, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12]).transpose(), y)
In [35]:
# modèle modifié pour pallier au problème de colinéarité
# valeur estiméee moyenne des 12 indicatrices
a = np.mean(reg.coef_[0][1:13])
# coef de t
b = reg.coef_[0][0]
# coef des 12 indicatrices corrigés
c = reg.coef_[0][1:13] - a
saisonnalite = (c[0]*s1+c[1]*s2+c[2]*s3+c[3]*s4+c[4]*s5+c[5]*s6+c[6]*s7+c[7]*s8+c[8]*s9+
c[9]*s10+c[10]*s11+c[11]*s12)
dft['saisonnalite'] = saisonnalite
dft['residus'] = dft['double_compo'] - dft['saisonnalite']
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 4))
dft['saisonnalite'].plot(ax=ax1, color="#ffa600")
dft['residus'].plot(ax=ax2, color="#ef5675")
y_lab=["Saisonnalité", "Résidus"]
for i,x in enumerate([ax1, ax2]):
x.set_xlim("2012-04-01", "2019-09-01")
x.set_ylabel(y_lab[i], color="gray", labelpad=15)
x.set_xlabel("", color="white", fontsize=0)
plt.suptitle("Décomposition de la double composante de la série corrigée de consommation d'électricité (en GWh)", y=0.95)
plt.tight_layout()
plt.savefig(path_graph+"11-decomposition_double_composante_saison_residus_conso_corrigee_dju.png", dpi=200)
plt.show();
3.2. Décomposition automatisée (seasonal_decompose)¶
3.2.1. Graphique des composantes de la série de consommation corrigée¶
In [36]:
pylab.rcParams['figure.figsize'] = (12, 6)
from statsmodels.tsa.seasonal import seasonal_decompose
def decompo_conso(s, model, num_graph):
result = seasonal_decompose(s, model=model, extrapolate_trend='freq')
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(15, 8))
zz = [ax1, ax2, ax3, ax4]
colors = ["#003f5c", "#7a5195", "#ef5675", "#ffa600"]
# Valeurs observées
result.observed.plot(ax=ax1, color=colors[0])
# Tendance
result.trend.plot(ax=ax2, color=colors[1])
# Saisonnalité
result.seasonal.plot(ax=ax3, color=colors[2])
for z in range(9):
year = 2012+z
ax3.axvline(x=str(year)+"-01-01", color='gray', linestyle='--', zorder=3)
# Résidus
result.resid.plot(ax=ax4, color=colors[3])
y_lab = ["Observations", "Tendance", "Saisonnalité", "Résidus"]
for i,z in enumerate(zz):
z.set_xlim('2011-11-01', '2020-07-01')
z.set_ylabel(y_lab[i], color="gray", labelpad=15)
z.set_xlabel("", color="white", fontsize=0)
plt.suptitle("Décomposition "+model+" de la consommation corrigée de l'effet température", y=0.97, color="darkred")
plt.tight_layout()
plt.savefig(path_graph+num_graph+"-decomposition_"+model+"_conso_corrigee_dju.png", dpi=200)
plt.show();
return result
In [37]:
result_add=decompo_conso(s=df_fr['conso_corr_dju'], model='additive', num_graph="12")
3.2.2. Récupération des valeurs des composantes¶
In [38]:
print(dir(result_add))
- Création d'un df avec les valeurs extraites de
seasonal_decompose
In [39]:
df_verif = pd.DataFrame(list(zip(result_add.trend, result_add.seasonal, result_add.resid)),
columns = ['sd_trend','sd_seasonal', 'sd_resid'],
index = df_fr.index)
df_verif
Out[39]:
3.3 Comparaison des résultats¶
- On va comparer les résultats obtenus !
- par le calcul
- par la fonction
seasonal_decomposede statsmodels
3.3.1 Comparaison graphique¶
In [40]:
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(15, 8))
zz = [ax1, ax2, ax3, ax4]
colors = ["#003f5c", "#7a5195", "#ef5675", "#ffa600"]
colors2 = ["#49c1f8", "#50fba3", "#302bda", "#085dbf"]
# Valeurs observées
result_add.observed.plot(ax=ax1, color=colors[0])
# Tendance
result_add.trend.plot(ax=ax2, color=colors[1], alpha=0.8)
dft['conso_corr_lissee'].plot(ax=ax2, color=colors2[1], linestyle="--", linewidth=1)
# Saisonnalité
result_add.seasonal.plot(ax=ax3, color=colors[2], alpha=0.8)
dft['saisonnalite'].plot(ax=ax3, color=colors2[2], linestyle="--", linewidth=1)
for z in range(9):
year = 2012+z
ax3.axvline(x=str(year)+"-01-01", color='gray', linestyle='--', zorder=3)
# Résidus
result_add.resid.plot(ax=ax4, color=colors[3], alpha=0.8, label="seasonal_decompose")
dft['residus'].plot(ax=ax4, color=colors2[3], linestyle="--", linewidth=1, label="composantes calculées")
y_lab = ["Observations", "Tendance", "Saisonnalité", "Résidus"]
for i,z in enumerate(zz):
z.set_xlim('2011-11-01', '2020-07-01')
z.set_ylabel(y_lab[i], color="gray", labelpad=15)
z.set_xlabel("", color="white", fontsize=0)
plt.legend(bbox_to_anchor=(1, 5), fontsize=8)
plt.suptitle("Comparaison de la décomposition calculée et de la décomposition de 'seasonal_decompose'", x=0.43, y=0.93, color="darkred")
plt.savefig(path_graph+"13-comparaison_decomposition_calculee_et_seasonal_decompose.png", dpi=200)
plt.show();
- On obtient la même tendance avec les 2 méthodes
- On constate quelques différences pour les composantes "saisonnalité" et "résidus" (différences complémentaires) mais elles sont marginales
On peut quantifier ces différences en comparant les valeurs de ces 2 composantes
3.3.1 Comparaison des valeurs¶
- On intégre les données calculées précédemment dans le df df_verif
In [41]:
df_fr_data = dft[['conso_corr_lissee', 'saisonnalite', 'residus']]
df_fr_data.columns = ['calc_trend', 'calc_seasonal', 'calc_resid']
df_verif = df_verif.merge(df_fr_data, left_index=True, right_index=True, how='inner')
df_verif
Out[41]:
In [42]:
comp = ['trend', 'seasonal', 'resid']
for c in comp:
df_verif['diff_'+c] = df_verif['sd_'+c] - df_verif['calc_'+c]
# affichage aléatoire de 10 valeurs
df_verif.sample(10)
Out[42]:
3.4. Consommation corrigée désaisonnalisée¶
- Effectuons la désaisonnalisation de la consommation corrigée de l'effet température
In [43]:
df_fr['conso_corr_deseason'] = df_fr['conso_corr_dju'] - result_add.seasonal
df_fr.sample(2)
Out[43]:
In [44]:
plt.figure(figsize=(12,4))
plt.plot(df_fr['conso_corr_dju'], color="#003f5c", linestyle='dashed', label="Consommation corrigée de l'effet température", zorder=2, alpha=0.5)
plt.plot(df_fr['conso_corr_deseason'], color="#ef5675", label="Consommation corrigée et désaisonnalisée", zorder=3)
plt.ylabel("Consommation (en GWh)")
plt.title("Consommation corrigée et consommation corrigée désaisonnalisée")
plt.legend(loc="best", fontsize=8)
plt.tight_layout()
plt.savefig(path_graph+"14-serie_temporelle_conso_corrigee_et_conso_corrigee_desaisonnalisee.png", dpi=200)
plt.show();
On constate graphiquement que la désaisonnalisation a permis de lisser la série :
-> la moyenne semble équivalente, avec une variance plus faible
In [45]:
df_fr[['conso_corr_dju', 'conso_corr_deseason']].describe().T.loc[:,['count', 'mean', 'std', 'min', 'max']].apply(lambda x: round(x,2))
Out[45]:
In [46]:
# df avec seulement la consommation corrigée et la consommation desaisonnalisée
df=df_fr.copy()
df=df_fr.drop(columns=['conso_corr_lissee', 'double_compo'])
df=df.rename(columns={'conso_corr_dju': 'conso_corr'})
df=df.rename(columns={'conso_corr_deseason': 'conso_desaison'})
df.head(2)
df.tail(12)
Out[46]:
Out[46]:
4.1. Méthode déterministe : méthode de Holt Winters¶
4.1.1. Le principe¶
Principe général :
La méthode de prévision de Holt Winters est une méthode empirique qui repose sur le lissage exponentiel :
-> Cette technique se base sur les moyennes pondérées d’observations antérieures, en accordant un poids d'autant plus important que l'observation est récente.La méthode de Hots Winters utilise un "triple lissage" : elle permet de gérer des séries comportant des niveaux, une tendance mais aussi une saisonnalité.
Principe d'application :
- Afin de pourvoir mesurer la performance des prévisions réalisées, nous allons diviser notre série initiale en deux jeux différents :
- un jeu d'entrainement x_train :
- données de janv-2012 à mai-2019
- on va utiliser ce jeu de données pour construire notre modèle afin de prédire la consommation d'électricité pour la période juin-19 à mai-20
- un jeu de test x_test
- données juin-19 à mai-20
- on comparera les prévisions et les valeurs réelles pour cette période
- un jeu d'entrainement x_train :
- Afin de pourvoir mesurer la performance des prévisions réalisées, nous allons diviser notre série initiale en deux jeux différents :
4.1.2. Exécution du modèle et prévisions obtenues (visualisation graphique)¶
In [47]:
# taille de x_train et x_test
df[df.index < '2019-01-01'].shape
df[df.index >= '2019-01-01'].shape
Out[47]:
Out[47]:
In [48]:
# preapration des données
df_x = df[['conso_corr']]
df_x_train = df[['conso_corr']][:84]
x_train=df_x_train['conso_corr']
In [49]:
from statsmodels.tsa.api import ExponentialSmoothing
import warnings
from statsmodels.tools.sm_exceptions import ConvergenceWarning
warnings.simplefilter('ignore', ConvergenceWarning)
warnings.simplefilter('ignore', FutureWarning)
# methode de Holt Winters
hw = ExponentialSmoothing(np.asarray(x_train), seasonal_periods=12, trend='add', seasonal='add').fit()
# prévisions
hw_pred = hw.forecast(17)
In [50]:
# ajout d'une colonne 'prevision' dans df
df['prevision_hw']=np.nan
df.loc[df.index > '2018-12-31', 'prevision_hw'] = hw_pred
df_prevision = df[['prevision_hw']]
- On affiche graphiquement les prévisions obtenues
In [51]:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 8))
zz = [ax1, ax2]
# x_train + x_test
df_x[:85].plot(ax=ax1, color='teal')
df_prevision[84:].plot(ax=ax1, color='red')
# x + x_test
df_x.plot(ax=ax2, color='teal')
df_prevision[84:].plot(ax=ax2, color='red', linestyle='dashed', alpha=0.8)
for z in zz:
z.set_xlim('2011-11-01', '2020-07-01')
z.set_ylabel("Consommation (en GWh)", color="gray", labelpad=15)
z.set_xlabel("", color="white", fontsize=0)
z.legend(['Consommation réelles', 'Consommation prédites'], loc="best", fontsize=9, frameon=True, edgecolor='black')
plt.suptitle("Prévision de la consommation corrigée à partir de 2019 - Méthode Holt Winters", y= 0.92, fontsize=14, color="darkred")
plt.savefig(path_graph+"15-Prevision_consommation_2019_Holt_Winters.png", dpi=200)
plt.show();
On fait un focus graphique sur la période prédite
Remarque :
Du fait du contexte particulier de la crise sanitaire débutée début 2020, on constate un "décrochage" des prévisions par rapport aux données réelles.
On présente donc 2 graphiques différents des prévisions :
- un graphique qui couvre la totalité de la période prédite (janv-19 à mai-20)
- un graphique qui se limite à l'année 2019
In [52]:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
zz = [ax1, ax2]
# x_train + x_test (19-mai20)
df_x[84:].plot(ax=ax1, color='teal')
df_prevision[84:].plot(ax=ax1, color='red', linestyle='dashed', alpha=0.8)
ax1.set_xlim('2018-12-01', '2020-06-01')
# x_train + x_test (2019)
df_x[84:-4].plot(ax=ax2, color='teal')
df_prevision[84:-4].plot(ax=ax2, color='red', linestyle='dashed', alpha=0.8)
ax2.set_xlim('2018-12-01', '2020-02-01')
for z in zz:
z.set_ylabel("Consommation (en GWh)", color="gray", labelpad=15)
z.set_xlabel("", color="white", fontsize=0)
z.legend(['Consommation réelles', 'Consommation prédites'], loc="best", fontsize=9, frameon=True, edgecolor='black', shadow=True)
plt.suptitle("Prévision de la consommation corrigée à partir de 2019 - Méthode Holt Winters", y= 0.96, fontsize=14, color="darkred")
plt.tight_layout()
plt.savefig(path_graph+"16-Prevision_consommation_2019_Holt_Winters_focus_periode_predite.png", dpi=200)
plt.show();
4.1.3. Analyse du modèle¶
- Calcul des résidus (et du pourcentage d'erreur pour chaque observation prédite)
In [53]:
dfr = df.copy()
dfr['residus'] = dfr['conso_corr'] - dfr['prevision_hw']
dfr['erreur_pct'] = (dfr['residus'] / dfr['conso_corr']).apply(lambda x: "{:.2%}".format(x))
for col in ['conso_corr', 'prevision_hw', 'residus']:
dfr[col] = dfr[col].apply(lambda x: round(x,2))
analys_res_hw = dfr[84:][['conso_corr', 'prevision_hw', 'residus', 'erreur_pct']]
analys_res_hw.columns=['y_test', 'pred_hw', 'res_hw', 'res_hw_pct']
analys_res_hw=analys_res_hw[:-5]
- Caclcul de la performance du modèle
In [54]:
print(dir(hw))
Calcul de :
- l'erreur quadratique moyenne (RMSE)
- l'erreur relative absolue moyenne (MAPE)
On récupère également les valeurs des critères d'informations AIC et BIC
In [55]:
from sklearn.metrics import mean_squared_error
def RMSE(y, predictions):
RMSE = mean_squared_error(y, predictions, squared=False)
return RMSE
def MAPE(y, predictions):
MAPE = np.mean(np.abs((y - predictions) / y)) * 100
return MAPE
res_pred_hw={}
res_pred_hw['res_pred_hw'] = analys_res_hw
y1 = df_x[84:]['conso_corr']
predictions1 = df_prevision[84:]['prevision_hw']
RMSE1=RMSE(y1, predictions1)
MAPE1=MAPE(y1, predictions1)
y2 = df_x[84:-4]['conso_corr']
predictions2 = df_prevision[84:-4]['prevision_hw']
RMSE2=RMSE(y2, predictions2)
MAPE2=MAPE(y2, predictions2)
print("#################################################")
print("#### Performances du modèle de Holt Winters ####")
print("#################################################")
print("#")
print(f"# AIC : {round(hw.aic,3)}")
res_pred_hw['aic']=round(hw.aic,3)
print(f"# BIC : {round(hw.bic,3)}")
res_pred_hw['bic']=round(hw.bic,3)
print("#")
print("######################################")
print("####### Période janv19-mai20 #######")
print("######################################")
print(f"# RMSE : {round(RMSE1, 2)}")
print(f"# MAPE : {round(MAPE1, 2)}")
print("#")
print("######################################")
print("####### Période janv19-dec19 #######")
print("######################################")
print(f"# RMSE : {round(RMSE2, 2)}")
res_pred_hw['rmse']=round(RMSE2, 2)
print(f"# MAPE : {round(MAPE2, 2)}")
res_pred_hw['mape']=round(MAPE2, 2)
4.2. Méthode stochastique : modèle SARIMA¶
4.2.1. Analyse de la stationnarité¶
Un processus temporel $X_t$ est dit stationnaire au sens faible (ou "de second ordre") s'il remplit les 3 conditions suivantes :
- l'espérance est constante au cours du temps
- la variance est constante au cours du temps
- l'autocorrélation entre 2 instants $t$ et $t - h$ ne dépend que de l'écart $h$ entre ces 2 instants (elle ne dépend donc pas du temps $t$)
Pour tester la stationnarité d'une série temporelle, on va utiliser différentes méthodes :
- méthodes graphiques : sortie ACF et PACF
- test de stationnarité : Test de racine unitaire tel que le test augmenté de Dickey-Fuller (test ADF)
- On commence par tester la stationnarité de la série avec un test ADF
Remarque concernant notre série temporelle initiale :
La crise sanitaire actuelle a eu pour conséquence une forte baisse de la consommation d'électricité corrigée.Ce choc exogène non prévisible est d'une ampleur non négligeable, et "perturbre" fortement les données.
L'objectif étant ici d'étudier les modèles ARIMA et SAMIRA, le choix a été fait de supprimer les données de 2020 du dataset initial :
-> ainsi, nos données s'arrêtent désormais au 31/12/2019.
In [56]:
df = df[df.index < '2020-01-01']
df
Out[56]:
In [57]:
y=df[['conso_corr']]
from fonctions_perso import test_stationnarite_adf
res=test_stationnarite_adf(y, "y", details=False)
Le test nous confirme que la série de consommation corrigée de l'effet température n'est pas stationnaire
Nous allons déstationnariser la série grâce à une ou plusieurs différenciations.
Dans un modèle SARIMA, 2 types de différenciation sont possibles :
une différenciation "en tendance", c'est-à-dire d'ordre 1
-> Concernant notre série temporelle, nous avons vu précédemment que la tendance est stable au cours du temps (une consommation corrigée mensuelle de l'ordre de 31500 GWh environ)
---> nous n'avons donc pas besoin de réaliser de différenciation en tendanceune différenciation "en saisonnalité", c'est-à-dire d'ordre équivalent au nombre de périodes d'une saison
-> Concernant notre série temporelle, nous avons déterminé une saisonnalité de 12 mois
---> nous allons donc a priori réaliser une ou plusieurs différenciations d'ordre 12
Analysons maintenant le graphique de la sortie ACF pour vérifier les autocorrélations
In [58]:
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plot_acf(y, zero=True, lags=40)
plt.xticks(np.arange(0, 42, step=3))
plt.show();
Nous allons commencer par la partie saisonnière, puis voir ensuite s'il est nécessaire ou non de différencier en tendance.
- Partie saisonnière:
La lecture de la sortie ACF indique des pics toutes les 6 périodes.
En regardant plus précisément la sortie ACF, on constate qu'il ne s'agit pas d'une saisonnalité de 6 mais bien d'une saisonnalité de 12 :
-> en effet, on constate le même schéma qui se répète :
- pour chaque période de 12, présence d'un pic à lag=6 et d'un pic plus important à lag=12 ;
- on constate également une décroissance lente de ces autocorélations simples de saisonnalité de 12 :
- pic 1 : lag = 6 + 12$k$
- pic 2 : lag = 12 + 12$k$
--> on va donc effectuer une différenciation d'ordre 12: $\nabla^{12} = (I - B^{12})$
On regarde la sortie ACF correspondante
In [59]:
ydiff_12 = y - y.shift(12)
plot_acf(ydiff_12, zero=True, lags=40, missing='drop')
plt.xticks(np.arange(0, 42, step=3))
plt.show();
La sortie ACF semble désormais correspondre à un autocorrélogramme simple (on ne constate plus de pics liés à une fréquence saisonnière).
Regardons la sortie PACF correspondante
In [60]:
plot_pacf(ydiff_12.dropna(), zero=True, lags=30)
plt.xticks(np.arange(0, 32, step=3))
plt.show();
Ainsi, cette différenciation d'ordre 12 semble suffisante pour rendre le processus startionnaire.
Vérifions la stationnarité avec un test de racine unitaire
In [61]:
res_y12 = test_stationnarite_adf(ydiff_12, 'ydiff_12', details=False)
La série ainsi différenciée parait donc stationnaire au seuil de 5% d'erreur.
Le modèle SARIMA sera donc de la forme :
$$\text SARIMA\; (p,\, 0,\, q)\,(P,\, 1,\, Q)\,_{12}$$
En effet :
- nous avons effectué une différenciation saisonnière d'ordre 12 1 fois : donc D = 1
- nous n'avons pas eu besoin de différencier la tendance pour stationnariser la série : donc d = 0
4.2.1. Détermination des paramètres du modèle SARIMA¶
Il faut maintenant déterminer les valeurs des paramètres p, q, P et Q :
-> on s'aide pour cela des graphiques des sorties ACF et PACF, à savoir :
- p est l'ordre de la partie autorégressive du modèle (ordre AR)
- P est l'ordre de la partie autorégressive saisonnière du modèle (ordre AR)
-> on utilise la sortie PACF pour déterminer le nombre optimal de termes à utiliser dans le modèle AR.
--> le nombre de termes détermine l’ordre du modèle. - q est l'ordre de la partie moyenne mobile du modèle (ordre MA)
- Q est l'ordre de la partie moyenne mobile saisonnière du modèle (ordre MA)
-> on utilise la sortie ACF pour déterminer le nombre optimal de termes à utiliser dans le modèle MA.
--> le nombre de termes détermine l’ordre du modèle.
Remarque : Pour les 2 graphiques ACF et PACF
-> seules les lignes verticales qui dépassent les lignes horizontales sont considérées comme significatives
In [62]:
# rappel des graphiques des sorties ACF et PACF
from fonctions_perso import graph_acf_pacf
graph_acf_pacf(ydiff_12, "Différenciation saisonnière d'ordre 12 effectuée 1 fois : ydiff_12", lags=30, zero=False, name_plt="18-ACF_PACF_ydiff12" )
In [63]:
ydiff_12_2 = ydiff_12 - ydiff_12.shift(12)
graph_acf_pacf(ydiff_12_2, "Différenciation saisonnière d'ordre 12 effectuée 1 fois : ydiff_12", lags=20, zero=False, name_plt="19-ACF_PACF_ydiff_12")
Essayons de déterminer les paramètres p,q, P et Q
- les autocorrélations simples et partielles au premier retard ne sont pas significatives :
--> on n'ajoute pas de terme AR ou MA à la partie tendance - l'autocorrélation simple au retard 12 est négative (situation normale quand une différenciation saisonnière a été effectuée) :
--> on ajoute un terme SMA au modèle (un seul car au retard 24 l'autocorrélation simple n'est pas significative) - les autocorrélations partielles aux retards 12 et 24 sont significatives :
--> on ajoute 2 termes SAR au modèle
Ainsi, un modèle possible serait le suivant :
- SARIMA: (0, 0, 0) x (2, 1, 1, 12)
La lecture des sorties ACF et PACF n'est pas évidente :
on va tester différentes valeurs pour les paramètres p,q, P et Q
---> on se limite, pour chaque paramètre, aux valeurs 0, 1 et 2 (des valeurs supérieures risquant d'entraîner un ajustement excessif des données et/ou des problèmes d'estimations)pour l'ordre des différenciations, nous gardons les valeurs qui nous ont permis de stationnariser la série, à savoir d=0 et D=1
--> déterminons toutes les combinaisons possibles
Déterminations de toutes les combinaisons envisageables¶
In [64]:
from fonctions_perso import combi_possibles
combi_sarima=combi_possibles(d_start=0, d_end=0, sd_start=1, sd_end=1, m=12, details=True)
Détermination du modèle optimal¶
Méthode utilisée :
a) Pour chacune des combinaisons déterminées précédemment :
1. on estime le modèle avec `sm.tsa.statespace.SARIMAX`
2. on vérifie la significativité des paramètres du modèle :
Cf. test de Student sur chaque paramètre :
-> si p_value < 0.05, alors le paramètre est significatif
3. si tous les paramètres ne sont pas significatifs, on passe au modèle suivant ;
si tous les paramètres sont significatifs :
* Test de blancheur du résidu (résidu assimilable à un bruit blanc ou pas) : test de Ljung-Box
-> si p_value > 0.05 (pour chaque retard indiqué), alors le résidu est assimilable à un bruit blanc
* Test de normalité du résidu : test de Shapiro
-> si p_value > 0.05 alors la distribution du résidu suit une loi normale
4. si tous les paramètres sont significatifs, et si le résidu est assimilable à un bruit blanc et suit une distribution normale, alors le modèle est valide :
-> le modèle fait partie des modèles potentiels -> on récupère les valeurs de ses critères d'information AIC et BIC
2. on vérifie la significativité des paramètres du modèle :
Cf. test de Student sur chaque paramètre :
-> si p_value < 0.05, alors le paramètre est significatif
3. si tous les paramètres ne sont pas significatifs, on passe au modèle suivant ;
si tous les paramètres sont significatifs :
* Test de blancheur du résidu (résidu assimilable à un bruit blanc ou pas) : test de Ljung-Box
-> si p_value > 0.05 (pour chaque retard indiqué), alors le résidu est assimilable à un bruit blanc
* Test de normalité du résidu : test de Shapiro
-> si p_value > 0.05 alors la distribution du résidu suit une loi normale
4. si tous les paramètres sont significatifs, et si le résidu est assimilable à un bruit blanc et suit une distribution normale, alors le modèle est valide :
-> le modèle fait partie des modèles potentiels -> on récupère les valeurs de ses critères d'information AIC et BIC
b) choix du modèle définitif parmi les les modèles potentiels
* on choisit le modèle qui possède le critère d'information le plus faible
Les modèles potentiels doivent satisfaire les conditions suivantes :
- La série est stationnaire
- Tous les paramètres SARIMA sont significatifs
- Le résidu :
- est assmilable à un bruit blanc
- et sa distribution suit une loi normale
S'il existe plusieurs modèles potentiels, on choisit le modèle dont le critère d'information (AIC ou BIC) est le plus faible
Nous allons procéder en 2 temps :
utilisation d'une fonction créée pour ce projet (
select_model_sarima, disponible dans le dossier "fonctions_perso") qui se base sur la méthode décrite précédemment ;
--> à noter que nous distinguerons 2 cas différents :- les modèles potentiels pour lesquels les 3 tests sont valides (test de significativité des paramètres, tests de blancheur et de normalité du résidu)
- les modèles possibles pour lesquels seuls les tests de significativité des paramètres et de blancheur des résidus sont valides ;
--> nous verrons pour ces derniers si l'hypothèse de normalité des résidus est acceptable via un Q-Q plot
utilisation de la fonction
auto-SARIMAdu packagepyramid-arimaqui utilise une méthode similaire (documentation)
In [65]:
# import des fonctions
from fonctions_perso import select_model_sarima
from fonctions_perso import affichage_graph_diagnostics
In [66]:
y = df['conso_corr']
pdq = combi_sarima['pdq']
seasonal_pdq = combi_sarima['s_pdq']
select_best_mod = select_model_sarima(y, pdq, seasonal_pdq, version_name="test_1_modele_optimal", periode_saison=12, perf="aic", details=False)
=> Aucun modèle potentiel ou possible;
Testons les mêmes paramètres avec la fonction auto-SARIMA du package pyramid-arima
In [67]:
import pmdarima as pm
y = df[['conso_corr']]
smodel = pm.auto_arima(y, start_p=0, start_q=0,
test='adf',
max_p=2, max_q=2, m=12,
start_P=0, seasonal=True,
d=0, D=1, trace=True,
error_action='trace',
suppress_warnings=True,
stepwise=True)
smodel.summary()
Out[67]:
Le modèle retourné n'est pas valide
--> la p_value d'un coefficient est supérieur à 0,05 ==> ce coefficient n'est donc pas significatif
Faisons une différenciation supplémentaire, à savoir une différenciation de la tendance, soit d=1, et testons la stationnarité
In [68]:
ydiff_12_1 = ydiff_12 - ydiff_12.shift(1)
res_y12_1 = test_stationnarite_adf(ydiff_12_1, 'ydiff_12_1', details=False)
La série est stationnaire avec cette nouvelle différenciation
=> regardons les sorties ACF et PACF correspondant à cette double différenciation
In [69]:
graph_acf_pacf(ydiff_12_1, "Double différenciation d=1 et D=12 : ydiff_12_1", lags=26, zero=False, name_plt="20-ACF_PACF_ydiff_12_1" )
La lecture des sorties ACF et PACF n'est pas évidente
=> on va maintenant tester différentes valeurs pour les paramètres p,q, P et Q (de 0 à 2), en fixant les paramètres de différenciation d=1 et D=12
In [70]:
combi2_sarima=combi_possibles(d_start=1, d_end=1, sd_start=1, sd_end=1, m=12, details=True)
In [71]:
y = df['conso_corr']
pdq = combi2_sarima['pdq']
seasonal_pdq = combi2_sarima['s_pdq']
select_best_mod2 = select_model_sarima(y, pdq, seasonal_pdq, version_name="test_2_modele_optimal", periode_saison=12, perf="aic", details=False)
Ainsi, le modèle optimal, c'est-à-dire celui qui possède les critères d'information AIC le plus faible, est le suivant :
SARIMA(0, 1, 2)x(1, 1, 2, 12)
On constate graphiquement que le Q-Q plot semble rendre compte d'une distribution normale des résidus
In [72]:
affichage_graph_diagnostics(select_best_mod2, 'mod_1', size_img='medium')
On va vérifier ce résultat en le comparant à celui obtenu avec la fonction auto-SARIMA du package pyramid-arima
In [73]:
y = df[['conso_corr']]
smodel = pm.auto_arima(y, start_p=0, start_q=0,
test='adf',
max_p=2, max_q=2, m=12,
start_P=0, seasonal=True,
d=1, D=1, trace=True,
error_action='trace',
suppress_warnings=True,
stepwise=True)
smodel.summary()
Out[73]:
Le modèle renvoyé par auto_arima est différent;
Nous allons donc tester ces 2 modèles :
- model1 (fonction perso) : SARIMA(0,1,2)x(1,1,2,12) AIC=1447.61
- model2 (fonction auto-arima) : SARIMA(2,1,0)x(0,1,2,12) AIC=1453.16
Vérification de la validité des modèles sur la série tronquée¶
On répartit la série temporelle de la consommation corrigée en 2 :
- y_train (base d'apprentissage)
--> du 01-01-2012 au 31-12-2018 (7 années) - y_test (base de test) :
--> du 01-01-2019 au 31-12-2019 (1 année)
- y_train (base d'apprentissage)
On vérifie que les 2 modèles sont toujours valides
--> on vérifie les éléments suivants :- stationnarité
- significativité des coefficients
- blancheur du résidu
- normalité du résidu
Si les modèles sont toujours valides, alors on peut réaliser des prévisions :
- on estime le modèle à partir de
y_train - on realise les prévions
y_prevpour une année supplémentaire (cf la période couverte par y_test) - on compare les valeurs prédites
y_prevet les valeurs réellesy_test - on mesure la performance des prévisions grâce au calcul des indicateurs RMSE et MAPE
- on estime le modèle à partir de
Test de stationnarité¶
--> on réalise un test ADF sur la série tronquée
In [74]:
ydiff_12_1_t=ydiff_12_1.dropna()[:-12]
res_y12_1_t = test_stationnarite_adf(ydiff_12_1_t, 'ydiff_12_1', details=False)
Ainsi, la p_value du test de stationnarité est légèrement supérieure à 0.05, soit 0.057
--> même si la valeur est limite, on va considérer malgré tout la série tronquée comme stationnaire
Tests de validité des modèles sur la période y_train¶
In [75]:
from fonctions_perso import prediction_sarima
y=df['conso_corr']
In [76]:
# Tests de validité du modèle 1
res_model1 = prediction_sarima(y, order=(0,1,2), s_order=(1,1,2,12), model_name="model_1", m=12, nb_graph=21, details=False, graph=True, res_only=False)
=> le modèle 1 est bien valide pour la période couverte par y_train (2012-2018)
In [77]:
# Test de validité du modèle 2
res_model2 = prediction_sarima(y, order=(2,1,0), s_order=(0,1,2,12), model_name="model_2", m=12, nb_graph=22, details=True, graph=True, res_only=False)
Le modèle 2 n'est plus un modèle potentiel sur la période 2012-2018 :
- en effet, un coefficient n'est plus significatif et un autre est très proche du seuil de 5%
de plus, le QQ Plot laisse penser que la distribution des résidus n'est pas forcément gaussienne
==> ainsi, on ne retient que le modèle 1, soit SARIMA(0,1,2)x(1,1,2,12)
4.3. Comparaison des prévisions obtenues¶
Nous allons comparer les résultats des prévisions obtenues par les deux types de modèles :
- Modèle Holt Winters
- Modèle SARIMA
4.3.1. Comparaison des performances¶
- Modèle Holt Winters
In [78]:
print(*res_model1)
print(*res_pred_hw)
In [79]:
df_res = res_pred_hw['res_pred_hw']
df_res
Out[79]:
- Modèle SARIMA
In [80]:
df_res['pred_sarima'] = res_model1['previsions_xtest']
df_res['res_sarima'] = df_res['y_test'] - df_res['pred_sarima']
df_res['res_sarima_pct'] = (df_res['res_sarima'] / df_res['y_test']).apply(lambda x: "{:.2%}".format(x))
for col in ['y_test', 'pred_sarima', 'res_sarima']:
df_res[col] = df_res[col].apply(lambda x: round(x,2))
df_res
Out[80]:
- Performances des modèles
In [81]:
print("Modèle Holt Winters :")
print("#####################")
print(f"AIC = {res_pred_hw['aic']}")
print(f"BIC = {res_pred_hw['bic']}")
print(f"RMSE = {res_pred_hw['rmse']}")
print(f"MAPE = {res_pred_hw['mape']}")
print()
print("Modèle SARIMA :")
print("#####################")
print(f"AIC = {res_model1['aic']}")
print(f"BIC = {res_model1['bic']}")
print(f"RMSE = {round(res_model1['rmse'],2)}")
print(f"MAPE = {round(res_model1['mape'],2)}")
4.3.2. Comparaison Graphique¶
In [82]:
plt.figure(figsize=(10,4))
plt.plot(df_res['y_test'], color='blue', label='Consommation réelle', alpha=0.35)
plt.plot(df_res['pred_hw'], color='red', label='Consommation prédite Holt Winters')
plt.plot(df_res['pred_sarima'], color='green', label='Consommation prédite SARIMA')
plt.ylabel("Consommation en GWh")
plt.title("Prédictions de la consommation corrigée pour l'année 2019 par le modèle Holt Winters et SARIMA")
plt.legend(loc="best", fontsize=9, frameon=True, edgecolor='black')
plt.tight_layout()
plt.savefig(path_graph+"23-Prediction_Conso_Corrigee_2019_HW_et_SARIMA_1_graph.png", dpi=200)
plt.show();
In [83]:
df_prevision_final = df_prevision[:-5]
df_prevision_final['prevision_sarima'] = df_res['pred_sarima']
df_x_final = df_x[:-5]
df_x_prev = df_x.merge(df_prevision_final, left_index=True, right_index=True)
df_x_prev
Out[83]:
In [84]:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 8))
zz = [ax1, ax2]
# x_train + x_test
df_x_prev['conso_corr'][:85].plot(ax=ax1, color='teal', alpha=0.4)
df_x_prev['prevision_hw'][84:].plot(ax=ax1, color='red')
df_x_prev['prevision_sarima'][84:].plot(ax=ax1, color='green')
# x + x_test
df_x_prev['conso_corr'].plot(ax=ax2, color='teal', alpha=0.4)
df_x_prev['prevision_hw'][84:].plot(ax=ax2, color='red', linestyle='dashed', alpha=0.8)
df_x_prev['prevision_sarima'][84:].plot(ax=ax2, color='green', linestyle='dashed', alpha=0.8)
for z in zz:
z.set_xlim('2011-11-01', '2020-02-01')
z.set_ylabel("Consommation (en GWh)", color="gray", labelpad=15)
z.set_xlabel("", color="white", fontsize=0)
z.legend(['Consommation réelles', 'Consommation prédites hw', 'Consommation prédites SARIMA'], loc="best", fontsize=9, frameon=True, edgecolor='black')
plt.suptitle("Prévision de la consommation corrigée à partir de 2019 - Méthode Holt Winters et SARIMA", y= 0.92, fontsize=14, color="darkred")
plt.savefig(path_graph+"24-Prediction_Conso_Corrigee_2019_HW_et_SARIMA_2_graphs.png", dpi=200)
plt.show();
Conclusion :
L'analyse des critères d'erreurs (cf. la RMSE et la MAPE) nous permet de conclure que le modèle SARIMA est plus performant que le modèle de Holt Winters :
* 2,56% d'erreurs pour le modèle déterministe
* 2,28% d'erreur pour le modèle stochastique
Cette conclusion est confirmée par l'analyse graphique
Mais nous pouvons néammoins souligner que ce pourcentage d'erreur est relativement important pour notre secteur d'activité.