﻿## Multinomial Logit model estimation

### Create the data array from the existing dataset
Data <- read.csv("C:/Data_Clean_English.csv",header=TRUE)
## Count the number of data rows
hh <- nrow(Data)

## Set the initial values of the parameters (the number in parenthesis corresponds to the number of parameters to be estimated)
b0 <- numeric(5)

## Recall that we estimate parameters by maximizing the log-likelihood function, hence we first need to define that function

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

fr <- function(x) {
  ### declare the parameters###
  ## Alternative specific constants
  b1 <- x[1]
  b2 <- x[2]
  b3 <- x[3]
  b4 <- x[4]
  
  ## Travel time to destination
  d1 <- x[5]
  
  ## declare the log-likelihood variable, set value to 0
  LL = 0
  
  ### For this choice problem, we consider the following 5 modes???
  ## 鉄道(train)
  ## バス(bus)
  ## 自動車(car)
  ## 自転車(bike)
  ## 徒歩(walk)
  
  ## calculate the utility function: :introduce the desired explanatory variables in the function
                                        # time           # fare      # constant
  train  <- Data$ModeAvailableTrain*exp(d1*Data$TotalTimeTrain/100  +b1*matrix(1,nrow =hh,ncol=1))
  bus    <- Data$ModeAvailableBus  *exp(d1*Data$TotalTimeBus/100    +b2*matrix(1,nrow =hh,ncol=1))
  car    <- Data$ModeAvailableCar  *exp(d1*Data$TimeCar/100         +b3*matrix(1,nrow =hh,ncol=1))
  bike   <- Data$ModeAvailableBike *exp(d1*Data$TimeBike/100        +b4*matrix(1,nrow =hh,ncol=1))
  walk   <- Data$ModeAvailableWalk *exp(d1*Data$TimeWalk/100                                                     )
  
  ### Calculate the choice probabilities
  ## calculate the Inclusive Value (the denominator of the choice probabilities equation)
  deno <- (car + train + bus + bike + walk)
  
  ## Calculate indiviudal choice probabilities
  Ptrain <- Data$ModeAvailableTrain*(train / deno)
  Pbus   <- Data$ModeAvailableTrain  *(bus   / deno)
  Pcar   <- Data$ModeAvailableCar  *(car   / deno)
  Pbike  <- Data$ModeAvailableBike *(bike  / deno)
  Pwalk  <- Data$ModeAvailableWalk *(walk  / deno)
  
  ## Avoid problems stemming from choice probabilities becoming zero.選択確??????0になってしまった???合に起こる問題???回避
  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)
  
  
  
  ## Choice results
  Ctrain   <- Data$MainModeENG =="Rail"
  Cbus     <- Data$MainModeENG =="Bus"
  Ccar     <- Data$MainModeENG =="Car"
  Cbike    <- Data$MainModeENG =="Bicycle"
  Cwalk    <- Data$MainModeENG =="Walk"
  
  ## Calculate the Log-likelihood function
  LL <- colSums(Ctrain*log(Ptrain) + Cbus*log(Pbus) +
                  Ccar  *log(Pcar)   + Cbike  *log(Pbike) +Cwalk *log(Pwalk))
  
}

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

##Parameter optimization 
res <- optim(b0,fr,gr=NULL, method = "Nelder-Mead", hessian = TRUE, control=list(fnscale=-1))

## Parmeter estimation、Hessian matrix calculation 
b   <- res$par
hhh <- res$hessian

## Calculate the t-statistic
tval <- b/sqrt(-diag(solve(hhh)))

## L(0), Log-Likelihood when all parameters are 0
L0 <- fr(b0)
## LL, maximium likelihood
LL <- res$value

##### Output #####
print(res)
## L(0)
print(L0)
## LL
print(LL)
##rho-square
print((L0-LL)/L0)
## adjusted rho-square
print((L0-(LL-length(b)))/L0)
##estimated parameter values
print(b)
## t-statistic 
print(tval)
#utility equations are needed to be written first followed by prob. formulas
 train  <- Data$ModeAvailableTrain*exp(b[5]*Data$TotalTimeTrain/100  +b[1]*matrix(1,nrow =hh,ncol=1))
   bus    <- Data$ModeAvailableBus  *exp(b[5]*Data$TotalTimeBus/100    +b[2]*matrix(1,nrow =hh,ncol=1))
   car    <- Data$ModeAvailableCar  *exp(b[5]*Data$TimeCar/100         +b[3]*matrix(1,nrow =hh,ncol=1))
   bike   <- Data$ModeAvailableBike *exp(b[5]*Data$TimeBike/100        +b[4]*matrix(1,nrow =hh,ncol=1))
  walk   <- Data$ModeAvailableWalk *exp(b[5]*Data$TimeWalk/100+b[6]*Data$Age/10 )
deno <- (car + train + bus + bike + walk)
Ptrain <- Data$ModeAvailableTrain*(train / deno)
  Pbus   <- Data$ModeAvailableBus  *(bus   / deno)
   Pcar   <- Data$ModeAvailableCar  *(car   / deno)
   Pbike  <- Data$ModeAvailableBike *(bike  / deno)
 Pwalk  <- Data$ModeAvailableWalk *(walk  / deno)
 A=matrix(c(Ptrain,Pbus,Pcar,Pbike,Pwalk),nrow=1522,ncol=5,byrow=FALSE)
Q=matrix(1,nrow=1522,ncol=1,byrow=FALSE)
 K=matrix(1,nrow=1522,ncol=1,byrow=FALSE)
 for(i in 1:1522){Q[i]=which.max(A[i,])}
for(i in 1:1522){if( Q[i]==1){K[i]="TRAIN"}
else if( Q[i]==2){K[i]="bus"}
 else if( Q[i]==3){K[i]="car"}
 else if( Q[i]==4){K[i]="bike"}
else  K[i]="walk"}
 n=table(Data$MainModeENG,K)
 n
table(K)

 i=table(K)
slices=c(i)
 lbls=c("Bike","BUS","Car","Rail","Walk")
 pct = round(slices/sum(slices)*100)
lbls <- paste(lbls, pct) # add percents to labels 
lbls <- paste(lbls,"%",sep="") # ad % to labels 
pie3D(slices,labels=lbls,explode=.1,main="Predicted Mode Share")

         K
          bike bus car TRAIN walk
  Bicycle   38   0  21   
  Bus        1   0  34     1    5
  Car       52   0 329    78   53
  Rail       1  15  62   442    8
  Walk       6   0   5    19  200
 round(prop.table(n,1),2)
         K
          bike  bus  car TRAIN walk
  Bicycle 0.18 0.00 0.10  0.26 0.46
  Bus     0.02 0.00 0.83  0.02 0.12
  Car     0.10 0.00 0.64  0.15 0.10
  Rail    0.00 0.03 0.12  0.84 0.02
  Walk    0.03 0.00 0.02  0.08 0.87
round(prop.table(n,2),2)
         K
          bike  bus  car TRAIN walk
  Bicycle 0.39 0.00 0.05  0.09 0.27
  Bus     0.01 0.00 0.08  0.00 0.01
  Car     0.53 0.00 0.73  0.13 0.15
  Rail    0.01 1.00 0.14  0.74 0.02
  Walk    0.06 0.00 0.01  0.03 0.55

#elaticity direct for Bus
elasticity=b[5]*(Data$TotalTimeBus/100)*(1-Pbus)
 elasticity
#cross elasticity
 elasticitycross=-b[5]*(Data$TotalTimeBus/100)*(Pbus)
 elasticitycross
#aggregate elasticitiespbus
 y=colSums(Pbus*elasticity)
 r=colSums(Pbus)
 w=y/r
 w
#aggregate cross elasticities
i=colSums(Pbus*elasticitycross)
 o=colSums(Pbus)
 p=i/o
 p
#elaticity direct for Train
elasticity=b[5]*(Data$TotalTimeTrain/100)*(1-Ptrain)
 elasticity
#cross elasticitytrain
 elasticitycross=-b[5]*(Data$TotalTimeTrain/100)*(Ptarin)
 elasticitycross
#aggregate elasticitiespbus
 y=colSums(Ptrain*elasticity)
 r=colSums(Ptrain)
 w=y/r
 w
#aggregate cross elasticities
i=colSums(Ptrain*elasticitycross)
 o=colSums(Ptrain)
 p=i/o
 p

#elaticity direct for Car
elasticity=b[5]*(Data$TimeCar/100)*(1-Pcar)
 elasticity
#cross elasticity
 elasticitycross=-b[5]*(Data$TimeCar/100)*(Pcar)
 elasticitycross
#aggregate elasticitiespbus
 y=colSums(Pcar*elasticity)
 r=colSums(Pcar)
 w=y/r
 w
#aggregate cross elasticities
i=colSums(Pcar*elasticitycross)
 o=colSums(Pcar)
 p=i/o
 p

#elaticity direct for bike
elasticity=b[5]*(Data$TimeBike/100)*(1-Pbike)
 elasticity
#cross elasticity
 elasticitycross=-b[5]*(Data$TimeBike/100)*(Pbike)
 elasticitycross
#aggregate elasticitiespbus
 y=colSums(Pbike*elasticity)
 r=colSums(Pbike)
 w=y/r
 w
#aggregate cross elasticities
i=colSums(Pbike*elasticitycross)
 o=colSums(Pbike)
 p=i/o
 p
#elaticity direct for walk
elasticity=b[5]*(Data$TimeWalk/100)*(1-Pwalk)
 elasticity
#cross elasticity
 elasticitycross=-b[5]*(Data$TimeWalk/100)*(Pwalk)
 elasticitycross
#aggregate elasticitieswalk
 y=colSums(Pwalk*elasticity)
 r=colSums(Pwalk)
 w=y/r
 w
#aggregate cross elasticities
i=colSums(Pwalk*elasticitycross)
 o=colSums(Pwalk)
 p=i/o
 p



