1 はじめに

2021年8月・9月号の「実証ビジネス・エコノミクス」の第3回「プライシングの真髄は代替性にあり：消費者需要モデルの推定［基礎編2］」に関するRプログラムの紹介である。本プログラムによってざ誌面で紹介した分析を再現しつつ、実際に手を動かしながら学ぶことを目的とする。

1.1 留意事項

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

1.2 全体の流れ

下準備：パッケージの導入
データの読み込み、クリーニング
記述統計の作成
ロジットモデルの推定
推定結果に基づくエクササイズ

1.3 謝辞

前回に引き続いて藤田光明さん(サイバーエージェント)・哥丸連太朗さん(東京大学大学院経済学研究科)、また今回のプログラムのチェックに島本幸典さん・牧野圭吾さん・山田雅広さん(東京大学大学院経済学研究科)にご尽力頂きました。

2 Rに関する下準備

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

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

## Loading required package: tidyverse

## -- Attaching packages --------------------------------------------------- tidyverse 1.3.0 --

## √ ggplot2 3.3.2     √ purrr   0.3.4
## √ tibble  3.0.3     √ dplyr   1.0.2
## √ tidyr   1.1.2     √ stringr 1.4.0
## √ readr   1.4.0     √ forcats 0.5.0

## -- Conflicts ------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

require(fixest)

## Loading required package: fixest

## Warning: package 'fixest' was built under R version 4.0.5

## fixest 0.9.0, BREAKING changes! (Permanently remove this message with fixest_startup_msg(FALSE).) 
## - In i():
##     + the first two arguments have been swapped! Now it's i(factor_var, continuous_var) for interactions. 
##     + argument 'drop' has been removed (put everything in 'ref' now).
## - In feglm(): 
##     + the default family becomes 'gaussian' to be in line with glm(). Hence, for Poisson estimations, please use fepois() instead.

require(summarytools)

## Loading required package: summarytools

## Warning: package 'summarytools' was built under R version 4.0.5

## Registered S3 method overwritten by 'pryr':
##   method      from
##   print.bytes Rcpp

## 
## Attaching package: 'summarytools'

## The following object is masked from 'package:fixest':
## 
##     ctable

## The following object is masked from 'package:tibble':
## 
##     view

3 データの準備

一気に前回の作業を読み込む。

#まず自動車データを読み込みます。
data <- readr::read_csv(file = "source/CleanData_20180222.csv", locale = locale(encoding = "shift-jis"))

## 
## -- Column specification --------------------------------------------------------------------
## cols(
##   .default = col_double(),
##   Maker = col_character(),
##   Type = col_character(),
##   Name = col_character(),
##   comment = col_character(),
##   Model = col_character(),
##   kudou = col_character(),
##   kata = col_character(),
##   TransMission = col_character(),
##   base_color = col_character(),
##   option_color = col_character(),
##   種類 = col_character(),
##   engine = col_character(),
##   過給器 = col_character(),
##   FuelType = col_character()
## )
## i Use `spec()` for the full column specifications.

#続いて、家計数データ（潜在的な市場規模）を読み込む。ソースは以下。 <https://www.airia.or.jp/publish/file/r5c6pv0000006s9v-att/r5c6pv0000006saa.pdf>
dataHH <-
  readr::read_csv(file = "source/HHsize.csv")

## 
## -- Column specification --------------------------------------------------------------------
## cols(
##   year = col_double(),
##   HH = col_number()
## )

#価格を実質化するべく、消費者物価指数を読み込む。
#ソースは以下：<https://www.e-stat.go.jp/stat-search/file-download?statInfId=000031431770&fileKind=1>

dataCPI <-
  readr::read_csv(file = "source/zni2015s.csv", locale = locale(encoding = "shift-jis"))

## 
## -- Column specification --------------------------------------------------------------------
## cols(
##   .default = col_character()
## )
## i Use `spec()` for the full column specifications.

dataCPI <- dataCPI[6:56,　]
dataCPI <- 
  dataCPI %>% 
  rename( year = '類・品目',
          cpi = '総合') %>% 
  select(year, cpi) %>% 
  mutate( year = as.numeric(year), 
          cpi = as.numeric(cpi))

## データクリーニング
#まずは、自動車データを整理しつつ、家計データをマージする。

# 必要な変数のみをキープ。
data <-
  data %>%
  dplyr::select(
    Maker, Type, Name, Year, Sales, Model, price, kata,
    weight, FuelEfficiency,　HorsePower,
    overall_length, overall_width, overall_height
  ) %>%
  dplyr::rename(year = Year)

# 家計サイズをマージする。
data <-
  data %>%
  dplyr::left_join(dataHH)

## Joining, by = "year"

# CPIをマージする
data <- 
  data %>% 
  dplyr::left_join(dataCPI)

## Joining, by = "year"

# 燃費が欠損しているデータがある。今回は簡便な処理として観測から落とす。
data <-
  data %>%
  dplyr::filter(is.na(FuelEfficiency) == 0)

# 価格の実質化を行う。ここでは、2016年を基準年とする。
# また、価格の単位を100万円にする。元のデータは1万円
cpi2016 <- 
  dataCPI %>% 
  filter(year == 2016) %>% 
  select(cpi) %>% 
  as.double()

data <- 
  data %>% 
  mutate( price = price / (cpi / cpi2016) ) %>% 
  mutate( price = price / 100) %>%
  select(-cpi) 

# サイズ(高さ＊幅＊長さ)、燃費の重量に対する比率を定義する。
data <-
  data %>%
  dplyr::mutate(size = (overall_length / 1000) * (overall_width / 1000) * (overall_height / 1000)) %>%
  dplyr::mutate(hppw = HorsePower / weight) %>%
  dplyr::select(-HorsePower, -weight) %>% 
  dplyr::select(-starts_with("overall")) 
  
# 自動車の車種IDを作成する。
data <-
  data %>%
  dplyr::group_by(Name) %>%
  dplyr::mutate(NameID = dplyr::cur_group_id()) %>%
  dplyr::ungroup() %>%
  dplyr::relocate(NameID, year)

# マーケットシェアとOutside option shareを定義する。
data <-
  data %>%
  dplyr::group_by(year) %>%
  dplyr::mutate(inside_total = sum(Sales)) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(outside_total = HH - inside_total) %>%
  dplyr::mutate(share = Sales / HH) %>%
  dplyr::mutate(share0 = outside_total / HH) %>%
  dplyr::select(-inside_total, -outside_total)

# 操作変数の構築

# まず、マーケット・企業レベルにおける、各製品属性の和と自乗和を計算する。
# ここでacross関数は、最初に文字列ベクトルで指定した変数について、後ろにリスト内で定義した操作を適用している。
data <-
  data %>%
  dplyr::group_by(year, Maker) %>%
  dplyr::mutate(
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sum_own = ~ sum(.x, na.rm = TRUE) )   ),
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sqr_sum_own = ~ sum(.x^2, na.rm = TRUE) ) ),
    group_n = n()
  ) %>%
  dplyr::ungroup()

# 次に、マーケットレベルでの、各製品属性の和を計算する。
data <- 
  data %>% 
  dplyr::group_by(year) %>%
  dplyr::mutate( 
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sum_mkt = ~ sum(.x, na.rm = TRUE) )  ),
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sqr_sum_mkt = ~ sum(.x^2, na.rm = TRUE) )    ),
    mkt_n = n()
  ) %>%
  dplyr::ungroup() 

# 以上で定義した変数を利用して、まずBLP操作変数を構築する。
data <- 
  data %>% 
  dplyr::mutate(
    iv_BLP_own_hppw = hppw_sum_own - hppw,
    iv_BLP_own_FuelEfficiency = FuelEfficiency_sum_own - FuelEfficiency,
    iv_BLP_own_size = size_sum_own - size,
    iv_BLP_other_hppw = hppw_sum_mkt - hppw_sum_own,
    iv_BLP_other_FuelEfficiency = FuelEfficiency_sum_mkt - FuelEfficiency_sum_own,
    iv_BLP_other_size = size_sum_mkt - size_sum_own, 
    iv_BLP_own_num = group_n - 1, 
    iv_BLP_other_num = mkt_n - group_n) 

# 続いて、Differentiation IVを構築する。
data <- 
  data %>% 
  mutate(
    iv_GH_own_hppw = (group_n - 1) * hppw^2 + (hppw_sqr_sum_own - hppw^2) - 2 * hppw * (hppw_sum_own - hppw),
    iv_GH_own_FuelEfficiency = (group_n - 1) * FuelEfficiency^2 + (FuelEfficiency_sqr_sum_own - FuelEfficiency^2) - 2 * FuelEfficiency * (FuelEfficiency_sum_own - FuelEfficiency),
    iv_GH_own_size = (group_n - 1) * size^2 + (size_sqr_sum_own - size^2) - 2 * size * (size_sum_own - size),
    iv_GH_other_hppw = (mkt_n - group_n) * hppw^2 + (hppw_sqr_sum_mkt - hppw_sqr_sum_own) - 2 * hppw * (hppw_sum_mkt - hppw_sum_own),
    iv_GH_other_FuelEfficiency = (mkt_n - group_n) * FuelEfficiency^2 + (FuelEfficiency_sqr_sum_mkt - FuelEfficiency_sqr_sum_own) - 2 * FuelEfficiency * (FuelEfficiency_sum_mkt - FuelEfficiency_sum_own),
    iv_GH_other_size = (mkt_n - group_n) * size^2 + (size_sqr_sum_mkt - size_sqr_sum_own) - 2 * size * (size_sum_mkt - size_sum_own),
  ) %>%
  dplyr::select(
    -starts_with("sum_own"),
    -starts_with("sum_mkt"),
    -starts_with("sqr_sum_own"),
    -starts_with("sqr_sum_mkt"),
    -mkt_n,
    -group_n
  )

4 Logitモデルにおける価格弾力性

4.1 イントロダクションの日評自動車

まず前回同様、日評自動車のデータセットを用意する。

# 車種のID
IDvec <- sort(unique(data$NameID))
J <- length(IDvec)

# 乱数のシードを固定する。これを行うことで、毎回同じ乱数を得ることができる。
set.seed(125)

# sample関数を使って、先に用意した車種IDベクトルから、重複無しで30車種を取得する。
# なお、sample関数は内部で乱数を利用しているため、先の乱数シード固定が重要となる。
NIPPYOautoIDvec <- sample(IDvec, size = 30, replace = FALSE)
NIPPYOautoIDvec <- sort(NIPPYOautoIDvec)

# 日評自動車のデータセットを作成する。
data_NIPPYO <-
  data %>% 
  filter( NameID %in% NIPPYOautoIDvec ) %>% 
  select( year, share, NameID, Sales, price, hppw, FuelEfficiency, size, Name ) %>% 
  mutate( log_sales = log(Sales),
          log_price = log(price) )

4.2 価格弾力性行列を作成する

data <-
  data %>%
  dplyr::mutate(logit_share = log(share) - log(share0))

# Differentiation IVを用いた結果
iv_GH <-
  fixest::feols(
    logit_share ~ hppw + FuelEfficiency + size | 0  |
      price ~ iv_GH_own_hppw + iv_GH_own_FuelEfficiency + iv_GH_own_size +
      iv_GH_other_hppw + iv_GH_other_FuelEfficiency + iv_GH_other_size,
    data = data
  )

data_NIPPYO %>% 
  filter(year == 2016) %>% 
  select(price, share, NameID) %>% 
  arrange(NameID) -> dt2016

price <- dt2016$price
share <- dt2016$share
NameID <- dt2016$NameID

own_elas <- iv_GH$coefficients["fit_price"]*price*(1-share)

cross_elas <- (-1)*iv_GH$coefficients["fit_price"]*price*share
J <- length(own_elas)

elas_mat <- matrix( rep(cross_elas, J), nrow = J, ncol = J )
diag(elas_mat) <- own_elas

表１の価格弾力性行列。以下の製品について見る。

87: ベータード
122: セダン(A)
155: SUV(B)
178: 軽自動車(C)

elas_mat[ NameID %in% c(87, 122, 155, 178),  NameID %in% c(87, 122, 155, 178)  ]

##               [,1]          [,2]         [,3]          [,4]
## [1,] -1.7645536904  0.0011492884  0.001149288  0.0011492884
## [2,]  0.0012204172 -0.8186885645  0.001220417  0.0012204172
## [3,]  0.0001481753  0.0001481753 -0.959449003  0.0001481753
## [4,]  0.0018450942  0.0018450942  0.001845094 -0.6717501635

5 入れ子型ロジット推定のためのBLP操作変数の作成

# まず、マーケット・企業レベルにおける、各製品属性の和と自乗和を計算する。
# ここでacross関数は、最初に文字列ベクトルで指定した変数について、後ろにリスト内で定義した操作を適用している。
data <-
  data %>%
  dplyr::group_by(year, Maker,Type) %>%
  dplyr::mutate(
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sum_own = ~ sum(.x, na.rm = TRUE) )   ),
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sqr_sum_own = ~ sum(.x^2, na.rm = TRUE) ) ),
    group_n = n()
  ) %>%
  dplyr::ungroup()

# 次に、マーケットレベルでの、各製品属性の和を計算する。
data <- 
  data %>% 
  dplyr::group_by(year, Type) %>%
  dplyr::mutate( 
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sum_mkt = ~ sum(.x, na.rm = TRUE) )  ),
    dplyr::across( c("hppw", "FuelEfficiency", "size"),
                   list( sqr_sum_mkt = ~ sum(.x^2, na.rm = TRUE) )    ),
    mkt_n = n()
  ) %>%
  dplyr::ungroup() 

# 以上で定義した変数を利用して、まずBLP操作変数を構築する。
data <- 
  data %>% 
  dplyr::mutate(
    iv_BLP_own_hppw_nest = hppw_sum_own - hppw,
    iv_BLP_own_FuelEfficiency_nest = FuelEfficiency_sum_own - FuelEfficiency,
    iv_BLP_own_size_nest = size_sum_own - size,
    iv_BLP_other_hppw_nest = hppw_sum_mkt - hppw_sum_own,
    iv_BLP_other_FuelEfficiency_nest = FuelEfficiency_sum_mkt - FuelEfficiency_sum_own,
    iv_BLP_other_size_nest = size_sum_mkt - size_sum_own, 
    iv_BLP_own_num_nest = group_n - 1, 
    iv_BLP_other_num_nest = mkt_n - group_n) 

# 続いて、Differentiation IVを構築する。
data <- 
  data %>% 
  mutate(
    iv_GH_own_hppw_nest = (group_n - 1) * hppw^2 + (hppw_sqr_sum_own - hppw^2) - 2 * hppw * (hppw_sum_own - hppw),
    iv_GH_own_FuelEfficiency_nest = (group_n - 1) * FuelEfficiency^2 + (FuelEfficiency_sqr_sum_own - FuelEfficiency^2) - 2 * FuelEfficiency * (FuelEfficiency_sum_own - FuelEfficiency),
    iv_GH_own_size_nest = (group_n - 1) * size^2 + (size_sqr_sum_own - size^2) - 2 * size * (size_sum_own - size),
    iv_GH_other_hppw_nest = (mkt_n - group_n) * hppw^2 + (hppw_sqr_sum_mkt - hppw_sqr_sum_own) - 2 * hppw * (hppw_sum_mkt - hppw_sum_own),
    iv_GH_other_FuelEfficiency_nest = (mkt_n - group_n) * FuelEfficiency^2 + (FuelEfficiency_sqr_sum_mkt - FuelEfficiency_sqr_sum_own) - 2 * FuelEfficiency * (FuelEfficiency_sum_mkt - FuelEfficiency_sum_own),
    iv_GH_other_size_nest = (mkt_n - group_n) * size^2 + (size_sqr_sum_mkt - size_sqr_sum_own) - 2 * size * (size_sum_mkt - size_sum_own),
  ) %>%
  dplyr::select(
    -starts_with("sum_own"),
    -starts_with("sum_mkt"),
    -starts_with("sqr_sum_own"),
    -starts_with("sqr_sum_mkt"),
    -mkt_n,
    -group_n
  )

6 入れ子型ロジットモデルの推定とその応用

本誌で導入した、以下のモデルを推定しよう。

\[ \log S_{jt} - \log S_{0t} = -\alpha p_{jt} + \beta'\mathbf{X}_{jt} + \xi_{jt} \]

推定に際しては、OLS、BLP操作変数による結果、そしてDifferentiation IVによる結果の3通りについて比較する。

# まず被説明変数を定義する。
data <-
  data %>%
  dplyr::mutate(logit_share = log(share) - log(share0))


# インサイドシェアを定義する。

data <- 
  data %>% 
  dplyr::group_by(year, Type) %>% 
  dplyr::mutate( sum_year_body = sum(Sales)) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate( inside_share = Sales / sum_year_body, 
    log_inside_share = log(Sales / sum_year_body) )


# OLSの結果
ols <-
  fixest::feols(logit_share ~ price + log_inside_share + hppw + FuelEfficiency + size | 0 , data = data)

# BLP操作変数を用いた結果
iv_BLP2 <-
  fixest::feols(
    logit_share ~ hppw + FuelEfficiency + size | 0 |
      price + log_inside_share ~ iv_BLP_own_hppw_nest + iv_BLP_own_FuelEfficiency_nest + iv_BLP_own_size_nest +
      iv_BLP_other_hppw_nest + iv_BLP_other_FuelEfficiency_nest + iv_BLP_other_size_nest + iv_BLP_own_num_nest +  iv_BLP_other_num_nest,
    data = data
  )



fixest::etable( list(ols, iv_BLP2),  
                se = "hetero",
                fitstat = c("r2", "n", "ivf" ) , 
                signifCode = NA, 
                dict = c(price = "自動車価格",
                         hppw = "馬力／重量",
                         FuelEfficiency = "燃費(キロメートル／ 1 リットル)",
                         size = "サイズ",
                         `(Intercept)` = "定数項"),
                digits = 2,
                digits.stats = 2,
                depvar = FALSE)

6.1 入れ子型ロジットモデルにおける自己価格弾力性の計算

各推定結果について、自己価格弾力性を計算しよう。

alpha1 <- ols$coefficients["price"]
sigma1 <- ols$coefficients["log_inside_share"]

alpha2 <- iv_BLP2$coefficients["fit_price"]
sigma2 <- iv_BLP2$coefficients["fit_log_inside_share"]


data %>% 
  dplyr::mutate( 
    own_elas_ols   = alpha1*price*(1-sigma1*inside_share - (1-sigma1)*share)/(1-sigma1),
    own_elas_ivblp = alpha2*price*(1-sigma2*inside_share - (1-sigma2)*share)/(1-sigma2)  ) -> data_elas

data_elas %>% 
  dplyr::select( starts_with("own_elas") ) %>% 
  summarytools::descr( transpose = TRUE, 
                    stats = c("mean", "sd", "med", "min", "max") )

価格弾力性行列を作る。

data_NIPPYO <-
  data %>% 
  filter( NameID %in% NIPPYOautoIDvec ) %>% 
  select( year, share,  Type, inside_share, NameID, Sales, price, hppw, FuelEfficiency, size ) %>% 
  mutate( log_sales = log(Sales),
          log_price = log(price) ) 

data_NIPPYO %>% 
  filter(year == 2016) %>% 
  select(price, Type, share,inside_share, NameID) %>% 
  arrange(NameID) -> dt2016

price <- dt2016$price
share <- dt2016$share
NameID <- dt2016$NameID
inside_share <- dt2016$inside_share
group <- dt2016$Type

# Own elasticity
own_elas <- alpha2*price*(1-sigma2*inside_share - (1-sigma2)*share)/(1-sigma2)

# Cross elasticity outside the group
cross_elas_othergroup <- (-1)*alpha2*price*share
J <- length(own_elas)
cross_elas_othergroup <- matrix( rep(cross_elas_othergroup, J), nrow = J, ncol = J )

# Cross elasticity within the group (matrix representation)

price_l_mat <- matrix(rep(price, J), nrow = J, ncol = J)
share_l_mat <- matrix(rep(share, J), nrow = J, ncol = J)
insideshare_l_mat <- matrix(rep(inside_share, J), nrow = J, ncol = J)
cross_elas_samegroup <- (-1)*alpha2*price_l_mat*(sigma2*insideshare_l_mat + (1-sigma2)* share_l_mat )/(1-sigma2)
                                             

# Indicator of group
temp_mat1 <- matrix(rep(group, J), nrow = J, ncol = J)
temp_mat2 <- t(temp_mat1)
ind_same_group <- (temp_mat1 == temp_mat2)
ind_other_group <- (temp_mat1 != temp_mat2)

# Construct Elasticity Matrix
elas_mat_nl <- cross_elas_samegroup*ind_same_group + cross_elas_othergroup*ind_other_group
diag(elas_mat_nl) <- own_elas


elas_mat_nl[ NameID %in% c(87, 122, 155, 178),  NameID %in% c(87, 122, 155, 178)  ]

##              [,1]         [,2]         [,3]         [,4]
## [1,] -5.119695908  0.047757294  0.047757294  0.001361722
## [2,]  0.050712966 -2.348807986  0.050712966  0.001445999
## [3,]  0.006157248  0.006157248 -2.802170978  0.000175564
## [4,]  0.002186141  0.002186141  0.002186141 -1.832963467

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

7.1 定式化

以下の定式化を考える。

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

ここで、\(\beta_i^{size}, \beta_i^0, \alpha_i\)は正規分布に従うランダム係数である。

7.2 推定のOverview

Step 1はデータに関する下準備である。

Step 2から4では、GMM推定のために必要な関数を順番に用意していく。

市場シェアを計算する関数 (Step 2)
縮小写像によるBerryインバージョンを行う関数 (Step 3)
GMMの目的関数 (Step 4)

市場シェアを計算する関数はStep 3・4両方で使用する。また、Berryインバージョンを行う関数はGMM目的関数を評価する際に必要となる。 GMM目的関数を定義した後、数値最適化でGMM推定を行う。

Step 5では推定したパラメタの標準誤差の計算を行う。最後に推定結果をテーブルでまとめる。

7.3 Step 1:下準備

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

# まずデータをソートする。マーケット順かつモデル順

data %>% 
    arrange( year, NameID) -> data

# マーケットとモデルの情報をGET
data %>% 
    select(year, NameID) -> Info

N <- length(Info$year)
T <- length(unique(Info$year))

# 平均効用に入ってくる部分。内生性がある価格を含んでいる。
X1 <-
  data %>%
  mutate(cons = 1) %>%
  select(cons, price, FuelEfficiency, hppw, size) %>%
  as.matrix()

# ランダム係数とInteractする部分。
X2 <-
  data %>%
  mutate(cons = 1) %>%
  select(price, cons, size) %>%
  as.matrix()


# 操作変数行列。ここは外生変数と追加的な操作変数を含む
Z <-
  data %>%
  mutate(cons = 1) %>%
  select(
    cons, FuelEfficiency, hppw, size,
    starts_with("iv_GH"),
    -ends_with("nest")
  ) %>%
  as.matrix()

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

datalist = list()
datalist$X1 <- X1
datalist$X2 <- X2
datalist$Z <- Z
datalist$ShareVec <- ShareVec
datalist$marketindex <- data$year
datalist$logitshare <- data$logit_share

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

set.seed(111)

Nsim = 500

draw_vec = rnorm(Nsim*ncol(X2))

datalist$draw_vec <- draw_vec

parameter = list()
parameter$Nsim = Nsim
parameter$T = T 
parameter$N = N
parameter$theta2 = c(0.001, 0.001)

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

7.4 Step 2: 市場シェアの計算コード

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

f_mktshare <- function(datalist, parameter, delta){
    
    # Unpack data
    X2 <- datalist$X2
    ns <- parameter$Nsim
    T <- parameter$T
    N <- parameter$N
    
    draw_vec <- datalist$draw_vec
    
    # 20210331 theta2とK2をunpackする位置がずれて回らなかったので修正しました．
    theta2 <- parameter$theta2
    K2 <- length(theta2)
    
    # 20210331 X2が複数の場合，nsのみが認識するので括弧(ns*K2)を付けました．
    draw_vec <- datalist$draw_vec[1:(ns * K2)]
    
    marketindex <- datalist$marketindex
    
    # Nonlinear component. mu (N, ns) matrix.
    mu <- X2 %*% diag(x = theta2, nrow = length(theta2)) %*% matrix(draw_vec, nrow = K2)
    
    
    denom_outside <-  exp(matrix(0, T, ns) )
    
    #delta = matrix(0, N,1)
    
    delta_mu <-  delta %*% matrix(1, nrow = 1, ncol = ns ) + mu
    exp_delta_mu <- exp(delta_mu)
    
    # ロジット確率の分母を計算する。これは、市場ごとに和を計算する。
    tempmat <- datalist$tempmat
    
    denom_temp <- t( t(exp_delta_mu) %*% tempmat )
    denom_temp <- denom_temp + denom_outside
    
    # denom is (N * ns) matrix that captures the denominator of the logit prob.
    denom <- tempmat %*% 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)
    
}

7.5 Step 3: 縮小写像によるBerryインバージョンのコード

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

f_contraction <- function(datalist, parameter, delta_ini){
    
    
    tol <- 1e-11
    share_obs <- datalist$ShareVec
    norm <- 1e+10
    
    delta_old <- delta_ini
    exp_delta_old <- exp(delta_old)
    
    iter <- 0
    
    while (norm > tol && iter < 1000){
        
        pred_mkt_share <- f_mktshare(datalist, parameter, delta_old)
        
        exp_delta <- exp_delta_old * share_obs / pred_mkt_share
        
        t <- abs(exp_delta - exp_delta_old)
        norm <- max(t)
        
        #print(norm)
        
        # Update
        #exp_delta_old <- 0.95*exp_delta + 0.05*exp_delta_old
        exp_delta_old <- exp_delta 
        
        delta_old <- log(exp_delta)
        iter <- iter + 1
        
    }
    
    return( log(exp_delta) )
    
}

7.6 Step 4: GMM 目的関数を定義し、数値最適化する

ここでは、非線形パラメタtheta2に対してGMM目的関数を返す関数f_GMMobjを定義する。以下、関数の留意点をいくつか。

以下では変数delta_globalをグローバル変数として扱っており、関数f_GMMobjを評価するたびに、delta_globalもアップデートするようになっている。この変数は、f_GMMobjを評価する際に必要となる縮小写像アルゴリズムの初期値である。最適化が進むにつれてf_GMMobjが評価するパラメタが徐々に変化するが、その評価の際に必要な縮小写像アルゴリズムの初期値もなるべく近いパラメタの元で得られた平均効用を用いる方が、計算速度上望ましい。なおここでは、永続代入演算子delta_global <<- deltaを利用しているが、その詳細については後述する。
GMM目的関数のインプットは非線形パラメタtheta2のみである。これは、平均効用に線形の形で入ってくるパラメタ\(\theta_1\)は、非線形パラメタを所与とすると線形GMMで解析的に求まるためである。この点については本連載ウェブ付録・第3回のページにアップされた補足資料を参照されたい。

f_GMMobj <- function(theta2, parameter, datalist, option) {

  # print(theta2)
  
  # 引数option
  # 0:最適化の際に指定
  # 1:推定値やmean utility(delta)を得る際に指定

  param <- parameter
  param$theta2 <- theta2

  # Contractionの初期値として、グローバル変数のdelta_globalを用いる。
  delta_ini <- delta_global

  # Contraction mapping
  delta <- f_contraction(datalist, param, delta_ini)
  
  # 平均効用を永続代入演算子でアップデートしている。 
  delta_global <<- delta

  # Obj function
  X1 <- datalist$X1

  # IV matrix
  Z <- datalist$Z

  # Z <- Z[, 1:8]/1000

  # 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
  # beta_hat2 = solve( t(X1)%*%X1 )%*%t(X1)%*%delta

  beta_hat <- solve(t(X1) %*% Z %*% W %*% t(Z) %*% X1) %*% t(X1) %*% Z %*% W %*% t(Z) %*% delta

  Xi <- delta - X1 %*% beta_hat


  # 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, beta_hat = beta_hat, delta = delta))
  } 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)"
    )
  }
}

定義した関数を用いて、実際にGMM推定をしてみよう。

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

# グローバル変数としてアップデートする平均効用
delta_global <- data$logit_share

# 最適化
result <-
  optim(
    par = c(0.3, 18, 0.01),
    fn = f_GMMobj,
    method = "L-BFGS-B",
    lower = c(0,0,0 ),
    parameter = parameter,
    datalist = datalist,
    option = 0
  )

7.6.1 技術的な補足：永続代入演算子とグローバル変数

通常、関数内部における変数はローカル変数と呼ばれ、関数内で処理が完結し、関数の外（すなわちワークスペース）の変数には影響を与えない。一方、永続代入演算子<<-を用いることで、関数内部での処理によって、関数の外の変数を変更することが可能である。このような変数をグローバル変数と呼ぶ。通常、永続代入演算子やグローバル変数を用いることは、変数間の関係が煩雑になったり、予期しない変数のアップデートが起きうる点から望ましくない。今回は計算速度を早めるべく利用しているが、そのような場合にも永続代入演算子を利用していることを明示的することが、コードをクリアに書く上で重要である。

7.7 Step 5: 標準誤差の計算

ここでは標準誤差を計算する関数を定義し、その上で計算する。

f_SE <- function(theta2, parameter, datalist) {
  param <- parameter
  param$theta2 <- theta2

  delta_ini <- datalist$logitshare

  # Contraction mapping
  delta <- f_contraction(datalist, param, delta_ini)

  # 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
  beta_hat <- solve(t(X1) %*% Z %*% W %*% t(Z) %*% X1) %*% t(X1) %*% Z %*% W %*% t(Z) %*% delta

  # ResidualをGetする。
  Xi <- delta - X1 %*% beta_hat

  # 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
    tempparam$theta2[k] <- theta2[k] + 1e-6

    delta_add <- f_contraction(datalist, tempparam, delta_ini)

    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)
}

以上の関数を用いて標準誤差を計算する。

result2 <- f_GMMobj( result$par, parameter, datalist, option = 1)

delta <- result2$delta

se <- f_SE(theta2 = result$par, parameter = parameter, datalist = datalist)

7.8 ランダム係数ロジットの推定結果

beta_hat <- rbind(result2$beta_hat, result$par[1], result$par[2], result$par[3])
rownames(beta_hat)[6] <- "random_price"
rownames(beta_hat)[7] <- "random_constant"
rownames(beta_hat)[8] <- "random_size"
colnames(beta_hat) <- "coeff"
names(se)[6] <- "random_price"
names(se)[7] <- "random_constant"
names(se)[8] <- "random_size"

result_table <-
  as.data.frame(beta_hat) %>%
  mutate(se = se)

print(result_table)

##          coeff          se
## 1 -32.98386528 16.32257486
## 2  -0.88319656  0.34237767
## 3   0.10535101  0.01416849
## 4   7.92300169  3.54262053
## 5   0.24477852  0.14386844
## 6   0.30024626  0.18926904
## 7  18.00016113 16.02090600
## 8   0.01023621  0.83677330

8 価格弾力性行列の計算

まず価格弾力性を計算する関数を定義する。

f_elasticity <- function(datalist, parameter, beta_hat, delta) {

  # Unpack data
  X2 <- datalist$X2
  ns <- parameter$Nsim
  N <- parameter$N
  price <- datalist$X1[, colnames(datalist$X1) == "price"] %>% as.matrix()
  theta2 <- parameter$theta2
  K2 <- length(theta2)
  draw_vec <- datalist$draw_vec[1:(ns * K2)]
  uniquemarketindex <- datalist$uniquemarketindex

  # 製品を示す行列を取り出す
  # tempmat2 <- datalist$tempmat2


  # 以下マーケットシェアの関数とほぼ同じ
  # Nonlinear component. mu (N, ns) matrix.
  mu <- X2 %*% diag(x = theta2, nrow = length(theta2)) %*% matrix(draw_vec, nrow = K2)

  denom_outside <- exp(matrix(0, T, ns))
  delta_mu <- delta %*% matrix(1, nrow = 1, ncol = ns) + mu
  exp_delta_mu <- exp(delta_mu)
  tempmat <- datalist$tempmat
  denom_temp <- t(t(exp_delta_mu) %*% tempmat)
  denom_temp <- denom_temp + denom_outside

  denom <- tempmat %*% 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, ]
  alpha_i <- beta_hat[rownames(beta_hat) == "price"] + beta_hat[rownames(beta_hat) == "random_price"]*draw_for_price

  # 弾力性行列を保存する空のリストを作成する
  elaslist <- as.list(rep(NA, T))

  # 各マーケットごとに弾力性行列を作成する．
  for (t in 1:T) {

    # マーケットの表示
    #cat("calculate the elasticity for Market", uniquemarketindex[t], "\n")
    
    # Number of products 
    J_t <- sum(marketindex == 2005 + t)

    # それぞれのマーケットに注目した(J*ns)の選択確率の行列を作成する
    # あるマーケットに対して，製品ごとに集計した選択確率の行列(J*ns)を作成する
    ag_model_s_i <- s_jt_i[ marketindex == 2005 + 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 == 2005 + 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+2005)
  }
  return(elaslist)
}

早速、価格弾力性行列を示そう。

# 弾力性の表示の前処理
parameter$theta2 <- result$par

elasticity <- f_elasticity(datalist = datalist, parameter = parameter, beta_hat = beta_hat, delta = delta)

NameID2016 <- data %>% filter(year == 2016) %>% select(NameID) %>% pull()
elasticity$elasmat_2016[ NameID2016 %in% c(87, 122, 155, 178),  NameID2016 %in% c(87, 122, 155, 178)  ]

##             [,1]         [,2]        [,3]         [,4]
## [1,] -2.16720791  0.018628907  0.01863600  0.018583792
## [2,]  0.01978241 -1.196164320  0.02403604  0.025724955
## [3,]  0.00240276  0.002918292 -1.38792236  0.003003349
## [4,]  0.02983586  0.038892599  0.03739835 -0.991414389

9 応用：プライシング

9.1 需要曲線と収入曲線を書く

需要関数の推定結果にもとづいて、需要曲線を書こう。前回と手順は同様である。

まず、推定結果から\(\xi_{jt}\)をゲットする。
関数を作成する。
どこか市場を固定する。
製品を一つ固定する。
その製品の価格を変えたときに、どの製品のSalesがどう変わるかをPredictする。その際、他の製品価格はデータのものに固定しておく。

まず関数を定義する。

f_revenue <- function(price_cand, data , datalist , datalist_temp, result2, parameter, option ){
  
  # ベータードの限界費用。今回はラーナーの公式から簡便的に定義する。より詳しい内容は次回の連載第４回を参照。
  mc_betado <- 3.198*(1- 1/abs(-2.16720791))

  tempprice <- data$price
  tempprice[data$NameID == 87 & data$year == 2016] <- price_cand
  
  datalist_temp <- datalist
  datalist_temp$X1[,"price"] <- tempprice
  datalist_temp$X2[,"price"] <- tempprice
  
  # Re-calculate 平均効用 with new price
  org_xi <- result2$delta - datalist$X1 %*% result2$beta_hat
  new_delta <- datalist_temp$X1 %*% result2$beta_hat + org_xi
  
  mktshare <- f_mktshare(datalist = datalist_temp  , parameter = parameter, delta = new_delta)
  
  quant <- mktshare*data$HH
  revenue <- tempprice*quant
  
  # ベータードのみの収入
  revenuevec <-　revenue[data$NameID == 87 & data$year == 2016]

  # 日評自動車全体の収入
  revenuevec2 <-　sum(revenue[data$NameID %in% NIPPYOautoIDvec & data$year == 2016] )
  
  
  # ベータードのみの利潤
  pivec <- revenuevec-mc_betado*quant[data$NameID == 87 & data$year == 2016]
  
  # ベータードのみの利潤+日評自動車の他の車種の収入
  pivec2 <- revenuevec2 - mc_betado*quant[data$NameID == 87 & data$year == 2016]

  if (option == "own"){
    return(revenuevec)
  } else if (option == "total"){
    return(revenuevec2)
  } else if (option == "ownpi"){
    return(pivec)
  } else if (option == "totalpi"){
    return(pivec2)
  }
  
}

定義した関数を用いて、まず図を作成する。

NameID_target <- 87

# 価格の範囲として、0.5(50万円)から5(500万円)を考える。なお、もとの価格は3.198 (319.8万円)
pricevec <- seq(from = 1.73, to = 4,　by = 0.01)
pivec <- numeric(length(pricevec))
pivec2 <- numeric(length(pricevec))


for ( i in 1:length(pricevec)){

  # ベータードのみの利潤
  pivec[i] <- f_revenue(pricevec[i], data , datalist , datalist_temp, result2, parameter, option= "ownpi")
  
  # ベータードのみの利潤+日評自動車の他の車種の収入
  pivec2[i] <- f_revenue(pricevec[i], data , datalist , datalist_temp, result2, parameter, option= "totalpi")
  
}


# ベータードのみ
dt_fig1 <- tibble( price = pricevec*100, 
                  pi1 = pivec*100/10000) %>% 
  pivot_longer( -price, names_to = "name")

fig_pi1 <- ggplot(dt_fig1, aes(price, value) ) + 
  geom_point() 
plot(fig_pi1)

# ベータードのみの利潤+日評自動車の他の車種の収入
dt_fig2 <- tibble( price = pricevec*100, 
                  pi2 = pivec2*100/10000) %>% 
  pivot_longer( -price, names_to = "name")

fig_pi2 <- ggplot(dt_fig2, aes(price, value) ) + 
  geom_point() 
plot(fig_pi2)

最後に、数値最適化を使って、利潤を最大にする価格を求める。

まずベータードのみ

# ベータードのみ
optimprice3 <- optimise( f_revenue, interval = c(0.3, 5), maximum = TRUE, 
          data = data , datalist = datalist , result2 = result2, parameter = parameter,
          option = "ownpi")
print(optimprice3)

## $maximum
## [1] 3.197997
## 
## $objective
## [1] 54601.47

つづいてベータードのみの利潤+日評自動車の他の車種の収入

optimprice4 <- optimise( f_revenue, interval = c(0.3, 5), maximum = TRUE, 
          data = data , datalist = datalist , result2 = result2, parameter = parameter,
          option = "totalpi")
print(optimprice4)

## $maximum
## [1] 3.382552
## 
## $objective
## [1] 741237.1

ベータードのみの利潤を最大にするような価格における、「ベータードのみの利潤+日評自動車の他の車種の収入」はどうなるか？

f_revenue( price = optimprice3$maximum,  data = data , datalist = datalist , result2 = result2, parameter = parameter,    option = "totalpi")

## [1] 740893.7

経済セミナー連載「実証ビジネス・エコノミクス」
第３回「プライシングの真髄は代替性にあり：消費者需要モデルの推定［基礎編2］」

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

最終更新: 2021-07-15

1 はじめに

1.1 留意事項

1.2 全体の流れ

1.3 謝辞

2 Rに関する下準備

3 データの準備

4 Logitモデルにおける価格弾力性

4.1 イントロダクションの日評自動車

4.2 価格弾力性行列を作成する

5 入れ子型ロジット推定のためのBLP操作変数の作成

6 入れ子型ロジットモデルの推定とその応用

6.1 入れ子型ロジットモデルにおける自己価格弾力性の計算

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

7.1 定式化

7.2 推定のOverview

7.3 Step 1:下準備

7.4 Step 2: 市場シェアの計算コード

7.5 Step 3: 縮小写像によるBerryインバージョンのコード

7.6 Step 4: GMM 目的関数を定義し、数値最適化する

7.6.1 技術的な補足：永続代入演算子とグローバル変数

7.7 Step 5: 標準誤差の計算

7.8 ランダム係数ロジットの推定結果

8 価格弾力性行列の計算

9 応用：プライシング

9.1 需要曲線と収入曲線を書く

経済セミナー 連載「実証ビジネス・エコノミクス」 第３回 「プライシングの真髄は代替性にあり：消費者需要モデルの推定［基礎編2］」

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

最終更新: 2021-07-15

1 はじめに

1.1 留意事項

1.2 全体の流れ

1.3 謝辞

2 Rに関する下準備

3 データの準備

4 Logitモデルにおける価格弾力性

4.1 イントロダクションの日評自動車

4.2 価格弾力性行列を作成する

5 入れ子型ロジット推定のためのBLP操作変数の作成

6 入れ子型ロジットモデルの推定とその応用

6.1 入れ子型ロジットモデルにおける自己価格弾力性の計算

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

7.1 定式化

7.2 推定のOverview

7.3 Step 1:下準備

7.4 Step 2: 市場シェアの計算コード

7.5 Step 3: 縮小写像によるBerryインバージョンのコード

7.6 Step 4: GMM 目的関数を定義し、数値最適化する

7.6.1 技術的な補足：永続代入演算子とグローバル変数

7.7 Step 5: 標準誤差の計算

7.8 ランダム係数ロジットの推定結果

8 価格弾力性行列の計算

9 応用：プライシング

9.1 需要曲線と収入曲線を書く

経済セミナー連載「実証ビジネス・エコノミクス」
第３回「プライシングの真髄は代替性にあり：消費者需要モデルの推定［基礎編2］」