.. _seance4projectionpopulationcorrectionrst:

========================================
Evoluation d’une population (correction)
========================================


.. only:: html

    **Links:** :download:`notebook <seance4_projection_population_correction.ipynb>`, :downloadlink:`html <seance4_projection_population_correction2html.html>`, :download:`PDF <seance4_projection_population_correction.pdf>`, :download:`python <seance4_projection_population_correction.py>`, :downloadlink:`slides <seance4_projection_population_correction.slides.html>`, :githublink:`GitHub|_doc/notebooks/sessions/seance4_projection_population_correction.ipynb|*`


Evolution d’une population à partir des tables de mortalités et d’une
situation initiale.

.. code:: ipython3

    %matplotlib inline
    import matplotlib.pyplot as plt

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





Exercice 1 : pyramides des âges
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    from actuariat_python.data import population_france_year
    population = population_france_year()
    df = population
    df.head(n=3)






.. 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>naissance</th>
          <th>age</th>
          <th>hommes</th>
          <th>femmes</th>
          <th>ensemble</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>2019</td>
          <td>0</td>
          <td>360058</td>
          <td>346324</td>
          <td>706382</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2018</td>
          <td>1</td>
          <td>365656</td>
          <td>350503</td>
          <td>716159</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2017</td>
          <td>2</td>
          <td>371835</td>
          <td>357304</td>
          <td>729139</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    hommes = df["hommes"]
    femmes = df["femmes"]
    somme = hommes - femmes

Je reprends ici le code exposé à `Damien Vergnaud’s
Homepage <http://www.di.ens.fr/~vergnaud/teaching1415/TPc_20150304_Python.html>`__
en l’adaptant un peu avec les fonctions de matplotlib via l’interface
`pyplot <http://matplotlib.org/api/pyplot_api.html>`__. Puis j’ajoute la
différence par âge. On commence souvent par la
`gallerie <http://matplotlib.org/gallery.html>`__ pour voir si un graphe
ou juste une partie est similaire à ce qu’on veut obtenir.

.. code:: ipython3

    from matplotlib import pyplot as plt
    from numpy import arange
    fig, ax = plt.subplots(figsize=(8,8))
    ValH = ax.barh(arange(len(hommes)), hommes, 1.0, label="Hommes",
                   color='b', linewidth=0, align='center')
    ValF = ax.barh(arange(len(femmes)), -femmes, 1.0, label="Femmes",
                   color='r', linewidth=0, align='center')
    diff, = ax.plot(somme, arange(len(femmes)), 'y', linewidth=2)
    ax.set_title("Pyramide des âges")
    ax.set_ylabel("Ages")
    ax.set_xlabel("Habitants")
    ax.set_ylim([0, 110])
    ax.legend((ValH[0], ValF[0], diff), ('Hommes', 'Femmes', 'différence'));



.. image:: seance4_projection_population_correction_7_0.png


Le même en utilisant la fonction insérée dans le module
`actuariat_python <http://www.xavierdupre.fr/app/actuariat_python/helpsphinx/index.html>`__.

.. code:: ipython3

    from actuariat_python.plots import plot_population_pyramid
    plot_population_pyramid(df["hommes"], df["femmes"], figsize=(8, 4));



.. image:: seance4_projection_population_correction_9_0.png


Exercice 2 : calcul de l’espérance de vie
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Le premier objectif est de calculer l’\ `espérance de
vie <https://en.wikipedia.org/wiki/Life_expectancy>`__ à l’âge :math:`t`
à partir de la table de mortalité. On récupère cette table.

.. code:: ipython3

    from actuariat_python.data import table_mortalite_france_00_02
    df=table_mortalite_france_00_02()
    import pandas
    pandas.concat([df.head(n=3), df.tail(n=3)])






.. 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>Age</th>
          <th>Homme</th>
          <th>Femme</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0</td>
          <td>100000.0</td>
          <td>100000.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>1</td>
          <td>99511.0</td>
          <td>99616.0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2</td>
          <td>99473.0</td>
          <td>99583.0</td>
        </tr>
        <tr>
          <th>115</th>
          <td>117</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>116</th>
          <td>118</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>117</th>
          <td>120</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
      </tbody>
    </table>
    </div>



On note :math:`P_t` la population l’âge :math:`t`. La probabilité de
mourir à la date :math:`t+d` lorsqu’on a l’âge :math:`t` correspond à la
probabilité de rester en vie à jusqu’à l’âge :math:`t+d` puis de mourir
dans l’année qui suit :

.. math:: m_{t+d} = \frac{P_{t+d}}{P_t}\frac{P_{t+d} - P_{t+d+1}}{P_{t+d}} 

.

L’espérance de vie s’exprime :

.. math:: \mathbb{E}(t) = \sum_{d=1}^\infty d m_{t+d} = \sum_{d=1}^\infty d \frac{P_{t+d}}{P_t}\frac{P_{t+d} - P_{t+d+1}}{P_{t+d}} = \sum_{d=1}^\infty d \frac{P_{t+d} - P_{t+d+1}}{P_{t}} 

On crée une matrice allant de 0 à 120 ans et on pose
:math:`\mathbb{E}(120)=0`. On utilise le module
`numpy <http://www.numpy.org/>`__.

.. code:: ipython3

    import numpy
    hf = df[["Homme", "Femme"]].values
    hf = numpy.vstack([hf, numpy.zeros((8, 2))])
    hf.shape




.. parsed-literal::
    (126, 2)



.. code:: ipython3

    nb = hf.shape[0]
    esp = numpy.zeros ((nb,2))
    for t in range(0,nb):
        for i in (0,1):
            if hf[t,i] == 0:
                esp[t,i] = 0
            else:
                somme  = 0.0
                for d in range(1,nb-t):
                    if hf[t+d,i] > 0:
                        somme += d * (hf[t+d,i] - hf[t+d+1,i]) / hf[t,i]
                esp[t,i] = somme
    esp[:1]            




.. parsed-literal::
    array([[75.00752, 82.48832]])



Enfin, on dessine le résultat avec matplotlib :

.. code:: ipython3

    h = plt.plot(esp)
    plt.legend(h, ["Homme", "Femme"])
    plt.title("Espérance de vie");



.. image:: seance4_projection_population_correction_17_0.png


Le calcul implémenté ci-dessus n’est pas le plus efficace. On fait deux
boucles imbriquées dont le coût global est en :math:`O(n^2)` mais
surtout on effectue les mêmes calculs plusieurs fois. Pour le réduire à
un coût linéaire :math:`O(n)`, il faut s’intéresser à la quantité :

.. math:: P_{t+1} \mathbb{E}(t+1) - P_t \mathbb{E}(t) = \sum_{d=1}^\infty d (P_{t+d+1} - P_{t+d+2}) - \sum_{d=1}^\infty d (P_{t+d} - P_{t+d+1})

L’implémentation devra utiliser la fonction
`numpy.cumsum <http://docs.scipy.org/doc/numpy/reference/generated/numpy.cumsum.html>`__
et cette astuce `Pandas Dataframe cumsum by row in reverse column
order? <http://stackoverflow.com/questions/26389959/pandas-dataframe-cumsum-by-row-in-reverse-column-order>`__.

.. code:: ipython3

    # à suivre

`numpy <http://www.numpy.org/>`__ et
`pandas <http://pandas.pydata.org/>`__ ont plusieurs fonction en commun
dès qu’il s’agit de parcourir les données. Il existe aussi la fonction
`DataFrame.cumsum <http://pandas.pydata.org/pandas-docs/dev/generated/pandas.DataFrame.cumsum.html>`__.

Exercice 3 : simulation de la pyramide l’année suivante
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

L’objectif est d’estimer la population française à l’année :math:`T`. Si
:math:`P(a,T)` désigne le nombre de personnes d’âge :math:`a` en
:math:`T`, on peut estimer :math:`P(a,T+1)` en utilisant la probabilité
de mourir :math:`m(a)` :

.. math::  P(a+1, T+1) = P(a,T) * (1 - m(a))

On commence par calculer les coefficients :math:`m(a)` avec la table
``hf`` obtenue lors de l’exercice précédent tout en gardant la même
dimension (on aura besoin de la fonction
`nan_to_num <http://docs.scipy.org/doc/numpy/reference/generated/numpy.nan_to_num.html>`__
:

.. code:: ipython3

    mortalite = (hf[:-1] - hf[1:]) / hf[:-1]
    mortalite = numpy.nan_to_num(mortalite)  # les divisions nulles deviennent nan, on les remplace par 0
    mortalite = numpy.vstack([mortalite, numpy.zeros((1, 2))])
    m = mortalite


.. parsed-literal::
    c:\python372_x64\lib\site-packages\ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in true_divide
      """Entry point for launching an IPython kernel.


La population a été obtenue lors de l’exercice 1, on la convertit en un
objet `numpy <http://www.numpy.org/>`__ :

.. code:: ipython3

    pop = population[["hommes","femmes"]].values
    pop = numpy.vstack( [pop, numpy.zeros((m.shape[0] - pop.shape[0], 2))])
    pop0 = pop.copy()
    pop.shape




.. parsed-literal::
    (126, 2)



Ensuite on calcule la population en 2020 :

.. code:: ipython3

    pop_next = pop * (1-m)
    pop_next = numpy.vstack([numpy.zeros((1, 2)), pop_next[:-1]])
    pop_next[:5]




.. parsed-literal::
    array([[     0.        ,      0.        ],
           [358297.31638   , 344994.11584   ],
           [365516.36791912, 350386.88814046],
           [371734.07266293, 357228.65195867],
           [382450.37346902, 366544.40263353]])



.. code:: ipython3

    pop[:5]




.. parsed-literal::
    array([[360058., 346324.],
           [365656., 350503.],
           [371835., 357304.],
           [382535., 366607.],
           [393693., 377204.]])



.. code:: ipython3

    from actuariat_python.plots import plot_population_pyramid
    plot_population_pyramid(pop_next[:, 0], pop_next[:, 1]);



.. image:: seance4_projection_population_correction_29_0.png


Exercice 4 : simulation jusqu’en 2100
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Il s’agit de répéter l’itération effectuée lors de l’exercice précédent.
Le plus est de recopier le code dans une fonction et de l’appeler un
grand nombre de fois.

.. code:: ipython3

    def iteration(pop, mortalite):
        pop_next = pop * (1-mortalite)
        pop_next = numpy.vstack([numpy.zeros((1, 2)), 
                                 pop_next[:-1]]) # aucune naissance
        return pop_next
    
    popt = pop
    for year in range(2020, 2051):
        popt = iteration(popt, mortalite)
                                 
    plot_population_pyramid(popt[:,0], popt[:,1], title="Pyramide des âges en 2050");



.. image:: seance4_projection_population_correction_31_0.png


Exercice 5 : simulation avec les naissances
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Dans l’exercice précédent, la seconde ligne de la fonction *iteration*
correspond à cas où il n’y a pas de naissance. On veut remplacer cette
ligne par quelque chose proche de la réalité :


-  les naissances sont calculées à partir de la population féminines et
   de la table de fécondité
-  on garde la même proportion homme/femme que celle actuellement
   observée

.. code:: ipython3

    ratio = pop[0, 0] / (pop[0, 1] + pop[0, 0])
    ratio




.. parsed-literal::
    0.5097213688910532



Il y a un peu plus de garçons qui naissent chaque année.

.. code:: ipython3

    from actuariat_python.data import fecondite_france
    df = fecondite_france()
    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>age</th>
          <th>2005</th>
          <th>2015</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>4</th>
          <td>16.0</td>
          <td>0.2</td>
          <td>0.2</td>
        </tr>
        <tr>
          <th>5</th>
          <td>17.0</td>
          <td>0.5</td>
          <td>0.5</td>
        </tr>
        <tr>
          <th>6</th>
          <td>18.0</td>
          <td>1.0</td>
          <td>1.1</td>
        </tr>
        <tr>
          <th>7</th>
          <td>19.0</td>
          <td>2.0</td>
          <td>2.1</td>
        </tr>
        <tr>
          <th>8</th>
          <td>20.0</td>
          <td>3.0</td>
          <td>3.2</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    from matplotlib import pyplot as plt
    df.plot(x="age", y=["2005", "2015"]);



.. image:: seance4_projection_population_correction_36_0.png


On convertit ces données en une matrice numpy sur 120 lignes comme les
précédentes. On se sert des méthodes
`fillna <http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.fillna.html>`__
et
`merge <http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.merge.html?highlight=merge#pandas.DataFrame.merge>`__.

.. code:: ipython3

    ages = pandas.DataFrame(dict(age=range(0,120)))
    merge = ages.merge(df, left_on="age", right_on="age", how="outer")
    fecondite = merge.fillna(0.0)
    fecondite[13:17]






.. 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>age</th>
          <th>2005</th>
          <th>2015</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>13</th>
          <td>13</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>14</th>
          <td>14</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>15</th>
          <td>15</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>16</th>
          <td>16</td>
          <td>0.2</td>
          <td>0.2</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    mat_fec = fecondite[["2015"]].values / 100  # les chiffres sont pour 100 femmes
    mat_fec.shape




.. parsed-literal::
    (120, 1)



.. code:: ipython3

    mat_fec.sum()




.. parsed-literal::
    1.899



Si la matrice ``pop`` a plus de ligne que la matrice ``mat_fec``, on
doit compléter la seconde avec autant de lignes nulles que la
précédente.

.. code:: ipython3

    if mat_fec.shape[0] < pop.shape[0]:
        zeros = numpy.zeros((pop.shape[0] - mat_fec.shape[0], mat_fec.shape[1]))
        mat_fec = numpy.vstack([mat_fec, zeros])
    mat_fec.sum()




.. parsed-literal::
    1.899



Il faut maintenant coder une fonction qui calcule le naissances pour
l’année suivantes.

.. code:: ipython3

    def naissances(pop, fec):
        # on suppose que pop est une matrice avec deux colonnes homme, femme
        # et que fec est une matrice avec une colonne fécondité
        n = pop[:, 1] * fec[:, 0]
        return n.sum()
    
    nais = naissances(pop, mat_fec)
    nais




.. parsed-literal::
    756839.2760000001



Et on reprend la fonction *iteration* et le code de l’exercice précédent
:

.. code:: ipython3

    def iteration(pop, mortalite, fec, ratio):
        pop_next = pop * (1 - mortalite)
        nais = naissances(pop, fec)
        row = numpy.array([[nais * ratio, nais * (1 - ratio)]])
        pop_next = numpy.vstack([row, pop_next[:-1]]) # aucune naissance
        return pop_next
    
    popt = pop
    for year in range(2020, 2101):
        popt = iteration(popt, m, mat_fec, ratio)
        if year == 2050:
            popt2050 = popt.copy()
        if year == 2100:
            popt2100 = popt.copy()
                                 
    fig, ax = plt.subplots(1, 2, figsize=(10, 4))
    plot_population_pyramid(popt2050[:, 0], popt2050[:, 1], ax=ax[0],
                            title="Pyramide des âges en 2050")
    plot_population_pyramid(popt2100[:, 0], popt2100[:, 1], ax=ax[1],
                            title="Pyramide des âges en 2100");



.. image:: seance4_projection_population_correction_46_0.png


On va plus loin et on stocke la population dans un vecteur :

.. code:: ipython3

    total = [[2020, pop[:,0].sum(),pop[:,1].sum()]]
    popt = pop
    for year in range(2020, 2101):
        popt = iteration(popt, m, mat_fec, ratio)
        total.append([year, popt[:,0].sum(),popt[:,1].sum()])
                                 
    plot_population_pyramid(popt[:, 0], popt[:, 1], title="Pyramide des âges en 2101");



.. image:: seance4_projection_population_correction_48_0.png


.. code:: ipython3

    df = pandas.DataFrame(data=total, columns=["année","hommes","femmes"])
    df.plot(x="année", y=["hommes", "femmes"], title="projection population française");



.. image:: seance4_projection_population_correction_49_0.png


Le code suivant permet de combiner les deux graphes sur la même ligne
avec la fonction
`subplots <http://matplotlib.org/api/pyplot_api.html?highlight=subplots#matplotlib.pyplot.subplots>`__
:

.. code:: ipython3

    from matplotlib import pyplot as plt
    fig, ax = plt.subplots(1, 2, figsize=(14,6))
    plot_population_pyramid(popt[:,0], popt[:,1], title="Pyramide des âges en 2050", ax=ax[0])
    df.plot(x="année", y=["hommes", "femmes"], title="projection population française", ax=ax[1]);



.. image:: seance4_projection_population_correction_51_0.png


Retraites
---------

La réforme de la retraite en 2019 laissera un souvenir indélébile. Une
des mesures far est celle de l’âge pivot. Le modèle développé jusqu’à
présent ne permet pas de simuler ce que chaque français pourrait obtenir
avec la nouvelle réforme mais il permet de donner un ordre de grandeur
sur le nombre de retraités selon que l’âge pivot est de 62 ans ou plus
vieux. Si on en croit cette page `1,7 actif cotisant par
retraité <https://reforme-retraite.gouv.fr/le-saviez-vous/les-chiffres-cles/article/1-7-actif-cotisant-par-retraite>`__,
il y a 1,7 cotisant par retraité.

Trois populations sont considérées : les jeunes de moins de 21 ans, à la
charge de leurs aînés, les actifs, les retraités dont le nombre dépend
de l’âge pivot. on veut donc calculer l’évolution de ces trois
populations.

.. code:: ipython3

    from tqdm import tqdm
    
    evol = []
    popt = pop0.copy()
    age_etude = 23
    pivot = list(range(60, 68))
    for year in tqdm(range(2020, 2101)):
        popt = iteration(popt, m, mat_fec, ratio)
        jeune = popt[:age_etude + 1].sum()
        row = dict(year=year)
        for p in pivot:
            actif = popt[age_etude + 1:p].sum()
            retraite = popt[p:].sum()
            rt = actif / retraite
            row['r%d' % p] = rt
        evol.append(row)
    
    df_evol = pandas.DataFrame(evol)
    df_evol.head()


.. parsed-literal::
    100%|████████████████████████████████████████████████████████████████████████████████| 81/81 [00:00<00:00, 6769.86it/s]






.. 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>year</th>
          <th>r60</th>
          <th>r61</th>
          <th>r62</th>
          <th>r63</th>
          <th>r64</th>
          <th>r65</th>
          <th>r66</th>
          <th>r67</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>2020</td>
          <td>1.661675</td>
          <td>1.792986</td>
          <td>1.936632</td>
          <td>2.092034</td>
          <td>2.262744</td>
          <td>2.450994</td>
          <td>2.659420</td>
          <td>2.893034</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2021</td>
          <td>1.637597</td>
          <td>1.766554</td>
          <td>1.907657</td>
          <td>2.062432</td>
          <td>2.230352</td>
          <td>2.415367</td>
          <td>2.620045</td>
          <td>2.847438</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2022</td>
          <td>1.616615</td>
          <td>1.742623</td>
          <td>1.881228</td>
          <td>2.033293</td>
          <td>2.200564</td>
          <td>2.382585</td>
          <td>2.583765</td>
          <td>2.807077</td>
        </tr>
        <tr>
          <th>3</th>
          <td>2023</td>
          <td>1.593194</td>
          <td>1.720874</td>
          <td>1.856290</td>
          <td>2.005630</td>
          <td>2.169929</td>
          <td>2.351192</td>
          <td>2.549055</td>
          <td>2.768466</td>
        </tr>
        <tr>
          <th>4</th>
          <td>2024</td>
          <td>1.572343</td>
          <td>1.698854</td>
          <td>1.836163</td>
          <td>1.982169</td>
          <td>2.143622</td>
          <td>2.321765</td>
          <td>2.518906</td>
          <td>2.734809</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    ax = df_evol.plot(x="year", y=["r%d" % p for p in pivot], figsize=(10, 4))
    years = df_evol.year
    ax.plot(years, [1.7 for _ in years], 'b--')
    ax.set_title("Ratio actifs / retraités en fonction de l'âge pivot");



.. image:: seance4_projection_population_correction_54_0.png


D’après ces simulations, il faudrait reculer l’âge de la retraite de
deux ou trois ans pour garder le même ratio entre actifs et retraités
tel qu’il est aujourd’hui.