<NME> Nested Logit Example--Python Biogeme benchmark--09NestedLogit.ipynb
<BEF> {
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this notebook is twofold. First, it demonstrates the basic functionality of PyLogit for estimating nested logit models. Secondly, it compares the nested logit capabilities of PyLogit with Python Biogeme. The dataset used is the SwissMetro dataset from <a href=\"http://biogeme.epfl.ch/examples_swissmetro.html\">http://biogeme.epfl.ch/examples_swissmetro.html</a>. For an explanation of the variables in the dataset, see http://www.strc.ch/conferences/2001/bierlaire1.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict    # For recording the model specification \n",
    "\n",
    "import pandas as pd                    # For file input/output\n",
    "import numpy as np                     # For vectorized math operations\n",
    "import statsmodels.tools.numdiff as numdiff       # For numeric hessian\n",
    "import scipy.linalg                    # For matrix inversion\n",
    "\n",
    "import pylogit as pl                   # For choice model estimation\n",
    "from pylogit import nested_logit as nl # For nested logit convenience funcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Load the Swissmetro Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the raw swiss metro data\n",
    "# Note the .dat files are tab delimited text files\n",
    "swissmetro_wide = pd.read_table(\"../data/swissmetro.dat\", sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Clean the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the 09NestedLogit.py file provided is an example from Python Biogeme (see: <a href=\"http://biogeme.epfl.ch/examples_swissmetro.html\">http://biogeme.epfl.ch/examples_swissmetro.html</a>). The 09NestedLogit.py file excludes observations meeting the following critera:\n",
    "<pre>\n",
    "exclude = (( PURPOSE != 1 ) * (  PURPOSE   !=  3  ) + ( CHOICE == 0 )) > 0\n",
    "</pre>\n",
    "As a result, their dataset has 6,768 observations. Below, I make the same exclusions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The cleaned number of observations is 6,768.\n"
     ]
    }
   ],
   "source": [
    "# Select obervations whose choice is known (i.e. CHOICE != 0)\n",
    "# **AND** whose PURPOSE is either 1 or 3\n",
    "include_criteria = (swissmetro_wide.PURPOSE.isin([1, 3]) &\n",
    "                    (swissmetro_wide.CHOICE != 0))\n",
    "\n",
    "# Use \".copy()\" so that later on, we avoid performing operations \n",
    "# on a view of a dataframe as opposed to on an actual dataframe\n",
    "clean_sm_wide = swissmetro_wide.loc[include_criteria].copy()\n",
    "\n",
    "# Look at how many observations we have after removing unwanted\n",
    "# observations\n",
    "final_num_obs = clean_sm_wide.shape[0]\n",
    "num_obs_statement = \"The cleaned number of observations is {:,.0f}.\"\n",
    "print(num_obs_statement.format(final_num_obs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Create an id column that ignores the repeat observations per individual\n",
    "\n",
    "In the simple example given on the Python Biogeme website for 09NestedLogit.py, the repeated observations per individual are treated as separate and independent observations. We will do the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a custom id column that ignores the fact that this is a \n",
    "# panel/repeated-observations dataset, and start the \"custom_id\" from 1\n",
    "clean_sm_wide[\"custom_id\"] = np.arange(clean_sm_wide.shape[0], dtype=int) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Convert the data from 'wide' format to 'long' format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4a. Determine the 'type' of each column in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['GROUP', 'SURVEY', 'SP', 'ID', 'PURPOSE', 'FIRST', 'TICKET', 'WHO',\n",
       "       'LUGGAGE', 'AGE', 'MALE', 'INCOME', 'GA', 'ORIGIN', 'DEST', 'TRAIN_AV',\n",
       "       'CAR_AV', 'SM_AV', 'TRAIN_TT', 'TRAIN_CO', 'TRAIN_HE', 'SM_TT', 'SM_CO',\n",
       "       'SM_HE', 'SM_SEATS', 'CAR_TT', 'CAR_CO', 'CHOICE', 'custom_id'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the columns of the swissmetro data\n",
    "clean_sm_wide.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the list of individual specific variables\n",
    "ind_variables = clean_sm_wide.columns.tolist()[:15]\n",
    "\n",
    "# Specify the variables that vary across individuals **AND** \n",
    "# across some or all alternatives\n",
    "alt_varying_variables = {u'travel_time': dict([(1, 'TRAIN_TT'),\n",
    "                                               (2, 'SM_TT'),\n",
    "                                               (3, 'CAR_TT')]),\n",
    "                          u'travel_cost': dict([(1, 'TRAIN_CO'),\n",
    "                                                (2, 'SM_CO'),\n",
    "                                                (3, 'CAR_CO')]),\n",
    "                          u'headway': dict([(1, 'TRAIN_HE'),\n",
    "                                            (2, 'SM_HE')]),\n",
    "                          u'seat_configuration': dict([(2, \"SM_SEATS\")])}\n",
    "\n",
    "# Specify the availability variables\n",
    "availability_variables = dict(zip(range(1, 4), ['TRAIN_AV', 'SM_AV', 'CAR_AV']))\n",
    "\n",
    "# Determine the columns that will denote the\n",
    "# new column of alternative ids, and the columns\n",
    "# that denote the custom observation ids and the \n",
    "# choice column\n",
    "new_alt_id = \"mode_id\"\n",
    "obs_id_column = \"custom_id\"\n",
    "choice_column = \"CHOICE\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4b. Actually perform the conversion from wide to long formats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>custom_id</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mode_id</th>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHOICE</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GROUP</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SURVEY</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SP</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PURPOSE</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FIRST</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TICKET</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>WHO</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LUGGAGE</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AGE</th>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MALE</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>INCOME</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GA</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ORIGIN</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DEST</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>travel_time</th>\n",
       "      <td>112.0</td>\n",
       "      <td>63.0</td>\n",
       "      <td>117.0</td>\n",
       "      <td>103.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>117.0</td>\n",
       "      <td>130.0</td>\n",
       "      <td>67.0</td>\n",
       "      <td>117.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>seat_configuration</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>headway</th>\n",
       "      <td>120.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>travel_cost</th>\n",
       "      <td>48.0</td>\n",
       "      <td>52.0</td>\n",
       "      <td>65.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>84.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>58.0</td>\n",
       "      <td>52.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        0     1      2      3     4      5      6     7      8\n",
       "custom_id             1.0   1.0    1.0    2.0   2.0    2.0    3.0   3.0    3.0\n",
       "mode_id               1.0   2.0    3.0    1.0   2.0    3.0    1.0   2.0    3.0\n",
       "CHOICE                0.0   1.0    0.0    0.0   1.0    0.0    0.0   1.0    0.0\n",
       "GROUP                 2.0   2.0    2.0    2.0   2.0    2.0    2.0   2.0    2.0\n",
       "SURVEY                0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "SP                    1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "ID                    1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "PURPOSE               1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "FIRST                 0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "TICKET                1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "WHO                   1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "LUGGAGE               0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "AGE                   3.0   3.0    3.0    3.0   3.0    3.0    3.0   3.0    3.0\n",
       "MALE                  0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "INCOME                2.0   2.0    2.0    2.0   2.0    2.0    2.0   2.0    2.0\n",
       "GA                    0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "ORIGIN                2.0   2.0    2.0    2.0   2.0    2.0    2.0   2.0    2.0\n",
       "DEST                  1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "travel_time         112.0  63.0  117.0  103.0  60.0  117.0  130.0  67.0  117.0\n",
       "seat_configuration    0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "headway             120.0  20.0    0.0   30.0  10.0    0.0   60.0  30.0    0.0\n",
       "travel_cost          48.0  52.0   65.0   48.0  49.0   84.0   48.0  58.0   52.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Perform the desired conversion\n",
    "long_swiss_metro = pl.convert_wide_to_long(clean_sm_wide, \n",
    "                                           ind_variables, \n",
    "                                           alt_varying_variables, \n",
    "                                           availability_variables, \n",
    "                                           obs_id_column, \n",
    "                                           choice_column,\n",
    "                                           new_alt_id_name=new_alt_id)\n",
    "\n",
    "# Look at the first 9 rows of the long-format dataframe\n",
    "long_swiss_metro.head(9).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Create the variables used in the Python Biogeme Nested Logit Model Example\n",
    "\n",
    "In 09NestedLogit.py, the travel time and travel cost variables are scaled for ease of numeric optimization. We will do the same such that our estimated coefficients are comparable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/timothyb0912/anaconda/envs/debug/lib/python3.5/site-packages/pandas/core/computation/expressions.py:180: UserWarning: evaluating in Python space because the '*' operator is not supported by numexpr for the bool dtype, use '&' instead\n",
      "  .format(op=op_str, alt_op=unsupported[op_str]))\n"
     ]
    }
   ],
   "source": [
    "# Scale both the travel time and travel cost by 100\n",
    "long_swiss_metro[\"travel_time_hundredth\"] = (long_swiss_metro[\"travel_time\"] /\n",
    "                                             100.0)\n",
    "\n",
    "# Figure out which rows correspond to train or swiss metro \n",
    "# alternatives for individuals with GA passes. These individuals face no \n",
    "# marginal costs for a trip\n",
    "train_pass_train_alt = ((long_swiss_metro[\"GA\"] == 1) *\n",
    "                        (long_swiss_metro[\"mode_id\"].isin([1, 2]))).astype(int)\n",
    "# Note that the (train_pass_train_alt == 0) term accounts for the\n",
    "# fact that those with a GA pass have no marginal cost for the trip\n",
    "long_swiss_metro[\"travel_cost_hundredth\"] = (long_swiss_metro[\"travel_cost\"] *\n",
    "                                             (train_pass_train_alt == 0) /\n",
    "                                             100.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Specify and Estimate the Python Biogeme Nested Logit Model Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6a. Specify the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the nesting values\n",
    "nest_membership = OrderedDict()\n",
    "nest_membership[\"Future Modes\"] = [2]\n",
    "nest_membership[\"Existing Modes\"] = [1, 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the model's specification dictionary and variable names dictionary\n",
    "# NOTE: - Keys should be variables within the long format dataframe.\n",
    "#         The sole exception to this is the \"intercept\" key.\n",
    "#       - For the specification dictionary, the values should be lists\n",
    "#         or lists of lists. Within a list, or within the inner-most\n",
    "#         list should be the alternative ID's of the alternative whose\n",
    "#         utility specification the explanatory variable is entering.\n",
    "\n",
    "example_specification = OrderedDict()\n",
    "example_names = OrderedDict()\n",
    "\n",
    "# Note that 1 is the id for the Train and 3 is the id for the Car.\n",
    "# The next two lines are placing alternative specific constants in\n",
    "# the utility equations for the Train and for the Car. The order\n",
    "# in which these variables are placed is chosen so the summary\n",
    "# dataframe which is returned will match that shown in the HTML\n",
    "# file of the python biogeme example.\n",
    "example_specification[\"intercept\"] = [3, 1]\n",
    "example_names[\"intercept\"] = ['ASC Car', 'ASC Train']\n",
    "\n",
    "# Note that the names used below are simply for consistency with\n",
    "# the coefficient names given in the Python Biogeme example.\n",
    "# example_specification[\"travel_cost_hundredth\"] = [[1, 2, 3]]\n",
    "# example_names[\"travel_cost_hundredth\"] = ['B_COST']\n",
    "\n",
    "example_specification[\"travel_cost_hundredth\"] = [[1, 2, 3]]\n",
    "example_names[\"travel_cost_hundredth\"] = ['B_COST']\n",
    "\n",
    "example_specification[\"travel_time_hundredth\"] = [[1, 2, 3]]\n",
    "example_names[\"travel_time_hundredth\"] = ['B_TIME']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6b. Estimate the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One main difference between the nested logit implementation in PyLogit and in Python Biogeme or mLogit in R is that PyLogit reparameterizes the 'standard' nested logit model. In particular, one standard reperesntation of the nested logit model is in terms of the inverse of the 'scale' parameter for each nest (see for example the representation given by Kenneth Train in section 4.2 <a href=\"http://eml.berkeley.edu/books/choice2nd/Ch04_p76-96.pdf\">here</a>). The 'scale' parameter has domain from zero to infinity, therefore the inverse of the scale parameter has the same domain.\n",
    "\n",
    "However, for econometric purposes (such as conforming to the assumptions that individuals are making choices through a utility maximizing decision protocol), the scale parameter of a 'lower level nest' is constrained to be greater than or equal to 1 (assuming that the 'upper level nest' is constrained to 1.0 for identification purposes). The inverse of the scale parameter would then be constrained to be between 0.0 and 1.0 in this case. In order to make use of unconstrained optimization algorithms, we therefore estimate the logit ( i.e. $\\ln \\left[ \\frac{\\textrm{scale}^{-1}}{1.0 - \\textrm{scale}^{-1}} \\right]$) of the inverse of the scale parameter, assuming that the inverse of the scale parameter will lie between zero and one (and accordingly that the scale parameter be greater than or equal to one)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function that calculates the \"logit\" transformation of values\n",
    "# between 0.0 and 1.0.\n",
    "def logit(x):\n",
    "    \"\"\"\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : int, float, or 1D ndarray.\n",
    "        If an array, all elements should be ints or floats. All\n",
    "        elements should be between zero and one, exclusive of 1.0.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    The logit of x:  `np.log(x / (1.0 - x))`.\n",
    "    \"\"\"\n",
    "    return np.log(x/(1.0 - x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log-likelihood at zero: -6,520.1578\n",
      "Initial Log-likelihood: -5,239.0782\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/timothyb0912/anaconda/envs/debug/lib/python3.5/site-packages/scipy/optimize/_minimize.py:507: RuntimeWarning: Method BFGS does not use Hessian information (hess).\n",
      "  RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation Time for Point Estimation: 22.88 seconds.\n",
      "Final log-likelihood: -5,236.9000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/timothyb0912/anaconda/envs/debug/lib/python3.5/site-packages/pylogit-0.2.2-py3.5.egg/pylogit/base_multinomial_cm_v2.py:1259: RuntimeWarning: invalid value encountered in sqrt\n",
      "/Users/timothyb0912/anaconda/envs/debug/lib/python3.5/site-packages/scipy/stats/_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in greater\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/Users/timothyb0912/anaconda/envs/debug/lib/python3.5/site-packages/scipy/stats/_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in less\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/Users/timothyb0912/anaconda/envs/debug/lib/python3.5/site-packages/scipy/stats/_distn_infrastructure.py:1821: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= self.a)\n"
     ]
    }
   ],
   "source": [
    "# Provide the module with the needed input arguments to create\n",
    "# an instance of the MNL model class\n",
    "example_nested = pl.create_choice_model(data=long_swiss_metro,\n",
    "                                        alt_id_col=new_alt_id,\n",
    "                                        obs_id_col=obs_id_column,\n",
    "                                        choice_col=choice_column,\n",
    "                                        specification=example_specification,\n",
    "                                        model_type=\"Nested Logit\",\n",
    "                                        names=example_names,\n",
    "                                        nest_spec=nest_membership)\n",
    "\n",
    "# Specify the initial nesting parameter values\n",
    "# Note: This should be in terms of the reparameterized values used\n",
    "# by PyLogit.\n",
    "\n",
    "# Note: The '40' corresponds to scale parameter that is numerically\n",
    "# indistinguishable from 1.0\n",
    "\n",
    "# Note: 2.05 is the scale parameter that is estimated by PythonBiogeme\n",
    "# so we invert it, then take the logit of this inverse to get the\n",
    "# corresponding starting value to be used by PyLogit.\n",
    "# Note the first value corresponds to the first nest in 'nest_spec'\n",
    "                       constrained_pos=[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
    "# alternative, and the nest parameter of degenerate nests is not\n",
    "# identified.\n",
    "example_nested.fit_mle(init_values,\n",
    "                       constrained_pos=[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, note that the functionality of using parameter constraints is restriced to the Mixed Logit and Nested Logit models at the moment. Moreover, this functionality is only relevant when using optimization method that make use of gradient information. Gradient-free estimation methods such as 'powell's' method or 'nelder-mead' will not make use of the constrained_pos keyword argument."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.c Compare the model output with that of Python Biogeme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Nested Logit Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>CHOICE</td>       <th>  No. Observations:  </th>    <td>6,768</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>         <td>Nested Logit Model</td> <th>  Df Residuals:      </th>    <td>6,762</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>MLE</td>        <th>  Df Model:          </th>      <td>6</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Fri, 02 Nov 2018</td>  <th>  Pseudo R-squ.:     </th>    <td>0.197</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>15:22:50</td>      <th>  Pseudo R-bar-squ.: </th>    <td>0.196</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AIC:</th>               <td>10,485.800</td>     <th>  Log-Likelihood:    </th> <td>-5,236.900</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>BIC:</th>               <td>10,526.720</td>     <th>  LL-Null:           </th> <td>-6,520.158</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>           <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Future Modes</th>   <td>   40.0000</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Existing Modes</th> <td>   -0.0527</td> <td>    0.082</td> <td>   -0.646</td> <td> 0.519</td> <td>   -0.212</td> <td>    0.107</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ASC Car</th>        <td>   -0.1672</td> <td>    0.032</td> <td>   -5.243</td> <td> 0.000</td> <td>   -0.230</td> <td>   -0.105</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ASC Train</th>      <td>   -0.5119</td> <td>    0.035</td> <td>  -14.781</td> <td> 0.000</td> <td>   -0.580</td> <td>   -0.444</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>B_COST</th>         <td>   -0.8567</td> <td>    0.036</td> <td>  -23.578</td> <td> 0.000</td> <td>   -0.928</td> <td>   -0.785</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>B_TIME</th>         <td>   -0.8987</td> <td>    0.034</td> <td>  -26.228</td> <td> 0.000</td> <td>   -0.966</td> <td>   -0.832</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                    Nested Logit Model Regression Results                     \n",
       "==============================================================================\n",
       "Dep. Variable:                 CHOICE   No. Observations:                6,768\n",
       "Model:             Nested Logit Model   Df Residuals:                    6,762\n",
       "Method:                           MLE   Df Model:                            6\n",
       "Date:                Fri, 02 Nov 2018   Pseudo R-squ.:                   0.197\n",
       "Time:                        15:22:50   Pseudo R-bar-squ.:               0.196\n",
       "AIC:                       10,485.800   Log-Likelihood:             -5,236.900\n",
       "BIC:                       10,526.720   LL-Null:                    -6,520.158\n",
       "==================================================================================\n",
       "                     coef    std err          z      P>|z|      [0.025      0.975]\n",
       "----------------------------------------------------------------------------------\n",
       "Future Modes      40.0000        nan        nan        nan         nan         nan\n",
       "Existing Modes    -0.0527      0.082     -0.646      0.519      -0.212       0.107\n",
       "ASC Car           -0.1672      0.032     -5.243      0.000      -0.230      -0.105\n",
       "ASC Train         -0.5119      0.035    -14.781      0.000      -0.580      -0.444\n",
       "B_COST            -0.8567      0.036    -23.578      0.000      -0.928      -0.785\n",
       "B_TIME            -0.8987      0.034    -26.228      0.000      -0.966      -0.832\n",
       "==================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the estimated coefficients and goodness-of-fit statistics\n",
    "example_nested.get_statsmodels_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare with PythonBiogeme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyLogit's estimated Mu is: 2.0541\n"
     ]
    }
   ],
   "source": [
    "# Note that the Mu (i.e the scale parameter) estimated by python biogeme is \n",
    "# 1.0 / nest_coefficient where\n",
    "# nest_coefficient = 1.0 / (1.0 + exp[-1 * estimated_nest_param])\n",
    "pylogit_mu = 1.0 + np.exp(-1 * example_nested.params[\"Existing Modes\"])\n",
    "print(\"PyLogit's estimated Mu is: {:,.4f}\".format(pylogit_mu))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "My parameter estimates match those of Python Biogeme. <br>\n",
    "The Python Biogeme log-likelihood is -5,236.900 and their estimated parameters are:\n",
    "<pre>\n",
    "ASC Car:    -0.167\n",
    "ASC Train:  -0.512\n",
    "B_COST:     -0.857\n",
    "B_TIME:     -0.899\n",
    "Mu:          2.05\n",
    "</pre>\n",
    "\n",
    "As shown above, my log-likelihood is -5,236.900, and my estimated parameters are:\n",
    "<pre>\n",
    "ASC Car:                   -0.1672\n",
    "ASC Train:                 -0.5119\n",
    "B_COST:                    -0.8567\n",
    "B_TIME:                    -0.8987\n",
    "Existing Modes Nest Param:  2.0541\n",
    "</pre>\n",
    "\n",
    "PyLogit's covariance estimates for the Nested Logit model are currently based on the BHHH approximation to the Fisher Information Matrix. This is the same procedure used by mlogit. However, based on the disaggreement between PyLogit's standard errors and those of Python Biogeme, Python Biogeme is clearly not using the BHHH approximation to the Fisher Information Matrix to calculate its standard errors. How does Python Biogeme calculate its standard errors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigate the use of numeric approximations to the Hessian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ASC Car           0.037136\n",
       "ASC Train         0.045180\n",
       "B_COST            0.046273\n",
       "B_TIME            0.056990\n",
       "Future Modes           NaN\n",
       "Existing Modes    0.111667\n",
       "dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create objects for all of the necessary arguments that are\n",
    "# needed to compute the log-likelihood of the nested logit model\n",
    "# given the data used in this example\n",
    "nested_design = example_nested.design\n",
    "mapping_res = example_nested.get_mappings_for_fit()\n",
    "choice_array = long_swiss_metro[\"CHOICE\"].values\n",
    "\n",
    "# Create a nested logit estimation object\n",
    "est_object = nl.NestedEstimator(example_nested, \n",
    "                                mapping_res, \n",
    "                                None, \n",
    "                                np.zeros(example_nested.params.size),\n",
    "                                nl.split_param_vec,\n",
    "                                constrained_pos=[0])\n",
    "\n",
    "# Create a 'convenience' function that simply returns the log-likelihood\n",
    "# given a vector of coefficients\n",
    "def convenient_log_likelihood(all_coefs):\n",
    "    return est_object.convenience_calc_log_likelihood(all_coefs)\n",
    "\n",
    "# Calculate the numeric hessian\n",
    "numeric_hess = numdiff.approx_hess(example_nested.params.values,\n",
    "                                   convenient_log_likelihood)\n",
    "# Account for the fact that the first param is constrained\n",
    "numeric_hess[0, :] = 0\n",
    "numeric_hess[:, 0] = 0\n",
    "numeric_hess[0, 0] = -1\n",
    "# Calculate the asymptotic covariance with the numeric hessian\n",
    "numeric_cov = -1 * scipy.linalg.inv(numeric_hess)\n",
    "# Get the numeric standard errors\n",
    "numeric_std_errs = pd.Series(np.sqrt(np.diag(numeric_cov)),\n",
    "                             index=example_nested.params.index)\n",
    "# Make sure the Future Modes Nest param has a standard error of np.nan\n",
    "numeric_std_errs.loc[\"Future Modes\"] = np.nan\n",
    "# Order the numeric standard errors according to the Python Biogeme\n",
    "# output\n",
    "numeric_std_errs = pd.concat([numeric_std_errs[example_nested.params.index[2:]],\n",
    "                              numeric_std_errs[example_nested.params.index[:2]]],\n",
    "                             axis=0)\n",
    "# Display the numeric standard errors\n",
    "numeric_std_errs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python Biogeme Output\n",
    "\n",
    "<pre>\n",
    "Name \t   Value\tStd err\tt-test\tp-value\n",
    "ASC_CAR\t   -0.167\t0.0371\t -4.50\t0.00\n",
    "ASC_TRAIN     -0.512\t0.0452\t-11.33\t0.00\n",
    "B_COST\t    -0.857\t0.0463\t-18.51\t0.00\n",
    "B_TIME\t    -0.899\t0.0570\t-15.77\t0.00\n",
    "MU\t         2.05\t 0.118\t  17.45\t0.00\n",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From above, we see that for the index coefficients, the standard errors that are calculated using the numeric approximation of the hessian match the standard errors returned by Python Biogeme. This suggests that the standard errors of Python Biogeme, for the nested logit model, are based on a numeric differentiation approximation to the Hessian.\n",
    "\n",
    "Below, we investigate whether the numeric approximation of the gradient via numeric differentiation is a close approximation to the analytic gradient. The premise is that if the numeric gradient does not adequately approximate the analytic gradient, then what chance does the numeric hessian have of adequately approximating the analytic hessian? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Numeric Differentiation</th>\n",
       "      <th>Analytic</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Future Modes</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Existing Modes</th>\n",
       "      <td>0.000610</td>\n",
       "      <td>-5.840467e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC Car</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-1.857643e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC Train</th>\n",
       "      <td>0.000119</td>\n",
       "      <td>-1.453355e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_COST</th>\n",
       "      <td>0.000071</td>\n",
       "      <td>3.240487e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_TIME</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-1.626236e-08</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Numeric Differentiation      Analytic\n",
       "Future Modes                   0.000000 -0.000000e+00\n",
       "Existing Modes                 0.000610 -5.840467e-09\n",
       "ASC Car                        0.000000 -1.857643e-10\n",
       "ASC Train                      0.000119 -1.453355e-08\n",
       "B_COST                         0.000071  3.240487e-09\n",
       "B_TIME                         0.000000 -1.626236e-08"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Approximate the gradient using numeric differentiation\n",
    "numeric_grad = numdiff.approx_fprime(example_nested.params.values,\n",
    "                                     convenient_log_likelihood)\n",
    "pd.DataFrame([numeric_grad,\n",
    "              example_nested.gradient.values],\n",
    "             index=[\"Numeric Differentiation\", \"Analytic\"],\n",
    "             columns=example_nested.params.index).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the dataframe above, we see that the numeric gradient does not adequately approximate the analytic gradient. The numeric gradient is incorrect by a factor of 2 to 4 depending on what variable is being examined (this excludes the Future Modes Nest Param that constrained to 40).\n",
    "\n",
    "Given that the numeric gradient does not provide a good approximation of the analytic gradient, I do not expect the numeric hessian (nor the BHHH approximation used by PyLogit) to provide a good approximation to the analytic hessian. Since the standard errors calculated from the numeric hessian matches PythonBioeme's standard errors, this suggests that Python Biogeme is calculating its hessian (and therefore standard errors) via numeric differentiation, and as suggested above, this is likely to be poor approximation to the analytic hessian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:debug]",
   "language": "python",
   "name": "conda-env-debug-py"
  },
  "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

<MSG> Updated nested logit example.

<DFF> @@ -754,6 +754,13 @@
     "                       constrained_pos=[0])"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Also, note that the functionality of using parameter constraints is restriced to the Mixed Logit and Nested Logit models at the moment. Moreover, this functionality is only relevant when using optimization method that make use of gradient information. Gradient-free estimation methods such as 'powell's' method or 'nelder-mead' will not make use of the constrained_pos keyword argument."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

Updated nested logit example.