Note
Go to the end to download the full example code.
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
- author:
Michel Bierlaire, EPFL
- date:
Wed Apr 12 20:57:00 2023
import sys
import numpy as np
import pandas as pd
from IPython.core.display_functions import display
try:
import matplotlib.pyplot as plt
can_plot = True
except ModuleNotFoundError:
can_plot = False
import biogeme.biogeme as bio
from biogeme.exceptions import BiogemeError
import biogeme.results as res
from biogeme.expressions import Derive
from biogeme.data.optima import read_data, normalized_weight
from scenarios import scenario
Obtain the specification for the default scenario The definition of the scenarios is available in Specification of a nested logit model.
V, _, _, _ = scenario()
V_PT = V[0]
V_CAR = V[1]
Calculation of the willingness to pay using derivatives.
WTP_PT_TIME = Derive(V_PT, 'TimePT') / Derive(V_PT, 'MarginalCostPT')
WTP_CAR_TIME = Derive(V_CAR, 'TimeCar') / Derive(V_CAR, 'CostCarCHF')
Formulas to simulate.
simulate = {
'weight': normalized_weight,
'WTP PT time': WTP_PT_TIME,
'WTP CAR time': WTP_CAR_TIME,
}
Read the data
database = read_data()
Create the Biogeme object.
the_biogeme = bio.BIOGEME(database, simulate)
Read the estimation results from the file.
try:
results = res.bioResults(pickle_file='saved_results/b02estimation.pickle')
except BiogemeError:
sys.exit(
'Run first the script b02estimation.py in order to generate '
'the file b02estimation.pickle.'
)
simulated_values is a Pandas dataframe with the same number of rows as the database, and as many columns as formulas to simulate.
simulated_values = the_biogeme.simulate(results.get_beta_values())
display(simulated_values)
weight WTP PT time WTP CAR time
0 0.886023 0.026513 0.026513
2 0.861136 0.026513 0.026513
3 0.861136 0.026513 0.026513
4 0.957386 0.086226 0.086226
5 0.861136 0.026513 0.026513
... ... ... ...
2259 2.036009 0.026513 0.026513
2261 0.861136 0.086226 0.086226
2262 0.861136 0.086226 0.086226
2263 0.957386 0.026513 0.026513
2264 0.957386 0.026513 0.026513
[1906 rows x 3 columns]
Note the multiplication by 60 to have the valus of time per hour and not per minute.
wtpcar = (60 * simulated_values['WTP CAR time'] * simulated_values['weight']).mean()
Calculate confidence intervals
b = results.get_betas_for_sensitivity_analysis(the_biogeme.free_beta_names)
Returns data frame containing, for each simulated value, the left and right bounds of the confidence interval calculated by simulation.
left, right = the_biogeme.confidence_intervals(b, 0.9)
Lower bounds of the confidence intervals
display(left)
weight WTP PT time WTP CAR time
0 0.886023 0.009167 0.009167
2 0.861136 0.009167 0.009167
3 0.861136 0.009167 0.009167
4 0.957386 0.071430 0.071430
5 0.861136 0.009167 0.009167
... ... ... ...
2259 2.036009 0.009167 0.009167
2261 0.861136 0.071430 0.071430
2262 0.861136 0.071430 0.071430
2263 0.957386 0.009167 0.009167
2264 0.957386 0.009167 0.009167
[1906 rows x 3 columns]
Upper bounds of the confidence intervals
display(right)
weight WTP PT time WTP CAR time
0 0.886023 0.066435 0.066435
2 0.861136 0.066435 0.066435
3 0.861136 0.066435 0.066435
4 0.957386 0.146445 0.146445
5 0.861136 0.066435 0.066435
... ... ... ...
2259 2.036009 0.066435 0.066435
2261 0.861136 0.146445 0.146445
2262 0.861136 0.146445 0.146445
2263 0.957386 0.066435 0.066435
2264 0.957386 0.066435 0.066435
[1906 rows x 3 columns]
Lower and upper bounds of the willingness to pay.
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}]'
)
Average WTP for car: 2.88 CI:[1.9, 5.72]
In this specific case, there are only two distinct values in the population: for workers and non workers
print(
'Unique values: ',
[f'{i:.3g}' for i in 60 * simulated_values['WTP CAR time'].unique()],
)
Unique values: ['1.59', '5.17']
Function calculating the willingness to pay for a group.
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
Full time workers.
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}]')
WTP car for workers: 5.17 CI:[4.29, 8.79]
Females.
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}]')
WTP car for females: 2.22 CI:[1.21, 4.83]
Males.
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}]')
WTP car for males : 3.72 CI:[2.77, 6.84]
We plot the distribution of WTP in the population. In this case, there are only two values
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()
Total running time of the script: (0 minutes 2.687 seconds)