######################################################################
#                                                                    #
# 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,
               LogsumIntra,
               LogsumInter
               )
{

  ## 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.PT1   = 2
  V.NM1   = 2
  V.PT2   = 2
  V.NM2   = 2
  V.Intra = 2
  V.Inter = 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.PT1 = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
  LS.NM1 = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))
  LS.PT2 = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
  LS.NM2 = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))

  LS.Intra = log(exp(LS.PT1/LogsumIntra) + exp(LS.NM1/LogsumIntra))
  LS.Inter = log(exp(LS.PT2/LogsumInter) + exp(LS.NM2/LogsumInter))

  ## Denominators
  DenoPT1  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
  DenoNM1  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)

  DenoPT2  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
  DenoNM2  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)

  DenoIntra = exp(Vcar) + exp(V.PT1 + LogsumPT*LS.PT1) + exp(V.NM1 + LogsumNM*LS.NM1)
  DenoInter = exp(Vcar) + exp(V.PT2 + LogsumPT*LS.PT2) + exp(V.NM2 + LogsumNM*LS.NM2)

  Deno = exp(V.Intra + LogsumIntra*LS.Intra) + exp(V.Inter + LogsumInter*LS.Inter)

  ## Probabilites
  # Conditional probabilits
  Ptrain.PT1  = exp(Vtrain/LogsumPT)/DenoPT1
  Pbus.PT1    = exp(Vbus/LogsumPT)/DenoPT1
  Pbike.NM1   = exp(Vbike/LogsumNM)/DenoNM1
  Pwalk.NM1   = exp(Vwalk/LogsumNM)/DenoNM1

  Ptrain.PT2  = exp(Vtrain/LogsumPT)/DenoPT2
  Pbus.PT2    = exp(Vbus/LogsumPT)/DenoPT2
  Pbike.NM2   = exp(Vbike/LogsumNM)/DenoNM2
  Pwalk.NM2   = exp(Vwalk/LogsumNM)/DenoNM2

  Pcar.Intra  = exp(Vcar)/DenoIntra
  P.PT1       = exp(V.PT1 + LogsumPT*LS.PT1)/DenoIntra
  P.NM1       = exp(V.NM1 + LogsumNM*LS.NM1)/DenoIntra

  Pcar.Inter  = exp(Vcar)/DenoInter
  P.PT2       = exp(V.PT2 + LogsumPT*LS.PT2)/DenoInter
  P.NM2       = exp(V.NM2 + LogsumNM*LS.NM2)/DenoInter

  P.Intra     = exp(V.Intra + LogsumIntra*LS.Intra)/Deno
  P.Inter     = exp(V.Inter + LogsumInter*LS.Inter)/Deno

  ##Joint probabilities
  Ptrain1     = Ptrain.PT1*P.PT1*P.Intra
  Pbus1       = Pbus.PT1*P.PT1*P.Intra
  Pbike1      = Pbike.NM1*P.NM1*P.Intra
  Pwalk1      = Pwalk.NM1*P.NM1*P.Intra
  Pcar1       = Pcar.Intra*P.Intra

  Ptrain2     = Ptrain.PT2*P.PT2*P.Inter
  Pbus2       = Pbus.PT2*P.PT2*P.Inter
  Pbike2      = Pbike.NM2*P.NM2*P.Inter
  Pwalk2      = Pwalk.NM2*P.NM2*P.Inter
  Pcar2       = Pcar.Inter*P.Inter

  Ptrain1 = (Ptrain1 != 0)*Ptrain1 + (Ptrain1 == 0)
  Pbus1   = (Pbus1   != 0)*Pbus1   + (Pbus1   == 0)
  Pcar1   = (Pcar1   != 0)*Pcar1   + (Pcar1   == 0)
  Pbike1  = (Pbike1  != 0)*Pbike1  + (Pbike1  == 0)
  Pwalk1  = (Pwalk1  != 0)*Pwalk1  + (Pwalk1  == 0)

  Ptrain2 = (Ptrain2 != 0)*Ptrain2 + (Ptrain2 == 0)
  Pbus2   = (Pbus2   != 0)*Pbus2   + (Pbus2   == 0)
  Pcar2   = (Pcar2   != 0)*Pcar2   + (Pcar2   == 0)
  Pbike2  = (Pbike2  != 0)*Pbike2  + (Pbike2  == 0)
  Pwalk2  = (Pwalk2  != 0)*Pwalk2  + (Pwalk2  == 0)


  ## Dummies for 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")
  DIntra   = ((CBD_Data$DistrictCodeO == 103&CBD_Data$DistrictCodeD == 103)|(CBD_Data$DistrictCodeO == 104&CBD_Data$DistrictCodeD == 104))
  DInter   = !DIntra

  ## Calculate the Log-likelihood function
  LL = -colSums(Dtrain*log(Ptrain1) + Dbus*log(Pbus1) + Dcar*log(Pcar1) + Dbike*log(Pbike1) + Dwalk*log(Pwalk1) +
                 Dtrain*log(Ptrain2) + Dbus*log(Pbus2) + Dcar*log(Pcar2) + Dbike*log(Pbike2) + Dwalk*log(Pwalk2)
                 )

}

##### 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,
                  LogsumNM = 1,
                  LogsumPT = 1,
                  LogsumInter = 1,
                  LogsumIntra = 1
                  )


# Clear the Console in R
cat("\014")

## Estimate parameters
NL = mle2(NL_LL, start = NL_initValue,  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)


