<jupyter_start><jupyter_text>TimeseriesCe notebook présente quelques étapes simples pour une série temporelle. La plupart utilise le module [statsmodels.tsa](https://www.statsmodels.org/stable/tsa.htmlmodule-statsmodels.tsa).<jupyter_code>from jyquickhelper import add_notebook_menu
add_notebook_menu()
%matplotlib inline<jupyter_output><empty_output><jupyter_text>DonnéesLes données sont artificielles mais simulent ce que pourraient être le chiffre d'affaires d'un magasin de quartier, des samedi très forts, une semaine morne, un Noël chargé, un été plat.<jupyter_code>from ensae_teaching_cs.data import generate_sells
import pandas
df = pandas.DataFrame(generate_sells())
df.head()<jupyter_output><empty_output><jupyter_text>Premiers graphiquesLa série a deux saisonnalités, hebdomadaire, mensuelle.<jupyter_code>import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(14, 4))
df.iloc[-30:].set_index('date').plot(ax=ax[0])
df.set_index('date').plot(ax=ax[1])
ax[0].set_title("chiffre d'affaire sur le dernier mois")
ax[1].set_title("chiffre d'affaire sur deux ans");<jupyter_output><empty_output><jupyter_text>Elle a une vague tendance, on peut calculer un tendance à l'ordre 1, 2, ...<jupyter_code>from statsmodels.tsa.tsatools import detrend
notrend = detrend(df.value, order=1)
df["notrend"] = notrend
df["trend"] = df['value'] - notrend
ax = df.plot(x="date", y=["value", "trend"], figsize=(14,4))
ax.set_title('tendance');<jupyter_output><empty_output><jupyter_text>Autocorrélations...<jupyter_code>from statsmodels.tsa.stattools import acf
cor = acf(df.value)
cor
fig, ax = plt.subplots(1, 1, figsize=(14,2))
ax.plot(cor)
ax.set_title("Autocorrélogramme");<jupyter_output><empty_output><jupyter_text>La première saisonalité apparaît, 7, 14, 21... Les autocorrélations partielles confirment cela, plutôt 7 jours.<jupyter_code>from statsmodels.tsa.stattools import pacf
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(df.value, lags=50);<jupyter_output><empty_output><jupyter_text>Comme il n'y a rien le dimanche, il vaut mieux les enlever. Garder des zéros nous priverait de modèles multiplicatifs.<jupyter_code>df["weekday"] = df.date.dt.weekday
df.head()
df_nosunday = df[df.weekday != 6]
df_nosunday.head(n=10)
fig, ax = plt.subplots(1, 1, figsize=(14,2))
cor = acf(df_nosunday.value)
ax.plot(cor)
ax.set_title("Autocorrélogramme");
plot_pacf(df_nosunday.value, lags=50);<jupyter_output><empty_output><jupyter_text>On décompose la série en tendance + saisonnalité. Les étés et Noël apparaissent.<jupyter_code>from statsmodels.tsa.seasonal import seasonal_decompose
res = seasonal_decompose(df_nosunday.value, freq=7)
res.plot();
plt.plot(res.seasonal[-30:])
plt.title("Saisonnalité");
cor = acf(res.trend[5:-5]);
plt.plot(cor);<jupyter_output><empty_output><jupyter_text>On cherche maintenant la saisonnalité de la série débarrassée de sa tendance herbdomadaire. On retrouve la saisonnalité mensuelle.<jupyter_code>res_year = seasonal_decompose(res.trend[5:-5], freq=25)
res_year.plot();<jupyter_output><empty_output><jupyter_text>Test de stationnaritéLe test [KPSS](https://en.wikipedia.org/wiki/KPSS_test) permet de tester la stationnarité d'une série.<jupyter_code>from statsmodels.tsa.stattools import kpss
kpss(res.trend[5:-5])<jupyter_output>c:\python372_x64\lib\site-packages\statsmodels\tsa\stattools.py:1276: InterpolationWarning: p-value is smaller than the indicated p-value
  warn("p-value is smaller than the indicated p-value", InterpolationWarning)<jupyter_text>Comme ce n'est pas toujours facile à interpréter, on simule une variable aléatoire gaussienne donc sans tendance.<jupyter_code>from numpy.random import randn
bruit = randn(1000)
kpss(bruit)<jupyter_output>c:\python372_x64\lib\site-packages\statsmodels\tsa\stattools.py:1278: InterpolationWarning: p-value is greater than the indicated p-value
  warn("p-value is greater than the indicated p-value", InterpolationWarning)<jupyter_text>Et puis une série avec une tendance forte.<jupyter_code>from numpy.random import randn
from numpy import arange
bruit = randn(1000) * 100 + arange(1000) / 10
kpss(bruit)<jupyter_output><empty_output><jupyter_text>Une valeur forte indique une tendance et la série en a clairement une.  PrédictionLes modèles *AR*, *ARMA*, *ARIMA* se concentrent sur une série à une dimension. En machine learning, il y a la série et plein d'autres informations. On construit une matrice avec des séries décalées.<jupyter_code>from statsmodels.tsa.tsatools import lagmat
lag = 8
X = lagmat(df_nosunday["value"], lag)
lagged = df_nosunday.copy()
for c in range(1,lag+1):
    lagged["lag%d" % c] = X[:, c-1]
lagged.tail()<jupyter_output><empty_output><jupyter_text>On ajoute ou on réécrit le jour de la semaine qu'on utilise comme variable supplémentaire.<jupyter_code>lagged["weekday"] = lagged.date.dt.weekday
X = lagged.drop(["date", "value", "notrend", "trend"], axis=1)
Y = lagged["value"]
X.shape, Y.shape
from numpy import corrcoef
corrcoef(X)<jupyter_output><empty_output><jupyter_text>Etrange autant de grandes valeurs, cela veut dire que la tendance est trop forte pour calculer des corrélations, il vaudrait mieux tout recommencer avec la série $\Delta Y_t = Y_t - Y_{t-1}$. Bref, passons...<jupyter_code>X.columns<jupyter_output><empty_output><jupyter_text>Une régression linéaire car les modèles linéaires sont toujours de bonnes baseline et pour connaître le modèle simulé, on ne fera pas beaucoup mieux.<jupyter_code>from sklearn.linear_model import LinearRegression
clr = LinearRegression()
clr.fit(X, Y)
from sklearn.metrics import r2_score
r2_score(Y, clr.predict(X))
clr.coef_<jupyter_output><empty_output><jupyter_text>On retrouve la saisonnalité, $Y_t$ et $Y_{t-6}$ sont de mèches.<jupyter_code>for i in range(1, X.shape[1]):
    print("X(t-%d)" % (i), r2_score(Y, X.iloc[:, i]))<jupyter_output>X(t-1) -0.5799344874123231
X(t-2) -1.0271995693021703
X(t-3) -1.2451181610331497
X(t-4) -1.0497781280554928
X(t-5) -0.6300696727754616
X(t-6) 0.8212380818165589
X(t-7) -0.643903114304551
X(t-8) -1.09689904382877<jupyter_text>Auparavant (l'année dernière en fait), je construisais deux bases, apprentissage et tests, comme ceci :<jupyter_code>n = X.shape[0]
X_train = X.iloc[:n * 2//3]
X_test = X.iloc[n * 2//3:]
Y_train = Y[:n * 2//3]
Y_test = Y[n * 2//3:]<jupyter_output><empty_output><jupyter_text>Et puis *scikit-learn* est arrivée avec [TimeSeriesSplit](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.TimeSeriesSplit.html).<jupyter_code>from sklearn.model_selection import TimeSeriesSplit
tscv = TimeSeriesSplit(n_splits=5)
for train_index, test_index in tscv.split(lagged):
    data_train, data_test = lagged.iloc[train_index, :], lagged.iloc[test_index, :]
    print("TRAIN:", data_train.shape, "TEST:", data_test.shape)<jupyter_output>TRAIN: (107, 13) TEST: (104, 13)
TRAIN: (211, 13) TEST: (104, 13)
TRAIN: (315, 13) TEST: (104, 13)
TRAIN: (419, 13) TEST: (104, 13)
TRAIN: (523, 13) TEST: (104, 13)<jupyter_text>Et on calé une forêt aléatoire...<jupyter_code>import warnings
from sklearn.ensemble import RandomForestRegressor
clr = RandomForestRegressor()

def train_test(clr, train_index, test_index):
    data_train = lagged.iloc[train_index, :]
    data_test = lagged.iloc[test_index, :]
    clr.fit(data_train.drop(["value", "date", "notrend", "trend"], 
                            axis=1), 
            data_train.value)
    r2 = r2_score(data_test.value,
                  clr.predict(data_test.drop(["value", "date", "notrend",
                                              "trend"], axis=1).as_matrix()))
    return r2

warnings.simplefilter("ignore")
last_test_index = None
for train_index, test_index in tscv.split(lagged):
    r2 = train_test(clr, train_index, test_index)        
    if last_test_index is not None:
        r2_prime = train_test(clr, last_test_index, test_index)        
        print(r2, r2_prime)
    else:
        print(r2)
    last_test_index = test_index<jupyter_output>0.6795844303954097
0.8643514472109561 0.8438674497139955
0.889424015517347 0.8566080912412136
0.8108804128996021 0.652299214153693
0.9082580468326908 0.9146950904017541<jupyter_text>2 ans coupé en 5, soit tous les 5 mois, ça veut dire que ce découpage inclut parfois Noël, parfois l'été et que les performances y seront très sensibles.<jupyter_code>from sklearn.metrics import r2_score
r2 = r2_score(data_test.value,
              clr.predict(data_test.drop(["value", "date", "notrend",
                                          "trend"], axis=1).as_matrix()))
r2<jupyter_output><empty_output><jupyter_text>On compare avec le $r_2$ avec le même $r_2$ obtenu en utilisant $Y_{t-1}$, $Y_{t-2}$, ... $Y_{t-d}$ comme prédiction.<jupyter_code>for i in range(1, 9):
    print(i, ":", r2_score(data_test.value, data_test["lag%d" % i]))
lagged[:5]<jupyter_output><empty_output><jupyter_text>En fait le jour de la semaine est une variable catégorielle, on crée une colonne par jour.<jupyter_code>from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
cols = ['lag1', 'lag2', 'lag3',
        'lag4', 'lag5', 'lag6', 'lag7', 'lag8']
ct = ColumnTransformer(
        [('pass', "passthrough", cols),
         ("dummies", OneHotEncoder(), ["weekday"])])
pred = ct.fit(lagged).transform(lagged[:5])
pred<jupyter_output><empty_output><jupyter_text>On met tout dans un pipeline parce que c'est plus joli, plus pratique aussi.<jupyter_code>from sklearn.pipeline import make_pipeline
from sklearn.decomposition import PCA, TruncatedSVD 
cols = ['lag1', 'lag2', 'lag3',
        'lag4', 'lag5', 'lag6', 'lag7', 'lag8']
model = make_pipeline(
         make_pipeline(
                 ColumnTransformer(
                    [('pass', "passthrough", cols),
                     ("dummies", make_pipeline(OneHotEncoder(), 
                                               TruncatedSVD(n_components=2)), ["weekday"])]),
        LinearRegression()))
model.fit(lagged, lagged["value"])<jupyter_output><empty_output><jupyter_text>C'est plus facile à voir visuellement.<jupyter_code>from mlinsights.plotting import pipeline2dot
dot = pipeline2dot(model, lagged)
from jyquickhelper import RenderJsDot
RenderJsDot(dot)
r2_score(lagged['value'], model.predict(lagged))<jupyter_output><empty_output><jupyter_text>TemplatingComplètement hors sujet mais utile.<jupyter_code>from jinja2 import Template
template = Template('Hello {{ name }}!')
template.render(name='<NAME>')
template = Template("""
{{ name }}
{{ "-" * len(name) }}
Possède :
{% for i in range(len(meubles)) %}
- {{meubles[i]}}{% endfor %}
""")
meubles = ['table', "tabouret"]
print(template.render(name='<NAME>', len=len,
                     meubles=meubles))<jupyter_output><NAME>
------------
Possède :

- table
- tabouret
StarcoderdataJupyterStructuredCleanDedup