class: title-slide-section, center, middle
name: logistics
# イントロダクション


---
## イントロダクション

- 内容：
  1. 事例：きのこの山 VS たけのこの里
  2. コード：`mlogit`パッケージを使った推定


---
class: title-slide-section, center, middle
name: logistics
# 事例：きのこの山 VS たけのこの里

---
## きのこの山 VS たけのこの里

.middle[
.center[
<img src="figuretable/Kinoko.jpg" width="350">
<img src="figuretable/Takenoko.jpg" width="350">
]
]


---
## きのこの山 VS たけのこの里


- アンケートデータ＆離散選択モデルを用いて、「きのこの山」と「たけのこの里」の需要関数を推定しよう。


- 何が知りたいか？
  - 「きのこの山」と「たけのこの里」はどちらが人気か？
  - 「ブランド価値」と「価格」のトレードオフ -> 需要関数・価格弾力性
  - 価格を変えたときに、需要量・収入はどのように変化するか？

---
## アンケート (メインクエスチョン)

- あなたはスーパーマーケットのお菓子コーナーに来ています。目の前に「きのこの山」と「たけのこの里」が並んでいます。「きのこの山」と「たけのこの里」の一箱あたり税込価格が、以下の組み合わせで与えられるとき、あなたにとって最も望ましい選択肢を一つ選んでください。
- 価格の組み合わせ：
  1.	（きのこの山、たけのこの里）＝（200円、200円）
  2.	（きのこの山、たけのこの里）＝（180円、200円）
  3.	（きのこの山、たけのこの里）＝（200円、170円）
  4.	（きのこの山、たけのこの里）＝（220円、200円）
  5.	（きのこの山、たけのこの里）＝（190円、210円）
- 選択肢：
  1.	きのこの山を買う。
  2.	たけのこの里を買う。
  3.	どちらも買わない。

---
## レクチャーで使うデータ

- 2024年春学期「産業組織論」＠早稲田大学政治経済学部で行ったアンケート結果を利用。


- アンケートでは上述の選択肢問題に加えて、消費者属性についての設問も。


---
## 結果を見る前に：選択型コンジョイント分析

- 今から見ていく分析を一般に、**選択型コンジョイント分析**と呼ぶ。

- 仮想状況における選択に関する質問から、消費者の需要・選好を推定することが可能。

- マーケティング、環境経済学などにおいて幅広く使われている手法。

- 限界点：あくまで仮想の状況であるので「実際の行動」を反映するか否かは要注意。

- マーケティングにおけるコンジョイント分析の解説として、[Allenby et al "Economic foundations of conjoint analysis" *Handbook of the Economics of Marketing*](https://www.sciencedirect.com/science/article/abs/pii/S2452261919300024)

---
## データの種類と限界点

- 今回用いたサーベイデータを表明選好データ(Stated preference) と呼ぶ。
  - 仮想のアンケートを用いて、人々の選好に関する情報を収集している。
  
  
- 一方、現実の購買行動データを、顕示選好データ(Revealed preference) と呼ぶ。
  - POS データ、購買履歴データ、など。


- 表明選好データにおいては、常に「その行動が真の行動を反映しているのか？」を考えなければならない。


- 同時に、サーベイ一般の問題として、サーベイ対象者が興味ある母集団全体を反映しているのかについても要検討。


---
class: title-slide-section, center, middle
name: logistics
# 分析の下準備


---
## 【R分析】下準備

```{r setup, include=FALSE, warning=FALSE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
```


```{r}
rm(list = ls())
library("tidyverse")
library("knitr")
library("mlogit") # ロジットモデル推定のためのパッケージ
library("stargazer") # 推定結果の表作成
```

---
## 【R分析】データの読み込み

```{r}
data <- readr::read_csv("data/KinokoTakenokoSurvey_IOSpring2024.csv")

# Cleaning
data <- data %>%
  filter(experience != "4 : 食べたことがない" ) %>% 
  drop_na() # NA(未回答？)の行を除外

#  filter(experience != "4 : 食べたことがない" & gender != "3 : 回答したくない") %>% 


```

```{css, echo = F}
.smalltable table {
  font-size: 18px;     
}
```

---
class: smalltable
## 変数一覧

| 変数          | 説明                         |
|-------------|----------------------------|
| ID          | 回答者のID                     |
| experience  | きのこの山・たけのこの里を最後にいつ食べたか？    |
| Q1          | （きのこの山、たけのこの里）＝（200円、200円） |
| Q2          | （きのこの山、たけのこの里）＝（180円、200円） |
| Q3          | （きのこの山、たけのこの里）＝（200円、170円） |
| Q4          | （きのこの山、たけのこの里）＝（220円、200円） |
| Q5          | （きのこの山、たけのこの里）＝（190円、210円） |
| age         | 年齢                         |
| gender      | 性別                         |
| region      | 出身地                        |
| familyhouse | 実家暮らしか一人暮らしか               |


---
## 【R分析】アンケートのメインの結果

各問における選択肢のシェアを計算する。

```{r}
# 回答者数
N <- length(data$ID)

# 必要なデータを抽出
data %>%
  select(ID, Q1, Q2, Q3, Q4, Q5)  %>% 
  gather(key = Q, value  = choice, Q1, Q2, Q3, Q4, Q5) -> datafig

# データ整形
datafig %>% 
    mutate( Q = ifelse(Q == "Q1", "Q1: (200円, 200円)" ,Q),  
            Q = ifelse(Q == "Q2", "Q2: (180円, 200円)" ,Q),  
            Q = ifelse(Q == "Q3", "Q3: (200円, 170円)" ,Q),  
            Q = ifelse(Q == "Q4", "Q4: (220円, 200円)" ,Q),  
            Q = ifelse(Q == "Q5", "Q5: (190円, 210円)" ,Q) ) -> datafig
```

---
## 【R分析】シェアの計算

```{r}
datafig %>% 
    group_by(Q, choice) %>% 
    tally() %>% 
    mutate( n = n/N) %>% 
    pivot_wider( id_cols = "Q", names_from = "choice", values_from = "n") %>%
    knitr::kable( digits = 2) -> tab

# Latex用出力
datafig %>% 
    group_by(Q, choice) %>% 
    tally() %>% 
    mutate( n = n/N) %>% 
    pivot_wider( id_cols = "Q", names_from = "choice", values_from = "n") %>%
    knitr::kable( digits = 2, "latex") %>%
    kableExtra::save_kable(tab, file = "output/tab_choice_share.tex")

```

---
## 選択肢のシェア

```{r}
tab
```

---
## 【R分析】グラフで出力

帯グラフで結果を出力する。

```{r, eval = TRUE, include = TRUE }
p <- ggplot() +
  geom_bar(data = datafig, aes(x = Q,  fill = as.factor(choice) ),  position = "fill" ) + 
    coord_flip() + 
    theme(legend.position="bottom", axis.title.x=element_blank())+
    scale_fill_discrete(name = "Choice") + theme_bw()
```

```{r, eval = FALSE}
plot(p)
```

---
## アンケートのメインの結果

```{r, echo=FALSE, out.width="40%"}
plot(p)
```

---
## アンケート結果からのポイント

```{r, echo = FALSE}
tab
```

- ポイント１：価格が上がるとシェアが下がる→価格に反応


- ポイント２：価格差が大きくても選ぶ人もいる→きのこ／たけのこのブランド価値


- 離散選択モデルを使って、より詳しく分析してみよう！


---
class: title-slide-section, center, middle
name: logistics
# 多項ロジットモデル


---
## 離散選択モデル：多項ロジットモデル

- 消費者 $i$ が設問 $k$ において、選択肢 $j \in \{Kinoko, Takenoko, outside\}$ から得る効用 
$$
\begin{aligned}
U_{i, k, Kinoko} &=\beta_{Kinoko} - \alpha \cdot p_{k, Kinoko} +\epsilon_{i, k, Kinoko} \\
U_{i, k, Takenoko} &=\beta_{Takenoko}- \alpha \cdot p_{k, Takenoko}+\epsilon_{i, k,Takenoko} \\
U_{i, k, outside} &= \epsilon_{i, k, outside}
\end{aligned}
$$
  - $\beta_{Kinoko}, \beta_{Takenoko}$ : きのこ・たけのこ自体への好み
  - $\alpha$ : 価格への反応度
  - $p_{k,j}$ 設問 $k$ における選択肢 $j$ の価格
  - $\epsilon_{i,k,j}$: ランダムな選好ショック


- 補足：「何も買わない」オプションを「アウトサイド・グッズ」と呼ぶ。


- 推定するパラメタは $\beta_{Kinoko}, \beta_{Takenoko}, \alpha$ の３つ。
  - 全消費者で共通のパラメタとする。後ほど拡張（ランダム係数ロジット）

---
## 多項ロジットモデルにおける選択確率

- 各選択肢の選択確率は
$$P_k(j|\theta ) = \frac{\exp(\beta_{j} - \alpha \cdot p_{j, k})}{1 + \exp(\beta_{Kinoko} - \alpha \cdot p_{Kinoko, k}) + \exp(\beta_{Takenoko} - \alpha \cdot p_{Takenoko, k})}$$

- 尤度関数は
$$L(\theta)= \prod_{i=1}^{N} \prod_{k=1}^{5}P_{k}\left(j=y_{i,k}|\theta\right)$$ 


- 対数尤度関数を、パラメタ $\theta$ に関して最大化する。
  - `mlogit` パッケージを用いる。
  - 自分でプログラムを書く。

---
## 【R分析】推定：`mlogit`パッケージ

- `mlogit`パッケージ：Rにおける離散選択モデル推定のためのパッケージ


- ホームページは[こちら](https://cran.r-project.org/web/packages/mlogit/index.html)
  
  
- 様々な手法をカバー
  - 多項ロジット:[[モデル]](https://cran.r-project.org/web/packages/mlogit/vignettes/c3.rum.html) [[例]](https://cran.r-project.org/web/packages/mlogit/vignettes/e1mlogit.html)
  - ランダム係数ロジット: [[モデル]](https://cran.r-project.org/web/packages/mlogit/vignettes/c5.mxl.html) [[例]](https://cran.r-project.org/web/packages/mlogit/vignettes/e3mxlogit.html)


- (個人的な感想として) パッケージで用いるためのデータ加工が若干独特。
  - データ加工の詳細は[こちら](https://cran.r-project.org/web/packages/mlogit/vignettes/c2.formula.data.html)

---
## 【R分析】`mlogit`パッケージのためのデータ加工

まずは推定用にデータを加工する。

```{r}
# きのこ・たけのこを食べたことない人をDropする。
data %>% 
  filter( experience != "4 : 食べたことがない") -> data_for_estimation

# データをLong形式に加工する。
# 各行は 「回答者ID-設問」という単位になる。
data_for_estimation %>% 
    gather(key = "occasion", value = choice, starts_with("Q")) -> data_for_estimation
```


---
## 【R分析】価格データの作成

- データにマージする、各設問・各選択肢の価格情報を準備する。

```{r}
pricedata <- 
    data.frame( occasion = c("Q1","Q2", "Q3", "Q4", "Q5" ), 
                price_0 = numeric(5), 
                price_1 = c(200, 180, 200, 220, 190 ), 
                price_2 = c(200, 200, 170, 200, 210 ) )
pricedata <- as_tibble(pricedata)
```


- 注意点
  - `price_0` という形で、`変数名`_`選択肢番号` という定義になっている。
  - 0 はアウトサイドグッズ、1はきのこ、2はたけのこ
  - 後ほど`mlogit`で利用するため（後述）


---
## 【R分析】データをマージ

価格データをアンケート結果データに結合する。

```{r}
data_for_estimation %>% 
  left_join(pricedata) %>%
  arrange(ID, occasion) -> data_for_estimation
```

---
## 【R分析】きのこの山・たけのこの里ダミー変数

- ダミー変数を作成する。

```{r}
data_for_estimation %>% 
    mutate( Kinoko_0 = 0, 
            Kinoko_1 = 1, 
            Kinoko_2 = 0, 
            Takenoko_0 = 0, 
            Takenoko_1 = 0, 
            Takenoko_2 = 1   ) %>% 
    arrange(ID, occasion) -> data_for_estimation
```

- 1は「きのこ」、２は「たけのこ」、０は「どちらも買わない」なので、kinoko_1はきのこを買うことから1、takenoko_2はたけのこを買うことから1となる。

---
## 【R分析】選択に関する変数

- 各設問での選択に関する変数`choice`を再定義する。

```{r}
data_for_estimation %>% 
  mutate( choice = case_when( choice == "1 : きのこの山を買う" ~ 1, 
                              choice == "2 : たけのこの里を買う" ~ 2, 
                              choice == "3 : どちらも買わない" ~ 0 )) -> data_for_estimation
```

---
## 【R分析】選択の状況に関する変数の定義

- 現在のデータの各行は、各個人が各設問においてどう行動したかを示している。

- 変数`choiceid`は、このような「選択の状況」を示すIDとなっている。
  - 構築としては、単純に上から下まで連番を振ればよい。

```{r}
data_for_estimation$choiceid <- 1:nrow(data_for_estimation)
# slide_4_code_scratch.RMDのために保存
write_csv(file = "intermediate/data_for_estimation.csv", x = data_for_estimation)

```


---
## 【R分析】`mlogit`用のデータセット

- `mlogit`専用のデータ形式に変更する。

```{r}
datalogit <- dfidx(data = as.data.frame(data_for_estimation),
                   choice = "choice", 
                   varying = 9:17, 
                   sep = "_",
                   idx = list(c("choiceid", "ID")), 
                   idnames = c("chid", "alt"),
                   opposite = c("price"))
```

---
## 【R分析】各引数の説明

- `data`: データセット。
- `choice`: 各選択状況（各行）における選択を示す変数
- `varying`: データセットにおいて、選択肢ごとに異なる値をとる変数。今回は9-17列目。
- `sep`: 上のvaryingな変数は、`変数名`_`選択肢`となっており、この変数名と選択肢を分離する記号を指定する。
- `idx`: 各選択状況を示す変数(`choiceid`)と、各個人を示す変数(`ID`)
を指定する。
- `idnames`: 新しいデータにおいて、各選択状況を示す変数と選択肢をしめす変数の名前
- `opposite`: 効用関数において符号がマイナスにはいるような変数の指定。今回は`price`.

---
## 【R分析】新しいデータセットにおけるインデックス情報

```{r}
head(datalogit$idx)
```


---
## 【R分析】多項ロジットモデルの推定

- 多項ロジットモデルの推定を行う。
```{r}
multilogit <- mlogit( formula = choice ~ price + Kinoko + Takenoko  | 0, 
                      data = datalogit)
```


- 引数：
  - `formula`: 効用の定式化を指定。今回は、価格、きのこダミー、たけのこダミー。なお、`|`の後の`0`はひとまず無視してOK.
  - `data`:上で作成したデータセット


- `stargazer`で推定結果をレポート
```{r, eval = FALSE}
stargazer::stargazer(multilogit, type = "text")
```

---
## 多項ロジットモデルの推定結果

```{r, echo = FALSE}
stargazer::stargazer(multilogit, type = "text")
```


---
class: title-slide-section, center, middle
name: logistics
# 拡張１：消費者属性の追加


---
## 消費者属性

- サーベイでは消費者属性についても質問している。
  - 年齢、性別、出身地、などなど


- きのこの山・たけのこの里への選好や価格の反応度はこれらにも依存するかもしれない。


- これらの要素を取り込んだ多項ロジットモデルの推定を行おう。


---
## 【R分析】消費者属性の記述統計準備

```{r}
# 回答者の属性の比率を表にするためのデータ整形
data %>%
  mutate(gender =　case_when(gender == "1 : 男性" ~ "男性", 
                            gender == "2 : 女性" ~ "女性",
                            gender == "3 : 回答したくない" ~ "未回答"),
         exp    =  case_when(experience == "1 : 過去半年以内" ~ "過去半年以内", 
                             experience == "2 : 過去半年から１年以内" ~ "過去半年から１年以内", 
                             experience == "3 : １年以上前" ~ "１年以上前",
                             experience == "4 : 食べたことがない" ~ "食べたことがない"),
         adult = if_else(age>=20, "成人", "未成年") ) -> temp
```


---

```{r}
temp %>% 
  mutate(region =  case_when(region == "1 : 北海道地方" ~ "関東", 
                             region == "2 : 東北地方" ~ "関東", 
                             region == "3 : 関東地方" ~ "関東" ,
                             region == "4 : 中部地方" ~ "関東",
                             region == "5 : 近畿地方" ~ "関西",
                             region == "6 : 中国地方" ~ "関西",
                             region == "7 : 四国地方" ~ "関西",
                             region == "8 : 九州地方(沖縄含む)" ~ "関西",
                             region == "9 : 海外" ~ "海外"),
         familyhouse = if_else(familyhouse==1,"実家暮らし","実家暮らしでない") )  -> data_atr
```


---

```{r}
# 地域
data_atr %>%
  group_by(region) %>%
  tally() %>%
  mutate(割合 = n/N) %>% 
  select (-n) %>%
  mutate(変数 = c("出身地方","","")) %>%
  rename(属性 = region) %>%
  select(変数, everything()) ->region_fig

# 性別
data_atr %>%
  group_by(gender) %>%
  tally() %>%
  mutate( 割合 = n/N) %>% 
  select(-n) %>%
  mutate(変数 = c("性別","", "")) %>%
  rename(属性 = gender) %>%
  select(変数, everything()) ->gender_fig

```

---

```{r}
# 実家
data_atr %>%
  group_by(familyhouse) %>%
  tally() %>%
  mutate( "割合" = n/N) %>% 
  select(-n) %>%
  mutate("変数" = c("実家暮らしかどうか","")) %>%
  rename(属性 = familyhouse) %>%
  select(変数, everything(),) ->　fam_fig

# 成人
data_atr %>%
  group_by(adult) %>%
  tally() %>%
  mutate( 割合 = n/N) %>% 
  select(-n) %>%
  mutate(変数 = c("成人かどうか","")) %>%
  rename(属性 = adult) %>%
  select(変数, everything()) ->adult_fig

```


---
class: smalltable

```{r}
rbind(gender_fig,adult_fig, fam_fig, region_fig)  %>%
  knitr::kable( digits = 2) -> tab_atr

tab_atr
```

---
## 推定のためのデータセット

```{r}
data_for_estimation %>% 
  mutate( gender =　case_when(gender == "1 : 男性" ~ "0", 
                             gender == "2 : 女性" ~ "female",
                             gender == "3 : 回答したくない" ~ "0"),
          exp    =  case_when(experience == "1 : 過去半年以内" ~ "halfyear", 
                              experience == "2 : 過去半年から１年以内" ~ "half_to_1year", 
                              experience == "3 : １年以上前" ~ "0" ),
          adult = if_else(age>=20, 1, 0),
          region =  case_when(region == "1 : 北海道地方" ~ "_Kanto", 
                              region == "2 : 東北地方" ~ "_Kanto", 
                              region == "3 : 関東地方" ~ "_Kanto" ,
                              region == "4 : 中部地方" ~ "_Kanto",
                              region == "5 : 近畿地方" ~ "Kansai",
                              region == "6 : 中国地方" ~ "Kansai",
                              region == "7 : 四国地方" ~ "Kansai",
                              region == "8 : 九州地方(沖縄含む)" ~ "Kansai",
                              region == "9 : 海外" ~ "Oversea",
          )
  ) -> data_for_estimation
```

---
```{r}
datalogit <- dfidx(data = as.data.frame(data_for_estimation),
                   choice = "choice", 
                   varying = 9:17, 
                   sep = "_",
                   idx = list(c("choiceid", "ID")), 
                   idnames = c("chid", "alt"),
                   opposite = c("price"))
```



---
## 消費者属性入りの定式化

- 価格の係数及びきのこ・たけのこ係数が、消費者属性 $z_i$ に依存すると考える。

$$
\begin{aligned}
U_{i, k, Kinoko} &=\beta_{Kinoko }(z_i) - \alpha(z_i) p_{k, Kinoko} +\epsilon_{i, k, Kinoko} \\
U_{i, k, Takenoko} &=\beta_{Takenoko }(z_i)- \alpha(z_i) p_{k, Takenoko}+\epsilon_{i, k,Takenoko} \\
U_{i, k, outside} &= \epsilon_{i, k, other}
\end{aligned}
$$
- ここで、係数は以下のように与える
$$
\begin{aligned}
\alpha(z_{i})	&=	\alpha_{0}+\sum_{d=1}^{D}\alpha_{d}z_{di} \\
\beta_{j}(z_{i})	&=	\beta_{0,j}+\sum_{d=1}^{D}\beta_{d,j}z_{di}, j\in \{Kinoko, Takenoko \} 
\end{aligned}
$$
  - $z_{di}$ は消費者 $i$ の属性 $d$ (成人ダミーなど)

---
## 【R分析】消費者属性入りの多項ロジットモデル推定

```{r}
with_atr <- formula(choice ~price+Kinoko+Takenoko+
                      price:gender+Takenoko:gender+Kinoko:gender+ 
                      price:familyhouse+Takenoko:familyhouse+Kinoko:familyhouse+
                      price:adult+Takenoko:adult+Kinoko:adult+
                      price:region+Takenoko:region+Kinoko:region | 0 | 0) # 式が長いので事前に格納

ml_consumer_atr <- mlogit(formula = with_atr, 
                          data = datalogit,
                          reflevel = '0')
```

- `gender`: 女性ダミー、 `familyhouse` : 実家ダミー、`adult`: 成人ダミー、 `region`: 地域（関西、海外）

---
## 推定結果の出力

```{r, eval=FALSE}
stargazer::stargazer(multilogit,ml_consumer_atr,
                     single.row=TRUE,type="text",align = TRUE)
```

---

```{r, eval=TRUE, echo=FALSE}
stargazer::stargazer(multilogit,ml_consumer_atr,
                     single.row=TRUE,type="text",align = TRUE) 
```

---
class: title-slide-section, center, middle
name: logistics
# ランダム係数ロジットモデル


---
## 離散選択モデル：ランダム係数ロジットモデル

- もし消費者間で異なった選好パラメタを持っていたらどうなるか？
$$
\begin{aligned}
U_{i, k, Kinoko} &=\beta_{i, Kinoko } - \alpha_{i} \cdot p_{k, Kinoko} +\epsilon_{i, k, Kinoko} \\
U_{i, k, Takenoko} &=\beta_{i, Takenoko }- \alpha_{i} \cdot p_{k, Takenoko}+\epsilon_{i, k,Takenoko} \\
U_{i, k, outside} &= \epsilon_{i, k, other}
\end{aligned}
$$


- ここで、
  -   $\beta_{i, j} \sim N\left(\theta_{j}, \sigma_{j}^{2}\right)$ for $j \in\{$ Kinoko, Takenoko $\}$
  -   $\alpha_{i} \sim  N\left(\theta_{\alpha}, \sigma_{\alpha}^{2}\right) .$


- これらの分布のパラメタを推定する。
  - 技術的な詳細は補足資料を参照(少し難しめ)


---
## 【R分析】ランダム係数ロジットモデルの推定

```{r, eval = FALSE}
# 乱数固定
set.seed(100)

# 推定
rcdclogit <- mlogit(choice ~ price + Kinoko + Takenoko | 0, 
                    data = datalogit, 
                    panel = TRUE, 
                    rpar = c(price = "n", Kinoko = "n", Takenoko = "n") , 
                    R = 100, 
                    correlation = FALSE)
```

- 引数：
  - `panel`: 今回は各消費者について複数回の選択を観察しているのでパネルデータ構造となっている。
  - `rpar`: ランダム係数の定式化。`ln`は対数正規分布、`n`は正規分布
  - `R`: モンテカルロ積分における乱数のドロー数。
  - `correlation`: ランダム係数の間の相関構造。多くの場合は独立(つまり`FALSE`)

---
## 推定時間
```{r, eval = TRUE}
# 乱数固定
set.seed(100)

system.time( 
rcdclogit <- mlogit(choice ~ price + Kinoko + Takenoko | 0, 
                    data = datalogit, 
                    panel = TRUE, 
                    rpar = c(price = "n", Kinoko = "n", Takenoko = "n") , 
                    R = 100, 
                    correlation = FALSE)
)
```

---
## 推定結果

```{r, eval = FALSE}
stargazer::stargazer(multilogit, rcdclogit, type="text")
```

---
## 【参考】乱数のドロー数

- 乱数のドロー数を便宜上100としているが、これは少なすぎる。


- ドロー数が多いとより正確。ただし、計算時間(推定時間)が長くなる。


- 推定結果がドロー数によってあまり変化しなくなるくらい大きい値にするのが重要。


- なお、ガウス求積など、正確な値を得つつ、計算負荷を下げる方法もある。(興味ある人は個人的に)





---
## 推定結果(ロジットとランダム係数ロジット)

```{r, echo = FALSE}
stargazer::stargazer(multilogit, rcdclogit, type="text", dep.var.caption = "")
```


---
## (おまけ) 推定結果の外部ファイル保存

```{r, echo = TRUE}
writeLines(capture.output(
    stargazer::stargazer(multilogit, rcdclogit, type="latex", dep.var.caption = "")
), "output/tab_estimated_parameter.tex")
```

---
## 推定したランダム係数の分布

- まず、 $\alpha_{i,Kinoko} - \alpha_{i,Takenoko}$ という選好パラメタの差をプロットする。

- 平均及び分散は、
$$
\begin{aligned}
\mathrm{E}[\beta_{i,Kinoko} - \beta_{i,Takenoko}] &= \mathrm{E}[\beta_{i,Kinoko}]  - \mathrm{E}[\beta_{i,Takenoko}] \\
\mathrm{Var}[\beta_{i,Kinoko} - \beta_{i,Takenoko}] &= \mathrm{Var}[\beta_{i,Kinoko}]  + \mathrm{Var}[\beta_{i,Takenoko}]
\end{aligned}
$$
なお、 $\beta_{i,Kinoko}$ と $\beta_{i,Takenoko}$ は独立という仮定を利用している。

---
## 【R分析】ランダム係数の分布

- `mlogit::rpar`関数を使うと、ランダム係数の分布パラメタを取得できる。

```{r}
dist_Kinoko = rpar(rcdclogit, 'Kinoko')
dist_Takenoko = rpar(rcdclogit, 'Takenoko')
dist_price = rpar(rcdclogit, 'price')
```


---
## 【R分析】ランダム係数の分布

- 推定値から、 $\beta_{i,Kinoko} - \beta_{i, Takenoko}$ の分布プロットを作成する。

```{r}
p1 <- ggplot(data = data.frame(x = c(-25, 25)), aes(x)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = dist_Kinoko$mean -dist_Takenoko$mean , 
                            sd = sqrt(dist_Kinoko$sigma^2 + dist_Takenoko$sigma^2) ) ) +
  ylab("") + theme_bw() +  ggtitle("Preference Kinoko over Takenoko")
ggsave(file = "output/fig_dist_K_over_T.pdf", plot = p1)
```

```{r, eval = FALSE}
plot(p1)
```


---
## ランダム係数の分布： $\beta_{i,Kinoko} - \beta_{i, Takenoko}$

```{r, echo = FALSE, fig.align = "center"}
plot(p1)
```

---
## 【R分析】価格係数の分布

- 価格係数 $\alpha_{i} \sim N\left(\theta_{\alpha}, \sigma_{\alpha}^{2}\right)$ についてプロットする。

```{r}
p2 <- ggplot(data = data.frame(x = c(0, 0.4) ), aes(x)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = dist_price$mean, sd = abs(dist_price$sigma)), colour = "red" ) + 
  ylab("") + theme_bw() + 
  scale_y_continuous(breaks = NULL) + ggtitle("Alpha (price coefficient)")
ggsave(file = "output/fig_dist_beta.pdf", plot = p2)
```

```{r, eval = FALSE}
plot(p2)
```

---
## 価格係数の分布
```{r, echo = FALSE, fig.align = "center"}
plot(p2)
```


---
## 支払い意思額(Willingness-to-pay, WTP)

- 各消費者のきのこ・たけのこへの支払い意志額は $ \beta_{i,Kinoko}/ \alpha_i$ で表される。


- 支払い意思額の分布をプロットしよう。


---
## WTPの計算

```{r}
# 乱数固定
set.seed(100)

# パラメタについての乱数を発生
dt_brand = 
  tibble( beta_kinoko_draw = rnorm(30000, mean = dist_Kinoko$mean, sd = abs(dist_Kinoko$sigma) ), 
          beta_takenoko_draw = rnorm(30000, mean = dist_Takenoko$mean, sd = abs(dist_Takenoko$sigma) ),  
          alpha_draw =rnorm(30000, mean = dist_price$mean, sd = abs(dist_price$sigma) ) )

# 支払い意思額の計算
dt_brand %>% 
  mutate( brand_Kinoko = beta_kinoko_draw / alpha_draw, 
          brand_Takenoko = beta_takenoko_draw / alpha_draw,
          love_Kinoko = ifelse(beta_kinoko_draw > beta_takenoko_draw, 1, 0 )) -> dt_brand 

```

---
## 続・WTPの計算

```{r}
dt_brand %>% 
  select( brand_Kinoko, brand_Takenoko) %>% 
  rename( WTP_Kinoko = brand_Kinoko, 
         WTP_Takenoko = brand_Takenoko) %>% 
  pivot_longer( cols = everything(), names_to = "Name", values_to = "WTP") -> dt

dt %>% 
  ggplot( aes(x = WTP, colour  = Name)) + theme_bw() + 
  geom_density( ) -> fig_brand

ggsave(filename = "output/brand_value.pdf", plot = fig_brand, width = 12, height = 6)

```


---
## WTPのプロット

```{r}
plot(fig_brand)
```


---
class: title-slide-section, center, middle
name: logistics
# 需要曲線の構築

---
## 需要曲線のプロット

- ランダム係数ロジットモデル(消費者属性なし)の推定結果にもとづいて、需要曲線をプロットする。


- 解釈：推定しているものは消費者の選択確率のモデル。しかしながら、選択確率＝市場シェアと解釈可能。


- 選択確率(市場シェア)に潜在的な消費者数をかけ合わせると、市場全体の需要が得られる。


---
## 設定

- 1000人の消費者がいると考える。

- 各消費者の選好パラメタは推定された分布に従っており、ランダムに決まっている。

- 各消費者の購買確率を計算し、それを足し合わせることで、全体の需要を予測する。
  
---
## 【R分析】下準備

- 推定した分布パラメタにもとづいて、個々の消費者のパラメタをドローする。

```{r}
# 以下では乱数を用いるので乱数のシードを固定
set.seed(101)

# 消費者数
R = 1000

# 消費者のパラメタを、乱数によって発生させる。
beta_Kinoko_vec = rnorm(n = R, mean = dist_Kinoko$mean, sd = abs(dist_Kinoko$sigma) )
beta_Takenoko_vec = rnorm(n = R, mean = dist_Takenoko$mean, sd = abs(dist_Takenoko$sigma) )
alpha_vec  = rnorm(n = R, mean = dist_price$mean, sd = abs(dist_price$sigma) )
```

---
## 乱数について

- 今、1000人の消費者が異なったパラメタをもっており、そのパラメタは推定した分布に従う。


- 1000人のパラメタについて、その分布に従う**乱数**として発生させる。
  - `rnorm` と `rlnorm` が正規乱数・対数正規乱数
  
  
- 【重要】乱数を発生させる際には、乱数のシードを固定する。
  - `set.seed(101)` に相当する。
  - これをしないと、分析を回す度に異なった乱数が生成される。(再現性の観点から重要)


---
## 【R分析】ロジット確率計算のための関数

- パラメタ $\beta, \beta_{Kinoko}, \beta_{Takenoko}$ と 価格を与えることで、選択確率を計算する。

```{r}
f_logit_prob <- function(beta_Kinoko, beta_Takenoko, beta, price_Kinoko, price_Takenoko){

  util_Kinoko <- beta_Kinoko - beta*price_Kinoko
  util_Takenoko <- beta_Takenoko - beta*price_Takenoko

  prob_Kinoko <- exp( util_Kinoko ) / ( 1 + exp( util_Kinoko ) + exp( util_Takenoko ))
  prob_Takenoko <- exp( util_Takenoko ) / ( 1 + exp( util_Kinoko ) + exp( util_Takenoko ))
  prob_Other <- 1 - (prob_Kinoko + prob_Takenoko)

  return( cbind(prob_Kinoko, prob_Takenoko, prob_Other))
}
```


---
## 【R分析】需要関数

```{r}
f_demand <- function( beta_Kinoko_vec, beta_Takenoko_vec, alpha_vec, price_Kinoko, price_Takenoko ){
    
    R = length(beta_Kinoko_vec)  # Number of consumers
    
    # 結果を保存するベクトルを事前に準備
    prob_Kinoko = numeric(R)
    prob_Takenoko = numeric(R)
    prob_Other = numeric(R)

    for (r in 1:R){
      result = f_logit_prob( beta_Kinoko_vec[r], beta_Takenoko_vec[r], alpha_vec[r], price_Kinoko, price_Takenoko)
      prob_Kinoko[r] = result[1]
      prob_Takenoko[r] = result[2]
      prob_Other[r] = result[3]
    }

    # 選択確率をすべての消費者について足し合わせて、各オプションの需要を得る。
    return( c(sum(prob_Kinoko), sum(prob_Takenoko), sum(prob_Other)))
}
```

---
## 【R分析】きのこの需要関数

- たけのこの価格を２００円に固定したときの、きのこの需要関数を求めてみる。

```{r}
price_Takenoko = 200 

# きのこの価格を100円から250円まで動かす。
price_vec = seq(from = 100, to = 250, by = 5)
kinoko_vec = numeric(length(price_vec))

# ループで需要を計算
for ( i in 1:length(price_vec)){
    result <- f_demand( beta_Kinoko_vec, beta_Takenoko_vec, 
                        alpha_vec, price_vec[i], price_Takenoko = 200 )
    kinoko_vec[i] <- result[1]
}

# データフレームに結果を保存
data_demand_kinoko = tibble( price = price_vec, 
                             demand = kinoko_vec,
                             revenue = price_vec*kinoko_vec)
```

---
## 【R分析】プロット

```{r}
fig_demand <- ggplot(data = data_demand_kinoko, aes(x = demand, y = price) ) +
  geom_line() + theme_bw() + 
  xlab("Demand (Share of Kinoko)") + 
  ylab("Price of Kinoko in JPY") + 
  ggtitle("Demand Curve of Kinoko when Takenoko's price = JPY 200")
ggsave(filename = "output/demand_kinoko.pdf", plot = fig_demand)
```

```{r, eval = FALSE}
plot(fig_demand)
```

---
## 推定された需要関数

```{r, echo = FALSE, fig.align = "center"}
plot(fig_demand)
```

---
## 【R分析】収入関数のプロット

```{r}
fig_rev <- ggplot(data = data_demand_kinoko, aes(x = price, y = revenue) ) +
  geom_line() + theme_bw() + 
  xlab("Price of Kinoko in JPY") + 
  ylab("Revenue of Kinoko in JPY") + 
  ggtitle("Revenue Curve of Kinoko when Takenoko's price = JPY 200")
ggsave(filename = "output/revenue_kinoko.pdf", plot = fig_rev)
```

```{r, eval = FALSE, fig.align = "center"}
plot(fig_rev)
```

---
## 推定された収入関数

```{r, echo = FALSE, fig.align = "center"}
plot(fig_rev)
```


