{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7d2f605d",
   "metadata": {},
   "source": [
    "\n",
    "File: 05-wtp_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Sun Aug 03 2025, 17:30:36\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aeb9e49",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.expressions import Beta, Derive, exp, log, logzero\n",
    "from biogeme.models import loglogit\n",
    "from biogeme.results_processing import (\n",
    "    EstimationResults,\n",
    "    get_pandas_estimated_parameters,\n",
    ")\n",
    "from matplotlib import pyplot as plt\n",
    "from pandas import DataFrame\n",
    "\n",
    "from optima_variables import (\n",
    "    Choice,\n",
    "    CostCarCHF,\n",
    "    MarginalCostPT,\n",
    "    NbTransf,\n",
    "    TimeCar,\n",
    "    TimePT,\n",
    "    database,\n",
    "    distance_km,\n",
    "    female,\n",
    "    fulltime,\n",
    "    male,\n",
    "    normalizedWeight,\n",
    "    notfulltime,\n",
    "    unreportedGender,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe98f965",
   "metadata": {},
   "source": [
    "The objective of this laboratory is to calculate the willingness to pay for the improvement of some attributes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99554107",
   "metadata": {},
   "source": [
    "# First model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d9291ba",
   "metadata": {},
   "source": [
    "In the first model, the variables involved for the calculation of the willingness to pay appear linearly in the\n",
    "utility functions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d5c798",
   "metadata": {},
   "source": [
    "Parameters to be estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c07cd4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "asc_car = Beta('asc_car', 0, None, None, 0)\n",
    "asc_sm = Beta('asc_sm', 0, None, None, 0)\n",
    "beta_time_fulltime_pt = Beta('beta_time_fulltime_pt', 0, None, None, 0)\n",
    "beta_time_other_pt = Beta('beta_time_other_pt', 0, None, None, 0)\n",
    "beta_time_fulltime_car = Beta('beta_time_fulltime_car', 0, None, None, 0)\n",
    "beta_time_other_car = Beta('beta_time_other_car', 0, None, None, 0)\n",
    "\n",
    "beta_dist_male = Beta('beta_dist_male', 0, None, None, 0)\n",
    "beta_dist_female = Beta('beta_dist_female', 0, None, None, 0)\n",
    "beta_dist_unreported = Beta('beta_dist_unreported', 0, None, None, 0)\n",
    "beta_cost = Beta('beta_cost', 0, None, None, 0)\n",
    "beta_transf = Beta('beta_transf', 0, None, None, 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47c0ad86",
   "metadata": {},
   "source": [
    "Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44a8e273",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_pt = (\n",
    "    beta_time_fulltime_pt * TimePT * fulltime\n",
    "    + beta_time_other_pt * TimePT * notfulltime\n",
    "    + beta_cost * MarginalCostPT\n",
    "    + beta_transf * NbTransf\n",
    ")\n",
    "v_car = (\n",
    "    asc_car\n",
    "    + beta_time_fulltime_car * TimeCar * fulltime\n",
    "    + beta_time_other_car * TimeCar * notfulltime\n",
    "    + beta_cost * CostCarCHF\n",
    ")\n",
    "v_sm = (\n",
    "    asc_sm\n",
    "    + beta_dist_male * distance_km * male\n",
    "    + beta_dist_female * distance_km * female\n",
    "    + beta_dist_unreported * distance_km * unreportedGender\n",
    ")\n",
    "v = {0: v_pt, 1: v_car, 2: v_sm}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "379f0998",
   "metadata": {},
   "source": [
    "Estimation of the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59bcbc0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v, None, Choice)\n",
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = 'wtp_first_model'\n",
    "results: EstimationResults = biogeme.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd53320a",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c818d3bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46ff4133",
   "metadata": {},
   "source": [
    "Estimated parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52f78d4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = get_pandas_estimated_parameters(estimation_results=results)\n",
    "display(parameters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca3e75e1",
   "metadata": {},
   "source": [
    "## Value of time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aefbc5d4",
   "metadata": {},
   "source": [
    "We  calculate the value of time, that is the willingness to pay to save travel time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53faf209",
   "metadata": {},
   "source": [
    "Public transportation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fcc662b",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_pt = Derive(v_pt, 'TimePT') / Derive(v_pt, 'MarginalCostPT')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6f7d33d",
   "metadata": {},
   "source": [
    "Car"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bc847f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_car = Derive(v_car, 'TimeCar') / Derive(v_car, 'CostCarCHF')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bc6456b",
   "metadata": {},
   "source": [
    "Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d01dff8",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = {\n",
    "    'weight': normalizedWeight,\n",
    "    'WTP PT time': vot_pt,\n",
    "    'WTP CAR time': vot_car,\n",
    "}\n",
    "biosim = BIOGEME(database, simulate)\n",
    "simulated_values = biosim.simulate(results.get_beta_values())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e621dcb",
   "metadata": {},
   "source": [
    "We first look at the value of time for the public transportation alternative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5452d538",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values['WTP PT time'].hist()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe9dd337",
   "metadata": {},
   "source": [
    "There are two values for the willingness to pay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03984a4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(simulated_values['WTP PT time'].unique())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6141cd85",
   "metadata": {},
   "source": [
    "They correspond to the two segments of the population characterized by `fulltime` and `notfulltime`. As the\n",
    "specification of the utility function is linear both in `TimePT` and `MarginalCostPT`, the value of time is simply\n",
    "the ratio of the coefficients."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8153085f",
   "metadata": {},
   "source": [
    "For the segment `fulltime`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f57599a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator = results.get_parameter_value('beta_time_fulltime_pt')\n",
    "denominator = results.get_parameter_value('beta_cost')\n",
    "vot_pt_fulltime_min = numerator / denominator\n",
    "print(f'VOT public transportation / full time: {vot_pt_fulltime_min:.3g} CHF/min.')\n",
    "print(f'VOT public transportation / full time: {vot_pt_fulltime_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33b8b034",
   "metadata": {},
   "source": [
    "According to this specification, individuals in this segment are willing to pay 0.208 CHF to save one minute of\n",
    "travel time, that is 12.5 CHF per hour."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502de426",
   "metadata": {},
   "source": [
    "For the segment `notfulltime`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be57d18e",
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator = results.get_parameter_value('beta_time_other_pt')\n",
    "denominator = results.get_parameter_value('beta_cost')\n",
    "vot_pt_not_fulltime_min = numerator / denominator\n",
    "print(f'VOT public transportation / part time: {vot_pt_not_fulltime_min:.3g} CHF/min.')\n",
    "print(\n",
    "    f'VOT public transportation / part time: {vot_pt_not_fulltime_min * 60:.3g} CHF/h.'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f29ceb1f",
   "metadata": {},
   "source": [
    "According to this specification, individuals in this segment are willing to pay 0.144 CHF to save one minute of\n",
    "travel time, that is 8.67 CHF per hour."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77214daf",
   "metadata": {},
   "source": [
    "We can do the same for the car alternative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58acb5e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values['WTP CAR time'].hist()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac22254b",
   "metadata": {},
   "source": [
    "There are also two values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f09a648d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(simulated_values['WTP CAR time'].unique())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0a4952d",
   "metadata": {},
   "source": [
    "For the segment `fulltime`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdbb0daf",
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator = results.get_parameter_value('beta_time_fulltime_car')\n",
    "denominator = results.get_parameter_value('beta_cost')\n",
    "vot_car_fulltime_min = numerator / denominator\n",
    "print(f'VOT car / full time: {vot_car_fulltime_min:.3g} CHF/min.')\n",
    "print(f'VOT car / full time: {vot_car_fulltime_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0530c9ce",
   "metadata": {},
   "source": [
    "For the segment `notfulltime`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "815be950",
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator = results.get_parameter_value('beta_time_other_car')\n",
    "denominator = results.get_parameter_value('beta_cost')\n",
    "vot_car_not_fulltime_min = numerator / denominator\n",
    "print(f'VOT car / part time: {vot_car_not_fulltime_min:.3g} CHF/min.')\n",
    "print(f'VOT car / part time: {vot_car_not_fulltime_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26a409f",
   "metadata": {},
   "source": [
    "The value of time for the car alternative is 27.9 CHF/hour for individuals in the segment `fulltime` and\n",
    "28.3 CHF/hour for individuals in the segment `nofulltime`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cfe396c",
   "metadata": {},
   "source": [
    "Finally, we can calculate the average value of time in the population."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "613c697d",
   "metadata": {},
   "source": [
    "For public transportation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5963e7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_vot_pt = (\n",
    "    simulated_values['WTP PT time'] * simulated_values['weight']\n",
    ").sum() / simulated_values['weight'].sum()\n",
    "print(f'Average value of time for PT: {60 * avg_vot_pt:.3g} CHF/hour')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1496301e",
   "metadata": {},
   "source": [
    "For car:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0029372",
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_vot_car = (\n",
    "    simulated_values['WTP CAR time'] * simulated_values['weight']\n",
    ").sum() / simulated_values['weight'].sum()\n",
    "print(f'Average value of time for car: {60 * avg_vot_car:.3g} CHF/hour')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4765d7c",
   "metadata": {},
   "source": [
    "## Willingness to pay to reduce the number of transfers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ae8d44d",
   "metadata": {},
   "source": [
    "In this case, as the variable is discrete, the derivatives cannot be used. The wtp to decrease the number of\n",
    "transfers is the value $population_shares$ that solves the equation:\n",
    "$$\\beta_\\text{transfer} \\text{NbTransf} + \\beta_\\text{cost} \\text{MarginalCostPT} =\n",
    "\\beta_\\text{transfer} (\\text{NbTransf}-1) + \\beta_\\text{cost} (\\text{MarginalCostPT} + population_shares).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "903ce1f0",
   "metadata": {},
   "source": [
    "In this case, the solution happens to be also the ratio of coefficients:\n",
    "$$population_shares=\\frac{\\beta_\\text{transfer}}{\\beta_\\text{cost}}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a12861c",
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator = results.get_parameter_value('beta_transf')\n",
    "denominator = results.get_parameter_value('beta_cost')\n",
    "wtp_transfer = numerator / denominator\n",
    "print(\n",
    "    f'Willingness to pay to reduce the number of transfers: {wtp_transfer:.3g} CHF/transfer'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33d279c8",
   "metadata": {},
   "source": [
    "# Second model: moneymetric specification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f98fd894",
   "metadata": {},
   "source": [
    "In the moneymetric specification, the cost parameter is normalized to -1, and the scale parameter is estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59e551f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "scale_parameter = Beta('scale_parameter', 1, 1.0e-5, None, 0)\n",
    "\n",
    "v_pt_moneymetric = scale_parameter * (\n",
    "    beta_time_fulltime_pt * TimePT * fulltime\n",
    "    + beta_time_other_pt * TimePT * notfulltime\n",
    "    - MarginalCostPT\n",
    "    + beta_transf * NbTransf\n",
    ")\n",
    "v_car_moneymetric = scale_parameter * (\n",
    "    asc_car\n",
    "    + beta_time_fulltime_car * TimeCar * fulltime\n",
    "    + beta_time_other_car * TimeCar * notfulltime\n",
    "    - CostCarCHF\n",
    ")\n",
    "v_sm_moneymetric = scale_parameter * (\n",
    "    asc_sm\n",
    "    + beta_dist_male * distance_km * male\n",
    "    + beta_dist_female * distance_km * female\n",
    "    + beta_dist_unreported * distance_km * unreportedGender\n",
    ")\n",
    "v_moneymetric = {0: v_pt_moneymetric, 1: v_car_moneymetric, 2: v_sm_moneymetric}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e183185",
   "metadata": {},
   "source": [
    "Estimation of the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f590e9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v_moneymetric, None, Choice)\n",
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = 'wtp_second_model'\n",
    "results_second: EstimationResults = biogeme.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4fcbd3a",
   "metadata": {},
   "source": [
    "General statistics. Note that the final log likelihood is the same as for the first model. It makes sense as\n",
    "the two models are equivalent, up to a normalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea3a3b29",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_second.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b6b53a",
   "metadata": {},
   "source": [
    "Estimated parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4966e90b",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_second = get_pandas_estimated_parameters(estimation_results=results_second)\n",
    "display(parameters_second)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e03b97c",
   "metadata": {},
   "source": [
    "As the utility function is expressed in CHF units, the parameters provide us directly with the WTP,\n",
    "with the opposite sign. You can verify that the values are  the same as for the first model (up to possible\n",
    "rounding errors)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c9b944",
   "metadata": {},
   "source": [
    "Value of time, public transportation, full time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "156feff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_pt_full_time_min = -results_second.get_parameter_value('beta_time_fulltime_pt')\n",
    "print(f'VOT public transportation / full time: {vot_pt_full_time_min:.3g} CHF/min.')\n",
    "print(f'VOT public transportation / full time: {vot_pt_full_time_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a83d19b6",
   "metadata": {},
   "source": [
    "Value of time, public transportation, part time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "044ac608",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_pt_part_time_min = -results_second.get_parameter_value('beta_time_other_pt')\n",
    "print(f'VOT public transportation / part time: {vot_pt_part_time_min:.3g} CHF/min.')\n",
    "print(f'VOT public transportation / part time: {vot_pt_part_time_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7e38587",
   "metadata": {},
   "source": [
    "Value of time, car, full time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37264133",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_car_full_time_min = -results_second.get_parameter_value('beta_time_fulltime_car')\n",
    "print(f'VOT car / full time: {vot_car_full_time_min:.3g} CHF/min.')\n",
    "print(f'VOT car / full time: {vot_car_full_time_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5db88f36",
   "metadata": {},
   "source": [
    "Value of time, car, part time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e951844",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_car_part_time_min = -results_second.get_parameter_value('beta_time_other_car')\n",
    "print(f'VOT car / part time: {vot_car_part_time_min:.3g} CHF/min.')\n",
    "print(f'VOT car / part time: {vot_car_part_time_min * 60:.3g} CHF/h.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d560795",
   "metadata": {},
   "source": [
    "Willingness to pay to reduce the transfers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c456920a",
   "metadata": {},
   "outputs": [],
   "source": [
    "wtp_transfer = -results_second.get_parameter_value('beta_transf')\n",
    "print(\n",
    "    f'Willingness to pay to reduce the number of transfers: {wtp_transfer:.3g} CHF/transfer'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6acaf69",
   "metadata": {},
   "source": [
    "Clearly, the aggregate values for the population can be calculated exactly like for the first model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fe0f76f",
   "metadata": {},
   "source": [
    "# Third model: nonlinear transformation of the variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8ec38f5",
   "metadata": {},
   "source": [
    "We use a Biogeme expression that returns the log of a variable if it is nonzero, and 0 otherwise.\n",
    "The variable is assumed to be non-negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "177888a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_pt = (\n",
    "    beta_time_fulltime_pt * logzero(TimePT) * fulltime\n",
    "    + beta_time_other_pt * logzero(TimePT) * notfulltime\n",
    "    + beta_cost * logzero(MarginalCostPT)\n",
    "    + beta_transf * logzero(NbTransf)\n",
    ")\n",
    "v_car = (\n",
    "    asc_car\n",
    "    + beta_time_fulltime_car * logzero(TimeCar) * fulltime\n",
    "    + beta_time_other_car * logzero(TimeCar) * notfulltime\n",
    "    + beta_cost * logzero(CostCarCHF)\n",
    ")\n",
    "v_sm = (\n",
    "    asc_sm\n",
    "    + beta_dist_male * distance_km * male\n",
    "    + beta_dist_female * distance_km * female\n",
    "    + beta_dist_unreported * distance_km * unreportedGender\n",
    ")\n",
    "v = {0: v_pt, 1: v_car, 2: v_sm}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d00a5282",
   "metadata": {},
   "source": [
    "Estimation of the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ad31751",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v, None, Choice)\n",
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = 'wtp_third_model'\n",
    "results_third: EstimationResults = biogeme.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65010596",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a836e9dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_third.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00a1a09f",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "230e0707",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_third = get_pandas_estimated_parameters(estimation_results=results_third)\n",
    "display(parameters_third)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95b30352",
   "metadata": {},
   "source": [
    "## Value of time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b398224",
   "metadata": {},
   "source": [
    "We calculate the value of time, that is the willingness to pay to save travel time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9273821d",
   "metadata": {},
   "outputs": [],
   "source": [
    "vot_pt = Derive(v_pt, 'TimePT') / Derive(v_pt, 'MarginalCostPT')\n",
    "vot_car = Derive(v_car, 'TimeCar') / Derive(v_car, 'CostCarCHF')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3ecc225",
   "metadata": {},
   "source": [
    "Simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8878f47",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = {\n",
    "    'weight': normalizedWeight,\n",
    "    'WTP PT time': vot_pt,\n",
    "    'WTP CAR time': vot_car,\n",
    "    'MarginalCostPT': MarginalCostPT,\n",
    "    'CostCarCHF': CostCarCHF,\n",
    "}\n",
    "biosim = BIOGEME(database, simulate)\n",
    "simulated_values = DataFrame(biosim.simulate(results_third.get_beta_values()))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a21b1404",
   "metadata": {},
   "source": [
    "Note that the formula are not valid when the cost is 0. In that case, Biogeme generates `inf` or `-inf`.\n",
    "For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd889b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(simulated_values.iloc[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c8abb5f",
   "metadata": {},
   "source": [
    "If the cost is zero, we set the value of time to zero as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58347d86",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_values['WTP PT time'] = simulated_values['WTP PT time'].where(\n",
    "    simulated_values['MarginalCostPT'] != 0, 0\n",
    ")\n",
    "simulated_values['WTP CAR time'] = simulated_values['WTP CAR time'].where(\n",
    "    simulated_values['CostCarCHF'] != 0, 0\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5f5cdba",
   "metadata": {},
   "source": [
    "We check that the values have been set correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ce49775",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(simulated_values.iloc[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2595a015",
   "metadata": {},
   "source": [
    "We first look at the value of time for the public transportation alternative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30a482d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values['WTP PT time'].hist(bins=200)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53b3e3ae",
   "metadata": {},
   "source": [
    "There are now many different values for the value of time. Indeed, because of the nonlinear specification, the\n",
    "value of time depends also on the value of the variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7172c146",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(simulated_values['WTP PT time'].value_counts())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e3b6c13",
   "metadata": {},
   "source": [
    "Consider only the nonzero values for a clearer plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1978241f",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values[simulated_values['WTP PT time'] != 0]['WTP PT time'].hist(bins=200)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "789dc67a",
   "metadata": {},
   "source": [
    "We can do the same for the car alternative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2695f0b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values['WTP CAR time'].hist(bins=200)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cba3838",
   "metadata": {},
   "source": [
    "Removing the zeros..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "455b3040",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values[simulated_values['WTP CAR time'] != 0]['WTP CAR time'].hist(\n",
    "    bins=200\n",
    ")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "602b3d65",
   "metadata": {},
   "source": [
    "Finally, we can calculate the average value of time in the population."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a968ea4f",
   "metadata": {},
   "source": [
    "For public transportation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41263f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_vot_pt = (\n",
    "    simulated_values['WTP PT time'] * simulated_values['weight']\n",
    ").sum() / simulated_values['weight'].sum()\n",
    "print(f'Average value of time for PT: {60 * avg_vot_pt:.3g} CHF/hour')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92d09dc1",
   "metadata": {},
   "source": [
    "For car:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b25c7337",
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_vot_car = (\n",
    "    simulated_values['WTP CAR time'] * simulated_values['weight']\n",
    ").sum() / simulated_values['weight'].sum()\n",
    "print(f'Average value of time for car: {60 * avg_vot_car:.3g} CHF/hour')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c17a6c24",
   "metadata": {},
   "source": [
    "Those values are abnormally low, triggering suspicion on the validity of the model specification."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6718e56b",
   "metadata": {},
   "source": [
    "## Willingness to pay to reduce the number of transfers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64e1a9ec",
   "metadata": {},
   "source": [
    "Again, as the variable is discrete, the derivatives cannot be used. The wtp to decrease the number of transfers is\n",
    "the value $population_shares$ that solves the equation:\n",
    "$$\\beta_\\text{transfer} \\log(\\text{NbTransf}) + \\beta_\\text{cost} \\log(\\text{MarginalCostPT}) =\n",
    "\\beta_\\text{transfer} \\log(\\text{NbTransf}-1) + \\beta_\\text{cost} \\log(\\text{MarginalCostPT} + population_shares).$$\n",
    "\n",
    "The solution is: $$population_shares=\\exp\\left(\\frac{\\beta_\\text{transfer}}{\\beta_\\text{cost}} [\\log(\\text{NbTransf}) -\n",
    "\\log(\\text{NbTransf}-1)] + \\log(\\text{MarginalCostPT})\\right)-\\text{MarginalCostPT}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29293f8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "wtp_transfer = (\n",
    "    exp(\n",
    "        (log(NbTransf) - log(NbTransf - 1)) * beta_transf / beta_cost\n",
    "        + log(MarginalCostPT)\n",
    "    )\n",
    "    - MarginalCostPT\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0991c3b",
   "metadata": {},
   "source": [
    "Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d88420c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = {\n",
    "    'weight': normalizedWeight,\n",
    "    'WTP Transfers': wtp_transfer,\n",
    "    'Nbr. of transfers': NbTransf,\n",
    "    'Cost PT': MarginalCostPT,\n",
    "}\n",
    "\n",
    "\n",
    "biosim = BIOGEME(database, simulate)\n",
    "simulated_values = biosim.simulate(results_third.get_beta_values())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "932398c9",
   "metadata": {},
   "source": [
    "Note that we do not verify if the arguments of the `log` is involved in the formula are valid, that is, positive.\n",
    "When they are not, Biogeme generates either a `NaN` or `inf`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "163c360b",
   "metadata": {},
   "source": [
    "For example, the following entry is invalid as the number of transfers is 1, and it is not possible to take the\n",
    "log of 1-1=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb066731",
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_index_1 = 5\n",
    "print(f'Simulated values: {simulated_values.loc[bad_index_1]}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f46cc55",
   "metadata": {},
   "source": [
    "The following entry is invalid as the number of transfers is 0, and it is not possible to take the log of 0-1=-1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff4b49d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_index_0 = 1903\n",
    "print(simulated_values.loc[bad_index_0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b214a8f0",
   "metadata": {},
   "source": [
    "We set the willingness to pay to zero in these cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e718c26a",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_values['WTP Transfers'] = simulated_values['WTP Transfers'].replace(\n",
    "    [np.nan, np.inf, -np.inf], 0\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8faf95ff",
   "metadata": {},
   "source": [
    "The willingness to pay for problematic cases is now set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8c4d10e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(simulated_values.loc[bad_index_1])\n",
    "print(simulated_values.loc[bad_index_0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c24eaea3",
   "metadata": {},
   "source": [
    "We can now plot the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c4ede2",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values['WTP Transfers'].hist(bins=100)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15be189d",
   "metadata": {},
   "source": [
    "As well as the distribution of nonzero values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d90d27fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = simulated_values[simulated_values['WTP Transfers'] > 0]['WTP Transfers'].hist(\n",
    "    bins=100\n",
    ")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}