.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/indicators/plot_b09wtp.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_indicators_plot_b09wtp.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_indicators_plot_b09wtp.py:


Calculation of willingness to pay
=================================

We calculate and plot willingness to pay.
Details about this example are available in Section 4 of `Bierlaire (2018)
Calculating indicators with PandasBiogeme
<http://transp-or.epfl.ch/documents/technicalReports/Bier18a.pdf>`_

:author: Michel Bierlaire, EPFL
:date: Wed Apr 12 20:57:00 2023

.. GENERATED FROM PYTHON SOURCE LINES 16-33

.. code-block:: default


    import sys

    try:
        import matplotlib.pyplot as plt

        can_plot = True
    except ModuleNotFoundError:
        can_plot = False
    import biogeme.biogeme as bio
    import biogeme.exceptions as excep
    import biogeme.results as res

    from biogeme.expressions import Derive
    from optima_data import database, normalized_weight
    from scenarios import scenario








.. GENERATED FROM PYTHON SOURCE LINES 34-36

Obtain the specification for the default scenario
The definition of the scenarios is available in :ref:`scenarios`.

.. GENERATED FROM PYTHON SOURCE LINES 36-41

.. code-block:: default

    V, _, _, _ = scenario()

    V_PT = V[0]
    V_CAR = V[1]








.. GENERATED FROM PYTHON SOURCE LINES 42-43

Calculation of the willingness to pay using derivatives.

.. GENERATED FROM PYTHON SOURCE LINES 43-46

.. code-block:: default

    WTP_PT_TIME = Derive(V_PT, 'TimePT') / Derive(V_PT, 'MarginalCostPT')
    WTP_CAR_TIME = Derive(V_CAR, 'TimeCar') / Derive(V_CAR, 'CostCarCHF')








.. GENERATED FROM PYTHON SOURCE LINES 47-48

Formulas to simulate.

.. GENERATED FROM PYTHON SOURCE LINES 48-54

.. code-block:: default

    simulate = {
        'weight': normalized_weight,
        'WTP PT time': WTP_PT_TIME,
        'WTP CAR time': WTP_CAR_TIME,
    }








.. GENERATED FROM PYTHON SOURCE LINES 55-56

Create the Biogeme object.

.. GENERATED FROM PYTHON SOURCE LINES 56-58

.. code-block:: default

    the_biogeme = bio.BIOGEME(database, simulate)








.. GENERATED FROM PYTHON SOURCE LINES 59-60

Read the estimation results from the file.

.. GENERATED FROM PYTHON SOURCE LINES 60-68

.. code-block:: default

    try:
        results = res.bioResults(pickleFile='saved_results/b02estimation.pickle')
    except excep.BiogemeError:
        sys.exit(
            'Run first the script b02estimation.py in order to generate '
            'the file b02estimation.pickle.'
        )








.. GENERATED FROM PYTHON SOURCE LINES 69-71

`simulated_values` is a Pandas dataframe with the same number of rows as
the database, and as many columns as formulas to simulate.

.. GENERATED FROM PYTHON SOURCE LINES 71-74

.. code-block:: default

    simulated_values = the_biogeme.simulate(results.getBetaValues())
    simulated_values






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <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>weight</th>
          <th>WTP PT time</th>
          <th>WTP CAR time</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.886023</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.861136</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.861136</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.957386</td>
          <td>0.111499</td>
          <td>0.111499</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.861136</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>2259</th>
          <td>2.036009</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
        <tr>
          <th>2261</th>
          <td>0.861136</td>
          <td>0.111499</td>
          <td>0.111499</td>
        </tr>
        <tr>
          <th>2262</th>
          <td>0.861136</td>
          <td>0.111499</td>
          <td>0.111499</td>
        </tr>
        <tr>
          <th>2263</th>
          <td>0.957386</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
        <tr>
          <th>2264</th>
          <td>0.957386</td>
          <td>0.040308</td>
          <td>0.040308</td>
        </tr>
      </tbody>
    </table>
    <p>1906 rows × 3 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 75-76

Note the multiplication by 60 to have the valus of time per hour and not per minute.

.. GENERATED FROM PYTHON SOURCE LINES 76-78

.. code-block:: default

    wtpcar = (60 * simulated_values['WTP CAR time'] * simulated_values['weight']).mean()








.. GENERATED FROM PYTHON SOURCE LINES 79-80

Calculate confidence intervals

.. GENERATED FROM PYTHON SOURCE LINES 80-82

.. code-block:: default

    b = results.getBetasForSensitivityAnalysis(the_biogeme.free_beta_names())








.. GENERATED FROM PYTHON SOURCE LINES 83-85

Returns data frame containing, for each simulated value, the left
and right bounds of the confidence interval calculated by simulation.

.. GENERATED FROM PYTHON SOURCE LINES 85-87

.. code-block:: default

    left, right = the_biogeme.confidenceIntervals(b, 0.9)








.. GENERATED FROM PYTHON SOURCE LINES 88-89

Lower bounds of the confidence intervals

.. GENERATED FROM PYTHON SOURCE LINES 89-90

.. code-block:: default

    left





.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <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>weight</th>
          <th>WTP PT time</th>
          <th>WTP CAR time</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.886023</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.861136</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.861136</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.957386</td>
          <td>0.069673</td>
          <td>0.069673</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.861136</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>2259</th>
          <td>2.036009</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
        <tr>
          <th>2261</th>
          <td>0.861136</td>
          <td>0.069673</td>
          <td>0.069673</td>
        </tr>
        <tr>
          <th>2262</th>
          <td>0.861136</td>
          <td>0.069673</td>
          <td>0.069673</td>
        </tr>
        <tr>
          <th>2263</th>
          <td>0.957386</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
        <tr>
          <th>2264</th>
          <td>0.957386</td>
          <td>0.011110</td>
          <td>0.011110</td>
        </tr>
      </tbody>
    </table>
    <p>1906 rows × 3 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 91-92

Upper bounds of the confidence intervals

.. GENERATED FROM PYTHON SOURCE LINES 92-94

.. code-block:: default

    right






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <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>weight</th>
          <th>WTP PT time</th>
          <th>WTP CAR time</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.886023</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.861136</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.861136</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.957386</td>
          <td>0.17001</td>
          <td>0.17001</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.861136</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>2259</th>
          <td>2.036009</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
        <tr>
          <th>2261</th>
          <td>0.861136</td>
          <td>0.17001</td>
          <td>0.17001</td>
        </tr>
        <tr>
          <th>2262</th>
          <td>0.861136</td>
          <td>0.17001</td>
          <td>0.17001</td>
        </tr>
        <tr>
          <th>2263</th>
          <td>0.957386</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
        <tr>
          <th>2264</th>
          <td>0.957386</td>
          <td>0.08798</td>
          <td>0.08798</td>
        </tr>
      </tbody>
    </table>
    <p>1906 rows × 3 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 95-96

Lower and upper bounds of the willingness to pay.

.. GENERATED FROM PYTHON SOURCE LINES 96-102

.. code-block:: default

    wtpcar_left = (60 * left['WTP CAR time'] * left['weight']).mean()
    wtpcar_right = (60 * right['WTP CAR time'] * right['weight']).mean()
    print(
        f'Average WTP for car: {wtpcar:.3g} ' f'CI:[{wtpcar_left:.3g}, {wtpcar_right:.3g}]'
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Average WTP for car: 3.96 CI:[1.93, 7.05]




.. GENERATED FROM PYTHON SOURCE LINES 103-105

In this specific case, there are only two distinct values in the
population: for workers and non workers

.. GENERATED FROM PYTHON SOURCE LINES 105-111

.. code-block:: default

    print(
        'Unique values:      ',
        [f'{i:.3g}' for i in 60 * simulated_values['WTP CAR time'].unique()],
    )






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Unique values:       ['2.42', '6.69']




.. GENERATED FROM PYTHON SOURCE LINES 112-113

Function calculating the willingness to pay for a group.

.. GENERATED FROM PYTHON SOURCE LINES 113-132

.. code-block:: default

    def wtp_for_subgroup(the_filter: 'pd.Series[np.bool_]') -> tuple[float, float, float]:
        """
        Check the value for groups of the population. Define a function that
        works for any filter to avoid repeating code.

        :param the_filter: pandas filter

        :return: willingness-to-pay for car and confidence interval
        """
        size = the_filter.sum()
        sim = simulated_values[the_filter]
        total_weight = sim['weight'].sum()
        weight = sim['weight'] * size / total_weight
        _wtpcar = (60 * sim['WTP CAR time'] * weight).mean()
        _wtpcar_left = (60 * left[the_filter]['WTP CAR time'] * weight).mean()
        _wtpcar_right = (60 * right[the_filter]['WTP CAR time'] * weight).mean()
        return _wtpcar, _wtpcar_left, _wtpcar_right









.. GENERATED FROM PYTHON SOURCE LINES 133-134

Full time workers.

.. GENERATED FROM PYTHON SOURCE LINES 134-138

.. code-block:: default

    aFilter = database.data['OccupStat'] == 1
    w, l, r = wtp_for_subgroup(aFilter)
    print(f'WTP car for workers: {w:.3g} CI:[{l:.3g}, {r:.3g}]')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    WTP car for workers: 6.69 CI:[4.18, 10.2]




.. GENERATED FROM PYTHON SOURCE LINES 139-140

Females.

.. GENERATED FROM PYTHON SOURCE LINES 140-144

.. code-block:: default

    aFilter = database.data['Gender'] == 2
    w, l, r = wtp_for_subgroup(aFilter)
    print(f'WTP car for females: {w:.3g} CI:[{l:.3g}, {r:.3g}]')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    WTP car for females: 3.17 CI:[1.29, 6.15]




.. GENERATED FROM PYTHON SOURCE LINES 145-146

Males.

.. GENERATED FROM PYTHON SOURCE LINES 146-150

.. code-block:: default

    aFilter = database.data['Gender'] == 1
    w, l, r = wtp_for_subgroup(aFilter)
    print(f'WTP car for males  : {w:.3g} CI:[{l:.3g}, {r:.3g}]')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    WTP car for males  : 4.96 CI:[2.75, 8.2]




.. GENERATED FROM PYTHON SOURCE LINES 151-153

We plot the distribution of WTP in the population. In this case,
there are only two values

.. GENERATED FROM PYTHON SOURCE LINES 153-161

.. code-block:: default

    if can_plot:
        plt.hist(
            60 * simulated_values['WTP CAR time'],
            weights=simulated_values['weight'],
        )
        plt.xlabel('WTP (CHF/hour)')
        plt.ylabel('Individuals')
        plt.show()



.. image-sg:: /auto_examples/indicators/images/sphx_glr_plot_b09wtp_001.png
   :alt: plot b09wtp
   :srcset: /auto_examples/indicators/images/sphx_glr_plot_b09wtp_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.353 seconds)


.. _sphx_glr_download_auto_examples_indicators_plot_b09wtp.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_b09wtp.py <plot_b09wtp.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_b09wtp.ipynb <plot_b09wtp.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_