######################################################################
#                                                                    #
# Name  :  15th Travel Behavior and Network Model Course             #
# Author:  Team C                                                    #
#          (Canh Do, Hong Nguyen, Aktar Hossain, Perez Barbosa David)#
#          Transportation Study Group (TSG),Hiroshima University     #
# Time  :  23th-25th September, 2016                                 #
#                                                                    #
######################################################################

#Run MNL Code
rm(list = ls())
source("BIN2016-TeamC-MNL.R")

##### Define the log-likelihood function of the logit model#####

NL_LL <- function(Ctrain,
                 Cbus,
                 Ccar,
                 Cbike,
                 TTime,
                 TCost,
                 WDayTrain,
                 WDayBus,
                 WDayCar,
                 WDayBike,
                 WDayWalk,
                 WPurTrain,
                 WPurBus,
                 WPurCar,
                 WPurBike,
                 WPurWalk,
                 LogsumPT,
                 LogsumNM
                 )
{

  ## Declare the log-likelihood variable, set value to 0
  NL_LL = 0

  ## Utility functions
  Vtrain  = Ctrain*matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeTrain  + TCost *CBD_Data$FareTrain                        + WDayTrain *CBD_Data$DWeekday + WPurTrain *CBD_Data$DWPur
  Vbus    = Cbus  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeBus    + TCost *CBD_Data$FareBus                          + WDayBus *CBD_Data$DWeekday   + WPurBus *CBD_Data$DWPur
  Vcar    = Ccar  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeCar         + TCost *(CBD_Data$ParkingCost + CBD_Data$FareCar) + WDayCar *CBD_Data$DWeekday   + WPurCar *CBD_Data$DWPur
  Vbike   = Cbike *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeBike                                                           + WDayBike *CBD_Data$DWeekday  + WPurBike *CBD_Data$DWPur
  Vwalk   =        matrix(0,nrow = N, ncol = 1) + TTime *CBD_Data$TimeWalk                                                           + WDayWalk *CBD_Data$DWeekday  + WPurWalk *CBD_Data$DWPur
  V.PT    = 2
  V.NM    = 2

  Vtrain = CBD_Data$ModeAvailableTrain * Vtrain
  Vbus   = CBD_Data$ModeAvailableBus * Vbus
  Vcar   = CBD_Data$ModeAvailableCar * Vcar
  Vbike  = CBD_Data$ModeAvailableBike * Vbike
  Vwalk  = CBD_Data$ModeAvailableWalk * Vwalk

  ## Log-sum variables
  LS.PT = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
  LS.NM = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))

  ## Denominators
  DenoPT  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
  DenoNM  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)

  Deno = exp(Vcar) + exp(V.PT + LogsumPT*LS.PT) + exp(V.NM + LogsumNM*LS.NM)

  ## Probabilites
  # Conditional probabilits
  Ptrain.PT  = exp(Vtrain/LogsumPT)/DenoPT
  Pbus.PT    = exp(Vbus  /LogsumPT)/DenoPT
  Pbike.NM   = exp(Vbike /LogsumNM)/DenoNM
  Pwalk.NM   = exp(Vwalk /LogsumNM)/DenoNM

  P.PT       = exp(V.PT + LogsumPT*LS.PT)/Deno
  P.NM       = exp(V.NM + LogsumNM*LS.NM)/Deno

  Pcar       = exp(Vcar)/Deno

  ##Joint probabilities
  Ptrain     = Ptrain.PT*P.PT
  Pbus       = Pbus.PT*P.PT
  Pbike      = Pbike.NM*P.NM
  Pwalk      = Pwalk.NM*P.NM

  ## Avoid problems stemming from choice probabilities becoming zero.
  Ptrain = (Ptrain != 0)*Ptrain + (Ptrain == 0)
  Pbus   = (Pbus   != 0)*Pbus   + (Pbus   == 0)
  Pcar   = (Pcar   != 0)*Pcar   + (Pcar   == 0)
  Pbike  = (Pbike  != 0)*Pbike  + (Pbike  == 0)
  Pwalk  = (Pwalk  != 0)*Pwalk  + (Pwalk  == 0)

  ## Dummies for mode choice
  Dtrain   = (CBD_Data$MainModeENG == "Rail")
  Dbus     = (CBD_Data$MainModeENG == "Bus")
  Dcar     = (CBD_Data$MainModeENG == "Car")
  Dbike    = (CBD_Data$MainModeENG == "Bicycle")
  Dwalk    = (CBD_Data$MainModeENG == "Walk")
  DPT      = (CBD_Data$MainModeENG == "Rail" | CBD_Data$MainModeENG == "Bus")
  DNM      = (CBD_Data$MainModeENG == "Bicycle" | CBD_Data$MainModeENG == "Walk")

  ## Calculate the Log-likelihood function
  NL_LL = -colSums(Dtrain*log(Ptrain) + Dbus*log(Pbus) + Dcar*log(Pcar) + Dbike*log(Pbike) + Dwalk*log(Pwalk)
                            )
  return(NL_LL)
}

##### Maximize the Log-likelihood function#####

## Setting up initValaue

NL_initValue = list(Ctrain = 0,
                  Cbus   = 0,
                  Ccar   = 0,
                  Cbike  = 0,
                  TTime  = 0,
                  TCost  = 0,
                  WDayTrain = 0,
                  WDayBus = 0,
                  WDayCar = 0,
                  WDayBike = 0,
                  WDayWalk = 0,
                  WPurTrain = 0,
                  WPurBus = 0,
                  WPurCar = 0,
                  WPurBike = 0,
                  WPurWalk = 0,
                  LogsumPT = 1,
                  LogsumNM = 1)



## Estimate parameters
NL = mle2(NL_LL, start = NL_initValue, method = "BFGS", control = list(trace = T, maxit = 100000000))

NL_LL0 = N*log(1/5)     ## Initial log-likelihood
NL_LL1 = logLik(NL)             ## Converged log-likelihood

## McFadden's Rho square
NL_Rho     = (1 - (NL_LL1/NL_LL0))
NL_Rho.adj = (1 - ((NL_LL1 - length(coef(NL)))/NL_LL0))
NL_GoF = cbind(NL_LL0,NL_LL1, NL_Rho, NL_Rho.adj)  ## Goodness of fit

## Preview NL model result in console
print(summary(NL))
print(NL_GoF)


