1 はじめに

経済セミナー2021年10月・11月号の「実証ビジネス・エコノミクス」の第3回「合併の効果は需要次第：消費者需要モデルの推定【応用編】」に関するRプログラムの紹介である。

1.1 留意事項

本連載の内容、およびサンプルコード等の資料は、情報の提供のみを目的としていますので、運用につきましては十分にご確認をいただき、お客様ご自身の責任とご判断によって行ってください。これらの情報の運用結果により損害等が生じた場合でも、日本評論社および著者はいかなる責任を負うことはできませんので、ご留意ください。

1.2 謝辞

今回のプログラムの作成に際して、島本幸典さん・牧野圭吾さん・山田雅広さん(東京大学大学院経済学研究科)にご尽力頂きました。

2 Rに関する下準備

今回はkableExtra,sjmisc,tictocを新たに利用する。なお、以下のrequireで読み込んでいるパッケージについてまだインストールしていない場合には、最初にインストールすること。

# ワークスペースを初期化
rm(list = ls())

# (もしまだなら) パッケージのインストール
# install.packages("tidyverse")
# install.packages("kableExtra")
# install.packages("sjmisc")
# install.packages("knitr")
# install.packages("tictoc")

# パッケージを読み込む
require(tidyverse)

## Warning: package 'tidyverse' was built under R version 3.6.3

## Warning: package 'tibble' was built under R version 3.6.3

## Warning: package 'tidyr' was built under R version 3.6.3

## Warning: package 'readr' was built under R version 3.6.3

## Warning: package 'purrr' was built under R version 3.6.3

## Warning: package 'dplyr' was built under R version 3.6.3

## Warning: package 'stringr' was built under R version 3.6.3

## Warning: package 'forcats' was built under R version 3.6.3

require(kableExtra)

## Warning: package 'kableExtra' was built under R version 3.6.3

require(sjmisc)
require(knitr)

## Warning: package 'knitr' was built under R version 3.6.3

require(tictoc)

## Warning: package 'tictoc' was built under R version 3.6.3

また、tictocパッケージを用いて時間計測を行う。tictoc::tic()で時間計測を開始し、tictoc::toc()で計測を終了し、そこまでの時間がレポートされる。推定に比較的時間がかかるような場合においては、コードの時間計測を適宜行うと作業の見通しが立ちやすくなる。

tictoc::tic()

3 データの準備

連載第三回の「6 入れ子型ロジットモデルの推定とその応用」までで作成したデータに、二つの変数（‘FuelType’, ‘capacity’）を加えたデータを読み込む

data <- readr::read_csv(file = "chap3_data.csv")

つづいて、各種ダミー変数の作成を行う。なおここで使っているto_dummyはsjmiscパッケージに含まれる関数

# ダミー変数の作成
# 外国車ダミー
data <- data %>%
  mutate(Foreign_d = if_else(Type=="Foreign", 1, 0))
# FuelType=レギュラーのダミー
data <- data %>%
  mutate(FuelRegular_d = if_else(FuelType=="レギュラー", 1, 0))
# 年ダミー
data <- data %>%
  mutate(year = as.factor(year)) %>%
  to_dummy(year, suffix = "label") %>%
  bind_cols(data) %>%
  mutate(year = as.numeric(year))
# year 2006を基準に
data <- data %>%
  select(-year_2006)
# ダミー変数を最後の列に移動
data <- data %>%
  relocate(starts_with("Type_"), starts_with("year_"), starts_with("Maker_"), .after = last_col())
# capacity ダミー(4以下と5以上で分ける)
data <- data %>%
  mutate(capacity_d = if_else(capacity > 4, 1, 0))

4 ランダム係数ロジットモデルの推定

前回の連載に従って、ランダム係数ロジットモデルを推定する。なお、コード自体は前回とほぼ同様であるため、詳細な説明は省く。詳しくは前回の連載のコードを参照されたい。

4.1 定式化

以下の定式化を考える。なお、前回から定式化が変化していることに留意されたい。

\[ \begin{eqnarray} u_{i j t}=& \beta^{0}+\beta^{size}{size}_{j t}+\beta^{\textit{fuelefficiency}}\textit{fuelefficiency}_{jt} +\beta^{h p p w} h p p w_{j t}+\beta^{capacity\_d}{capacity\_d_{j t}}\\ &+\beta^{FuelRegular\_d}{FuelRegular\_d_{jt}}+\beta^{Foreign\_d}{Foreign\_d_{jt}} -\alpha_{i} p_{j t}+\xi_{j t}+\xi_t+\epsilon_{i j t} \end{eqnarray} \]

ここで、\(\alpha_i\)は正規分布に従うランダム係数である。

4.2 データの下準備

まずは、データセットから、行列を抽出する。

# まずデータをソートする。マーケット順、メーカー順かつ価格順
data %>% 
  arrange(year, Maker, price) -> data
# データセットより抽出した行列等をdatalistに格納する
datalist = list()
# 平均効用に入ってくる部分。内生性がある価格を含んでいる。
datalist$X1 <-
  data %>%
  mutate(cons = 1) %>%
  # 価格を必ず2つ目にする
  select(cons, price, FuelEfficiency, hppw, size, 
         capacity_d, FuelRegular_d, Foreign_d,
         starts_with("year_")
  ) %>% 
  as.matrix()
# ランダム係数とInteractする部分。
datalist$X2 <-
  data %>%
  mutate(cons = 1) %>%
  # 価格を必ず1つ目にする
  select(price) %>% 
  as.matrix()
# 操作変数行列。ここは外生変数と追加的な操作変数を含む
datalist$Z <-
  data %>%
  mutate(cons = 1) %>%
  select(
    cons, FuelEfficiency, hppw, size,
    capacity_d, FuelRegular_d, Foreign_d,
    starts_with("year_"),
    starts_with("iv_GH"),
    -ends_with("nest")
  ) %>%
  as.matrix()
# 市場シェア
datalist$ShareVec <-
  data %>%
  select(share) %>%
  as.matrix()
# 
datalist$marketindex <- data$year
datalist$logitshare <- data$logit_share
# 観察数、年数を取得
N <- length(datalist$marketindex)
T <- length(unique(datalist$marketindex))

次に、乱数を用意しておく。

set.seed(111)
Nsim <- 1000
datalist$draw_vec <- rnorm(Nsim*ncol(datalist$X2))

parameter = list()
parameter$Nsim <- Nsim
parameter$T = T 
parameter$N = N
marketindex <- datalist$marketindex
uniquemarketindex = sort(unique(marketindex))
temp1 = matrix( rep(uniquemarketindex, N), T, N) %>% t()
temp2 = matrix( rep(marketindex, T), N, T)
mkt_denom_d = (temp1 == temp2)*1
datalist$mkt_denom_d <- mkt_denom_d

4.3 各種関数の定義

前回と同様にGMM推定に必要な関数を定義していく。これらは前回と同じである。

市場シェアを計算する関数 f_marketshare
縮小写像によるBerryインバージョンを行う関数 f_contraction
GMMの目的関数 f_GMMobj
標準誤差の計算 f_SE
価格弾力性の計算 f_elasticity

4.3.1 市場シェアの計算コード

ここでは、市場シェアを計算する関数f_mktshareを定義しよう。

f_mktshare <- function(datalist, parameter, delta, theta2){
  # datalist, parameterの中身を読み込む
  X2 <- datalist$X2
  N <- length(datalist$marketindex)
  T <- length(unique(datalist$marketindex))
  draw_vec <- datalist$draw_vec
  marketindex <- datalist$marketindex
  mkt_denom_d <- datalist$mkt_denom_d
  Nsim <- parameter$Nsim
  K2 <- length(theta2)
  # Nonlinear component. mu (N, ns) matrix.
  mu <- X2 %*% diag(x = theta2, nrow = K2) %*% matrix(draw_vec, nrow = K2)
  delta_mu <- delta %*% matrix(1, nrow = 1, ncol = Nsim) + mu
  exp_delta_mu <- exp(delta_mu)
  # ロジット確率の分母を計算する。これは、市場ごとに和を計算する。
  denom_temp <- t(t(exp_delta_mu) %*% mkt_denom_d)
  denom_outside <- exp(matrix(0, T, Nsim))
  denom_temp <- denom_temp + denom_outside
  # denom is (N * ns) matrix that captures the denominator of the logit prob.
  denom <- mkt_denom_d %*% denom_temp
  # Choice prob for individual (N * ns) matrix
  s_jt_i <- exp_delta_mu / denom
  # Market share (N * 1) vector
  s_jt <- apply(s_jt_i, 1, mean)
  return(s_jt)
}

4.3.2 縮小写像によるBerryインバージョンのコード

ここでは、Berryインバージョンによる平均効用の計算を行う関数を定義する。なお、関数の内部でStep2で定義したf_mktshareを利用している。

f_contraction <- function(datalist, parameter, delta_ini, theta2){
  # 初期値の代入
  ShareVec <- datalist$ShareVec
  exp_delta_old <- exp(delta_ini)
  # 収束の設定
  tol <- 1e-11
  norm <- 1e+10
  iter <- 0
  while (norm > tol && iter < 1000){
    # マーケットシェアの計算
    pred_mkt_share <- f_mktshare(datalist, parameter, log(exp_delta_old), theta2)
    # deltaを更新
    exp_delta <- exp_delta_old * ShareVec / pred_mkt_share
    # 距離を評価
    norm <- max(abs(exp_delta - exp_delta_old))
    # 計算結果を次のループに渡す
    exp_delta_old <- exp_delta 
    iter <- iter + 1
  }
  return(log(exp_delta))
}

4.3.3 GMM 目的関数を定義

f_GMMobj <- function(theta2, parameter, datalist, option) {
  # 引数option
  # 0:最適化の際に指定
  # 1:推定値やmean utility(delta)を得る際に指定
  
  
  
  # Contractionの初期値として、グローバル変数のdelta_globalを用いる。
  delta_ini <- delta_global
  
  # Contraction mapping
  delta <- f_contraction(datalist, parameter, delta_ini, theta2)
  
  # 平均効用を永続代入演算子でアップデートしている。 
  delta_global <<- delta
  
  # 目的関数に必要な行列を準備
  X1 <- datalist$X1
  Z <- datalist$Z
  # Weight matrix W
  if (datalist$weight_mat_option == "2SLS") {
    W <- solve(t(Z) %*% Z)
  } else if (datalist$weight_mat_option == "Ident") {
    W <- eye(dim(Z)[2])
  } else if (is.null(datalist$weight_mat_option)) {
    stop("SET 'datalist$weight_mat_option' !!")
  }
  # Obtain linear parameter by regression
  theta1 <- solve(t(X1) %*% Z %*% W %*% t(Z) %*% X1) %*% t(X1) %*% Z %*% W %*% t(Z) %*% delta
  # xiを計算
  Xi <- delta - X1 %*% theta1
  
  # Objective function
  output <- t(Xi) %*% Z %*% W %*% t(Z) %*% Xi
  # optionによって戻り値を変える
  if (option == 0) {
    return(output = output)
  } else if (option == 1) {
    return(list(output = output, theta1 = theta1, delta = delta, Xi = Xi))
  } else {
    stop(
      "SET the argument 'option' to 0 or 1 !! \n 
      0: Specify when optimizing \n
      1: Specify when getting coefficient and mean utility (delta)"
    )
  }
}

4.3.4 標準誤差の計算

f_SE <- function(theta2, parameter, datalist) {
  param <- parameter
  
  delta_ini <- datalist$logitshare
  
  # Contraction mapping
  delta <- f_contraction(datalist, param, delta_ini, theta2)
  
  # Obj function
  X1 <- datalist$X1
  
  # IV matrix
  Z <- datalist$Z
  
  # Weight matrix W
  if (datalist$weight_mat_option == "2SLS") {
    W <- solve(t(Z) %*% Z)
  } else if (datalist$weight_mat_option == "Ident") {
    W <- eye(dim(Z)[2])
  } else if (is.null(datalist$weight_mat_option)) {
    STOP("SET 'datalist$weight_mat_option' !!")
  }
  
  # Obtain linear parameter by regression
  theta1 <- solve(t(X1) %*% Z %*% W %*% t(Z) %*% X1) %*% t(X1) %*% Z %*% W %*% t(Z) %*% delta
  
  # ResidualをGetする。
  Xi <- delta - X1 %*% theta1
  
  # Omega
  Omega_hat <- matrix(0, nrow = ncol(Z), ncol = ncol(Z))
  for (ii in 1:param$N){
    Omega_hat <- Omega_hat + ( matrix(Z[ii, ]) %*% t(matrix(Z[ii, ])) * Xi[ii]^2 )/parameter$N
  }
  
  # Gradient of delta. 単純化のために数値微分で計算する。
  Ddelta <- matrix(0, nrow = parameter$N, ncol = length(theta2))
  
  for (k in 1:length(theta2)) {
    tempparam <- param
    temptheta2 <- theta2
    temptheta2[k] <- theta2[k] + 1e-6
    
    delta_add <- f_contraction(datalist, tempparam, delta_ini, temptheta2)
    
    Ddelta[, k] <- (delta_add - delta) / 1e-6
  }
  
  G <- (parameter$N^(-1)) * t(Z) %*% cbind(-X1, Ddelta)
  #G <-  t(Z) %*% cbind(-X1, Ddelta)
  
  # Asymptotic Var-Cov Matrix
  AsyVarMat <- solve(t(G) %*% W %*% G) %*% t(G) %*% W %*% Omega_hat %*% W %*% G %*% solve(t(G) %*% W %*% G)
  # Asymptotic Standard Errors.
  Ase <- sqrt( diag(AsyVarMat) / parameter$N )
  
  return(Ase)
}

4.3.5 価格弾力性を計算する関数

このコードも基本的に前回と同様である。

f_elasticity <- function(datalist, parameter, theta1, theta2, delta) {
  # Unpack data
  X2 <- datalist$X2
  Nsim <- parameter$Nsim
  N <- parameter$N
  price <- datalist$X1[, colnames(datalist$X1) == "price"] %>% as.matrix()
  K2 <- length(theta2)
  draw_vec <- datalist$draw_vec
  uniquemarketindex <- datalist$uniquemarketindex
  
  # 以下マーケットシェアの関数とほぼ同じ
  # Nonlinear component. mu (N, ns) matrix.
  mu <- X2 %*% diag(x = theta2, nrow = K2) %*% matrix(draw_vec, nrow = K2)
  
  denom_outside <- exp(matrix(0, T, Nsim))
  delta_mu <- delta %*% matrix(1, nrow = 1, ncol = Nsim) + mu
  exp_delta_mu <- exp(delta_mu)
  mkt_denom_d <- datalist$mkt_denom_d
  denom_temp <- t(t(exp_delta_mu) %*% mkt_denom_d)
  denom_temp <- denom_temp + denom_outside
  denom <- mkt_denom_d %*% denom_temp
  # Choice prob for individual (N * ns) matrix
  s_jt_i <- exp_delta_mu / denom
  
  # Price parameter alpha_i
  draw_for_price <- matrix(draw_vec, nrow = K2)[1, ]
  # theta1のなかで2つ目、theta2のなかで1つ目に価格係数がある
  alpha_i <- theta1[2] + theta2[1]*draw_for_price
  
  # 弾力性行列を保存する空のリストを作成する
  elaslist <- as.list(rep(NA, T))
  
  # 各マーケットごとに弾力性行列を作成する．
  for (t in 1:T) {
    year_beg <- 2016 - T 
    # Number of products 
    J_t <- sum(marketindex == year_beg + t)
    
    # それぞれのマーケットに注目した(J*ns)の選択確率の行列を作成する
    # あるマーケットに対して，製品ごとに集計した選択確率の行列(J*ns)を作成する
    ag_model_s_i <- s_jt_i[ marketindex == year_beg + t,]
    
    # 製品ごとの選択確率のベクトル(J * 1)
    ag_model_s <- apply(ag_model_s_i, 1, mean) %>% as.matrix()
    
    # 空の弾力性行列(J * J)を定義
    elasmat <- matrix(0, nrow = J_t, ncol = J_t)
    
    # あるマーケットに注目した，製品ごとの価格のベクトル(J*1)を作成する
    price_t <- price[marketindex == year_beg + t]
    
    # 弾力性を計算し，弾力性行列を作成する
    for (k in 1:J_t) {
      for (j in 1:J_t) {
        if (k != j) {
          # cross elasticity
          elasmat[k, j] <- (-1)*price_t[k] * ag_model_s[j, ]^(-1) * mean(alpha_i * ag_model_s_i[j, ] * ag_model_s_i[k, ])
        } else if (k == j) {
          # Own elasticity
          elasmat[k, j] <- price_t[j] * ag_model_s[j, ]^(-1) * mean(alpha_i * ag_model_s_i[j, ] * (1 - ag_model_s_i[j, ]) )
        } 
      }
    }
    # 弾力性行列をマーケットごとに名前を付けて保存する．
    elaslist[[t]] <- elasmat
    names(elaslist)[[t]] <- sprintf("elasmat_%d", t+year_beg)
  }
  return(elaslist)
}

4.4 パラメタ推定及び標準誤差の計算

上で定義した関数を用いてGMM推定を行う。

# GMMの荷重行列のオプション
datalist$weight_mat_option <- "2SLS"

# グローバル変数としてアップデートする平均効用
delta_global <- data$logit_share
# 最適化
GMM_nonlinear <-
  optim(
    par = c(0.7), 
    fn = f_GMMobj,
    method = "L-BFGS-B",
    lower = c(0), 
    parameter = parameter,
    datalist = datalist,
    option = 0
  )

# nonlinear parameter estimates
theta2_hat <- GMM_nonlinear$par

推定値に基づいて標準誤差の計算を行う。

GMM_linear <- f_GMMobj(theta2_hat, parameter, datalist, option = 1)

# deltaと観測されないXiをdataに保存
delta <- GMM_linear$delta
Xi <- GMM_linear$Xi

data <- data %>%
  mutate(Xi = Xi, delta = delta)

# linear parameter estimates
theta1_hatmat <- GMM_linear$theta1
theta1_hat <- as.vector(theta1_hatmat)

se <- f_SE(theta2 = theta2_hat, parameter = parameter, datalist = datalist)

以上に基づいて、推定値のテーブルを出力する。

theta2_hatmat <- as.matrix(theta2_hat)
rownames(theta2_hatmat)[1] <- "random_price"
beta_hat <- rbind(theta1_hatmat, theta2_hatmat) 
colnames(beta_hat) <- "coeff"
names(se)[length(GMM_linear$theta1)+1] <- "random_price"

est_table <- as.data.frame(cbind(beta_hat, se))
kable(est_table,  align=rep('c', 3), booktabs = T, 
      caption = "Estimates of The Parameters")

Estimates of The Parameters
	coeff	se
cons	-13.7359205	0.6295519
price	-2.2459258	0.7410923
FuelEfficiency	0.1936160	0.0125437
hppw	13.3473417	3.1461977
size	0.5294779	0.0682720
capacity_d	-0.2964342	0.1353892
FuelRegular_d	-1.0261309	0.2486721
Foreign_d	0.9835681	0.1734038
year_2007	-0.0991314	0.1507647
year_2008	-0.2489589	0.1477847
year_2009	-0.4799054	0.1471831
year_2010	-0.6238478	0.1645326
year_2011	-0.8349815	0.1675212
year_2012	-0.7000877	0.1691600
year_2013	-0.8510638	0.1714373
year_2014	-1.0413253	0.1834505
year_2015	-1.1387600	0.1809877
year_2016	-1.2454347	0.1833420
random_price	0.7001090	0.2376130

5 限界費用の推定

需要推定の結果に基づいて、限界費用を推定する。

差別化財ベルトラン競争モデルのナッシュ均衡条件に基づき、限界費用は以下の式で求める。 \[\hat{mc} = p^* - \Delta^{pre}(p^*)^{-1}q(p^*)\] ちなみに

\[\Delta^{pre}_{jr}(p) = \left\{ \begin{array}{ll} - \frac{\partial q_j(p)}{\partial p_r} & {\rm if} \ (r,j) \in \mathcal{F}_f \\ 0 & その他 \end{array} \right.\]

これは需要推定値から求めることができる。詳しくは本誌を参照されたい。

以上の考え方に基づいて、限界費用を推定する関数を作成する。

f_Merginal_Cost_Est <- function(data, marketindex, elasticity) {
  # マーケットの数を取り出す
  T <- length(unique(data$year))
  
  # 限界費用を保存する空のリストを作成する
  MC_list <- as.list(rep(NA, T))
  
  year_beg <- min(data$year) - 1 
  
  # 各マーケットごとに限界費用計算する
  for (t in 1:T) {
    # あるマーケットの弾力性行列を取り出す
    elasmat_t <- elasticity[[t]] 
    
    # あるマーケットの大きさ(車種の数)
    J_t <- sum(marketindex == year_beg + t)
    
    # あるマーケットのデータを抽出
    data_t <- data[ data$year == year_beg + t, ]
    
    # あるマーケットの価格とマーケットシェアを取り出す
    Pricevec_t <- data_t$price
    Sharevec_t <- data_t$share
    
    # 空の所有構造行列と弾力性の微分部分を用意する
    Ownership_t <- matrix(0, J_t, J_t)
    # 各弾力性に対して計算をする
    for (j in 1:J_t) {
      for (k in 1:J_t) {
        
        # 車種jと車種kが同じ企業で生産されていれば1を入れる
        if (data_t$Maker[j] == data_t$Maker[k]){
          Ownership_t[j,k] = 1
        }
      }
    }
    Derivative_t = - (elasmat_t)*(kronecker(matrix(1,J_t,1), t(Sharevec_t)))/(kronecker(matrix(1,1,J_t),Pricevec_t))
    # 所有構造行列に弾力性の微分部分を掛けることでDeltaを作成する
    Delta_t <- Derivative_t * Ownership_t
    
    # 限界費用を計算する
    Marginal_Cost_t <- Pricevec_t - (solve(Delta_t) %*% Sharevec_t)
    
    # データフレーム化してNameと企業名を追加する(NameIDはleft_joinに必要)
    Marginal_Cost_t_D <- data.frame(Name = data_t$Name, NameID = data_t$NameID, 
                                    Maker = data_t$Maker, 
                                    price = Pricevec_t, mc = Marginal_Cost_t) %>%
      mutate(margin = (price-mc)/price)
    # 限界費用をマーケットごとに名前を付けて保存する．
    MC_list[[t]] <- Marginal_Cost_t_D
    names(MC_list)[[t]] <- sprintf("Marginal_Cost_%d", t+year_beg)
  }
  return(MC_list)
}

作成した関数を使って限界費用を計算する。加えて価格費用マージンを計算してその分布を確認する。

# 価格弾力性行列の計算。限界費用推定において用いる。
elasticity <- f_elasticity(datalist = datalist, parameter = parameter, theta1 = theta1_hat, theta2 = theta2_hat, delta = delta)

# 作成した関数を実行する
marketindex <- datalist$marketindex
Marginal_Cost_Est <- f_Merginal_Cost_Est(data, marketindex, elasticity)

# 2016年の3つの企業のそれぞれの車種の限界費用を取り出す
# まず適当に3つの企業の名前を選ぶ
Maker <- unique(data$Maker)
Choice_Maker <- Maker[c(8, 10, 11)]

#2016年の限界費用を抽出
Marginal_Cost_2016_Est <- Marginal_Cost_Est$Marginal_Cost_2016

# Nevo (2000) Table 4
Table_4dt <-  Marginal_Cost_2016_Est %>%
  select(Maker, Name, price, mc, margin) %>%
  # この後の表記に合わせて企業名を、Honda=Nippyo, Nissan=Brand A, Subaru=Brand B, Toyota=Brand Cとする
  mutate(Maker = 
           recode(Maker, Honda = "Nippyo", Nissan = "Brand_A", Subaru = "Brand_B", Toyota = "Brand_C"))
  
colnames(Table_4dt) <- c("Maker", "Name", "Price", "Marginal Cost", "Margin (p-mc)/p")


kable(Table_4dt, align=rep('c', 5), booktabs = T, caption = "Predicted Marginal Costs") %>%
  kable_styling(full_width = F)

Predicted Marginal Costs
Maker	Name	Price	Marginal Cost	Margin (p-mc)/p
Audi	A3シリーズ	3.280	2.3045725	0.2973864
Audi	A4シリーズ	5.180	3.5186546	0.3207230
BMW	ミニ	2.400	1.6424919	0.3156284
BMW	1シリーズ	3.100	2.1721540	0.2993052
BMW	X1	3.670	2.5684160	0.3001591
BMW	2シリーズ	3.810	2.6612243	0.3015159
BMW	3シリーズ	4.490	3.0919859	0.3113617
Daihatsu	ミラ	0.885	0.3583615	0.5950717
Daihatsu	ムーヴ	1.134	0.5779224	0.4903683
Daihatsu	ブーン	1.150	0.5919350	0.4852739
Daihatsu	キャスト	1.220	0.6530946	0.4646766
Daihatsu	タント	1.220	0.6530946	0.4646766
Daihatsu	ウェイク	1.350	0.7660195	0.4325782
Daihatsu	アトレーワゴン	1.404	0.8126625	0.4211806
Daihatsu	トール	1.463	0.8634383	0.4098166
Daihatsu	ビーゴ	1.738	1.0973459	0.3686157
Daihatsu	コペン	1.852	1.1928669	0.3559034
Daihatsu	メビウス	2.571	1.7721395	0.3107198
Daihatsu	アルティス	2.745	1.9055523	0.3058097
Fiat	500	1.998	1.3204514	0.3391134
Nippyo	N-WGN	1.164	0.6001353	0.4844198
Nippyo	N-ONE	1.185	0.6184386	0.4781109
Nippyo	N-BOX	1.270	0.6922947	0.4548861
Nippyo	フィット	1.300	0.7182716	0.4474834
Nippyo	バモス	1.374	0.7821395	0.4307573
Nippyo	シャトル	1.695	1.0554603	0.3773096
Nippyo	グレイス	1.750	1.1016307	0.3704967
Nippyo	フリード	1.880	1.2099237	0.3564236
Nippyo	ヴェゼル	1.920	1.2429991	0.3526047
Nippyo	S660	1.980	1.2923877	0.3472789
Nippyo	ステップワゴン	2.288	1.5414222	0.3263015
Nippyo	CR-V	2.470	1.6847911	0.3178984
Nippyo	ジェイド	2.530	1.7314019	0.3156514
Nippyo	CR-Z	2.700	1.8616451	0.3105018
Nippyo	オデッセイ	2.760	1.9069581	0.3090731
Nippyo	アコード	3.850	2.6700234	0.3064874
Nippyo	シビック	4.280	2.9434205	0.3122849
Nippyo	レジェンド	6.800	4.5028643	0.3378141
Lexas	CT200h	3.662	2.5550328	0.3022849
Lexas	NX	4.280	2.9488113	0.3110254
Lexas	HS250h	4.347	2.9899057	0.3121910
Lexas	IS	4.700	3.2028687	0.3185386
Lexas	RX	4.950	3.3512526	0.3229793
Lexas	RC	5.650	3.7657682	0.3334924
Lexas	GS	5.770	3.8377497	0.3348787
Lexas	LS	8.548	5.7126187	0.3317011
Lexas	RCF	9.670	6.5805658	0.3194865
Lexas	LX570	11.000	7.6618522	0.3034680
Lexas	GSF	11.110	7.7531982	0.3021424
Matsuda	キャロル	0.916	0.3884554	0.5759220
Matsuda	フレア	1.166	0.6086610	0.4779923
Matsuda	フレアクロスオーバー	1.338	0.7584250	0.4331652
Matsuda	デミオ	1.350	0.7688160	0.4305067
Matsuda	フレアワゴン	1.361	0.7783341	0.4281160
Matsuda	スクラムワゴン	1.436	0.8430522	0.4129163
Matsuda	アクセラ	1.760	1.1187613	0.3643402
Matsuda	プレマシー	1.850	1.1941276	0.3545256
Matsuda	ビアンテ	2.344	1.5968801	0.3187371
Matsuda	CX-3	2.376	1.6222765	0.3172237
Matsuda	CX-5	2.446	1.6775164	0.3141797
Matsuda	ロードスター	2.495	1.7159235	0.3122551
Matsuda	アテンザ	2.765	1.9235782	0.3043117
Matsuda	MPV	2.777	1.9326476	0.3040520
Mercedes	Aクラス	2.980	2.0819540	0.3013577
Mercedes	CLA	3.590	2.5092412	0.3010470
Mercedes	Cクラス	4.360	3.0014903	0.3115848
Mercedes	GLC	6.280	4.1583187	0.3378473
Mercedes	Eクラス	6.750	4.4588582	0.3394284
Mercedes	Sクラス	9.980	6.8432880	0.3142998
Mitsubishi	eKワゴン	1.080	0.5347760	0.5048371
Mitsubishi	デリカD:2	1.787	1.1436579	0.3600123
Mitsubishi	RVR	1.912	1.2480793	0.3472389
Mitsubishi	デリカD:5	2.408	1.6505581	0.3145523
Mitsubishi	アウトランダー	2.520	1.7385596	0.3100954
Mitsubishi	パジェロワゴン	2.927	2.0485430	0.3001220
Brand_A	クリッパーリオ	1.068	0.5196689	0.5134186
Brand_A	デイズ	1.150	0.5915111	0.4856425
Brand_A	モコ	1.193	0.6290542	0.4727123
Brand_A	ノート	1.480	0.8771352	0.4073411
Brand_A	キューブ	1.620	0.9964221	0.3849246
Brand_A	ウイングロード	1.676	1.0437898	0.3772137
Brand_A	ジューク	1.738	1.0959904	0.3693956
Brand_A	NV200バネット	1.873	1.2087339	0.3546536
Brand_A	シルフィ	1.993	1.3078368	0.3437849
Brand_A	キャラバンコーチ	2.035	1.3422637	0.3404110
Brand_A	エクストレイル	2.239	1.5074675	0.3267229
Brand_A	ラフェスタ	2.304	1.5593761	0.3231874
Brand_A	セレナワゴン	2.317	1.5697141	0.3225231
Brand_A	セドリック	2.344	1.5911381	0.3211868
Brand_A	ティアナ	2.564	1.7632683	0.3122978
Brand_A	エルグランド	3.213	2.2442067	0.3015230
Brand_A	スカイライン	3.893	2.7050769	0.3051434
Brand_A	フェアレディZ	3.910	2.7160910	0.3053476
Brand_A	フーガ	4.776	3.2543191	0.3186099
Brand_A	シーマ	7.560	5.0316863	0.3344330
Brand_A	GT-R	9.961	6.8694786	0.3103626
Brand_B	プレオ	0.807	0.2919479	0.6382305
Brand_B	ステラ	1.134	0.5809766	0.4876749
Brand_B	シフォン	1.285	0.7127770	0.4453097
Brand_B	ディアスワゴン	1.409	0.8201211	0.4179410
Brand_B	トレジア	1.430	0.8382162	0.4138348
Brand_B	ジャスティ	1.528	0.9223216	0.3963864
Brand_B	インプレッサ	1.598	0.9820437	0.3854545
Brand_B	フォレスター	2.149	1.4404556	0.3297089
Brand_B	エクシーガ	2.484	1.7074501	0.3126207
Brand_B	BRZ	2.738	1.9031568	0.3049099
Brand_B	レヴォーグ	2.776	1.9319115	0.3040665
Brand_B	レガシィ	2.938	2.0529313	0.3012487
Brand_B	WRX	3.348	2.3477896	0.2987486
Suzuki	アルト	0.778	0.2629785	0.6619813
Suzuki	ハスラー	1.079	0.5294083	0.5093529
Suzuki	ワゴンR	1.079	0.5294083	0.5093529
Suzuki	MRワゴン	1.207	0.6414934	0.4685224
Suzuki	スペーシア	1.274	0.6998433	0.4506725
Suzuki	スイフト	1.317	0.7371702	0.4402656
Suzuki	イグニス	1.382	0.7934079	0.4258987
Suzuki	ジムニー	1.407	0.8149764	0.4207701
Suzuki	エブリイワゴン	1.426	0.8313451	0.4170090
Suzuki	ソリオ	1.455	0.8562896	0.4114848
Suzuki	ジムニーシエラ	1.747	1.1046264	0.3677010
Suzuki	エクシード	2.128	1.4198991	0.3327542
Suzuki	ランディ	2.306	1.5633224	0.3220631
Suzuki	キザシ	2.867	1.9970190	0.3034465
Brand_C	ピクシスエポック	0.766	0.2399710	0.6867219
Brand_C	パッソ	1.150	0.5752309	0.4997993
Brand_C	ヴィッツ	1.154	0.5786836	0.4985411
Brand_C	bB	1.419	0.8053614	0.4324444
Brand_C	iQ	1.440	0.8231406	0.4283746
Brand_C	タンク	1.463	0.8425802	0.4240737
Brand_C	ルーミー	1.463	0.8425802	0.4240737
Brand_C	カローラ	1.485	0.8611422	0.4201063
Brand_C	ラクティス	1.615	0.9701570	0.3992836
Brand_C	シエンタ	1.690	1.0325082	0.3890484
Brand_C	ラッシュ	1.738	1.0721957	0.3830865
Brand_C	アクア	1.761	1.0911511	0.3803798
Brand_C	スペイド	1.778	1.1051357	0.3784389
Brand_C	ポルテ	1.778	1.1051357	0.3784389
Brand_C	オーリス	1.790	1.1149937	0.3770985
Brand_C	イスト	1.826	1.1445004	0.3732199
Brand_C	アリオン	1.898	1.2032055	0.3660667
Brand_C	ウィッシュ	1.906	1.2097024	0.3653188
Brand_C	プレミオ	1.909	1.2121374	0.3650406
Brand_C	アイシス	1.990	1.2776008	0.3579896
Brand_C	RAV4	2.211	1.4533281	0.3426829
Brand_C	ヴォクシー	2.254	1.4870072	0.3402808
Brand_C	ノア	2.254	1.4870072	0.3402808
Brand_C	コンフォート	2.274	1.5026131	0.3392203
Brand_C	プリウス	2.429	1.6222650	0.3321264
Brand_C	C-HR	2.516	1.6883952	0.3289367
Brand_C	トヨタ86	2.647	1.7865316	0.3250731
Brand_C	マークX	2.657	1.7939508	0.3248209
Brand_C	エスクァイア	2.658	1.7946921	0.3247960
Brand_C	ハイエースワゴン	2.706	1.8301553	0.3236677
Brand_C	ハリアー	2.798	1.8974556	0.3218529
Brand_C	アルファード	3.198	2.1796845	0.3184226
Brand_C	ヴェルファイア	3.198	2.1796845	0.3184226
Brand_C	カムリ	3.208	2.1865240	0.3184152
Brand_C	FJクルーザー	3.240	2.2083403	0.3184135
Brand_C	エスティマ	3.271	2.2293733	0.3184429
Brand_C	SAI	3.304	2.2516542	0.3185066
Brand_C	クラウン	3.812	2.5811266	0.3228944
Brand_C	ランドクルーザー	4.728	3.1260318	0.3388258
Brand_C	センチュリー	12.538	8.8289634	0.2958236
Volkswagen	ポロ	1.999	1.3201259	0.3396068
Volkswagen	ゴルフ	2.499	1.7220515	0.3109038
Volkswagen	ビートル	2.699	1.8766610	0.3046828
Volkswagen	パサート	3.486	2.4478178	0.2978147
Volvo	40シリーズ	3.240	2.2766244	0.2973382
Volvo	60シリーズ	4.540	3.1309450	0.3103645

# distribution of margin
g <- ggplot(Marginal_Cost_2016_Est, aes(x = margin)) + 
  geom_histogram(binwidth = 0.01)
g

6 合併シミュレーション

需要推定結果を用いて、限界費用関数の推定及び合併シミュレーションを行う。

留意点: 連載第二回・第三回においては、日評自動車はランダムに選んだ車種を生産・販売していると考えたが、今回第４回ではホンダの車種を生産・販売する企業として扱う。

6.1 準備

2016年のみを取り出し、日評自動車（ホンダを日評自動車とする）と国内ブランドA社（日産を国内ブランドA社とする）、日評自動車と国内ブランドB社（スバルを国内ブランドB社とする）を合併した場合を考える。加えて比較のため国内ブランドC社を用意する（トヨタを国内ブランドC社とする）。

# 最新の2016年のみ使用
data_2016 <-
  data %>%
  filter(year == 2016) %>%
  left_join(Marginal_Cost_2016_Est) %>%
  # 企業名を変更する
  mutate(Maker = 
           recode(Maker, Honda = "Nippyo", Nissan = "Brand_A", Subaru = "Brand_B", Toyota = "Brand_C")) %>%
  # 合併後の企業
  mutate(MakerNippyoA = 
           recode(Maker, Nippyo = "Nippyo_A", Brand_A = "Nippyo_A"), 
         MakerNippyoB = 
           recode(Maker, Nippyo = "Nippyo_B", Brand_B = "Nippyo_B"))

## Joining, by = c("NameID", "Maker", "Name", "price")

6.2 所有構造行列を作成する

日評自動車とブランドA社、日評自動車とブランドB社が合併した場合の所有構造行列を作成する。

# ownership matrix
J = sum(marketindex == 2016)
Ownership_NippyoA <- matrix(0, J, J)
Ownership_NippyoB <- matrix(0, J, J)
Ownership_true <- matrix(0, J, J)

for (j in 1:J) {
  for (k in 1:J) {
    if (data_2016$MakerNippyoA[j] == data_2016$MakerNippyoA[k]){
      Ownership_NippyoA[j,k] = 1
    }
    if (data_2016$MakerNippyoB[j] == data_2016$MakerNippyoB[k]){
      Ownership_NippyoB[j,k] = 1
    }
    if (data_2016$Maker[j] == data_2016$Maker[k]){
      Ownership_true[j,k] = 1
    }
  }
}

6.3 BLP推定同様の下準備

# 2016年に絞った状態で推定に必要な行列・リストを作成
datalist_2016 = list()

# マーケットとモデルの情報をGET
N <- length(data_2016$year)
T <- length(unique(data_2016$year))

# 平均効用に入ってくる部分。内生性がある価格を含んでいる。
datalist_2016$X1 <- as.data.frame(datalist$X1) %>%
  mutate(year = datalist$marketindex) %>%
  filter(year == 2016) %>%
  select(-year) %>%
  as.matrix()

# ランダム係数とInteractする部分。
datalist_2016$X2 <- as.data.frame(datalist$X2) %>%
  mutate(year = datalist$marketindex) %>%
  filter(year == 2016) %>%
  select(-year) %>%
  as.matrix()
# 操作変数行列。ここは外生変数と追加的な操作変数を含む
datalist_2016$Z <-  as.data.frame(datalist$Z) %>%
  mutate(year = datalist$marketindex) %>%
  filter(year == 2016) %>%
  select(-year) %>%
  as.matrix()

# 市場シェア
datalist_2016$ShareVec <-
  data_2016 %>%
  select(share) %>%
  as.matrix()

datalist_2016$marketindex <- data_2016$year
datalist_2016$logitshare <- data_2016$logit_share

# 乱数を用意
set.seed(111)

Nsim = 1000

datalist_2016$draw_vec <- rnorm(Nsim*ncol(datalist_2016$X2))

parameter = list()
parameter$Nsim = Nsim

marketindex <- datalist_2016$marketindex
uniquemarketindex = sort(unique(marketindex))
temp1 = matrix(rep(uniquemarketindex, N), T, N) %>% t()
temp2 = matrix(rep(marketindex, T), N, T)
mkt_denom_d = (temp1 == temp2)*1
datalist_2016$mkt_denom_d <- mkt_denom_d

6.4 反実仮想シミュレーションにおける均衡計算コード

合併後の価格を計算する。

# 推定した費用を用いる
mc = as.matrix(data_2016$mc)
# 残差として得られた観測されないXiは固定されているものとする
Xi = as.matrix(data_2016$Xi)

ここでは合併後の所有構造に基づくベルトラン競争のナッシュ均衡を計算する。方法については本誌連載4.3節を参照されたい。大枠の考え方は以下である。

ある価格ベクトル(p_oldと呼ぶ)を所与とする。
その価格のもとでの、各製品の最適反応価格を求める(この作業が以下のf_update関数)
その最適反応価格を新たな価格(p_new)とする。そのp_newに基づいて上のプロセスを再び行う。
上記手順２と３を繰り返し、価格のアップデートがほとんど起きなくなったら、プロセスを止める。この作業が以下のf_eqpriceで行われている。

# 均衡を求める

# 価格をupdateする関数
f_update <- function(datalist,p_old,Ownership,parameter,theta1,theta2,mc,Xi){
  
  # datalist内の価格を更新
  datalist$X1[,'price'] = as.matrix(p_old)
  datalist$X2 = as.matrix(p_old)
  
  X1 = datalist$X1
  
  # Market Shareを求める
  delta = (X1 %*% theta1) + Xi
  Sharevec = f_mktshare(datalist, parameter, delta, theta2)
  
  # elasticityを求める
  elas <- f_elasticity(datalist, parameter, theta1, theta2, delta)
  elasmat <- elas$elasmat_2016
  
  J = length(p_old)
  Derivative = - (elasmat)*(kronecker(matrix(1,J,1), t(Sharevec)))/(kronecker(matrix(1,1,J), p_old))
  Delta = Ownership*Derivative
  p_new = mc + (solve(Delta) %*% Sharevec)
  return(p_new)
}

# iterationで均衡価格を求める関数
f_eqprice = function(datalist,p_ini,Ownership,parameter,theta1,theta2,mc,Xi){
  
  # set the threshold
  lambda <- 1e-6
  p_old <- p_ini
  distance <- 10000
  while (distance > lambda) {
    p_new <- f_update(datalist,p_old,Ownership,parameter,theta1,theta2,mc,Xi)
    distance <- max(abs(p_new - p_old))
    p_old <- p_new
  }
  return(p_new)
}

6.5 合併シミュレーション

まず、日評自動車とA社が合併したケースのシミュレーションを行う。

# Nippyo-Brand A Simulation

# set the initial price
p_ini <- data_2016$price
p_NippyoA <- f_eqprice(datalist_2016,p_ini,Ownership_NippyoA,parameter,theta1_hat,theta2_hat,mc,Xi)
data_2016 <- data_2016 %>%
  mutate(p_NippyoA = p_NippyoA)

つづいて、日評自動車とB社が合併したケースのシミュレーションを行う

# Nippyo-Brand B Simulation
# set the initial price
p_ini <- data_2016$price
p_NippyoB <- f_eqprice(datalist_2016,p_ini,Ownership_NippyoB,parameter,theta1_hat,theta2_hat,mc,Xi)
data_2016 <- data_2016 %>%
  mutate(p_NippyoB = p_NippyoB)

# シミュレートされた価格でのシェアを与える関数
f_mktshare_sim <- function(datalist,p,parameter,theta1,theta2, Xi){
  
  # datalist内の価格を更新
  datalist$X1[,'price'] = p
  datalist$X2 = p
  X1 = datalist$X1
  
  # Market Shareを求める
  delta = (X1 %*% theta1) + Xi
  Sharevec = f_mktshare(datalist, parameter, delta, theta2)
  
  return(Sharevec)
}

合併後の価格からそれぞれの車種のマーケットシェアを計算する。

# 合併時のシェアを計算
share_NippyoA <- f_mktshare_sim(datalist_2016,data_2016$p_NippyoA,parameter,theta1_hat,theta2_hat, Xi)
share_NippyoB <- f_mktshare_sim(datalist_2016,data_2016$p_NippyoB,parameter,theta1_hat,theta2_hat, Xi)

data_2016 <- data_2016 %>%
  mutate(share_NippyoA = share_NippyoA, share_NippyoB = share_NippyoB)

7 シミュレーション結果のまとめ

7.1 合併シミュレーションによる価格・販売台数変化

なお、この表はNevo(2000, RAND)のTABLE 5と同様。（比較対象として国内ブランドC社（トヨタを国内ブランドC社とする）を表に含む）

# 表を作るために専用のデータフレームを作る
Table_5dt <- data_2016 %>%
  mutate(p_NippyoA_PPC = (p_NippyoA - price)/price*100,
         share_NippyoA_PPC = (share_NippyoA - share)/share*100,
         p_NippyoB_PPC = (p_NippyoB - price)/price*100,
         share_NippyoB_PPC = (share_NippyoB - share)/share*100) %>%
  select(Maker, Name, p_NippyoA_PPC, share_NippyoA_PPC, p_NippyoB_PPC,
         share_NippyoB_PPC) %>%
  filter(Maker %in% c("Nippyo","Brand_A","Brand_B","Brand_C"))

# 名前を付けなおす
colnames(Table_5dt) <- c("Maker", "Name", "p", "q", "p", "q")

# 表を作る
kable(Table_5dt,  align=rep('c', 5), booktabs = T, caption = "Predicted Percent Change in Prices and Quantities as a Result of Mergers") %>%
  add_header_above(c(" ", " ", "Nippyo and Brand A" = 2, "Nippyo and Brand B" = 2)) %>%
  kable_styling(full_width = F)

Predicted Percent Change in Prices and Quantities as a Result of Mergers
		Nippyo and Brand A		Nippyo and Brand B
Maker	Name	p	q	p	q
Nippyo	N-WGN	0.5650051090	-1.1305656	2.147882e-01	-0.4301405
Nippyo	N-ONE	0.5619195185	-1.1393250	2.139901e-01	-0.4342648
Nippyo	N-BOX	0.5512648510	-1.1751910	2.114143e-01	-0.4512262
Nippyo	フィット	0.5481416775	-1.1880070	2.107336e-01	-0.4573151
Nippyo	バモス	0.5416890287	-1.2199697	2.095063e-01	-0.4725626
Nippyo	シャトル	0.5299517358	-1.3643188	2.101514e-01	-0.5424223
Nippyo	グレイス	0.5300531854	-1.3899629	2.110545e-01	-0.5549865
Nippyo	フリード	0.5322736158	-1.4515992	2.139518e-01	-0.5853509
Nippyo	ヴェゼル	0.5334750206	-1.4708471	2.150451e-01	-0.5948782
Nippyo	S660	0.5357023972	-1.4999615	2.168525e-01	-0.6093275
Nippyo	ステップワゴン	0.5543577778	-1.6536920	2.290282e-01	-0.6862994
Nippyo	CR-V	0.5703964399	-1.7474591	2.382808e-01	-0.7337201
Nippyo	ジェイド	0.5764125642	-1.7787539	2.416361e-01	-0.7496138
Nippyo	CR-Z	0.5952660306	-1.8681980	2.519084e-01	-0.7952041
Nippyo	オデッセイ	0.6025202014	-1.8999644	2.557896e-01	-0.8114498
Nippyo	アコード	0.7680315969	-2.4447195	3.407087e-01	-1.0943499
Nippyo	シビック	0.8318621352	-2.5980688	3.731327e-01	-1.1769169
Nippyo	レジェンド	0.8083368861	-2.3283768	3.703539e-01	-1.0759421
Brand_A	クリッパーリオ	0.9548086364	-1.8127766	2.333679e-03	0.0143496
Brand_A	デイズ	0.9262780206	-1.8591021	2.275363e-03	0.0147666
Brand_A	モコ	0.9134980357	-1.8835542	2.249356e-03	0.0149868
Brand_A	ノート	0.8566126899	-2.0493505	2.132949e-03	0.0164795
Brand_A	キューブ	0.8420739912	-2.1316881	2.100765e-03	0.0172193
Brand_A	ウイングロード	0.8380774711	-2.1648544	2.090844e-03	0.0175167
Brand_A	ジューク	0.8347184766	-2.2017136	2.081415e-03	0.0178469
Brand_A	NV200バネット	0.8308680085	-2.2824193	2.065292e-03	0.0185683
Brand_A	シルフィ	0.8309510270	-2.3545796	2.054452e-03	0.0192114
Brand_A	キャラバンコーチ	0.8316773853	-2.3799082	2.051139e-03	0.0194367
Brand_A	エクストレイル	0.8397634872	-2.5032415	2.036203e-03	0.0205316
Brand_A	ラフェスタ	0.8437842454	-2.5425708	2.031193e-03	0.0208802
Brand_A	セレナワゴン	0.8446654685	-2.5504342	2.030144e-03	0.0209499
Brand_A	セドリック	0.8465753968	-2.5667610	2.027901e-03	0.0210946
Brand_A	ティアナ	0.8658721956	-2.6992971	2.004941e-03	0.0222710
Brand_A	エルグランド	0.9524601030	-3.0723311	1.838942e-03	0.0257051
Brand_A	スカイライン	1.0617429658	-3.3801901	1.405422e-03	0.0292232
Brand_A	フェアレディZ	1.0643065159	-3.3859929	1.390698e-03	0.0293091
Brand_A	フーガ	1.1506568873	-3.5054543	4.664819e-04	0.0333048
Brand_A	シーマ	0.8564932533	-2.4889525	-1.159927e-03	0.0358299
Brand_A	GT-R	0.6174407974	-1.9303320	-9.647124e-04	0.0327698
Brand_B	プレオ	0.0054652730	0.0341278	1.094062e+00	-1.6867412
Brand_B	ステラ	0.0046228092	0.0379398	9.266954e-01	-1.8681757
Brand_B	シフォン	0.0044007132	0.0397265	8.854287e-01	-1.9540019
Brand_B	ディアスワゴン	0.0042619366	0.0412014	8.616408e-01	-2.0253316
Brand_B	トレジア	0.0042414254	0.0414517	8.583441e-01	-2.0374819
Brand_B	ジャスティ	0.0041549928	0.0426205	8.453678e-01	-2.0944332
Brand_B	インプレッサ	0.0041012101	0.0434560	8.382803e-01	-2.1353487
Brand_B	フォレスター	0.0037976896	0.0499942	8.281672e-01	-2.4621735
Brand_B	エクシーガ	0.0036278399	0.0538840	8.498397e-01	-2.6617994
Brand_B	BRZ	0.0034662857	0.0567689	8.760657e-01	-2.8108881
Brand_B	レヴォーグ	0.0034383552	0.0571953	8.805906e-01	-2.8328679
Brand_B	レガシィ	0.0033058434	0.0589972	9.013808e-01	-2.9252063
Brand_B	WRX	0.0028541204	0.0634384	9.621136e-01	-3.1438898
Brand_C	ピクシスエポック	0.0075525989	0.0307629	3.405698e-03	0.0117937
Brand_C	パッソ	0.0063939416	0.0343873	2.951100e-03	0.0133609
Brand_C	ヴィッツ	0.0063872622	0.0344250	2.948713e-03	0.0133774
Brand_C	bB	0.0060788181	0.0369090	2.849757e-03	0.0144723
Brand_C	iQ	0.0060633597	0.0371044	2.845905e-03	0.0145592
Brand_C	タンク	0.0060476115	0.0373182	2.842216e-03	0.0146543
Brand_C	ルーミー	0.0060476115	0.0373182	2.842216e-03	0.0146543
Brand_C	カローラ	0.0060336496	0.0375223	2.839181e-03	0.0147453
Brand_C	ラクティス	0.0059705357	0.0387219	2.829981e-03	0.0152822
Brand_C	シエンタ	0.0059469205	0.0394083	2.830475e-03	0.0155913
Brand_C	ラッシュ	0.0059358703	0.0398451	2.832646e-03	0.0157887
Brand_C	アクア	0.0059315806	0.0400538	2.834148e-03	0.0158831
Brand_C	スペイド	0.0059288007	0.0402077	2.835437e-03	0.0159529
Brand_C	ポルテ	0.0059288007	0.0402077	2.835437e-03	0.0159529
Brand_C	オーリス	0.0059270303	0.0403161	2.836436e-03	0.0160021
Brand_C	イスト	0.0059226241	0.0406408	2.839850e-03	0.0161496
Brand_C	アリオン	0.0059175016	0.0412864	2.848387e-03	0.0164439
Brand_C	ウィッシュ	0.0059172071	0.0413579	2.849463e-03	0.0164765
Brand_C	プレミオ	0.0059171098	0.0413847	2.849873e-03	0.0164888
Brand_C	アイシス	0.0059169689	0.0421042	2.862097e-03	0.0168184
Brand_C	RAV4	0.0059342942	0.0440335	2.903792e-03	0.0177101
Brand_C	ヴォクシー	0.0059394929	0.0444031	2.912774e-03	0.0178822
Brand_C	ノア	0.0059394929	0.0444031	2.912774e-03	0.0178822
Brand_C	コンフォート	0.0059420132	0.0445744	2.917000e-03	0.0179622
Brand_C	プリウス	0.0059621892	0.0458880	2.950056e-03	0.0185785
Brand_C	C-HR	0.0059726134	0.0466153	2.968158e-03	0.0189223
Brand_C	トヨタ86	0.0059844839	0.0476980	2.993501e-03	0.0194376
Brand_C	マークX	0.0059851251	0.0477800	2.995302e-03	0.0194768
Brand_C	エスクァイア	0.0059851868	0.0477882	2.995481e-03	0.0194807
Brand_C	ハイエースワゴン	0.0059875960	0.0481811	3.003789e-03	0.0196690
Brand_C	ハリアー	0.0059887816	0.0489298	3.017968e-03	0.0200293
Brand_C	アルファード	0.0059185694	0.0521396	3.040585e-03	0.0216007
Brand_C	ヴェルファイア	0.0059185694	0.0521396	3.040585e-03	0.0216007
Brand_C	カムリ	0.0059148287	0.0522194	3.040106e-03	0.0216404
Brand_C	FJクルーザー	0.0059020982	0.0524746	3.038167e-03	0.0217674
Brand_C	エスティマ	0.0058886371	0.0527219	3.035685e-03	0.0218908
Brand_C	SAI	0.0058730551	0.0529851	3.032371e-03	0.0220224
Brand_C	クラウン	0.0054511909	0.0570609	2.881507e-03	0.0241066
Brand_C	ランドクルーザー	0.0038487105	0.0642479	2.114705e-03	0.0280192
Brand_C	センチュリー	0.0001387364	0.0604210	-4.030637e-05	0.0275315

7.2 日評自動車とブランドAの合併後に価格が変化しないような限界費用を計算する

ここでは、需要推定値・シェア・価格はデータ上のものに固定した上で、「日評自動車とA社が合併していて共同利潤最大化していた場合」における限界費用を計算している。

このようにして得られた限界費用の元では、「日評自動車とA社」が合併した場合に利潤最大化するような価格というものは、データ上の価格と一致（すなわち価格が変化しない）することとなる。

# 合併前の弾力性行列
elasmat_2016 = list(elasticity$elasmat_2016)

# 日評自動車とブランドA社の名前を同じにする（所有構造を変更）
data_2016_NippyoA <- data_2016 %>%
  mutate(Maker = MakerNippyoA) %>%
  select(year, Maker, Name, NameID, price, share)

# 限界費用を計算（シェアと価格、弾力性は合併以前のものを使用）
mc_NippyoA_pfix <- f_Merginal_Cost_Est(data_2016_NippyoA, marketindex, elasmat_2016)
mc_NippyoA_pfix_2016 <- mc_NippyoA_pfix$Marginal_Cost_2016
data_2016 <- data_2016 %>%
  mutate(mc_NippyoA_pfix = mc_NippyoA_pfix_2016$mc)

7.3 日評自動車とブランドBの合併後に価格が変化しないような限界費用を計算する

考え方は上と同様。

# 日評自動車とブランドB社の名前を同じにする（所有構造を変更）
data_2016_NippyoB <- data_2016 %>%
  mutate(Maker = MakerNippyoB) %>%
  select(year, Maker, Name, NameID, price, share)

# 限界費用を計算（シェアと価格、弾力性は合併以前のものを使用）
mc_NippyoB_pfix <- f_Merginal_Cost_Est(data_2016_NippyoB, marketindex, elasmat_2016)
mc_NippyoB_pfix_2016 <- mc_NippyoB_pfix$Marginal_Cost_2016
data_2016 <- data_2016 %>%
  mutate(mc_NippyoB_pfix = mc_NippyoB_pfix_2016$mc)

以上の結果に基づき、合併による価格上昇を打ち消すのに必要な限界費用変化を求めることができる。

# 表を作るために専用のデータフレームを作る
Table_6dt <- data_2016 %>%
  mutate(mc_NippyoA_PPC = (mc_NippyoA_pfix - mc)/mc*100,
         mc_NippyoB_PPC = (mc_NippyoB_pfix - mc)/mc*100) %>%
  select(Maker, Name, mc_NippyoA_PPC, mc_NippyoB_PPC) %>%
  filter(Maker %in% c("Nippyo","Brand_A","Brand_B","Brand_C"))

# 名前を付けなおす
colnames(Table_6dt) <- c("Origin Maker", "Name", "Nippyo and Brand A",
                         "Nippyo and Brand B")

# 表を作る
kable(Table_6dt,  align=rep('c', 5), booktabs = T, 
      caption = "Percent Reduction in Marginal Costs Required for No Change in Predicted Postmerger Prices")

Percent Reduction in Marginal Costs Required for No Change in Predicted Postmerger Prices
Origin Maker	Name	Nippyo and Brand A	Nippyo and Brand B
Nippyo	N-WGN	-0.9800274	-0.3735398
Nippyo	N-ONE	-0.9615378	-0.3671281
Nippyo	N-BOX	-0.8978292	-0.3451965
Nippyo	フィット	-0.8788993	-0.3387413
Nippyo	バモス	-0.8384657	-0.3250852
Nippyo	シャトル	-0.7304423	-0.2903292
Nippyo	グレイス	-0.7190931	-0.2869900
Nippyo	フリード	-0.6976706	-0.2810869
Nippyo	ヴェゼル	-0.6923789	-0.2797510
Nippyo	S660	-0.6854216	-0.2781109
Nippyo	ステップワゴン	-0.6642920	-0.2751378
Nippyo	CR-V	-0.6603191	-0.2765856
Nippyo	ジェイド	-0.6600507	-0.2774577
Nippyo	CR-Z	-0.6616074	-0.2808090
Nippyo	オデッセイ	-0.6628753	-0.2822648
Nippyo	アコード	-0.7232485	-0.3223134
Nippyo	シビック	-0.7541154	-0.3398870
Nippyo	レジェンド	-0.8175309	-0.3749426
Brand_A	クリッパーリオ	-1.7455320	0.0000000
Brand_A	デイズ	-1.5939844	0.0000000
Brand_A	モコ	-1.5292960	0.0000000
Brand_A	ノート	-1.2506234	0.0000000
Brand_A	キューブ	-1.1715657	0.0000000
Brand_A	ウイングロード	-1.1461951	0.0000000
Brand_A	ジューク	-1.1214295	0.0000000
Brand_A	NV200バネット	-1.0774240	0.0000000
Brand_A	シルフィ	-1.0473718	0.0000000
Brand_A	キャラバンコーチ	-1.0384917	0.0000000
Brand_A	エクストレイル	-1.0049800	0.0000000
Brand_A	ラフェスタ	-0.9970935	0.0000000
Brand_A	セレナワゴン	-0.9956537	0.0000000
Brand_A	セドリック	-0.9928014	0.0000000
Brand_A	ティアナ	-0.9756494	0.0000000
Brand_A	エルグランド	-0.9667036	0.0000000
Brand_A	スカイライン	-0.9880065	0.0000000
Brand_A	フェアレディZ	-0.9886698	0.0000000
Brand_A	フーガ	-1.0176416	0.0000000
Brand_A	シーマ	-0.9189211	0.0000000
Brand_A	GT-R	-0.7309373	0.0000000
Brand_B	プレオ	0.0000000	-2.7015019
Brand_B	ステラ	0.0000000	-1.5883454
Brand_B	シフォン	0.0000000	-1.3888272
Brand_B	ディアスワゴン	0.0000000	-1.2773524
Brand_B	トレジア	0.0000000	-1.2616996
Brand_B	ジャスティ	0.0000000	-1.1981654
Brand_B	インプレッサ	0.0000000	-1.1607702
Brand_B	フォレスター	0.0000000	-0.9997983
Brand_B	エクシーガ	0.0000000	-0.9632454
Brand_B	BRZ	0.0000000	-0.9509910
Brand_B	レヴォーグ	0.0000000	-0.9499827
Brand_B	レガシィ	0.0000000	-0.9476066
Brand_B	WRX	0.0000000	-0.9519346
Brand_C	ピクシスエポック	0.0000000	0.0000000
Brand_C	パッソ	0.0000000	0.0000000
Brand_C	ヴィッツ	0.0000000	0.0000000
Brand_C	bB	0.0000000	0.0000000
Brand_C	iQ	0.0000000	0.0000000
Brand_C	タンク	0.0000000	0.0000000
Brand_C	ルーミー	0.0000000	0.0000000
Brand_C	カローラ	0.0000000	0.0000000
Brand_C	ラクティス	0.0000000	0.0000000
Brand_C	シエンタ	0.0000000	0.0000000
Brand_C	ラッシュ	0.0000000	0.0000000
Brand_C	アクア	0.0000000	0.0000000
Brand_C	スペイド	0.0000000	0.0000000
Brand_C	ポルテ	0.0000000	0.0000000
Brand_C	オーリス	0.0000000	0.0000000
Brand_C	イスト	0.0000000	0.0000000
Brand_C	アリオン	0.0000000	0.0000000
Brand_C	ウィッシュ	0.0000000	0.0000000
Brand_C	プレミオ	0.0000000	0.0000000
Brand_C	アイシス	0.0000000	0.0000000
Brand_C	RAV4	0.0000000	0.0000000
Brand_C	ヴォクシー	0.0000000	0.0000000
Brand_C	ノア	0.0000000	0.0000000
Brand_C	コンフォート	0.0000000	0.0000000
Brand_C	プリウス	0.0000000	0.0000000
Brand_C	C-HR	0.0000000	0.0000000
Brand_C	トヨタ86	0.0000000	0.0000000
Brand_C	マークX	0.0000000	0.0000000
Brand_C	エスクァイア	0.0000000	0.0000000
Brand_C	ハイエースワゴン	0.0000000	0.0000000
Brand_C	ハリアー	0.0000000	0.0000000
Brand_C	アルファード	0.0000000	0.0000000
Brand_C	ヴェルファイア	0.0000000	0.0000000
Brand_C	カムリ	0.0000000	0.0000000
Brand_C	FJクルーザー	0.0000000	0.0000000
Brand_C	エスティマ	0.0000000	0.0000000
Brand_C	SAI	0.0000000	0.0000000
Brand_C	クラウン	0.0000000	0.0000000
Brand_C	ランドクルーザー	0.0000000	0.0000000
Brand_C	センチュリー	0.0000000	0.0000000

7.4 合併シミュレーションの厚生分析

ある個人\(i\)の補償変分（CV）は以下の式で求められる。

\[ CV_i = \frac{\ln(\sum^J_{j=0}\exp(V^{post}_{ij})) - \ln(\sum^J_{j=0}\exp(V^{pre}_{ij}))}{\alpha_i} \] まず、消費者余剰を推定する関数を用意する。

# Nevo (2000) 式(6), (7)参照
# p.414に従いdrawの中でCV_iの平均を取ってから消費者数を乗する

# consumer surplusを計算する関数
# see Small and Rosen (1981), ignoring an unknown constant term
# see also Tran (2009) "Discrete Choice Methods with Simulation" Chapter 3
f_CS = function(datalist, p, parameter, theta1, theta2, Xi, HH){
  # datalist内の価格を更新
  datalist$X1[,'price'] = p
  datalist$X2 = as.matrix(p)
  X1 <- datalist$X1
  
  # deltaを求める
  delta = (X1 %*% theta1) + Xi
  
  # f_mktshareの分子と同じ計算
  X2 <- datalist$X2
  N <- length(datalist$marketindex)
  draw_vec <- datalist$draw_vec
  marketindex <- datalist$marketindex
  mkt_denom_d <- datalist$mkt_denom_d
  Nsim <- parameter$Nsim
  K2 <- length(theta2)
  # Nonlinear component. mu (N, ns) matrix.
  mu <- X2 %*% diag(x = theta2, nrow = K2) %*% matrix(draw_vec, nrow = K2)
  
  # prepare exp_V
  V <- delta %*% matrix(1, nrow = 1, ncol = Nsim) + mu
  exp_V <- exp(V)
  # numerator
  # the second term is for the outside good
  numerator <- log(colSums(exp_V) + exp(matrix(0, 1, Nsim)))
  
  # Price parameter alpha_i
  draw_for_price <- matrix(draw_vec, nrow = K2)[1, ]
  # theta1のなかで2つ目、theta2のなかで1つ目に価格係数がある
  alpha_i <- - (theta1[2] + theta2[1]*draw_for_price)
  
  # Consumer Surplus
  CS <- mean(numerator/alpha_i) * HH
  return(CS)
}

消費者余剰と補償変分を求める。

# compensating variation
HH_2016 <- unique(data_2016$HH)
CS_2016 <- f_CS(datalist_2016, data_2016$price, parameter, theta1_hat, theta2_hat, Xi, HH_2016)
CS_NippyoA <- f_CS(datalist_2016, data_2016$p_NippyoA, parameter, theta1_hat, theta2_hat, Xi, HH_2016)
CS_NippyoB <- f_CS(datalist_2016, data_2016$p_NippyoB, parameter, theta1_hat, theta2_hat, Xi, HH_2016)

CV_NippyoA <- CS_NippyoA - CS_2016
CV_NippyoB <- CS_NippyoB - CS_2016

次に企業の利潤を計算する。ある企業\(f\)の利潤の式は

\[ \sum_{j \in \mathcal{F}_f} (p_j - mc_j)q_j(\mathbf{p}) \] である。ここでpは価格、mcは限界費用、qは数量（ここではシェアに家計の数を掛けることで求める）である。

# Profit and Revenue
f_profit <-  function(Maker, price, mc, share, HH){
  
  pro_rev_dt <- data.frame(Maker, 
                           profit_each = (price - mc)*share*HH,
                           revenue_each = price*share*HH)
  pro_rev <- pro_rev_dt %>%
    group_by(Maker) %>%
    summarise(profit = sum(profit_each), revenue = sum(revenue_each)) %>%
    ungroup() %>%
    as.data.frame()
  return(pro_rev)
}

総余剰の変化を導出する。

pro_rev_2016 <- f_profit(data_2016$Maker, data_2016$price, data_2016$mc, data_2016$share, HH_2016)  
pro_rev_NippyoA <- f_profit(data_2016$Maker, data_2016$p_NippyoA, data_2016$mc,
                            data_2016$share_NippyoA, HH_2016) 
pro_rev_NippyoB <- f_profit(data_2016$Maker, data_2016$p_NippyoB, data_2016$mc,
                            data_2016$share_NippyoB, HH_2016) 

# Total Surplus Change
TS_Chage_NippyoA <- CV_NippyoA + sum(pro_rev_NippyoA$profit - pro_rev_2016$profit)
TS_Chage_NippyoB <- CV_NippyoB + sum(pro_rev_NippyoB$profit - pro_rev_2016$profit)

消費者余剰と総余剰の変化の表を作成

Table_7dt1 <- tibble(
  clo1 = c("Consumer surplus", "Total Welfare"),
  col2 = c(CV_NippyoA, TS_Chage_NippyoA),
  col3 = c(CV_NippyoB, TS_Chage_NippyoB)
)

colnames(Table_7dt1) <- c("Measure", "Nippyo and Brand A",
                          "Nippyo and Brand B")

kable(Table_7dt1,  align=rep('c', 3), booktabs = T, 
      caption = "Consumer Surplus and Welfare")

Consumer Surplus and Welfare
Measure	Nippyo and Brand A	Nippyo and Brand B
Consumer surplus	-12527.56	-5106.183
Total Welfare	-11170.69	-4497.038

表を作る前の下準備をする関数を作成

gen_pro_rev_table <- function(pro_rev_2016, pro_rev_NippyoA, pro_rev_NippyoB){
  result_df <- tibble(
    Maker = pro_rev_2016$Maker,
    col2 = pro_rev_NippyoA$profit - pro_rev_2016$profit,
    col3 = pro_rev_NippyoA$revenue - pro_rev_2016$revenue,
    col4 = pro_rev_NippyoB$profit - pro_rev_2016$profit,
    col5 = pro_rev_NippyoB$revenue - pro_rev_2016$revenue,
  ) 
  Total <- tibble(Maker = "Total",
                  col2 = sum(result_df$col2),
                  col3 = sum(result_df$col3),
                  col4 = sum(result_df$col4),
                  col5 = sum(result_df$col5))
  result_df <- rbind(result_df, Total)
  
  colnames(result_df) <- c("Maker", "Profits", "Revenues", "Profits", "Revenues")
  
  return(result_df)
}

合併前と合併後の収入を比較する表を作成（表中の企業名は日評自動車としているホンダ、ブランドA社としている日産、ブランドB社としているスバル、ブランドC社としているトヨタ以外はそのままの社名を使っている。）

Table_7dt2 <- gen_pro_rev_table(pro_rev_2016, pro_rev_NippyoA, pro_rev_NippyoB)

# 表を作る
kable(Table_7dt2,  align=rep('c', 5), booktabs = T, caption = "Change in Profits and Revenues") %>%
  add_header_above(c(" ", "Nippyo and Brand A" = 2, "Nippyo and Brand B" = 2)) %>%
  kable_styling(full_width = F)

Change in Profits and Revenues
	Nippyo and Brand A		Nippyo and Brand B
Maker	Profits	Revenues	Profits	Revenues
Audi	13.059879	41.092100	5.8799765	18.152078
BMW	43.267742	130.362381	19.1943987	56.435463
Brand_A	16.155150	-10640.014691	64.7961967	153.373590
Brand_B	62.024083	163.371990	-13.3922039	-5019.738372
Brand_C	640.418088	1501.995588	274.5824215	625.040629
Daihatsu	119.923532	223.596837	49.0816713	89.986153
Fiat	1.708678	4.434631	0.7248364	1.834694
Lexas	66.344666	208.280302	30.1277171	93.484230
Matsuda	66.515353	162.321187	28.0901545	67.108239
Mercedes	57.825532	180.788053	26.1795950	80.967982
Mitsubishi	22.875957	53.839981	9.6180951	22.235082
Nippyo	88.338898	-9404.628796	48.5061687	-3779.899786
Suzuki	125.547334	231.517649	51.3571936	93.212884
Volkswagen	22.106758	62.784761	9.5756474	26.486001
Volvo	10.759187	33.699215	4.8230457	14.759988
Total	1356.870838	-17046.558811	609.1449142	-7456.561145

7.5 シミュレーションに関する計算時間

最後に、ここまでの計算にかかった合計時間をレポートする。

tictoc::toc()

## 65.53 sec elapsed

8 Appendix 1: 追加的なシミュレーション

ここでは「総余剰が変化しない」ために要求される合併企業の限界費用削減について計算しよう。

なお、Appendix 1の分析については比較的計算時間がかかる(パソコンのスペックにもよるが10-20分程度)ので、その点留意されたい。ここでもtictocで計算時間を計測しよう。

tictoc::tic()

8.1 関数の定義：総余剰が変化しないような限界費用削減

# 特定の企業のみ限界費用を定数倍した時の、総余剰の変化を関数
f_effect_cost_reduction <- function(cost_red, cost_red_firm, Ownership, data, datalist, 
                                    parameter, theta1, theta2, HH, p_pre, pro_rev_pre, CS_pre){
  # 特定の企業のみ、限界費用を一律に係数倍
  data <- data %>%
    mutate(mc = if_else(Maker %in% cost_red_firm, mc * cost_red, mc))
  mc = data$mc
  
  # 均衡価格
  p_post <- f_eqprice(datalist,p_pre,Ownership,parameter,theta1,theta2,mc,Xi)
  
  # CVの計算
  CV <- (f_CS(datalist, p_post, parameter, theta1, theta2, Xi, HH) - CS_pre)
  # 利益と売上の計算
  share_post <- f_mktshare_sim(datalist,p_post,parameter,theta1,theta2, Xi)
  pro_rev_post <- f_profit(data$Maker, p_post, data$mc, share_post, HH) 
  # 総余剰の変化
  TS_Chage <- CV + sum(pro_rev_post$profit - pro_rev_pre$profit)
  obj <- TS_Chage
  return(obj)
}

8.2 日評自動車とブランドAについて限界費用削減を計算

# Nippyo-Brand A 限界費用削減
# cost_redについて総余剰変化は単調減少するため、二分法を適用

# 両企業が限界費用削減 （日評自動車とbrand Aを取り出す）
cost_red_firm = c("Nippyo","Brand_A")

distance <- 100
lambda <- 10^(-6)
max_cost_red <- 1
min_cost_red <- 0

# ここのループにおいて、総余剰を合併前の水準に保つために必要な限界費用削減率を計算している。
while (distance > lambda){
  mid_cost_red = (max_cost_red + min_cost_red)/2

  mid_eval = f_effect_cost_reduction(
    cost_red = mid_cost_red, 
    cost_red_firm = cost_red_firm, Ownership = Ownership_NippyoA, data = data_2016, 
    datalist = datalist_2016,  
    parameter = parameter, theta1 = theta1_hat, theta2 = theta2_hat, HH = HH_2016, 
    p_pre = data_2016$p_NippyoA, pro_rev_pre = pro_rev_2016, CS_pre = CS_2016)

  if (mid_eval > 0){
    min_cost_red <- mid_cost_red
  }else{
    max_cost_red <- mid_cost_red
  }
  
  distance = abs(mid_eval-0)
  print(distance)
}

## [1] 2067673
## [1] 743127.3
## [1] 316809.1
## [1] 142416.9
## [1] 63221.5
## [1] 25447.89
## [1] 6996.972
## [1] 2121.935
## [1] 2428.707
## [1] 151.1891
## [1] 985.9214
## [1] 417.5034
## [1] 133.1914
## [1] 8.990259
## [1] 62.10274
## [1] 26.55677
## [1] 8.783391
## [1] 0.1034005
## [1] 4.340004
## [1] 2.118304
## [1] 1.007452
## [1] 0.452026
## [1] 0.1743128
## [1] 0.03545614
## [1] 0.03397217
## [1] 0.000741987
## [1] 0.01661509
## [1] 0.007936553
## [1] 0.00359728
## [1] 0.001427646
## [1] 0.0003428299
## [1] 0.0001995791
## [1] 7.162669e-05
## [1] 6.397752e-05
## [1] 3.824866e-06
## [1] 3.007801e-05
## [1] 1.312723e-05
## [1] 4.651383e-06
## [1] 4.135545e-07

cost_red_NippyoA <- mid_cost_red

data_2016 <- data_2016 %>%
  mutate(mc_NippyoA_TSfix = if_else(Maker %in% cost_red_firm, mc * cost_red_NippyoA, mc))

p_NippyoA_TSfix <- f_eqprice(datalist_2016,data_2016$p_NippyoA,Ownership_NippyoA,
                             parameter,theta1_hat,theta2_hat,data_2016$mc_NippyoA_TSfix,Xi)
share_NippyoA_TSfix <- f_mktshare_sim(datalist_2016,p_NippyoA_TSfix,parameter,theta1_hat,theta2_hat, Xi)
data_2016 <- data_2016 %>%
  mutate(p_NippyoA_TSfix = p_NippyoA_TSfix, share_NippyoA_TSfix = share_NippyoA_TSfix)

8.3 日評自動車とブランドBについて限界費用削減を計算

# Nippyo-Brand B 限界費用削減
# cost_redについて総余剰変化は単調減少するため、二分法を適用

# 両企業が限界費用削減 （日評自動車とbrand Bを取り出す）
cost_red_firm = c("Nippyo","Brand_B")

distance <- 100
lambda <- 10^(-6)
max_cost_red <- 1
min_cost_red <- 0

# ここのループにおいて、総余剰を合併前の水準に保つために必要な限界費用削減率を計算している。
while (distance > lambda){
  mid_cost_red = (max_cost_red + min_cost_red)/2
  mid_eval = f_effect_cost_reduction(
    cost_red = mid_cost_red, 
    cost_red_firm = cost_red_firm, Ownership = Ownership_NippyoB, data = data_2016, 
    datalist = datalist_2016,  
    parameter = parameter, theta1 = theta1_hat, theta2 = theta2_hat, HH = HH_2016, 
    p_pre = data_2016$p_NippyoB, pro_rev_pre = pro_rev_2016, CS_pre = CS_2016)
  if (mid_eval > 0){
    min_cost_red <- mid_cost_red
  }else{
    max_cost_red <- mid_cost_red
  }
  
  distance = abs(mid_eval-0)
  print(distance)
}

## [1] 1633622
## [1] 584455
## [1] 250588.8
## [1] 114752.8
## [1] 53217.74
## [1] 23901.56
## [1] 9589.808
## [1] 2518.543
## [1] 996.1794
## [1] 759.4478
## [1] 118.7974
## [1] 320.2167
## [1] 100.6826
## [1] 9.064212
## [1] 45.80748
## [1] 18.37121
## [1] 4.653394
## [1] 2.205436
## [1] 1.223972
## [1] 0.4907334
## [1] 0.366619
## [1] 0.06205731
## [1] 0.1522808
## [1] 0.04511175
## [1] 0.008472777
## [1] 0.01831949
## [1] 0.004923353
## [1] 0.00177471
## [1] 0.001574322
## [1] 0.0001001958
## [1] 0.0007370655
## [1] 0.0003184348
## [1] 0.0001091182
## [1] 4.462312e-06
## [1] 4.786623e-05
## [1] 2.170176e-05
## [1] 8.619282e-06
## [1] 2.078166e-06
## [1] 1.190073e-06
## [1] 4.442077e-07

cost_red_NippyoB <- mid_cost_red

data_2016 <- data_2016 %>%
  mutate(mc_NippyoB_TSfix = if_else(Maker %in% cost_red_firm, mc * cost_red_NippyoB, mc))

p_NippyoB_TSfix <- f_eqprice(datalist_2016,data_2016$p_NippyoB,Ownership_NippyoB,
                             parameter,theta1_hat,theta2_hat,data_2016$mc_NippyoB_TSfix,Xi)
share_NippyoB_TSfix <- f_mktshare_sim(datalist_2016,p_NippyoB_TSfix,parameter,theta1_hat,theta2_hat, Xi)
data_2016 <- data_2016 %>%
  mutate(p_NippyoB_TSfix = p_NippyoB_TSfix, share_NippyoB_TSfix= share_NippyoB_TSfix)
print(cost_red_NippyoB)

## [1] 0.9974925

8.4 限界費用削減の表を作成

Table_7dt3 <- tibble(
  col1 = c("Cost reduction"),
  col2 = (1 - cost_red_NippyoA) * 100,
  col3 = (1 - cost_red_NippyoB) * 100
)

colnames(Table_7dt3) <- c("Measure", "Nippyo and Brand A",
                          "Nippyo and Brand B")

kable(Table_7dt3,  align=rep('c', 3), booktabs = T, 
      caption = "Cost Reduction (Percent)")

Cost Reduction (Percent)
Measure	Nippyo and Brand A	Nippyo and Brand B
Cost reduction	0.4817918	0.2507483

8.5 限界費用削減をした場合の利潤と収入の変化

（表中の企業名は日評自動車としているホンダ、ブランドA社としている日産、ブランドB社としているスバル、ブランドC社としているトヨタ以外はそのままの社名を使っている。）

# 日評自動車とブランドA社が合併した場合の利益の導出
pro_rev_NippyoA_rc <- f_profit(data_2016$Maker, data_2016$p_NippyoA_TSfix, data_2016$mc_NippyoA_TSfix,
                               data_2016$share_NippyoA_TSfix, HH_2016) 

# 日評自動車とブランドB社が合併した場合の利益の導出
pro_rev_NippyoB_rc <- f_profit(data_2016$Maker, data_2016$p_NippyoB_TSfix, data_2016$mc_NippyoB_TSfix,
                               data_2016$share_NippyoB_TSfix, HH_2016) 

Table_7dt4 <- gen_pro_rev_table(pro_rev_2016, pro_rev_NippyoA_rc, pro_rev_NippyoB_rc)



# 表を作る
kable(Table_7dt4,  align=rep('c', 5), booktabs = T, caption = "Change in Profits and Revenues") %>%
  add_header_above(c(" ", "Nippyo and Brand A" = 2, "Nippyo and Brand B" = 2)) %>%
  kable_styling(full_width = F)

Change in Profits and Revenues
	Nippyo and Brand A		Nippyo and Brand B
Maker	Profits	Revenues	Profits	Revenues
Audi	5.9526234	18.979679	3.0075925	9.1100361
BMW	19.8982433	60.768864	9.6965900	27.9428483
Brand_A	2145.9781319	-5997.986468	31.7352231	74.1663144
Brand_B	29.0182365	77.222281	463.4234740	-3810.6605630
Brand_C	299.4337152	711.792297	135.5356403	302.2609392
Daihatsu	57.2135126	107.420965	23.5440181	42.7555185
Fiat	0.8028854	2.106037	0.3549816	0.8847146
Lexas	29.9674014	95.229809	15.5931903	47.5925500
Matsuda	31.3146894	77.120020	13.7284848	32.3561249
Mercedes	26.1588416	82.711499	13.5240619	41.2041190
Mitsubishi	10.7905482	25.592893	4.6895681	10.7156568
Nippyo	3254.6842912	-3256.890357	1707.2075360	-525.3996557
Suzuki	59.9106183	111.208310	24.6300595	44.2993308
Volkswagen	10.2906841	29.581151	4.7534686	12.9134976
Volvo	4.9204873	15.629353	2.4552115	7.3605475
Total	5986.3349098	-7839.513667	2453.8791003	-3682.4980210

8.6 限界費用削減をした場合の全体の利潤と収入の変化率

print((sum(pro_rev_NippyoA_rc$profit) - sum(pro_rev_2016$profit))/sum(pro_rev_2016$profit)*100)

## [1] 0.2122612

print((sum(pro_rev_NippyoA_rc$revenue) - sum(pro_rev_2016$revenue))/sum(pro_rev_2016$revenue)*100)

## [1] -0.1050482

print((sum(pro_rev_NippyoB_rc$profit) - sum(pro_rev_2016$profit))/sum(pro_rev_2016$profit)*100)

## [1] 0.0870087

print((sum(pro_rev_NippyoB_rc$revenue) - sum(pro_rev_2016$revenue))/sum(pro_rev_2016$revenue)*100)

## [1] -0.04934489

8.7 計算時間レポート

Appendix 1にかかった合計計算時間をレポートする。

tictoc::toc()

## 638.85 sec elapsed

経済セミナー連載「実証ビジネス・エコノミクス」
第４回「合併の効果は需要次第：消費者需要モデルの推定【応用編】」

上武康亮・遠山祐太・若森直樹・渡辺安虎

最終更新: 2021-09-14

1 はじめに

1.1 留意事項

1.2 謝辞

2 Rに関する下準備

3 データの準備

4 ランダム係数ロジットモデルの推定

4.1 定式化

4.2 データの下準備

4.3 各種関数の定義

4.3.1 市場シェアの計算コード

4.3.2 縮小写像によるBerryインバージョンのコード

4.3.3 GMM 目的関数を定義

4.3.4 標準誤差の計算

4.3.5 価格弾力性を計算する関数

4.4 パラメタ推定及び標準誤差の計算

5 限界費用の推定

6 合併シミュレーション

6.1 準備

6.2 所有構造行列を作成する

6.3 BLP推定同様の下準備

6.4 反実仮想シミュレーションにおける均衡計算コード

6.5 合併シミュレーション

7 シミュレーション結果のまとめ

7.1 合併シミュレーションによる価格・販売台数変化

7.2 日評自動車とブランドAの合併後に価格が変化しないような限界費用を計算する

7.3 日評自動車とブランドBの合併後に価格が変化しないような限界費用を計算する

7.4 合併シミュレーションの厚生分析

7.5 シミュレーションに関する計算時間

8 Appendix 1: 追加的なシミュレーション

8.1 関数の定義：総余剰が変化しないような限界費用削減

8.2 日評自動車とブランドAについて限界費用削減を計算

8.3 日評自動車とブランドBについて限界費用削減を計算

8.4 限界費用削減の表を作成

8.5 限界費用削減をした場合の利潤と収入の変化

8.6 限界費用削減をした場合の全体の利潤と収入の変化率

8.7 計算時間レポート

経済セミナー 連載「実証ビジネス・エコノミクス」 第４回 「合併の効果は需要次第：消費者需要モデルの推定【応用編】」

上武康亮・遠山祐太・若森直樹・渡辺安虎

最終更新: 2021-09-14

1 はじめに

1.1 留意事項

1.2 謝辞

2 Rに関する下準備

3 データの準備

4 ランダム係数ロジットモデルの推定

4.1 定式化

4.2 データの下準備

4.3 各種関数の定義

4.3.1 市場シェアの計算コード

4.3.2 縮小写像によるBerryインバージョンのコード

4.3.3 GMM 目的関数を定義

4.3.4 標準誤差の計算

4.3.5 価格弾力性を計算する関数

4.4 パラメタ推定及び標準誤差の計算

5 限界費用の推定

6 合併シミュレーション

6.1 準備

6.2 所有構造行列を作成する

6.3 BLP推定同様の下準備

6.4 反実仮想シミュレーションにおける均衡計算コード

6.5 合併シミュレーション

7 シミュレーション結果のまとめ

7.1 合併シミュレーションによる価格・販売台数変化

7.2 日評自動車とブランドAの合併後に価格が変化しないような限界費用を計算する

7.3 日評自動車とブランドBの合併後に価格が変化しないような限界費用を計算する

7.4 合併シミュレーションの厚生分析

7.5 シミュレーションに関する計算時間

8 Appendix 1: 追加的なシミュレーション

8.1 関数の定義：総余剰が変化しないような限界費用削減

8.2 日評自動車とブランドAについて限界費用削減を計算

8.3 日評自動車とブランドBについて限界費用削減を計算

8.4 限界費用削減の表を作成

8.5 限界費用削減をした場合の利潤と収入の変化

8.6 限界費用削減をした場合の全体の利潤と収入の変化率

8.7 計算時間レポート

経済セミナー連載「実証ビジネス・エコノミクス」
第４回「合併の効果は需要次第：消費者需要モデルの推定【応用編】」