﻿### Multinomial Logit model estimation

###

# データファイルの読み込み

Data <- read.csv("model2015_5.csv",header=T)



## データ数:Dataの行数を数える

hh <- nrow(Data)



## パラメータの初期値の設定

b0 <- numeric(10) 



##### Logit model の対数尤度関数の定義#####


fr <- function(x) {

### パラメータの宣言：
  

## 定数項
  
b1 <- x[1]
  
b2 <- x[2]
  
b3 <- x[3]
  
b4 <- x[4]

  

## 目的地までの所要時間
  
d1 <- x[5]
  
  

## 料金
  
f1 <- x[6]
  


#a1 <- x[7]
#i1 <- x[8]
w1 <- x[7]  
y1 <- x[8]
s1 <- x[9]
h1 <- x[10]

## 対数尤度のための変数を宣言
  
LL = 0

  

### 今回用いる目的地は以下の5つ．
  
## 鉄道(train)
  
## バス(bus)
  
## 自動車(car)
  
## 自転車(bike)
  
## 徒歩(walk)
  
  

## 効用の計算:説明変数にしたい列を入れる．
                                         # 時間                     # 料金                    # 定数項
  
train  <- Data$代替手段生成可否train*exp(d1*Data$総所要時間train/100 +f1*Data$費用train/100 + b1*matrix(1,nrow =hh,ncol=1))
  
bus    <- Data$代替手段生成可否bus  *exp(d1*Data$総所要時間bus/100   +f1*Data$費用bus/100   + b2*matrix(1,nrow =hh,ncol=1))
  
car    <- Data$代替手段生成可否car  *exp(d1*Data$所要時間car/100                        + b3*matrix(1,nrow =hh,ncol=1))
  
bike   <- Data$代替手段生成可否bike *exp(h1*Data$hight3 + s1*Data$female + w1*Data$work + y1*Data$young + d1*Data$所要時間bike/100 + b4*matrix(1,nrow =hh,ncol=1))
  
walk   <- Data$代替手段生成可否walk *exp(d1*Data$所要時間walk/100                                                     )

  

### 選択確率の計算
  
## 分母となる，各々のexp(V)の和をつくる
  
deno <- (car + train + bus + bike + walk)

  

## それぞれ計算する
  
Ptrain <- Data$代替手段生成可否train*(train / deno)
  
Pbus   <- Data$代替手段生成可否bus  *(bus   / deno)
  
Pcar   <- Data$代替手段生成可否car  *(car   / deno)
  
Pbike  <- Data$代替手段生成可否bike *(bike  / deno)
  
Pwalk  <- Data$代替手段生成可否walk *(walk  / deno)

  

## 選択確率が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)
    


  

## 選択結果
  
Ctrain   <- Data$代表交通手段 =="鉄道"
  
Cbus     <- Data$代表交通手段 =="バス"
  
Ccar     <- Data$代表交通手段 =="自動車"
  
Cbike    <- Data$代表交通手段 =="自転車"
  
Cwalk    <- Data$代表交通手段 =="徒歩"

  

## 対数尤度の計算
  
LL <- colSums(Ctrain*log(Ptrain) + Cbus*log(Pbus) +

                Ccar  *log(Pcar)   + Cbike  *log(Pbike) +Cwalk *log(Pwalk))


}



##### 対数尤度関数frの最大化####
#

##パラメータ値の最適化 

res <- optim(b0,fr, method = "Nelder-Mead", hessian = TRUE, control=list(fnscale=-1,maxit=1000000000))

## パラメータ推定値、ヘッセ行列 

b   <- res$par

hhh <- res$hessian



## t値の計算

tval <- b/sqrt(-diag(solve(hhh)))



## 初期尤度

L0 <- fr(b0)

## 最終尤度

LL <- res$value


##### 結果の出力 #####

print(res)

## 初期尤度

print(L0)

## 最終尤度

print(LL)

##ρ^2値

print((L0-LL)/L0)

## 修正済ρ^2値

print((L0-(LL-length(b)))/L0)

##パラメータ推定値 

print(b)

## t値 

print(tval)




# Display all estimation results
results <- cbind(b,tval)
#rownames(results) <- c("Param_cost","Param_time","Cbus","Ctrain")
Goodness_of_fit <- cbind(L0,LL,rho,rho.adj)
rownames(Goodness_of_fit) <- c("GoF")


res <- function (a,b) {
  cat("-------------Estimation results--------------", "\n")
  print(a)
  cat("--------------Goodness of fit----------------", "\n")
  print(b)
  cat("---------------------------------------------", "\n")
}

res(results,Goodness_of_fit)

