### Nested Logit model estimation on Peak-segregated data 
 ### Create the data array from the existing dataset
Data <- read.csv("F:/Tokyo Trip 2016/Data/R/Data_PP5.csv",header=T, sep=",")
 ## 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(12) 
   ##### 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]
             
   ## Fare
  ##f1 <- x[6]
 ##Scale parameter
 pa <- x[6] 
 a1 <- x[7]
 a2 <- x[8]
 p1<-x[9]
 p2<-x[10]
 p3<-x[11]
 p4<-x[12]
 #a3 <- x[13]
#a4 <- x[10]
   ## declare the log-likelihood variable, set value to 0
   LL = 0
             
## calculate the utility function: :introduce the desired explanatory variables in the function
   train  <- Data$ModeAvailableTrain*exp(d1*Data$TotalTimeTrain/100+p1*Data$Peak  +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+p2*Data$Peak          +b3*matrix(1,nrow =hh,ncol=1))
   bike   <- Data$ModeAvailableBike *exp(d1*Data$TimeBike/100+p3*Data$Peak        +b4*matrix(1,nrow =hh,ncol=1))
   walk   <- Data$ModeAvailableWalk *exp(d1*Data$TimeWalk/100+p4*Data$Peak )
   ### Construct nests.First, calculate public or private choice probabilityP(LV1)
   ## setting of log-sum variables
  logsum.act  <- log( ( (bike +walk)!=0) * (bike + walk) + ((bike+walk)==0))
  logsum.nonact <- log( ( (car + bus+ train)!=0) * (car + bus+ train) + ((car + bus+ train)==0))
            
  ##Calculate public or private choice probability with these log-sum variables
  nume.act   <- (logsum.act != 0)*exp(pa*logsum.act+a2*Data$Leisure) + (logsum.act == 0) 
  nume.nonact  <- (logsum.nonact != 0)*exp(pa*logsum.nonact + a1*Data$AgeFemale45) + (logsum.nonact == 0) 
                  
  deno <- nume.act + nume.nonact
  P.act <- nume.act / deno
  P.nonact <- nume.nonact / deno
                      
  ### Next, calculate conditional probability.P(LV2|LV1)
  ## conditional probability under public choice
 deno.act    <- bike + walk
 P.bike.act  <- bike / ((deno.act!=0)*deno.act + (deno.act==0))
 P.walk.act   <- walk   / ((deno.act!=0)*deno.act+ (deno.act==0))
                    
## conditional probability under private choice
 deno.nonact   <- car + bus + train
 P.car.nonact  <- car   / ((deno.nonact!=0)*deno.nonact + (deno.nonact==0))
 P.bus.nonact  <- bus  / ((deno.nonact!=0)*deno.nonact + (deno.nonact==0))
 P.train.nonact <- train / ((deno.nonact!=0)*deno.nonact + (deno.nonact==0))
                    
 ### Finally, calculate joint probability. P(LV1, LV2) = P(LV2|LV1)*P(LV1)
   P.bike <- P.bike.act * P.act
   P.walk   <- P.walk.act * P.act
   P.car   <- P.car.nonact  * P.nonact
   P.bus <- P.bus.nonact * P.nonact
   P.train  <- P.train.nonact * P.nonact
   ## 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)
   ## 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(P.train) + Cbus*log(P.bus) + Ccar  *log(P.car)   + Cbike  *log(P.bike) +Cwalk *log(P.walk))
                                     
 }
     ##### Maximize the Log-likelihood function#####
     
         ##Parameter optimization 
         res <- optim(b0,fr, upper=c(Inf, Inf, Inf, Inf, Inf, Inf, Inf), method = "L-BFGS-B", 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)
                
                  