######################################################################
#                                                                    #
# 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                                 #
#                                                                    #
######################################################################

#Load library
library(bbmle)
library(xlsx)

### Create the CBD_Data array from the existing CBD_Dataset
rm(list = ls())
Yokohama_Data = read.csv("Data_Clean_English.csv", header = T)

CBD_Data = Yokohama_Data[ which((Yokohama_Data$DistrictCodeO == 103 | Yokohama_Data$DistrictCodeO == 104) & (Yokohama_Data$DistrictCodeD  == 103 | Yokohama_Data$DistrictCodeD  == 104)),]
rm(Yokohama_Data)

## Count the number of CBD_Data rows
N = nrow(CBD_Data)

## CBD_Data cleaning
# Calculating travel cost for car with fuel price and fuel consumption of spassenger cars
CBD_Data$FareCar = 120/15.1*(CBD_Data$DistanceCar/1000)

# Dummy variable for work-related purpose
CBD_Data$DWPur = ifelse(CBD_Data$PurposeCode < 400 & CBD_Data$PurposeCode != 200, 1,0)
#Data$WPur_Dummy = (CBD_Data$PurposeCode < 400 & CBD_Data$PurposeCode != 200)

# Dummy variable for destination choice
CBD_Data$DesChoice = ifelse(CBD_Data$DistrictCodeO == 103 & CBD_Data$DistrictCodeD == 103,1,
                            ifelse(CBD_Data$DistrictCodeO == 104 & CBD_Data$DistrictCodeD == 104,2,
                                   ifelse(CBD_Data$DistrictCodeO == 103 & CBD_Data$DistrictCodeD == 104,3,4)))
## Convert to Day of week
CBD_Data$DayOfWeek = weekdays(as.Date(CBD_Data$ArrivalTime,origin = "1899-12-30"), abbreviate = F)

## Dummy for weekdays
CBD_Data$DWeekday = !(CBD_Data$DayOfWeek == "Saturday" | CBD_Data$DayOfWeek == "Sunday")

## Parking cost for an hour
CBD_Data$ParkingCost = ifelse(CBD_Data$DistrictCodeD == 103, ifelse(CBD_Data$DWeekday == T,1350,1850), ifelse(CBD_Data$DWeekday == T,1400, 1700))

## Modify CBD_Data
CBD_Data$TotalTimeTrain = CBD_Data$TotalTimeTrain/10
CBD_Data$TotalTimeBus   = CBD_Data$TotalTimeBus/10
CBD_Data$TimeCar        = CBD_Data$TimeCar/10
CBD_Data$TimeBike       = CBD_Data$TimeBike/10
CBD_Data$TimeWalk       = CBD_Data$TimeWalk/10

######## MNL Model #######

# Maximizeing the Log-likelihood function
MNL_LL <- function(Ctrain,
                   Cbus,
                   Ccar,
                   Cbike,
                   TTime,
                   TCost,
                   WDayTrain,
                   WDayBus,
                   WDayCar,
                   WDayBike,
                   WDayWalk,
                   WPurTrain,
                   WPurBus,
                   WPurCar,
                   WPurBike,
                   WPurWalk
                   )
{
  ## Calculate the utility function

  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

  ### Calculate the choice probabilities
  ## calculate the Inclusive Value (the denominator of the choice probabilities equation)
  deno = (CBD_Data$ModeAvailableTrain * exp(Vtrain) + CBD_Data$ModeAvailableBus * exp(Vbus) + CBD_Data$ModeAvailableCar * exp(Vcar) + CBD_Data$ModeAvailableBike  * exp(Vbike) + CBD_Data$ModeAvailableWalk  * exp(Vwalk))

  ## Calculate indiviudal choice probabilities
  Ptrain = (CBD_Data$ModeAvailableTrain * exp(Vtrain) / deno)
  Pbus   = (CBD_Data$ModeAvailableBus   * exp(Vbus)   / deno)
  Pcar   = (CBD_Data$ModeAvailableCar   * exp(Vcar)   / deno)
  Pbike  = (CBD_Data$ModeAvailableBike  * exp(Vbike)  / deno)
  Pwalk  = (CBD_Data$ModeAvailableWalk  * exp(Vwalk)  / deno)

  ## 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)

  ## Dummy 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")

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

## Setting up initValaue
#initValue = vector("list", 16)

MNL_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
                  )

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

## Estimate parameters
MNL = mle2(MNL_LL, start = MNL_initValue, method = "BFGS", control = list(trace = T, maxit = 10000))

## Loglikelihood
MNL_LL0 = N*log(1/5)     ## Initial log-likelihood
MNL_LL1 = logLik(MNL)             ## Converged log-likelihood

## McFadden's Rho square
MNL_Rho     = 1 - (MNL_LL1/MNL_LL0)
MNL_Rho.adj = 1 - ((MNL_LL1 - length(coef(MNL)))/MNL_LL0)

## Goodness of fit
MNL_GoF = cbind(MNL_LL0, MNL_LL1, MNL_Rho, MNL_Rho.adj)

## Preview MNL model result in console
print(summary(MNL))
print(MNL_GoF)



