{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#%%\n",
        "\"\"\"File 02nestedPlot.py\n",
        "\n",
        ":author: Michel Bierlaire, EPFL\n",
        ":date: Wed Sep 11 10:15:18 2019\n",
        "\n",
        " We use a previously estimated nested logit model.\n",
        " Three alternatives: public transporation, car and slow modes.\n",
        " RP data.\n",
        " We simulate pricing scenarios and their impact on the revenues.\n",
        "\"\"\"\n",
        "\n",
        "import sys\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import biogeme.database as db\n",
        "import biogeme.biogeme as bio\n",
        "from biogeme import models\n",
        "import biogeme.results as res\n",
        "from biogeme.expressions import Beta\n",
        "\n",
        "# Library for plotting\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Read the data\n",
        "df = pd.read_csv('optima.dat', sep='\\t')\n",
        "database = db.Database('optima', df)\n",
        "\n",
        "# The following statement allows you to use the names of the variable\n",
        "# as Python variable.\n",
        "globals().update(database.variables)\n",
        "\n",
        "# Exclude observations such that the chosen alternative is -1\n",
        "database.remove(Choice == -1.0)\n",
        "\n",
        "# Define new variables. Must be consistent with estimation results.\n",
        "TimePT_scaled = TimePT / 200\n",
        "TimeCar_scaled = TimeCar / 200\n",
        "CostCarCHF_scaled = CostCarCHF / 10\n",
        "distance_km_scaled = distance_km / 5\n",
        "male = Gender == 1\n",
        "female = Gender == 2\n",
        "unreportedGender = Gender == -1\n",
        "fulltime = OccupStat == 1\n",
        "notfulltime = OccupStat != 1\n",
        "\n",
        "# Normalize the weights\n",
        "sumWeight = database.data['Weight'].sum()\n",
        "numberOfRows = database.data.shape[0]\n",
        "normalizedWeight = Weight * numberOfRows / sumWeight\n",
        "\n",
        "\n",
        "def scenario(scale):\n",
        "    \"\"\"Simulate a scenarios modifying the price of public transportation\n",
        "\n",
        "    :param scale: price multiplier.\n",
        "    :type scale: float\n",
        "\n",
        "    :return: simulated revenues\n",
        "    :rtype: float\n",
        "    \"\"\"\n",
        "    # This is the only variable that depends on scale\n",
        "    MarginalCostScenario = MarginalCostPT * scale\n",
        "    MarginalCostPT_scaled = MarginalCostScenario / 10\n",
        "\n",
        "    # The rest of the model is the same for all scenarios\n",
        "    ASC_CAR = Beta('ASC_CAR', 0, None, None, 0)\n",
        "    ASC_PT = Beta('ASC_PT', 0, None, None, 1)\n",
        "    ASC_SM = Beta('ASC_SM', 0, None, None, 0)\n",
        "    BETA_TIME_FULLTIME = Beta('BETA_TIME_FULLTIME', 0, None, None, 0)\n",
        "    BETA_TIME_OTHER = Beta('BETA_TIME_OTHER', 0, None, None, 0)\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",
        "    # Utility functions\n",
        "    V_PT = (\n",
        "        ASC_PT\n",
        "        + BETA_TIME_FULLTIME * TimePT_scaled * fulltime\n",
        "        + BETA_TIME_OTHER * TimePT_scaled * notfulltime\n",
        "        + BETA_COST * MarginalCostPT_scaled\n",
        "    )\n",
        "    V_CAR = (\n",
        "        ASC_CAR\n",
        "        + BETA_TIME_FULLTIME * TimeCar_scaled * fulltime\n",
        "        + BETA_TIME_OTHER * TimeCar_scaled * notfulltime\n",
        "        + BETA_COST * CostCarCHF_scaled\n",
        "    )\n",
        "    V_SM = (\n",
        "        ASC_SM\n",
        "        + BETA_DIST_MALE * distance_km_scaled * male\n",
        "        + BETA_DIST_FEMALE * distance_km_scaled * female\n",
        "        + BETA_DIST_UNREPORTED * distance_km_scaled * unreportedGender\n",
        "    )\n",
        "    V = {0: V_PT, 1: V_CAR, 2: V_SM}\n",
        "    MU_NOCAR = Beta('MU_NOCAR', 1.0, 1.0, None, 0)\n",
        "    CAR_NEST = 1.0, [1]\n",
        "    NO_CAR_NEST = MU_NOCAR, [0, 2]\n",
        "    nests = CAR_NEST, NO_CAR_NEST\n",
        "    prob_pt = models.nested(V, None, nests, 0)\n",
        "    simulate = {\n",
        "        'weight': normalizedWeight,\n",
        "        'Revenue public transportation': prob_pt * MarginalCostScenario,\n",
        "    }\n",
        "\n",
        "    biogeme = bio.BIOGEME(database, simulate)\n",
        "    biogeme.modelName = '02nestedPlot'\n",
        "\n",
        "    # Read the estimation results from the file\n",
        "    try:\n",
        "        results = res.bioResults(pickleFile='01nestedEstimation.pickle')\n",
        "    except FileNotFoundError:\n",
        "        sys.exit(\n",
        "            'Run first the script 01nestedEstimation.py in order to generate '\n",
        "            'the file 01nestedEstimation.pickle.'\n",
        "        )\n",
        "    # Simulation\n",
        "    simulatedValues = biogeme.simulate(results.getBetaValues())\n",
        "\n",
        "    # We calculate the sum for all individuals of the generated revenues.\n",
        "    revenues_pt = (\n",
        "        simulatedValues['Revenue public transportation']\n",
        "        * simulatedValues['weight']\n",
        "    ).sum()\n",
        "    return revenues_pt\n",
        "\n",
        "\n",
        "# We now plot the relationship between the cost and the revenues\n",
        "scales = np.arange(0.0, 5.0, 0.1)\n",
        "revenues = [scenario(s) for s in scales]\n",
        "plt.plot(scales, revenues)\n",
        "plt.xlabel('Modification of the price of public transportation (%)')\n",
        "plt.ylabel('Revenues')\n",
        "plt.show()\n"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}