Calculation of revenues

We use an estimated model to calculate revenues.

author:

Michel Bierlaire, EPFL

date:

Wed Apr 12 21:02:19 2023

import sys
import numpy as np

try:
    import matplotlib.pyplot as plt

    can_plot = True
except ModuleNotFoundError:
    can_plot = False
from tqdm import tqdm
from biogeme import models
from biogeme.exceptions import BiogemeError
import biogeme.biogeme as bio
import biogeme.results as res
from biogeme.data.optima import read_data, normalized_weight
from scenarios import scenario

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 b02simulation.py '
        'in order to generate the '
        'file b02estimation.pickle.'
    )

Read the data

database = read_data()

Function calculating the revenues

def revenues(factor: float) -> tuple[float, float, float]:
    """Calculate the total revenues generated by public transportation,
        when the price is multiplied by a factor.

    :param factor: factor that multiplies the current cost of public
        transportation

    :return: total revenues, followed by the lower and upper bound of
        the confidence interval.

    """
    # Obtain the specification for the default scenario
    V, nests, _, marginal_cost_scenario = scenario(factor=factor)

    # Obtain the expression for the choice probability of each alternative
    prob_pt = models.nested(V, None, nests, 0)

    # We now simulate the choice probabilities,the weight and the
    # price variable

    simulate = {
        'weight': normalized_weight,
        'Revenue public transportation': prob_pt * marginal_cost_scenario,
    }

    the_biogeme = bio.BIOGEME(database, simulate)
    simulated_values = the_biogeme.simulate(results.get_beta_values())

    # We also calculate confidence intervals for the calculated quantities

    betas = the_biogeme.free_beta_names
    beta_bootstrap = results.get_betas_for_sensitivity_analysis(betas)
    left, right = the_biogeme.confidence_intervals(beta_bootstrap, 0.9)

    revenues_pt = (
        simulated_values['Revenue public transportation'] * simulated_values['weight']
    ).sum()
    revenues_pt_left = (left['Revenue public transportation'] * left['weight']).sum()
    revenues_pt_right = (right['Revenue public transportation'] * right['weight']).sum()
    return revenues_pt, revenues_pt_left, revenues_pt_right

Current revenues for public transportation

r, r_left, r_right = revenues(factor=1.0)
print(
    f'Total revenues for public transportation (for the sample): {r:.1f} CHF '
    f'[{r_left:.1f} CHF, '
    f'{r_right:.1f} CHF]'
)

Total revenues for public transportation (for the sample): 2878.3 CHF [2360.6 CHF, 3424.8 CHF]

We now investigate how the revenues vary with the multiplicative factor

factors = np.arange(0.0, 5.0, 0.1)
plot_revenues = [revenues(s) for s in tqdm(factors)]
zipped = zip(*plot_revenues)
rev = next(zipped)
lower = next(zipped)
upper = next(zipped)

largest_revenue = max(rev)
max_index = rev.index(largest_revenue)

  0%|          | 0/50 [00:00<?, ?it/s]
  2%|▏         | 1/50 [00:11<09:25, 11.54s/it]
  4%|▍         | 2/50 [00:22<08:50, 11.06s/it]
  6%|▌         | 3/50 [00:33<08:37, 11.00s/it]
  8%|▊         | 4/50 [00:41<07:28,  9.75s/it]
 10%|█         | 5/50 [00:47<06:23,  8.52s/it]
 12%|█▏        | 6/50 [00:53<05:44,  7.82s/it]
 14%|█▍        | 7/50 [01:00<05:16,  7.35s/it]
 16%|█▌        | 8/50 [01:06<04:56,  7.06s/it]
 18%|█▊        | 9/50 [01:13<04:41,  6.86s/it]
 20%|██        | 10/50 [01:19<04:28,  6.71s/it]
 22%|██▏       | 11/50 [01:25<04:19,  6.65s/it]
 24%|██▍       | 12/50 [01:32<04:10,  6.59s/it]
 26%|██▌       | 13/50 [01:38<04:02,  6.55s/it]
 28%|██▊       | 14/50 [01:45<03:55,  6.53s/it]
 30%|███       | 15/50 [01:52<03:55,  6.74s/it]
 32%|███▏      | 16/50 [02:05<04:55,  8.69s/it]
 34%|███▍      | 17/50 [02:18<05:23,  9.79s/it]
 36%|███▌      | 18/50 [02:30<05:41, 10.68s/it]
 38%|███▊      | 19/50 [02:42<05:40, 10.99s/it]
 40%|████      | 20/50 [02:53<05:31, 11.04s/it]
 42%|████▏     | 21/50 [03:05<05:22, 11.13s/it]
 44%|████▍     | 22/50 [03:16<05:09, 11.07s/it]
 46%|████▌     | 23/50 [03:26<04:50, 10.75s/it]
 48%|████▊     | 24/50 [03:35<04:26, 10.23s/it]
 50%|█████     | 25/50 [03:45<04:17, 10.30s/it]
 52%|█████▏    | 26/50 [03:57<04:17, 10.71s/it]
 54%|█████▍    | 27/50 [04:08<04:10, 10.91s/it]
 56%|█████▌    | 28/50 [04:17<03:49, 10.43s/it]
 58%|█████▊    | 29/50 [04:24<03:13,  9.20s/it]
 60%|██████    | 30/50 [04:30<02:47,  8.35s/it]
 62%|██████▏   | 31/50 [04:36<02:27,  7.74s/it]
 64%|██████▍   | 32/50 [04:43<02:12,  7.34s/it]
 66%|██████▌   | 33/50 [04:49<02:00,  7.08s/it]
 68%|██████▊   | 34/50 [04:56<01:50,  6.89s/it]
 70%|███████   | 35/50 [05:02<01:41,  6.77s/it]
 72%|███████▏  | 36/50 [05:09<01:33,  6.67s/it]
 74%|███████▍  | 37/50 [05:15<01:26,  6.63s/it]
 76%|███████▌  | 38/50 [05:22<01:19,  6.59s/it]
 78%|███████▊  | 39/50 [05:28<01:12,  6.56s/it]
 80%|████████  | 40/50 [05:37<01:11,  7.11s/it]
 82%|████████▏ | 41/50 [05:49<01:18,  8.74s/it]
 84%|████████▍ | 42/50 [06:01<01:18,  9.79s/it]
 86%|████████▌ | 43/50 [06:15<01:16, 10.96s/it]
 88%|████████▊ | 44/50 [06:26<01:06, 11.11s/it]
 90%|█████████ | 45/50 [06:37<00:53, 10.79s/it]
 92%|█████████▏| 46/50 [06:47<00:42, 10.72s/it]
 94%|█████████▍| 47/50 [06:57<00:31, 10.59s/it]
 96%|█████████▌| 48/50 [07:08<00:21, 10.53s/it]
 98%|█████████▊| 49/50 [07:17<00:10, 10.11s/it]
100%|██████████| 50/50 [07:28<00:00, 10.28s/it]
100%|██████████| 50/50 [07:28<00:00,  8.96s/it]

print(
    f'Largest revenue: {largest_revenue:.1f} obtained with '
    f'factor {factors[max_index]:.1f}'
)

Largest revenue: 2879.6 obtained with factor 0.9

if can_plot:
    # We plot the results
    ax = plt.gca()
    ax.plot(factors, rev, label="Revenues")
    ax.plot(factors, lower, label="Lower bound of the CI")
    ax.plot(factors, upper, label="Upper bound of the CI")
    ax.legend()

    plt.show()