<jupyter_start><jupyter_text>Calibration Estimation of the model parametersquetzal_germany is being estimated using [PandasBiogeme](https://biogeme.epfl.ch/). This notebook estimates calibration parameters for the model's utility functions.- Documentation and reference: [<NAME>. (2020). A short introduction to PandasBiogeme. Technical report TRANSP-OR 200605. Transport and Mobility Laboratory, ENAC, EPFL.](https://transp-or.epfl.ch/documents/technicalReports/Bier20.pdf)- Tutorial: https://www.youtube.com/watch?v=OiM94B8WayA  Model formulationThe model consists of systematic utility functions, one for each mode. They comprise an alternaive-specific constant (ASC), a distance-dependent part with travel time and cost summarised as generalised cost (GC), and the accessibility> V_i = ASC + GC(T, C, b_gc_i) + AC • b_ac_iIndex i marks the demand group. I = {'commuting' (1), 'education' (2), 'shopping/medical' (3), 'business' (4), 'private' (6)}Note: The cost variable already includes subscriptions<jupyter_code>input_path = '../input/'
model_path = '../model/'<jupyter_output><empty_output><jupyter_text>Prepare the database<jupyter_code>df = pd.read_csv(input_path + 'transport_demand/calibration_inter-cellular_trips_MiD2017.csv')
df = df[['cost_rail', 'cost_car', 'cost_coach', 'cost_bus', 'cost_walk', 'cost_air',
         'time_rail', 'time_car', 'time_coach', 'time_bus', 'time_walk', 'time_air',
#         'accessibility_rail', 'accessibility_car', 'accessibility_coach', 'accessibility_bus',
#         'accessibility_walk', 'accessibility_air',
         'mode_model', 'purpose_vp', 'car_avail', 'distance']]
df.columns = ['C_RAIL', 'C_CAR', 'C_COACH', 'C_BUS', 'C_NON_MOTOR', 'C_AIR',
              'T_RAIL', 'T_CAR', 'T_COACH', 'T_BUS', 'T_NON_MOTOR', 'T_AIR',
#              'AC_RAIL', 'AC_CAR', 'AC_COACH', 'AC_BUS', 'AC_NON_MOTOR', 'AC_AIR',
              'MODE', 'PURPOSE', 'CAR_AV', 'DIST']
inf = 1000
df = df.replace({np.inf:inf})
df.describe()
# Scale time to hours
df['T_RAIL_S'] = df['T_RAIL'] / 60
df['T_CAR_S'] = df['T_CAR'] / 60
df['T_COACH_S'] = df['T_COACH'] / 60
df['T_BUS_S'] = df['T_BUS'] / 60
df['T_AIR_S'] = df['T_AIR'] / 60
df['T_NON_MOTOR_S'] = df['T_NON_MOTOR'] / 60
# Make car availability binary
df['CAR_AV'] = df['CAR_AV'].replace({9:0})
# Remove trips where mode is car but the car availability is zero
# because it irritates the MLE algorithm
mask = ((df['MODE']==6) & (df['CAR_AV']==0))
print('Share of car trips dropped: {}'.format(
    len(df.loc[mask])/len(df.loc[df['MODE']==6])))
df = df.loc[~mask]
# Remove trips where cost are infinity
# Share of drops per mode must be equal, otherwise the calibration is skewed
# Ignore air trips because this dataset has too few observations anyways
max_drop_ratio = 0
lengths = []
for mode, col in zip(range(1,5), ['C_RAIL', 'C_RAIL', 'C_COACH', 'C_BUS']):
    drops = df.loc[((df['MODE']==mode) & (df[col]==inf))].index
    lengths = lengths + list(df.loc[drops, 'DIST'])
    if len(drops) > 0: print('mode ' + str(mode) + ': ' + str(len(drops)) + ' drops')
    max_drop_ratio = max(len(drops) / len(df.loc[df['MODE']==mode]), max_drop_ratio)
print('max_drop_ratio: ' + str(max_drop_ratio))
if max_drop_ratio > 0:
    print('Dropped trips length mean: {}; min: {}; max: {}'.format(
        sum(lengths)/len(lengths), min(lengths), max(lengths)))
# Drop trips with infinite cost
if max_drop_ratio > 0:
    df = df.loc[~((df['MODE']==2) & (df['C_RAIL']==inf))]
    # Don't drop air trips because there are only a few
    for m in [1, 3, 4, 6]:
        n_drops = int(max_drop_ratio * len(df.loc[df['MODE']==m]))
        df = df.drop(df.loc[(df['MODE']==m) & (df['DIST']<max(lengths)) &
                            (df['DIST']>min(lengths))].sample(n_drops).index)
        print('mode ' + str(m) + ': ' + str(n_drops) + ' drops')
    print('New number of observations: ' + str(len(df.index)))
# Merge rail short and long because people percieve it the same
df['MODE'] = df['MODE'].replace({2:1})<jupyter_output><empty_output><jupyter_text>Model variablesAll columns are variables. DefineVariable creates a new column in the database.  Cost dampingMany modelling studies have shown that cost damping is required in order to flatten the tail of time/cost elasticities, i.e. decrease the impact of long distances on choice results to prevent from overestimation of time/cost parameters. Cost damping represents the property of decreasing marginal utility. It is commonly approached with Box-Cox transformations of generalised cost (usually defined as the sum of travel time and travel expenditures divided by the value of time). Daly (2010) proposes a hybrid function as a sum of the linear term and the narural logarithm of the same. However, the linear term still dominates cost on long distances. Rich (2015), main developer of the Danish National Transport Model, proposes a more complex spline function which successfully manages cost damping and even outperforms the Box-Cox transformation in terms of stability of elasticities.<jupyter_code># Using the hybrid linear-logarithmic version requires a new column
# for the logarithm of time in minutes
'''df['T_RAIL_S'] = np.log(df['T_RAIL'])
df['T_CAR_S'] = np.log(df['T_CAR'])
df['T_COACH_S'] = np.log(df['T_COACH'])
df['T_BUS_S'] = np.log(df['T_BUS'])
df['T_AIR_S'] = np.log(df['T_AIR'])
df['T_NON_MOTOR_S'] = np.log(df['T_NON_MOTOR'])'''
print('This functional form performes poorly')<jupyter_output>This functional form performes poorly<jupyter_text>Generalised costBoth, travel time and monetary cost should be affected by cost damping measures. It is logical to define a generalised cost term `GC` with dimension of time units. This requires definition or estimation of values of time, in order to rescale monetary units, for all segments. Usually, the value of time (VoT) is distance-dependent. In the case of Germany, VoT can be taken from research undertaken within the Federal Government's transport study "Bundesverkehrswegeplan 2030": Axhausen et al. 2015. Ermittlung von Bewertungsansätzen für Reisezeiten und Zuverlässigkeit auf der Basis eines Modells für modale Verlagerungen im nicht-gewerblichen und gewerblichen Personenverkehr für die Bundesverkehrswegeplanung<jupyter_code># VoT from literature, distance-dependent, see cal19
VoT = pd.read_csv(input_path + 'vot.csv', header=[0,1], index_col=0)
VoT.sample(2)
# Make distance integer
df['DIST'] = df['DIST'].apply(int)
# Generate generalised cost
VoT = VoT.to_dict()
df['GC_RAIL'] = df['T_RAIL_S'] + [c / VoT[('Fz'+str(p), 'PT')][d]
                                  for c,d,p in zip(df['C_RAIL'], df['DIST'], df['PURPOSE'])]
df['GC_COACH'] = df['T_COACH_S'] + [c / VoT[('Fz'+str(p), 'PT')][d]
                                    for c,d,p in zip(df['C_COACH'], df['DIST'], df['PURPOSE'])]
df['GC_BUS'] = df['T_BUS_S'] + [c / VoT[('Fz'+str(p), 'PT')][d]
                                for c,d,p in zip(df['C_BUS'], df['DIST'], df['PURPOSE'])]
df['GC_AIR'] = df['T_AIR_S'] + [c / VoT[('Fz'+str(p), 'air')][d]
                                for c,d,p in zip(df['C_AIR'], df['DIST'], df['PURPOSE'])]
df['GC_CAR'] = df['T_CAR_S'] + [c / VoT[('Fz'+str(p), 'car')][d]
                                for c,d,p in zip(df['C_CAR'], df['DIST'], df['PURPOSE'])]
df['GC_NON_MOTOR'] = df['T_NON_MOTOR_S']
# Create the initial database and make columns global variables
database = db.Database('MiD2017', df.copy())
globals().update(database.variables)
database.getSampleSize()<jupyter_output><empty_output><jupyter_text>Estimation parameters<jupyter_code>asc_rail = ex.Beta('asc_rail', 0, None, None, 0)
asc_coach = ex.Beta('asc_coach', 0, None, None, 0)
asc_bus = ex.Beta('asc_bus', 0, None, None, 0)
asc_air = ex.Beta('asc_air', 0, None, None, 0)
asc_car = ex.Beta('asc_car', 0, None, None, 1)
asc_non_motor = ex.Beta('asc_non_motor', 0, None, None, 0)
b_t = ex.Beta('b_t', 0, None, None, 0)
b_c = ex.Beta('b_c', 0, None, None, 0)
b_ac = ex.Beta('b_ac', 0, None, None, 0)
# non-linear time component
b_t2 = ex.Beta('b_t2', 0, None, None, 0)
# generalised cost function parameter
b_gc = ex.Beta('b_gc', 0, None, None, 0)
# Parameters for the nested logit structure
mu_pt = ex.Beta('mu_pt', 1, 1, 10, 0)<jupyter_output><empty_output><jupyter_text>Utility functions<jupyter_code># Aggregated formulation
V_RAIL = asc_rail + b_t * T_RAIL_S + b_c * C_RAIL + b_ac * AC_RAIL
V_COACH = asc_coach + b_t * T_COACH_S + b_c * C_COACH + b_ac * AC_COACH
V_BUS = asc_bus + b_t * T_BUS_S + b_c * C_BUS + b_ac * AC_BUS
V_AIR = asc_air + b_t * T_AIR_S + b_c * C_AIR + b_ac * AC_AIR
V_CAR = asc_car + b_t * T_CAR_S + b_c * C_CAR
V_NON_MOTOR = asc_non_motor + b_t * T_NON_MOTOR_S
# Aggregated formulation with non-linear perception of travel time
V_RAIL = asc_rail + b_t * T_RAIL + b_t2 * T_RAIL_S + b_c * C_RAIL + b_ac * AC_RAIL
V_COACH = asc_coach + b_t * T_COACH + b_t2 * T_COACH_S + b_c * C_COACH + b_ac * AC_COACH
V_BUS = asc_bus + b_t * T_BUS + b_t2 * T_BUS_S + b_c * C_BUS + b_ac * AC_BUS
V_AIR = asc_air + b_t * T_AIR + b_t2 * T_AIR_S + b_c * C_AIR + b_ac * AC_AIR
V_CAR = asc_car + b_t * T_CAR + b_t2 * T_CAR_S + b_c * C_CAR
V_NON_MOTOR = asc_non_motor + b_t * T_NON_MOTOR + b_t2 * T_NON_MOTOR_S
# The cost damping function
def spline(x, beta=b_gc, Q=3, c0=0, c1=20, c2=40, c3=np.inf):
    alpha = [0, -beta/2*ex.Power(ex.log(c1),3),
             -beta/2*ex.log(c1)*(3*ex.Power(ex.log(c2),2)+ex.Power(ex.log(c1),2))] # for Q=3
    theta = [1, 3/2*ex.log(c1), 3*ex.log(c1)*ex.log(c2)] # for Q=3
    def component(x, q):
        return beta*theta[q-1]*ex.Power(ex.log(x),Q-q+1) + alpha[q-1]
    if x < c1:
        return component(x,1)
    elif x < c2:
        return component(x,2)
    elif x >= c2:
        return component(x,3)
    #return sum([component(x, q) if x>=c[q-1] and x<c[q] else 0 for q in range(1,Q+1)])
# Aggregated formulation with damped generalised cost
V_RAIL = asc_rail + spline(GC_RAIL)# + b_ac * AC_RAIL
V_COACH = asc_coach + spline(GC_COACH)# + b_ac * AC_COACH
V_BUS = asc_bus + spline(GC_BUS)# + b_ac * AC_BUS
V_AIR = asc_air + spline(GC_AIR)# + b_ac * AC_AIR
V_CAR = asc_car + spline(GC_CAR)
V_NON_MOTOR = asc_non_motor + spline(GC_NON_MOTOR)<jupyter_output><empty_output><jupyter_text>Run the estimation<jupyter_code># Define level of verbosity
logger = message.bioMessage()
#logger.setSilent()
logger.setWarning()
#logger.setGeneral()
#logger.setDetailed()
# Map modes to utility functions
V = {1:V_RAIL,
     2:V_RAIL,
     3:V_COACH,
     4:V_BUS,
     5:V_AIR,
     6:V_CAR,
     7:V_NON_MOTOR}
# Map the availability of alternatives with MODE as key
# PT is always available
av = {1:1,
      2:1,
      3:1,
      4:1,
      5:1,
      6:CAR_AV,
      7:1}
# Mode nests as tuples with nest name and dictionary where
# alternative IDs are mapped to alpha values. Missing ID's alpha is zero
# Alternatively use lists with mode ID without alpha
pt = mu_pt, [1, 2, 3, 4]
air = 1, [5]
car = 1, [6]
walk = 1, [7]
nests = pt, air, car, walk
# Choose the multinomial logit model
mnl = models.loglogit(V, av, MODE)
model_mnl = bio.BIOGEME(database, mnl)
model_mnl.modelName = 'MNL'
results_mnl = model_mnl.estimate()
results = results_mnl.getEstimatedParameters()
for key, val in results_mnl.getGeneralStatistics().items():
    results.loc[key] = [val[0], val[1]] + ['' for i in range(len(results.columns)-2)]
results
# Write results to a file
writer = pd.ExcelWriter(input_path + 'estimation_results_splineGC_distVoT_knots2040_2.xlsx', engine='xlsxwriter')
results.to_excel(writer, sheet_name=model_mnl.modelName)
# Choose the logarithmic nested logit model
nl = models.lognested(V, av, nests, MODE)
# Nested Logit
model_nl = bio.BIOGEME(database, nl)
model_nl.modelName = 'NL'
results_nl = model_nl.estimate()
results = results_nl.getEstimatedParameters()
for key, val in results_nl.getGeneralStatistics().items():
    results.loc[key] = [val[0], val[1]] + ['' for i in range(len(results.columns)-2)]
results
results.to_excel(writer, sheet_name=model_nl.modelName)
# Run the estimation by purpose
results = []
for p in [1,2,3,4,6]:
    database = db.Database('MiD2017', df.copy())
    database.remove(PURPOSE!=p)
    print('Sample size for purpose {}: {}'.format(p, database.getSampleSize()))
    model = bio.BIOGEME(database, nl) # Choose the model formulation
    model.modelName = 'NL_Fz' + str(p) # Name it
    results.append(model.estimate()) # Estimation
    output = results[-1].getEstimatedParameters()
    # Add results to the Excel file
    for key, val in results[-1].getGeneralStatistics().items():
        output.loc[key] = [val[0], val[1]] + ['' for i in range(len(output.columns)-2)]
    output.to_excel(writer, sheet_name=model.modelName)
writer.save()
# Generate LaTeX code
file = open(input_path + 'estimation_results_LaTeX_code.txt', 'w')
for r in results:
    file.write(r.getLaTeX())
file.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1477794197713308, 'fraction_numerical': 0.016877965310725676, 'mean_word_length': 4.551813471502591, 'pattern_counts': {'":': 1, '<': 33, '<?xml version=': 0, '>': 36, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '312457', 'license': 'Apache-2.0', 'n_tokens_mistral': 5005, 'n_tokens_neox': 4784, 'n_words': 1261, 'path': 'quetzal_germany-master/notebooks/cal20_estimation.ipynb', 'repo_name': 'avi-shake/Flensburg_Modeling', 'size': '12863'}