経済セミナー 連載「実証ビジネス・エコノミクス」
第２回 「需要を制すものはプライシングを制す：消費者需要モデルの推定(基礎編1)」

1 はじめに

2021年6月・7月号の「実証ビジネス・エコノミクス」の第2回「需要を制すものはプライシングを制す：消費者需要モデルの推定(基礎編1)」に関するRプログラムの紹介である。本プログラムによって紙面で紹介した分析を再現しつつ、実際に手を動かしながら学ぶことを目的とする。

1.1 留意事項

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

1.2 全体の流れ

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

1.3 謝辞

消費者需要の推定(連載第2・3・4回)において使用する日本の自動車市場のデータは藤田光明さん(現・サイバーエージェント)から提供頂きました。また、連載第２・3回の消費者需要モデルの推定プログラムの準備に際して哥丸連太朗さん(現・東京大学大学院経済学研究科)に手助け頂きました。この場を借りて御礼申し上げます。

2 Rに関する下準備

2.1 R及びRstudioのインストールに関して

Rの初心者向けのオンラインリソースとして、宋財泫・矢内勇生「私たちのR: ベストプラクティスの探究」を推奨する。。 インストールに関しては、同ウェブサイトにおける、「3. Rのインストール」を参照されたい。OSに応じたインストール方法が詳細に解説されている。

2.2 必要なパッケージの読み込み

以下のコードを回すのに必要なパッケージ一式を準備します。 もしインストールされていないものがありましたら、install.package() コマンドを利用しインストールして下さい。

簡単な紹介：

• tidyverse: データ加工やグラフ作成などのパッケージが一式入っています。特に、データ加工のdplyr、グラフ作成のggplot2、そしてパイプ演算のためのmagrittrを使います。
• fixest: 線形モデル推定のためのパッケージ。今回は利用しませんが、ダミー変数を多数入れた場合に計算が非常に早いパッケージになります。
• summarytools: 記述統計作成用のパッケージ。

## Loading required package: tidyverse

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

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

## √ ggplot2 3.3.3     √ purrr   0.3.4
## √ tibble  3.0.3     √ dplyr   1.0.5
## √ tidyr   1.1.2     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0

## Warning: package 'ggplot2' was built under R version 4.0.3

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

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

## Loading required package: fixest

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

## Loading required package: summarytools

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

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

## For best results, restart R session and update pander using devtools:: or remotes::install_github('rapporter/pander')

## 
## Attaching package: 'summarytools'

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

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

2.2.1 細かい留意点

Rで外部ライブラリを利用する際には、以下の二通りがあります。

1. require()やlibrary()で事前に読み込み、パッケージ内の関数を呼び出す。例：require(dplyr)の後にmutate().

2. パッケージ全体を事前に読み込みはせず、特定のパッケージの関数を使う際に[package_name]::[function_name]として呼び出す。例：dplyr::mutate().

今回は、どの関数がどのライブラリから来ているをわかりやすくするために、後者の方法でプログラミングしていきます。

3 データの準備

3.1 データ読み込み

まず自動車データを読み込みます。

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

## Parsed with 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()
## )

## See spec(...) for full column specifications.

続いて、家計数データ（潜在的な市場規模）を読み込む。ソースは以下。 https://www.airia.or.jp/publish/file/r5c6pv0000006s9v-att/r5c6pv0000006saa.pdf

dataHH <-
  readr::read_csv(file = "source/HHsize.csv")

## Parsed with 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"))

## Parsed with column specification:
## cols(
##   .default = col_character()
## )

## See spec(...) for full column specifications.

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

3.2 データクリーニング

まずは、自動車データを整理しつつ、家計データをマージする。

# 必要な変数のみをキープ。
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)

続いて、市場シェアとアウトサイドオプションシェアを定義する。ここでは、マーケット(年)ごとに総販売台数(変数inside_total)を定義し、それを用いてアウトサイドオプションのシェアを計算している。

# マーケットシェアと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)

3.3 操作変数の構築

操作変数を構築しよう。 ここでは、BLP型操作変数と、Gandhi-HoudeのDifferentiation IVを定義する。 コードが若干複雑なので、構築のアイデアを先に説明しよう。

まず、Gandhi-HoudeのDifferentiation IVは以下のように定義される。

\[ \begin{split}Z_{jtk}^{\text{Quad,Other}}(X)=\sum_{k\in J_{ft}\setminus\{j\}}d_{jkt\ell}^{2},\\ Z_{jtk}^{\text{Quad,Rival}}(X)=\sum_{k\notin J_{ft}}d_{jkt\ell}^{2}. \end{split} \]

まず\(Z_{jtk}^{\text{Quad,Other}}\)の作成方法を考える． 具体的なイメージを持つべく、以下のように，市場と企業でグループ化したうえで，同じ企業が生産する財\(1,2,3,4\)と，財の特質\(X\)を持ったデータ(data)を考える．

市場	企業	財	特性\(X\)	特性の差の二乗\(d^{2}\)
1	1	1	A	\(\left(A-B\right)^{2}+\left(A-C\right)^{2}+\left(A-D\right)^{2}\)
1	1	2	B	\(\left(B-A\right)^{2}+\left(B-C\right)^{2}+\left(B-D\right)^{2}\)
1	1	3	C	\(\left(C-A\right)^{2}+\left(C-B\right)^{2}+\left(C-D\right)^{2}\)
1	1	4	D	\(\left(D-A\right)^{2}+\left(D-B\right)^{2}+\left(D-C\right)^{2}\)

ここで，財\(1\)についての\(d^{2}\)を考えると，以下のように書ける．

\[ \left(A-B\right)^{2}+\left(A-C\right)^{2}+\left(A-D\right)^{2}=3A^{2}+\left(B^{2}+C^{2}+D^{2}\right)-2A\left(B+C+D\right) \]

これを一般化した形で書くと、

\[ d^{2}=\left(\mathrm{nrow(data)}-1\right)*X^{2}+\left(\mathrm{sum}(X^{2})-X^{2}\right)-2*X*\left(\mathrm{sum}(X)-X\right) \]

となる。

以下では、このアイデアをコードに落とし込んだものを見ていこう。

# まず、マーケット・企業レベルにおける、各製品属性の和と自乗和を計算する。
# ここで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 ) 

# 続いて、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 記述統計と基礎的な分析

まずは、連載のケースに載せたものを再現しよう。

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

イントロダクションにおける、日評自動車が社内で行った分析についてみていこう。 日評自動車はあくまで仮想の自動車会社であるため、ここでは、全車種からランダム選ばれた30車種 を保有する自動車メーカーと考えよう。

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

# 車種の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( Sales, price, hppw, FuelEfficiency, size ) %>% 
  mutate( log_sales = log(Sales),
          log_price = log(price) ) 

その上で、イントロダクションにおける以下の回帰式をOLSで推定しよう。

\[ \log q_j\ =\alpha_0+\alpha_1\ \log\ p_j\ +\alpha_2 size_j+\alpha_3 fuelefficiency_j+\alpha_4 hppw_j \]

ols_intro <-
  fixest::feols(log_sales ~ log_price + hppw + FuelEfficiency + size, data = data_NIPPYO)

# 結果のアウトプット
fixest::etable( ols_intro, 
                se = "hetero",
                signifCode = NA, 
                fitstat = c("r2", "n" ) , 
                dict = c(log_price = "log(価格)",
                         hppw = "馬力／重量",
                         FuelEfficiency = "燃費(キロメートル／ 1 リットル)",
                         size = "サイズ",
                         `(Intercept)` = "定数項"),
                                digits = 2,
                digits.stats = 2,
                depvar = FALSE)

参考：価格と販売台数の散布図

#日評自動車のデータのみ取り出す。
data_graph <-
  data_NIPPYO %>% dplyr::select(price, Sales)

g1 <-
  ggplot2::ggplot(data_graph, aes(price, Sales)) +
  scale_x_continuous(trans = "log10") +
  scale_y_continuous(trans = "log10") +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x)

plot(g1)

4.2 記述統計

データセット全体の記述統計を見てみよう。

data %>% 
  dplyr::select( Sales, price, FuelEfficiency, size, hppw ) %>% 
  summarytools::descr( transpose = TRUE, stats = c("mean", "sd", "q1", "med", "q3"), order = "preserve")

5 ロジットモデルの推定とその応用

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

\[ \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))

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

# BLP操作変数を用いた結果
iv_BLP <-
  fixest::feols(
    logit_share ~ hppw + FuelEfficiency + size | 0 |
      price ~ iv_BLP_own_hppw + iv_BLP_own_FuelEfficiency + iv_BLP_own_size +
      iv_BLP_other_hppw + iv_BLP_other_FuelEfficiency + iv_BLP_other_size,
    data = data
  )

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

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

5.1 1st stage

# BLP操作変数を用いた結果
iv1st_BLP <-
  fixest::feols(
    price ~ hppw + FuelEfficiency + size +
      iv_BLP_own_hppw + iv_BLP_own_FuelEfficiency + iv_BLP_own_size +
      iv_BLP_other_hppw + iv_BLP_other_FuelEfficiency + iv_BLP_other_size,
    data = data
  )

# Differentiation IVを用いた結果
iv1st_GH <-
  fixest::feols(
    price ~ hppw + FuelEfficiency + size + 
      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
  )

# 推定結果をレポートする。
fixest::etable( list( iv1st_BLP, iv1st_GH),  
                se = "hetero",
                fitstat = c("r2", "n") , 
                signifCode = NA, 
                dict = c(price = "自動車価格",
                         hppw = "馬力／重量",
                         FuelEfficiency = "燃費(キロメートル／ 1 リットル)",
                         size = "サイズ",
                         `(Intercept)` = "定数項"),
                digits = 2,
                digits.stats = 2,
                depvar = FALSE)

5.2 自己価格弾力性の計算

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

data %>% 
  dplyr::mutate( own_elas_ols = ols$coefficients["price"]*price*(1-share), 
          own_elas_ivblp = iv_BLP$coefficients["fit_price"]*price*(1-share), 
          own_elas_ivgh  = iv_GH$coefficients["fit_price"]*price*(1-share) ) -> data_elas

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

6 推定結果の応用

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

需要関数の推定結果にもとづいて、需要曲線を書こう。

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

data %>% 
  select( NameID, year, Sales, price, FuelEfficiency, size, hppw, HH, share ) %>% 
  mutate( xi_fit = resid(iv_GH)) -> dt_application

まず、日評自動車が販売する中での特定の車種に着目しよう。ここではNameID 87について考えてみよう。

NameID_target <- 87

dt_application %>% 
  filter( NameID == NameID_target & year == 2016) %>% 
  head()

価格をインプットとして、販売台数を返す関数を考える。

f_share <- function(price_cand, year, modelID_target, dt, estparam){
  
  dt %>% 
    filter( year == 2016) %>% 
    mutate( temp_price = ifelse( NameID == modelID_target, price_cand, price) ) %>% 
    mutate( delta = estparam[1] + estparam[2]*temp_price + estparam[3]*hppw + 
              estparam[4]*FuelEfficiency + estparam[5]*size + xi_fit ) %>% 
    mutate( denom = 1 + sum(exp(delta))) %>% 
    mutate( pred_sales = exp(delta)/denom*HH) %>% 
    filter( NameID == modelID_target ) -> dt_result
  
  return(dt_result$pred_sales)  
  
}

estparam <- iv_GH$coefficients
NameID_target <- 87

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

for ( i in 1:length(pricevec)){
  
  quantvec[i] <- f_share(price_cand = pricevec[i], year = 2016, modelID_target = NameID_target, dt = dt_application, estparam = estparam )
  revenuevec[i] <- pricevec[i]*quantvec[i]
  
}

まず、需要曲線をプロットする。 なお、以下ではggplot内で日本語を利用するためにフォントの設定を行っている。もしエラーが出る場合には、ggplot内のラベルを英語にしてもらいたい。 (なおWindows環境で作業を行っている。)

# フォントの設定。
windowsFonts(
  `Yu Mincho` = windowsFont("Yu Mincho"),
  `Yu Gothic` = windowsFont("Yu Gothic")
)
theme_set(theme(text = element_text(family = "Yu Gothic")))

# 既にテーマに設定したフォントファミリをデフォルトに設定
update_geom_defaults("text", list(family = theme_get()$text$family))
update_geom_defaults("label", list(family = theme_get()$text$family))

# fig_demand <- qplot(quantvec, pricevec*100, xlab = "Sales", ylab = "Price(Unit:10,000JPY)", geom = "line" )
fig_demand <- qplot(quantvec, pricevec*100, xlab = "販売台数", ylab = "価格(単位：万円)", geom = "line" )
plot(fig_demand)

ggsave(file = "demand.png", plot = fig_demand, width = 9, height = 6)
ggsave(file = "demand.pdf", plot = fig_demand, width = 9, height = 6, device = cairo_pdf )

続いて、収入と価格の関係をプロットする。

# pricevecの価格の単位を1万円にする。
fig_rev <- qplot(pricevec*100, revenuevec*100/10000, ylab = "収入(単位：億円)", xlab = "価格(単位：万円)", geom = c("line"), ylim = c(500,1500) )
plot(fig_rev)

ggsave(file = "revenue.png", plot = fig_rev, width = 9, height = 6)
ggsave(file = "revenue.pdf", plot = fig_rev, width = 9, height = 6, device = cairo_pdf)

6.2 収入を最大化にする価格

実際に収入が最大になる価格を、最適化を使って求めてみよう。

# 最適化で利用する関数を定義する。
f_revenue <- function(price_cand, year , modelID_target, dt , estparam){
  
  quantity <- f_share(price_cand, year , modelID_target, dt , estparam  )
  revenue <- price_cand*quantity
  return(revenue)
  
}

# 最適化アルゴリズムを使って、収入を最大にする価格を求める。
result <- optimise( f_revenue, interval = c(0.3, 3), maximum = TRUE, 
          year = 2016, 
          modelID_target = NameID_target, dt = dt_application, estparam = estparam )

print(result)

## $maximum
## [1] 1.813709
## 
## $objective
## [1] 144273.8