.. _timeseriesssarst:

==============================
Single Spectrum Analysis (SSA)
==============================


.. only:: html

    **Links:** :download:`notebook <timeseries_ssa.ipynb>`, :downloadlink:`html <timeseries_ssa2html.html>`, :download:`python <timeseries_ssa.py>`, :downloadlink:`slides <timeseries_ssa.slides.html>`, :githublink:`GitHub|_doc/notebooks/2a/timeseries_ssa.ipynb|*`


Illustration de la méthode
`SSA <https://en.wikipedia.org/wiki/Singular_spectrum_analysis>`__ pour
les séries temporelles appliquée à la détection de points aberrants. La
méthode est décrite dans `Singular Spectrum Analysis: Methodology and
Comparison <http://perso.ens-lyon.fr/patrick.flandrin/hassani.pdf>`__.
Voir aussi `Automated outlier detection in Singular Spectrum
Analysis <https://www.statistics.gov.hk/wsc/CPS204-P35-S.pdf>`__.

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





.. code:: ipython3

    %matplotlib inline

Une série articielle
--------------------

On introduit quelques points aberrants, pour le reste, elle suit le
modèle :math:`y_t = \frac{9}{10} y_{t-2} + \epsilon_t + a_t` où
:math:`a_t` est le bruit aberrant qui survient quelques fois.

.. code:: ipython3

    import numpy.random as rnd
    import numpy
    N = 2000
    bruit1 = rnd.normal(size=(N,))
    temps = numpy.arange(N)
    bruit1[:5], temps[:5]




.. parsed-literal::
    (array([-0.42127366, -0.80614375,  1.87947411,  0.13890672,  0.55054172]),
     array([0, 1, 2, 3, 4]))



On crée un bruit aberrant.

.. code:: ipython3

    import random
    bruit2 = numpy.zeros((N,))
    for i in range(0, 10):
        h = random.randint(0, N-1)
        bruit2[h] = rnd.normal() + 10

.. code:: ipython3

    serie = []
    y = 10
    for i in range(N//2+100):
        serie.append(y + bruit1[i] + 0.0004 * temps[i] + bruit2[i])
        if i > 30:
            y = 0.9 * serie[-2]
    Y = numpy.array(serie[-1000:])
    Y[:5]




.. parsed-literal::
    array([ 0.35328897,  0.43350731, -0.25385586,  1.62514162,  0.09987641])



.. code:: ipython3

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(1, 2, figsize=(12, 4))
    ax[0].plot(numpy.arange(len(Y)), Y)
    ax[1].plot(numpy.arange(800, 900), Y[800:900])
    ax[0].set_title("Série temporelle simulée")
    ax[1].set_title("Même série temporelle simulée");



.. image:: timeseries_ssa_8_0.png


Autocorrélations
----------------

L’autocorrélogramme est définie par la série
:math:`cor(Y_t, Y_{t-d})_d`. On le calcule sur la série nettoyée de sa
tendance.

.. code:: ipython3

    from pandas import DataFrame
    df = DataFrame(dict(Y=Y))
    df.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Y</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.353289</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.433507</td>
        </tr>
        <tr>
          <th>2</th>
          <td>-0.253856</td>
        </tr>
        <tr>
          <th>3</th>
          <td>1.625142</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.099876</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    from statsmodels.tsa.tsatools import detrend
    df["notrend"] = detrend(df.Y)
    df.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Y</th>
          <th>notrend</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.353289</td>
          <td>0.216265</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.433507</td>
          <td>0.292403</td>
        </tr>
        <tr>
          <th>2</th>
          <td>-0.253856</td>
          <td>-0.399040</td>
        </tr>
        <tr>
          <th>3</th>
          <td>1.625142</td>
          <td>1.475877</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.099876</td>
          <td>-0.053468</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    ax = df.plot()
    ax.set_title("Sans tendance");



.. image:: timeseries_ssa_12_0.png


L’autocorrélogramme à proprement parler.

.. code:: ipython3

    from statsmodels.tsa.stattools import acf
    cor = acf(df.notrend)
    fig, ax = plt.subplots(1, 1)
    ax.plot(cor)
    ax.set_title("Autocorrélogramme");



.. image:: timeseries_ssa_14_0.png


Etant donnée que la série :math:`Y_t` dépend de :math:`Y_{t-2}`, on
observe un pic pour :math:`cor(Y_t, Y_{t-2})_d` et pour tous les
:math:`d` pairs. :math:`cor(Y_t, Y_{t-4}) \sim cor(Y_t, Y_{t-2})^2`. On
enlève ces effets récursifs en calculant l’autocorrélogramme partiel qui
correspond à l’estimation des coefficients d’un `modèle
autorégressif <https://fr.wikipedia.org/wiki/Processus_autor%C3%A9gressif>`__
infini.

.. code:: ipython3

    from statsmodels.tsa.stattools import pacf
    pcor = pacf(df.notrend)
    fig, ax = plt.subplots(1, 1)
    ax.plot(pcor[1:])
    ax.set_title("Autocorrélogramme partiel");



.. image:: timeseries_ssa_16_0.png


SSA
---

Ou Singular Spectrum Analysis. La méthode part de la matrice des séries
décalées qu’on décompose avec la méthode SVD ou `Singular Value
Decomposition <https://fr.wikipedia.org/wiki/D%C3%A9composition_en_valeurs_singuli%C3%A8res>`__.

.. code:: ipython3

    def lagged_ts(serie, lag):
        dim = serie.shape[0]
        res = numpy.zeros((dim - lag + 1, lag))
        for i in range(lag):
            res[:, i] = serie[i:dim-lag+i+1]
        return res
    
    lagged_ts(Y, 3)




.. parsed-literal::
    array([[ 0.35328897,  0.43350731, -0.25385586],
           [ 0.43350731, -0.25385586,  1.62514162],
           [-0.25385586,  1.62514162,  0.09987641],
           ...,
           [ 5.59900199,  4.51154279,  6.48711916],
           [ 4.51154279,  6.48711916,  4.94491822],
           [ 6.48711916,  4.94491822,  6.90403987]])



.. code:: ipython3

    lag = lagged_ts(Y, 60)
    lag.shape




.. parsed-literal::
    (941, 60)



.. code:: ipython3

    from numpy.linalg import svd
    # u @ numpy.diag(s) @ vh
    u, s, vh = svd(lag)

.. code:: ipython3

    u.shape, s.shape, vh.shape




.. parsed-literal::
    ((941, 941), (60,), (60, 60))



.. code:: ipython3

    d = numpy.zeros((941, 60))
    d[:60,:60] = numpy.diag(s)

.. code:: ipython3

    (u @ d @ vh).shape




.. parsed-literal::
    (941, 60)



.. code:: ipython3

    fig, ax = plt.subplots(1,2, figsize=(12,4))
    DataFrame(dict(valeur=s)).plot(kind="bar", ax=ax[0])
    DataFrame(dict(valeur=s[1:15])).plot(kind="bar", ax=ax[1])
    ax[0].set_title("Valeurs propres")
    ax[1].set_title("Valeurs propres sans la première");



.. image:: timeseries_ssa_24_0.png


Je me représente la méthode SVD comme une façon de projeter des vecteurs
sur l’espace vectoriel constitué des premiers vecteurs propres, à chaque
dimension supplémentaire, c’est comme une pièce du puzzle qui s’assemble
jusqu’à recomposer l’ensemble. Ce qu’on peut voir aussi comme ceci :

.. code:: ipython3

    np = 12
    fig, ax = plt.subplots(np, 3, figsize=(14, np*2))
    for n in range(np):
        i = n if n < 5 else n * 5 - 15
        d = numpy.zeros((941, 60))
        d[i, i] = s[i]
        X2 = u @ d @ vh
        pos = 0 #X2.shape[1] - 1
        
        # série reconstruites avec un axe
        ax[n, 0].plot(X2[:,pos])
        ax[n, 1].set_title("i=%d" % i)    
        # série reconstruites avec un axe
        ax[n, 1].plot(X2[800:850,pos])
        ax[n, 1].set_title("i=%d" % i)
        
        d = numpy.zeros((941, 60))
        d[:i+1, :i+1] = numpy.diag(s[:i+1])
        X2 = u @ d @ vh
        ax[n, 2].plot(X2[800:850,pos])    
        ax[n, 2].plot(Y[800:850])
        ax[n, 2].set_title("-->i=%d + 1" % i)



.. image:: timeseries_ssa_26_0.png


La prédiction
-------------

On veut prédire :math:`Y_{t+1}`. L’idée consiste à appliquer la méthode
en considérent :math:`Y_{t+1}` égale à :math:`Y_t` puis à remplacer
cette prédiction par la valeur de la série reconstruite. On peut même
prédire à un horizon plus grand que la valeur suivante.

.. code:: ipython3

    fig, ax = plt.subplots(1, 1, figsize=(5,5))
    for i in range(0, 8):
        ax.plot([0, 5], [i, i], 'k-')
        if i < 6:
            ax.plot([i, i], [0, 7], 'k-')
        if i < 4:
            ax.text(i + 0.1, 1.5, "Y(t-%d)" % (4-i))
            ax.text(i + 0.1, 0.5, "Y(t-%d)" % (3-i))
    ax.text(4.1, 1.5, "Y(t)")
    ax.text(4.05, 0.5, "Y(t+1)=?")
    plt.axis('off');



.. image:: timeseries_ssa_28_0.png


Les points aberrants
--------------------

On repère les points aberrants avec l’une méthode de son choix sur la
série reconstruite.

.. code:: ipython3

    d = numpy.zeros((941, 60))
    for i in range(0, 30):
        d[i, i] = s[i]
    X2 = u @ d @ vh

.. code:: ipython3

    from sklearn.covariance import EllipticEnvelope
    env = EllipticEnvelope(support_fraction=0.9)
    env.fit(X2[:,:30])




.. parsed-literal::
    EllipticEnvelope(assume_centered=False, contamination=0.1, random_state=None,
             store_precision=True, support_fraction=0.9)



L’idéal serait d’utiliser une méthode basée sur une ACP. Le plus proche
reste le modèle gaussien avec
`EllipticEnvelope <http://scikit-learn.org/dev/modules/generated/sklearn.covariance.EllipticEnvelope.html#sklearn.covariance.EllipticEnvelope>`__.

.. code:: ipython3

    out = env.predict(X2[:,:30])
    
    fig, ax = plt.subplots(1, 1, figsize=(12,2))
    ax.plot((1 - out)/2, "-")
    ax.set_title("Points aberrants d'une série temporelles.");



.. image:: timeseries_ssa_33_0.png


On observe des plages contiguës. Cela signifie que d’une valeur
aberrante contamine des vecteurs décalées consécutifs de la série
:math:`Y`. Il ne reste plus qu’à repérer la valeur incriminée.

.. code:: ipython3

    fig, ax = plt.subplots(1, 1, figsize=(14,2))
    ax.plot(X2[:,0], label="serie")
    ax.plot((1 - out)*5, "-", label="outlier")
    ax.set_title("Points aberrants sur la série reconstruite")
    ax.legend();



.. image:: timeseries_ssa_35_0.png


.. code:: ipython3

    fig, ax = plt.subplots(1, 1, figsize=(14,2))
    ax.plot(Y, label="serie")
    ax.plot((1 - out)*5, "-", label="outlier")
    ax.set_title("Points aberrants sur la série brute")
    ax.legend();



.. image:: timeseries_ssa_36_0.png


Ce qui a l’air de correspondre à la fin des grandes plages. On
recommence avec la probabilité d’être un outlier.

.. code:: ipython3

    fig, ax = plt.subplots(1, 1, figsize=(14,2))
    outp = env.decision_function(X2[:,:30])
    ax.plot(Y, label="serie")
    ax.plot(outp, "-", label="Proba not outlier")
    ax.set_title("Points aberrants sur la série brute")
    ax.legend();



.. image:: timeseries_ssa_38_0.png