2021年12月・2022年1月号の「実証ビジネス・エコノミクス」の第5回「競争の激しさをデータで読み解く：参入ゲームの推定［基礎編］」に関するRプログラムの紹介である。本プログラムによって誌面で紹介した分析を再現しつつ、実際に手を動かしながら学ぶことを目的とする。
本連載の内容、およびサンプルコード等の資料は、情報の提供のみを目的としていますので、運用につきましては十分にご確認をいただき、お客様ご自身の責任とご判断によって行ってください。これらの情報の運用結果により損害等が生じた場合でも、日本評論社および著者はいかなる責任を負うことはできませんので、ご留意ください。
今回の病院レベルのデータは、橋本千代さん(Stanford University)に提供して頂きました。 また、楠井俊朗さん(東京大学工学部)には追加的なデータの収集に協力して頂き、プログラムの作成に際して塩田凌平さん(東京大学大学院経済学研究科)にご尽力頂きました。 この場を借りて、御礼申し上げます。
# ワークスペースを初期化
rm(list = ls())
# パッケージを読み込む。
require(tidyverse)## Loading required package: tidyverse## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.1.1 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'purrr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'stringr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()require(optimx)## Loading required package: optimxまず病院レベルのデータを読み込みます。
data <- readr::read_csv(file = "source/MRIData.csv", locale = locale(encoding = "UTF-8"))##
## -- Column specification --------------------------------------------------------
## cols(
## HospitalID = col_double(),
## CityCode = col_double(),
## Kyukyu = col_double(),
## Kinou = col_double(),
## Sien = col_double(),
## Hyoka = col_double(),
## DepNeurology = col_double(),
## DepNeurosurgery = col_double(),
## NumBeds = col_double(),
## MRINumOwn = col_double(),
## Tesla = col_double(),
## Population = col_double(),
## Menseki = col_double(),
## PopDensity = col_double(),
## TaxableIncome = col_double(),
## MRIOwnDum = col_double()
## )data <- data.frame(data)今回のデータは病院レベルのデータなので、推定のために市区町村レベルの集計データに変換します。 そのために、まず市区町村コードのリストを最初に作成し、それに続いて（長さが市区町村数の）各変数のベクトルを準備します。
listCode <- unique(data$CityCode)
NumMRI <- numeric(length(listCode))
Pop <- numeric(length(listCode))
Menseki <- numeric(length(listCode))
PopDen <- numeric(length(listCode))
Income <- numeric(length(listCode))準備した各変数のベクトルに、各市区町村の数を逐次代入します。
for (i in 1:length(listCode)){
subdata <- data[data$CityCode==listCode[i],]
NumMRI[i] <- sum(subdata$MRIOwnDum)
Pop[i] <- as.integer(unique(subdata$Population))
Menseki[i] <- as.integer(unique(subdata$Menseki)[1])
PopDen[i] <- unique(subdata$PopDensity)
Income[i] <- as.integer(unique(subdata$TaxableIncome))
}上で作成した市区町村レベルの変数のデータをdatasetとして格納します。
dataset <- data.frame(Code = listCode,
NumMRI = NumMRI,
Pop = Pop,
Menseki = Menseki,
PopDen = PopDen,
Income = Income)変数に欠損がある市区町村があるため、今回は簡便な処理として観測から落とし、最後に人口を百万で除します。最終的に推定で用いるサンプルサイズ（市区町村数）は1459となります。
# cleaning
dataset <- na.omit(dataset)
dataset$Pop <- dataset$Pop/1000000
M <- nrow(dataset)
print(M)## [1] 1459本文の表3からもわかる通り、MRIを導入している病院数が6を超える市区町村は全体の10%程度と、サンプルサイズはそこまで多くないため、ここでは\(N_{max}=6\)として、\(N_{max}\)以上の病院がMRIスキャナーを購入している市区町村の参入企業数を\(N_{max}\)に置換します。（また、本文では頑健性を試すために、\(N_{max}=7,\dots,10\)も試しています。興味のある読者は、各自でチェックをお願いします。）
N_max <- 6
# MRIを導入している病院数がN_maxより大きい市区町村のMRI導入済み病院数をN_maxに置換する
dataset$NumMRI[dataset$NumMRI>N_max]<- N_maxBresnahan and Reiss (1991b)に従って、目的関数を設定します。ただし、以下で出てくる\(V\)や\(VV\)は、\(N_{max} = 6\)でれば、 \(V=\begin{bmatrix} \alpha_1 & 0 & 0 & 0 & 0 & 0 \\ \alpha_1 & -\alpha_2 & 0 & 0 & 0 & 0 \\ \alpha_1 & -\alpha_2 & -\alpha_3 & 0 & 0 & 0 \\ \alpha_1 & -\alpha_2 & -\alpha_3 & -\alpha_4 & 0 & 0 \\ \alpha_1 & -\alpha_2 & -\alpha_3 & -\alpha_4 & -\alpha_5 & 0 \\ \alpha_1 & -\alpha_2 & -\alpha_3 & -\alpha_4 & -\alpha_5 & -\alpha_6 \end{bmatrix}\)、 \(VV\)は\(\begin{bmatrix} \alpha_1 \\ \alpha_1 - \alpha_2 \\ \alpha_1 - \alpha_2 - \alpha_3 \\ \alpha_1 - \alpha_2 - \alpha_3 -\alpha_4 \\ \alpha_1 - \alpha_2 - \alpha_3 - \alpha_4 - \alpha_5 \\ \alpha_1 - \alpha_2 - \alpha_3 - \alpha_4 - \alpha_5 - \alpha_6 \end{bmatrix}\)がM回（市場の個数分）繰り返される行列であり、病院がMRIスキャナーを購入した際の利潤を定義することに役立つ行列です。
obj <- function(params, dataset, N_max){
# 本文に沿ってパラメターの名前をつけ直す
# ただし、alpha2はマイナス1を乗じていることに注意
alpha1 <- params[1]
alpha2 <- -params[2:N_max]
alpha <- c(alpha1, alpha2)
gamma <- params[N_max+1]
# 推定に用いる変数をデータから抽出する
NumMRI <- dataset$NumMRI
M <- nrow(dataset)
pop <- as.matrix(dataset$Pop)
pop <- matrix(pop, M, N_max)
# 可変利潤を定義する際に便利な行列（上で説明済）を定義する
V <- matrix(0, N_max, N_max)
for (i in 1:length(alpha)){
V[i:N_max, i] <- alpha[i]
}
# 可変利潤部分の行列（ただしPop_mを乗じる前）を定義する
VV <- t(V %*% matrix(1, N_max, M))
# 固定費用部分の行列を定義する
F <- matrix(gamma, M, N_max)
# 各市場mの企業の参入時の利潤を定義する
pi <- pop*VV-F
# 標準正規分布のCDFに各値を代入する
phi <- pnorm(pi, mean=0, sd=1)
# 端点（n=0 および n=N_max）に注意しながら、各企業数になる確率を計算する
mat <- cbind(matrix(1, M, 1)-phi[,1], phi[,1:N_max-1]-phi[,2:N_max], phi[,N_max])
# データで観測される企業数の確率を抽出する
ml <- rep(0, M)
for (i in 0:N_max){
ml[NumMRI==i] <- log(mat[, i+1][NumMRI==i])
}
# 対数尤度が定義できない場合（確率が小さすぎて対数をとれない場合）の対処
ml <- replace(ml, which(is.infinite(ml)), -10000)
# 対数尤度の和をとることで目的関数の値を計算する
val <- sum(ml)/M
return(val)
}目的関数の最大化を行います。ここではoptimxを用いていて、初期値を全て1として始めています。
# パラメターの初期値を設定
initial <- rep(1, N_max+1)
# 最適化
res <- optimx(par = initial, obj, dataset = dataset, N_max = N_max, hessian=T, method='L-BFGS-B', lower=0, control = list(fnscale=-1))標準誤差を計算します。（その際に用いるHessian行列はoptimxが必要であることに注意）
Hessian <- gHgen(par = as.numeric(res[0:N_max+1]), obj, dataset = dataset, N_max = N_max)
se <- sqrt(diag(solve(-Hessian$Hn)/M)) 推定値と標準誤差を表示します。上から\(\alpha_0, \alpha_1, \cdots\)で、最終行が\(\gamma_0\)の推定値・標準誤差です。
estimates <- as.numeric(res[0:N_max+1])
cbind(estimates, se)## estimates se
## [1,] 54.626157 1.5288645
## [2,] 33.546017 1.3621141
## [3,] 8.083133 0.4822667
## [4,] 4.236733 0.3265227
## [5,] 2.512332 0.2602603
## [6,] 1.632938 0.2151856
## [7,] 1.390409 0.0462287ここでは、上の分析で得られた推定値を用いて、\(S_n\)と\(s_n=S_n/n\)を計算し、表示します（１列目が\(S_n\)、２列目が\(s_n\)です）。
# 推定値からalphaとgammaの値を定義する
alpha <- estimates[1:N_max]
gamma <- estimates[N_max+1]
# Entry Threshold の行列を定義する
EntryThreshold <- matrix(0, N_max, 2)
# Entry Threshold(S_1とs_1)を計算して、行列に代入する
deno <- alpha[1]
EntryThreshold[1,] <- c(as.integer(gamma/deno*10^6), as.integer(gamma/deno*10^6))
# Entry Threshold(S_n, s_n, n>1)を計算して、行列に代入する
for (i in 2:N_max){
deno <- deno-alpha[i]
EntryThreshold[i,] <- c(as.integer(gamma/deno*10^6), as.integer(gamma/deno*10^6/i))
}
# 行列を表示する
print(EntryThreshold)## [,1] [,2]
## [1,] 25453 25453
## [2,] 65958 32979
## [3,] 106979 35659
## [4,] 158717 39679
## [5,] 222538 44507
## [6,] 301280 50213