Etude prévisionnelle des ventes
Abdelghafour Ait bennassar, Salma Mounji
Table des matières
- 1 Introduction
- 2 Importation des outils nécessaires
- 3 Analyse exploratoire des données
- 4 Modélisation statistique univarié : ARIMA
- 5 Modélisation univarié par l’approche d’apprentissage automatique
- 6 Modélisation multivariée : Vector AutoRegressif Moving Average (VARMA)
Introduction¶
En raison de la rude concurrence qui existe aujourd'hui, les sociétés s'efforcent continuellement d'augmenter leurs bénéfices, de réduire leurs coûts et d'obtenir une flexibilité face aux changements. Une prévision précise des ventes est un outil de gestion précieux pour atteindre les objectifs mentionnés, car cela conduit à un service client amélioré en évitant les ruptures de stock, ainsi qu'à une réduction des coûts perdus causés par les excès de stock.
En d'autres termes, des prĂ©visions de ventes exactes sont utilisĂ©es pour saisir le compromis entre la satisfaction de la demande des clients et les coĂ»ts d'inventaire. En particulier,pour l'industrie pharmaceutique, des systèmes de prĂ©vision des ventes rĂ©ussis peuvent Ăªtre très bĂ©nĂ©fiques, en raison de la courte durĂ©e de conservation de nombreux produits pharmaceutiques et de l'importance de la qualitĂ© du produit qui est Ă©troitement liĂ©e Ă la santĂ© humaine.
Le notebook suivant représente l'étude de cas réalisée dans le cadre de notre projet de fin d'année basée sur les données historiques des ventes de l'une des principales sociétés de distribution pharmaceutique au MAROC (SOREMED).
Afin de construire et comparer des modèles de prévision des ventes selon les approches d’intelligence artificielle et économétrique, nous avons travaillé sur les ventes mensuelles d'un échantillon de 70 produits sur une durée de cinq ans.
Importation des outils nécessaires¶
from IPython.display import display
from math import ceil, sqrt
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import seaborn as sns
import statsmodels.api as sm
import joblib
import random
# colors generator
from itertools import cycle
color_pal = plt.rcParams['axes.prop_cycle'].by_key()['color']
color_cycle = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
# Preprocessing
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from statsmodels.tsa.seasonal import seasonal_decompose
# Clustering Algorithms
from tslearn.barycenters import dtw_barycenter_averaging
from sklearn.decomposition import PCA
from sklearn.metrics import pairwise_distances_argmin_min
from sklearn.cluster import KMeans
# ARIMA
from pmdarima import auto_arima
from pmdarima.utils import diff_inv
from pmdarima.arima.utils import ndiffs, nsdiffs
# Machine learning
from xgboost import XGBRegressor
from sklearn.linear_model import LinearRegression, LassoCV, RidgeCV
from sklearn.model_selection import TimeSeriesSplit
# VARMAX
from autots import AutoTS
# metrics
from sklearn.model_selection import cross_val_score
from sklearn.metrics import accuracy_score
from sklearn.metrics import mean_absolute_percentage_error, mean_absolute_error, mean_squared_error
# assert reproducibility
random.seed(7)
np.random.seed(7)
print("setup complete")
import warnings
warnings.filterwarnings('ignore')
import inspect
def get_default_args(func):
signature = inspect.signature(func)
return {
k: v.default
for k, v in signature.parameters.items()
if v.default is not inspect.Parameter.empty
}
from functools import wraps
class ForceRunException(Exception): pass
import joblib
def meta_decorator(joblib_Filename, compress=False):
def real_decorator(function):
@wraps(function)
def wrapper(*args, **kwargs):
force_run_arg = kwargs.get('force_run', None)
default_args = get_default_args(function)
default_force_run = default_args.get('force_run', False)
if force_run_arg is None:
force_run = default_force_run
else :
force_run = force_run_arg
assert force_run in {True, False}
filename = joblib_Filename
if compress:
filename += '.z'
try:
if force_run:
raise ForceRunException
My_Model = joblib.load(filename)
warnings.warn("function call and arguments are not evaluated !!")
print("data found and successfully imported,",
"No need to call the function")
except (ForceRunException, FileNotFoundError) as e:
if e.__class__ == ForceRunException:
print("Explicit execution requested,", end=' ')
else:
print("data not found,", end=' ')
print("Executing calculations...")
My_Model = function(*args, **kwargs)
print("calculations are done.", "saving data to " + filename)
joblib.dump(My_Model, filename)
print("the data is now saved")
return My_Model
return wrapper
return real_decorator
from time import sleep
@meta_decorator("joblib_test.joblib")
def test_save(number, times, *, force_run=False):
"""repeat number in list for number of times"""
sleep(3) # simulate waiting
to_save = [number for _ in range(times)]
return to_save
def pretty_df(data_f):
df = data_f.copy()
df["Series"] = df.index
cols = df.columns.tolist()
cols = cols[-1:] + cols[:-1]
df = df[cols]
return df.style.hide_index().set_properties(
**{'color': 'black !important',
'border': '1px black solid !important', "text-align": "left"}
).set_table_styles([{
'selector': 'th',
'props': [('border', '1px black solid !important'), ("text-align", "center")]
}]).set_properties(
subset=[cols[0]],
**{'font-weight': 'bold'}
)
def remove_axes(n_used_axes, all_axs):
if n_used_axes % 2:
all_axs = all_axs.flatten()
size = len(all_axs)
for useless_ax in all_axs.flatten()[-(size - n_used_axes):]:
useless_ax.axis('off')
sales = pd.read_csv('70prod_data.csv', index_col=0, parse_dates=True)
display(sales.head())
display(sales.describe().T)
Description de données¶
Initialement, on disposait des donnĂ©es sous forme d’un classeur Excel contenant cinq feuilles oĂ¹ chaque feuille contient les ventes mensuelles des produits dans un an, pour obtenir les donnĂ©es sous la forme requise pour l’étude, on a joint les donnĂ©es des cinq annĂ©es pour obtenir un jeu de donnĂ©es dans lequel les lignes reprĂ©sentent les dates (63 mois) et les colonnes reprĂ©sentent les produits (70 produits).
Le corrélogramme : visualisation de la matrice de corrélation¶
#ax = sns.heatmap(sales.corr())
cg = sns.clustermap(sales.corr(), figsize=(8, 8),
cbar_pos=(.1, .1, .03, .65), cmap="coolwarm")
cg.ax_row_dendrogram.set_visible(False)
cg.ax_col_dendrogram.set_visible(False)
plt.show()
On observe une forte correlation entre certains groupes de produits, ce qui suggère qu'on peut les classifier selon leurs comportements de ventes.
clustering : Regroupement des produits¶
Normalisation des variables¶
L'un des problèmes qu'on doit traiter avant de commencer le partitionnement (clustering) est l'échelle de la série. Sans normaliser les données, les séries qui se ressemblent seront vues si différentes les unes des autres et affecteront la précision du processus de clustering. Nous pouvons voir l'effet de la normalisation dans les images suivantes.
a = [[2], [7], [11], [14], [19], [23], [26]]
b = [[20000000], [40000000], [60000000], [80000000],
[100000000], [120000000], [140000000]]
fig, axs = plt.subplots(1, 3, figsize=(25, 5))
axs[0].plot(a)
axs[0].set_title("Series 1")
axs[1].plot(b)
axs[1].set_title("Series 2")
axs[2].plot(a)
axs[2].plot(b)
axs[2].set_title("Series 1 & 2")
plt.figure(figsize=(25, 5))
plt.plot(MinMaxScaler().fit_transform(a))
plt.plot(MinMaxScaler().fit_transform(b))
plt.title("Normalized Series 1 & Series 2")
plt.show()
mySeries = sales.copy()
scaler = MinMaxScaler()
mySeries = MinMaxScaler().fit_transform(mySeries).T
Réduction de la dimension : Analyse en composantes principales (ACP)¶
Afin de regrouper nos séries avec des k-moyennes, les indices temporels des séries chronologiques seront considérés comme des caractéristiques différentes et seront les dimensions des points de données(les séries).
Puisque la longueur des série temporelles en pratique est souvent importante, un autre problème auquel on doit faire face est la malédiction de la dimensionnalité,ce terme a été inventé pour la première fois par Richard E. Bellman lors de l'examen des problèmes de programmation dynamique. Cela signifie essentiellement que lorsque la dimensionnalité des données augmente, la distance entre les points de données augmente également. Ainsi, ce changement de mesure de la distance affecte gravement les algorithmes basés sur la distance.
Pour résoudre ce problème, on va effectuer une analyse en composantes principales ACP avant le partitionnement.
pca = PCA(n_components=2)
pca_res = pca.fit_transform(mySeries)
Désormais avec moins de dimensions qu'avant, nous pouvons voir comment nos séries se répartissent en 2 dimensions.
plt.figure(figsize=(25, 10))
plt.scatter(pca_res[:, 0], pca_res[:, 1], s=125)
plt.show()
Clustering : K-means¶
Le partitionnement en k-moyennes (ou k-means en anglais) est une méthode de partitionnement de données et un problème d'optimisation combinatoire. Étant donnés des points et un entier k, le problème est de diviser les points en k groupes, souvent appelés clusters, de façon à minimiser une certaine fonction. On considère la distance d'un point à la moyenne des points de son cluster (centroid) ; la fonction à minimiser est la somme des carrés de ces distances,la distance euclidienne est souvent utilisée.
Initialement on va éstimer le nombre de clusters requis par la racine du nombre des points (séries), on obtient 9 clusters.
cluster_count = ceil(sqrt(len(mySeries)))
cluster_count
et on va appliquer la méthode de Kmeans sur le résultat de l'ACP.
@meta_decorator("kmeans_object.joblib")
def create_fit_kmeans(pca_res, cluster_count):
kmeans = KMeans(n_clusters=cluster_count, max_iter=50000)
kmeans.fit(pca_res)
return kmeans
kmeans = create_fit_kmeans(pca_res, cluster_count)
labels = kmeans.predict(pca_res)
Après l'entrainement du modèle, on a prĂ©sentĂ© les rĂ©sultats. Pour chaque cluster, on a tracĂ© les sĂ©ries en gris, et afin de visualiser le mouvement et le comportement des sĂ©ries appartenant au mĂªme cluster, on a pris la moyenne du cluster et tracĂ© cette sĂ©rie moyenne en rouge.
labels
clusters = [[] for _ in range(cluster_count)]
for i, lab in enumerate(labels):
clusters[lab].append(i)
clusters
plot_count = ceil(sqrt(cluster_count))
fig, axs = plt.subplots(plot_count, plot_count, figsize=(20, 20))
axs = axs.flatten()
for i, cluster in enumerate(clusters):
pd.DataFrame(mySeries).T.iloc[:, cluster].plot(
ax=axs[i], c="gray", alpha=0.4, legend=False)
if len(cluster) > 0:
axs[i].plot(dtw_barycenter_averaging(
pd.DataFrame(mySeries).iloc[cluster]), c="red")
axs[i].set_title("Cluster "+str(i))
plt.show()
Nous pouvons voir la distribution des séries chronologiques en clusters dans le graphique suivant.
les séries semblent uniformement distribuées aux clusters.
clusters_c = list(map(len, clusters))
clusters_n = ["cluster_"+str(i) for i in range(cluster_count)]
plt.figure(figsize=(12, 4))
plt.title("Cluster Distribution for KMeans")
plt.bar(clusters_n, clusters_c)
plt.show()
Pour bénéficier de partitionnement qu'on a fait, il est nécessaire de lister pour chaque cluster, les séries qui lui appartiennent.
for cluster, cluster_n in zip(clusters, clusters_n):
# print(cluster_n, ':', [sales.columns[p] for p in cluster])
print(cluster_n+':', ', '.join(['P_'+str(p) for p in cluster]))
L’un des avantages que nous offre le regroupement effectuĂ© est l’exploitation de l’information disponible non seulement sur les ventes du produit qu’on souhaite prĂ©voir, mais Ă©galement l’information sur les produits ayant un comportement de ventes similaire (produits appartenant au mĂªme groupe).
On a selectionnĂ© le produit le plus proche du centroide de chaque cluster, autrement dit le plus proche aux autres produits au mĂªme groupe.
closest, _ = pairwise_distances_argmin_min(kmeans.cluster_centers_, pca_res)
closest
Pour visualiser ce résultat, on a tracé l'ensemble des points ainsi que les centres des groupes (étoiles), le point le plus proche de chaque centre est présenté en gris.
fig, ax = plt.subplots(figsize=(20, 10))
ax.scatter(pca_res[:, 0], pca_res[:, 1], c=labels, s=300)
sc = ax.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[
:, 1], marker='*', c=np.unique(labels), edgecolors='black', s=500, label=np.unique(labels))
legend1 = ax.legend(*sc.legend_elements(),
loc="upper left", title="Clusters", fontsize="x-large")
legend1.get_title().set_fontsize('15')
ax.add_artist(legend1)
ax.scatter(pca_res[closest, 0], pca_res[closest, 1], c='grey',
edgecolors=[sc.to_rgba(lab) for lab in np.unique(labels)], s=300)
for cluster_lab, cluster_x, cluster_y in zip(np.unique(labels), kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1]):
for pt_lab, pt_x, pt_y in zip(labels, pca_res[:, 0], pca_res[:, 1]):
if cluster_lab == pt_lab:
x, y = [cluster_x, pt_x], [cluster_y, pt_y]
ax.plot(x, y, c=sc.to_rgba(cluster_lab), alpha=0.5)
plt.show()
pour re vérifier ces résultats, on peut voir que les séries choisis sont centrées pour chaque groupe.
fig, axs = plt.subplots(plot_count, plot_count, figsize=(20, 20))
axs = axs.flatten()
for i, cluster in enumerate(clusters):
pd.DataFrame(mySeries).T.iloc[:, cluster].plot(
ax=axs[i], c="gray", alpha=0.4, legend=False)
if len(cluster) > 0:
axs[i].plot(dtw_barycenter_averaging(
pd.DataFrame(mySeries).iloc[cluster]), c="red")
pd.DataFrame(mySeries).T.iloc[:, [p for p in cluster if p in closest]].plot(
ax=axs[i], c="blue", legend=False)
axs[i].set_title("Cluster "+str(i))
plt.show()
Analyse de données¶
On a extrait les produits représentatifs des groupes et on a construit un nouveau jeu de données de séries temporelles sur lequel on va travailler par la suite.
sales_s = sales.iloc[:, sorted(closest)]
sales_s.head()
Statistiques descriptives :¶
Dans cette section, nous allons trouver quelques mesures de tendance et de dispersion comme la moyenne, médiane, la variance et l’écart type pour mieux percevoir les caractéristiques des médicaments.
display(sales_s.describe().T.round(1))
plt.figure(figsize=(9, 7))
sns.boxplot(data=sales_s)
plt.show()
D'après le tableau ci-dessus et le graphique, il est clair que ces médicaments ont des natures diverses et nous pouvons étendre cette remarque à tous les médicaments. Leurs moyennes et leur somme de ventes sont très différentes. Certains d'entre eux ont été vendus fortement comme le médicament 30, et d’autres sont moins vendus comme les médicaments 47 et 64. Certains médicaments ont une variance élevée comme les médicaments 30 et 22. Par conséquent, nous pouvons conclure que la majorité des médicaments ont des caractéristiques différentes en plus de la dissemblance dans le comportement de vente.
Visualisation de la tendance et saisonnalité¶
sales_s.plot(subplots=True, layout=(3, 3), figsize=(20, 15), marker='.')
plt.show()
On peut observer une saisonnalitĂ© dans certains produits notamment P_64 et p_47, avec une variĂ©tĂ© de tendances, certaines gardent la mĂªme monotonie au cours du temps tel que celles des P_64 et d’autre produits possèdent une tendance qui alterne entre haussière et baissière. L’existence de hausses soudaine et significative dans le produit P_22 pourra biaiser la prĂ©vision
Pour mieux visualiser la saisonnalité et tendance des ventes des médicaments, on a tracé pour chaque produit, une moyenne mobile 4-mois (pour détecter la saisonnalité) ainsi qu’une moyenne mobile 12-mois (pour capter la tendance)
fig, axs = plt.subplots(5, 2, figsize=(20, 20))
axs[-1, -1].axis('off')
axs = axs.flatten()
cols = [next(color_cycle), next(color_cycle), next(color_cycle)]
for i, item in enumerate(sales_s.columns):
sales_s[item].plot(title=item,
color=cols[0],
ax=axs[i],
marker='.', label="monthly")
sales_s[item].rolling(3).mean().plot(
color=cols[1],
ax=axs[i], label="quarter rolling mean", alpha=0.75)
sales_s[item].rolling(12).mean().plot(
color=cols[2],
ax=axs[i], linestyle='--', label="12months rolling mean")
axs[i].legend()
plt.show()
Les graphiques ci-dessous représentent la décomposition des séries temporelles
def plotseasonal(res, axes, name):
res.observed.plot(ax=axes[0], legend=False, title=name)
axes[0].set_ylabel('Observed')
res.trend.plot(ax=axes[1], legend=False)
axes[1].set_ylabel('Trend')
res.seasonal.plot(ax=axes[2], legend=False)
axes[2].set_ylabel('Seasonal')
res.resid.plot(ax=axes[3], legend=False, marker='o',
linestyle='', markersize=5)
axes[3].set_ylabel('Residual')
plotseasonal.counter += 1
def grid_plotseasonal(df, ncol):
nrow = ceil(len(df.columns)/ncol)*4
fig, axes = plt.subplots(ncols=ncol, nrows=nrow,
sharex=True, figsize=(15, 1.5*nrow))
plotseasonal.counter, k = 0, 0
for i, col in enumerate(df.columns):
res = seasonal_decompose(df[col].dropna(), period=12, filt=None)
plotseasonal(res, axes[k:k+4, i % ncol], col)
if (i % ncol) == (ncol-1):
k += 4
plt.tight_layout()
plt.show()
grid_plotseasonal(sales_s, ncol=3)
L'analyse des portions de chaque composant de la série chronologique est particulièrement utile pour déterminer l'absorption des résidus dans les données, sur la base des données décomposées. Le volume de cette absorption implique la prévisibilité de la série chronologique - plus les résidus sont élevés, moins la prévisibilité.
Le tableau ci-dessous prĂ©sente la moyenne des observations et la moyenne des rĂ©sidus pour chaque produit, ainsi que la proportion de l’erreur par rapport Ă la sĂ©rie chronologique observĂ©e, pour l’ensemble des mĂ©dicaments l’erreur ne dĂ©passe pas 7% des ventes observĂ©es, chose qui indique qu’une quantitĂ© importante d’information pourra Ăªtre expliquĂ©e.
d = pd.DataFrame(0, index=sales_s.columns, columns=[
"RESMEAN", "OBSMEAN", "PERC"], dtype=float)
for col in sales_s.columns:
result = seasonal_decompose(sales_s[col], period=12, model='additive')
res, obs = result.resid, result.observed
d.loc[col][:2] = list(map(lambda x: np.mean(
np.abs(x)), (res, obs[~np.isnan(res)])))
d.PERC = d.RESMEAN*100/d.OBSMEAN
d.round(2)
sales_s.to_pickle("sales_s.pkl")
def choose_d(l):
return list(sorted(l))[1]
def Oui_non(x):
return "Oui" if x else "Non"
def test_all(y):
y = y.dropna()
n_ch = nsdiffs(remove_trend(
y)["serie"].dropna(), m=12, max_D=12, test='ch')
n_ocsb = nsdiffs(remove_trend(
y)["serie"].dropna(), m=12, max_D=12, test='ocsb')
n_adf = ndiffs(y, test='adf')
n_kpss = ndiffs(y, test='kpss')
n_pp = ndiffs(y, test='pp')
return (Oui_non(n_ch), Oui_non(n_ocsb), max(n_ch, n_ocsb), Oui_non(n_adf),
Oui_non(n_kpss), Oui_non(n_pp), choose_d((n_adf, n_kpss, n_pp)))
def multiple_test_all(df):
res = pd.DataFrame(df.columns, columns=["Series"])
res[1], res[2], res[3], res[4], res[5], res[6], res[7] = [
[test_all(df[i])[j] for i in df.columns] for j in range(7)]
res = res.set_index("Series")
res.columns = ["Canova-Hansen", "OCSB", "Seasonality Order",
"ADF", "KPSS", "PP", "differencing order"]
return res
def remove_trend(y):
n_adf = ndiffs(y, test='adf')
n_kpss = ndiffs(y, test='kpss')
n_pp = ndiffs(y, test='pp')
n_diff = choose_d((n_adf, n_kpss, n_pp))
for i in range(n_diff):
y = y.diff(1)
return {"serie": y, "n_diff": n_diff}
def stationnarise(y):
y, n_diff = remove_trend(y).values()
y = y.dropna()
n_ch = nsdiffs(y, m=12, max_D=12, test='ch')
n_ocsb = nsdiffs(y, m=12, max_D=12, test='ocsb')
n_sdiff = max(n_ch, n_ocsb)
for i in range(n_sdiff):
y = y.diff(12)
return {"serie": y, "n_diff": n_diff, "n_sdiff": n_sdiff}
def series_from_summary(summary):
return pd.DataFrame([dic['serie'] for dic in summary.values()]).T
def excess_short(forecast, actual):
epsilon = forecast - actual
mape = np.mean(np.abs(epsilon)/np.abs(actual))
exces_mask = epsilon > 0
short_mask = epsilon < 0
mape_exces = np.sum(np.abs(
forecast[exces_mask] - actual[exces_mask])/np.abs(actual[exces_mask]))/len(forecast)
mape_short = np.sum(np.abs(
forecast[short_mask] - actual[short_mask])/np.abs(actual[short_mask]))/len(forecast)
return 100*mape_exces/mape, 100*mape_short/mape
def forecast_accuracy(forecast, actual):
assert len(forecast) == len(actual)
mape = mean_absolute_percentage_error(actual, forecast)*100
mae = mean_absolute_error(actual, forecast) # MAE
rmse = mean_squared_error(actual, forecast, squared=True) # RMSE
mape_exces, mape_short = excess_short(forecast, actual)
return pd.Series({'mape': mape, 'mae': mae, 'rmse': rmse, 'mape_exces': mape_exces, 'mape_short': mape_short})
\ L'analyse d'autocorrélation illustre le potentiel de prédiction des données de séries chronologiques. Les graphiques d'autocorrélation résument graphiquement la force d'une relation avec une observation dans une série chronologique avec des observations à des pas de temps antérieurs. Le coefficient de Pearson est utilisé pour mesurer l'autocorrélation. Ainsi, l'analyse suivante n'est pertinente que pour les données avec une distribution gaussienne normale.
Un tracé de l'autocorrélation d'une série chronologique par décalage est appelé la fonction d'autocorrélation (ACF). Ce graphique est parfois appelé corrélogramme ou graphique d'autocorrélation. Le graphique montre la valeur de décalage le long de l'axe des x et la corrélation sur l'axe des y entre -1 et 1. Les intervalles de confiance sont dessinés sous forme de cône. Par défaut, il est défini sur un intervalle de confiance de 95%, ce qui suggère que les valeurs de corrélation en dehors de ce code sont très probablement une corrélation.
La figure ci-dessous présente la fonction d’autocorrélation (ACF) pour chacun des médicaments avec un intervalle de confiance de 95%.
fig, axes = plt.subplots(5, 2, figsize=(20, 30), sharex=False)
axes[-1, -1].axis('off')
axs = axes.flatten()
alpha = .05
for i, col in enumerate(sales_s.columns):
sm.graphics.tsa.plot_acf(sales_s[col].values.squeeze(
), lags=40, ax=axs[i], title=f"{col} ACF", alpha=alpha)
plt.show()
D'après les correlogrames des séries temporelles ci dessus, on peut constater que la fonction d'autocorrelation se dégrade ou diminue lentement, ce qui indique la non stationnarité des séries. De plus, dans les deux séries des produits p_47 et p_64 l'ACF montre une oscillation et des pics qui se reproduisent avec des décalages de 12 mois indicant un effet de saisonnalité.
Tests de stationnarité¶
Afin d'estimer le terme de différenciation saisonnière D, on a effectué deux tests de saisonnalité pour différents niveaux de D pour estimer le nombre de différences saisonnières nécessaires pour rendre la série chronologique stationnaire et on a sélectionné la valeur maximale de D pour laquelle la série chronologique est jugée saisonnièrement stationnaire par le test statistique.
De mĂªme, pour estimer le terme de diffĂ©renciation d, on a effectuĂ© trois test statistiques de racine unitaire : Augmented Dickey Fuller (ADF), Phillips–Perron (PP) et Kwiatkowski–Phillips–Schmidt–Shin (KPSS) pour diffĂ©rents niveaux de d afin d'estimer le nombre de diffĂ©rences nĂ©cessaires pour rendre chaque sĂ©rie chronologique stationnaire, et on a sĂ©lectionnĂ© la valeur maximale de d pour laquelle la sĂ©rie temporelle est jugĂ©e stationnaire par le test statistique.
Le tableau ci dessous synthĂ©tise les rĂ©sultats, chaque test de saisonnalitĂ© indique la prĂ©sence ou absence d'un effet de saisonnaliĂ©, et de mĂªme pour les tests de racine unitaire, ainsi que les ordres de differenciation nĂ©caissaires pour stationnarier la sĂ©rie.
multiple_test_all(sales_s)
Stationnarisation des séries¶
On a procedĂ© par la suite Ă la stationnarisation des sĂ©ries en appliquant la differenciation adĂ©quate, et on a reĂ©ffectuĂ© les mĂªmes tests pour confirmer qu'on a atteint la stationnaritĂ© des sĂ©ries.
stationnarisation_summary = {key: 0 for key in sales_s.columns}
for col in sales_s.columns:
stationnarisation_summary[col] = stationnarise(sales_s[col])
stationnaries = series_from_summary(stationnarisation_summary)
fig, axs = plt.subplots(5, 2, figsize=(18, 18))
remove_axes(len(sales_s.columns), axs)
axs = axs.flatten()
cols = [next(color_cycle)]*3
for i, item in enumerate(sales_s.columns):
stationnaries[item].plot(title=f"stationnarised {item}",
color=next(color_cycle),
ax=axs[i])
stationnaries[item].rolling(12).mean().plot(
color=cols[1],
ax=axs[i], label="rolling mean")
stationnaries[item].rolling(12).std().plot(
color=cols[2],
ax=axs[i], linestyle='--', label="rolling std")
axs[i].legend()
plt.show()
multiple_test_all(stationnaries)
Les tests finals confirment qu’on a bien atteint la stationnarité des séries.
fig, axes = plt.subplots(9, 2, figsize=(20, 35), sharex=False)
alpha = .05
for i, col in enumerate(sales_s.columns):
sm.graphics.tsa.plot_acf(stationnaries[col].dropna().values.squeeze(
), lags=40, ax=axes[i, 0], title=f"{col} ACF", alpha=alpha)
sm.graphics.tsa.plot_pacf(stationnaries[col].dropna().values.squeeze(
), lags=15, ax=axes[i, 1], title=f"{col} PACF", alpha=alpha)
plt.show()
A l’addition de la stationnaritĂ© des sĂ©ries qui peut Ăªtre confirmĂ© Ă partir de ces corrĂ©logrammes, ces derniers nous permettent aussi de dĂ©terminer l'ordre d'ARIMA, mais la dĂ©tection visuelle ne permet pas toujours d'atteindre la performance maximale, c'est pourquoi on a utilisĂ© la fonction Auto_Arima qui permet d'ajuste le meilleur modèle ARIMA Ă une sĂ©rie temporelle univariĂ©e selon un critère d'information fourni (soit AIC, AICc, BIC...). La fonction effectue une recherche (soit pas Ă pas, soit parallĂ©lisĂ©e) sur des ordres de modèle et saisonniers possibles dans les contraintes fournies, et sĂ©lectionne les paramètres qui minimisent la mĂ©trique donnĂ©e.# L'approche parallèle est une mĂ©thode de brute force naĂ¯ve, il s'agit d'une recherche de grille (grid search) sur diverses combinaisons d'hyperparamètres. Cela prendra gĂ©nĂ©ralement plus de temps pour plusieurs raisons. Tout d'abord, il n'y a pas de procĂ©dure intelligente sur la façon dont les commandes de modèles sont testĂ©es ; ils sont tous testĂ©s (pas de court-circuit), ce qui peut prendre un certain temps.
L'approche par étapes (stepwise) suit la stratégie définie par Hyndman et Khandakar dans leur article de 2008, « Automatic Time Series Forecasting : The Forecast Package for R ».
Étape 1 : Essayez quatre modèles possibles pour commencer :
ARIMA(2,d,2) si m = 1 et ARIMA(2,d,2)(1,D,1) si m > 1\ ARIMA(0,d,0) si m = 1 et ARIMA(0,d,0)(0,D,0) si m > 1\ ARIMA(1,d,0) si m = 1 et ARIMA(1,d,0)(1,D,0) si m > 1\ ARIMA(0,d,1) si m = 1 et ARIMA(0,d,1)(0,D,1) si m > 1\ Le modèle avec le plus petit AIC (ou BIC, ou AICc, etc., selon les critères de minimisation) est sélectionné. C'est le "meilleur" modèle actuel.
Étape 2 : Considérez un certain nombre d'autres modèles :
OĂ¹ l'un des p, q, P et Q peut varier de ±1 par rapport au meilleur modèle actuel\ OĂ¹ p et q varient tous les deux de ±1 par rapport au meilleur modèle actuel\ OĂ¹ P et Q varient tous les deux de ±1 par rapport au meilleur modèle actuel\ Chaque fois qu'un modèle avec un critère d'information infĂ©rieur est trouvĂ©, il devient le nouveau meilleur modèle actuel, et la procĂ©dure est rĂ©pĂ©tĂ©e jusqu'Ă ce qu'il ne puisse pas trouver un modèle proche du meilleur modèle actuel avec un critère d'information infĂ©rieur ou que le processus dĂ©passe l'un des seuils d'exĂ©cution dĂ©fini via arima.StepwiseContext.
@meta_decorator("joblib_ARIMA_sales_s_Models.joblib")
def auto_arima_for_df(df, trace: bool):
"""find best arima model for each series in a data frame, returns dict of models"""
models = {}
for col in df.columns:
if trace:
print('\n', col, ':')
models[col] = auto_arima(df[col], start_p=1, start_q=1,
test='adf', # use adftest to find optimal 'd'
max_p=3, max_q=3, # maximum p and q
m=12, # frequency of series
# d=None, # let model determine 'd'
# seasonal=True, # No Seasonality
# start_P=0,
# D=0,
trace=trace,
# error_action='ignore',
# suppress_warnings=True,
stepwise=True)
return models
models = auto_arima_for_df(sales_s, trace=False)
Les meilleurs modèles sélectionnés sont :
res = pd.DataFrame(columns=['Best model', 'AIC', 'BIC'])
for col in sales_s.columns:
tbl = models[col].summary().tables[0].data
res = res.append({"Best model": tbl[1][1], "AIC": float(
tbl[2][3]), "BIC": float(tbl[3][3])}, ignore_index=True)
res.set_index(sales_s.columns)
Prévisions statiques et dynamiques¶
On a effectué la prévision de l'ensemble des séries par les meilleurs modèles ARIMA séléctionnés selon deux approches.
La méthode de prévision statique calcule une séquence de prévisions à un pas d'avance, en utilisant les valeurs réelles plutôt que prévues pour les variables dépendantes décalées, si elles sont disponibles. En effet, pour effectuer la prévision d'un nouveau mois, le modèle est mis a jour en incluant la valeur réelle du mois précédant, cette méthode est plus performante et donne des prévisions moins biaisées mais ne peut etre utilisée que pour prévoir un seul mois à l'avance.
La méthode de prévision dynamique calcule des prévisions dynamiques en plusieurs étapes à partir de la première période de l'échantillon de prévision, les valeurs précedemment prévues pour les variables dépendantes décalées sont utilisées pour former la prévision de la valeur actuelle. Cette méthode est utile lorsque l'entreprise s'interesse à effectuer une prévision à long terme, c'est dire prévoir les ventes de plusieurs mois à l'avance.
train_arima = sales_s.loc[:'2018-12-31']
test_arima = sales_s.loc['2019-01-31':]
def empty_df_like(df):
return pd.DataFrame(columns=df.columns, index=df.index)
def empty_df_to_date(df, n_periods):
return pd.DataFrame(columns=df.columns,
index=pd.date_range(df.index[-1], freq="M", periods=n_periods).shift(1, freq='M'))
def arima_dynamic_out_of_sample_forecast(models, train_df, n_periods):
"""enerate future dynamic forecast to the given periods"""
# fitted_dynamic = [0]*len(train_df.columns)
fc_df = empty_df_to_date(train_df, n_periods)
lower_df = empty_df_to_date(train_df, n_periods)
upper_df = empty_df_to_date(train_df, n_periods)
for col in train_df.columns:
fitted_dynamic = models[col].fit(train_df[col])
fc, confint = fitted_dynamic.predict(
n_periods=n_periods, return_conf_int=True)
fc_df[col] = fc
lower_df[col] = confint[:, 0]
upper_df[col] = confint[:, 1]
return {'fc_df': fc_df, 'lower_df': lower_df, 'upper_df': upper_df}
def arima_dynamic_in_sample_forecast(models, train_df, test_df):
"""generate dynamic forecast with confidence intervales and accuracy"""
n_series = len(train_df.columns)
n_periods = len(test_df)
# fitted_dynamic = [0]*n_series
accuracies_dynamic = pd.DataFrame(columns=test_df.columns, index=[
'mape', 'mae', 'rmse', 'mape_exces', 'mape_short'])
fc_df = empty_df_like(test_df)
lower_df = empty_df_like(test_df)
upper_df = empty_df_like(test_df)
for col in train_df.columns:
fitted_dynamic = models[col].fit(train_df[col])
fc, confint = fitted_dynamic.predict(
n_periods=n_periods, return_conf_int=True)
fc_df[col] = fc
lower_df[col] = confint[:, 0]
upper_df[col] = confint[:, 1]
accuracies_dynamic[col] = forecast_accuracy(fc_df[col], test_df[col])
return {'fc_df': fc_df, 'lower_df': lower_df, 'upper_df': upper_df, 'accuracies_dynamic': accuracies_dynamic}
def arima_static_in_sample_forecast(models, train_df, test_df):
"""generate dynamic forecast with confidence intervales and accuracy"""
n_series = len(train_df.columns)
n_periods = len(test_df)
# fitted_static = [0]*n_series
accuracies_static = pd.DataFrame(columns=test_df.columns, index=[
'mape', 'mae', 'rmse', 'mape_exces', 'mape_short'])
fc_df = empty_df_like(test_df)
lower_df = empty_df_like(test_df)
upper_df = empty_df_like(test_df)
for col in train_df.columns:
fitted_static = models[col].fit(train_df[col])
preds = np.zeros((len(test_df)))
conf_int = np.zeros((len(test_df), 2))
for j in range(len(test_df)):
new_preds, new_conf_int = fitted_static.predict(
n_periods=1, return_conf_int=True)
fitted_static.update([test_df[col].iloc[j]])
preds[j] = new_preds
conf_int[j] = new_conf_int
fc_df[col] = preds
lower_df[col] = conf_int[:, 0]
upper_df[col] = conf_int[:, 1]
accuracies_static[col] = forecast_accuracy(fc_df[col], test_df[col])
return {'fc_df': fc_df, 'lower_df': lower_df, 'upper_df': upper_df, 'accuracies_static': accuracies_static}
def plot_forecast_ARIMA(forecast_data, df, plot_anomalies=False):
n_series = len(df.columns)
fig_n_lines = ceil(n_series/2)
# fig, axs = plt.subplots(fig_n_lines, 2 if n_series>1 else 1, figsize=(20, fig_n_lines*4))
if n_series > 1:
fig, axs = plt.subplots(fig_n_lines, 2, figsize=(20, fig_n_lines*4))
remove_axes(n_series, axs)
axs = axs.flatten()
else:
fig, axs = plt.subplots(fig_n_lines, 1, figsize=(15, 7))
axs = [axs]
for i, item in enumerate(df.columns):
axs[i].plot(df[item], label="actual")
# axs[i].plot(pd.concat([df[item], forecast_data['fc_df'][item]]), label="actual")
axs[i].plot(forecast_data['fc_df'][item],
color='darkgreen', label="forecast")
axs[i].set_title(f"{item} forecast")
lower = forecast_data['lower_df'][item]
upper = forecast_data['upper_df'][item]
axs[i].plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
axs[i].plot(upper, "r--", alpha=0.5)
axs[i].axvspan(lower.index[0], lower.index[-1],
color=sns.xkcd_rgb['grey'], alpha=0.2)
if plot_anomalies:
index_of_pred = forecast_data['fc_df'].index
y_test = df[item][index_of_pred]
anomalies = np.array([np.NaN]*len(y_test))
anomalies[y_test < lower] = y_test[y_test < lower]
anomalies[y_test > upper] = y_test[y_test > upper]
anomalies_series = pd.Series(anomalies, index=index_of_pred)
if ~np.isnan(anomalies).all():
axs[i].plot(anomalies_series, "o",
markersize=10, label="Anomalies")
axs[i].legend(loc="best")
plt.show()
Les graphiques ci-dessous reprĂ©sentent les prĂ©visions effectuĂ©es par les deux mĂ©thodes ainsi que les intervalles de confiance (95%). On peut constater que les sĂ©ries des deux produits P_22 et P_30 correspondent Ă des marches alĂ©atoires, la meilleure prĂ©vision dans ce cas est la valeur prĂ©cendente. Les intervalles de prĂ©vision associĂ©s Ă la mĂ©thode dynamique s'allongent Ă mesure que l'horizon de prĂ©vision s'allonge, plus nous prĂ©voyons Ă l'avance, plus l'incertitude est associĂ©e Ă la prĂ©vision, et donc plus les intervalles de prĂ©vision sont larges. Dans les prĂ©visions par la mĂ©thode statique, la variation de la distribution des prĂ©visions est presque la mĂªme que la variation des rĂ©sidus.
Prévisions dynamiques dans l'échantillon¶
arima_dynamic_in_sample_forecast_res = arima_dynamic_in_sample_forecast(
models, train_arima, test_arima)
plot_forecast_ARIMA(arima_dynamic_in_sample_forecast_res, sales_s, True)
pretty_df(arima_dynamic_in_sample_forecast_res['accuracies_dynamic'])
Prévisions dynamiques hors échantillon¶
arima_dynamic_out_of_sample_forecast_res = arima_dynamic_out_of_sample_forecast(
models, sales_s, 12)
plot_forecast_ARIMA(arima_dynamic_out_of_sample_forecast_res, sales_s)
Prévisions statiques¶
arima_static_in_sample_forecast_res = arima_static_in_sample_forecast(
models, train_arima, test_arima)
errors_ARIMA = arima_static_in_sample_forecast_res['accuracies_static']
plot_forecast_ARIMA(arima_static_in_sample_forecast_res, sales_s, True)
pretty_df(errors_ARIMA)
Diagnostic des résidus¶
Pour la validation des modèles obtenus, on a effectué un diagnostic des résidus. le tableau ci-dessous synthétise les résultats des tests de bruit blanc, homoscedasticité et normalité pour le résidus de chacune des séries
def Verif(x):
return "Vérifié" if x else "Non Vérifié"
res = pd.DataFrame(columns=['Prob(Q)', 'White noise',
'Prob(H)', 'Homoscedasticity', 'Prob(JB)', 'Normality'])
for col in sales_s.columns:
tbl = models[col].summary().tables[2].data
res = res.append({"Prob(Q)": float(tbl[1][1]), "White noise": Verif(float(tbl[1][1]) > 0.05),
"Prob(H)": float(tbl[3][1]), "Homoscedasticity": Verif(float(tbl[3][1]) > 0.05),
"Prob(JB)": float(tbl[1][3]), "Normality": Verif(float(tbl[1][3]) > 0.05)}, ignore_index=True)
res.set_index(sales_s.columns)
Un diagnostic visuel des résidus est appliqué sur le produit P_2 à titre d'exemple.
models['P_2'].plot_diagnostics(figsize=(16, 8))
plt.show()
display(models[col].summary())
def double_exponential_smoothing(series, alpha=0.9, beta=0.9):
"""
series - dataset with timeseries
alpha - float [0.0, 1.0], smoothing parameter for level
beta - float [0.0, 1.0], smoothing parameter for trend
"""
# first value is same as series
result = [series[0]]
for n in range(1, len(series)+1):
if n == 1:
level, trend = series[0], series[1] - series[0]
if n >= len(series): # forecasting
value = result[-1]
else:
value = series[n]
last_level, level = level, alpha*value + (1-alpha)*(level+trend)
trend = beta*(level-last_level) + (1-beta)*trend
result.append(level+trend)
return result
month_in_year = 12
month_quarter = 4
window1 = 3
window2 = 6
window3 = 12
def add_vars(y):
df = y.copy()
y_name = df.columns[0]
df['year'] = df.index.year
df['month'] = df.index.month
df['quarter'] = df.index.quarter
df['month_sin'] = np.sin(2*np.pi*df['month']/month_in_year)
df['month_cos'] = np.cos(2*np.pi*df['month']/month_in_year)
df['quarter_sin'] = np.sin(2*np.pi*df['quarter']/month_quarter)
df['quarter_cos'] = np.cos(2*np.pi*df['quarter']/month_quarter)
for i in range(1, 13):
df[f"lag{i}"] = df[y_name].shift(i)
df["Date"] = y[[y_name]].index
df.reset_index(drop=True, inplace=True)
df_rolled_3d = df[[y_name]].rolling(window=window1, min_periods=0)
df_rolled_7d = df[[y_name]].rolling(window=window2, min_periods=0)
df_rolled_30d = df[[y_name]].rolling(window=window3, min_periods=0)
df_mean_3d = df_rolled_3d.mean().shift(1).reset_index().astype(np.float32)
df_mean_7d = df_rolled_7d.mean().shift(1).reset_index().astype(np.float32)
df_mean_30d = df_rolled_30d.mean().shift(1).reset_index().astype(np.float32)
df_std_3d = df_rolled_3d.std().shift(1).reset_index().astype(np.float32)
df_std_7d = df_rolled_7d.std().shift(1).reset_index().astype(np.float32)
df_std_30d = df_rolled_30d.std().shift(1).reset_index().astype(np.float32)
df[f"{y_name}_mean_lag{window1}"] = df_mean_3d[y_name]
df[f"{y_name}_mean_lag{window2}"] = df_mean_7d[y_name]
df[f"{y_name}_mean_lag{window3}"] = df_mean_30d[y_name]
df[f"{y_name}_std_lag{window1}"] = df_std_3d[y_name]
df[f"{y_name}_std_lag{window2}"] = df_std_7d[y_name]
df[f"{y_name}_std_lag{window3}"] = df_std_30d[y_name]
df[f"{y_name}_dbl_exp{'0.9, 0.9'}"] = pd.Series(
double_exponential_smoothing(df[y_name])).shift(1)
df.fillna(df.mean(), inplace=True)
df.set_index("Date", drop=False, inplace=True)
return df
def timeseries_train_test_split(X, y, test_size):
"""
Perform train-test split with respect to time series structure
"""
# get the index after which test set starts
test_index = int(len(X)*(1-test_size))
X_train = X.iloc[:test_index]
y_train = y.iloc[:test_index]
X_test = X.iloc[test_index:]
y_test = y.iloc[test_index:]
return X_train, X_test, y_train, y_test
def add_vars_plus(df, test_size):
"""
Add vars
"""
ts = df.copy()
ts = add_vars(ts)
y_name = ts.columns[0]
y = ts.dropna()[y_name]
X = ts.dropna().drop([y_name, 'Date'], axis=1)
return X, y
scaler = StandardScaler()
def prepare_data(df, test_size=0.15):
"""
Add vars & return X_train_scaled, X_test_scaled, y_train, y_test
"""
X, y = add_vars_plus(df, test_size)
X_train, X_test, y_train, y_test = timeseries_train_test_split(
X, y, test_size)
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
X_train_scaled = pd.DataFrame(X_train_scaled)
X_train_scaled.columns = X_train.columns
X_test_scaled = pd.DataFrame(X_test_scaled)
X_test_scaled.columns = X_train.columns
return X_train_scaled, X_test_scaled, y_train, y_test
# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)
from hyperopt import hp, fmin, tpe, STATUS_OK, Trials
def objective(space):
# Instantiate the classifier
clf = XGBRegressor(n_estimators=1000, colsample_bytree=space['colsample_bytree'],
learning_rate=space['learning_rate'],
max_depth=int(space['max_depth']),
min_child_weight=space['min_child_weight'],
subsample=space['subsample'],
gamma=space['gamma'])
# Fit the classsifier
clf.fit(X_train_scaled, y_train, verbose=False)
# Predict on Cross Validation data
pred = clf.predict(X_test_scaled)
# Calculate our Metric - mape
error = mean_absolute_percentage_error(y_test, pred)
return {'loss': error, 'status': STATUS_OK}
space = {'learning_rate': hp.uniform('learning_rate', 0.01, 0.3),
'max_depth': hp.quniform("max_depth", 4, 16, 1),
'min_child_weight': hp.quniform('min_child_weight', 1, 10, 1),
'subsample': hp.uniform('subsample', 0.7, 1),
'gamma': hp.uniform('gamma', 0.1, 0.5),
'colsample_bytree': hp.uniform('colsample_bytree', 0.7, 1)
}
def empty_df_like(df):
return pd.DataFrame(columns=df.columns, index=df.index)
@meta_decorator("joblib_XGB_sales_s_Models.joblib")
def get_tuned_xgb_for_df(df):
"""find best xgb model for each series in a data frame, returns dict of models"""
models = {}
for col in df.columns:
global X_train_scaled, y_train, X_test_scaled, y_test
X_train_scaled, X_test_scaled, y_train, y_test = prepare_data(df[[col]])
trials = Trials()
best = fmin(fn=objective,
space=space,
algo=tpe.suggest,
max_evals=100,
trials=trials)
best['max_depth'] = int(best['max_depth'])
models[col] = XGBRegressor(n_estimators=1000, **best)
return models
xgb_models = get_tuned_xgb_for_df(sales_s)
def ML_forecast(series, models, test_size=0.15):
test_index = int(len(series)*(1-test_size))
n_models = len(models)
models_errors = [pd.DataFrame(columns=series.columns, index=[
'mape', 'mae', 'rmse', 'mape_exces', 'mape_short']) for x in range(n_models)]
models_predictions = [empty_df_like(series) for x in range(n_models)]
models_lower = [empty_df_like(series) for x in range(n_models)]
models_upper = [empty_df_like(series) for x in range(n_models)]
test = [empty_df_like(sales_s) for x in range(2)]
for i, col in enumerate(series.columns):
global X_train_scaled, y_train, X_test_scaled, y_test
X_train_scaled, X_test_scaled, y_train, y_test = prepare_data(
series[[col]], test_size)
if i == 0:
models_coeffs = [pd.DataFrame(
columns=series.columns, index=X_train_scaled.columns) for x in range(n_models)]
for j in range(n_models):
if "xgboost" in models[j][0]:
if "tune" in models[j][0]:
models[j][1] = xgb_models[col]
models[j][1].fit(X_train_scaled, y_train, verbose=False)
models_coeffs[j][col] = models[j][1].feature_importances_
else:
models[j][1].fit(X_train_scaled, y_train)
models_coeffs[j][col] = models[j][1].coef_
models_predictions[j][col] = models[j][1].predict(
pd.concat([X_train_scaled, X_test_scaled]))
models_errors[j][col] = forecast_accuracy(
models_predictions[j][col][test_index:], y_test)
cv = cross_val_score(models[j][1], X_train_scaled, y_train,
cv=tscv,
scoring="neg_mean_absolute_error")
mae = cv.mean() * (-1)
deviation = cv.std()
scale = 1.96
models_lower[j][col] = models_predictions[j][col] - \
(mae + scale * deviation)
models_upper[j][col] = models_predictions[j][col] + \
(mae + scale * deviation)
models_predictions = {models[i][0]: pred for i,
pred in enumerate(models_predictions)}
models_lower = {models[i][0]: pred for i, pred in enumerate(models_lower)}
models_upper = {models[i][0]: pred for i, pred in enumerate(models_upper)}
models_errors = {models[i][0]: pred for i,
pred in enumerate(models_errors)}
models_coeffs = {models[i][0]: pred for i,
pred in enumerate(models_coeffs)}
return {'models_predictions': models_predictions, 'models_lower': models_lower, 'models_upper': models_upper,
'models_errors': models_errors, 'models_coeffs': models_coeffs, 'test_index': test_index}
def best_models(ML_forecast_res, criterion):
df_template = lambda model_data: empty_df_like(list(ML_forecast_res[model_data].values())[0])
best_predictions, models_lower, models_upper, models_errors, models_coeffs = map(
df_template, ('models_predictions', 'models_lower', 'models_upper', 'models_errors', 'models_coeffs'))
ML_forecast_error = {model_name: df.loc[criterion]
for model_name, df in ML_forecast_res['models_errors'].items()}
ML_forecast_error_df = pd.DataFrame(ML_forecast_error)
best_models = ML_forecast_error_df.idxmin(axis=1)
for col in best_predictions.columns:
best_model = best_models[col]
best_predictions[col] = ML_forecast_res['models_predictions'][best_model][col]
models_lower[col] = ML_forecast_res['models_lower'][best_model][col]
models_upper[col] = ML_forecast_res['models_upper'][best_model][col]
models_errors[col] = ML_forecast_res['models_errors'][best_model][col]
models_coeffs[col] = ML_forecast_res['models_coeffs'][best_model][col]
return {'best_predictions': best_predictions, 'models_lower': models_lower, 'models_upper': models_upper,
'models_errors': models_errors, 'models_coeffs': models_coeffs, 'best_models': best_models,
'criterion': criterion, 'test_index': ML_forecast_res['test_index']}
def plotCoefficients(model_coeff, features_names, model_name, item, axs):
"""
Plots sorted coefficient values of the model
"""
coefs = pd.DataFrame(model_coeff, features_names)
coefs.columns = ["coef"]
coefs["abs"] = coefs.coef.apply(np.abs)
coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)
axs.bar(coefs.index, coefs.coef)
axs.tick_params(axis='x', labelrotation=90)
axs.grid(True, axis='y')
axs.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')
axs.set_title(
f"{item}, Features importance of {'xgboost' if 'xgboost' in model_name else model_name}")
def plot_forecast_ML(forecast_data, df, plot_train_fc=False, plot_anomalies=True, plot_coeff=False):
n_series = len(df.columns)
test_index = forecast_data['test_index']
fig_n_lines = ceil(n_series/2)
if n_series > 1:
fig, axs = plt.subplots(fig_n_lines, 2, figsize=(20, fig_n_lines*4))
remove_axes(n_series, axs)
axs = axs.flatten()
else:
fig, axs = plt.subplots(fig_n_lines, 1, figsize=(15, 7))
axs = [axs]
for i, item in enumerate(df.columns):
axs[i].plot(df[item], label="actual")
axs[i].plot(forecast_data['best_predictions'][item]
[0 if plot_train_fc else test_index:], color='darkgreen', label="forecast")
axs[i].set_title(f"{item}, model : {'xgboost' if 'xgboost' in forecast_data['best_models'][item] else forecast_data['best_models'][item]}, {forecast_data['criterion']} : {forecast_data['models_errors'][item][forecast_data['criterion']]:.2f} { '%' if forecast_data['criterion']=='mape' else ''}")
lower = forecast_data['models_lower'][item][test_index:]
upper = forecast_data['models_upper'][item][test_index:]
axs[i].plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
axs[i].plot(upper, "r--", alpha=0.5)
y_test = df[item][test_index:]
if plot_anomalies:
anomalies = np.array([np.NaN]*len(y_test))
anomalies[y_test < lower] = y_test[y_test < lower]
anomalies[y_test > upper] = y_test[y_test > upper]
anomalies_series = pd.Series(anomalies, index=y_test.index)
if ~np.isnan(anomalies).all():
axs[i].plot(anomalies_series[-len(y_test):], "o",
markersize=10, label="Anomalies")
axs[i].axvspan(lower.index[0], lower.index[-1],
color=sns.xkcd_rgb['grey'], alpha=0.2)
axs[i].legend(loc="best")
# plot features importance
if plot_coeff:
features_names = forecast_data['models_coeffs'].index
if n_series > 1:
fig, axs = plt.subplots(
fig_n_lines, 2, figsize=(16, fig_n_lines*5))
remove_axes(n_series, axs)
axs = axs.flatten()
else:
fig, axs = plt.subplots(fig_n_lines, 1, figsize=(13, 7))
axs = [axs]
for i, item in enumerate(df.columns):
plotCoefficients(forecast_data['models_coeffs'][item], features_names,
forecast_data['best_models'][item], item, axs[i])
fig.tight_layout()
plt.show()
def full_forecast(series, models, test_size=0.15, criterion='mape', plot_train_fc=False, plot_anomalies=True, plot_coeff=False, return_errors=False):
ML_forecast_res = ML_forecast(sales_s, models, test_size)
res = best_models(ML_forecast_res, criterion)
plot_forecast_ML(res, series, plot_train_fc, plot_anomalies, plot_coeff)
if return_errors:
return ML_forecast_res['models_errors']
lr = LinearRegression()
ridge = RidgeCV(cv=tscv)
lasso = LassoCV(cv=tscv)
Pour la démonstration de la demarche suivie dans l'apprentissage supervisé, on a selectionné aléatoirement un produit (P_64).
sales_s_p64 = sales_s[['P_64']]
sales_s_p64.head()
On a commencé par implémenter le concept de "feature engineering" ou on a généré les variables nécessaires pour entrainer les modèles d'apprentissage.
- Variables de temps : année, mois, trimestre
- Valeurs retardées des séries temporelles
Statistiques : la moyenne et l'écart type des valeurs retardées.
Les variables temporelles tels que le mois et le trimestre sont cycliques. Par exemple, le mois oscille entre 1 et 12 pour chaque année. Alors que la différence entre chaque mois augmente de 1 au cours de l'année, entre deux ans, la caractéristique du mois passe de 12 (décembre) à 1 (janvier). Il en résulte une différence de -11, ce qui peut dérouter beaucoup de modèles.
demo = sales_s_p64.copy()
demo['year'] = demo.index.year
demo['month'] = demo.index.month
demo['quarter'] = demo.index.quarter
demo.tail()
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 3))
sns.lineplot(x=demo.index, y=demo['quarter'], color='dodgerblue')
# ax.set_xlim(["2014-07-31", "2019-09-30"])
plt.show()
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 3))
sns.lineplot(x=demo.index, y=demo['month'], color='dodgerblue')
# ax.set_xlim(["2014-07-31", "2019-09-30"])
plt.show()
Pour rĂ©soudre ce problème, on a effectuĂ© une transformation par la fonction cosinus après avoir normalisĂ© les variables entre 0 et 2Ï€, ce qui correspond Ă un cycle de cosinus, cette solution ne rĂ©gle pas le problème complĂ©tement car deux valeurs diffĂ©rentes peuvent avoir la mĂªme image par la fonction cosinus. La meilleure façon de rĂ©soudre ce nouveau problème Ă©tait d'ajouter une autre information cyclique pour distinguer deux temps avec des valeurs de cosinus identiques, il s'agit de la fonction sinus. Nous pourrions le considĂ©rer comme un système de coordonnĂ©es Ă deux axes.
month_in_year = 12
demo['month_sin'] = np.sin(2*np.pi*demo['month']/month_in_year)
demo['month_cos'] = np.cos(2*np.pi*demo['month']/month_in_year)
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 6))
sns.scatterplot(x=demo.month_sin, y=demo.month_cos, color='dodgerblue')
plt.show()
month_quarter = 4
demo['quarter_sin'] = np.sin(2*np.pi*demo['quarter']/month_quarter)
demo['quarter_cos'] = np.cos(2*np.pi*demo['quarter']/month_quarter)
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 6))
sns.scatterplot(x=demo.quarter_sin, y=demo.quarter_cos, color='dodgerblue')
plt.show()
demo.head()
Comme indiqué précédemment, on a utilisé les variables de la série décalées pour beneficier de l'autocorrélation existante dans la série. Pour prendre en compte l'effet de la saisonnalité, on a considéré 12 retards.
for i in range(1, 13):
demo[f"lag{i}"] = demo.P_64.shift(i)
demo2 = demo.copy()
demo2["Date"] = demo[['P_64']].index
On a généré les moyennes mobiles et ecart-types mobiles (trimestrielle, semestrielle et annuelle) des valeurs décalées afin de détecter la tendance globale des données.
demo2.reset_index(drop=True, inplace=True)
window1 = 3
window2 = 6
window3 = 12
df_rolled_3d = demo2[["P_64"]].rolling(window=window1, min_periods=0)
df_rolled_7d = demo2[["P_64"]].rolling(window=window2, min_periods=0)
df_rolled_30d = demo2[["P_64"]].rolling(window=window3, min_periods=0)
df_mean_3d = df_rolled_3d.mean().shift(1).reset_index().astype(np.float32)
df_mean_7d = df_rolled_7d.mean().shift(1).reset_index().astype(np.float32)
df_mean_30d = df_rolled_30d.mean().shift(1).reset_index().astype(np.float32)
df_std_3d = df_rolled_3d.std().shift(1).reset_index().astype(np.float32)
df_std_7d = df_rolled_7d.std().shift(1).reset_index().astype(np.float32)
df_std_30d = df_rolled_30d.std().shift(1).reset_index().astype(np.float32)
demo2[f"{'P_64'}_mean_lag{window1}"] = df_mean_3d["P_64"]
demo2[f"{'P_64'}_mean_lag{window2}"] = df_mean_7d["P_64"]
demo2[f"{'P_64'}_mean_lag{window3}"] = df_mean_30d["P_64"]
demo2[f"{'P_64'}_std_lag{window1}"] = df_std_3d["P_64"]
demo2[f"{'P_64'}_std_lag{window2}"] = df_std_7d["P_64"]
demo2[f"{'P_64'}_std_lag{window3}"] = df_std_30d["P_64"]
demo2[f"{'P_64'}_dbl_exp{'0.9, 0.9'}"] = pd.Series(
double_exponential_smoothing(demo2["P_64"])).shift(1)
demo2.fillna(demo2.mean(), inplace=True)
demo2.set_index("Date", drop=True, inplace=True)
demo2.head()
full_forecast(sales_s_p64, [["Linear Regression", lr]],
plot_train_fc=True, plot_coeff=True)
On peut observer que la courbe de régression est collée à la série temporelle du produit, ce qui indique qu’il s’agit d’un modèle complexe et flexible qui surapprend les données d’entrainement (un biais faible) mais conduit à des erreurs plus importantes sur les nouvelles données (une variance élevée). La figure de l’importance des variables confirme le surapprentissage. En effet, le modèle concentre toutes ses forces sur les variables du temps et néglige l’information inclue dans les autres variables. Lorsque les variables indépendantes sont corrélées, cela indique que les changements dans une variable sont associés à des changements dans une autre variable. Plus la corrélation est forte, plus il est difficile de changer une variable sans en changer une autre. Il devient difficile pour le modèle d'estimer la relation entre chaque variable indépendante et la variable dépendante indépendamment parce que les variables indépendantes ont tendance à changer à l'unisson. Dans notre cas, plusieurs variables sont fortement corrélées comme l’indique le corrélogramme suivant :
cg = sns.clustermap(demo2.corr(), figsize=(
10, 10), cbar_pos=(.1, .1, .03, .65))
cg.ax_row_dendrogram.set_visible(False)
cg.ax_col_dendrogram.set_visible(False)
plt.show()
D'autre part, puisqu’on a diffĂ©rentes Ă©chelles dans les variables, celles-ci doivent Ăªtre transformĂ©es en une mĂªme Ă©chelle pour explorer l'importance des caractĂ©ristiques et Ăªtre utilisĂ©s par la suite dans la rĂ©gularisation.
modèles de régularisation : Ridge et Lasso¶
Pour une meilleure optimisation des variables explicatives, on a appliqué la régularisation. Deux des modèles de régression avec régularisation les plus populaires sont les régressions Ridge et Lasso. Ils ajoutent tous deux des contraintes supplémentaires à notre fonction de perte.
Dans le cas de la régression Ridge, ces contraintes sont la somme des carrés des coefficients multipliée par le coefficient de régularisation. Plus le coefficient d'une caractéristique est grand, plus notre perte sera importante. Par conséquent, nous essaierons d'optimiser le modèle tout en gardant les coefficients assez bas.
À la suite de cette régularisation L2, nous aurons un biais plus élevé et une variance plus faible, donc le modèle se généralisera mieux.
Le deuxième modèle de régression, la régression Lasso, ajoute à la fonction de perte, non pas des carrés, mais des valeurs absolues des coefficients. Par conséquent, au cours du processus d'optimisation, les coefficients des caractéristiques sans importance peuvent devenir des zéros, ce qui permet une sélection automatisée des caractéristiques. Ce type de régularisation est appelé L1. Les prévisions par les deux régularisations ainsi que le graphique de l’importance des variables sont présentées dans les figures suivantes :
full_forecast(sales_s_p64, [["ridge", ridge]],
plot_train_fc=True, plot_coeff=True)
On peut remarquer que la qualitĂ© d’ajustement s’est amĂ©liorĂ©e par la rĂ©gression Ridge, l’erreur sur le train set est la mĂªme que l’erreur sur le test set, donc le modèle ne surapprend pas les donnĂ©es d’entrainement et le problème de surapprentissage est rĂ©glĂ©. On peut Ă©galement remarquer que les intervalles de prĂ©diction sont moins larges que dans la rĂ©gression linĂ©aire. Nous pouvons clairement voir que certains coefficients se rapprochent de plus en plus de zĂ©ro (bien qu'ils ne l'atteignent jamais rĂ©ellement) Ă mesure que leur importance dans le modèle diminue.
full_forecast(sales_s_p64, [["lasso", lasso]],
plot_train_fc=True, plot_coeff=True)
De mĂªme que la rĂ©gression Ridge, la rĂ©gression Lasso permet de rĂ©gler le problème de surapprentissage. La rĂ©gression au lasso s'est avĂ©rĂ©e plus conservatrice, le modèle elle a supprimĂ© plusieurs variables.
La figure ci-dessous représente les prévisions pour chaque produit par le modèle qui minimise l'erreur (MAPE) parmis les deux modèles de régularisation.
full_forecast(sales_s, [["Lasso", lasso], ["Ridge", ridge]])
Gradient Boosting : XGBoost¶
On a effectuĂ© les prĂ©visions par un autre modèle d’apprentissage, Ă savoir XGBoost. Ce modèle est très performant sur beaucoup de problèmes de prĂ©diction, Ă©valuons sa performance sur notre problème de prĂ©vision des sĂ©ries temporelles. L’implĂ©mentation XGBoost du grandient Boosting propose plusieurs hyperparamètres tel que : -Max_depth : Profondeur maximale d'un arbre. L'augmentation de cette valeur rendra le modèle plus complexe et plus susceptible de surapprentissage. -colsample_bytree est le sous-Ă©chantillon des colonnes considĂ©rĂ© lors de la construction de chaque arbre. Le sous-Ă©chantillonnage a lieu une fois pour chaque arbre construit. -max_leaves :Nombre maximum de nÅ“uds Ă ajouter. On considère encore une fois le produit P_64, on a essayĂ© plusieurs combinaisons des hyperparamètres jusqu’à trouver la combinaison optimale qui minimise l’erreur (MAPE), on a entrainĂ© le modèle en utilisant cette dernière, les prĂ©visions donnĂ©es sont prĂ©sentĂ©es dans le graphique ci-dessous.
best = {'colsample_bytree': 0.7855436372368847,
'gamma': 0.2715026337470021,
'max_depth': 11,
'min_child_weight': 1.0,
'reg_lambda': 0.358974472027333,
'subsample': 0.8129759325416941}
reg = XGBRegressor(n_estimators=1000, **best)
full_forecast(sales_s_p64, [["default_xgboost", reg]], plot_coeff=True)
On a réeffectué le tunning pour chacun des produits, choisi la combinaison des hyperparamètres qui minimise l’erreur et ajusté un modèle adapté. Le tableau ci-dessous présente l’erreur des modèles entrainés et testé sur le test set pour chaque produit.
full_forecast(sales_s, [["tuned_xgboost", 0], [
"default_xgboost", XGBRegressor(n_estimators=1000)]])
Finalement, pour chaque produit on a retenu le modèle qui minimise l’erreur sur le test set parmi les modèles d’apprentissage supervisé testés : RLM, Ridge, Lasso, XGBoost. Les prévisions données par le modèle qui minimise l’erreur pour chaque produit ainsi que son erreur sur le test set sont présentées dans le graphique suivant.
errors_ML = full_forecast(sales_s, [["tuned_xgboost", 0], ["default_xgboost", XGBRegressor(
n_estimators=1000)], ["Lasso", lasso], ["Ridge", ridge]], return_errors=True)
errors_LASSO = errors_ML['Lasso']
errors_RIDGE = errors_ML['Ridge']
errors_XGB = errors_ML['tuned_xgboost'].combine(
errors_ML['default_xgboost'], np.minimum)
Comparaison des modèles univariés¶
print("Performance des modèles ARIMA")
display(pretty_df(errors_ARIMA))
print("Performance des modèles Lasso")
display(pretty_df(errors_LASSO))
print("Performance des modèles Ridge")
display(pretty_df(errors_RIDGE))
print("Performance des modèles XGBoost")
display(pretty_df(errors_XGB))
Modélisation multivariée : Vector AutoRegressif Moving Average (VARMA)¶
Dans cette section, on a exploitĂ© le clustering effectuĂ© auparavant en s’intĂ©ressant Ă l’ensemble des produits d’un mĂªme groupe (Cluster). Ayant des comportements de ventes similaires, ces produits peuvent s'influencer mutuellement et par la suite l’information supplĂ©mentaire sur les Ă©volutions des autres produits peut amĂ©liorer la qualitĂ© du modèle. On dispose de 9 clusters, Pour la dĂ©monstration on considère le cluster qui comporte les produits :P_3, P_21, P_23, P_24, P_25, P_27. Visuellement, on peut constater que les produits ont des Ă©volutions similaires avec des effets de saisonnalitĂ© semblables.
def choose_d(l):
return list(sorted(l))[1]
def Oui_non(x):
return "Oui" if x else "Non"
def test_all(y):
y = y.dropna()
n_ch = nsdiffs(remove_trend(
y)["serie"].dropna(), m=12, max_D=12, test='ch')
n_ocsb = nsdiffs(remove_trend(
y)["serie"].dropna(), m=12, max_D=12, test='ocsb')
n_adf = ndiffs(y, test='adf')
n_kpss = ndiffs(y, test='kpss')
n_pp = ndiffs(y, test='pp')
return (Oui_non(n_ch), Oui_non(n_ocsb), max(n_ch, n_ocsb), Oui_non(n_adf), Oui_non(n_kpss), Oui_non(n_pp), choose_d((n_adf, n_kpss, n_pp)))
def multiple_test_all(df):
res = pd.DataFrame(df.columns, columns=["Series"])
res[1], res[2], res[3], res[4], res[5], res[6], res[7] = [
[test_all(df[i])[j] for i in df.columns] for j in range(7)]
res = res.set_index("Series")
res.columns = ["Canova-Hansen", "OCSB", "Seasonality Order",
"ADF", "KPSS", "PP", "differencing order"]
return res
def remove_trend(y):
n_adf = ndiffs(y, test='adf')
n_kpss = ndiffs(y, test='kpss')
n_pp = ndiffs(y, test='pp')
n_diff = max(n_adf, n_kpss, n_pp)
for i in range(n_diff):
y = y.diff(1)
return {"serie": y, "n_diff": n_diff}
def stationnarise(y):
y, n_diff = remove_trend(y).values()
y = y.dropna()
n_ch = nsdiffs(y, m=12, max_D=12, test='ch')
n_ocsb = nsdiffs(y, m=12, max_D=12, test='ocsb')
n_sdiff = max(n_ch, n_ocsb)
for i in range(n_sdiff):
y = y.diff(12)
return {"serie": y, "n_diff": n_diff, "n_sdiff": n_sdiff}
def series_from_summary(summary):
return pd.DataFrame([dic['serie'] for dic in summary.values()]).T
def orders_from_summary(summary):
return [(dic["n_diff"], dic["n_sdiff"]) for dic in summary.values()]
def diff_inv_wrapper(train, diff_test, lags=1, differences=1):
last_train = train[-lags:]
added_orig = pd.concat([last_train, diff_test])
X_diff_inv = diff_test.copy()
X_diff_inv.iloc[:, 0] = diff_inv(
added_orig, lag=lags, differences=differences)[-len(diff_test):]
return X_diff_inv
def undiff_d1_D1(train, diff_test):
diff_train_s = train.diff(1)[-12:]
unseasonal_test = diff_inv_wrapper(diff_train_s, diff_test, lags=12)
return diff_inv_wrapper(train, unseasonal_test, lags=1)
def invert_transformation(train, test, d=1, D=0):
if d == 1 and D == 1:
return undiff_d1_D1(train, test)
elif d == 1 and D == 0:
return diff_inv_wrapper(train, test)
elif d == 0 and D == 1:
return diff_inv_wrapper(train, test, lags=12)
elif d == 0 and D == 0:
return test
else:
raise ValueError("unexpected case")
def remove_axes(n_used_axes, all_axs):
if n_used_axes % 2:
all_axs = all_axs.flatten()
size = len(all_axs)
for useless_ax in all_axs.flatten()[-(size - n_used_axes):]:
useless_ax.axis('off')
cluster_count = len(np.unique(labels))
clusters = [[] for _ in range(cluster_count)]
for i, lab in enumerate(labels):
clusters[lab].append(i)
clusters_n = ["cluster_"+str(i) for i in range(cluster_count)]
cluster_p = {}
for cluster, cluster_n in zip(clusters, clusters_n):
cluster_p[cluster_n] = [sales.columns[p] for p in cluster]
cluster_7 = sales[cluster_p['cluster_7']]
cluster_7.head()
cluster_7.plot(subplots=True, layout=(3, 2), figsize=(20, 15), marker='.')
plt.show()
Le corrélogramme suivant confirme la forte corrélation entre les produits, plus particulièrement les produits P_21, P_23 et P_24.
cg = sns.clustermap(cluster_7.corr(), figsize=(
8, 8), cbar_pos=(.1, .1, .03, .65), annot=True)
cg.ax_row_dendrogram.set_visible(False)
cg.ax_col_dendrogram.set_visible(False)
plt.show()
fig, axes = plt.subplots(3, 2, figsize=(15, 12), sharex=False)
remove_axes(len(cluster_7.columns), axes)
axs = axes.flatten()
alpha = .05
for i, col in enumerate(cluster_7.columns):
sm.graphics.tsa.plot_acf(cluster_7[col].values.squeeze(
), lags=40, ax=axs[i], title=f"{col} ACF", alpha=alpha)
plt.show()
D’une manière analogue Ă ARIMA, le modèle VARMA nĂ©cessite la stationnaritĂ© des sĂ©ries temporelle. Pour vĂ©rifier la stationnaritĂ©, on a effectuĂ© les mĂªmes tests de racine unitaire et de saisonnalitĂ©, le tableau suivant synthĂ©tise les rĂ©sultats.
multiple_test_all(cluster_7)
Après les différenciations nécessaires, on obtient les séries suivantes :
stationnarisation_summary = {}
for col in cluster_7.columns:
stationnarisation_summary[col] = stationnarise(cluster_7[col])
cluster_7_stat = series_from_summary(stationnarisation_summary)
orders_7 = orders_from_summary(stationnarisation_summary)
cluster_7_stat.plot(subplots=True, layout=(3, 3), figsize=(20, 15), marker='.')
plt.show()
train_var = cluster_7.loc[:'2018-12-31']
test_var = cluster_7.loc['2019-01-31':]
Pour trouver le meilleur modèle VARMA, on a utilisé le module AutoTS, qui fonctionne actuellement de la manière suivante :
Le processus commence par le remodelage et le prĂ©-traitement de donnĂ©es selon le besoin Un premier fractionnement train/test est gĂ©nĂ©rĂ© oĂ¹ le test est de longueur 'Forecast_length' Le modèle initial est une combinaison d'apprentissage par transfert et de modèles gĂ©nĂ©rĂ©s alĂ©atoirement. Ceci est testĂ© sur le train/test initial Les modèles se composent d'une Ă©tape de prĂ©-transformation (options de remplissage, options de suppression des valeurs aberrantes, etc.) Les meilleurs modèles (sĂ©lectionnĂ©s par une combinaison de mĂ©triques) sont recombinĂ©s avec des mutations alĂ©atoires pour n_gĂ©nĂ©rations Un pourcentage des meilleurs modèles issus de ce processus est soumis Ă une validation croisĂ©e, oĂ¹ ils sont rĂ©Ă©valuĂ©s sur de nouvelles divisions train/test. S'il est utilisĂ©, l'assemblage horizontal utilise les donnĂ©es de validation pour choisir le meilleur modèle pour chaque sĂ©rie. Le meilleur modèle ou ensemble en cours de validation est sĂ©lectionnĂ© comme best_model et utilisĂ© dans la mĂ©thode .predict() pour gĂ©nĂ©rer des prĂ©visions.
cluster_7_filename = "cluster_7_export.json" # .csv/.json
best_model = AutoTS(frequency='M', max_generations=0,
num_validations=0, verbose=0)
try:
best_model = best_model.import_template(
cluster_7_filename, method='only', enforce_model_list=False)
print("template found and imported")
except Exception:
print("template not found", "Fitting new models ...", sep='\n')
model = AutoTS(
forecast_length=len(test_var),
frequency='M',
prediction_interval=0.95,
model_list=["VARMAX"],
ensemble="all",
# transformer_list="fast",
max_generations=8,
num_validations=2,
# models_to_validate=0.2,
validation_method="backwards"
)
model.fit(train_var)
model.export_template(cluster_7_filename, models='best', n=1)
best_model = best_model.import_template(
cluster_7_filename, method='only', enforce_model_list=False)
print("models fitted and ready to predict",
"template is saved to "+cluster_7_filename, sep='\n')
best_model.fit(train_var)
prediction = best_model.predict(forecast_length=len(test_var))
forecast_cluster_7 = prediction.forecast
lower_forecast_cluster_7 = prediction.lower_forecast
upper_forecast_cluster_7 = prediction.upper_forecast
Les prévisions par le modèle selectionné sur les séries donne les résultats suivants :
fig, axes = plt.subplots(3, 2, figsize=(15, 12), sharex=False)
remove_axes(len(forecast_cluster_7.columns), axes)
axs = axes.flatten()
alpha = .05
accuracies = pd.DataFrame(columns=cluster_7.columns, index=[
'mape', 'mae', 'rmse', 'mape_exces', 'mape_short'])
for i, col in enumerate(forecast_cluster_7.columns):
axs[i].plot(cluster_7[col], label="actual")
axs[i].plot(forecast_cluster_7[col], 'g', label="prediction")
axs[i].plot(lower_forecast_cluster_7[col], "r--",
label="upper bond / lower bond", alpha=0.5)
axs[i].plot(upper_forecast_cluster_7[col], "r--", alpha=0.5)
accuracies[col] = forecast_accuracy(forecast_cluster_7[col], test_var[col])
axs[i].set_title(
f"{col}, Mean absolute percentage error {accuracies[col]['mape']:.2f}%")
axs[i].axvspan(forecast_cluster_7.index[0], forecast_cluster_7.index[-1],
color=sns.xkcd_rgb['grey'], alpha=0.2)
anomalies = np.array([np.NaN]*len(test_var))
anomalies[test_var[col] < lower_forecast_cluster_7[col]
] = test_var[col][test_var[col] < lower_forecast_cluster_7[col]]
anomalies[test_var[col] > upper_forecast_cluster_7[col]
] = test_var[col][test_var[col] > upper_forecast_cluster_7[col]]
anomalies_series = pd.Series(anomalies, index=forecast_cluster_7.index)
if ~np.isnan(anomalies).all():
axs[i].plot(anomalies_series, "o", markersize=10, label="Anomalies")
axs[i].legend(loc="best")
plt.show()
accuracies