### Nested Logit model estimation

### Create the data array from the existing dataset
setwd("C:/Users/Guest/Desktop/kkk/")
Data <- read.csv("summer.csv",header=TRUE)
## Count the number of data rows
hh <- nrow(Data)
print(hh)

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

##### 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]
  b5 <- x[5]
  b6 <- x[6]
  b7 <- x[7]
  b8 <- x[8]

  ## Travel time to destination
  d1 <- x[9] #所要時間
  d2 <- x[10] #徒歩とチャリの所要時間
  d3 <- x[11] #費用 
  d4 <- x[12] #アクセス距離
  d5 <- x[13] #イグレス距離
#  d6 <- x[10] #駐車場料金

  d7 <- x[14] #駐車場料金 ##+d7*Data[,177]/10+d8*Data[,180]
  d8 <- x[15] #駐車場料金 ##+d7*Data[,177]+d8*Data[,180]+d9*Data[,181]
  d9 <- x[16] #駐車場料金 ##+d9*Data[,182]/1000
   
  ##Scale parameter
  pa <- x[17] 

  ## declare the log-likelihood variable, set value to 0
  LL = 0

  ### For this choice problem, we consider the following 5 modes
  ## 鉄道(train)  - (public)
  ## バス(bus)    - (public)
  ## 自動車(car)  - (private)
  ## 自転車(bike) - (private)
  ## 徒歩(walk)   - (private)
  
  ## calculate the utility function: :introduce the desired explanatory variables in the function
                                        # time           # fare      # constant
  train  <- Data[,44] *exp(d1*Data[,55]/10 +d3*Data[,47]/1000   +d4*Data[,51]/1000+d5*Data[,53]/1000	+b1*matrix(1,nrow =hh,ncol=1))
  bus    <- Data[,62] *exp(d1*Data[,70]/10 +d3*Data[,65]/1000   +d4*Data[,66]/1000+d5*Data[,68]/1000	+b2*matrix(1,nrow =hh,ncol=1))
  car    <- Data[,37] *exp(d1*Data[,43]/10 +d3*Data[,39]/1000  						+b3*matrix(1,nrow =hh,ncol=1))
  bike   <- Data[,56] *exp(d2*Data[,58]/10 								+b4*matrix(1,nrow =hh,ncol=1))
  walk   <- Data[,59] *exp(d2*Data[,61]/10                                                     )

  train2  <- Data[,190] *exp(d1*Data[,200]/10 +d3*Data[,192]/1000   +d4*Data[,196]/1000+d5*Data[,198]/1000	+b5*matrix(1,nrow =hh,ncol=1))
  bus2    <- Data[,207] *exp(d1*Data[,213]/10 +d3*Data[,210]/1000   +d4*Data[,211]/1000+d5*Data[,212]/1000	+b6*matrix(1,nrow =hh,ncol=1))
  car2    <- Data[,186] *exp(d1*Data[,189]/10 +d3*Data[,188]/1000 						+b7*matrix(1,nrow =hh,ncol=1))
  bike2   <- Data[,201] *exp(d2*Data[,203]/10 									+b8*matrix(1,nrow =hh,ncol=1))
  walk2   <- Data[,204] *exp(d2*Data[,206]/10                                                     )

  ### Construct nests.First, calculate public or private choice probabilityP(LV1)
  ## setting of log-sum variables
  ##logsum.public→西区中区
  ##logsum.private→自分の市区町村

  logsum.public  <- log( ( (train+bus+car+bike+walk)!=0) * (train + bus + car + bike + walk) + ((train+bus+car+bike+walk)==0))
  logsum.private <- log( ( (train2+bus2+car2+bike2+walk2)!=0) * (train2 + bus2 + car2 + bike2 + walk2) + ((train2+bus2+car2+bike2+walk2)==0))

  ##Calculate public or private choice probability with these log-sum variables
  nume.public   <- (logsum.public != 0)*exp(pa*logsum.public+d7*Data[,177]+d8*Data[,180]+d9*Data[,182]/1000) + (logsum.public == 0)
  nume.private  <- (logsum.private != 0)*exp(pa*logsum.private) + (logsum.private == 0)

  deno <- nume.public + nume.private
  P.public  <- nume.public  / deno
  P.private <- 1-P.public #nume.private / deno

  ### Next, calculate conditional probability.P(LV2|LV1)
  ## conditional probability under public choice
  deno.public     <- train+bus+car+bike+walk
  P.train.public  <- train / ((deno.public!=0)*deno.public + (deno.public==0))
  P.bus.public    <- bus   / ((deno.public!=0)*deno.public + (deno.public==0))
  P.car.public   <- car   / ((deno.public!=0)*deno.public + (deno.public==0))
  P.bike.public  <- bike  / ((deno.public!=0)*deno.public + (deno.public==0))
  P.walk.public  <- walk  / ((deno.public!=0)*deno.public + (deno.public==0))

  ## conditional probability under private choice
  deno.private    <- train2+bus2+car2+bike2+walk2
  P.train.private  <- train2 / ((deno.private!=0)*deno.private + (deno.private==0))
  P.bus.private    <- bus2   / ((deno.private!=0)*deno.private + (deno.private==0))
  P.car.private   <- car2   / ((deno.private!=0)*deno.private + (deno.private==0))
  P.bike.private  <- bike2  / ((deno.private!=0)*deno.private + (deno.private==0))
  P.walk.private  <- walk2  / ((deno.private!=0)*deno.private + (deno.private==0))

  ### Finally, calculate joint probability. P(LV1, LV2) = P(LV2|LV1)*P(LV1)
  P.train <- P.train.public * P.public
  PP.train <- P.train.private * P.private 
  P.bus   <- P.bus.public   * P.public
  PP.bus <- P.bus.private * P.private 
  P.car   <- P.car.public   * P.public
  PP.car <- P.car.private  * P.private
  P.bike  <- P.bike.public   * P.public
  PP.bike <- P.bike.private * P.private
  P.walk  <- P.walk.public   * P.public
  PP.walk <- P.walk.private * P.private

  ## Avoid problems stemming from choice probabilities becoming zero.
  P.train <- (P.train!=0)*P.train + (P.train==0)
  P.bus   <- (P.bus!=0)  *P.bus   + (P.bus==0)
  P.car   <- (P.car!=0)  *P.car   + (P.car==0)
  P.bike  <- (P.bike!=0) *P.bike  + (P.bike==0)
  P.walk  <- (P.walk!=0) *P.walk  + (P.walk==0)
  PP.train <- (PP.train!=0)*PP.train + (PP.train==0)
  PP.bus   <- (PP.bus!=0)  *PP.bus   + (PP.bus==0)
  PP.car   <- (PP.car!=0)  *PP.car   + (PP.car==0)
  PP.bike  <- (PP.bike!=0) *PP.bike  + (PP.bike==0)
  PP.walk  <- (PP.walk!=0) *PP.walk  + (PP.walk==0)

  ## Choice results
  C.train   <- Data[,185] =="2001"
  C.bus     <- Data[,185] =="2401"
  C.car     <- Data[,185] =="1001"
  C.bike    <- Data[,185] =="4101"
  C.walk    <- Data[,185] =="4201"
  CC.train   <- Data[,185] =="2000"
  CC.bus     <- Data[,185] =="2400"
  CC.car     <- Data[,185] =="1000"
  CC.bike    <- Data[,185] =="4100"
  CC.walk    <- Data[,185] =="4200"

  ## Calculate the Log-likelihood function
  LL <- colSums(C.train*log(P.train) + C.bus*log(P.bus) + C.car  *log(P.car)   + C.bike  *log(P.bike) +C.walk *log(P.walk)
		+CC.train*log(PP.train) + CC.bus*log(PP.bus) + CC.car  *log(PP.car)   + CC.bike  *log(PP.bike) +CC.walk *log(PP.walk))

}

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

##Parameter optimization 
res <- optim(b0,fr, method = "BFGS", 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)
